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
	// 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) {
David Rohr's avatar
David Rohr committed
230
		cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length";
David Rohr's avatar
David Rohr committed
231
232
233
234
		exit(1);
	}
	// Envelop param should be positive
	if( startGridCTF_amp < 0 || startGridCTF_amp > 1){
David Rohr's avatar
David Rohr committed
235
		cout << "PSF Amplitud should be between 0 and 1\n";
David Rohr's avatar
David Rohr committed
236
237
238
239
		exit(1);
	}

	if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1){
David Rohr's avatar
David Rohr committed
240
		cout << "Value should be between 0 and 1\n" ;
David Rohr's avatar
David Rohr committed
241
242
243
		exit(1);
	}

244
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
245
246
247
248
249
250

	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);
251
	Alignment = 64;
252

253
254
255
256
	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);
257

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

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

	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";*/

278
279
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";

280
	return(0);
qon's avatar
qon committed
281
282
}

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

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

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

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

	// Euler Angle Array
310
	delete[] angles;
David Rohr's avatar
David Rohr committed
311
	angles = new myfloat3_t[angleGridPointsAlpha * angleGridPointsBeta * angleGridPointsAlpha];
312
313
314
315
316
317
	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
	{
		for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
		{
			for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
			{
318
319
320
				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
321
322
323
324
				n++;
			}
		}
	}
325
	nTotGridAngles = n;
326

David Rohr's avatar
David Rohr committed
327
	// ********** Calculating normalized volumen element *********
328

329
	param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
330
331
						* 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);
332

David Rohr's avatar
David Rohr committed
333
	// *** Number of total pixels***
334

335
336
	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);
337
338

	return(0);
qon's avatar
qon committed
339
340
341
342
343

}

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

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

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

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

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

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

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

375
				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
376
377

				//Assigning CTF values
378
				for(int i = 0; i < nctfmax; i++)
379
				{
380
					for(int j = 0; j < nctfmax; j++)
381
					{
382
383
						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));
384

385
386
387
388
						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;
389
390
391
					}
				}
				//Calling FFT_Forward
392
				myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
393
394

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

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

	return(0);
qon's avatar
qon committed
418
419
420
421
}

bioem_param::~bioem_param()
{
422
	releaseFFTPlans();
423
	param_device.NumberPixels = 0;
424
425
426
427
428
429
430
	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
431
432
433
	if (refCTF) delete[] refCTF;
	if (CtfParam) delete[] CtfParam;
	if (angles) delete[] angles;
qon's avatar
qon committed
434
}