map.cpp 4.57 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
14

using namespace std;

15
int bioem_RefMap::readRefMaps(bioem_param& param)
16
{
17
18
19
	/**************************************************************************************/
	/***********************Reading reference Particle Maps************************/
	/**************************************************************************************/
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
	if (loadMap)
	{
		FILE* fp = fopen("maps.dump", "rb");
		if (fp == NULL)
		{
			cout << "Error opening dump file\n";
			exit(1);
		}
		fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
		if (ntotRefMap > BIOEM_MAX_MAPS)
		{
			cout << "BIOEM_MAX_MAPS too small\n";
			exit(1);
		}
		fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
		fclose(fp);
36

37
		cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
38
		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
39
40
41
42
	}
	else
	{
		int nummap=-1;
43
44
		int lasti=0;
		int lastj=0;
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
		ifstream input(filemap);
		if (!input.good())
		{
			cout << "Particle Maps Failed to open file" << endl ;
			exit(1);
		}

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

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

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

			if (strcmp(token,"PARTICLE")==0) // to count the number of maps
			{
				nummap++;
				if (nummap % 128 == 0)
				{
					cout << "..." << nummap << "\n";
				}
				if (nummap == BIOEM_MAX_MAPS)
				{
					cout << "BIOEM_MAX_MAPS too small\n";
					exit(1);
				}
74
75
76
77
78
				if(lasti!=param.param_device.NumberPixels && lastj!=param.param_device.NumberPixels && nummap>0)
				{
					cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
					exit(1);
				}
79
80
81
82
83
			}
			else
			{
				int i,j;
				myfloat_t z;
84

85
86
87
88
89
90
91
92
93
94
95
96
97
98
				char tmpVals[36]  = {' '};

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

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

				strncpy (tmpVals,line+16,12);
				sscanf (tmpVals,"%f",&z);
				//checking for Map limits
				if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
				{
					Ref[nummap].points[i-1][j-1]=z;
99
100
					lasti=i;
					lastj=j;
101
102
103
104
105
106
107
108
109
110
111
112
				}
				else
				{
					cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
					exit(1);
				}
			}
		}
		cout << ".";
		ntotRefMap=nummap+1;
		cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
113

114
115
116
117
118
119
120
121
122
123
124
125
126
127
		if (dumpMap)
		{
			FILE* fp = fopen("maps.dump", "w+b");
			if (fp == NULL)
			{
				cout << "Error opening dump file\n";
				exit(1);
			}
			fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
			fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
			fclose(fp);
		}
	}

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
	return(0);
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
	/**************************************************************************************/
	/********** Routine that pre-calculates Kernels for Convolution **********************/
	/************************************************************************************/

	mycomplex_t* localMap;

	localMap= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
	RefMapFFT = new bioem_map_forFFT[ntotRefMap];
	mycomplex_t* localout;
	localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);

	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
		//Assigning localMap values to padded Map
		for(int i=0; i< param.param_device.NumberPixels; i++)
		{
			for(int j=0; j< param.param_device.NumberPixels; j++)
			{
				localMap[i*param.param_device.NumberPixels+j][0]=Ref[iRefMap].points[i][j];
				localMap[i*param.param_device.NumberPixels+j][1]=0.0;
			}
		}

		//Calling FFT_Forward
157
		myfftw_execute_dft(param.fft_plan_c2c_forward,localMap,localout);
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

		// Normalizing and saving the Reference CTFs
		for(int i=0; i < param.param_device.NumberPixels ; i++ )
		{
			for(int j=0; j < param.param_device.NumberPixels ; j++ )
			{
				RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]=localout[i*param.param_device.NumberPixels+j][0];
				RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1]=localout[i*param.param_device.NumberPixels+j][1];
			}
		}


	}

	myfftw_free(localMap);
	myfftw_free(localout);


	return(0);
177
}