param.cpp 11.5 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
217
218
219
220
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
221
222
int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
{
David Rohr's avatar
David Rohr committed
223
224
225
	// **************************************************************************************
	// **************** Routine that pre-calculates Euler angle grids **********************
	// ************************************************************************************
226
227
	myfloat_t grid_alpha, cos_grid_beta;
	int n = 0;
228
229

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

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

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

David Rohr's avatar
David Rohr committed
251
	// ********** Calculating normalized volumen element *********
252

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

David Rohr's avatar
David Rohr committed
257
	// *** Number of total pixels***
258

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

	return(0);
qon's avatar
qon committed
263
264
265
266
267

}

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

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

277
	localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
278
279
280

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

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

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

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

297
				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
298
299

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

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

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

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

	return(0);
qon's avatar
qon committed
343
344
345
346
347
}


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