map.cpp 5.24 KB
Newer Older
1
2
3
4
5
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstring>
6
7
#include <math.h>
#include <fftw3.h>
8
9

#include "map.h"
10
11
#include "param.h"

12
13
using namespace std;

14
int bioem_RefMap::readRefMaps(bioem_param& param)
15
{
16
17
	numPixels = param.param_device.NumberPixels;
	refMapSize = param.param_device.NumberPixels * param.param_device.NumberPixels;
David Rohr's avatar
David Rohr committed
18
19
20
	// **************************************************************************************
	// ***********************Reading reference Particle Maps************************
	// **************************************************************************************
21
22
	int allocsize = 0;
	if (param.loadMap)
23
24
25
26
27
28
29
30
	{
		FILE* fp = fopen("maps.dump", "rb");
		if (fp == NULL)
		{
			cout << "Error opening dump file\n";
			exit(1);
		}
		fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
31
32
		maps = (myfloat_t*) mallocchk(ntotRefMap * refMapSize * sizeof(myfloat_t));
		fread(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
33
		fclose(fp);
34

35
		cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
36
		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
37
38
39
	}
	else
	{
40
41
42
		int nummap = -1;
		int lasti = 0;
		int lastj = 0;
43
		ifstream input(param.filemap);
44
45
46
47
48
49
50
51
52
53
54
		if (!input.good())
		{
			cout << "Particle Maps Failed to open file" << endl ;
			exit(1);
		}

		char line[512] = {' '};
		char tmpLine[512] = {' '};

		while (!input.eof())
		{
55
			input.getline(line, 512);
56

57
58
			strncpy(tmpLine, line, strlen(line));
			char *token = strtok(tmpLine, " ");
59

60
			if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
61
62
			{
				nummap++;
63
64
65
66
67
68
69
70
71
72
73
				if (allocsize == 0)
				{
					allocsize = 64;
					maps = (myfloat_t*) mallocchk(refMapSize * sizeof(myfloat_t) * allocsize);

				}
				else if (nummap + 1 >= allocsize)
				{
					allocsize *= 2;
					maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * allocsize);
				}
74
75
76
77
78
79
80
81
82
				if (nummap % 128 == 0)
				{
					cout << "..." << nummap << "\n";
				}
				if (nummap == BIOEM_MAX_MAPS)
				{
					cout << "BIOEM_MAX_MAPS too small\n";
					exit(1);
				}
83
				if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
84
85
86
87
				{
					cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
					exit(1);
				}
88
89
90
			}
			else
			{
91
				int i, j;
92
				float z;
93

94
95
				char tmpVals[36]  = {' '};

96
97
				strncpy (tmpVals, line, 8);
				sscanf (tmpVals, "%d", &i);
98

99
100
				strncpy (tmpVals, line + 8, 8);
				sscanf (tmpVals, "%d", &j);
101

102
103
				strncpy (tmpVals, line + 16, 12);
				sscanf (tmpVals, "%f", &z);
104
				//checking for Map limits
105
				if(i > 0 && i - 1 < BIOEM_MAP_SIZE_X && j > 0 && j - 1 < BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
106
				{
107
108
109
					maps[nummap * refMapSize + (i - 1) * numPixels + (j - 1)] = (myfloat_t) z;
					lasti = i;
					lastj = j;
110
111
112
113
114
115
116
117
118
				}
				else
				{
					cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
					exit(1);
				}
			}
		}
		cout << ".";
119
		ntotRefMap = nummap + 1;
120
		maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
121
122
		cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
123

124
		if (param.dumpMap)
125
126
127
128
129
130
131
132
		{
			FILE* fp = fopen("maps.dump", "w+b");
			if (fp == NULL)
			{
				cout << "Error opening dump file\n";
				exit(1);
			}
			fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
133
			fwrite(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
134
135
136
137
			fclose(fp);
		}
	}

138
139
140
	sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
	sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);

141
142
143
144
145
	return(0);
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
David Rohr's avatar
David Rohr committed
146
147
148
	// **************************************************************************************
	// ********** Routine that pre-calculates Kernels for Convolution **********************
	// ************************************************************************************
149

150
151
	myfloat_t* localMap;

152
	localMap = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
153

154
	RefMapsFFT = new mycomplex_t[ntotRefMap * param.FFTMapSize];
155

156
	mycomplex_t* localout;
157
	localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
158
159
160
161

	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
		//Assigning localMap values to padded Map
162
		for(int i = 0; i < param.param_device.NumberPixels; i++)
163
		{
164
			for(int j = 0; j < param.param_device.NumberPixels; j++)
165
			{
166
				localMap[i * param.param_device.NumberPixels + j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
167
168
169
170
			}
		}

		//Calling FFT_Forward
171
		myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
172
173

		// Normalizing and saving the Reference CTFs
174
		mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.FFTMapSize];
175
		for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
176
		{
177
178
			RefMap[i][0] = localout[i][0];
			RefMap[i][1] = localout[i][1];
179
180
181
182
183
184
185
		}
	}

	myfftw_free(localMap);
	myfftw_free(localout);

	return(0);
186
}