param.cpp 15.3 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
17
	//Number of Pixels
	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;

Pilar Cossio's avatar
Pilar Cossio 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;
qon's avatar
qon committed
36
37
38
39
}

int bioem_param::readParameters()
{
40
41
42
	//**************************************************************************************
	//***************************** Reading Input Parameters ******************************
	//**************************************************************************************
43

Pilar Cossio's avatar
Pilar Cossio committed
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        // 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 ;	
	
	
	
63
64
65
66
67
68
69
70
71
72
	ifstream input(fileinput);
	if (!input.good())
	{
		cout << "Failed to open file: " << fileinput << "\n";
		exit(0);
	}

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

Pilar Cossio's avatar
Pilar Cossio committed
73
74
	cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
	cout << "\n	      READING PARAMETERS  \n\n";
75

Pilar Cossio's avatar
Pilar Cossio committed
76
	cout << "+++++++++++++++++++++++++++++++++++++++++++ \n";
77
78
79
80
81
82
83
84
85
86
87
88
89
90
	while (!input.eof())
	{
		input.getline(line,512);
		strcpy(saveline,line);
		char *token = strtok(line," ");

		if (token==NULL || line[0] == '#' || strlen(token)==0)
		{
			// comment or blank line
		}
		else if (strcmp(token,"PIXEL_SIZE")==0)
		{
			token = strtok(NULL," ");
			pixelSize=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
91
			if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);}
92
			cout << "Pixel Sixe " << pixelSize << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
93
                        yesPixSi= true;
94
95
96
97
98
99

		}
		else if (strcmp(token,"NUMBER_PIXELS")==0)
		{
			token = strtok(NULL," ");
			param_device.NumberPixels=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
100
101
102
			if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);}
			cout << "Number of Pixels " << param_device.NumberPixels << "\n";			
                        yesNumPix= true ;
103
104
105
106
107
		}
		else if (strcmp(token,"GRIDPOINTS_ALPHA")==0)
		{
			token = strtok(NULL," ");
			angleGridPointsAlpha=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
108
			if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);}
109
			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
110
                        yesGPal= true;
111
112
113
114
		}
		else if (strcmp(token,"GRIDPOINTS_BETA")==0)
		{
			token = strtok(NULL," ");
Pilar Cossio's avatar
Pilar Cossio committed
115
116
			angleGridPointsBeta=int(atoi(token));			
			if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);}
117
			cout << "Grid points beta " << angleGridPointsBeta << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
118
                        yesGPbe= true;
119
120
121
122
123
124
		}

		else if (strcmp(token,"GRIDPOINTS_ENVELOPE")==0)
		{
			token = strtok(NULL," ");
			numberGridPointsEnvelop=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
125
			if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);}
126
			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
127
			yesGPEnv = true;
128
129
130
131
132
		}
		else if (strcmp(token,"START_ENVELOPE")==0)
		{
			token = strtok(NULL," ");
			startGridEnvelop=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
133
			if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START_ENVELOPE "; exit(1);}
134
			cout << "Start Envelope " << startGridEnvelop << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
135
			yesSTEnv = true ;
136
137
138
139
140
141

		}
		else if (strcmp(token,"GRIDSPACE_ENVELOPE")==0)
		{
			token = strtok(NULL," ");
			gridEnvelop=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
142
			if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);}
143
			cout << "Grid spacing Envelope " << gridEnvelop << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
144
			yesGSPEnv = true ;
145
146

		}
Pilar Cossio's avatar
Pilar Cossio committed
147
		else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0)
148
149
150
		{
			token = strtok(NULL," ");
			numberGridPointsCTF_phase=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
151
152
			cout << "Grid points PSF " << numberGridPointsCTF_phase << "\n";
			yesGPpha = true;
153
154

		}
Pilar Cossio's avatar
Pilar Cossio committed
155
		else if (strcmp(token,"START_PSF_PHASE")==0)
156
157
158
		{
			token = strtok(NULL," ");
			startGridCTF_phase=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
159
160
			cout << "Start PSF " << startGridCTF_phase << "\n";
			yesSTpha = true ;
161
162

		}
Pilar Cossio's avatar
Pilar Cossio committed
163
		else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0)
164
165
166
		{
			token = strtok(NULL," ");
			gridCTF_phase=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
167
168
			cout << "Grid Space PSF " << gridCTF_phase << "\n";
			yesGSPpha = true ;
169

Pilar Cossio's avatar
Pilar Cossio committed
170
		} else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0)
171
172
173
		{
			token = strtok(NULL," ");
			numberGridPointsCTF_amp=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
174
175
			cout << "Grid points PSF " << numberGridPointsCTF_amp << "\n";
			yesGPamp = true ;
176
177

		}
Pilar Cossio's avatar
Pilar Cossio committed
178
		else if (strcmp(token,"START_PSF_AMP")==0)
179
180
181
		{
			token = strtok(NULL," ");
			startGridCTF_amp=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
182
183
			cout << "Start PSF " << startGridCTF_amp << "\n";
			yesSTamp = true ;
184
185

		}
Pilar Cossio's avatar
Pilar Cossio committed
186
		else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0)
187
188
189
		{
			token = strtok(NULL," ");
			gridCTF_amp=atof(token);
Pilar Cossio's avatar
Pilar Cossio committed
190
191
			cout << "Grid Space PSF " << gridCTF_amp << "\n";
			yesGSPamp = true ;
192
193
194
195
196
197

		}
		else if (strcmp(token,"MAX_D_CENTER")==0)
		{
			token = strtok(NULL," ");
			param_device.maxDisplaceCenter=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
198
			if (param_device.maxDisplaceCenter < 0 ) { cout << "*** Error: Negative MAX_D_CENTER "; exit(1);}
199
			cout << "Maximum displacement Center " <<  param_device.maxDisplaceCenter << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
200
			yesMDC = true;
201
202
203
204
205
206

		}
		else if (strcmp(token,"PIXEL_GRID_CENTER")==0)
		{
			token = strtok(NULL," ");
			param_device.GridSpaceCenter=int(atoi(token));
Pilar Cossio's avatar
Pilar Cossio committed
207
			if (param_device.GridSpaceCenter < 0 ) { cout << "*** Error: Negative PIXEL_GRID_CENTER "; exit(1);}
208
			cout << "Grid space displacement center " <<   param_device.GridSpaceCenter << "\n";
Pilar Cossio's avatar
Pilar Cossio committed
209
			yesGCen = true;
210
211
212
213
214
215
216
217
218
219
220

		}
		else if (strcmp(token,"WRITE_PROB_ANGLES")==0)
		{
			writeAngles=true;
			cout << "Writing Probabilies of each angle \n";

		}

	}
	input.close();
Pilar Cossio's avatar
Pilar Cossio committed
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
246
247
248
249
250
251
252
253
254
255
256
257
258
	
	// 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);
	}
	
	
	  
259
	param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
260
	RefMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;
Pilar Cossio's avatar
Pilar Cossio committed
261
	cout << "+++++++++++++++++++++++++++++++++++++++++++ \n";
262

Pilar Cossio's avatar
Pilar Cossio committed
263
	cout << "                Preparing FFTs\n";
264
265
266
267
268
	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);

269
270
271
272
	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);
273

274
	if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
275
276
277
278
279
280
281
282
	{
		cout << "Error planing FFTs\n";
		exit(1);
	}

	myfftw_free(tmp_map);
	myfftw_free(tmp_map2);
	fft_plans_created = 1;
Pilar Cossio's avatar
Pilar Cossio committed
283
	cout << "+++++++++++++++++++++++++++++++++++++++++++ \n";
284

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

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

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

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

	//cosine beta is uniform in -1,1
310
	cos_grid_beta=2.f / (myfloat_t) angleGridPointsBeta;
311
312
313
314
315
316
317
318

	// Euler Angle Array
	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
	{
		for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
		{
			for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
			{
319
320
321
				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
322
323
324
325
326
327
				n++;
			}
		}
	}
	nTotGridAngles=n;

Pilar Cossio's avatar
Pilar Cossio committed
328
	//********** Calculating normalized volumen element *********
329

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

Pilar Cossio's avatar
Pilar Cossio committed
334
	//*** Number of total pixels***
335

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

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

}

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

	myfloat_t amp,env,phase,ctf,radsq;
350
	myfloat_t* localCTF;
351
	int nctfmax= param_device.NumberPixels / 2;
352
353
	int n=0;

354
	localCTF= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels*param_device.NumberPixels);
355
356
357
358
359
360

	nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
	delete[] refCTF;
	refCTF = new mycomplex_t[nTotCTFs * RefMapSize];
	delete[] CtfParam;
	CtfParam = new myfloat3_t[nTotCTFs];
361
362
363

	for (int iamp = 0; iamp <  numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
	{
364
		amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
Pilar Cossio's avatar
Pilar Cossio committed
365
		
366
367
		for (int iphase = 0; iphase <  numberGridPointsCTF_phase ; iphase++)//Loop over phase
		{
368
			phase = (myfloat_t) iphase * gridCTF_phase + startGridCTF_phase;
369

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

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

				//Assigning CTF values
				for(int i=0; i< nctfmax; i++)
				{
					for(int j=0; j< nctfmax; j++)
					{
381
						radsq=(myfloat_t) (i*i+j*j) * pixelSize * pixelSize;
382
						ctf=exp(-radsq*env/2.0)*(amp*cos(radsq*phase/2.0)-sqrt((1-amp*amp))*sin(radsq*phase/2.0));
Pilar Cossio's avatar
Pilar Cossio committed
383
						
384
385
386
387
						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;
388
389
390
391
					}
				}
				//Creatting FFT_Forward of Kernel to store
				mycomplex_t* localout;
392
				localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param_device.NumberPixels*param_device.NumberFFTPixels1D);
393
				//Calling FFT_Forward
394
				myfftw_execute_dft_r2c(fft_plan_r2c_forward,localCTF,localout);
395
396

				// Normalizing and saving the Reference CTFs
397
398
				mycomplex_t* curRef = &refCTF[n * RefMapSize];
				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
406
407
408
409
410
411
412
				}
				CtfParam[n].pos[0]=amp;
				CtfParam[n].pos[1]=phase;
				CtfParam[n].pos[2]=env;
				n++;
				myfftw_free(localout);
			}
		}
	}

	myfftw_free(localCTF);
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
424
}


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