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);
}