From e7310deaf61fb85b058657a0c11ec99c9378caa8 Mon Sep 17 00:00:00 2001
From: David Rohr <drohr@jwdt.org>
Date: Mon, 7 Jul 2014 18:59:05 +0200
Subject: [PATCH] add printModel flag, get ride of unnecessary fourier
 transform for FFTALGO

---
 bioem.cpp       | 57 +++++++++++++++++++++++++++++++++++--------------
 include/param.h |  4 ++++
 param.cpp       |  2 ++
 3 files changed, 47 insertions(+), 16 deletions(-)

diff --git a/bioem.cpp b/bioem.cpp
index 51156b1..0f268a2 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -283,6 +283,13 @@ int bioem::run()
 	// **************************************************************************************
 	// **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
 	// **************************************************************************************
+	
+	if (param.printModel)
+	{
+		
+		//....
+		return(0);
+	}
 
 	// **** If we want to control the number of threads -> omp_set_num_threads(XX); ******
 	// ****************** Declarying class of Probability Pointer  *************************
@@ -720,27 +727,45 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
 		// cout << "GG " << i << " " << j << " " << refCTF[i][0] << " " << refCTF[i][1] <<" " <<lproj[i][0] <<" " <<lproj[i][1] << "\n";
 	}
 
-	//FFTW_C2R will destroy the input array, so we have to work on a copy here
-	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, Mapconv);
-
 	// *** Calculating Cross-correlations of cal-convoluted map with its self *****
-	sumC = 0;
+	sumC = localmultFFT[0][0];
+	
 	sumsquareC = 0;
-	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
+	if (FFTAlgo)
 	{
-		sumC += Mapconv[i];
-		sumsquareC += Mapconv[i] * Mapconv[i];
+		for(int i = 0; i < param.param_device.NumberPixels; i++)
+		{
+			for (int j = 1;j < param.param_device.NumberFFTPixels1D;j++)
+			{
+				int k = i * param.param_device.NumberFFTPixels1D + j;
+				sumsquareC += (localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1]) * 2;
+			}
+			int k = i * param.param_device.NumberFFTPixels1D;
+			sumsquareC += localmultFFT[k][0] * localmultFFT[k][0] + localmultFFT[k][1] * localmultFFT[k][1];
+		}
+		
+		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
+		sumsquareC = sumsquareC / norm2;
 	}
-	// *** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
-	// Normalizing
-	myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
-	myfloat_t norm4 = norm2 * norm2;
-	sumC = sumC / norm2;
-	sumsquareC = sumsquareC / norm4;
+	else
+	{
+		//FFTW_C2R will destroy the input array, so we have to work on a copy here
+		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, Mapconv);
 
+		for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
+		{
+			sumsquareC += Mapconv[i] * Mapconv[i];
+			sumC += Mapconv[i];
+		}
+
+		myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
+		myfloat_t norm4 = norm2 * norm2;
+		sumsquareC = sumsquareC / norm4;
+	}
+	
 	cuda_custom_timeslot_end;
 
 	return(0);
diff --git a/include/param.h b/include/param.h
index e395831..b251b81 100644
--- a/include/param.h
+++ b/include/param.h
@@ -64,6 +64,10 @@ public:
 	myfloat3_t* angles;
 	int nTotGridAngles;
 	int nTotCTFs;
+	
+	bool printModel;
+	int printModelOrientation;
+	int printModelConvolution;
 
 	int fft_plans_created;
 	myfftw_plan fft_plan_c2c_forward, fft_plan_c2c_backward, fft_plan_r2c_forward, fft_plan_c2r_backward;
diff --git a/param.cpp b/param.cpp
index db86fd0..a2f8796 100644
--- a/param.cpp
+++ b/param.cpp
@@ -43,6 +43,8 @@ bioem_param::bioem_param()
 	refCTF = NULL;
 	CtfParam = NULL;
 	angles = NULL;
+	
+	printModel = false;
 }
 
 int bioem_param::readParameters(const char* fileinput)
-- 
GitLab