From 65ca85371255dcdf720660f2d8e3844e43388d33 Mon Sep 17 00:00:00 2001
From: David Rohr <drohr@jwdt.org>
Date: Sat, 19 Apr 2014 16:55:36 +0200
Subject: [PATCH] convert indent spaces to tabs, unify indentation
---
bioem.cpp | 160 +++++++++++++++++++++----------------------
bioem_algorithm.h | 50 +++++++-------
bioem_cuda.cu | 32 ++++-----
include/bioem.h | 4 +-
include/map.h | 8 +--
include/model.h | 2 +-
main.cpp | 4 +-
map.cpp | 58 ++++++++--------
model.cpp | 160 +++++++++++++++++++++----------------------
param.cpp | 170 +++++++++++++++++++++++-----------------------
10 files changed, 324 insertions(+), 324 deletions(-)
diff --git a/bioem.cpp b/bioem.cpp
index 9a1e3db..7c50a53 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -54,9 +54,9 @@ 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;
param.dumpMap = false;
param.loadMap = false;
@@ -104,20 +104,20 @@ 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("DumpMaps"))
@@ -136,7 +136,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)
@@ -184,15 +184,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.getmap(iRefMap),sum,sumsquare);
+ calcross_cor(RefMap.getmap(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
@@ -219,13 +219,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;
}
}
/**************************************************************************************/
@@ -237,11 +237,11 @@ int bioem::run()
mycomplex_t* proj_mapFFT;
myfloat_t* conv_map = new myfloat_t[param.param_device.NumberPixels * param.param_device.NumberPixels];
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;
@@ -304,7 +304,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: ";
@@ -326,7 +326,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";
}
}
@@ -346,7 +346,7 @@ int bioem::run()
if (param.refCTF)
{
delete[] param.refCTF;
- param.refCTF =NULL;
+ param.refCTF = NULL;
}
RefMap.freePointers();
@@ -357,12 +357,12 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* conv_m
{
if (FFTAlgo)
{
-#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();
@@ -372,7 +372,7 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* 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);
@@ -380,30 +380,30 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* 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)
{
const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
- 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****************
@@ -412,43 +412,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];
}
}
}
@@ -456,41 +456,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,myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
+int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
{
/**************************************************************************************/
/**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
@@ -499,14 +499,14 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,m
/**************************************************************************************/
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 = ¶m.refCTF[iConv * param.FFTMapSize];
- 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];
@@ -517,20 +517,20 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,m
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[i*param.param_device.NumberPixels+j] = localconvFFT[i*param.param_device.NumberPixels+j];
+ Mapconv[i * param.param_device.NumberPixels + 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];
@@ -550,20 +550,20 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,m
return(0);
}
-int bioem::calcross_cor(myfloat_t* localmap,myfloat_t& sum,myfloat_t& sumsquare)
+int bioem::calcross_cor(myfloat_t* 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++)
{
// Calculate Sum of pixels
- sum += localmap[i*param.param_device.NumberPixels+j];
+ sum += localmap[i * param.param_device.NumberPixels + j];
// Calculate Sum of pixels squared
- sumsquare += localmap[i*param.param_device.NumberPixels+j] * localmap[i*param.param_device.NumberPixels+j];
+ sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
}
}
return(0);
diff --git a/bioem_algorithm.h b/bioem_algorithm.h
index f62127f..c2754e5 100644
--- a/bioem_algorithm.h
+++ b/bioem_algorithm.h
@@ -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, 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_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), 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), 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), 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), param.NumberPixels - cent_x, param.NumberPixels - cent_y, pProb, param, RefMap);
}
}
}
template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* 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 myfloat_t* 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)
{
- for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
+ for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
{
compareRefMap<GPUAlgo>(iRefMap, iOrient, iConv, Mapconv, pProb, param, RefMap, cent_x, cent_y);
}
diff --git a/bioem_cuda.cu b/bioem_cuda.cu
index 60c34e5..26b6c6a 100644
--- a/bioem_cuda.cu
+++ b/bioem_cuda.cu
@@ -37,7 +37,7 @@ __global__ void compareRefMap_kernel(const int iOrient, const int iConv, const m
const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
if (iRefMap < maxRef)
{
- compareRefMap<0>(iRefMap,iOrient,iConv, pMap, pProb, param, *RefMap, cent_x, cent_y);
+ compareRefMap<0>(iRefMap, iOrient, iConv, pMap, pProb, param, *RefMap, cent_x, cent_y);
}
}
@@ -46,7 +46,7 @@ __global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv,
const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
if (iRefMap < maxRef)
{
- compareRefMapShifted<1>(iRefMap,iOrient,iConv, pMap, pProb, param, *RefMap);
+ compareRefMapShifted<1>(iRefMap, iOrient, iConv, pMap, pProb, param, *RefMap);
}
}
@@ -56,7 +56,7 @@ __global__ void cudaZeroMem(void* ptr, size_t size)
int mysize = size / sizeof(int);
int myid = myBlockDimX * myBlockIdxX + myThreadIdxX;
int mygrid = myBlockDimX * myGridDimX;
- for (int i = myid;i < mysize;i += mygrid) myptr[i] = 0;
+ for (int i = myid; i < mysize; i += mygrid) myptr[i] = 0;
}
__global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iConv, const myfloat_t* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap RefMap, const int blockoffset, const int nShifts, const int nShiftBits, const int maxRef)
@@ -72,7 +72,7 @@ __global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iCon
const bool threadActive = myShiftIdx < nShifts && myShiftIdy < nShifts && iRefMap < maxRef;
- compareRefMap<2>(iRefMap,iOrient,iConv, pMap, pProb, param, RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef, threadActive);
+ compareRefMap<2>(iRefMap, iOrient, iConv, pMap, pProb, param, RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef, threadActive);
}
__global__ void multComplexMap(const mycomplex_t* convmap, const mycomplex_t* refmap, mycuComplex_t* out, const int NumberPixelsTotal, const int MapSize, const int NumberMaps, const int Offset)
@@ -80,7 +80,7 @@ __global__ void multComplexMap(const mycomplex_t* convmap, const mycomplex_t* re
if (myBlockIdxX >= NumberMaps) return;
const mycomplex_t* myin = &refmap[myBlockIdxX * MapSize + Offset];
mycuComplex_t* myout = &out[(myBlockIdxX * MapSize)];
- for(int i = myThreadIdxX;i < NumberPixelsTotal;i += myBlockDimX)
+ for(int i = myThreadIdxX; i < NumberPixelsTotal; i += myBlockDimX)
{
myout[i].x = convmap[i][0] * myin[i][0] + convmap[i][1] * myin[i][1];
myout[i].y = convmap[i][1] * myin[i][0] - convmap[i][0] * myin[i][1];
@@ -123,7 +123,7 @@ int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* c
if (FFTAlgo)
{
checkCudaErrors(cudaMemcpyAsync(&pConvMapFFT[(iConv & 1) * param.FFTMapSize], localmultFFT, param.FFTMapSize * sizeof(mycomplex_t), cudaMemcpyHostToDevice, cudaStream));
- for (int i = 0;i < maxRef;i += CUDA_FFTS_AT_ONCE)
+ for (int i = 0; i < maxRef; i += CUDA_FFTS_AT_ONCE)
{
const int num = min(CUDA_FFTS_AT_ONCE, maxRef - i);
multComplexMap<<<num, CUDA_THREAD_COUNT, 0, cudaStream>>>(&pConvMapFFT[(iConv & 1) * param.FFTMapSize], pRefMapsFFT, pFFTtmp2, param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D, param.FFTMapSize, num, i);
@@ -161,24 +161,24 @@ int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* c
const int nShiftBits = ilog2(nShifts);
size_t totalBlocks = divup((size_t) maxRef * (size_t) nShifts * (size_t) nShifts, (size_t) CUDA_THREAD_COUNT);
size_t nBlocks = CUDA_BLOCK_COUNT;
- for (size_t i = 0;i < totalBlocks;i += nBlocks)
+ for (size_t i = 0; i < totalBlocks; i += nBlocks)
{
- compareRefMapLoopShifts_kernel <<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, *gpumap, i, nShifts, nShiftBits, maxRef);
+ compareRefMapLoopShifts_kernel<<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream >>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, *gpumap, i, nShifts, nShiftBits, maxRef);
}
}
else if (GPUAlgo == 1) //Split shifts in multiple kernels
{
- for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
+ for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x = cent_x + param.param_device.GridSpaceCenter)
{
- for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
+ for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y = cent_y + param.param_device.GridSpaceCenter)
{
- compareRefMap_kernel <<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y, maxRef);
+ compareRefMap_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y, maxRef);
}
}
}
else if (GPUAlgo == 0) //All shifts in one kernel
{
- compareRefMapShifted_kernel <<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, maxRef);
+ compareRefMapShifted_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, maxRef);
}
else
{
@@ -224,7 +224,7 @@ int bioem_cuda::deviceInit()
checkCudaErrors(cudaStreamCreate(&cudaStream));
checkCudaErrors(cudaMalloc(&pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap));
- for (int i = 0;i < 2;i++)
+ for (int i = 0; i < 2; i++)
{
checkCudaErrors(cudaEventCreate(&cudaEvent[i]));
checkCudaErrors(cudaMalloc(&pConvMap_device[i], sizeof(myfloat_t) * RefMap.refMapSize));
@@ -259,7 +259,7 @@ int bioem_cuda::deviceExit()
cudaFree(pProb_device);
cudaFree(sum);
cudaFree(sumsquare);
- for (int i = 0;i < 2;i++)
+ for (int i = 0; i < 2; i++)
{
cudaEventDestroy(cudaEvent[i]);
cudaFree(pConvMap_device[i]);
@@ -294,7 +294,7 @@ int bioem_cuda::deviceStartRun()
if (FFTAlgo)
{
- for (int i = 0;i < 2;i++)
+ for (int i = 0; i < 2; i++)
{
int n[2] = {param.param_device.NumberPixels, param.param_device.NumberPixels};
if (cufftPlanMany(&plan[i], 2, n, NULL, 1, 0, NULL, 1, 0, CUFFT_C2R, i ? (maxRef % CUDA_FFTS_AT_ONCE) : CUDA_FFTS_AT_ONCE) != CUFFT_SUCCESS)
@@ -324,7 +324,7 @@ int bioem_cuda::deviceFinishRun()
if (FFTAlgo)
{
- for (int i = 0;i < 2;i++) cufftDestroy(plan[i]);
+ for (int i = 0; i < 2; i++) cufftDestroy(plan[i]);
}
return(0);
diff --git a/include/bioem.h b/include/bioem.h
index 4295794..4fda36d 100644
--- a/include/bioem.h
+++ b/include/bioem.h
@@ -21,13 +21,13 @@ public:
int dopreCalCrossCorrelation(int iRefMap, int iRefMapLocal);
int run();
int doProjections(int iMap);
- int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj, myfloat_t* Mapconv,mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC);
+ int createConvolutedProjectionMap(int iOreint, int iMap, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC);
virtual int compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
int createProjection(int iMap, mycomplex_t* map);
int calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare);
- void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC);
+ void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC);
bioem_Probability* pProb;
diff --git a/include/map.h b/include/map.h
index f97a650..7c30479 100644
--- a/include/map.h
+++ b/include/map.h
@@ -63,12 +63,12 @@ public:
ntotRefMap = map.ntotRefMap;
memcpy(sum_RefMap, map.sum_RefMap, sizeof(sum_RefMap));
memcpy(sumsquare_RefMap, map.sumsquare_RefMap, sizeof(sumsquare_RefMap));
-#pragma omp parallel for
- for (int i = 0;i < ntotRefMap;i++)
+ #pragma omp parallel for
+ for (int i = 0; i < ntotRefMap; i++)
{
- for (int j = 0;j < BIOEM_MAP_SIZE_X;j++)
+ for (int j = 0; j < BIOEM_MAP_SIZE_X; j++)
{
- for (int k = 0;k < BIOEM_MAP_SIZE_Y;k++)
+ for (int k = 0; k < BIOEM_MAP_SIZE_Y; k++)
{
Ref[j][k][i] = map.get(i, j, k);
}
diff --git a/include/model.h b/include/model.h
index bc08678..f325e14 100644
--- a/include/model.h
+++ b/include/model.h
@@ -15,7 +15,7 @@ public:
myfloat_t getAminoAcidRad(char *name);
myfloat_t getAminoAcidDensity(char *name);
- myfloat_t NormDen;
+ myfloat_t NormDen;
const char* filemodel;
diff --git a/main.cpp b/main.cpp
index edcb7cb..c9b6ad3 100644
--- a/main.cpp
+++ b/main.cpp
@@ -25,7 +25,7 @@ int main(int argc, char* argv[])
/********************************* Main BioEM code **********************************/
/************************************************************************************/
-#pragma omp parallel
+ #pragma omp parallel
{
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
}
@@ -45,7 +45,7 @@ int main(int argc, char* argv[])
/************ Configuration and Pre-calculating necessary objects *****************/
printf("Configuring\n");
- if (bio->configure(argc,argv)) return(1);
+ if (bio->configure(argc, argv)) return(1);
/******************************* Run BioEM routine ******************************/
printf("Running\n");
diff --git a/map.cpp b/map.cpp
index 3ef06d2..a81bc62 100644
--- a/map.cpp
+++ b/map.cpp
@@ -37,9 +37,9 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
}
else
{
- int nummap=-1;
- int lasti=0;
- int lastj=0;
+ int nummap = -1;
+ int lasti = 0;
+ int lastj = 0;
ifstream input(param.filemap);
if (!input.good())
{
@@ -52,12 +52,12 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
while (!input.eof())
{
- input.getline(line,512);
+ input.getline(line, 512);
- strncpy(tmpLine,line,strlen(line));
- char *token = strtok(tmpLine," ");
+ strncpy(tmpLine, line, strlen(line));
+ char *token = strtok(tmpLine, " ");
- if (strcmp(token,"PARTICLE")==0) // to count the number of maps
+ if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
{
nummap++;
if (allocsize == 0)
@@ -80,7 +80,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
cout << "BIOEM_MAX_MAPS too small\n";
exit(1);
}
- if(lasti!=param.param_device.NumberPixels && lastj!=param.param_device.NumberPixels && nummap>0)
+ if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
exit(1);
@@ -88,25 +88,25 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
}
else
{
- int i,j;
+ int i, j;
float z;
char tmpVals[36] = {' '};
- strncpy (tmpVals,line,8);
- sscanf (tmpVals,"%d",&i);
+ strncpy (tmpVals, line, 8);
+ sscanf (tmpVals, "%d", &i);
- strncpy (tmpVals,line+8,8);
- sscanf (tmpVals,"%d",&j);
+ strncpy (tmpVals, line + 8, 8);
+ sscanf (tmpVals, "%d", &j);
- strncpy (tmpVals,line+16,12);
- sscanf (tmpVals,"%f",&z);
+ strncpy (tmpVals, line + 16, 12);
+ sscanf (tmpVals, "%f", &z);
//checking for Map limits
- if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
+ if(i > 0 && i - 1 < BIOEM_MAP_SIZE_X && j > 0 && j - 1 < BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
{
- maps[nummap * refMapSize + (i-1) * numPixels + (j-1)] = (myfloat_t) z;
- lasti=i;
- lastj=j;
+ maps[nummap * refMapSize + (i - 1) * numPixels + (j - 1)] = (myfloat_t) z;
+ lasti = i;
+ lastj = j;
}
else
{
@@ -116,7 +116,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
}
}
cout << ".";
- ntotRefMap=nummap+1;
+ ntotRefMap = nummap + 1;
maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
@@ -149,33 +149,33 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
myfloat_t* localMap;
- localMap= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
+ localMap = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
RefMapsFFT = new mycomplex_t[ntotRefMap * param.FFTMapSize];
mycomplex_t* localout;
- localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
+ localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
{
//Assigning localMap values to padded 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++)
{
- localMap[i * param.param_device.NumberPixels+j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
+ localMap[i * param.param_device.NumberPixels + j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
}
}
//Calling FFT_Forward
- myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localMap,localout);
+ myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
// Normalizing and saving the Reference CTFs
mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.FFTMapSize];
- 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++ )
{
- RefMap[i][0]=localout[i][0];
- RefMap[i][1]=localout[i][1];
+ RefMap[i][0] = localout[i][0];
+ RefMap[i][1] = localout[i][1];
}
}
diff --git a/model.cpp b/model.cpp
index a2151f9..0eaa1fd 100644
--- a/model.cpp
+++ b/model.cpp
@@ -30,19 +30,19 @@ int bioem_model::readModel()
char line[512] = {' '};
char tmpLine[512] = {' '};
- int numres=0;
- NormDen=0.0;
+ int numres = 0;
+ NormDen = 0.0;
// cout << " HERE " << filemodel ;
// for eachline in the file
while (!input.eof())
{
- input.getline(line,512);
+ input.getline(line, 512);
- strncpy(tmpLine,line,strlen(line));
- char *token = strtok(tmpLine," ");
+ strncpy(tmpLine, line, strlen(line));
+ char *token = strtok(tmpLine, " ");
- if (strcmp(token,"ATOM")==0) // Optional,Mandatory if standard residues exist
+ if (strcmp(token, "ATOM") == 0) // Optional,Mandatory if standard residues exist
{
/*
1-6 "ATOM "
@@ -66,26 +66,26 @@ int bioem_model::readModel()
char tmp[6] = {' '};
// parse name
- strncpy(tmp,line+12,4);
- sscanf (tmp,"%s",name);
+ strncpy(tmp, line + 12, 4);
+ sscanf (tmp, "%s", name);
// parse resName
- strncpy(tmp,line+17,3);
- sscanf (tmp,"%s",resName);
+ strncpy(tmp, line + 17, 3);
+ sscanf (tmp, "%s", resName);
// parse x, y, z
char tmpVals[36] = {' '};
- strncpy (tmpVals,line+30,8);
- sscanf (tmpVals,"%f",&x);
+ strncpy (tmpVals, line + 30, 8);
+ sscanf (tmpVals, "%f", &x);
- strncpy (tmpVals,line+38,8);
- sscanf (tmpVals,"%f",&y);
+ strncpy (tmpVals, line + 38, 8);
+ sscanf (tmpVals, "%f", &y);
- strncpy (tmpVals,line+46,8);
- sscanf (tmpVals,"%f",&z);
+ strncpy (tmpVals, line + 46, 8);
+ sscanf (tmpVals, "%f", &z);
- if (strcmp(name,"CA") == 0)
+ if (strcmp(name, "CA") == 0)
{
if (numres >= BIOEM_MODEL_SIZE)
{
@@ -93,30 +93,30 @@ int bioem_model::readModel()
exit(1);
}
//Getting residue Radius and electron density
- radiusPointsModel[numres]=getAminoAcidRad(resName);
- densityPointsModel[numres]=getAminoAcidDensity(resName);
- NormDen+=densityPointsModel[numres];
+ radiusPointsModel[numres] = getAminoAcidRad(resName);
+ densityPointsModel[numres] = getAminoAcidDensity(resName);
+ NormDen += densityPointsModel[numres];
//Getting the coordinates
PointsModel[numres].pos[0] = (myfloat_t) x;
PointsModel[numres].pos[1] = (myfloat_t) y;
PointsModel[numres].pos[2] = (myfloat_t) z;
- exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres] << "\n";
+ exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
numres++;
}
}
}
- nPointsModel=numres;
+ nPointsModel = numres;
cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
}
else //Reading model from FILE FORMAT x,y,z,rad,density
{
char line[128];
- int numres=0;
- NormDen=0.0;
+ int numres = 0;
+ NormDen = 0.0;
FILE *file = fopen ( filemodel , "r" );
if ( file != NULL )
{
@@ -135,11 +135,11 @@ int bioem_model::readModel()
radiusPointsModel[numres] = (myfloat_t) tmpval[3];
densityPointsModel[numres] = (myfloat_t) tmpval[4];
- exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]<< "\n";
- NormDen+=densityPointsModel[numres];
+ exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " << PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
+ NormDen += densityPointsModel[numres];
numres++;
}
- nPointsModel=numres;
+ nPointsModel = numres;
cout << "Protein structure read from Standard File \nTotal Number of Residues " << nPointsModel ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
}
@@ -149,22 +149,22 @@ int bioem_model::readModel()
//Moving to Model to its center of mass:
myfloat3_t r_cm;
- for(int n=0; n < 3; n++)r_cm.pos[n] = 0.0;
+ for(int n = 0; n < 3; n++)r_cm.pos[n] = 0.0;
- for(int n=0; n < nPointsModel; n++)
+ for(int n = 0; n < nPointsModel; n++)
{
r_cm.pos[0] += PointsModel[n].pos[0];
r_cm.pos[1] += PointsModel[n].pos[1];
r_cm.pos[2] += PointsModel[n].pos[2];
}
- r_cm.pos[0]=r_cm.pos[0] / (myfloat_t) nPointsModel;
- r_cm.pos[1]=r_cm.pos[1] / (myfloat_t) nPointsModel;
- r_cm.pos[2]=r_cm.pos[2] / (myfloat_t) nPointsModel;
- for(int n=0; n < nPointsModel; n++)
+ r_cm.pos[0] = r_cm.pos[0] / (myfloat_t) nPointsModel;
+ r_cm.pos[1] = r_cm.pos[1] / (myfloat_t) nPointsModel;
+ r_cm.pos[2] = r_cm.pos[2] / (myfloat_t) nPointsModel;
+ for(int n = 0; n < nPointsModel; n++)
{
- PointsModel[n].pos[0]-= r_cm.pos[0];
- PointsModel[n].pos[1]-= r_cm.pos[1];
- PointsModel[n].pos[2]-= r_cm.pos[2];
+ PointsModel[n].pos[0] -= r_cm.pos[0];
+ PointsModel[n].pos[1] -= r_cm.pos[1];
+ PointsModel[n].pos[2] -= r_cm.pos[2];
}
@@ -174,28 +174,28 @@ int bioem_model::readModel()
myfloat_t bioem_model::getAminoAcidRad(char *name)
{
/*************** Function that gets the radius for each amino acid ****************/
- myfloat_t iaa=0;
-
- if(std::strcmp(name,"CYS")==0)iaa=2.75;
- if(std::strcmp(name,"PHE")==0)iaa=3.2;
- if(std::strcmp(name,"LEU")==0)iaa=3.1;
- if(std::strcmp(name,"TRP")==0)iaa=3.4;
- if(std::strcmp(name,"VAL")==0)iaa=2.95;
- if(std::strcmp(name,"ILE")==0)iaa=3.1;
- if(std::strcmp(name,"MET")==0)iaa=3.1;
- if(std::strcmp(name,"HIS")==0)iaa=3.05;
- if(std::strcmp(name,"TYR")==0)iaa=3.25;
- if(std::strcmp(name,"ALA")==0)iaa=2.5;
- if(std::strcmp(name,"GLY")==0)iaa=2.25;
- if(std::strcmp(name,"PRO")==0)iaa=2.8;
- if(std::strcmp(name,"ASN")==0)iaa=2.85;
- if(std::strcmp(name,"THR")==0)iaa=2.8;
- if(std::strcmp(name,"SER")==0)iaa=2.6;
- if(std::strcmp(name,"ARG")==0)iaa=3.3;
- if(std::strcmp(name,"GLN")==0)iaa=3.0;
- if(std::strcmp(name,"ASP")==0)iaa=2.8;
- if(std::strcmp(name,"LYS")==0)iaa=3.2;
- if(std::strcmp(name,"GLU")==0)iaa=2.95;
+ myfloat_t iaa = 0;
+
+ if(std::strcmp(name, "CYS") == 0)iaa = 2.75;
+ if(std::strcmp(name, "PHE") == 0)iaa = 3.2;
+ if(std::strcmp(name, "LEU") == 0)iaa = 3.1;
+ if(std::strcmp(name, "TRP") == 0)iaa = 3.4;
+ if(std::strcmp(name, "VAL") == 0)iaa = 2.95;
+ if(std::strcmp(name, "ILE") == 0)iaa = 3.1;
+ if(std::strcmp(name, "MET") == 0)iaa = 3.1;
+ if(std::strcmp(name, "HIS") == 0)iaa = 3.05;
+ if(std::strcmp(name, "TYR") == 0)iaa = 3.25;
+ if(std::strcmp(name, "ALA") == 0)iaa = 2.5;
+ if(std::strcmp(name, "GLY") == 0)iaa = 2.25;
+ if(std::strcmp(name, "PRO") == 0)iaa = 2.8;
+ if(std::strcmp(name, "ASN") == 0)iaa = 2.85;
+ if(std::strcmp(name, "THR") == 0)iaa = 2.8;
+ if(std::strcmp(name, "SER") == 0)iaa = 2.6;
+ if(std::strcmp(name, "ARG") == 0)iaa = 3.3;
+ if(std::strcmp(name, "GLN") == 0)iaa = 3.0;
+ if(std::strcmp(name, "ASP") == 0)iaa = 2.8;
+ if(std::strcmp(name, "LYS") == 0)iaa = 3.2;
+ if(std::strcmp(name, "GLU") == 0)iaa = 2.95;
if(iaa == 0)
{
@@ -209,28 +209,28 @@ myfloat_t bioem_model::getAminoAcidRad(char *name)
myfloat_t bioem_model::getAminoAcidDensity(char *name)
{
/*************** Function that gets the number of electrons for each amino acid ****************/
- myfloat_t iaa=0.0;
-
- if(std::strcmp(name,"CYS")==0)iaa=64.0;
- if(std::strcmp(name,"PHE")==0)iaa=88.0;
- if(std::strcmp(name,"LEU")==0)iaa=72.0;
- if(std::strcmp(name,"TRP")==0)iaa=108.0;
- if(std::strcmp(name,"VAL")==0)iaa=64.0;
- if(std::strcmp(name,"ILE")==0)iaa=72.0;
- if(std::strcmp(name,"MET")==0)iaa=80.0;
- if(std::strcmp(name,"HIS")==0)iaa=82.0;
- if(std::strcmp(name,"TYR")==0)iaa=96.0;
- if(std::strcmp(name,"ALA")==0)iaa=48.0;
- if(std::strcmp(name,"GLY")==0)iaa=40.0;
- if(std::strcmp(name,"PRO")==0)iaa=62.0;
- if(std::strcmp(name,"ASN")==0)iaa=66.0;
- if(std::strcmp(name,"THR")==0)iaa=64.0;
- if(std::strcmp(name,"SER")==0)iaa=56.0;
- if(std::strcmp(name,"ARG")==0)iaa=93.0;
- if(std::strcmp(name,"GLN")==0)iaa=78.0;
- if(std::strcmp(name,"ASP")==0)iaa=59.0;
- if(std::strcmp(name,"LYS")==0)iaa=79.0;
- if(std::strcmp(name,"GLU")==0)iaa=53.0;
+ myfloat_t iaa = 0.0;
+
+ if(std::strcmp(name, "CYS") == 0)iaa = 64.0;
+ if(std::strcmp(name, "PHE") == 0)iaa = 88.0;
+ if(std::strcmp(name, "LEU") == 0)iaa = 72.0;
+ if(std::strcmp(name, "TRP") == 0)iaa = 108.0;
+ if(std::strcmp(name, "VAL") == 0)iaa = 64.0;
+ if(std::strcmp(name, "ILE") == 0)iaa = 72.0;
+ if(std::strcmp(name, "MET") == 0)iaa = 80.0;
+ if(std::strcmp(name, "HIS") == 0)iaa = 82.0;
+ if(std::strcmp(name, "TYR") == 0)iaa = 96.0;
+ if(std::strcmp(name, "ALA") == 0)iaa = 48.0;
+ if(std::strcmp(name, "GLY") == 0)iaa = 40.0;
+ if(std::strcmp(name, "PRO") == 0)iaa = 62.0;
+ if(std::strcmp(name, "ASN") == 0)iaa = 66.0;
+ if(std::strcmp(name, "THR") == 0)iaa = 64.0;
+ if(std::strcmp(name, "SER") == 0)iaa = 56.0;
+ if(std::strcmp(name, "ARG") == 0)iaa = 93.0;
+ if(std::strcmp(name, "GLN") == 0)iaa = 78.0;
+ if(std::strcmp(name, "ASP") == 0)iaa = 59.0;
+ if(std::strcmp(name, "LYS") == 0)iaa = 79.0;
+ if(std::strcmp(name, "GLU") == 0)iaa = 53.0;
if(iaa == 0.0)
{
diff --git a/param.cpp b/param.cpp
index 64e2178..342ea4b 100644
--- a/param.cpp
+++ b/param.cpp
@@ -14,7 +14,7 @@ using namespace std;
bioem_param::bioem_param()
{
//Number of Pixels
- param_device.NumberPixels=0;
+ param_device.NumberPixels = 0;
param_device.NumberFFTPixels1D = 0;
// Euler angle grid spacing
angleGridPointsAlpha = 0;
@@ -57,122 +57,122 @@ int bioem_param::readParameters()
cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
while (!input.eof())
{
- input.getline(line,512);
- strcpy(saveline,line);
- char *token = strtok(line," ");
+ input.getline(line, 512);
+ strcpy(saveline, line);
+ char *token = strtok(line, " ");
- if (token==NULL || line[0] == '#' || strlen(token)==0)
+ if (token == NULL || line[0] == '#' || strlen(token) == 0)
{
// comment or blank line
}
- else if (strcmp(token,"PIXEL_SIZE")==0)
+ else if (strcmp(token, "PIXEL_SIZE") == 0)
{
- token = strtok(NULL," ");
- pixelSize=atof(token);
+ token = strtok(NULL, " ");
+ pixelSize = atof(token);
cout << "Pixel Sixe " << pixelSize << "\n";
}
- else if (strcmp(token,"NUMBER_PIXELS")==0)
+ else if (strcmp(token, "NUMBER_PIXELS") == 0)
{
- token = strtok(NULL," ");
- param_device.NumberPixels=int(atoi(token));
+ token = strtok(NULL, " ");
+ param_device.NumberPixels = int(atoi(token));
cout << "Number of Pixels " << param_device.NumberPixels << "\n";
}
- else if (strcmp(token,"GRIDPOINTS_ALPHA")==0)
+ else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
{
- token = strtok(NULL," ");
- angleGridPointsAlpha=int(atoi(token));
+ token = strtok(NULL, " ");
+ angleGridPointsAlpha = int(atoi(token));
cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
}
- else if (strcmp(token,"GRIDPOINTS_BETA")==0)
+ else if (strcmp(token, "GRIDPOINTS_BETA") == 0)
{
- token = strtok(NULL," ");
- angleGridPointsBeta=int(atoi(token));
+ token = strtok(NULL, " ");
+ angleGridPointsBeta = int(atoi(token));
cout << "Grid points beta " << angleGridPointsBeta << "\n";
}
- else if (strcmp(token,"GRIDPOINTS_ENVELOPE")==0)
+ else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
{
- token = strtok(NULL," ");
- numberGridPointsEnvelop=int(atoi(token));
+ token = strtok(NULL, " ");
+ numberGridPointsEnvelop = int(atoi(token));
cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
}
- else if (strcmp(token,"START_ENVELOPE")==0)
+ else if (strcmp(token, "START_ENVELOPE") == 0)
{
- token = strtok(NULL," ");
- startGridEnvelop=atof(token);
+ token = strtok(NULL, " ");
+ startGridEnvelop = atof(token);
cout << "Start Envelope " << startGridEnvelop << "\n";
}
- else if (strcmp(token,"GRIDSPACE_ENVELOPE")==0)
+ else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0)
{
- token = strtok(NULL," ");
- gridEnvelop=atof(token);
+ token = strtok(NULL, " ");
+ gridEnvelop = atof(token);
cout << "Grid spacing Envelope " << gridEnvelop << "\n";
}
- else if (strcmp(token,"GRIDPOINTS_CTF_PHASE")==0)
+ else if (strcmp(token, "GRIDPOINTS_CTF_PHASE") == 0)
{
- token = strtok(NULL," ");
- numberGridPointsCTF_phase=int(atoi(token));
+ token = strtok(NULL, " ");
+ numberGridPointsCTF_phase = int(atoi(token));
cout << "Grid points CTF " << numberGridPointsCTF_phase << "\n";
}
- else if (strcmp(token,"START_CTF_PHASE")==0)
+ else if (strcmp(token, "START_CTF_PHASE") == 0)
{
- token = strtok(NULL," ");
- startGridCTF_phase=atof(token);
+ token = strtok(NULL, " ");
+ startGridCTF_phase = atof(token);
cout << "Start CTF " << startGridCTF_phase << "\n";
}
- else if (strcmp(token,"GRIDSPACE_CTF_PHASE")==0)
+ else if (strcmp(token, "GRIDSPACE_CTF_PHASE") == 0)
{
- token = strtok(NULL," ");
- gridCTF_phase=atof(token);
+ token = strtok(NULL, " ");
+ gridCTF_phase = atof(token);
cout << "Grid Space CTF " << gridCTF_phase << "\n";
- } else if (strcmp(token,"GRIDPOINTS_CTF_AMP")==0)
+ } else if (strcmp(token, "GRIDPOINTS_CTF_AMP") == 0)
{
- token = strtok(NULL," ");
- numberGridPointsCTF_amp=int(atoi(token));
+ token = strtok(NULL, " ");
+ numberGridPointsCTF_amp = int(atoi(token));
cout << "Grid points CTF " << numberGridPointsCTF_amp << "\n";
}
- else if (strcmp(token,"START_CTF_AMP")==0)
+ else if (strcmp(token, "START_CTF_AMP") == 0)
{
- token = strtok(NULL," ");
- startGridCTF_amp=atof(token);
+ token = strtok(NULL, " ");
+ startGridCTF_amp = atof(token);
cout << "Start CTF " << startGridCTF_amp << "\n";
}
- else if (strcmp(token,"GRIDSPACE_CTF_AMP")==0)
+ else if (strcmp(token, "GRIDSPACE_CTF_AMP") == 0)
{
- token = strtok(NULL," ");
- gridCTF_amp=atof(token);
+ token = strtok(NULL, " ");
+ gridCTF_amp = atof(token);
cout << "Grid Space CTF " << gridCTF_amp << "\n";
}
- else if (strcmp(token,"MAX_D_CENTER")==0)
+ else if (strcmp(token, "MAX_D_CENTER") == 0)
{
- token = strtok(NULL," ");
- param_device.maxDisplaceCenter=int(atoi(token));
+ token = strtok(NULL, " ");
+ param_device.maxDisplaceCenter = int(atoi(token));
cout << "Maximum displacement Center " << param_device.maxDisplaceCenter << "\n";
}
- else if (strcmp(token,"PIXEL_GRID_CENTER")==0)
+ else if (strcmp(token, "PIXEL_GRID_CENTER") == 0)
{
- token = strtok(NULL," ");
- param_device.GridSpaceCenter=int(atoi(token));
+ token = strtok(NULL, " ");
+ param_device.GridSpaceCenter = int(atoi(token));
cout << "Grid space displacement center " << param_device.GridSpaceCenter << "\n";
}
- else if (strcmp(token,"WRITE_PROB_ANGLES")==0)
+ else if (strcmp(token, "WRITE_PROB_ANGLES") == 0)
{
- writeAngles=true;
+ writeAngles = true;
cout << "Writing Probabilies of each angle \n";
}
@@ -223,14 +223,14 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
/**************************************************************************************/
/**************** Routine that pre-calculates Euler angle grids **********************/
/************************************************************************************/
- myfloat_t grid_alpha,cos_grid_beta;
- int n=0;
+ myfloat_t grid_alpha, cos_grid_beta;
+ int n = 0;
//alpha and gamma are uniform in -PI,PI
- grid_alpha=2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
+ grid_alpha = 2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
//cosine beta is uniform in -1,1
- cos_grid_beta=2.f / (myfloat_t) angleGridPointsBeta;
+ cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
// Euler Angle Array
for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
@@ -239,25 +239,25 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
{
for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
{
- angles[n].pos[0]= (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
- angles[n].pos[1]= acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
- angles[n].pos[2]= (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
+ angles[n].pos[0] = (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
+ angles[n].pos[1] = acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
+ angles[n].pos[2] = (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
n++;
}
}
}
- nTotGridAngles=n;
+ nTotGridAngles = n;
/********** Calculating normalized volumen element *********/
param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
- * gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
- / ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
+ * gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
+ / ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
/*** Number of total pixels***/
- param_device.Ntotpi= (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
- param_device.NtotDist=(2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1);
+ param_device.Ntotpi = (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
+ param_device.NtotDist = (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1);
return(0);
@@ -269,12 +269,12 @@ int bioem_param::CalculateRefCTF()
/********** Routine that pre-calculates Kernels for Convolution **********************/
/************************************************************************************/
- myfloat_t amp,env,phase,ctf,radsq;
+ myfloat_t amp, env, phase, ctf, radsq;
myfloat_t* localCTF;
- int nctfmax= param_device.NumberPixels / 2;
- int n=0;
+ int nctfmax = param_device.NumberPixels / 2;
+ int n = 0;
- localCTF= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels*param_device.NumberPixels);
+ localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
delete[] refCTF;
@@ -292,40 +292,40 @@ int bioem_param::CalculateRefCTF()
for (int ienv = 0; ienv < numberGridPointsEnvelop ; ienv++)//Loop over envelope
{
- env= (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
+ env = (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
- memset(localCTF,0,param_device.NumberPixels*param_device.NumberPixels*sizeof(myfloat_t));
+ memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
//Assigning CTF values
- for(int i=0; i< nctfmax; i++)
+ for(int i = 0; i < nctfmax; i++)
{
- for(int j=0; j< nctfmax; j++)
+ for(int j = 0; j < nctfmax; j++)
{
- radsq=(myfloat_t) (i*i+j*j) * pixelSize * pixelSize;
- ctf=exp(-radsq*env/2.0)*(amp*cos(radsq*phase/2.0)-sqrt((1-amp*amp))*sin(radsq*phase/2.0));
+ radsq = (myfloat_t) (i * i + j * j) * pixelSize * pixelSize;
+ ctf = exp(-radsq * env / 2.0) * (amp * cos(radsq * phase / 2.0) - sqrt((1 - amp * amp)) * sin(radsq * phase / 2.0));
- localCTF[i*param_device.NumberPixels+j]=(myfloat_t) ctf;
- localCTF[i*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
- localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+j]=(myfloat_t) ctf;
- localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
+ localCTF[i * param_device.NumberPixels + j] = (myfloat_t) ctf;
+ localCTF[i * param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
+ localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + j] = (myfloat_t) ctf;
+ localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
}
}
//Creatting FFT_Forward of Kernel to store
mycomplex_t* localout;
- localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param_device.NumberPixels*param_device.NumberFFTPixels1D);
+ localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberFFTPixels1D);
//Calling FFT_Forward
- myfftw_execute_dft_r2c(fft_plan_r2c_forward,localCTF,localout);
+ myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
// Normalizing and saving the Reference CTFs
mycomplex_t* curRef = &refCTF[n * FFTMapSize];
- for(int i=0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
+ for(int i = 0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
{
curRef[i][0] = localout[i][0];
curRef[i][1] = localout[i][1];
}
- CtfParam[n].pos[0]=amp;
- CtfParam[n].pos[1]=phase;
- CtfParam[n].pos[2]=env;
+ CtfParam[n].pos[0] = amp;
+ CtfParam[n].pos[1] = phase;
+ CtfParam[n].pos[2] = env;
n++;
myfftw_free(localout);
}
@@ -346,7 +346,7 @@ int bioem_param::CalculateRefCTF()
bioem_param::~bioem_param()
{
releaseFFTPlans();
- param_device.NumberPixels=0;
+ param_device.NumberPixels = 0;
angleGridPointsAlpha = 0;
angleGridPointsBeta = 0;
numberGridPointsEnvelop = 0;
--
GitLab