bioem.cpp 22.1 KB
Newer Older
Pilar Cossio's avatar
License  
Pilar Cossio committed
1
2
3
4
5
6
7
8
9
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        < BioEM software for Bayesian inference of Electron Microscopy images>
            Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
            Max Planck Institute of Biophysics, Frankfurt, Germany.
 
                See license statement for terms of distribution.

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#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)
{
41
42
	copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
	return os;
43
44
45
46
}

bioem::bioem()
{
47
	FFTAlgo = getenv("FFTALGO") == NULL ? 0 : atoi(getenv("FFTALGO"));
48
49
50
51
}

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

53
54
55
56
}

int bioem::configure(int ac, char* av[])
{
David Rohr's avatar
David Rohr committed
57
58
59
60
61
	// **************************************************************************************
	// **** Configuration Routine using boost for extracting parameters, models and maps ****
	// **************************************************************************************
	// ****** And Precalculating necessary grids, map crosscorrelations and kernels  ********
	// *************************************************************************************
62

David Rohr's avatar
David Rohr committed
63
	// *** Inizialzing default variables ***
64
65
66
	std::string infile, modelfile, mapfile;
	Model.readPDB = false;
	param.writeAngles = false;
67
68
	param.dumpMap = false;
	param.loadMap = false;
David Rohr's avatar
David Rohr committed
69
70
	RefMap.readMRC = false;
	RefMap.readMultMRC = false;
71

David Rohr's avatar
David Rohr committed
72
	// *************************************************************************************
73
	cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
David Rohr's avatar
David Rohr committed
74
	// *************************************************************************************
75

David Rohr's avatar
David Rohr committed
76
	// ********************* Command line reading input with BOOST ************************
77

78
79
	try {
		po::options_description desc("Command line inputs");
David Rohr's avatar
David Rohr committed
80
81
82
83
84
85
86
87
88
89
90
		desc.add_options()
		("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")
		("ReadPDB", "(Optional) If reading model file in PDB format")
		("ReadMRC", "(Optional) If reading particle file in MRC format")
		("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
		("help", "(Optional) Produce help message")
		;
David Rohr's avatar
David Rohr committed
91

92
93
94
95
96
97

		po::positional_options_description p;
		p.add("Inputfile", -1);
		p.add("Modelfile", -1);
		p.add("Particlesfile", -1);
		p.add("ReadPDB", -1);
David Rohr's avatar
David Rohr committed
98
99
		p.add("ReadMRC", -1);
		p.add("ReadMultipleMRC", -1);
100
101
102
		p.add("DumpMaps", -1);
		p.add("LoadMapDump", -1);

103
104
105
106
107
108
109
		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;
110
			return 1;
111
112
113
114
		}
		if (vm.count("help")) {
			cout << "Usage: options_description [options]\n";
			cout << desc;
115
			return 1;
116
117
118
119
120
		}

		if (vm.count("Inputfile"))
		{
			cout << "Input file is: ";
121
122
			cout << vm["Inputfile"].as< std::string >() << "\n";
			infile = vm["Inputfile"].as< std::string >();
123
124
125
126
127
		}
		if (vm.count("Modelfile"))
		{
			cout << "Model file is: "
				 << vm["Modelfile"].as<  std::string  >() << "\n";
128
			modelfile = vm["Modelfile"].as<  std::string  >();
129
130
131
132
133
		}

		if (vm.count("ReadPDB"))
		{
			cout << "Reading model file in PDB format.\n";
134
			Model.readPDB = true;
135
136
		}

David Rohr's avatar
David Rohr committed
137
138
139
140
141
142
143
144
145
146
147
		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;
		}
148
149
150
151

		if (vm.count("DumpMaps"))
		{
			cout << "Dumping Maps after reading from file.\n";
152
			param.dumpMap = true;
153
154
155
156
157
		}

		if (vm.count("LoadMapDump"))
		{
			cout << "Loading Map dump.\n";
158
			param.loadMap = true;
159
160
161
162
163
164
		}

		if (vm.count("Particlesfile"))
		{
			cout << "Paricle file is: "
				 << vm["Particlesfile"].as< std::string >() << "\n";
165
			mapfile = vm["Particlesfile"].as< std::string >();
166
167
168
169
170
171
172
		}
	}
	catch(std::exception& e)
	{
		cout << e.what() << "\n";
		return 1;
	}
Pilar Cossio's avatar
Pilar Cossio committed
173
174
175
176
177
        //check for consitency in multiple MRCs
        if(  RefMap.readMultMRC && not(RefMap.readMRC) ){
         cout << "For Multiple MRCs command --ReadMRC is necesary too";
         exit(1);
        }
David Rohr's avatar
David Rohr committed
178
	// ********************* Reading Parameter Input ***************************
179
180
181
182
	// copying inputfile to param class
	param.fileinput = infile.c_str();
	param.readParameters();

David Rohr's avatar
David Rohr committed
183
	// ********************* Reading Model Input ******************************
184
185
186
187
	// copying modelfile to model class
	Model.filemodel = modelfile.c_str();
	Model.readModel();

David Rohr's avatar
David Rohr committed
188
189
	// ********************* Reading Particle Maps Input **********************
	// ********* HERE: PROBLEM if maps dont fit on the memory!! ***************
190
	// copying mapfile to ref map class
191
	param.filemap = mapfile.c_str();
192
193
	RefMap.readRefMaps(param);

David Rohr's avatar
David Rohr committed
194
	// ****************** Precalculating Necessary Stuff *********************
195
	precalculate();
David Rohr's avatar
David Rohr committed
196

197
198
199
200
201
	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
202

203
	pProb.init(RefMap.ntotRefMap, param.nTotGridAngles, *this);
204

205
206
	deviceInit();

207
	return(0);
208
209
}

210
211
212
void bioem::cleanup()
{
	//Deleting allocated pointers
213
	free_device_host(pProb.ptr);
214
215
216
	RefMap.freePointers();
}

217
218
int bioem::precalculate()
{
David Rohr's avatar
David Rohr committed
219
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
220
	// **Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels**
David Rohr's avatar
David Rohr committed
221
	// **************************************************************************************
222

223
224
	// Generating Grids of orientations
	param.CalculateGridsParam();
225

226
227
	// Precalculating CTF Kernels stored in class Param
	param.CalculateRefCTF();
228

229
230
	//Precalculate Maps
	RefMap.precalculate(param, *this);
231

232
	return(0);
233
234
235
236
}

int bioem::run()
{
David Rohr's avatar
David Rohr committed
237
238
239
	// **************************************************************************************
	// **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
	// **************************************************************************************
240

David Rohr's avatar
David Rohr committed
241
242
	// **** If we want to control the number of threads -> omp_set_num_threads(XX); ******
	// ****************** Declarying class of Probability Pointer  *************************
243

244
	printf("\tInitializing Probabilities\n");
245
246
247
	// Inizialzing Probabilites to zero and constant to -Infinity
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
248
249
250
251
252
		bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

		pProbMap.Total = 0.0;
		pProbMap.Constoadd = -9999999;
		pProbMap.max_prob = -9999999;
253
		for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
254
		{
255
256
257
258
			bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);

			pProbAngle.forAngles = 0.0;
			pProbAngle.ConstAngle = -99999999;
259
260
		}
	}
David Rohr's avatar
David Rohr committed
261
	// **************************************************************************************
262
	deviceStartRun();
263
264
265
266
267
268
269
270
271
272
	{
		const int count = omp_get_max_threads();
		localCCT = new mycomplex_t*[count];
		lCC = new myfloat_t*[count];
		for (int i = 0;i < count;i++)
		{
			localCCT[i] = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
			lCC[i] = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
		}
	}
273

David Rohr's avatar
David Rohr committed
274
	// ******************************** MAIN CYCLE ******************************************
David Rohr's avatar
David Rohr committed
275

David Rohr's avatar
David Rohr committed
276
	// *** Declaring Private variables for each thread *****
277
	mycomplex_t* proj_mapFFT;
278
	myfloat_t* conv_map = new myfloat_t[param.param_device.NumberPixels * param.param_device.NumberPixels];
279
	mycomplex_t* conv_mapFFT;
280
	myfloat_t sumCONV, sumsquareCONV;
281
282

	//allocating fftw_complex vector
283
284
	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);
285
286
287
288

	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);
Pilar Cossio's avatar
Pilar Cossio committed
289
//	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)));
290
	for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
291
	{
David Rohr's avatar
David Rohr committed
292
293
		// ***************************************************************************************
		// ***** Creating Projection for given orientation and transforming to Fourier space *****
294
		timer.ResetStart();
295
296
		createProjection(iOrient, proj_mapFFT);
		printf("Time Projection %d: %f\n", iOrient, timer.GetCurrentElapsedTime());
297

David Rohr's avatar
David Rohr committed
298
299
		// ***************************************************************************************
		// ***** **** Internal Loop over convolutions **** *****
300
301
		for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
		{
302
			printf("\t\tConvolution %d %d\n", iOrient, iConv);
David Rohr's avatar
David Rohr committed
303
			// *** Calculating convolutions of projection map and crosscorrelations ***
304

305
			timer.ResetStart();
306
307
			createConvolutedProjectionMap(iOrient, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
			printf("Time Convolution %d %d: %f\n", iOrient, iConv, timer.GetCurrentElapsedTime());
308

David Rohr's avatar
David Rohr committed
309
310
			// ***************************************************************************************
			// *** Comparing each calculated convoluted map with all experimental maps ***
311
			timer.ResetStart();
312
			compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
313

314
315
316
			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 *
317
								  (((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;
318
			const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
319
								(((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;
320
321
			const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

322
			printf("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iOrient, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
323
324
325
326
327
		}
	}
	//deallocating fftw_complex vector
	myfftw_free(proj_mapFFT);
	myfftw_free(conv_mapFFT);
328
	delete[] conv_map;
David Rohr's avatar
David Rohr committed
329

330
	deviceFinishRun();
331
332
333
334
335
336
337
338
339
340
	{
		const int count = omp_get_max_threads();
		for (int i = 0;i < count;i++)
		{
			myfftw_free(localCCT[i]);
			myfftw_free(lCC[i]);
		}
		delete[] localCCT;
		delete[] lCC;
	}
341

David Rohr's avatar
David Rohr committed
342
	// ************* Writing Out Probabilities ***************
343

David Rohr's avatar
David Rohr committed
344
	// *** Angular Probability ***
345

346
347
348
349
	// if(param.writeAngles){
	ofstream angProbfile;
	angProbfile.open ("ANG_PROB_iRefMap");
	// }
350

351
352
353
354
	ofstream outputProbFile;
	outputProbFile.open ("Output_Probabilities");
	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
	{
David Rohr's avatar
David Rohr committed
355
		// **** Total Probability ***
356
357
358
		bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);

		outputProbFile << "RefMap " << iRefMap << " Probability  "  << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProbMap.Constoadd  << "\n";
359
360
361

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

David Rohr's avatar
David Rohr committed
362
		// *** Param that maximize probability****
363
364
365
366
367
368
369
370
371
		outputProbFile << (pProbMap.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[pProbMap.max_prob_orient].pos[0] << " ";
		outputProbFile << param.angles[pProbMap.max_prob_orient].pos[1] << " ";
		outputProbFile << param.angles[pProbMap.max_prob_orient].pos[2] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[0] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[1] << " ";
		outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[2] << " ";
		outputProbFile << pProbMap.max_prob_cent_x << " ";
		outputProbFile << pProbMap.max_prob_cent_y;
372
		outputProbFile << "\n";
373

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

376
		if(param.writeAngles)
377
		{
378
			for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
379
			{
380
				bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
381

382
				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";
383
384
385
			}
		}
	}
386

387
388
	angProbfile.close();
	outputProbFile.close();
389

390
	return(0);
391
392
}

393
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
394
{
David Rohr's avatar
David Rohr committed
395
396
	//***************************************************************************************
	//***** BioEM routine for comparing reference maps to convoluted maps *****
397
	if (FFTAlgo)
398
	{
David Rohr's avatar
David Rohr committed
399
		//With FFT Algorithm
400
401
		#pragma omp parallel for
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
402
		{
403
			const int num = omp_get_thread_num();
404
			calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, localCCT[num], lCC[num]);
405
406
407
		}
	}
	else
408
	{
David Rohr's avatar
David Rohr committed
409
		//Without FFT Algorithm
410
		#pragma omp parallel for
411
		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
412
		{
413
			compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap);
414
415
416
417
418
		}
	}
	return(0);
}

419
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)
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.FFTMapSize];
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
	// **************************************************************************************
David Rohr's avatar
David Rohr committed
439
440
	// ****  BioEM Create Projection routine in Euler angle predefined grid******************
	// ********************* and turns projection into Fourier space ************************
David Rohr's avatar
David Rohr committed
441
	// **************************************************************************************
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
	localproj = lCC[omp_get_thread_num()];
449
	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
	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++)
469
	{
470
471
472
		RotatedPointsModel[n].pos[0] = 0.0;
		RotatedPointsModel[n].pos[1] = 0.0;
		RotatedPointsModel[n].pos[2] = 0.0;
473
	}
474
	for(int n = 0; n < Model.nPointsModel; n++)
475
	{
476
		for(int k = 0; k < 3; k++)
477
		{
478
			for(int j = 0; j < 3; j++)
479
			{
480
				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.PointsModel[n].pos[j];
481
482
483
484
485
486
			}
		}
	}

	int i, j;

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

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

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

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

	return(0);
}

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

530
	mycomplex_t* tmp = localCCT[omp_get_thread_num()];
531

David Rohr's avatar
David Rohr committed
532
	// **** Multiplying FFTmap with corresponding kernel ****
533
	const mycomplex_t* refCTF = &param.refCTF[iConv * param.FFTMapSize];
534
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
535
	{
536
537
538
		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";
539
540
	}

541
542
543
	//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
544
	// **** Bringing convoluted Map to real Space ****
David Rohr's avatar
David Rohr committed
545
	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, Mapconv);
546

David Rohr's avatar
David Rohr committed
547
	// *** Calculating Cross-correlations of cal-convoluted map with its self *****
548
549
	sumC = 0;
	sumsquareC = 0;
550
	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
551
	{
David Rohr's avatar
David Rohr committed
552
553
		sumC += Mapconv[i];
		sumsquareC += Mapconv[i] * Mapconv[i];
554
	}
David Rohr's avatar
David Rohr committed
555
	// *** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
556
	// Normalizing
557
558
559
560
	myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
	myfloat_t norm4 = norm2 * norm2;
	sumC = sumC / norm2;
	sumsquareC = sumsquareC / norm4;
561
562

	return(0);
563
564
}

565
int bioem::calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare)
566
{
David Rohr's avatar
David Rohr committed
567
	// *********************** Routine to calculate Cross correlations***********************
568

569
570
	sum = 0.0;
	sumsquare = 0.0;
571
572
573
574
575
	for (int i = 0; i < param.param_device.NumberPixels; i++)
	{
		for (int j = 0; j < param.param_device.NumberPixels; j++)
		{
			// Calculate Sum of pixels
576
			sum += localmap[i * param.param_device.NumberPixels + j];
577
			// Calculate Sum of pixels squared
578
			sumsquare += localmap[i * param.param_device.NumberPixels + j] * localmap[i * param.param_device.NumberPixels + j];
579
580
581
		}
	}
	return(0);
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
}

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

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

int bioem::deviceFinishRun()
{
	return(0);
}
598
599
600
601
602
603
604
605
606
607

void* bioem::malloc_device_host(size_t size)
{
	return(mallocchk(size));
}

void bioem::free_device_host(void* ptr)
{
	free(ptr);
}