bioem.cpp 22.4 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[])
{
David Rohr's avatar
David Rohr committed
50
51
52
53
54
	// **************************************************************************************
	// **** Configuration Routine using boost for extracting parameters, models and maps ****
	// **************************************************************************************
	// ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
	// *************************************************************************************
55

David Rohr's avatar
David Rohr committed
56
	// *** Inizialzing default variables ***
57
58
59
	std::string infile, modelfile, mapfile;
	Model.readPDB = false;
	param.writeAngles = false;
60
61
	RefMap.dumpMap = false;
	RefMap.loadMap = false;
62
	RefMap.readMRC = false;
63
	RefMap.readMultMRC = false;
64

65

David Rohr's avatar
David Rohr committed
66
	// *************************************************************************************
Pilar Cossio's avatar
Pilar Cossio committed
67
	cout << "+++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
David Rohr's avatar
David Rohr committed
68
	// *************************************************************************************
69

David Rohr's avatar
David Rohr committed
70
	// ********************* Command line reading input with BOOST ************************
71

72
73
74
	try {
		po::options_description desc("Command line inputs");
		desc.add_options()
Pilar Cossio's avatar
Pilar Cossio committed
75
76
77
		("Inputfile", po::value<std::string>(), "(Mandatory) Name of input parameter file")
		("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
		("Particlesfile", po::value< std::string>(), "(Mandatory) Name of paricles file")
78
		("ReadPDB", "(Optional) If reading model file in PDB format")
79
80
		("ReadMRC", "(Optional) If reading particle file in MRC format")
		("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
81
82
		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
83
84
85
86
87
88
89
90
		("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);
91
92
		p.add("ReadMRC", -1);
		p.add("ReadMultipleMRC", -1);
93
94
95
		p.add("DumpMaps", -1);
		p.add("LoadMapDump", -1);

96
97
98
99
100
101
		po::variables_map vm;
		po::store(po::command_line_parser(ac, av).
				  options(desc).positional(p).run(), vm);
		po::notify(vm);

		if((ac < 6)) {
102
			std::cout << desc << std::endl;
Pilar Cossio's avatar
Pilar Cossio committed
103
			exit(1);
104
			return 0;
105
106
107
108
109
110
111
112
113
114
		}
		if (vm.count("help")) {
			cout << "Usage: options_description [options]\n";
			cout << desc;
			return 0;
		}

		if (vm.count("Inputfile"))
		{
			cout << "Input file is: ";
115
116
			cout << vm["Inputfile"].as< std::string >() << "\n";
			infile = vm["Inputfile"].as< std::string >();
117
118
119
120
121
		}
		if (vm.count("Modelfile"))
		{
			cout << "Model file is: "
				 << vm["Modelfile"].as<  std::string  >() << "\n";
122
			modelfile = vm["Modelfile"].as<  std::string  >();
123
124
125
126
127
		}

		if (vm.count("ReadPDB"))
		{
			cout << "Reading model file in PDB format.\n";
128
			Model.readPDB = true;
129
		}
130
		if (vm.count("ReadMRC"))
131
132
133
134
135
		{
			cout << "Reading particle file in MRC format.\n";
			RefMap.readMRC = true;
		}

136
		if (vm.count("ReadMultipleMRC"))
137
138
139
140
		{
			cout << "Reading Multiple MRCs.\n";
			RefMap.readMultMRC = true;
		}
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

		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";
158
			mapfile = vm["Particlesfile"].as< std::string >();
159
160
161
162
163
164
165
166
		}
	}
	catch(std::exception& e)
	{
		cout << e.what() << "\n";
		return 1;
	}

David Rohr's avatar
David Rohr committed
167
	// ********************* Reading Parameter Input ***************************
168
169
170
171
	// copying inputfile to param class
	param.fileinput = infile.c_str();
	param.readParameters();

David Rohr's avatar
David Rohr committed
172
	// ********************* Reading Model Input ******************************
173
174
175
176
	// copying modelfile to model class
	Model.filemodel = modelfile.c_str();
	Model.readModel();

David Rohr's avatar
David Rohr committed
177
178
	// ********************* Reading Particle Maps Input **********************
	// ********* HERE: PROBLEM if maps dont fit on the memory!! ***************
179
180
181
182
	// copying mapfile to ref map class
	RefMap.filemap = mapfile.c_str();
	RefMap.readRefMaps(param);

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

186
187
188
189
190
	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
191

192
193
	deviceInit();

194
	return(0);
195
196
197
198
}

int bioem::precalculate()
{
David Rohr's avatar
David Rohr committed
199
	// **************************************************************************************
Pilar Cossio's avatar
Pilar Cossio committed
200
	//* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels *
David Rohr's avatar
David Rohr committed
201
	// **************************************************************************************
202

203
204
	// Generating Grids of orientations
	param.CalculateGridsParam();
205

206
	myfloat_t sum, sumsquare;
207
208
209
210

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

217
218
	// Precalculating CTF Kernels stored in class Param
	param.CalculateRefCTF();
219

220
221
	// Precalculating Maps in Fourier space
	RefMap.PreCalculateMapsFFT(param);
222

223
	return(0);
224
225
226
227
228
}


int bioem::run()
{
David Rohr's avatar
David Rohr committed
229
230
231
	// **************************************************************************************
	// **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
	// **************************************************************************************
232

David Rohr's avatar
David Rohr committed
233
234
	// **** If we want to control the number of threads -> omp_set_num_threads(XX); ******
	// ****************** Declarying class of Probability Pointer  *************************
235
236
237
238
239
240
	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 ++)
	{
241
242
243
		pProb[iRefMap].Total = 0.0;
		pProb[iRefMap].Constoadd = -9999999;
		pProb[iRefMap].max_prob = -9999999;
244
		for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
245
		{
246
247
			pProb[iRefMap].forAngles[iOrient] = 0.0;
			pProb[iRefMap].ConstAngle[iOrient] = -99999999;
248
249
		}
	}
David Rohr's avatar
David Rohr committed
250
	// **************************************************************************************
251
252
	deviceStartRun();

David Rohr's avatar
David Rohr committed
253
	// ******************************** MAIN CYCLE ******************************************
David Rohr's avatar
David Rohr committed
254

David Rohr's avatar
David Rohr committed
255
	// *** Declaring Private variables for each thread *****
256
257
258
	mycomplex_t* proj_mapFFT;
	bioem_map conv_map;
	mycomplex_t* conv_mapFFT;
259
	myfloat_t sumCONV, sumsquareCONV;
260
261

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

	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)));
269
270
	for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
	{
David Rohr's avatar
David Rohr committed
271
272
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
273
		timer.ResetStart();
274
		createProjection(iProjectionOut, proj_mapFFT);
275
276
		printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());

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

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

David Rohr's avatar
David Rohr committed
288
289
			// ***************************************************************************************
			// *** Comparing each calculated convoluted map with all experimental maps ***
290
			timer.ResetStart();
291
			compareRefMaps(iProjectionOut, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
292

293
294
295
			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 *
296
								  (((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;
297
			const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
298
								(((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;
299
300
301
			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.);
302
303
304
305
306
		}
	}
	//deallocating fftw_complex vector
	myfftw_free(proj_mapFFT);
	myfftw_free(conv_mapFFT);
David Rohr's avatar
David Rohr committed
307

308
309
	deviceFinishRun();

David Rohr's avatar
David Rohr committed
310
	// ************* Writing Out Probabilities ***************
311

David Rohr's avatar
David Rohr committed
312
	// *** Angular Probability ***
313

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

319
320
	ofstream outputProbFile;
	outputProbFile.open ("Output_Probabilities");
321

322
323
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
David Rohr's avatar
David Rohr committed
324
		// **** Total Probability ***
325
		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";
326
327
328

		outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";

David Rohr's avatar
David Rohr committed
329
		// *** Param that maximize probability****
330
331
332
333
334
335
336
		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] << " ";
337
338
		outputProbFile << pProb[iRefMap].max_prob_cent_x << " ";
		outputProbFile << pProb[iRefMap].max_prob_cent_y;
339
		outputProbFile << "\n";
340

David Rohr's avatar
David Rohr committed
341
		// *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
342

343
		if(param.writeAngles)
344
		{
345
346
			for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
			{
347
				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";
348

349
350
351
			}
		}
	}
352

353
354
	angProbfile.close();
	outputProbFile.close();
355

356
	//Deleting allocated pointers
357

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

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

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

378
int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
379
{
Pilar Cossio's avatar
Pilar Cossio committed
380

David Rohr's avatar
David Rohr committed
381
382
	// ***************************************************************************************
	// ***** BioEM routine for comparing reference maps to convoluted maps *****
Pilar Cossio's avatar
Pilar Cossio committed
383
384

	if (FFTAlgo) //IF using the fft algorithm
385
	{
386
		#pragma omp parallel
387
388
389
		{
			mycomplex_t *localCCT;
			myfloat_t *lCC;
390
391
			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);
392
393
394
395
396
397
398
399
400

			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 ++)
			{
401
				calculateCCFFT(iRefMap, iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT, lCC);
402
403
404
405
406
407
			}
			myfftw_free(localCCT);
			myfftw_free(lCC);
		}
	}
	else
408
	{
409
		#pragma omp parallel for
410
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
411
		{
412
			compareRefMapShifted < -1 > (iRefMap, iProjectionOut, iConv, conv_map, pProb, param.param_device, RefMap);
413
414
415
416
417
		}
	}
	return(0);
}

418
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)
419
{
Pilar Cossio's avatar
Pilar Cossio committed
420

David Rohr's avatar
David Rohr committed
421
422
	// ***************************************************************************************
	// ***** Calculating cross correlation in FFTALGOrithm *****
Pilar Cossio's avatar
Pilar Cossio committed
423

424
	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize];
425
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
426
	{
427
428
		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];
429
430
	}

431
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
432

433
	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
434
}
435

436
int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
437
{
David Rohr's avatar
David Rohr committed
438
439
440
441
	// **************************************************************************************
	// ****  BioEM Create Projection routine in Euler angle predefined grid****************
	// ********************* and turns projection into Fourier space **********************
	// **************************************************************************************
442
443
444

	myfloat3_t RotatedPointsModel[Model.nPointsModel];
	myfloat_t rotmat[3][3];
445
	myfloat_t alpha, gam, beta;
446
	myfloat_t* localproj;
447

448
449
	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));
450

451
452
453
	alpha = param.angles[iMap].pos[0];
	beta = param.angles[iMap].pos[1];
	gam = param.angles[iMap].pos[2];
454

David Rohr's avatar
David Rohr committed
455
	// **** To see how things are going: cout << "Id " << omp_get_thread_num() <<  " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
456

David Rohr's avatar
David Rohr committed
457
	// ********** Creat Rotation with pre-defiend grid of orientations**********
458

459
460
461
462
463
464
465
466
467
468
469
	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++)
470
	{
471
472
473
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
474
	}
475
	for(int n = 0; n < Model.nPointsModel; n++)
476
	{
477
		for(int k = 0; k < 3; k++)
478
		{
479
			for(int j = 0; j < 3; j++)
480
			{
481
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.PointsModel[n].pos[j];
482
483
484
485
486
487
			}
		}
	}

	int i, j;

David Rohr's avatar
David Rohr committed
488
	// ************ Projection over the Z axis********************
489
	for(int n = 0; n < Model.nPointsModel; n++)
490
491
	{
		//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
492
493
		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);
494

495
		localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
496
497
	}

David Rohr's avatar
David Rohr committed
498
	// **** Output Just to check****
499
	if(iMap == 10)
500
501
502
503
504
505
	{
		ofstream myexamplemap;
		ofstream myexampleRot;
		myexamplemap.open ("MAP_i10");
		myexampleRot.open ("Rot_i10");
		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
506
		for(int k = 0; k < param.param_device.NumberPixels; k++)
507
		{
508
			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
509
510
		}
		myexamplemap << " \n";
511
		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
512
513
514
515
		myexamplemap.close();
		myexampleRot.close();
	}

David Rohr's avatar
David Rohr committed
516
517
	// ***** Converting projection to Fourier Space for Convolution later with kernel****
	// ********** Omp Critical is necessary with FFTW*******
518
	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
519
520
521
522

	return(0);
}

523
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
524
{
David Rohr's avatar
David Rohr committed
525
526
527
528
529
	// **************************************************************************************
	// ****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
	// **************** calculated Projection with convoluted precalculated Kernel**********
	// *************** and Backtransforming it to real Space ******************************
	// **************************************************************************************
530

531
	myfloat_t* localconvFFT;
532
	localconvFFT = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
533
534
	mycomplex_t* tmp;
	tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
535

David Rohr's avatar
David Rohr committed
536
	// **** Multiplying FFTmap with corresponding kernel ****
537

538
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize];
539
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
540
	{
541
542
543
		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";
544
545
	}

546
547
548
	//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);

David Rohr's avatar
David Rohr committed
549
	// **** Bringing convoluted Map to real Space ****
550
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, localconvFFT);
551

David Rohr's avatar
David Rohr committed
552
	// ****Asigning convolution fftw_complex to bioem_map ****
553
	for(int i = 0; i < param.param_device.NumberPixels ; i++ )
554
	{
555
		for(int j = 0; j < param.param_device.NumberPixels ; j++ )
556
		{
557
			Mapconv.points[i][j] = localconvFFT[i * param.param_device.NumberPixels + j];
558
559
560
		}
	}

David Rohr's avatar
David Rohr committed
561
	// *** Calculating Cross-correlations of cal-convoluted map with its self *****
562
563
	sumC = 0;
	sumsquareC = 0;
564
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
565
	{
566
567
		sumC += localconvFFT[i];
		sumsquareC += localconvFFT[i] * localconvFFT[i];
568
	}
David Rohr's avatar
David Rohr committed
569
	// *** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
570
	// Normalizing
571
572
573
574
	myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
	myfloat_t norm4 = norm2 * norm2;
	sumC = sumC / norm2;
	sumsquareC = sumsquareC / norm4;
575

David Rohr's avatar
David Rohr committed
576
	// **** Freeing fftw_complex created (dont know if omp critical is necessary) ****
577
	myfftw_free(localconvFFT);
578
	myfftw_free(tmp);
579
580

	return(0);
581
582
}

583
int bioem::calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare)
584
{
David Rohr's avatar
David Rohr committed
585
	// *********************** Routine to calculate Cross correlations***********************
586

587
588
	sum = 0.0;
	sumsquare = 0.0;
589
590
591
592
593
594
595
	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
596
			sumsquare += localmap.points[i][j] * localmap.points[i][j];
597
598
599
		}
	}
	return(0);
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
}

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

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

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