Commit 1e3203f7 authored by David Rohr's avatar David Rohr
Browse files

convert indent spaces to tabs, unify indentation

parent 87ad9163
......@@ -54,14 +54,14 @@ int bioem::configure(int ac, char* av[])
//*************************************************************************************
/*** Inizialzing default variables ***/
std::string infile,modelfile,mapfile;
Model.readPDB=false;
param.writeAngles=false;
std::string infile, modelfile, mapfile;
Model.readPDB = false;
param.writeAngles = false;
RefMap.dumpMap = false;
RefMap.loadMap = false;
RefMap.readMRC = false;
RefMap.readMRC = false;
RefMap.readMultMRC = false;
//*************************************************************************************
cout << "+++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
......@@ -76,8 +76,8 @@ int bioem::configure(int ac, char* av[])
("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
("Particlesfile", po::value< std::string>(), "(Mandatory) Name of paricles file")
("ReadPDB", "(Optional) If reading model file in PDB format")
("ReadMRC", "(Optional) If reading particle file in MRC format")
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
("ReadMRC", "(Optional) If reading particle file in MRC format")
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
("DumpMaps", "(Optional) Dump maps after they were red from maps file")
("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
("help", "(Optional) Produce help message")
......@@ -88,8 +88,8 @@ int bioem::configure(int ac, char* av[])
p.add("Modelfile", -1);
p.add("Particlesfile", -1);
p.add("ReadPDB", -1);
p.add("ReadMRC", -1);
p.add("ReadMultipleMRC", -1);
p.add("ReadMRC", -1);
p.add("ReadMultipleMRC", -1);
p.add("DumpMaps", -1);
p.add("LoadMapDump", -1);
......@@ -99,9 +99,9 @@ int bioem::configure(int ac, char* av[])
po::notify(vm);
if((ac < 6)) {
std::cout << desc << std::endl;
std::cout << desc << std::endl;
exit(1);
return 0;
return 0;
}
if (vm.count("help")) {
cout << "Usage: options_description [options]\n";
......@@ -112,32 +112,32 @@ int bioem::configure(int ac, char* av[])
if (vm.count("Inputfile"))
{
cout << "Input file is: ";
cout << vm["Inputfile"].as< std::string >()<< "\n";
infile=vm["Inputfile"].as< std::string >();
cout << vm["Inputfile"].as< std::string >() << "\n";
infile = vm["Inputfile"].as< std::string >();
}
if (vm.count("Modelfile"))
{
cout << "Model file is: "
<< vm["Modelfile"].as< std::string >() << "\n";
modelfile=vm["Modelfile"].as< std::string >();
modelfile = vm["Modelfile"].as< std::string >();
}
if (vm.count("ReadPDB"))
{
cout << "Reading model file in PDB format.\n";
Model.readPDB=true;
Model.readPDB = true;
}
if (vm.count("ReadMRC"))
{
cout << "Reading particle file in MRC format.\n";
RefMap.readMRC=true;
}
{
cout << "Reading particle file in MRC format.\n";
RefMap.readMRC = true;
}
if (vm.count("ReadMultipleMRC"))
{
cout << "Reading Multiple MRCs.\n";
RefMap.readMultMRC=true;
}
{
cout << "Reading Multiple MRCs.\n";
RefMap.readMultMRC = true;
}
if (vm.count("DumpMaps"))
{
......@@ -155,7 +155,7 @@ int bioem::configure(int ac, char* av[])
{
cout << "Paricle file is: "
<< vm["Particlesfile"].as< std::string >() << "\n";
mapfile=vm["Particlesfile"].as< std::string >();
mapfile = vm["Particlesfile"].as< std::string >();
}
}
catch(std::exception& e)
......@@ -203,15 +203,15 @@ int bioem::precalculate()
// Generating Grids of orientations
param.CalculateGridsParam();
myfloat_t sum,sumsquare;
myfloat_t sum, sumsquare;
//Precalculating cross-correlations of maps
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
{
calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
calcross_cor(RefMap.Ref[iRefMap], sum, sumsquare);
//Storing Crosscorrelations in Map class
RefMap.sum_RefMap[iRefMap]=sum;
RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
RefMap.sum_RefMap[iRefMap] = sum;
RefMap.sumsquare_RefMap[iRefMap] = sumsquare;
}
// Precalculating CTF Kernels stored in class Param
......@@ -238,13 +238,13 @@ int bioem::run()
// Inizialzing Probabilites to zero and constant to -Infinity
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
pProb[iRefMap].Total=0.0;
pProb[iRefMap].Constoadd=-9999999;
pProb[iRefMap].max_prob=-9999999;
pProb[iRefMap].Total = 0.0;
pProb[iRefMap].Constoadd = -9999999;
pProb[iRefMap].max_prob = -9999999;
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
{
pProb[iRefMap].forAngles[iOrient]=0.0;
pProb[iRefMap].ConstAngle[iOrient]=-99999999;
pProb[iRefMap].forAngles[iOrient] = 0.0;
pProb[iRefMap].ConstAngle[iOrient] = -99999999;
}
}
//**************************************************************************************
......@@ -256,11 +256,11 @@ int bioem::run()
mycomplex_t* proj_mapFFT;
bioem_map conv_map;
mycomplex_t* conv_mapFFT;
myfloat_t sumCONV,sumsquareCONV;
myfloat_t sumCONV, sumsquareCONV;
//allocating fftw_complex vector
proj_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
conv_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
proj_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
HighResTimer timer;
......@@ -322,7 +322,7 @@ int bioem::run()
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
//**** Total Probability ***
outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total)+pProb[iRefMap].Constoadd+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n";
outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total) + pProb[iRefMap].Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n";
outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
......@@ -344,7 +344,7 @@ int bioem::run()
{
for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
{
angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut])+pProb[iRefMap].ConstAngle[iProjectionOut]+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut]) + pProb[iRefMap].ConstAngle[iProjectionOut] + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
}
}
......@@ -364,7 +364,7 @@ int bioem::run()
if (param.refCTF)
{
delete[] param.refCTF;
param.refCTF =NULL;
param.refCTF = NULL;
}
if(RefMap.RefMapsFFT)
......@@ -378,17 +378,17 @@ int bioem::run()
int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
{
//***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps *****
//***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps *****
if (FFTAlgo) //IF using the fft algorithm
{
#pragma omp parallel
#pragma omp parallel
{
mycomplex_t *localCCT;
myfloat_t *lCC;
localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
lCC= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
localCCT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
lCC = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
const int num_threads = omp_get_num_threads();
const int thread_id = omp_get_thread_num();
......@@ -398,7 +398,7 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
for (int iRefMap = iStart; iRefMap < iEnd; iRefMap ++)
{
calculateCCFFT(iRefMap,iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT,lCC);
calculateCCFFT(iRefMap, iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT, lCC);
}
myfftw_free(localCCT);
myfftw_free(lCC);
......@@ -406,34 +406,34 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
}
else
{
#pragma omp parallel for
#pragma omp parallel for
for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
compareRefMapShifted < -1 > (iRefMap, iProjectionOut, iConv, conv_map, pProb, param.param_device, RefMap);
}
}
return(0);
}
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC)
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
{
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize];
for(int i = 0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
{
localCCT[i][0] = localConvFFT[i][0] * RefMapFFT[i][0] + localConvFFT[i][1] * RefMapFFT[i][1];
localCCT[i][1] = localConvFFT[i][1] * RefMapFFT[i][0] - localConvFFT[i][0] * RefMapFFT[i][1];
}
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,localCCT,lCC);
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
}
int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
{
//**************************************************************************************
//**** BioEM Create Projection routine in Euler angle predefined grid****************
......@@ -442,43 +442,43 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
myfloat3_t RotatedPointsModel[Model.nPointsModel];
myfloat_t rotmat[3][3];
myfloat_t alpha, gam,beta;
myfloat_t alpha, gam, beta;
myfloat_t* localproj;
localproj= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
memset(localproj,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
localproj = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
alpha=param.angles[iMap].pos[0];
beta=param.angles[iMap].pos[1];
gam=param.angles[iMap].pos[2];
alpha = param.angles[iMap].pos[0];
beta = param.angles[iMap].pos[1];
gam = param.angles[iMap].pos[2];
//**** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
//********** Creat Rotation with pre-defiend grid of orientations**********
rotmat[0][0]=cos(gam)*cos(alpha)-cos(beta)*sin(alpha)*sin(gam);
rotmat[0][1]=cos(gam)*sin(alpha)+cos(beta)*cos(alpha)*sin(gam);
rotmat[0][2]=sin(gam)*sin(beta);
rotmat[1][0]=-sin(gam)*cos(alpha)-cos(beta)*sin(alpha)*cos(gam);
rotmat[1][1]=-sin(gam)*sin(alpha)+cos(beta)*cos(alpha)*cos(gam);
rotmat[1][2]=cos(gam)*sin(beta);
rotmat[2][0]=sin(beta)*sin(alpha);
rotmat[2][1]=-sin(beta)*cos(alpha);
rotmat[2][2]=cos(beta);
for(int n=0; n< Model.nPointsModel; n++)
rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam);
rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam);
rotmat[0][2] = sin(gam) * sin(beta);
rotmat[1][0] = -sin(gam) * cos(alpha) - cos(beta) * sin(alpha) * cos(gam);
rotmat[1][1] = -sin(gam) * sin(alpha) + cos(beta) * cos(alpha) * cos(gam);
rotmat[1][2] = cos(gam) * sin(beta);
rotmat[2][0] = sin(beta) * sin(alpha);
rotmat[2][1] = -sin(beta) * cos(alpha);
rotmat[2][2] = cos(beta);
for(int n = 0; n < Model.nPointsModel; n++)
{
RotatedPointsModel[n].pos[0]=0.0;
RotatedPointsModel[n].pos[1]=0.0;
RotatedPointsModel[n].pos[2]=0.0;
RotatedPointsModel[n].pos[0] = 0.0;
RotatedPointsModel[n].pos[1] = 0.0;
RotatedPointsModel[n].pos[2] = 0.0;
}
for(int n=0; n< Model.nPointsModel; n++)
for(int n = 0; n < Model.nPointsModel; n++)
{
for(int k=0; k< 3; k++)
for(int k = 0; k < 3; k++)
{
for(int j=0; j< 3; j++)
for(int j = 0; j < 3; j++)
{
RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.PointsModel[n].pos[j];
}
}
}
......@@ -486,41 +486,41 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
int i, j;
//************ Projection over the Z axis********************
for(int n=0; n< Model.nPointsModel; n++)
for(int n = 0; n < Model.nPointsModel; n++)
{
//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
i=floor(RotatedPointsModel[n].pos[0]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
localproj[i*param.param_device.NumberPixels+j]+=Model.densityPointsModel[n]/Model.NormDen;
localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
}
//**** Output Just to check****
if(iMap==10)
if(iMap == 10)
{
ofstream myexamplemap;
ofstream myexampleRot;
myexamplemap.open ("MAP_i10");
myexampleRot.open ("Rot_i10");
myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
for(int k=0; k<param.param_device.NumberPixels; k++)
for(int k = 0; k < param.param_device.NumberPixels; k++)
{
for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j];
for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
}
myexamplemap << " \n";
for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
myexamplemap.close();
myexampleRot.close();
}
//***** Converting projection to Fourier Space for Convolution later with kernel****
//********** Omp Critical is necessary with FFTW*******
myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localproj,mapFFT);
myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
return(0);
}
int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
{
//**************************************************************************************
//**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
......@@ -529,14 +529,14 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
//**************************************************************************************
myfloat_t* localconvFFT;
localconvFFT= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t)*param.param_device.NumberPixels*param.param_device.NumberPixels);
localconvFFT = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
mycomplex_t* tmp;
tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
//**** Multiplying FFTmap with corresponding kernel ****
const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize];
for(int i=0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
{
localmultFFT[i][0] = lproj[i][0] * refCTF[i][0] + lproj[i][1] * refCTF[i][1];
localmultFFT[i][1] = lproj[i][1] * refCTF[i][0] - lproj[i][0] * refCTF[i][1];
......@@ -547,20 +547,20 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
//**** Bringing convoluted Map to real Space ****
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,tmp,localconvFFT);
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, localconvFFT);
//****Asigning convolution fftw_complex to bioem_map ****
for(int i=0; i < param.param_device.NumberPixels ; i++ )
for(int i = 0; i < param.param_device.NumberPixels ; i++ )
{
for(int j=0; j < param.param_device.NumberPixels ; j++ )
for(int j = 0; j < param.param_device.NumberPixels ; j++ )
{
Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j];
Mapconv.points[i][j] = localconvFFT[i * param.param_device.NumberPixels + j];
}
}
//*** Calculating Cross-correlations of cal-convoluted map with its self *****
sumC=0;
sumsquareC=0;
sumC = 0;
sumsquareC = 0;
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
{
sumC += localconvFFT[i];
......@@ -580,12 +580,12 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
return(0);
}
int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
int bioem::calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare)
{
//*********************** Routine to calculate Cross correlations***********************
sum=0.0;
sumsquare=0.0;
sum = 0.0;
sumsquare = 0.0;
for (int i = 0; i < param.param_device.NumberPixels; i++)
{
for (int j = 0; j < param.param_device.NumberPixels; j++)
......@@ -593,7 +593,7 @@ int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
// Calculate Sum of pixels
sum += localmap.points[i][j];
// Calculate Sum of pixels squared
sumsquare += localmap.points[i][j]*localmap.points[i][j];
sumsquare += localmap.points[i][j] * localmap.points[i][j];
}
}
return(0);
......
......@@ -61,8 +61,8 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
const myfloat_t ForLogProb = (sumsquare * param.Ntotpi - sum * sum);
// Products of different cross-correlations (first element in formula)
const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare-crossproMapConv * crossproMapConv) +
2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare;
const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare - crossproMapConv * crossproMapConv) +
2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare;
//******* Calculating log of Prob*********/
// As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1);
......@@ -70,7 +70,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
return(logpro);
}
__device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
{
/********************************************************/
/*********** Calculates the BioEM probability ***********/
......@@ -80,8 +80,8 @@ __device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat
//update_prob<-1>(logpro, iRefMap, iOrient, iConv, disx, disy, pProb);
//GCC is too stupid to inline properly, so the code is copied here
if(pProb[iRefMap].Constoadd < logpro)
{
if(pProb[iRefMap].Constoadd < logpro)
{
pProb[iRefMap].Total = pProb[iRefMap].Total * exp(-logpro + pProb[iRefMap].Constoadd);
pProb[iRefMap].Constoadd = logpro;
}
......@@ -106,33 +106,33 @@ __device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat
__device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
{
for (int cent_x = 0; cent_x <= param.maxDisplaceCenter; cent_x=cent_x+param.GridSpaceCenter)
for (int cent_x = 0; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
{
for (int cent_y = 0; cent_y <= param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
for (int cent_y = 0; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
{
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
}
for (int cent_y = param.NumberPixels-param.maxDisplaceCenter; cent_y < param.NumberPixels; cent_y=cent_y+param.GridSpaceCenter)
for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y < param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
{
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
}
}
for (int cent_x = param.NumberPixels-param.maxDisplaceCenter; cent_x < param.NumberPixels; cent_x=cent_x+param.GridSpaceCenter)
for (int cent_x = param.NumberPixels - param.maxDisplaceCenter; cent_x < param.NumberPixels; cent_x = cent_x + param.GridSpaceCenter)
{
for (int cent_y = 0; cent_y < param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
for (int cent_y = 0; cent_y < param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
{
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
}
for (int cent_y = param.NumberPixels-param.maxDisplaceCenter; cent_y <= param.NumberPixels; cent_y=cent_y+param.GridSpaceCenter)
for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y <= param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
{
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
}
}
}
template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap,
const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true)
const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true)
{
/**************************************************************************************/
/********************** Calculating BioEM Probability ********************************/
......@@ -145,9 +145,9 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
myfloat_t sum, sumsquare, crossproMapConv;
__m128 sum_v = _mm_setzero_ps(), sumsquare_v = _mm_setzero_ps(), cross_v = _mm_setzero_ps(), d1, d2;
#else
myfloat_t sum=0.0;
myfloat_t sumsquare=0.0;
myfloat_t crossproMapConv=0.0;
myfloat_t sum = 0.0;
myfloat_t sumsquare = 0.0;
myfloat_t crossproMapConv = 0.0;
#endif
/****** Loop over Pixels to calculate dot product and cross-correlations of displaced Ref Conv. Map***/
myfloat_t logpro;
......@@ -182,7 +182,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
const float* ptr2 = RefMap.getp(iRefMap, i, jStart);
int j;
const int count = jEnd - jStart;
for (j = 0;j <= count - 4;j += 4)
for (j = 0; j <= count - 4; j += 4)
{
d1 = _mm_loadu_ps(ptr1);
d2 = _mm_loadu_ps(ptr2);
......@@ -201,7 +201,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
// Crosscorrelation of calculated displaced map
sum += pointMap;
// Calculate Sum of pixels squared
sumsquare += pointMap*pointMap;
sumsquare += pointMap * pointMap;
}
#endif
}
......@@ -380,18 +380,18 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
else
#endif
/***** Summing & Storing total/Orientation Probabilites for each map ************/
/***** Summing & Storing total/Orientation Probabilites for each map ************/
{
update_prob<-1>(logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
}
}
template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMapShifted(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap)
{
for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x=cent_x+param.GridSpaceCenter)
for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
{