New Class for printing out Cross Correlations

bioem.cpp

@@ -546,6 +546,11 @@ int bioem::run()
 		{
 			angProbfile.open ("ANG_PROB");
 		}
+		ofstream ccProbfile;
+                if(param.param_device.writeCC)
+                {
+                        ccProbfile.open ("CROSS_CORRELATION");
+                }
 
 		ofstream outputProbFile;
 		outputProbFile.open ("Output_Probabilities");
@@ -583,12 +588,31 @@ int bioem::run()
 					angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 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";
 				}
 			}
+                        if(param.param_device.writeCC)
+                        {
+				 int  cc=0;
+		           for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
+                	   {
+	                        for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
+        	                {
+	      	                bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
+                                ccProbfile << " " << iRefMap << " " << cent_x << " " << cent_y << " " << log(pProbCC.forCC) <<"\n";
+                        	}
+                	   }
+                             
+			}
 		}
 
 		if(param.param_device.writeAngles)
 		{
 			angProbfile.close();
 		}
+                if(param.param_device.writeCC)
+                {
+                        ccProbfile.close();
+                }
+
+
 		outputProbFile.close();
 	}
 

bioem_algorithm.h

@@ -114,7 +114,8 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
 			pProbAngle.ConstAngle = logpro;
 		}
 		pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle);
-	}
+	}      
+
 }
 
 __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)
@@ -141,6 +142,32 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
 			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
 		}
 	}
+
+	if (param.writeCC)
+        {
+          int  cc=0;
+	   for (int cent_x = 0; cent_x < param.NumberPixels; cent_x = cent_x + param.CCdisplace)
+        	{
+                	for (int cent_y = 0; cent_y < param.NumberPixels; cent_y = cent_y + param.CCdisplace)
+                	{     
+       
+	                bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
+
+        /*      		  if(pProbAngle.ConstAngle < logpro)
+                	{
+                        	pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle);
+	                        pProbAngle.ConstAngle = logpro;
+        		 }*/
+                        myfloat_t ttmp;
+                        ttmp=lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels);
+                        if(ttmp>0)pProbCC.forCC += exp( - ttmp);
+                        cc++;  
+                	}
+             	}
+        }
+
+
+      
 }
 
 template <int GPUAlgo, class RefT>

include/map.h

@@ -102,17 +102,27 @@ public:
 	myfloat_t ConstAngle;
 };
 
+class bioem_Probability_cc
+{
+public:
+        myfloat_t forCC;
+};
+
 class bioem_Probability
 {
 public:
 	int nMaps;
 	int nAngles;
+        int nCC;
+
 	__device__ __host__ bioem_Probability_map& getProbMap(int map) {return(ptr_map[map]);}
 	__device__ __host__ bioem_Probability_angle& getProbAngle(int map, int angle) {return(ptr_angle[angle * nMaps + map]);}
+        __device__ __host__ bioem_Probability_cc& getProbCC(int map, int cc) {return(ptr_cc[cc * nMaps + map]);}
 
 	void* ptr;
 	bioem_Probability_map* ptr_map;
 	bioem_Probability_angle* ptr_angle;
+        bioem_Probability_cc* ptr_cc;
 
 	static size_t get_size(size_t maps, size_t angles, bool writeAngles)
 	{
@@ -120,14 +130,22 @@ public:
 		if (writeAngles) size += angles * sizeof(bioem_Probability_angle);
 		return(maps * size);
 	}
+         static size_t get_sizeCC(size_t maps, size_t cc, bool writeCC)
+        {
+                size_t size = sizeof(bioem_Probability_map);
+                if (writeCC) size += cc * sizeof(bioem_Probability_cc) ;
+                return(maps * size);
+        }
 
 	void init(size_t maps, size_t angles, bioem& bio);
+        void initCC(size_t maps, size_t cc, bioem& bio); 
 	void copyFrom(bioem_Probability* from, bioem& bio);
 
 	void set_pointers()
 	{
 		ptr_map = (bioem_Probability_map*) ptr;
 		ptr_angle = (bioem_Probability_angle*) (&ptr_map[nMaps]);
+                ptr_cc = (bioem_Probability_cc*) (&ptr_map[nMaps]);
 	}
 };
 

include/param.h

@@ -16,11 +16,14 @@ public:
 	int NumberPixels;
 	int NumberFFTPixels1D;
 	int NtotDist;
+      
 	myfloat_t Ntotpi;
 	myfloat_t volu;
 // If to write Probabilities of Angles from Model
 	bool writeAngles;
-        bool writeCrosscor;
+        bool writeCC;
+	 int CCdisplace;
+
 };
 
 class bioem_param
@@ -66,6 +69,7 @@ public:
 	myfloat3_t* angles;
 	int nTotGridAngles;
 	int nTotCTFs;
+	int nTotCC;
 
 	bool printModel;
 	int printModelOrientation;

map.cpp

@@ -288,6 +288,15 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio)
 	set_pointers();
 }
 
+void bioem_Probability::initCC(size_t maps, size_t cc, bioem& bio)
+{
+        nMaps = maps;
+        nCC = cc;
+        ptr = bio.malloc_device_host(get_sizeCC(maps, cc, bio.param.param_device.writeCC));
+        set_pointers();
+}
+
+
 void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio)
 {
 	bioem_Probability_map& pProbMap = getProbMap(0);

param.cpp

@@ -220,9 +220,14 @@ int bioem_param::readParameters(const char* fileinput)
 		}
                  else if (strcmp(token, "WRITE_CROSSCOR") == 0)
                 {
-                        param_device.writeCrosscor = true;
+                        param_device.writeCC = true;
                         cout << "Writing CrossCorrelations \n";
                 }
+		else if (strcmp(token, "CROSSCOR_DISPLACE") == 0)
+                {
+                        param_device.CCdisplace=int(atoi(token));
+                        cout << "Writing Cross Correlation Displacement " <<  param_device.CCdisplace << "  \n";
+                }
                 else if (strcmp(token, "DO_AA_RADIUS") == 0)
                 {
                         doaaradius = true;
@@ -393,6 +398,7 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
 	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);
 
+        nTotCC=(myfloat_t) (param_device.NumberPixels * param_device.NumberPixels) / (myfloat_t) ( param_device.CCdisplace * param_device.CCdisplace );
 	return(0);
 
 }