Openholo  v1.1
Open Source Digital Holographic Library
ophLUT.cpp
Go to the documentation of this file.
1 #define OPH_DM_EXPORT
2 
3 #include "ophLUT.h"
4 #include <fstream>
5 #include <iostream>
6 #include <string>
7 #include <iomanip>
8 
9 #include "sys.h"
10 
11 //CGHEnvironmentData CONF; // config
12 
13 using namespace std;
14 
16 {
17 }
18 
20 {
21  delete[] m_pHologram;
22 }
23 
24 int ophLUT::init(const char* _filename, CGHEnvironmentData* _CGHE)
25 {
26  cout << _filename << endl;
27 
28  loadConfig("ConfigInfo.ini", _CGHE);
29 
30  return 0;
31 }
32 
33 bool ophLUT::loadConfig(const char* filename, CGHEnvironmentData* conf)
34 {
35 #define MAX_SIZE 1000
36  char inputString[MAX_SIZE];
37 
38  ifstream inFile(filename);
39  if (!(inFile.is_open()))
40  {
41  cout << "파일을 찾을 수 없습니다." << endl;
42  return false;
43  }
44  else {
45  while (!inFile.eof())
46  {
47  inFile.getline(inputString, MAX_SIZE);
48 
49  // 주석 및 빈칸 제거 if (!(inputString[0] == NULL || (inputString[0] == '#' && inputString[1] == ' '))) { char* token = NULL; char* parameter = NULL; char* context = nullptr; token = strtok_s(inputString, "=", &context); // 데이터 받아서 구조체에 저장 if (strcmp(token, "CGH width ") == 0) { token = strtok_s(NULL, "=", &context); conf->CghWidth = atoi(trim(token)); } else if (strcmp(token, "CGH height ") == 0) { token = strtok_s(NULL, "=", &context); conf->CghHeight = atoi(trim(token)); } else if (strcmp(token, "Segmentation size ") == 0) { token = strtok_s(NULL, "=", &context); conf->SegmentationSize = atoi(trim(token)); } else if (strcmp(token, "FFT segmentation size ") == 0) { token = strtok_s(NULL, "=", &context); conf->fftSegmentationSize = atoi(trim(token)); } else if (strcmp(token, "Red wavelength ") == 0) { token = strtok_s(NULL, "=", &context); conf->rWaveLength = atof(trim(token)); } else if (strcmp(token, "Tilting angle on x axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->ThetaX = atof(trim(token)); } else if (strcmp(token, "Tilting angle on y axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->ThetaY = atof(trim(token)); } else if (strcmp(token, "Default depth ") == 0) { token = strtok_s(NULL, "=", &context); conf->DefaultDepth = atof(trim(token)); } else if (strcmp(token, "3D point interval on x axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->xInterval = atof(trim(token)); } else if (strcmp(token, "3D point interval on y axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->yInterval = atof(trim(token)); } else if (strcmp(token, "Hologram interval on xi axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->xiInterval = atof(trim(token)); } else if (strcmp(token, "Hologram interval on eta axis ") == 0) { token = strtok_s(NULL, "=", &context); conf->etaInterval = atof(trim(token)); } else if (strcmp(token, "CGH scale ") == 0) { token = strtok_s(NULL, "=", &context); conf->CGHScale = atof(trim(token)); } } } } inFile.close(); return false; } bool ophLUT::loadPoint(const char* _filename, VoxelStruct* h_vox) { #define MAX_SIZE 1000 char inputString[MAX_SIZE]; ifstream inFile(_filename); if (!(inFile.is_open())) { cout << "포인트 클라우드 파일을 찾을 수 없습니다." << endl; return false; } else { inFile.getline(inputString, MAX_SIZE); int no = 0; while (!inFile.eof()) { inFile.getline(inputString, MAX_SIZE); if (inputString[0] != NULL) { char* token = NULL; char* context; token = strtok_s(inputString, "\t", &context); h_vox[no].num = atoi(token); // 인덱스 token = strtok_s(NULL, "\t", &context); h_vox[no].x = atof(token); // x 좌표 token = strtok_s(NULL, "\t", &context); h_vox[no].y = atof(token); // y 좌표 token = strtok_s(NULL, "\t", &context); h_vox[no].z = atof(token); // z 좌표 token = strtok_s(NULL, "\t", &context); h_vox[no].ph = atof(token); // phase token = strtok_s(NULL, "\t", &context); h_vox[no].r = atof(token); // red //token = strtok_s(NULL, "\t"); //h_vox[no].g = atof(token); // green //token = strtok_s(NULL, "\t"); //h_vox[no].b = atof(token); // blue no++; } } } inFile.close(); return true; } bool ophLUT::load_Num_Point(const char* _filename, long* num_point) { #define MAX_SIZE 1000 char inputString[MAX_SIZE]; ifstream inFile(_filename); if (!(inFile.is_open())) { cout << "포인트 클라우드 파일을 찾을 수 없습니다." << endl; return false; } else { inFile.getline(inputString, MAX_SIZE); *num_point = atoi(trim(inputString)); //cout << *num_point << endl; } inFile.close(); return true; } void ophLUT::DataInit(CGHEnvironmentData* _CGHE) { m_pHologram = new double[_CGHE->CghHeight*_CGHE->CghWidth]; memset(m_pHologram, 0x00, sizeof(double)*_CGHE->CghHeight*_CGHE->CghWidth); for (int i = 0; i<NUMTBL; i++) { float theta = (float)M2_PI * (float)(i + i - 1) / (float)(2 * NUMTBL); m_COStbl[i] = (float)cos(theta); m_SINtbl[i] = (float)sin(theta); } } // APAS int ophLUT::LUTCalcuation(long voxnum, unsigned char *cghfringe, VoxelStruct* h_vox, CGHEnvironmentData* _CGHE) { long i, j; double Max = -1E9, Min = 1E9; double myBuffer; int cghwidth = _CGHE->CghWidth; int cghheight = _CGHE->CghHeight; DataInit(_CGHE); // APAS(voxnum, h_vox, _CGHE); // for (i = 0; i<cghheight; i++) { for (j = 0; j<cghwidth; j++) { if (Max < m_pHologram[i*cghwidth + j]) Max = m_pHologram[i*cghwidth + j]; if (Min > m_pHologram[i*cghwidth + j]) Min = m_pHologram[i*cghwidth + j]; } } for (i = 0; i<cghheight; i++) { for (j = 0; j<cghwidth; j++) { myBuffer = 1.0*(((m_pHologram[i*cghwidth + j] - Min) / (Max - Min))*255. + 0.5); if (myBuffer >= 255.0) cghfringe[i*cghwidth + j] = 255; else cghfringe[i*cghwidth + j] = (unsigned char)(myBuffer); } } return 0; } void ophLUT::APAS(long voxelnum, VoxelStruct * _h_vox, CGHEnvironmentData * _CGHE) { int no; // voxel Number int i, j; int segx, segy; // coordinate in a Segment float R; int cghwidth = _CGHE->CghWidth; int cghheight = _CGHE->CghHeight; float xiInterval = _CGHE->xiInterval; float etaInterval = _CGHE->etaInterval; float rLamda = _CGHE->rWaveLength; float rWaveNum = _CGHE->rWaveNumber; float thetaX = _CGHE->ThetaX; float thetaY = _CGHE->ThetaY; int segSize = _CGHE->SegmentationSize; int hsegSize = (int)(segSize / 2); int dsegSize = segSize*segSize; int segNumx = (int)(cghwidth / segSize); int segNumy = (int)(cghheight / segSize); int hsegNumx = (int)(segNumx / 2); int hsegNumy = (int)(segNumy / 2); int FFTsegSize = _CGHE->fftSegmentationSize; int FFThsegSize = (int)(FFTsegSize / 2); int FFTdsegSize = FFTsegSize*FFTsegSize; float X, Y, Z; ; // x, y, real distance float theta_cx, theta_cy; float *SFrequency_cx = new float[segNumx]; // Center phase float *SFrequency_cy = new float[segNumy]; // Center phase int *PickPoint_cx = new int[segNumx]; int *PickPoint_cy = new int[segNumy]; int *Coefficient_cx = new int[segNumx]; int *Coefficient_cy = new int[segNumy]; float *xc = new float[segNumx]; float *yc = new float[segNumy]; float Amplitude; float phase; float sf_base = 1.0 / (xiInterval*FFTsegSize); int segxx, segyy; float theta; float theta_s, theta_c; int dtheta_s, dtheta_c; int idx_c, idx_s; //CString mm; fftw_complex *in, *out; fftw_plan plan; double **inRe = new double *[segNumy * segNumx]; double **inIm = new double *[segNumy * segNumx]; for (i = 0; i<segNumy; i++) { for (j = 0; j<segNumx; j++) { inRe[i*segNumx + j] = new double[FFTsegSize * FFTsegSize]; inIm[i*segNumx + j] = new double[FFTsegSize * FFTsegSize]; memset(inRe[i*segNumx + j], 0x00, sizeof(double) * FFTsegSize * FFTsegSize); memset(inIm[i*segNumx + j], 0x00, sizeof(double) * FFTsegSize * FFTsegSize); } } in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTsegSize * FFTsegSize); out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTsegSize * FFTsegSize); memset(in, 0x00, sizeof(fftw_complex) * FFTsegSize * FFTsegSize); memset(m_pHologram, 0x00, sizeof(double)*cghwidth*cghheight); for (segy = 0; segy<segNumy; segy++) yc[segy] = ((segy - hsegNumy) * segSize + hsegSize) * etaInterval; for (segx = 0; segx<segNumx; segx++) xc[segx] = (((segx - hsegNumx) * segSize) + hsegSize) * xiInterval; clock_t start, finish; double duration; start = clock(); // Iteration according to the point number for (no = 0; no<voxelnum; no++) { // point coordinate X = (_h_vox[no].x) * _CGHE->CGHScale; Y = (_h_vox[no].y) * _CGHE->CGHScale; Z = _h_vox[no].z * _CGHE->CGHScale - _CGHE->DefaultDepth; Amplitude = _h_vox[no].r; phase = _h_vox[no].ph; for (segy = 0; segy<segNumy; segy++) { theta_cy = (yc[segy] - Y) / Z; SFrequency_cy[segy] = (theta_cy + thetaY) / rLamda; (SFrequency_cy[segy] >= 0) ? PickPoint_cy[segy] = (int)(SFrequency_cy[segy] / sf_base + 0.5) : PickPoint_cy[segy] = (int)(SFrequency_cy[segy] / sf_base - 0.5); (abs(PickPoint_cy[segy])<FFThsegSize) ? Coefficient_cy[segy] = ((FFTsegSize - PickPoint_cy[segy]) % FFTsegSize) : Coefficient_cy[segy] = 0; } for (segx = 0; segx<segNumx; segx++) { theta_cx = (xc[segx] - X) / Z; SFrequency_cx[segx] = (theta_cx + thetaX) / rLamda; (SFrequency_cx[segx] >= 0) ? PickPoint_cx[segx] = (int)(SFrequency_cx[segx] / sf_base + 0.5) : PickPoint_cx[segx] = (int)(SFrequency_cx[segx] / sf_base - 0.5); (abs(PickPoint_cx[segx])<FFThsegSize) ? Coefficient_cx[segx] = ((FFTsegSize - PickPoint_cx[segx]) % FFTsegSize) : Coefficient_cx[segx] = 0; } for (segy = 0; segy<segNumy; segy++) { for (segx = 0; segx<segNumx; segx++) { segyy = segy*segNumx + segx; segxx = Coefficient_cy[segy] * FFTsegSize + Coefficient_cx[segx]; R = sqrt((xc[segx] - X)*(xc[segx] - X) + (yc[segy] - Y)*(yc[segy] - Y) + Z*Z); theta = rWaveNum * R + phase; theta_c = theta; theta_s = theta + PI; dtheta_c = ((int)(theta_c*NUMTBL / M2_PI)); dtheta_s = ((int)(theta_s*NUMTBL / M2_PI)); idx_c = (dtheta_c) & (NUMTBL2); idx_s = (dtheta_s) & (NUMTBL2); inRe[segyy][segxx] += (double)(Amplitude * m_COStbl[idx_c]); inIm[segyy][segxx] += (double)(Amplitude * m_SINtbl[idx_s]); } } } plan = fftw_plan_dft_2d(FFTsegSize, FFTsegSize, in, out, FFTW_BACKWARD, FFTW_ESTIMATE); for (segy = 0; segy<segNumy; segy++) { for (segx = 0; segx<segNumx; segx++) { segyy = segy*segNumx + segx; memset(in, 0x00, sizeof(fftw_complex) * FFTsegSize * FFTsegSize); for (i = 0; i <FFTsegSize; i++) { for (j = 0; j < FFTsegSize; j++) { segxx = i*FFTsegSize + j; in[i*FFTsegSize + j][0] = inRe[segyy][segxx]; in[i*FFTsegSize + j][1] = inIm[segyy][segxx]; } } fftw_execute(plan); for (i = 0; i <segSize; i++) { for (j = 0; j < segSize; j++) { m_pHologram[(segy*segSize + i)*cghwidth + (segx*segSize + j)] += out[(i + FFThsegSize - hsegSize) * FFTsegSize + (j + FFThsegSize - hsegSize)][0];// - out[l * SEGSIZE + m][1]; // m_pHologram[(segy*segSize+i)*cghwidth+(segx*segSize+j)] += atan2( out[i * FFTsegSize + j][1], out[i * FFTsegSize + j][0]);// - out[l * SEGSIZE + m][1]; } } } } finish = clock(); duration = (double)(finish - start) / CLOCKS_PER_SEC; cout << duration << endl; //mm.Format("%f", duration); //AfxMessageBox(mm); fftw_destroy_plan(plan); fftw_free(in); fftw_free(out); delete[] SFrequency_cx; delete[] SFrequency_cy; delete[] PickPoint_cx; delete[] PickPoint_cy; delete[] Coefficient_cx; delete[] Coefficient_cy; delete[] xc; delete[] yc; for (i = 0; i<segNumy; i++) { for (j = 0; j<segNumx; j++) { delete[] inRe[i*segNumx + j]; delete[] inIm[i*segNumx + j]; } } delete[] inRe; delete[] inIm; } int ophLUT::saveAsImg(const char * fname, uint8_t bitsperpixel, void * src, int pic_width, int pic_height) { LOG("Saving...%s...", fname); auto start = CUR_TIME; int _width = pic_width, _height = pic_height; int _pixelbytesize = _height * _width * bitsperpixel / 8; int _filesize = _pixelbytesize + sizeof(bitmap8bit); FILE *fp; fopen_s(&fp, fname, "wb"); if (fp == nullptr) return -1; bitmap8bit *pbitmap = (bitmap8bit*)calloc(1, sizeof(bitmap8bit)); memset(pbitmap, 0x00, sizeof(bitmap8bit)); pbitmap->fileheader.signature[0] = 'B'; pbitmap->fileheader.signature[1] = 'M'; pbitmap->fileheader.filesize = _filesize; pbitmap->fileheader.fileoffset_to_pixelarray = sizeof(bitmap8bit); for (int i = 0; i < 256; i++) { pbitmap->rgbquad[i].rgbBlue = i; pbitmap->rgbquad[i].rgbGreen = i; pbitmap->rgbquad[i].rgbRed = i; } pbitmap->bitmapinfoheader.dibheadersize = sizeof(bitmapinfoheader); pbitmap->bitmapinfoheader.width = _width; pbitmap->bitmapinfoheader.height = _height; //pbitmap->bitmapinfoheader.planes = _planes; pbitmap->bitmapinfoheader.bitsperpixel = bitsperpixel; //pbitmap->bitmapinfoheader.compression = _compression; pbitmap->bitmapinfoheader.imagesize = _pixelbytesize; //pbitmap->bitmapinfoheader.ypixelpermeter = _ypixelpermeter; //pbitmap->bitmapinfoheader.xpixelpermeter = _xpixelpermeter; pbitmap->bitmapinfoheader.numcolorspallette = 256; fwrite(pbitmap, 1, sizeof(bitmap8bit), fp); fwrite(src, 1, _pixelbytesize, fp); fclose(fp); free(pbitmap); auto end = CUR_TIME; auto during = ((std::chrono::duration<Real>)(end - start)).count(); LOG("%.5lfsec...done\n", during); return 0; } // 문자열 우측 공백문자 삭제 함수 char* ophLUT::rtrim(char* s) { char t[MAX_STR_LEN]; char *end; // Visual C 2003 이하에서는 // strcpy(t, s); // 이렇게 해야 함 strcpy_s(t, s); // 이것은 Visual C 2005용 end = t + strlen(t) - 1; while (end != t && isspace(*end)) end--; *(end + 1) = '\0'; s = t; return s; } // 문자열 좌측 공백문자 삭제 함수 char* ophLUT::ltrim(char* s) { char* begin; begin = s; while (*begin != '\0') { if (isspace(*begin)) begin++; else { s = begin; break; } } return s; } // 문자열 앞뒤 공백 모두 삭제 함수 char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
50  if (!(inputString[0] == NULL || (inputString[0] == '#' && inputString[1] == ' ')))
51  {
52  char* token = NULL;
53  char* parameter = NULL;
54  char* context = nullptr;
55  token = strtok_s(inputString, "=", &context);
56 
57  // 데이터 받아서 구조체에 저장
58  if (strcmp(token, "CGH width ") == 0) {
59  token = strtok_s(NULL, "=", &context);
60  conf->CghWidth = atoi(trim(token));
61  }
62  else if (strcmp(token, "CGH height ") == 0) {
63  token = strtok_s(NULL, "=", &context);
64  conf->CghHeight = atoi(trim(token));
65  }
66  else if (strcmp(token, "Segmentation size ") == 0) {
67  token = strtok_s(NULL, "=", &context);
68  conf->SegmentationSize = atoi(trim(token));
69  }
70  else if (strcmp(token, "FFT segmentation size ") == 0) {
71  token = strtok_s(NULL, "=", &context);
72  conf->fftSegmentationSize = atoi(trim(token));
73  }
74  else if (strcmp(token, "Red wavelength ") == 0) {
75  token = strtok_s(NULL, "=", &context);
76  conf->rWaveLength = atof(trim(token));
77  }
78  else if (strcmp(token, "Tilting angle on x axis ") == 0) {
79  token = strtok_s(NULL, "=", &context);
80  conf->ThetaX = atof(trim(token));
81  }
82  else if (strcmp(token, "Tilting angle on y axis ") == 0) {
83  token = strtok_s(NULL, "=", &context);
84  conf->ThetaY = atof(trim(token));
85  }
86  else if (strcmp(token, "Default depth ") == 0) {
87  token = strtok_s(NULL, "=", &context);
88  conf->DefaultDepth = atof(trim(token));
89  }
90  else if (strcmp(token, "3D point interval on x axis ") == 0) {
91  token = strtok_s(NULL, "=", &context);
92  conf->xInterval = atof(trim(token));
93  }
94  else if (strcmp(token, "3D point interval on y axis ") == 0) {
95  token = strtok_s(NULL, "=", &context);
96  conf->yInterval = atof(trim(token));
97  }
98  else if (strcmp(token, "Hologram interval on xi axis ") == 0) {
99  token = strtok_s(NULL, "=", &context);
100  conf->xiInterval = atof(trim(token));
101  }
102  else if (strcmp(token, "Hologram interval on eta axis ") == 0) {
103  token = strtok_s(NULL, "=", &context);
104  conf->etaInterval = atof(trim(token));
105  }
106  else if (strcmp(token, "CGH scale ") == 0) {
107  token = strtok_s(NULL, "=", &context);
108  conf->CGHScale = atof(trim(token));
109  }
110  }
111  }
112  }
113  inFile.close();
114  return false;
115 }
116 
117 bool ophLUT::loadPoint(const char* _filename, VoxelStruct* h_vox)
118 {
119 #define MAX_SIZE 1000
120  char inputString[MAX_SIZE];
121 
122  ifstream inFile(_filename);
123  if (!(inFile.is_open()))
124  {
125  cout << "포인트 클라우드 파일을 찾을 수 없습니다." << endl;
126  return false;
127  }
128  else
129  {
130  inFile.getline(inputString, MAX_SIZE);
131  int no = 0;
132  while (!inFile.eof())
133  {
134  inFile.getline(inputString, MAX_SIZE);
135 
136  if (inputString[0] != NULL)
137  {
138  char* token = NULL;
139  char* context;
140  token = strtok_s(inputString, "\t", &context);
141  h_vox[no].num = atoi(token); // 인덱스
142 
143  token = strtok_s(NULL, "\t", &context);
144  h_vox[no].x = atof(token); // x 좌표
145 
146  token = strtok_s(NULL, "\t", &context);
147  h_vox[no].y = atof(token); // y 좌표
148 
149  token = strtok_s(NULL, "\t", &context);
150  h_vox[no].z = atof(token); // z 좌표
151 
152  token = strtok_s(NULL, "\t", &context);
153  h_vox[no].ph = atof(token); // phase
154 
155  token = strtok_s(NULL, "\t", &context);
156  h_vox[no].r = atof(token); // red
157 
158  //token = strtok_s(NULL, "\t");
159  //h_vox[no].g = atof(token); // green
160 
161  //token = strtok_s(NULL, "\t");
162  //h_vox[no].b = atof(token); // blue
163 
164  no++;
165  }
166  }
167  }
168  inFile.close();
169  return true;
170 }
171 
172 bool ophLUT::load_Num_Point(const char* _filename, long* num_point)
173 {
174 #define MAX_SIZE 1000
175  char inputString[MAX_SIZE];
176  ifstream inFile(_filename);
177  if (!(inFile.is_open()))
178  {
179  cout << "포인트 클라우드 파일을 찾을 수 없습니다." << endl;
180  return false;
181  }
182  else
183  {
184  inFile.getline(inputString, MAX_SIZE);
185  *num_point = atoi(trim(inputString));
186  //cout << *num_point << endl;
187  }
188  inFile.close();
189  return true;
190 }
191 
193 {
194  m_pHologram = new double[_CGHE->CghHeight*_CGHE->CghWidth];
195  memset(m_pHologram, 0x00, sizeof(double)*_CGHE->CghHeight*_CGHE->CghWidth);
196 
197  for (int i = 0; i<NUMTBL; i++) {
198  float theta = (float)M2_PI * (float)(i + i - 1) / (float)(2 * NUMTBL);
199  m_COStbl[i] = (float)cos(theta);
200  m_SINtbl[i] = (float)sin(theta);
201  }
202 }
203 
204 // APAS
205 int ophLUT::LUTCalcuation(long voxnum, unsigned char *cghfringe, VoxelStruct* h_vox, CGHEnvironmentData* _CGHE)
206 {
207  long i, j;
208 
209  double Max = -1E9, Min = 1E9;
210  double myBuffer;
211  int cghwidth = _CGHE->CghWidth;
212  int cghheight = _CGHE->CghHeight;
213 
214  DataInit(_CGHE);
215 
216  //
217  APAS(voxnum, h_vox, _CGHE);
218  //
219 
220  for (i = 0; i<cghheight; i++) {
221  for (j = 0; j<cghwidth; j++) {
222  if (Max < m_pHologram[i*cghwidth + j]) Max = m_pHologram[i*cghwidth + j];
223  if (Min > m_pHologram[i*cghwidth + j]) Min = m_pHologram[i*cghwidth + j];
224  }
225  }
226 
227  for (i = 0; i<cghheight; i++) {
228  for (j = 0; j<cghwidth; j++) {
229  myBuffer = 1.0*(((m_pHologram[i*cghwidth + j] - Min) / (Max - Min))*255. + 0.5);
230  if (myBuffer >= 255.0) cghfringe[i*cghwidth + j] = 255;
231  else cghfringe[i*cghwidth + j] = (unsigned char)(myBuffer);
232  }
233  }
234 
235  return 0;
236 }
237 
238 void ophLUT::APAS(long voxelnum, VoxelStruct * _h_vox, CGHEnvironmentData * _CGHE)
239 {
240  int no; // voxel Number
241  int i, j;
242  int segx, segy; // coordinate in a Segment
243  float R;
244 
245  int cghwidth = _CGHE->CghWidth;
246  int cghheight = _CGHE->CghHeight;
247  float xiInterval = _CGHE->xiInterval;
248  float etaInterval = _CGHE->etaInterval;
249  float rLamda = _CGHE->rWaveLength;
250  float rWaveNum = _CGHE->rWaveNumber;
251  float thetaX = _CGHE->ThetaX;
252  float thetaY = _CGHE->ThetaY;
253 
254  int segSize = _CGHE->SegmentationSize;
255  int hsegSize = (int)(segSize / 2);
256  int dsegSize = segSize*segSize;
257 
258  int segNumx = (int)(cghwidth / segSize);
259  int segNumy = (int)(cghheight / segSize);
260  int hsegNumx = (int)(segNumx / 2);
261  int hsegNumy = (int)(segNumy / 2);
262 
263  int FFTsegSize = _CGHE->fftSegmentationSize;
264  int FFThsegSize = (int)(FFTsegSize / 2);
265  int FFTdsegSize = FFTsegSize*FFTsegSize;
266 
267  float X, Y, Z; ; // x, y, real distance
268  float theta_cx, theta_cy;
269  float *SFrequency_cx = new float[segNumx]; // Center phase
270  float *SFrequency_cy = new float[segNumy]; // Center phase
271  int *PickPoint_cx = new int[segNumx];
272  int *PickPoint_cy = new int[segNumy];
273  int *Coefficient_cx = new int[segNumx];
274  int *Coefficient_cy = new int[segNumy];
275  float *xc = new float[segNumx];
276  float *yc = new float[segNumy];
277  float Amplitude;
278  float phase;
279  float sf_base = 1.0 / (xiInterval*FFTsegSize);
280  int segxx, segyy;
281  float theta;
282  float theta_s, theta_c;
283  int dtheta_s, dtheta_c;
284  int idx_c, idx_s;
285 
286  //CString mm;
287 
288  fftw_complex *in, *out;
289  fftw_plan plan;
290 
291  double **inRe = new double *[segNumy * segNumx];
292  double **inIm = new double *[segNumy * segNumx];
293  for (i = 0; i<segNumy; i++) {
294  for (j = 0; j<segNumx; j++) {
295  inRe[i*segNumx + j] = new double[FFTsegSize * FFTsegSize];
296  inIm[i*segNumx + j] = new double[FFTsegSize * FFTsegSize];
297  memset(inRe[i*segNumx + j], 0x00, sizeof(double) * FFTsegSize * FFTsegSize);
298  memset(inIm[i*segNumx + j], 0x00, sizeof(double) * FFTsegSize * FFTsegSize);
299  }
300  }
301 
302  in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTsegSize * FFTsegSize);
303  out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTsegSize * FFTsegSize);
304  memset(in, 0x00, sizeof(fftw_complex) * FFTsegSize * FFTsegSize);
305  memset(m_pHologram, 0x00, sizeof(double)*cghwidth*cghheight);
306 
307  for (segy = 0; segy<segNumy; segy++)
308  yc[segy] = ((segy - hsegNumy) * segSize + hsegSize) * etaInterval;
309  for (segx = 0; segx<segNumx; segx++)
310  xc[segx] = (((segx - hsegNumx) * segSize) + hsegSize) * xiInterval;
311 
312  clock_t start, finish;
313  double duration;
314  start = clock();
315 
316  // Iteration according to the point number
317  for (no = 0; no<voxelnum; no++)
318  {
319  // point coordinate
320  X = (_h_vox[no].x) * _CGHE->CGHScale;
321  Y = (_h_vox[no].y) * _CGHE->CGHScale;
322  Z = _h_vox[no].z * _CGHE->CGHScale - _CGHE->DefaultDepth;
323  Amplitude = _h_vox[no].r;
324  phase = _h_vox[no].ph;
325 
326  for (segy = 0; segy<segNumy; segy++)
327  {
328  theta_cy = (yc[segy] - Y) / Z;
329  SFrequency_cy[segy] = (theta_cy + thetaY) / rLamda;
330  (SFrequency_cy[segy] >= 0) ? PickPoint_cy[segy] = (int)(SFrequency_cy[segy] / sf_base + 0.5)
331  : PickPoint_cy[segy] = (int)(SFrequency_cy[segy] / sf_base - 0.5);
332  (abs(PickPoint_cy[segy])<FFThsegSize) ? Coefficient_cy[segy] = ((FFTsegSize - PickPoint_cy[segy]) % FFTsegSize)
333  : Coefficient_cy[segy] = 0;
334  }
335 
336  for (segx = 0; segx<segNumx; segx++)
337  {
338  theta_cx = (xc[segx] - X) / Z;
339  SFrequency_cx[segx] = (theta_cx + thetaX) / rLamda;
340  (SFrequency_cx[segx] >= 0) ? PickPoint_cx[segx] = (int)(SFrequency_cx[segx] / sf_base + 0.5)
341  : PickPoint_cx[segx] = (int)(SFrequency_cx[segx] / sf_base - 0.5);
342  (abs(PickPoint_cx[segx])<FFThsegSize) ? Coefficient_cx[segx] = ((FFTsegSize - PickPoint_cx[segx]) % FFTsegSize)
343  : Coefficient_cx[segx] = 0;
344  }
345  for (segy = 0; segy<segNumy; segy++) {
346  for (segx = 0; segx<segNumx; segx++) {
347  segyy = segy*segNumx + segx;
348  segxx = Coefficient_cy[segy] * FFTsegSize + Coefficient_cx[segx];
349  R = sqrt((xc[segx] - X)*(xc[segx] - X) + (yc[segy] - Y)*(yc[segy] - Y) + Z*Z);
350  theta = rWaveNum * R + phase;
351  theta_c = theta;
352  theta_s = theta + PI;
353  dtheta_c = ((int)(theta_c*NUMTBL / M2_PI));
354  dtheta_s = ((int)(theta_s*NUMTBL / M2_PI));
355  idx_c = (dtheta_c) & (NUMTBL2);
356  idx_s = (dtheta_s) & (NUMTBL2);
357  inRe[segyy][segxx] += (double)(Amplitude * m_COStbl[idx_c]);
358  inIm[segyy][segxx] += (double)(Amplitude * m_SINtbl[idx_s]);
359  }
360  }
361  }
362 
363  plan = fftw_plan_dft_2d(FFTsegSize, FFTsegSize, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
364 
365  for (segy = 0; segy<segNumy; segy++) {
366  for (segx = 0; segx<segNumx; segx++) {
367  segyy = segy*segNumx + segx;
368  memset(in, 0x00, sizeof(fftw_complex) * FFTsegSize * FFTsegSize);
369  for (i = 0; i <FFTsegSize; i++) {
370  for (j = 0; j < FFTsegSize; j++) {
371  segxx = i*FFTsegSize + j;
372  in[i*FFTsegSize + j][0] = inRe[segyy][segxx];
373  in[i*FFTsegSize + j][1] = inIm[segyy][segxx];
374  }
375  }
376  fftw_execute(plan);
377  for (i = 0; i <segSize; i++) {
378  for (j = 0; j < segSize; j++) {
379  m_pHologram[(segy*segSize + i)*cghwidth + (segx*segSize + j)] +=
380  out[(i + FFThsegSize - hsegSize) * FFTsegSize + (j + FFThsegSize - hsegSize)][0];// - out[l * SEGSIZE + m][1];
381  // m_pHologram[(segy*segSize+i)*cghwidth+(segx*segSize+j)] += atan2( out[i * FFTsegSize + j][1], out[i * FFTsegSize + j][0]);// - out[l * SEGSIZE + m][1];
382  }
383  }
384  }
385  }
386  finish = clock();
387 
388  duration = (double)(finish - start) / CLOCKS_PER_SEC;
389  cout << duration << endl;
390  //mm.Format("%f", duration);
391  //AfxMessageBox(mm);
392 
393  fftw_destroy_plan(plan);
394  fftw_free(in);
395  fftw_free(out);
396  delete[] SFrequency_cx;
397  delete[] SFrequency_cy;
398  delete[] PickPoint_cx;
399  delete[] PickPoint_cy;
400  delete[] Coefficient_cx;
401  delete[] Coefficient_cy;
402  delete[] xc;
403  delete[] yc;
404  for (i = 0; i<segNumy; i++) {
405  for (j = 0; j<segNumx; j++) {
406  delete[] inRe[i*segNumx + j];
407  delete[] inIm[i*segNumx + j];
408  }
409  }
410  delete[] inRe;
411  delete[] inIm;
412 }
413 
414 int ophLUT::saveAsImg(const char * fname, uint8_t bitsperpixel, void * src, int pic_width, int pic_height)
415 {
416  LOG("Saving...%s...", fname);
417  auto start = CUR_TIME;
418 
419  int _width = pic_width, _height = pic_height;
420 
421  int _pixelbytesize = _height * _width * bitsperpixel / 8;
422  int _filesize = _pixelbytesize + sizeof(bitmap8bit);
423 
424  FILE *fp;
425  fopen_s(&fp, fname, "wb");
426  if (fp == nullptr) return -1;
427 
428  bitmap8bit *pbitmap = (bitmap8bit*)calloc(1, sizeof(bitmap8bit));
429  memset(pbitmap, 0x00, sizeof(bitmap8bit));
430 
431  pbitmap->fileheader.signature[0] = 'B';
432  pbitmap->fileheader.signature[1] = 'M';
433  pbitmap->fileheader.filesize = _filesize;
435 
436  for (int i = 0; i < 256; i++) {
437  pbitmap->rgbquad[i].rgbBlue = i;
438  pbitmap->rgbquad[i].rgbGreen = i;
439  pbitmap->rgbquad[i].rgbRed = i;
440  }
441 
443  pbitmap->bitmapinfoheader.width = _width;
444  pbitmap->bitmapinfoheader.height = _height;
445  //pbitmap->bitmapinfoheader.planes = _planes;
446  pbitmap->bitmapinfoheader.bitsperpixel = bitsperpixel;
447  //pbitmap->bitmapinfoheader.compression = _compression;
448  pbitmap->bitmapinfoheader.imagesize = _pixelbytesize;
449  //pbitmap->bitmapinfoheader.ypixelpermeter = _ypixelpermeter;
450  //pbitmap->bitmapinfoheader.xpixelpermeter = _xpixelpermeter;
451  pbitmap->bitmapinfoheader.numcolorspallette = 256;
452  fwrite(pbitmap, 1, sizeof(bitmap8bit), fp);
453 
454  fwrite(src, 1, _pixelbytesize, fp);
455  fclose(fp);
456  free(pbitmap);
457 
458  auto end = CUR_TIME;
459 
460  auto during = ((std::chrono::duration<Real>)(end - start)).count();
461 
462  LOG("%.5lfsec...done\n", during);
463 
464  return 0;
465 }
466 
467 // 문자열 우측 공백문자 삭제 함수char* ophLUT::rtrim(char* s) { char t[MAX_STR_LEN]; char *end; // Visual C 2003 이하에서는 // strcpy(t, s); // 이렇게 해야 함 strcpy_s(t, s); // 이것은 Visual C 2005용 end = t + strlen(t) - 1; while (end != t && isspace(*end)) end--; *(end + 1) = '\0'; s = t; return s; } // 문자열 좌측 공백문자 삭제 함수 char* ophLUT::ltrim(char* s) { char* begin; begin = s; while (*begin != '\0') { if (isspace(*begin)) begin++; else { s = begin; break; } } return s; } // 문자열 앞뒤 공백 모두 삭제 함수 char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
468 char* ophLUT::rtrim(char* s)
469 {
470  char t[MAX_STR_LEN];
471  char *end;
472 
473  // Visual C 2003 이하에서는 // strcpy(t, s); // 이렇게 해야 함 strcpy_s(t, s); // 이것은 Visual C 2005용 end = t + strlen(t) - 1; while (end != t && isspace(*end)) end--; *(end + 1) = '\0'; s = t; return s; } // 문자열 좌측 공백문자 삭제 함수 char* ophLUT::ltrim(char* s) { char* begin; begin = s; while (*begin != '\0') { if (isspace(*begin)) begin++; else { s = begin; break; } } return s; } // 문자열 앞뒤 공백 모두 삭제 함수 char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
474  // strcpy(t, s);
475  // 이렇게 해야 함
476  strcpy_s(t, s); // 이것은 Visual C 2005용 end = t + strlen(t) - 1; while (end != t && isspace(*end)) end--; *(end + 1) = '\0'; s = t; return s; } // 문자열 좌측 공백문자 삭제 함수 char* ophLUT::ltrim(char* s) { char* begin; begin = s; while (*begin != '\0') { if (isspace(*begin)) begin++; else { s = begin; break; } } return s; } // 문자열 앞뒤 공백 모두 삭제 함수 char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
477  end = t + strlen(t) - 1;
478  while (end != t && isspace(*end))
479  end--;
480  *(end + 1) = '\0';
481  s = t;
482 
483  return s;
484 }
485 
486 // 문자열 좌측 공백문자 삭제 함수char* ophLUT::ltrim(char* s) { char* begin; begin = s; while (*begin != '\0') { if (isspace(*begin)) begin++; else { s = begin; break; } } return s; } // 문자열 앞뒤 공백 모두 삭제 함수 char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
487 char* ophLUT::ltrim(char* s)
488 {
489  char* begin;
490  begin = s;
491 
492  while (*begin != '\0') {
493  if (isspace(*begin))
494  begin++;
495  else {
496  s = begin;
497  break;
498  }
499  }
500 
501  return s;
502 }
503 
504 
505 // 문자열 앞뒤 공백 모두 삭제 함수char* ophLUT::trim(char* s) { return rtrim(ltrim(s)); }
506 char* ophLUT::trim(char* s)
507 {
508  return rtrim(ltrim(s));
509 }
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:113
void APAS(long voxelnum, VoxelStruct *_h_vox, CGHEnvironmentData *_CGHE)
Definition: ophLUT.cpp:238
uint32_t dibheadersize
Definition: struct.h:57
float ph
Definition: ophACPAS.h:99
bool loadConfig(const char *filename, CGHEnvironmentData *_CGHE)
Definition: ophLUT.cpp:33
char * ltrim(char *s)
Definition: ophLUT.cpp:487
#define MAX_SIZE
fileheader fileheader
Definition: struct.h:80
float m_SINtbl[NUMTBL]
Definition: ophLUT.h:53
#define PI
Definition: ophACPAS.h:9
crop phase
double * m_pHologram
Definition: ophLUT.h:50
bool load_Num_Point(const char *_filename, long *num_point)
Definition: ophLUT.cpp:172
uint32_t imagesize
Definition: struct.h:63
#define FFTW_BACKWARD
Definition: fftw3.h:380
void Amplitude(Complex< Real > *holo, Real *encoded, const int size)
Definition: ophGen.cpp:2154
int saveAsImg(const char *fname, uint8_t bitsperpixel, void *src, int pic_width, int pic_height)
Definition: ophLUT.cpp:414
uint8_t signature[2]
Definition: struct.h:51
bitmapinfoheader bitmapinfoheader
Definition: struct.h:81
#define NUMTBL
Definition: ophACPAS.h:15
uint8_t rgbGreen
Definition: struct.h:71
#define CUR_TIME
Definition: function.h:58
uint16_t bitsperpixel
Definition: struct.h:61
#define M2_PI
Definition: ophACPAS.h:10
uint8_t rgbBlue
Definition: struct.h:70
float r
Definition: ophACPAS.h:100
char * rtrim(char *s)
Definition: ophLUT.cpp:468
bool loadPoint(const char *_filename, VoxelStruct *h_vox)
Definition: ophLUT.cpp:117
fclose('all')
float m_COStbl[NUMTBL]
Definition: ophLUT.h:52
uint32_t numcolorspallette
Definition: struct.h:66
rgbquad rgbquad[256]
Definition: struct.h:82
#define MAX_STR_LEN
Definition: ophACPAS.h:17
uint32_t filesize
Definition: struct.h:52
virtual ~ophLUT()
Definition: ophLUT.cpp:19
uint32_t fileoffset_to_pixelarray
Definition: struct.h:54
uint8_t rgbRed
Definition: struct.h:72
uint32_t width
Definition: struct.h:58
char * trim(char *s)
Definition: ophLUT.cpp:506
uint32_t height
Definition: struct.h:59
ophLUT()
Definition: ophLUT.cpp:15
int LUTCalcuation(long voxnum, unsigned char *cghfringe, VoxelStruct *h_vox, CGHEnvironmentData *_CGHE)
Definition: ophLUT.cpp:205
float x
Definition: ophACPAS.h:96
int init(const char *_filename, CGHEnvironmentData *_CGHE)
Definition: ophLUT.cpp:24
#define FFTW_ESTIMATE
Definition: fftw3.h:392
float y
Definition: ophACPAS.h:97
void DataInit(CGHEnvironmentData *_CGHE)
Definition: ophLUT.cpp:192
Definition: ophGen.h:68
#define NUMTBL2
Definition: ophACPAS.h:16
float z
Definition: ophACPAS.h:98