param.cpp 15.8 KB
Newer Older
qon's avatar
qon committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <cstring>
#include <math.h>
#include <fftw3.h>

#include "param.h"
#include "map.h"

using namespace std;

bioem_param::bioem_param()
{
16
	//Number of Pixels
17
	param_device.NumberPixels = 0;
18
	param_device.NumberFFTPixels1D = 0;
19 20 21 22 23 24 25 26 27
	// Euler angle grid spacing
	angleGridPointsAlpha = 0;
	angleGridPointsBeta = 0;
	//Envelop function paramters
	numberGridPointsEnvelop = 0;
	//Contrast transfer function paramters
	numberGridPointsCTF_amp = 0;
	numberGridPointsCTF_phase = 0;

David Rohr's avatar
David Rohr committed
28
	// ****center displacement paramters Equal in both directions***
29 30
	param_device.maxDisplaceCenter = 0;
	numberGridPointsDisplaceCenter = 0;
31 32

	fft_plans_created = 0;
33 34 35

	refCTF = NULL;
	CtfParam = NULL;
David Rohr's avatar
David Rohr committed
36
	angles = NULL;
qon's avatar
qon committed
37 38 39 40
}

int bioem_param::readParameters()
{
David Rohr's avatar
David Rohr committed
41 42 43
	// **************************************************************************************
	// ***************************** Reading Input Parameters ******************************
	// **************************************************************************************
44

David Rohr's avatar
David Rohr committed
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
	// Control for Parameters
	bool yesPixSi = false;
	bool yesNumPix = false;
	bool yesGPal = false;
	bool yesGPbe = false;
	bool yesGPEnv = false;
	bool yesGPamp = false;
	bool yesGPpha = false;
	bool yesSTEnv = false;
	bool yesSTamp = false;
	bool yesSTpha = false;
	bool yesGSPamp = false ;
	bool yesGSPEnv = false ;
	bool yesGSPpha = false ;
	bool yesMDC = false ;
	bool yesGCen = false ;	
61

62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
	ifstream input(fileinput);
	if (!input.good())
	{
		cout << "Failed to open file: " << fileinput << "\n";
		exit(0);
	}

	char line[512] = {0};
	char saveline[512];

	cout << "\n +++++++++++++++++++++++++++++++++++++++++ \n";
	cout << "\n		   READING PARAMETERS  \n\n";

	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
	while (!input.eof())
	{
78 79 80
		input.getline(line, 512);
		strcpy(saveline, line);
		char *token = strtok(line, " ");
81

82
		if (token == NULL || line[0] == '#' || strlen(token) == 0)
83 84 85
		{
			// comment or blank line
		}
86
		else if (strcmp(token, "PIXEL_SIZE") == 0)
87
		{
88 89
			token = strtok(NULL, " ");
			pixelSize = atof(token);
David Rohr's avatar
David Rohr committed
90 91 92
			if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);}
			cout << "Pixel Sixe " << pixelSize << "\n";
			yesPixSi= true;
93
		}
94
		else if (strcmp(token, "NUMBER_PIXELS") == 0)
95
		{
96 97
			token = strtok(NULL, " ");
			param_device.NumberPixels = int(atoi(token));
David Rohr's avatar
David Rohr committed
98 99 100
			if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);}
			cout << "Number of Pixels " << param_device.NumberPixels << "\n";			
			yesNumPix= true ;
101
		}
102
		else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
103
		{
104 105
			token = strtok(NULL, " ");
			angleGridPointsAlpha = int(atoi(token));
David Rohr's avatar
David Rohr committed
106 107 108
			if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);}
			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
			yesGPal= true;
109
		}
110
		else if (strcmp(token, "GRIDPOINTS_BETA") == 0)
111
		{
112 113
			token = strtok(NULL, " ");
			angleGridPointsBeta = int(atoi(token));
David Rohr's avatar
David Rohr committed
114 115 116
			if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);}
			cout << "Grid points beta " << angleGridPointsBeta << "\n";
			yesGPbe= true;
117 118
		}

119
		else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
120
		{
121 122
			token = strtok(NULL, " ");
			numberGridPointsEnvelop = int(atoi(token));
David Rohr's avatar
David Rohr committed
123 124 125
			if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);}
			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
			yesGPEnv = true;
126
		}
127
		else if (strcmp(token, "START_ENVELOPE") == 0)
128
		{
129 130
			token = strtok(NULL, " ");
			startGridEnvelop = atof(token);
David Rohr's avatar
David Rohr committed
131 132 133
			if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START_ENVELOPE "; exit(1);}
			cout << "Start Envelope " << startGridEnvelop << "\n";
			yesSTEnv = true ;
134
		}
135
		else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0)
136
		{
137 138
			token = strtok(NULL, " ");
			gridEnvelop = atof(token);
David Rohr's avatar
David Rohr committed
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 183 184
			if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);}
			cout << "Grid spacing Envelope " << gridEnvelop << "\n";
			yesGSPEnv = true ;
		}
		else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0)
		{
			token = strtok(NULL," ");
			numberGridPointsCTF_phase=int(atoi(token));
			cout << "Grid points PSF " << numberGridPointsCTF_phase << "\n";
			yesGPpha = true;
		}
		else if (strcmp(token,"START_PSF_PHASE")==0)
		{
			token = strtok(NULL," ");
			startGridCTF_phase=atof(token);
			cout << "Start PSF " << startGridCTF_phase << "\n";
			yesSTpha = true ;
		}
		else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0)
		{
			token = strtok(NULL," ");
			gridCTF_phase=atof(token);
			cout << "Grid Space PSF " << gridCTF_phase << "\n";
			yesGSPpha = true ;
		}
		else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0)
		{
			token = strtok(NULL," ");
			numberGridPointsCTF_amp=int(atoi(token));
			cout << "Grid points PSF " << numberGridPointsCTF_amp << "\n";
			yesGPamp = true ;
		}
		else if (strcmp(token,"START_PSF_AMP")==0)
		{
			token = strtok(NULL," ");
			startGridCTF_amp=atof(token);
			cout << "Start PSF " << startGridCTF_amp << "\n";
			yesSTamp = true ;
		}
		else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0)
		{
			token = strtok(NULL," ");
			gridCTF_amp=atof(token);
			cout << "Grid Space PSF " << gridCTF_amp << "\n";
			yesGSPamp = true ;
		}
185
		else if (strcmp(token, "MAX_D_CENTER") == 0)
186
		{
187 188
			token = strtok(NULL, " ");
			param_device.maxDisplaceCenter = int(atoi(token));
David Rohr's avatar
David Rohr committed
189 190 191
			if (param_device.maxDisplaceCenter < 0 ) { cout << "*** Error: Negative MAX_D_CENTER "; exit(1);}
			cout << "Maximum displacement Center " <<  param_device.maxDisplaceCenter << "\n";
			yesMDC = true;
192
		}
193
		else if (strcmp(token, "PIXEL_GRID_CENTER") == 0)
194
		{
195 196
			token = strtok(NULL, " ");
			param_device.GridSpaceCenter = int(atoi(token));
David Rohr's avatar
David Rohr committed
197 198 199
			if (param_device.GridSpaceCenter < 0 ) { cout << "*** Error: Negative PIXEL_GRID_CENTER "; exit(1);}
			cout << "Grid space displacement center " <<   param_device.GridSpaceCenter << "\n";
			yesGCen = true;
200
		}
201
		else if (strcmp(token, "WRITE_PROB_ANGLES") == 0)
202
		{
203
			writeAngles = true;
204 205 206 207
			cout << "Writing Probabilies of each angle \n";
		}
	}
	input.close();
208

David Rohr's avatar
David Rohr committed
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
	// Checks for ALL INPUT 
	if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);};
	if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);};
	if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);};
	if( not ( yesGPbe ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);};
	if( not (  yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);};
	if( not (  yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);};
	if( not (  yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);};
	if( not (  yesSTEnv  ) ) { cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ; exit (1);};
	if( not (  yesSTamp  ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ; exit (1);};
	if( not (  yesSTpha  ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ; exit (1);};
	if( not (  yesGSPamp  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ; exit (1);};
	if( not ( yesGSPEnv  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ; exit (1);};
	if( not ( yesGSPpha  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ; exit (1);};
	if( not (  yesMDC  ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);};
	if( not (  yesGCen  ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);};

	//More checks with input parameters

	// Envelope should not have a standard deviation greater than Npix/2 
	if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter  * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) {
		cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length " ;
		exit(1);
	}
	// Envelop param should be positive
	if( startGridCTF_amp < 0 || startGridCTF_amp > 1){
		cout << "PSF Amplitud should be between 0 and 1  \n" ;	   
		exit(1);
	}

	if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1){
		cout << "PSF Amplitud should be between 0 and 1  \n" ;	   
		exit(1);
	}

David Rohr's avatar
David Rohr committed
244 245


246
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
247 248 249 250 251 252

	cout << "Preparing FFTs\n";
	releaseFFTPlans();
	mycomplex_t *tmp_map, *tmp_map2;
	tmp_map = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels);
	tmp_map2 = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels);
253
	Alignment = 64;
254

255 256 257 258
	fft_plan_c2c_forward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_FORWARD, FFTW_MEASURE | FFTW_DESTROY_INPUT);
	fft_plan_c2c_backward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_BACKWARD, FFTW_MEASURE | FFTW_DESTROY_INPUT);
	fft_plan_r2c_forward = myfftw_plan_dft_r2c_2d(param_device.NumberPixels, param_device.NumberPixels, (myfloat_t*) tmp_map, tmp_map2, FFTW_MEASURE | FFTW_DESTROY_INPUT);
	fft_plan_c2r_backward = myfftw_plan_dft_c2r_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, (myfloat_t*) tmp_map2, FFTW_MEASURE | FFTW_DESTROY_INPUT);
259

260
	if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
261 262 263 264 265 266 267 268
	{
		cout << "Error planing FFTs\n";
		exit(1);
	}

	myfftw_free(tmp_map);
	myfftw_free(tmp_map2);
	fft_plans_created = 1;
269 270 271 272 273 274 275 276 277 278 279

	param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
	FFTMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;

	//Does currently not work with custom alignment on GPU
	/*if (FFTMapSize % Alignment)
	{
		FFTMapSize += Alignment - FFTMapSize % Alignment;
	}
	cout << "Using MAP Size " << FFTMapSize << " (Alignment " << Alignment << ", Unaligned Size " << param_device.NumberPixels * param_device.NumberFFTPixels1D << ")\n";*/

280 281
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";

282
	return(0);
qon's avatar
qon committed
283 284
}

285 286 287 288 289 290
void bioem_param::releaseFFTPlans()
{
	if (fft_plans_created)
	{
		myfftw_destroy_plan(fft_plan_c2c_forward);
		myfftw_destroy_plan(fft_plan_c2c_backward);
291 292
		myfftw_destroy_plan(fft_plan_r2c_forward);
		myfftw_destroy_plan(fft_plan_c2r_backward);
293 294 295 296
	}
	fft_plans_created = 0;
}

qon's avatar
qon committed
297 298
int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
{
David Rohr's avatar
David Rohr committed
299 300 301
	// **************************************************************************************
	// **************** Routine that pre-calculates Euler angle grids **********************
	// ************************************************************************************
302 303
	myfloat_t grid_alpha, cos_grid_beta;
	int n = 0;
304 305

	//alpha and gamma are uniform in -PI,PI
306
	grid_alpha = 2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
307 308

	//cosine beta is uniform in -1,1
309
	cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
310 311

	// Euler Angle Array
312
	delete[] angles;
David Rohr's avatar
David Rohr committed
313
	angles = new myfloat3_t[angleGridPointsAlpha * angleGridPointsBeta * angleGridPointsAlpha];
314 315 316 317 318 319
	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
	{
		for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
		{
			for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
			{
320 321 322
				angles[n].pos[0] = (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
				angles[n].pos[1] = acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
				angles[n].pos[2] = (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
323 324 325 326
				n++;
			}
		}
	}
327
	nTotGridAngles = n;
328

David Rohr's avatar
David Rohr committed
329
	// ********** Calculating normalized volumen element *********
330

331
	param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
332 333
						* gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
						/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
334

David Rohr's avatar
David Rohr committed
335
	// *** Number of total pixels***
336

337 338
	param_device.Ntotpi = (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
	param_device.NtotDist = (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1);
339 340

	return(0);
qon's avatar
qon committed
341 342 343 344 345

}

int bioem_param::CalculateRefCTF()
{
David Rohr's avatar
David Rohr committed
346 347 348
	// **************************************************************************************
	// ********** Routine that pre-calculates Kernels for Convolution **********************
	// ************************************************************************************
349

350
	myfloat_t amp, env, phase, ctf, radsq;
351
	myfloat_t* localCTF;
David Rohr's avatar
David Rohr committed
352
	mycomplex_t* localout;
353 354
	int nctfmax = param_device.NumberPixels / 2;
	int n = 0;
355

356
	localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
David Rohr's avatar
David Rohr committed
357
	localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberFFTPixels1D);
358 359 360

	nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
	delete[] refCTF;
361
	refCTF = new mycomplex_t[nTotCTFs * FFTMapSize];
362 363
	delete[] CtfParam;
	CtfParam = new myfloat3_t[nTotCTFs];
364 365 366

	for (int iamp = 0; iamp <  numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
	{
367
		amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
368 369 370

		for (int iphase = 0; iphase <  numberGridPointsCTF_phase ; iphase++)//Loop over phase
		{
371
			phase = (myfloat_t) iphase * gridCTF_phase + startGridCTF_phase;
372

373
			for (int ienv = 0; ienv <  numberGridPointsEnvelop ; ienv++)//Loop over envelope
374
			{
375
				env = (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
376

377
				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
378 379

				//Assigning CTF values
380
				for(int i = 0; i < nctfmax; i++)
381
				{
382
					for(int j = 0; j < nctfmax; j++)
383
					{
384 385
						radsq = (myfloat_t) (i * i + j * j) * pixelSize * pixelSize;
						ctf = exp(-radsq * env / 2.0) * (amp * cos(radsq * phase / 2.0) - sqrt((1 - amp * amp)) * sin(radsq * phase / 2.0));
386

387 388 389 390
						localCTF[i * param_device.NumberPixels + j] = (myfloat_t) ctf;
						localCTF[i * param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
						localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + j] = (myfloat_t) ctf;
						localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
391 392 393
					}
				}
				//Calling FFT_Forward
394
				myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
395 396

				// Normalizing and saving the Reference CTFs
397
				mycomplex_t* curRef = &refCTF[n * FFTMapSize];
398
				for(int i = 0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
399
				{
400 401
					curRef[i][0] = localout[i][0];
					curRef[i][1] = localout[i][1];
402
				}
403 404 405
				CtfParam[n].pos[0] = amp;
				CtfParam[n].pos[1] = phase;
				CtfParam[n].pos[2] = env;
406 407 408 409 410 411
				n++;
			}
		}
	}

	myfftw_free(localCTF);
David Rohr's avatar
David Rohr committed
412
	myfftw_free(localout);
413 414 415 416 417
	if (nTotCTFs != n)
	{
		cout << "Internal error during CTF preparation\n";
		exit(1);
	}
418 419

	return(0);
qon's avatar
qon committed
420 421 422 423
}

bioem_param::~bioem_param()
{
424
	releaseFFTPlans();
425
	param_device.NumberPixels = 0;
426 427 428 429 430 431 432
	angleGridPointsAlpha = 0;
	angleGridPointsBeta = 0;
	numberGridPointsEnvelop = 0;
	numberGridPointsCTF_amp = 0;
	numberGridPointsCTF_phase = 0;
	param_device.maxDisplaceCenter = 0;
	numberGridPointsDisplaceCenter = 0;
David Rohr's avatar
David Rohr committed
433 434 435
	if (refCTF) delete[] refCTF;
	if (CtfParam) delete[] CtfParam;
	if (angles) delete[] angles;
qon's avatar
qon committed
436
}