Skip to content

GitLab

B
BioEM

There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

Commit 65ca8537 authored by David Rohr's avatar David Rohr
Browse files
Options

convert indent spaces to tabs, unify indentation

parent a30096fc
Hide whitespace changes
Inline Side-by-side
Showing with 324 additions and 324 deletions
bioem.cpp
View file @ 65ca8537
......@@ -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 = &param.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);
......
bioem_algorithm.h
View file @ 65ca8537
......@@ -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);
}
......
bioem_cuda.cu
View file @ 65ca8537
......@@ -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];
......