param.cpp 16 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
			if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);}
			cout << "Grid spacing Envelope " << gridEnvelop << "\n";
			yesGSPEnv = true ;
142
		}
David Rohr's avatar
David Rohr committed
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
		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
244
245
	// 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);
	}



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
}