Merge FFT Algorithm into BioEM version in git repository, cleanup, ensure...

Merge FFT Algorithm into BioEM version in git repository, cleanup, ensure consistent indentation style

bioem.cpp

@@ -37,7 +37,7 @@ ostream& operator<<(ostream& os, const vector<T>& v)
 
 bioem::bioem()
 {
-
+	FFTAlgo = getenv("FFTALGO") == NULL ? 0 : atoi(getenv("FFTALGO"));
 }
 
 bioem::~bioem()
@@ -159,7 +159,7 @@ int bioem::configure(int ac, char* av[])
 	/********* HERE: PROBLEM if maps dont fit on the memory!! ***************/
 	// copying mapfile to ref map class
 	RefMap.filemap = mapfile.c_str();
-    RefMap.readRefMaps();
+	RefMap.readRefMaps(param);
 
 	/****************** Precalculating Necessary Stuff *********************/
 	precalculate();
@@ -200,6 +200,9 @@ int bioem::precalculate()
 	// Precalculating CTF Kernels stored in class Param
 	param.CalculateRefCTF();
 
+	// Precalculating Maps in Fourier space
+	RefMap.PreCalculateMapsFFT(param);
+
 	return(0);
 }
 
@@ -213,6 +216,7 @@ int bioem::run()
 	/**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/
 	/****************** Declarying class of Probability Pointer  *************************/
 	pProb = new bioem_Probability[RefMap.ntotRefMap];
+	crossCor = new bioem_crossCor[RefMap.ntotRefMap];
 
 	printf("\tInitializing\n");
 	// Inizialzing Probabilites to zero and constant to -Infinity
@@ -235,9 +239,13 @@ int bioem::run()
 	/*** Declaring Private variables for each thread *****/
 	mycomplex_t* proj_mapFFT;
 	bioem_map conv_map;
+	mycomplex_t* conv_mapFFT;
+	myfloat_t sumCONV,sumsquareCONV;
 
 	//allocating fftw_complex vector
-    proj_mapFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
+	proj_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
+	conv_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t)*param.param_device.NumberPixels*param.param_device.NumberPixels);
+
 
 	HighResTimer timer;
 
@@ -258,14 +266,24 @@ int bioem::run()
 		{
 			printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
 			/*** Calculating convolutions of projection map and crosscorrelations ***/
+			memset(conv_mapFFT,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*conv_mapFFT));
+
 			timer.ResetStart();
-            createConvolutedProjectionMap(iProjectionOut,iConv,proj_mapFFT,conv_map);
+			createConvolutedProjectionMap(iProjectionOut,iConv,proj_mapFFT,conv_map,conv_mapFFT,sumCONV,sumsquareCONV);
 			printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());
 
 			/***************************************************************************************/
 			/*** Comparing each calculated convoluted map with all experimental maps ***/
 			timer.ResetStart();
+			if (FFTAlgo == 0)
+			{
 				compareRefMaps(iProjectionOut, iConv, conv_map);
+			}
+			else
+			{
+				compareRefMaps2(iProjectionOut, iConv,conv_mapFFT,sumCONV,sumsquareCONV);
+			}
+
 			const double compTime = timer.GetCurrentElapsedTime();
 			const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
 			const double nFlops = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
@@ -279,7 +297,8 @@ int bioem::run()
 
 	}
 	//deallocating fftw_complex vector
-    fftw_free(proj_mapFFT);
+	myfftw_free(proj_mapFFT);
+	myfftw_free(conv_mapFFT);
 
 	deviceFinishRun();
 
@@ -343,12 +362,23 @@ int bioem::run()
 		param.refCTF =NULL;
 	}
 
+	if (crossCor)
+	{
+		delete[] crossCor;
+		crossCor = NULL;
+	}
+
+	if(RefMap.RefMapFFT)
+	{
+		delete[] RefMap.RefMapFFT;
+		RefMap.RefMapFFT = NULL;
+	}
 	return(0);
 }
 
 int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, const int startMap)
 {
-#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);
@@ -356,6 +386,170 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
 	return(0);
 }
 
+int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myfloat_t sumC,myfloat_t sumsquareC)
+{
+	#pragma omp parallel for
+	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
+	{
+
+		mycomplex_t* localCCT;
+		localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
+
+		mycomplex_t* lCC;
+		lCC= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
+
+		//setting crossCor value to zero for each projection
+		for(int n=0; n < param.param_device.NtotDist ; n++)
+		{
+			crossCor[iRefMap].value[n]=0.0;
+			crossCor[iRefMap].disx[n]=-99999;
+			crossCor[iRefMap].disy[n]=-99999;
+		}
+
+// Before::		compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
+
+		calculateCCFFT(iRefMap,iConv, localConvFFT, localCCT,lCC);
+
+		myfftw_free(localCCT);
+		myfftw_free(lCC);
+	}
+
+
+//Not in openMP loop SUM OVER PROBABILITIES
+	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
+	{
+		calProb(iRefMap,iOrient,iConv,sumC,sumsquareC);
+	}
+
+	return(0);
+}
+
+/////////////NEW ROUTINE ////////////////
+int bioem::calculateCCFFT(int iRefMap, int iConv, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC)
+{
+	myfftw_plan plan;
+
+	memset(lCC,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*lCC));
+	memset(localCCT,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localCCT));
+
+	for(int i=0; i < param.param_device.NumberPixels ; i++ )
+	{
+		for(int j=0; j < param.param_device.NumberPixels ; j++ )
+		{
+			localCCT[i*param.param_device.NumberPixels+j][0]=localConvFFT[i*param.param_device.NumberPixels+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]+localConvFFT[i*param.param_device.NumberPixels+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1];
+			localCCT[i*param.param_device.NumberPixels+j][1]=localConvFFT[i*param.param_device.NumberPixels+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]-localConvFFT[i*param.param_device.NumberPixels+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1];
+		}
+	}
+
+	#pragma omp critical
+	{
+		plan = myfftw_plan_dft_2d(param.param_device.NumberPixels,param.param_device.NumberPixels,localCCT,lCC,FFTW_BACKWARD,FFTW_ESTIMATE);
+		myfftw_execute(plan);
+	}
+
+// Storing CORRELATIONS FOR CORRESPONDING DISPLACEMENTS & Normalizing after Backward FFT
+	int n=0;
+	for (int cent_x = 0; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
+	{
+		for (int cent_y = 0; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
+		{
+			crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels);
+			crossCor[iRefMap].disx[n]=cent_x;
+			crossCor[iRefMap].disy[n]=cent_y;
+			n++;
+		}
+		for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y < param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
+		{
+			crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels);
+			crossCor[iRefMap].disx[n]=cent_x;
+			crossCor[iRefMap].disy[n]=param.param_device.NumberPixels-cent_y;
+			n++;
+		}
+	}
+	for (int cent_x = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_x < param.param_device.NumberPixels; cent_x=cent_x+param.param_device.GridSpaceCenter)
+	{
+		for (int cent_y = 0; cent_y < param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
+		{
+			crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels);
+			crossCor[iRefMap].disx[n]=param.param_device.NumberPixels-cent_x;
+			crossCor[iRefMap].disy[n]=cent_y;
+			n++;
+		}
+		for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y <= param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
+		{
+			crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels);
+			crossCor[iRefMap].disx[n]=param.param_device.NumberPixels-cent_x;
+			crossCor[iRefMap].disy[n]=param.param_device.NumberPixels-cent_y;
+			n++;
+		}
+	}
+//#pragma omp critical
+	{
+		myfftw_destroy_plan(plan);
+	}
+	/* Controll but slows down the parallelisim
+	  if(n> MAX_DISPLACE && n> param.param_device.NtotDist)
+	 {
+		cout << "Problem with displace center and Max allowed" ;
+		exit(1);
+		}*/
+//
+	return (0);
+}
+
+int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC)
+{
+
+	/********************************************************/
+	/*********** Calculates the BioEM probability ***********/
+	/********************************************************/
+
+	const myfloat_t ForLogProb = (sumsquareC * param.param_device.Ntotpi - sumC * sumC);
+
+// Loop again over displacements
+	for(int n=0; n < param.param_device.NtotDist ; n++)
+	{
+		// Products of different cross-correlations (first element in formula)
+		const myfloat_t firstele = param.param_device.Ntotpi * (RefMap.sumsquare_RefMap[iRefMap] * sumsquareC -   crossCor[iRefMap].value[n]*  crossCor[iRefMap].value[n]) +
+								   2 * RefMap.sum_RefMap[iRefMap] * sumC *   crossCor[iRefMap].value[n] - RefMap.sumsquare_RefMap[iRefMap] * sumC * sumC - RefMap.sum_RefMap[iRefMap] * RefMap.sum_RefMap[iRefMap] * sumsquareC;
+
+		//******* 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);
+		const myfloat_t logpro = (3 - param.param_device.Ntotpi) * 0.5 * log(firstele) + (param.param_device.Ntotpi * 0.5 - 2) * log((param.param_device.Ntotpi - 2) * ForLogProb);
+//   cout << n <<" " << firstele << " "<< logpro << "\n";
+		{
+			/*******  Summing total Probabilities *************/
+			/******* Need a constant because of numerical divergence*****/
+			if(pProb[iRefMap].Constoadd < logpro)
+			{
+				pProb[iRefMap].Total = pProb[iRefMap].Total * exp(-logpro + pProb[iRefMap].Constoadd);
+				pProb[iRefMap].Constoadd = logpro;
+			}
+			pProb[iRefMap].Total += exp(logpro - pProb[iRefMap].Constoadd);
+
+			//Summing probabilities for each orientation
+			if(pProb[iRefMap].ConstAngle[iOrient] < logpro)
+			{
+				pProb[iRefMap].forAngles[iOrient] = pProb[iRefMap].forAngles[iOrient] * exp(-logpro + pProb[iRefMap].ConstAngle[iOrient]);
+				pProb[iRefMap].ConstAngle[iOrient] = logpro;
+			}
+			pProb[iRefMap].forAngles[iOrient] += exp(logpro - pProb[iRefMap].ConstAngle[iOrient]);
+
+			/********** Getting parameters that maximize the probability ***********/
+			if(pProb[iRefMap].max_prob < logpro)
+			{
+				pProb[iRefMap].max_prob = logpro;
+				pProb[iRefMap].max_prob_cent_x = crossCor[iRefMap].disx[n];
+				pProb[iRefMap].max_prob_cent_y = crossCor[iRefMap].disy[n];
+				pProb[iRefMap].max_prob_orient = iOrient;
+				pProb[iRefMap].max_prob_conv = iConv;
+			}
+		}
+	}
+	return (0);
+}
+
+
 int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 {
 
@@ -367,12 +561,11 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 	myfloat3_t RotatedPointsModel[Model.nPointsModel];
 	myfloat_t rotmat[3][3];
 	myfloat_t alpha, gam,beta;
-    fftw_plan plan;
+	myfftw_plan plan;
 	mycomplex_t* localproj;
-    int totnumPixFFT=2*param.param_device.NumberPixels;
 
-    localproj= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
-    memset(localproj,0,4*param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
+	localproj= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_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];
@@ -416,10 +609,11 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 	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)+float(param.param_device.NumberPixels)/2.0+0.5);
-        j=floor(RotatedPointsModel[n].pos[1]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
+		i=floor(RotatedPointsModel[n].pos[0]/param.pixelSize+float(param.param_device.NumberPixels)/2.0+0.5);
+		j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+float(param.param_device.NumberPixels)/2.0+0.5);
+
+		localproj[i*param.param_device.NumberPixels+j][0]+=Model.densityPointsModel[n]/Model.NormDen;
 
-        localproj[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0]+=Model.densityPointsModel[n];
 
 	}
 
@@ -431,9 +625,9 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 		myexamplemap.open ("MAP_i10");
 		myexampleRot.open ("Rot_i10");
 		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
-        for(int k=0; k<2*param.param_device.NumberPixels; k++)
+		for(int k=0; k<param.param_device.NumberPixels; k++)
 		{
-            for(int j=0; j<2*param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*2*param.param_device.NumberPixels+j][0];
+			for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j][0];
 		}
 		myexamplemap << " \n";
 		for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
@@ -445,17 +639,17 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 	/********** Omp Critical is necessary with FFTW*******/
 	//#pragma omp critical
 	{
-        plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localproj,mapFFT,FFTW_FORWARD,FFTW_ESTIMATE);
-        fftw_execute(plan);
-        fftw_destroy_plan(plan);
-        fftw_free(localproj);
+		plan = myfftw_plan_dft_2d(param.param_device.NumberPixels,param.param_device.NumberPixels,localproj,mapFFT,FFTW_FORWARD,FFTW_ESTIMATE);
+
+		myfftw_execute(plan);
+		myfftw_destroy_plan(plan);
+		myfftw_free(localproj);
 	}
 
 	return(0);
 }
 
-
-int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv)
+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 **********
@@ -463,31 +657,28 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
 	*************** and Backtransforming it to real Space ******************************/
 	/**************************************************************************************/
 
-    fftw_plan plan;
-    mycomplex_t* localmultFFT;
+	myfftw_plan plan;
 	mycomplex_t* localconvFFT;
-    int totnumPixFFT=2*param.param_device.NumberPixels;
-    localmultFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
-    localconvFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
+	localconvFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t)*param.param_device.NumberPixels*param.param_device.NumberPixels);
 
 
 	/**** Multiplying FFTmap with corresponding kernel ****/
 
-    for(int i=0; i < 2*param.param_device.NumberPixels ; i++ )
+	for(int i=0; i < param.param_device.NumberPixels ; i++ )
 	{
-        for(int j=0; j < 2*param.param_device.NumberPixels ; j++ )
+		for(int j=0; j < param.param_device.NumberPixels ; j++ )
 		{   //Projection*CONJ(KERNEL)
-            localmultFFT[i*2*param.param_device.NumberPixels+j][0]=lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]+lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
-            localmultFFT[i*2*param.param_device.NumberPixels+j][1]=lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]-lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
-            // cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][1] << "\n";
+			localmultFFT[i*param.param_device.NumberPixels+j][0]=lproj[i*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0]+lproj[i*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1];
+			localmultFFT[i*param.param_device.NumberPixels+j][1]=lproj[i*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0]-lproj[i*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1];
+			// cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*param.param_device.NumberPixels+j][1] << "\n";
 		}
 	}
 
 	/**** Bringing convoluted Map to real Space ****/
 	//#pragma omp critical
 	{
-        plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localmultFFT,localconvFFT,FFTW_BACKWARD,FFTW_ESTIMATE);
-        fftw_execute(plan);
+		plan = myfftw_plan_dft_2d(param.param_device.NumberPixels,param.param_device.NumberPixels,localmultFFT,localconvFFT,FFTW_BACKWARD,FFTW_ESTIMATE);
+		myfftw_execute(plan);
 	}
 
 
@@ -496,16 +687,31 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
 	{
 		for(int j=0; j < param.param_device.NumberPixels ; j++ )
 		{
-            Mapconv.points[i][j]=localconvFFT[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0];
+			Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j][0];
+		}
+	}
+
+	/*** Calculating Cross-correlations of cal-convoluted map with its self *****/
+	sumC=0;
+	sumsquareC=0;
+	for(int i=0; i < param.param_device.NumberPixels ; i++ )
+	{
+		for(int j=0; j < param.param_device.NumberPixels ; j++ )
+		{
+			sumC+=localconvFFT[i*param.param_device.NumberPixels+j][0];
+			sumsquareC+=localconvFFT[i*param.param_device.NumberPixels+j][0]*localconvFFT[i*param.param_device.NumberPixels+j][0];
 		}
 	}
+	/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***/
+	// Normalizing
+	sumC=sumC/float(param.param_device.NumberPixels*param.param_device.NumberPixels);
+	sumsquareC=sumsquareC/pow(float(param.param_device.NumberPixels),4);
 
 	/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
 	//#pragma omp critical
 	{
-        fftw_destroy_plan(plan);
-        fftw_free(localconvFFT);
-        fftw_free(localmultFFT);
+		myfftw_destroy_plan(plan);
+		myfftw_free(localconvFFT);
 	}
 
 	return(0);

include/bioem.h

@@ -21,13 +21,19 @@ public:
 	int dopreCalCrossCorrelation(int iRefMap, int iRefMapLocal);
 	int run();
 	int doProjections(int iMap);
-    int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj,bioem_map& Mapconv);
+	int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj,bioem_map& Mapconv,mycomplex_t* localmultFFT,myfloat_t& sumC,myfloat_t& sumsquareC);
+
 
 	virtual int compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, const int startMap = 0);
+	int compareRefMaps2(int iProjectionOut, int iConv,mycomplex_t* localmultFFT,myfloat_t sumC,myfloat_t sumsquareC);
+
 	int createProjection(int iMap, mycomplex_t* map);
 	int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
+	int calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC);
+	int calculateCCFFT(int iMap, int iConv, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC);
 
 	bioem_Probability* pProb;
+	bioem_crossCor* crossCor;
 
 protected:
 	virtual int deviceInit();
@@ -43,6 +49,8 @@ protected:
 
 	int nProjectionMaps;		//Maps in memory at a time
 	int nProjectionMapsTotal;	//Maps in total
+
+	int FFTAlgo;
 };
 
 #endif

include/defs.h

 #ifndef BIOEM_DEFS_H
 #define BIOEM_DEFS_H
 
+#ifndef BIOEM_USE_DOUBLE
 typedef float myfloat_t;
-typedef double mycomplex_t[2];
+#define myfftw_malloc fftwf_malloc
+#define myfftw_free fftwf_free
+#define myfftw_destroy_plan fftwf_destroy_plan
+#define myfftw_execute fftwf_execute
+#define myfftw_plan_dft_2d fftwf_plan_dft_2d
+#define myfftw_plan fftwf_plan
+#else
+typedef double myfloat_t;
+#define myfftw_malloc fftw_malloc
+#define myfftw_free fftw_free
+#define myfftw_destroy_plan fftw_destroy_plan
+#define myfftw_execute fftw_execute
+#define myfftw_plan_dft_2d fftw_plan_dft_2d
+#define myfftw_plan fftw_plan
+#endif
+typedef myfloat_t mycomplex_t[2];
 
 #define BIOEM_FLOAT_3_PHYSICAL_SIZE 3	//Possible set to 4 for GPU
 #define BIOEM_MAP_SIZE_X 224
@@ -11,6 +27,7 @@ typedef double mycomplex_t[2];
 #define BIOEM_MAX_MAPS 12000
 #define MAX_REF_CTF 200
 #define MAX_ORIENT 20000
+#define MAX_DISPLACE 224
 
 struct myfloat3_t
 {

include/map.h

@@ -16,14 +16,25 @@ public:
 class bioem_map_forFFT
 {
 public:
-    mycomplex_t cpoints[2*BIOEM_MAP_SIZE_X*2*BIOEM_MAP_SIZE_Y];
+	mycomplex_t cpoints[BIOEM_MAP_SIZE_X*BIOEM_MAP_SIZE_Y];
 
 };
 
+class bioem_convolutedMap
+{
+public:
+	bioem_map conv[MAX_REF_CTF];
+	myfloat_t sum_convMap[MAX_REF_CTF];
+	myfloat_t sumsquare_convMap[MAX_REF_CTF];
+	myfloat_t ForLogProbfromConv[MAX_REF_CTF];
+};
+
 class bioem_RefMap
 {
 public:
-	int readRefMaps();
+	int readRefMaps(bioem_param& param);
+	int PreCalculateMapsFFT(bioem_param& param);
+	bioem_map_forFFT* RefMapFFT;
 
 	const char* filemap;
 	int ntotRefMap;
@@ -74,6 +85,13 @@ public:
 	}
 };
 
+class bioem_crossCor
+{
+public:
+	int disx[MAX_DISPLACE];
+	int disy[MAX_DISPLACE];
+	myfloat_t value[MAX_DISPLACE];
+};
 
 class bioem_Probability
 {

include/model.h

@@ -15,6 +15,7 @@ public:
 
 	myfloat_t getAminoAcidRad(char *name);
 	myfloat_t getAminoAcidDensity(char *name);
+	myfloat_t NormDen;
 
 	const char* filemodel;
 

include/param.h

@@ -14,6 +14,7 @@ public:
 	int maxDisplaceCenter;
 	int GridSpaceCenter;
 	int NumberPixels;
+	int NtotDist;
 	myfloat_t Ntotpi;
 	myfloat_t volu;
 };
@@ -40,7 +41,7 @@ public:
 	bool writeAngles;
 // Pixel size &&  BIOEM_MAP_SIZE_X should be defined here too
 	//int NumberPixels; //in device class
-    myfloat_t* pixelSize;
+	myfloat_t pixelSize;
 // Grid Points in Euler angles, assuming uniform sampling d_alpha=d_gamma (in 2pi) & cos(beta)=-1,1
 	int angleGridPointsAlpha;
 	int angleGridPointsBeta;
@@ -50,16 +51,16 @@ public:
 	//int GridSpaceCenter; //in device class
 // Grid sampling for the convolution kernel
 	//ENVELOPE
-    myfloat_t* startGridEnvelop;
+	myfloat_t startGridEnvelop;
 	int numberGridPointsEnvelop;
-    myfloat_t* gridEnvelop;
+	myfloat_t gridEnvelop;
 	//CTF=Amp*cos(phase*x)-sqrt(1-Amp**2)*sin(phase*x)
-    myfloat_t* startGridCTF_phase;
+	myfloat_t startGridCTF_phase;
 	int numberGridPointsCTF_phase;
-    myfloat_t* gridCTF_phase;
-    myfloat_t* startGridCTF_amp;
+	myfloat_t gridCTF_phase;
+	myfloat_t startGridCTF_amp;
 	int numberGridPointsCTF_amp;
-    myfloat_t* gridCTF_amp;
+	myfloat_t gridCTF_amp;
 	// Others
 	//myfloat_t volu;//in device class
 	myfloat3_t angles[MAX_ORIENT];

map.cpp

@@ -3,12 +3,16 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <cstring>
+#include <math.h>
+#include <fftw3.h>
 
 #include "map.h"
+#include "param.h"
+
 
 using namespace std;
 
-int bioem_RefMap::readRefMaps()
+int bioem_RefMap::readRefMaps(bioem_param& param)
 {
 	/**************************************************************************************/
 	/***********************Reading reference Particle Maps************************/
@@ -36,6 +40,8 @@ int bioem_RefMap::readRefMaps()
 	else
 	{
 		int nummap=-1;
+		int lasti=0;
+		int lastj=0;
 		ifstream input(filemap);