param.cpp 11.6 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;
qon's avatar
qon committed
36
37
38
39
}

int bioem_param::readParameters()
{
David Rohr's avatar
David Rohr committed
40
41
42
	// **************************************************************************************
	// ***************************** Reading Input Parameters ******************************
	// **************************************************************************************
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

	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())
	{
60
61
62
		input.getline(line, 512);
		strcpy(saveline, line);
		char *token = strtok(line, " ");
63

64
		if (token == NULL || line[0] == '#' || strlen(token) == 0)
65
66
67
		{
			// comment or blank line
		}
68
		else if (strcmp(token, "PIXEL_SIZE") == 0)
69
		{
70
71
			token = strtok(NULL, " ");
			pixelSize = atof(token);
72
73
74
			cout << "Pixel Sixe " << pixelSize << "\n";

		}
75
		else if (strcmp(token, "NUMBER_PIXELS") == 0)
76
		{
77
78
			token = strtok(NULL, " ");
			param_device.NumberPixels = int(atoi(token));
79
80
81
			cout << "Number of Pixels " << param_device.NumberPixels << "\n";

		}
82
		else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
83
		{
84
85
			token = strtok(NULL, " ");
			angleGridPointsAlpha = int(atoi(token));
86
87
88
			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";

		}
89
		else if (strcmp(token, "GRIDPOINTS_BETA") == 0)
90
		{
91
92
			token = strtok(NULL, " ");
			angleGridPointsBeta = int(atoi(token));
93
94
95
96
			cout << "Grid points beta " << angleGridPointsBeta << "\n";

		}

97
		else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
98
		{
99
100
			token = strtok(NULL, " ");
			numberGridPointsEnvelop = int(atoi(token));
101
102
103
			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";

		}
104
		else if (strcmp(token, "START_ENVELOPE") == 0)
105
		{
106
107
			token = strtok(NULL, " ");
			startGridEnvelop = atof(token);
108
109
110
			cout << "Start Envelope " << startGridEnvelop << "\n";

		}
111
		else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0)
112
		{
113
114
			token = strtok(NULL, " ");
			gridEnvelop = atof(token);
115
116
117
			cout << "Grid spacing Envelope " << gridEnvelop << "\n";

		}
118
		else if (strcmp(token, "GRIDPOINTS_CTF_PHASE") == 0)
119
		{
120
121
			token = strtok(NULL, " ");
			numberGridPointsCTF_phase = int(atoi(token));
122
123
124
			cout << "Grid points CTF " << numberGridPointsCTF_phase << "\n";

		}
125
		else if (strcmp(token, "START_CTF_PHASE") == 0)
126
		{
127
128
			token = strtok(NULL, " ");
			startGridCTF_phase = atof(token);
129
130
131
			cout << "Start CTF " << startGridCTF_phase << "\n";

		}
132
		else if (strcmp(token, "GRIDSPACE_CTF_PHASE") == 0)
133
		{
134
135
			token = strtok(NULL, " ");
			gridCTF_phase = atof(token);
136
137
			cout << "Grid Space CTF " << gridCTF_phase << "\n";

138
		} else if (strcmp(token, "GRIDPOINTS_CTF_AMP") == 0)
139
		{
140
141
			token = strtok(NULL, " ");
			numberGridPointsCTF_amp = int(atoi(token));
142
143
144
			cout << "Grid points CTF " << numberGridPointsCTF_amp << "\n";

		}
145
		else if (strcmp(token, "START_CTF_AMP") == 0)
146
		{
147
148
			token = strtok(NULL, " ");
			startGridCTF_amp = atof(token);
149
150
151
			cout << "Start CTF " << startGridCTF_amp << "\n";

		}
152
		else if (strcmp(token, "GRIDSPACE_CTF_AMP") == 0)
153
		{
154
155
			token = strtok(NULL, " ");
			gridCTF_amp = atof(token);
156
157
158
			cout << "Grid Space CTF " << gridCTF_amp << "\n";

		}
159
		else if (strcmp(token, "MAX_D_CENTER") == 0)
160
		{
161
162
			token = strtok(NULL, " ");
			param_device.maxDisplaceCenter = int(atoi(token));
163
164
165
			cout << "Maximum displacement Center " <<  param_device.maxDisplaceCenter << "\n";

		}
166
		else if (strcmp(token, "PIXEL_GRID_CENTER") == 0)
167
		{
168
169
			token = strtok(NULL, " ");
			param_device.GridSpaceCenter = int(atoi(token));
170
171
172
			cout << "Grid space displacement center " <<   param_device.GridSpaceCenter << "\n";

		}
173
		else if (strcmp(token, "WRITE_PROB_ANGLES") == 0)
174
		{
175
			writeAngles = true;
176
177
178
179
180
181
			cout << "Writing Probabilies of each angle \n";

		}

	}
	input.close();
182
	param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
183
	FFTMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;
184
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
185
186
187
188
189
190
191

	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);

192
193
194
195
	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);
196

197
	if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
198
199
200
201
202
203
204
205
206
207
	{
		cout << "Error planing FFTs\n";
		exit(1);
	}

	myfftw_free(tmp_map);
	myfftw_free(tmp_map2);
	fft_plans_created = 1;
	cout << " +++++++++++++++++++++++++++++++++++++++++ \n";

208
	return(0);
qon's avatar
qon committed
209
210
}

211
212
213
214
215
216
void bioem_param::releaseFFTPlans()
{
	if (fft_plans_created)
	{
		myfftw_destroy_plan(fft_plan_c2c_forward);
		myfftw_destroy_plan(fft_plan_c2c_backward);
217
218
		myfftw_destroy_plan(fft_plan_r2c_forward);
		myfftw_destroy_plan(fft_plan_c2r_backward);
219
220
221
222
	}
	fft_plans_created = 0;
}

qon's avatar
qon committed
223
224
int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
{
David Rohr's avatar
David Rohr committed
225
226
227
	// **************************************************************************************
	// **************** Routine that pre-calculates Euler angle grids **********************
	// ************************************************************************************
228
229
	myfloat_t grid_alpha, cos_grid_beta;
	int n = 0;
230
231

	//alpha and gamma are uniform in -PI,PI
232
	grid_alpha = 2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
233
234

	//cosine beta is uniform in -1,1
235
	cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
236
237
238
239
240
241
242
243

	// Euler Angle Array
	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
	{
		for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
		{
			for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
			{
244
245
246
				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
247
248
249
250
				n++;
			}
		}
	}
251
	nTotGridAngles = n;
252

David Rohr's avatar
David Rohr committed
253
	// ********** Calculating normalized volumen element *********
254

255
	param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
256
257
						* 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);
258

David Rohr's avatar
David Rohr committed
259
	// *** Number of total pixels***
260

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

	return(0);
qon's avatar
qon committed
265
266
267
268
269

}

int bioem_param::CalculateRefCTF()
{
David Rohr's avatar
David Rohr committed
270
271
272
	// **************************************************************************************
	// ********** Routine that pre-calculates Kernels for Convolution **********************
	// ************************************************************************************
273

274
	myfloat_t amp, env, phase, ctf, radsq;
275
	myfloat_t* localCTF;
276
277
	int nctfmax = param_device.NumberPixels / 2;
	int n = 0;
278

279
	localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
280
281
282

	nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
	delete[] refCTF;
283
	refCTF = new mycomplex_t[nTotCTFs * FFTMapSize];
284
285
	delete[] CtfParam;
	CtfParam = new myfloat3_t[nTotCTFs];
286
287
288

	for (int iamp = 0; iamp <  numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
	{
289
		amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
290
291
292

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

295
			for (int ienv = 0; ienv <  numberGridPointsEnvelop ; ienv++)//Loop over envelope
296
			{
297
				env = (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
298

299
				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
300
301

				//Assigning CTF values
302
				for(int i = 0; i < nctfmax; i++)
303
				{
304
					for(int j = 0; j < nctfmax; j++)
305
					{
306
307
						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));
308

309
310
311
312
						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;
313
314
315
316
					}
				}
				//Creatting FFT_Forward of Kernel to store
				mycomplex_t* localout;
317
				localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberFFTPixels1D);
318
				//Calling FFT_Forward
319
				myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
320
321

				// Normalizing and saving the Reference CTFs
322
				mycomplex_t* curRef = &refCTF[n * FFTMapSize];
323
				for(int i = 0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
324
				{
325
326
					curRef[i][0] = localout[i][0];
					curRef[i][1] = localout[i][1];
327
				}
328
329
330
				CtfParam[n].pos[0] = amp;
				CtfParam[n].pos[1] = phase;
				CtfParam[n].pos[2] = env;
331
332
333
334
335
336
337
				n++;
				myfftw_free(localout);
			}
		}
	}

	myfftw_free(localCTF);
338
339
340
341
342
	if (nTotCTFs != n)
	{
		cout << "Internal error during CTF preparation\n";
		exit(1);
	}
343
344

	return(0);
qon's avatar
qon committed
345
346
347
348
349
}


bioem_param::~bioem_param()
{
350
	releaseFFTPlans();
351
	param_device.NumberPixels = 0;
352
353
354
355
356
357
358
	angleGridPointsAlpha = 0;
	angleGridPointsBeta = 0;
	numberGridPointsEnvelop = 0;
	numberGridPointsCTF_amp = 0;
	numberGridPointsCTF_phase = 0;
	param_device.maxDisplaceCenter = 0;
	numberGridPointsDisplaceCenter = 0;
359
360
	delete[] refCTF;
	delete[] CtfParam;
qon's avatar
qon committed
361
}