bioem.cpp 21.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include <fstream>
#include <boost/program_options.hpp>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
#include <omp.h>

#include <fftw3.h>
#include <math.h>
#include "cmodules/timer.h"

#include "param.h"
#include "bioem.h"
#include "model.h"
#include "map.h"


#include "bioem_algorithm.h"


using namespace boost;
namespace po = boost::program_options;

using namespace std;

// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
34
35
	copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
	return os;
36
37
38
39
}

bioem::bioem()
{
40
	FFTAlgo = getenv("FFTALGO") == NULL ? 0 : atoi(getenv("FFTALGO"));
41
42
43
44
}

bioem::~bioem()
{
David Rohr's avatar
David Rohr committed
45

46
47
48
49
}

int bioem::configure(int ac, char* av[])
{
50
51
52
53
54
55
56
57
58
59
	/**************************************************************************************/
	/**** Configuration Routine using boost for extracting parameters, models and maps ****/
	/**************************************************************************************/
	/****** And Precalculating necessary grids, map crosscorrelations and kernels  ********/
	/*************************************************************************************/

	/*** Inizialzing default variables ***/
	std::string infile,modelfile,mapfile;
	Model.readPDB=false;
	param.writeAngles=false;
60
61
	RefMap.dumpMap = false;
	RefMap.loadMap = false;
62
63
        RefMap.readMRC = false;
	RefMap.readMultMRC = false;
64

65
66
67
	/*************************************************************************************/
	cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
	/*************************************************************************************/
68

69
	/********************* Command line reading input with BOOST ************************/
70

71
72
73
74
75
76
77
	try {
		po::options_description desc("Command line inputs");
		desc.add_options()
		("Inputfile", po::value<std::string>(), "Name of input parameter file")
		("Modelfile", po::value< std::string>() , "Name of model file")
		("Particlesfile", po::value< std::string>(), "Name of paricles file")
		("ReadPDB", "(Optional) If reading model file in PDB format")
78
79
                ("ReadMRC", "(Optional) If reading particle file in MRC format")
		 ("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
80
81
		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
82
83
84
85
86
87
88
89
		("help", "(Optional) Produce help message")
		;

		po::positional_options_description p;
		p.add("Inputfile", -1);
		p.add("Modelfile", -1);
		p.add("Particlesfile", -1);
		p.add("ReadPDB", -1);
90
91
                p.add("ReadMRC", -1);
                p.add("ReadMultipleMRC", -1);
92
93
94
		p.add("DumpMaps", -1);
		p.add("LoadMapDump", -1);

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
		po::variables_map vm;
		po::store(po::command_line_parser(ac, av).
				  options(desc).positional(p).run(), vm);
		po::notify(vm);

		if((ac < 6)) {
			std::cout << desc << std::endl;
			return 0;
		}
		if (vm.count("help")) {
			cout << "Usage: options_description [options]\n";
			cout << desc;
			return 0;
		}

		if (vm.count("Inputfile"))
		{
			cout << "Input file is: ";
			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  >();
		}

		if (vm.count("ReadPDB"))
		{
			cout << "Reading model file in PDB format.\n";
			Model.readPDB=true;
		}
128
129
130
131
132
133
134
135
136
137
138
		if (vm.count("ReadMRC"))
                {
                        cout << "Reading particle file in MRC format.\n";
                        RefMap.readMRC=true;
                }
                
		if (vm.count("ReadMultipleMRC"))
                {
                        cout << "Reading Multiple MRCs.\n";
                        RefMap.readMultMRC=true;
                }
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

		if (vm.count("DumpMaps"))
		{
			cout << "Dumping Maps after reading from file.\n";
			RefMap.dumpMap = true;
		}

		if (vm.count("LoadMapDump"))
		{
			cout << "Loading Map dump.\n";
			RefMap.loadMap = true;
		}

		if (vm.count("Particlesfile"))
		{
			cout << "Paricle file is: "
				 << vm["Particlesfile"].as< std::string >() << "\n";
			mapfile=vm["Particlesfile"].as< std::string >();
		}
	}
	catch(std::exception& e)
	{
		cout << e.what() << "\n";
		return 1;
	}

	/********************* Reading Parameter Input ***************************/
	// copying inputfile to param class
	param.fileinput = infile.c_str();
	param.readParameters();

	/********************* Reading Model Input ******************************/
	// copying modelfile to model class
	Model.filemodel = modelfile.c_str();
	Model.readModel();

	/********************* Reading Particle Maps Input **********************/
	/********* HERE: PROBLEM if maps dont fit on the memory!! ***************/
	// copying mapfile to ref map class
	RefMap.filemap = mapfile.c_str();
	RefMap.readRefMaps(param);

	/****************** Precalculating Necessary Stuff *********************/
	precalculate();
David Rohr's avatar
David Rohr committed
183

184
185
186
187
188
	if (getenv("BIOEM_DEBUG_BREAK"))
	{
		param.nTotGridAngles = atoi(getenv("BIOEM_DEBUG_BREAK"));
		param.nTotCTFs = atoi(getenv("BIOEM_DEBUG_BREAK"));
	}
David Rohr's avatar
David Rohr committed
189

190
191
	deviceInit();

192
	return(0);
193
194
195
196
}

int bioem::precalculate()
{
197
198
199
	/**************************************************************************************/
	/* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels */
	/**************************************************************************************/
200

201
202
	// Generating Grids of orientations
	param.CalculateGridsParam();
203

204
205
206
207
208
209
210
211
212
213
	myfloat_t sum,sumsquare;

	//Precalculating cross-correlations of maps
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
	{
		calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
		//Storing Crosscorrelations in Map class
		RefMap.sum_RefMap[iRefMap]=sum;
		RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
	}
214

215
216
	// Precalculating CTF Kernels stored in class Param
	param.CalculateRefCTF();
217

218
219
	// Precalculating Maps in Fourier space
	RefMap.PreCalculateMapsFFT(param);
220

221
	return(0);
222
223
224
225
226
}


int bioem::run()
{
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
	/**************************************************************************************/
	/**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****/
	/**************************************************************************************/

	/**** 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];

	printf("\tInitializing\n");
	// 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;
		for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
243
		{
244
245
246
247
248
			pProb[iRefMap].forAngles[iOrient]=0.0;
			pProb[iRefMap].ConstAngle[iOrient]=-99999999;
		}
	}
	/**************************************************************************************/
249
250
	deviceStartRun();

251
	/******************************** MAIN CYCLE ******************************************/
David Rohr's avatar
David Rohr committed
252

253
254
255
256
257
	/*** Declaring Private variables for each thread *****/
	mycomplex_t* proj_mapFFT;
	bioem_map conv_map;
	mycomplex_t* conv_mapFFT;
	myfloat_t sumCONV,sumsquareCONV;
258
259

	//allocating fftw_complex vector
260
261
	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);
262
263
264
265
266

	HighResTimer timer;

	printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %d²\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels);
	printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1)));
267
268
269
270
	for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
	{
		/***************************************************************************************/
		/***** Creating Projection for given orientation and transforming to Fourier space *****/
271
		timer.ResetStart();
272
		createProjection(iProjectionOut, proj_mapFFT);
273
274
		printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());

275
276
277
278
		/***************************************************************************************/
		/***** **** Internal Loop over convolutions **** *****/
		for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
		{
279
			printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
280
281
			/*** Calculating convolutions of projection map and crosscorrelations ***/

282
			timer.ResetStart();
283
			createConvolutedProjectionMap(iProjectionOut, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
284
285
			printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());

286
287
			/***************************************************************************************/
			/*** Comparing each calculated convoluted map with all experimental maps ***/
288
			timer.ResetStart();
289
			compareRefMaps(iProjectionOut, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
290

291
292
293
			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 *
294
								  (((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 5. + 25.) / compTime;
295
			const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
296
								(((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 2. + 8.) * (double) sizeof(myfloat_t) / compTime;
297
298
299
			const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

			printf("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iProjectionOut, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
300
301
302
303
304
		}
	}
	//deallocating fftw_complex vector
	myfftw_free(proj_mapFFT);
	myfftw_free(conv_mapFFT);
David Rohr's avatar
David Rohr committed
305

306
307
	deviceFinishRun();

308
	/************* Writing Out Probabilities ***************/
309

310
	/*** Angular Probability ***/
311

312
313
314
315
	// if(param.writeAngles){
	ofstream angProbfile;
	angProbfile.open ("ANG_PROB_iRefMap");
	// }
316

317
318
	ofstream outputProbFile;
	outputProbFile.open ("Output_Probabilities");
319

320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
	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 << " Maximizing Param: ";

		/*** Param that maximize probability****/
		outputProbFile << (pProb[iRefMap].max_prob + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " ";
		outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[0] << " ";
		outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " ";
		outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[2] << " ";
		outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[0] << " ";
		outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[1] << " ";
		outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[2] << " ";
335
336
		outputProbFile << pProb[iRefMap].max_prob_cent_x << " ";
		outputProbFile << pProb[iRefMap].max_prob_cent_y;
337
		outputProbFile << "\n";
338

339
		/*** For individual files***/ //angProbfile.open ("ANG_PROB_"iRefMap);
340

341
		if(param.writeAngles)
342
		{
343
344
345
			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";
346

347
348
349
			}
		}
	}
350

351
352
	angProbfile.close();
	outputProbFile.close();
353

354
	//Deleting allocated pointers
355

356
357
358
359
360
361
362
363
364
365
366
	if (pProb)
	{
		delete[] pProb;
		pProb = NULL;
	}

	if (param.refCTF)
	{
		delete[] param.refCTF;
		param.refCTF =NULL;
	}
367

368
	if(RefMap.RefMapsFFT)
369
	{
370
371
		delete[] RefMap.RefMapsFFT;
		RefMap.RefMapsFFT = NULL;
372
373
	}
	return(0);
374
375
}

376
int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
377
{
378
	if (FFTAlgo)
379
	{
380
#pragma omp parallel
381
382
383
384
385
386
387
388
389
390
391
392
393
394
		{
			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);

			const int num_threads = omp_get_num_threads();
			const int thread_id = omp_get_thread_num();
			const int mapsPerThread = (RefMap.ntotRefMap - startMap + num_threads - 1) / num_threads;
			const int iStart = startMap + thread_id * mapsPerThread;
			const int iEnd = min(RefMap.ntotRefMap, startMap + (thread_id + 1) * mapsPerThread);

			for (int iRefMap = iStart; iRefMap < iEnd; iRefMap ++)
			{
395
				calculateCCFFT(iRefMap,iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT,lCC);
396
397
398
399
400
401
			}
			myfftw_free(localCCT);
			myfftw_free(lCC);
		}
	}
	else
402
	{
403
404
#pragma omp parallel for
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
405
		{
406
			compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
407
408
409
410
411
		}
	}
	return(0);
}

412
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)
413
{
414
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize];
415
	for(int i = 0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
416
	{
417
418
		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];
419
420
	}

421
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,localCCT,lCC);
422

423
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
424
}
425

426
int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
427
{
428
429
430
431
432
433
434
435
	/**************************************************************************************/
	/****  BioEM Create Projection routine in Euler angle predefined grid****************
	********************* and turns projection into Fourier space **********************/
	/**************************************************************************************/

	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
	myfloat_t alpha, gam,beta;
436
	myfloat_t* localproj;
437

438
	localproj= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
	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];

	/**** 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++)
	{
		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 k=0; k< 3; k++)
		{
			for(int j=0; j< 3; j++)
			{
				RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
			}
		}
	}

	int i, j;

	/************ Projection over the Z axis********************/
	for(int n=0; n< Model.nPointsModel; n++)
	{
		//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
482
483
		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);
484

485
		localproj[i*param.param_device.NumberPixels+j]+=Model.densityPointsModel[n]/Model.NormDen;
486
487
488
489
490
491
492
493
494
495
496
497
	}

	/**** Output Just to check****/
	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++)
		{
498
			for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j];
499
500
501
502
503
504
505
506
507
		}
		myexamplemap << " \n";
		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*******/
508
	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localproj,mapFFT);
509
510
511
512

	return(0);
}

513
int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
514
515
516
517
518
519
520
{
	/**************************************************************************************/
	/****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
	**************** calculated Projection with convoluted precalculated Kernel**********
	*************** and Backtransforming it to real Space ******************************/
	/**************************************************************************************/

521
522
523
524
	myfloat_t* localconvFFT;
	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);
525
526
527

	/**** Multiplying FFTmap with corresponding kernel ****/

528
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize];
529
	for(int i=0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
530
	{
531
532
533
		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];
		// cout << "GG " << i << " " << j << " " << refCTF[i][0] << " " << refCTF[i][1] <<" " <<lproj[i][0] <<" " <<lproj[i][1] << "\n";
534
535
	}

536
537
538
	//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);

539
	/**** Bringing convoluted Map to real Space ****/
540
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,tmp,localconvFFT);
541
542
543
544
545
546

	/****Asigning convolution fftw_complex to bioem_map ****/
	for(int i=0; i < param.param_device.NumberPixels ; i++ )
	{
		for(int j=0; j < param.param_device.NumberPixels ; j++ )
		{
547
			Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j];
548
549
550
551
552
553
		}
	}

	/*** Calculating Cross-correlations of cal-convoluted map with its self *****/
	sumC=0;
	sumsquareC=0;
554
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
555
	{
556
557
		sumC += localconvFFT[i];
		sumsquareC += localconvFFT[i] * localconvFFT[i];
558
559
560
	}
	/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***/
	// Normalizing
561
562
563
564
	myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
	myfloat_t norm4 = norm2 * norm2;
	sumC = sumC / norm2;
	sumsquareC = sumsquareC / norm4;
565
566

	/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
567
	myfftw_free(localconvFFT);
568
	myfftw_free(tmp);
569
570

	return(0);
571
572
573
574
}

int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
{
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
	/*********************** Routine to calculate Cross correlations***********************/

	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.points[i][j];
			// Calculate Sum of pixels squared
			sumsquare += localmap.points[i][j]*localmap.points[i][j];
		}
	}
	return(0);
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
}

int bioem::deviceInit()
{
	return(0);
}

int bioem::deviceStartRun()
{
	return(0);
}

int bioem::deviceFinishRun()
{
	return(0);
}