param.cpp 16.2 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
			if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);}
115
			cout << "Grid points in Cosine ( beta ) " << angleGridPointsBeta << "\n";
David Rohr's avatar
David Rohr committed
116
			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
	// 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);};


227
228
229
230
231
232
        //if only one grid point for PSF kernel:
        if( (myfloat_t) numberGridPointsCTF_amp == 1 ) gridCTF_amp = startGridCTF_amp;
        if( (myfloat_t) numberGridPointsCTF_phase == 1 ) gridCTF_phase = startGridCTF_phase;
        if( (myfloat_t) numberGridPointsEnvelop == 1 ) gridEnvelop = startGridEnvelop;

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

244
	if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) (numberGridPointsCTF_amp - 1) * gridCTF_amp + startGridCTF_amp > 1){
David Rohr's avatar
David Rohr committed
245
		cout << "Value should be between 0 and 1\n" ;
David Rohr's avatar
David Rohr committed
246
247
248
		exit(1);
	}

249
        
250
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
251
252
253
254
255
256

	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);
257
	Alignment = 64;
258

259
260
261
262
	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);
263

264
	if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
265
266
267
268
269
270
271
272
	{
		cout << "Error planing FFTs\n";
		exit(1);
	}

	myfftw_free(tmp_map);
	myfftw_free(tmp_map2);
	fft_plans_created = 1;
273
274
275
276
277
278
279
280
281
282
283

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

284
285
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";

286
	return(0);
qon's avatar
qon committed
287
288
}

289
290
291
292
293
294
void bioem_param::releaseFFTPlans()
{
	if (fft_plans_created)
	{
		myfftw_destroy_plan(fft_plan_c2c_forward);
		myfftw_destroy_plan(fft_plan_c2c_backward);
295
296
		myfftw_destroy_plan(fft_plan_r2c_forward);
		myfftw_destroy_plan(fft_plan_c2r_backward);
297
298
299
300
	}
	fft_plans_created = 0;
}

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

	//alpha and gamma are uniform in -PI,PI
310
	grid_alpha = 2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
311
312

	//cosine beta is uniform in -1,1
313
	cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
314
315

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

333
        
David Rohr's avatar
David Rohr committed
334
	// ********** Calculating normalized volumen element *********
335

336
	param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
337
338
						* 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 )
						/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop );
339

340
      //  cout << "VOLU " << param_device.volu  << " " << gridCTF_amp << "\n";
David Rohr's avatar
David Rohr committed
341
	// *** Number of total pixels***
342

343
344
	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);
345
346

	return(0);
qon's avatar
qon committed
347
348
349
350
351

}

int bioem_param::CalculateRefCTF()
{
David Rohr's avatar
David Rohr committed
352
353
354
	// **************************************************************************************
	// ********** Routine that pre-calculates Kernels for Convolution **********************
	// ************************************************************************************
355

356
	myfloat_t amp, env, phase, ctf, radsq;
357
	myfloat_t* localCTF;
David Rohr's avatar
David Rohr committed
358
	mycomplex_t* localout;
359
360
	int nctfmax = param_device.NumberPixels / 2;
	int n = 0;
361

362
	localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
David Rohr's avatar
David Rohr committed
363
	localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberFFTPixels1D);
364
365
366

	nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
	delete[] refCTF;
367
	refCTF = new mycomplex_t[nTotCTFs * FFTMapSize];
368
369
	delete[] CtfParam;
	CtfParam = new myfloat3_t[nTotCTFs];
370
371
372

	for (int iamp = 0; iamp <  numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
	{
373
		amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
374
375
376

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

379
			for (int ienv = 0; ienv <  numberGridPointsEnvelop ; ienv++)//Loop over envelope
380
			{
381
				env = (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
382

383
				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
384
385

				//Assigning CTF values
386
				for(int i = 0; i < nctfmax; i++)
387
				{
388
					for(int j = 0; j < nctfmax; j++)
389
					{
390
391
						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));
392

393
394
395
396
						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;
397
398
399
					}
				}
				//Calling FFT_Forward
400
				myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
401
402

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

	myfftw_free(localCTF);
David Rohr's avatar
David Rohr committed
418
	myfftw_free(localout);
419
420
421
422
423
	if (nTotCTFs != n)
	{
		cout << "Internal error during CTF preparation\n";
		exit(1);
	}
424
425

	return(0);
qon's avatar
qon committed
426
427
428
429
}

bioem_param::~bioem_param()
{
430
	releaseFFTPlans();
431
	param_device.NumberPixels = 0;
432
433
434
435
436
437
438
	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
439
440
441
	if (refCTF) delete[] refCTF;
	if (CtfParam) delete[] CtfParam;
	if (angles) delete[] angles;
qon's avatar
qon committed
442
}