Openholo  v5.0
Open Source Digital Holographic Library
ophSig.cpp
Go to the documentation of this file.
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install, copy or use the software.
7 //
8 //
9 // License Agreement
10 // For Open Source Digital Holographic Library
11 //
12 // Openholo library is free software;
13 // you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
14 //
15 // Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
16 // E-mail : contact.openholo@gmail.com
17 // Web : http://www.openholo.org
18 //
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
21 //
22 // 1. Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
24 //
25 // 2. Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other materials provided with the distribution.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the copyright holder or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 // This software contains opensource software released under GNU Generic Public License,
41 // NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
42 // Check whether software you use contains licensed software.
43 //
44 //M*/
45 
46 #include "ophSig.h"
47 #include "include.h"
48 
49 
50 ophSig::ophSig(void)
51  : is_CPU(true)
52 {
53  _foc[0] = _foc[1] = _foc[2] = 0.0f;
54 }
55 
56 void ophSig::cField2Buffer(matrix<Complex<Real>>& src, Complex<Real> **dst,int nx,int ny) {
57  ivec2 bufferSize(nx, ny); //= src.getSize();
58 
59  *dst = new oph::Complex<Real>[nx * ny];
60 
61  int idx = 0;
62 
63  for (int x = 0; x < bufferSize[_X]; x++) {
64  for (int y = 0; y < bufferSize[_Y]; y++) {
65  (*dst)[idx] = src[x][y];
66  idx++;
67  }
68  }
69 }
70 
71 void ophSig::ColorField2Buffer(matrix<Complex<Real>>& src, Complex<Real> **dst, int nx, int ny) {
72  ivec2 bufferSize(nx, ny); //= src.getSize();
73 
74  *dst = new oph::Complex<Real>[3*nx*ny];
75 
76  int idx = 0;
77  for (int x = 0; x < bufferSize[_X]; x++) {
78  for (int y = 0; y < bufferSize[_Y]; y++) {
79  (*dst)[idx] = src[x][y];
80  idx++;
81  }
82  }
83 }
84 void ophSig::setMode(bool is_CPU)
85 {
86  this->is_CPU = is_CPU;
87 }
88 
89 template<typename T>
90 void ophSig::linInterp(vector<T> &X, matrix<Complex<T>> &src, vector<T> &Xq, matrix<Complex<T>> &dst)
91 {
92  if (src.size != dst.size) {
93  dst.resize(src.size[_X], src.size[_Y]);
94  }
95  int size = src.size[_Y];
96 
97  for (int i = 0, j = 0; j < dst.size[_Y]; j++)
98  {
99  if ((Xq[j]) >= (X[size - 2]))
100  {
101  i = size - 2;
102  }
103  else
104  {
105  while ((Xq[j]) >(X[i + 1])) i++;
106  }
107  dst(0, j)[_RE] = src(0, i).real() + (src(0, i + 1).real() - src(0, i).real()) / (X[i + 1] - X[i]) * (Xq[j] - X[i]);
108  dst(0, j)[_IM] = src(0, i).imag() + (src(0, i + 1).imag() - src(0, i).imag()) / (X[i + 1] - X[i]) * (Xq[j] - X[i]);
109  }
110 }
111 
112 template<typename T>
113 void ophSig::fft1(matrix<Complex<T>> &src, matrix<Complex<T>> &dst, int sign, uint flag)
114 {
115  if (src.size != dst.size) {
116  dst.resize(src.size[_X], src.size[_Y]);
117  }
118  fftw_complex *fft_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * src.size[_Y]);
119  fftw_complex *fft_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * src.size[_Y]);
120 
121  for (int i = 0; i < src.size[_Y]; i++) {
122  fft_in[i][_RE] = src(0, i).real();
123  fft_in[i][_IM] = src(0, i).imag();
124  }
125 
126  fftw_plan plan = fftw_plan_dft_1d(src.size[_Y], fft_in, fft_out, sign, flag);
127 
128  fftw_execute(plan);
129  if (sign == OPH_FORWARD)
130  {
131  for (int i = 0; i < src.size[_Y]; i++) {
132  dst(0, i)[_RE] = fft_out[i][_RE];
133  dst(0, i)[_IM] = fft_out[i][_IM];
134  }
135  }
136  else if (sign == OPH_BACKWARD)
137  {
138  for (int i = 0; i < src.size[_Y]; i++) {
139  dst(0, i)[_RE] = fft_out[i][_RE] / src.size[_Y];
140  dst(0, i)[_IM] = fft_out[i][_IM] / src.size[_Y];
141  }
142  }
143 
144  fftw_destroy_plan(plan);
145  fftw_free(fft_in);
146  fftw_free(fft_out);
147 }
148 template<typename T>
149 void ophSig::fft2(matrix<Complex<T>> &src, matrix<Complex<T>> &dst, int sign, uint flag)
150 {
151  if (src.size != dst.size) {
152  dst.resize(src.size[_X], src.size[_Y]);
153  }
154 
155  fftw_complex *fft_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * src.size[_X] * src.size[_Y]);
156  fftw_complex *fft_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * src.size[_X] * src.size[_Y]);
157 
158  for (int i = 0; i < src.size[_X]; i++) {
159  for (int j = 0; j < src.size[_Y]; j++) {
160  fft_in[src.size[_Y] * i + j][_RE] = src(i, j).real();
161  fft_in[src.size[_Y] * i + j][_IM] = src(i, j).imag();
162  }
163  }
164 
165  fftw_plan plan = fftw_plan_dft_2d(src.size[_X], src.size[_Y], fft_in, fft_out, sign, flag);
166 
167  fftw_execute(plan);
168  if (sign == OPH_FORWARD)
169  {
170  for (int i = 0; i < src.size[_X]; i++) {
171  for (int j = 0; j < src.size[_Y]; j++) {
172  dst(i, j)[_RE] = fft_out[src.size[_Y] * i + j][_RE];
173  dst(i, j)[_IM] = fft_out[src.size[_Y] * i + j][_IM];
174  }
175  }
176  }
177  else if (sign == OPH_BACKWARD)
178  {
179  for (int i = 0; i < src.size[_X]; i++) {
180  for (int j = 0; j < src.size[_Y]; j++) {
181  dst(i, j)[_RE] = fft_out[src.size[_Y] * i + j][_RE] / (src.size[_X] * src.size[_Y]);
182  dst(i, j)[_IM] = fft_out[src.size[_Y] * i + j][_IM] / (src.size[_X] * src.size[_Y]);
183 
184  }
185  }
186  }
187 
188  fftw_destroy_plan(plan);
189  fftw_free(fft_in);
190  fftw_free(fft_out);
191 }
192 
193 
194 
195 
196 
197 
198 bool ophSig::loadAsOhc(const char *fname)
199 {
200  std::string fullname = fname;
201  if (!checkExtension(fname, ".ohc")) fullname.append(".ohc");
202  OHC_decoder->setFileName(fullname.c_str());
203 
204  if (!OHC_decoder->load()) return false;
205  vector<Real> wavelengthArray;
206  OHC_decoder->getWavelength(wavelengthArray);
207  _wavelength_num = OHC_decoder->getNumOfWavlen();
208  int wavelength_num = OHC_decoder->getNumOfWavlen();
209 
210  context_.pixel_number = OHC_decoder->getNumOfPixel();
211 
212  context_.wave_length = new Real[_wavelength_num];
213 
214  ComplexH = new OphComplexField[_wavelength_num];
215 
216  for (int i = 0; i < _wavelength_num; i++)
217  {
218  context_.wave_length[i] = wavelengthArray[(_wavelength_num - 1) - i];
219 
220  ComplexH[i].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
221  OHC_decoder->getComplexFieldData(ComplexH[i], (_wavelength_num - 1) - i);
222  //OHC_decoder->getComplexFieldData(&ComplexH);
223  }
224  return true;
225 }
226 
227 bool ophSig::saveAsOhc(const char *fname)
228 {
229  std::string fullname = fname;
230  if (!checkExtension(fname, ".ohc")) fullname.append(".ohc");
231  OHC_encoder->setFileName(fullname.c_str());
232 
233  ohcHeader header;
234 
235  OHC_encoder->getOHCheader(header);
236 
237  OHC_encoder->setNumOfPixel(context_.pixel_number[_X], context_.pixel_number[_Y]);
238 
239  OHC_encoder->setFieldEncoding(FldStore::Directly, FldCodeType::RI);
240 
241  OHC_encoder->setNumOfWavlen(_wavelength_num);
242 
243  for (int i = _wavelength_num - 1; i >= 0; i--)
244  {
245  //int wl = context_.wave_length[i] * 1000000000;
246  OHC_encoder->setWavelength(context_.wave_length[i], LenUnit::nm);
247 
248  OHC_encoder->addComplexFieldData(ComplexH[i]);
249  }
250 
251  if (!OHC_encoder->save()) return false;
252 
253 
254  return true;
255 }
256 
257 bool ophSig::load(const char *real, const char *imag)
258 {
259  string realname = real;
260  string imagname = imag;
261 
262  char* RGB_name[3] = { 0, };
263  const char* postfix = "";
264 
265  if (_wavelength_num > 1) {
266  postfix = "_B";
267  strcpy(RGB_name[0], postfix);
268  postfix = "_G";
269  strcpy(RGB_name[1], postfix);
270  postfix = "_R";
271  strcpy(RGB_name[2], postfix);
272  }
273  else
274  strcpy(RGB_name[0], postfix);
275 
276  int checktype = static_cast<int>(realname.rfind("."));
277 
278  OphRealField* realMat = new OphRealField[_wavelength_num];
279  OphRealField* imagMat = new OphRealField[_wavelength_num];
280 
281  std::string realtype = realname.substr(checktype + 1, realname.size());
282  std::string imgtype = imagname.substr(checktype + 1, realname.size());
283 
284  ComplexH = new OphComplexField[_wavelength_num];
285 
286  if (realtype != imgtype) {
287  LOG("failed : The data type between real and imaginary is different!\n");
288  return false;
289  }
290  if (realtype == "bmp")
291  {
292  realname = real;
293  imagname = imag;
294 
295  oph::uchar* realdata = loadAsImg(realname.c_str());
296  oph::uchar* imagdata = loadAsImg(imagname.c_str());
297 
298  if (realdata == 0 && imagdata == 0) {
299  cout << "failed : hologram data load was failed." << endl;
300  return false;
301  }
302 
303  for (int z = 0; z < _wavelength_num; z++)
304  {
305  realMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
306  imagMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
307  for (int i = context_.pixel_number[_X] - 1; i >= 0; i--)
308  {
309  for (int j = 0; j < context_.pixel_number[_Y]; j++)
310  {
311  realMat[z](context_.pixel_number[_X] - i - 1, j) = (double)realdata[(i * context_.pixel_number[_Y] + j)*_wavelength_num + z];
312  imagMat[z](context_.pixel_number[_X] - i - 1, j) = (double)imagdata[(i * context_.pixel_number[_Y] + j)*_wavelength_num + z];
313  }
314  }
315  }
316  delete[] realdata;
317  delete[] imagdata;
318  }
319  else if (realtype == "bin")
320  {
321  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
322  double *imagdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
323 
324  for (int z = 0; z < _wavelength_num; z++)
325  {
326  realname = real;
327  imagname = imag;
328 
329  realname.insert(checktype, RGB_name[z]);
330  imagname.insert(checktype, RGB_name[z]);
331 
332  ifstream freal(realname.c_str(), ifstream::binary);
333  ifstream fimag(imagname.c_str(), ifstream::binary);
334 
335  freal.read(reinterpret_cast<char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
336  fimag.read(reinterpret_cast<char*>(imagdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
337 
338  realMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
339  imagMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
340 
341  for (int i = 0; i < context_.pixel_number[_X]; i++)
342  {
343  for (int j = 0; j < context_.pixel_number[_Y]; j++)
344  {
345  realMat[z](i, j) = realdata[i + j * context_.pixel_number[_X]];
346  imagMat[z](i, j) = imagdata[i + j * context_.pixel_number[_X]];
347  }
348  }
349  freal.close();
350  fimag.close();
351  }
352  delete[] realdata;
353  delete[] imagdata;
354  }
355  else
356  {
357  LOG("Error: wrong type\n");
358  }
359  //nomalization
360  //double realout, imagout;
361  for (int z = 0; z < _wavelength_num; z++)
362  {
363  /*meanOfMat(realMat[z], realout); meanOfMat(imagMat[z], imagout);
364 
365  realMat[z] / realout; imagMat[z] / imagout;
366  absMat(realMat[z], realMat[z]);
367  absMat(imagMat[z], imagMat[z]);
368  realout = maxOfMat(realMat[z]); imagout = maxOfMat(imagMat[z]);
369  realMat[z] / realout; imagMat[z] / imagout;
370  realout = minOfMat(realMat[z]); imagout = minOfMat(imagMat[z]);
371  realMat[z] - realout; imagMat[z] - imagout;*/
372 
373  ComplexH[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
374 
375  for (int i = 0; i < context_.pixel_number[_X]; i++)
376  {
377  for (int j = 0; j < context_.pixel_number[_Y]; j++)
378  {
379  ComplexH[z](i, j)[_RE] = realMat[z](i, j);
380  ComplexH[z](i, j)[_IM] = imagMat[z](i, j);
381  }
382  }
383  }
384  LOG("Reading Openholo Complex Field File...%s, %s\n", realname.c_str(), imagname.c_str());
385 
386  return true;
387 }
388 
389 
390 
391 bool ophSig::save(const char *real, const char *imag)
392 {
393  string realname = real;
394  string imagname = imag;
395 
396  char* RGB_name[3] = { 0, };
397  const char* postfix = "";
398 
399  if (_wavelength_num > 1) {
400  postfix = "_B";
401  strcpy(RGB_name[0], postfix);
402  postfix = "_G";
403  strcpy(RGB_name[1], postfix);
404  postfix = "_R";
405  strcpy(RGB_name[2], postfix);
406  }
407  else
408  strcpy(RGB_name[0], postfix);
409 
410  int checktype = static_cast<int>(realname.rfind("."));
411  string type = realname.substr(checktype + 1, realname.size());
412  if (type == "bin")
413  {
414  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
415  double *imagdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
416 
417  for (int z = 0; z < _wavelength_num; z++)
418  {
419  realname = real;
420  imagname = imag;
421  realname.insert(checktype, RGB_name[z]);
422  imagname.insert(checktype, RGB_name[z]);
423  std::ofstream cos(realname.c_str(), std::ios::binary);
424  std::ofstream sin(imagname.c_str(), std::ios::binary);
425 
426  if (!cos.is_open()) {
427  LOG("real file not found.\n");
428  cos.close();
429  delete[] realdata;
430  delete[] imagdata;
431  return false;
432  }
433 
434  if (!sin.is_open()) {
435  LOG("imag file not found.\n");
436  sin.close();
437  delete[] realdata;
438  delete[] imagdata;
439  return false;
440  }
441 
442  for (int i = 0; i < context_.pixel_number[_X]; i++)
443  {
444  for (int j = 0; j < context_.pixel_number[_Y]; j++)
445  {
446  realdata[i + j * context_.pixel_number[_X]] = ComplexH[z](i, j)[_RE];
447  imagdata[i + j * context_.pixel_number[_X]] = ComplexH[z](i, j)[_IM];
448  }
449  }
450  cos.write(reinterpret_cast<const char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
451  sin.write(reinterpret_cast<const char*>(imagdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
452 
453  cos.close();
454  sin.close();
455  }
456  delete[]realdata;
457  delete[]imagdata;
458 
459  LOG("Writing Openholo Complex Field...%s, %s\n", realname.c_str(), imagname.c_str());
460  }
461  else if (type == "bmp")
462  {
463  oph::uchar* realdata;
464  oph::uchar* imagdata;
465  int _pixelbytesize = 0;
466  int _width = context_.pixel_number[_Y], _height = context_.pixel_number[_X];
467  int bitpixel = _wavelength_num * 8;
468 
469  if (bitpixel == 8)
470  {
471  _pixelbytesize = _height * _width;
472  }
473  else
474  {
475  _pixelbytesize = _height * _width * 3;
476  }
477  int _filesize = 0;
478 
479 
480  FILE *freal, *fimag;
481  freal = fopen(realname.c_str(), "wb");
482  fimag = fopen(imagname.c_str(), "wb");
483 
484  if ((freal == nullptr) || (fimag == nullptr))
485  {
486  LOG("file not found\n");
487  return false;
488  }
489 
490  if (bitpixel == 8)
491  {
492  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height);
493  imagdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height);
494  _filesize = _pixelbytesize + sizeof(bitmap8bit);
495 
496  bitmap8bit *pbitmap = (bitmap8bit*)calloc(1, sizeof(bitmap8bit));
497  memset(pbitmap, 0x00, sizeof(bitmap8bit));
498 
499  pbitmap->_fileheader.signature[0] = 'B';
500  pbitmap->_fileheader.signature[1] = 'M';
501  pbitmap->_fileheader.filesize = _filesize;
502  pbitmap->_fileheader.fileoffset_to_pixelarray = sizeof(bitmap8bit);
503 
504  for (int i = 0; i < 256; i++) {
505  pbitmap->_rgbquad[i].rgbBlue = i;
506  pbitmap->_rgbquad[i].rgbGreen = i;
507  pbitmap->_rgbquad[i].rgbRed = i;
508  }
509 
510 
512  for (int i = _height - 1; i >= 0; i--)
513  {
514  for (int j = 0; j < _width; j++)
515  {
516  if (ComplexH[0].mat[_height - i - 1][j][_RE] < 0)
517  {
518  ComplexH[0].mat[_height - i - 1][j][_RE] = 0;
519  }
520 
521  if (ComplexH[0].mat[_height - i - 1][j][_IM] < 0)
522  {
523  ComplexH[0].mat[_height - i - 1][j][_IM] = 0;
524  }
525  }
526  }
527 
528  double minVal, iminVal, maxVal, imaxVal;
529  maxVal = minVal = ComplexH[0](0, 0)[_RE];
530  imaxVal = iminVal = ComplexH[0](0, 0)[_IM];
531 
532  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
533  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
534  if (ComplexH[0](i, j)[_RE] < minVal) minVal = ComplexH[0](i, j).real();
535  if (ComplexH[0](i, j)[_RE] > maxVal) maxVal = ComplexH[0](i, j).real();
536  if (ComplexH[0](i, j)[_IM] < iminVal) iminVal = ComplexH[0](i, j)[_IM];
537  if (ComplexH[0](i, j)[_IM] > imaxVal) imaxVal = ComplexH[0](i, j)[_IM];
538  }
539  }
540  for (int i = _height - 1; i >= 0; i--)
541  {
542  for (int j = 0; j < _width; j++)
543  {
544  realdata[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255 + 0.5);
545  imagdata[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)[_IM] - iminVal) / (imaxVal - iminVal) * 255 + 0.5);
546  }
547  }
548 
549  pbitmap->_bitmapinfoheader.dibheadersize = sizeof(bitmapinfoheader);
550  pbitmap->_bitmapinfoheader.width = _width;
551  pbitmap->_bitmapinfoheader.height = _height;
552  pbitmap->_bitmapinfoheader.planes = OPH_PLANES;
553  pbitmap->_bitmapinfoheader.bitsperpixel = bitpixel;
554  pbitmap->_bitmapinfoheader.compression = OPH_COMPRESSION;
555  pbitmap->_bitmapinfoheader.imagesize = _pixelbytesize;
556  pbitmap->_bitmapinfoheader.ypixelpermeter = 0;
557  pbitmap->_bitmapinfoheader.xpixelpermeter = 0;
558  pbitmap->_bitmapinfoheader.numcolorspallette = 256;
559 
560  fwrite(pbitmap, 1, sizeof(bitmap8bit), freal);
561  fwrite(realdata, 1, _pixelbytesize, freal);
562 
563  fwrite(pbitmap, 1, sizeof(bitmap8bit), fimag);
564  fwrite(imagdata, 1, _pixelbytesize, fimag);
565 
566  fclose(freal);
567  fclose(fimag);
568  free(pbitmap);
569  }
570  else
571  {
572  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height * bitpixel / 3);
573  imagdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height * bitpixel / 3);
574  _filesize = _pixelbytesize + sizeof(fileheader) + sizeof(bitmapinfoheader);
575 
576  fileheader *hf = (fileheader*)calloc(1, sizeof(fileheader));
577  bitmapinfoheader *hInfo = (bitmapinfoheader*)calloc(1, sizeof(bitmapinfoheader));
578 
579  hf->signature[0] = 'B';
580  hf->signature[1] = 'M';
581  hf->filesize = _filesize;
582  hf->fileoffset_to_pixelarray = sizeof(fileheader) + sizeof(bitmapinfoheader);
583 
584  double minVal, iminVal, maxVal, imaxVal;
585 
586  for (int z = 0; z < 3; z++)
587  {
588  maxVal = minVal = ComplexH[z](0, 0)[_RE];
589  imaxVal = iminVal = ComplexH[z](0, 0)[_IM];
590 
591  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
592  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
593  if (ComplexH[z](i, j)[_RE] < minVal)
594  minVal = ComplexH[z](i, j)[_RE];
595  if (ComplexH[z](i, j)[_RE] > maxVal)
596  maxVal = ComplexH[z](i, j)[_RE];
597  if (ComplexH[z](i, j)[_IM] < iminVal)
598  iminVal = ComplexH[z](i, j)[_IM];
599  if (ComplexH[z](i, j)[_IM] > imaxVal)
600  imaxVal = ComplexH[z](i, j)[_IM];
601  }
602  }
603 
604  for (int i = _height - 1; i >= 0; i--)
605  {
606  for (int j = 0; j < _width; j++)
607  {
608  realdata[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255);
609  imagdata[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)[_IM] - iminVal) / (imaxVal - iminVal) * 255);
610 
611  }
612  }
613  }
614  hInfo->dibheadersize = sizeof(bitmapinfoheader);
615  hInfo->width = _width;
616  hInfo->height = _height;
617  hInfo->planes = OPH_PLANES;
618  hInfo->bitsperpixel = bitpixel;
619  hInfo->compression = OPH_COMPRESSION;
620  hInfo->imagesize = _pixelbytesize;
621  hInfo->ypixelpermeter = 0;
622  hInfo->xpixelpermeter = 0;
623 
624  fwrite(hf, 1, sizeof(fileheader), freal);
625  fwrite(hInfo, 1, sizeof(bitmapinfoheader), freal);
626  fwrite(realdata, 1, _pixelbytesize, freal);
627 
628  fwrite(hf, 1, sizeof(fileheader), fimag);
629  fwrite(hInfo, 1, sizeof(bitmapinfoheader), fimag);
630  fwrite(imagdata, 1, _pixelbytesize, fimag);
631 
632  fclose(freal);
633  fclose(fimag);
634  free(hf);
635  free(hInfo);
636  }
637 
638  free(realdata);
639  free(imagdata);
640  std::cout << "file save bmp complete\n" << endl;
641 
642  }
643  else {
644  LOG("failed : The Invalid data type! - %s\n", type.c_str());
645  }
646  return true;
647 }
648 
649 bool ophSig::save(const char *fname)
650 {
651  string fullname = fname;
652 
653  char* RGB_name[3] = { 0, };
654  const char* postfix = "";
655 
656  if (_wavelength_num > 1) {
657  postfix = "_B";
658  strcpy(RGB_name[0], postfix);
659  postfix = "_G";
660  strcpy(RGB_name[1], postfix);
661  postfix = "_R";
662  strcpy(RGB_name[2], postfix);
663  }
664  else
665  strcpy(RGB_name[0], postfix);
666 
667  int checktype = static_cast<int>(fullname.rfind("."));
668 
669  if (fullname.substr(checktype + 1, fullname.size()) == "bmp")
670  {
671  oph::uchar* realdata;
672  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * context_.pixel_number[_X] * context_.pixel_number[_Y] * _wavelength_num);
673 
674  double gamma = 0.5;
675  double maxIntensity = 0.0;
676  double realVal = 0.0;
677  double imagVal = 0.0;
678  double intensityVal = 0.0;
679 
680  for (int z = 0; z < _wavelength_num; z++)
681  {
682  for (int j = 0; j < context_.pixel_number[_Y]; j++) {
683  for (int i = 0; i < context_.pixel_number[_X]; i++) {
684  realVal = ComplexH[z](i, j)[_RE];
685  imagVal = ComplexH[z](i, j)[_RE];
686  intensityVal = realVal*realVal + imagVal*imagVal;
687  if (intensityVal > maxIntensity) {
688  maxIntensity = intensityVal;
689  }
690  }
691  }
692  for (int i = context_.pixel_number[_X] - 1; i >= 0; i--)
693  {
694  for (int j = 0; j < context_.pixel_number[_Y]; j++)
695  {
696  realVal = ComplexH[z](context_.pixel_number[_X] - i - 1, j)[_RE];
697  imagVal = ComplexH[z](context_.pixel_number[_X] - i - 1, j)[_IM];
698  intensityVal = realVal*realVal + imagVal*imagVal;
699  realdata[(i*context_.pixel_number[_Y] + j)* _wavelength_num + z] = (uchar)(pow(intensityVal / maxIntensity, gamma)*255.0);
700  }
701  }
702  //sprintf(str, "_%.2u", z);
703  //realname.insert(checktype, RGB_name[z]);
704  }
705  saveAsImg(fullname.c_str(), _wavelength_num * 8, realdata, context_.pixel_number[_X], context_.pixel_number[_Y]);
706 
707  delete[] realdata;
708  }
709  else if (fullname.substr(checktype + 1, fullname.size()) == "bin")
710  {
711  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
712 
713  for (int z = 0; z < _wavelength_num; z++)
714  {
715  fullname = fname;
716  fullname.insert(checktype, RGB_name[z]);
717  std::ofstream cos(fullname.c_str(), std::ios::binary);
718 
719  if (!cos.is_open()) {
720  LOG("Error: file name not found.\n");
721  cos.close();
722  delete[] realdata;
723  return false;
724  }
725 
726  for (int i = 0; i < context_.pixel_number[_X]; i++)
727  {
728  for (int j = 0; j < context_.pixel_number[_Y]; j++)
729  {
730  realdata[context_.pixel_number[_Y] * i + j] = ComplexH[z](i, j)[_RE];
731  }
732  }
733  cos.write(reinterpret_cast<const char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
734 
735  cos.close();
736  }
737  delete[] realdata;
738  }
739 
740  LOG("Writing Openholo Complex Field...%s\n", fullname.c_str());
741  return true;
742 }
743 
744 void ophSig::Data_output(uchar *data, int pos, int bitpixel)
745 {
746 
747  int _width = context_.pixel_number[_Y], _height = context_.pixel_number[_X];
748  OphComplexField* abs_data;
749  abs_data = new OphComplexField[1];
750  abs_data->resize(context_.pixel_number[_Y], _height = context_.pixel_number[_X]);
751  if(pos == 0)
752  {
753  if (bitpixel == 8)
754  {
755 
756  absMat(ComplexH[0], *abs_data);
757  for (int i = _height - 1; i >= 0; i--)
758  {
759  for (int j = 0; j < _width; j++)
760  {
761  if (abs_data[0].mat[_height - i - 1][j][_RE] < 0)
762  {
763  abs_data[0].mat[_height - i - 1][j][_RE] = 0;
764  }
765  }
766  }
767 
768  double minVal, iminVal, maxVal, imaxVal;
769  for (int j = 0; j < abs_data[0].size[_Y]; j++) {
770  for (int i = 0; i < abs_data[0].size[_X]; i++) {
771  if ((i == 0) && (j == 0))
772  {
773  minVal = abs_data[0](i, j)[_RE];
774  maxVal = abs_data[0](i, j)[_RE];
775  }
776  else {
777  if (abs_data[0](i, j)[_RE] < minVal)
778  {
779  minVal = abs_data[0](i, j).real();
780  }
781  if (abs_data[0](i, j)[_RE] > maxVal)
782  {
783  maxVal = abs_data[0](i, j).real();
784  }
785  }
786  }
787  }
788  for (int i = _height - 1; i >= 0; i--)
789  {
790  for (int j = 0; j < _width; j++)
791  {
792  data[i*_width + j] = (uchar)((abs_data[0](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255 + 0.5);
793  }
794  }
795  }
796 
797  else
798  {
799  double minVal, iminVal, maxVal, imaxVal;
800  for (int z = 0; z < 3; z++)
801  {
802  absMat(ComplexH[z], abs_data[0]);
803  for (int j = 0; j < abs_data[0].size[_Y]; j++) {
804  for (int i = 0; i < abs_data[0].size[_X]; i++) {
805  if ((i == 0) && (j == 0))
806  {
807  minVal = abs_data[0](i, j)[_RE];
808  maxVal = abs_data[0](i, j)[_RE];
809  }
810  else {
811  if (abs_data[0](i, j)[_RE] < minVal)
812  {
813  minVal = abs_data[0](i, j)[_RE];
814  }
815  if (abs_data[0](i, j)[_RE] > maxVal)
816  {
817  maxVal = abs_data[0](i, j)[_RE];
818  }
819  }
820  }
821  }
822 
823  for (int i = _height - 1; i >= 0; i--)
824  {
825  for (int j = 0; j < _width; j++)
826  {
827  data[3 * j + 3 * i * _width + z] = (uchar)((abs_data[0](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255);
828  }
829  }
830  }
831  }
832  }
833 
834  else if (pos == 2)
835  {
836  if (bitpixel == 8)
837  {
838  for (int i = _height - 1; i >= 0; i--)
839  {
840  for (int j = 0; j < _width; j++)
841  {
842  if (ComplexH[0].mat[_height - i - 1][j][_RE] < 0)
843  {
844  ComplexH[0].mat[_height - i - 1][j][_RE] = 0;
845  }
846  }
847  }
848 
849  double minVal, iminVal, maxVal, imaxVal;
850  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
851  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
852  if ((i == 0) && (j == 0))
853  {
854  minVal = ComplexH[0](i, j)[_RE];
855  maxVal = ComplexH[0](i, j)[_RE];
856  }
857  else {
858  if (ComplexH[0](i, j)[_RE] < minVal)
859  {
860  minVal = ComplexH[0](i, j).real();
861  }
862  if (ComplexH[0](i, j)[_RE] > maxVal)
863  {
864  maxVal = ComplexH[0](i, j).real();
865  }
866  }
867  }
868  }
869  for (int i = _height - 1; i >= 0; i--)
870  {
871  for (int j = 0; j < _width; j++)
872  {
873  data[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255 + 0.5);
874  }
875  }
876  }
877 
878  else
879  {
880  double minVal, iminVal, maxVal, imaxVal;
881  for (int z = 0; z < 3; z++)
882  {
883  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
884  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
885  if ((i == 0) && (j == 0))
886  {
887  minVal = ComplexH[z](i, j)[_RE];
888  maxVal = ComplexH[z](i, j)[_RE];
889  }
890  else {
891  if (ComplexH[z](i, j)[_RE] < minVal)
892  {
893  minVal = ComplexH[z](i, j)[_RE];
894  }
895  if (ComplexH[z](i, j)[_RE] > maxVal)
896  {
897  maxVal = ComplexH[z](i, j)[_RE];
898  }
899  }
900  }
901  }
902 
903  for (int i = _height - 1; i >= 0; i--)
904  {
905  for (int j = 0; j < _width; j++)
906  {
907  data[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)[_RE] - minVal) / (maxVal - minVal) * 255);
908  }
909  }
910  }
911  }
912  }
913 
914  else if (pos == 1)
915  {
916  if (bitpixel == 8)
917  {
918  for (int i = _height - 1; i >= 0; i--)
919  {
920  for (int j = 0; j < _width; j++)
921  {
922  if (ComplexH[0].mat[_height - i - 1][j][_IM] < 0)
923  {
924  ComplexH[0].mat[_height - i - 1][j][_IM] = 0;
925  }
926  }
927  }
928 
929  double minVal, iminVal, maxVal, imaxVal;
930  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
931  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
932  if ((i == 0) && (j == 0))
933  {
934  minVal = ComplexH[0](i, j)[_IM];
935  maxVal = ComplexH[0](i, j)[_IM];
936  }
937  else {
938  if (ComplexH[0](i, j)[_IM] < minVal)
939  {
940  minVal = ComplexH[0](i, j).imag();
941  }
942  if (ComplexH[0](i, j)[_IM] > maxVal)
943  {
944  maxVal = ComplexH[0](i, j).imag();
945  }
946  }
947  }
948  }
949  for (int i = _height - 1; i >= 0; i--)
950  {
951  for (int j = 0; j < _width; j++)
952  {
953  data[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)[_IM] - minVal) / (maxVal - minVal) * 255 + 0.5);
954  }
955  }
956  }
957 
958  else
959  {
960  double minVal, iminVal, maxVal, imaxVal;
961  for (int z = 0; z < 3; z++)
962  {
963  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
964  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
965  if ((i == 0) && (j == 0))
966  {
967  minVal = ComplexH[z](i, j)[_IM];
968  maxVal = ComplexH[z](i, j)[_IM];
969  }
970  else {
971  if (ComplexH[z](i, j)[_IM] < minVal)
972  {
973  minVal = ComplexH[z](i, j)[_IM];
974  }
975  if (ComplexH[z](i, j)[_IM] > maxVal)
976  {
977  maxVal = ComplexH[z](i, j)[_IM];
978  }
979  }
980  }
981  }
982 
983  for (int i = _height - 1; i >= 0; i--)
984  {
985  for (int j = 0; j < _width; j++)
986  {
987  data[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)[_IM] - minVal) / (maxVal - minVal) * 255);
988  }
989  }
990  }
991  }
992  }
993 }
994 
995 
996 bool ophSig::sigConvertOffaxis(Real angleX, Real angleY) {
997  auto start_time = CUR_TIME;
998 
999  if (is_CPU == true)
1000  {
1001  std::cout << "Start Single Core CPU" << endl;
1002  cvtOffaxis_CPU(angleX,angleY);
1003  }
1004  else {
1005  std::cout << "Start Multi Core GPU" << std::endl;
1006  cvtOffaxis_GPU(angleX, angleY);
1007  }
1008  auto end_time = CUR_TIME;
1009  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1010  LOG("Implement time : %.5lf sec\n", during_time);
1011  return true;
1012 }
1013 
1014 bool ophSig::sigConvertHPO(Real depth, Real_t redRate) {
1015  auto start_time = CUR_TIME;
1016  if (is_CPU == true)
1017  {
1018  std::cout << "Start Single Core CPU" << endl;
1019  sigConvertHPO_CPU(depth,redRate);
1020 
1021  }
1022  else {
1023  std::cout << "Start Multi Core GPU" << std::endl;
1024 
1025  sigConvertHPO_GPU(depth, redRate);
1026 
1027  }
1028 
1029  auto end_time = CUR_TIME;
1030  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1031  LOG("Implement time : %.5lf sec\n", during_time);
1032  return true;
1033 }
1034 
1035 bool ophSig::sigConvertCAC(double red, double green, double blue){
1036  auto start_time = CUR_TIME;
1037  if (is_CPU == true)
1038  {
1039  std::cout << "Start Single Core CPU" << endl;
1040  sigConvertCAC_CPU(red, green, blue);
1041 
1042  }
1043  else {
1044  std::cout << "Start Multi Core GPU" << std::endl;
1045  sigConvertCAC_GPU(red, green, blue);
1046 
1047  }
1048  auto end_time = CUR_TIME;
1049  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1050  LOG("Implement time : %.5lf sec\n", during_time);
1051  return true;
1052 }
1053 
1054 bool ophSig::propagationHolo(float depth) {
1055  auto start_time = CUR_TIME;
1056  if (is_CPU == true)
1057  {
1058  std::cout << "Start Single Core CPU" << endl;
1059  propagationHolo_CPU(depth);
1060 
1061  }
1062  else {
1063  std::cout << "Start Multi Core GPU" << std::endl;
1064  if (_wavelength_num == 1)
1065  {
1066  propagationHolo_GPU(depth);
1067  }
1068  else if (_wavelength_num == 3)
1069  {
1070  Color_propagationHolo_GPU(depth);
1071  }
1072 
1073  }
1074  auto end_time = CUR_TIME;
1075  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1076  LOG("Implement time : %.5lf sec\n", during_time);
1077  return true;
1078 }
1079 
1080 double ophSig::sigGetParamSF(float zMax, float zMin, int sampN, float th) {
1081  auto start_time = CUR_TIME;
1082  double out = 0;
1083  if (is_CPU == true)
1084  {
1085  std::cout << "Start Single Core CPU" << endl;
1086  out = sigGetParamSF_CPU(zMax,zMin,sampN,th);
1087 
1088  }
1089  else {
1090  std::cout << "Start Multi Core GPU" << std::endl;
1091  out = sigGetParamSF_GPU(zMax, zMin, sampN, th);
1092 
1093  }
1094  auto end_time = CUR_TIME;
1095  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1096  LOG("Implement time : %.5lf sec\n", during_time);
1097  return out;
1098 }
1099 
1100 
1101 bool ophSig::cvtOffaxis_CPU(Real angleX, Real angleY) {
1102 
1103  int nx = context_.pixel_number[_X];
1104  int ny = context_.pixel_number[_Y];
1105  OphRealField H1(nx,ny);
1106  Real x, y;
1107  Complex<Real> F;
1108  H1.resize(nx, ny);
1109 
1110  for (int i = 0; i < nx; i++)
1111  {
1112  y = (_cfgSig.height / (nx - 1)*i - _cfgSig.height / 2);
1113  for (int j = 0; j < ny; j++)
1114  {
1115  x = (_cfgSig.width / (ny - 1)*j - _cfgSig.width / 2);
1116 
1117  //(*ComplexH)(i, j)[_RE] = cos(((2 * M_PI) / *context_.wave_length)*(x*sin(angle[_X]) + y *sin(angle[_Y])));
1118  //(*ComplexH)(i, j)[_IM] = sin(((2 * M_PI) / *context_.wave_length)*(x*sin(angle[_X]) + y *sin(angle[_Y])));
1119  F[_RE] = cos(((2 * M_PI) / *context_.wave_length)*(x*sin(angleX) + y *sin(angleY)));
1120  F[_IM] = sin(((2 * M_PI) / *context_.wave_length)*(x*sin(angleX) + y *sin(angleY)));
1121  H1(i, j) = ((*ComplexH)(i, j) * F)[_RE];
1122  }
1123  }
1124  double out = minOfMat(H1);
1125  H1 - out;
1126  out = maxOfMat(H1);
1127  H1 / out;
1128  //normalizeMat(H1, H1);
1129 
1130 
1131  for (int i = 0; i < nx; i++)
1132  {
1133  for (int j = 0; j < ny; j++)
1134  {
1135  (*ComplexH)(i, j)[_RE] = H1(i, j);
1136  (*ComplexH)(i, j)[_IM] = 0;
1137  }
1138  }
1139 
1140 
1141 
1142  return true;
1143 }
1144 
1145 bool ophSig::sigConvertHPO_CPU(Real depth, Real_t redRate) {
1146 
1147 
1148  int nx = context_.pixel_number[_X];
1149  int ny = context_.pixel_number[_Y];
1150 
1151  Real wl = *context_.wave_length;
1152  Real NA = _cfgSig.width/(2*depth);
1153 
1154  int xshift = nx / 2;
1155  int yshift = ny / 2;
1156 
1157  Real y;
1158 
1159  Real_t NA_g = NA * redRate;
1160 
1161  OphComplexField F1(nx, ny);
1162  OphComplexField FH(nx, ny);
1163 
1164 
1165  Real Rephase = -(1 / (4 * M_PI)*pow((wl / NA_g), 2));
1166  Real Imphase = ((1 / (4 * M_PI))*depth*wl);
1167 
1168  for (int i = 0; i < ny; i++)
1169  {
1170  int ii = (i + yshift) % ny;
1171 
1172  for (int j = 0; j < nx; j++)
1173  {
1174  y = (2 * M_PI * (j) / _cfgSig.height - M_PI * (nx - 1) / _cfgSig.height);
1175  int jj = (j + xshift) % nx;
1176  F1(jj, ii)[_RE] = std::exp(Rephase*pow(y, 2))*cos(Imphase*pow(y, 2));
1177  F1(jj, ii)[_IM] = std::exp(Rephase*pow(y, 2))*sin(Imphase*pow(y, 2));
1178  }
1179  }
1180  fft2((*ComplexH), FH, OPH_FORWARD);
1181  F1.mulElem(FH);
1182  fft2(F1, (*ComplexH), OPH_BACKWARD);
1183 
1184 
1185 
1186 
1187  return true;
1188 
1189 }
1190 
1191 bool ophSig::sigConvertCAC_CPU(double red, double green, double blue) {
1192 
1193  Real x, y;
1194  //OphComplexField exp, conj, FH_CAC;
1195  int nx = context_.pixel_number[_X];
1196  int ny = context_.pixel_number[_Y];
1197 
1198  if (_wavelength_num != 3) {
1199  _wavelength_num = 3;
1200  delete[] context_.wave_length;
1201  context_.wave_length = new Real[_wavelength_num];
1202  }
1203 
1204  context_.wave_length[0] = blue;
1205  context_.wave_length[1] = green;
1206  context_.wave_length[2] = red;
1207 
1208  OphComplexField FFZP(nx, ny);
1209  OphComplexField FH(nx, ny);
1210 
1211  for (int z = 0; z < _wavelength_num; z++)
1212  {
1213  double sigmaf = ((_foc[2] - _foc[z]) * context_.wave_length[z]) / (4 * M_PI);
1214  int xshift = nx / 2;
1215  int yshift = ny / 2;
1216 
1217  for (int i = 0; i < ny; i++)
1218  {
1219  int ii = (i + yshift) % ny;
1220  y = (2 * M_PI * i) / _radius - (M_PI*(ny - 1)) / _radius;
1221  for (int j = 0; j < nx; j++)
1222  {
1223  x = (2 * M_PI * j) / _radius - (M_PI*(nx - 1)) / _radius;
1224 
1225  int jj = (j + xshift) % nx;
1226 
1227  FFZP(jj, ii)[_RE] = cos(sigmaf * (pow(x, 2) + pow(y, 2)));
1228  FFZP(jj, ii)[_IM] = -sin(sigmaf * (pow(x, 2) + pow(y, 2))); //conjugate ¶§¹®¿¡ -ºÙÀ½¿©
1229  }
1230  }
1231  fft2(ComplexH[z], FH, OPH_FORWARD);
1232  FH.mulElem(FFZP);
1233  fft2(FH, ComplexH[z], OPH_BACKWARD);
1234  }
1235 
1236  return true;
1237 }
1238 void ophSig::Parameter_Set(int nx, int ny, double width, double height, double NA)
1239 {
1240  context_.pixel_number[_X] = nx;
1241  context_.pixel_number[_Y] = ny;
1242  _cfgSig.width = width;
1243  _cfgSig.height = height;
1244  _cfgSig.NA = NA;
1245 }
1246 
1247 void ophSig::wavelength_Set(double wavelength)
1248 {
1249  *context_.wave_length = wavelength;
1250 }
1251 
1252 void ophSig::focal_length_Set(double red, double green, double blue,double rad)
1253 {
1254  _foc[2] = red;
1255  _foc[1] = green;
1256  _foc[0] = blue;
1257  _radius = rad;
1258 }
1259 
1260 void ophSig::Wavenumber_output(int &wavenumber)
1261 {
1262  wavenumber = _wavelength_num;
1263 }
1264 
1265 bool ophSig::readConfig(const char* fname)
1266 {
1267  //LOG("Reading....%s...\n", fname);
1268 
1269  /*XML parsing*/
1270  tinyxml2::XMLDocument xml_doc;
1271  tinyxml2::XMLNode* xml_node;
1272 
1273  if (!checkExtension(fname, ".xml"))
1274  {
1275  LOG("file's extension is not 'xml'\n");
1276  return false;
1277  }
1278  auto ret = xml_doc.LoadFile(fname);
1279  if (ret != tinyxml2::XML_SUCCESS)
1280  {
1281  LOG("Failed to load file \"%s\"\n", fname);
1282  return false;
1283  }
1284 
1285  xml_node = xml_doc.FirstChild();
1286 
1287  (xml_node->FirstChildElement("pixel_number_x"))->QueryIntText(&context_.pixel_number[_X]);
1288  (xml_node->FirstChildElement("pixel_number_y"))->QueryIntText(&context_.pixel_number[_Y]);
1289  (xml_node->FirstChildElement("width"))->QueryFloatText(&_cfgSig.width);
1290  (xml_node->FirstChildElement("height"))->QueryFloatText(&_cfgSig.height);
1291  (xml_node->FirstChildElement("wavelength_num"))->QueryIntText(&_wavelength_num);
1292 
1293  context_.wave_length = new Real[_wavelength_num];
1294 
1295  (xml_node->FirstChildElement("wavelength"))->QueryDoubleText(context_.wave_length);
1296  (xml_node->FirstChildElement("NA"))->QueryFloatText(&_cfgSig.NA);
1297  (xml_node->FirstChildElement("z"))->QueryFloatText(&_cfgSig.z);
1298  (xml_node->FirstChildElement("radius_of_lens"))->QueryFloatText(&_radius);
1299  (xml_node->FirstChildElement("focal_length_R"))->QueryFloatText(&_foc[2]);
1300  (xml_node->FirstChildElement("focal_length_G"))->QueryFloatText(&_foc[1]);
1301  (xml_node->FirstChildElement("focal_length_B"))->QueryFloatText(&_foc[0]);
1302  // 2024.04.23. mwnam
1303 // set variable for resolution
1304  resCfg = context_.pixel_number;
1305 
1306  return true;
1307 }
1308 
1309 
1310 
1311 bool ophSig::propagationHolo_CPU(float depth) {
1312  int i, j;
1313  Real x, y, sigmaf;
1314  int nx = context_.pixel_number[_X];
1315  int ny = context_.pixel_number[_Y];
1316 
1317  OphComplexField FH(nx, ny);
1318  //OphComplexField FFZP(nx, ny);
1319  int xshift = nx / 2;
1320  int yshift = ny / 2;
1321 
1322  for (int z = 0; z < _wavelength_num; z++) {
1323 
1324  sigmaf = (depth * context_.wave_length[z]) / (4 * M_PI);
1325 
1326  /*FH.resize(nx, ny);
1327  FFZP.resize(nx, ny);*/
1328 
1329  fft2(ComplexH[z], FH, OPH_FORWARD);
1330 
1331  for (i = 0; i < ny; i++)
1332  {
1333  int ii = (i + yshift) % ny;
1334  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1335 
1336  for (j = 0; j < nx; j++)
1337  {
1338  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1339  int jj = (j + xshift) % nx;
1340  double temp = FH(jj, ii)[_RE];
1341  FH(jj, ii)[_RE] = cos(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_RE] - sin(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_IM];
1342  FH(jj, ii)[_IM] = sin(sigmaf * (pow(x, 2) + pow(y, 2))) * temp + cos(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_IM];
1343 
1344  }
1345  }
1346 
1347  fft2(FH, ComplexH[z], OPH_BACKWARD);
1348  }
1349  return true;
1350 }
1351 
1352 OphComplexField ophSig::propagationHolo(OphComplexField complexH, float depth) {
1353  int i, j;
1354  Real x, y, sigmaf;
1355  int nx = context_.pixel_number[_X];
1356  int ny = context_.pixel_number[_Y];
1357 
1358  OphComplexField FH(nx, ny);
1359  int xshift = nx / 2;
1360  int yshift = ny / 2;
1361 
1362 
1363  sigmaf = (depth * (*context_.wave_length)) / (4 * M_PI);
1364 
1365  fft2(complexH, FH);
1366 
1367  for (i = 0; i < ny; i++)
1368  {
1369  int ii = (i + yshift) % ny;
1370  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1371 
1372  for (j = 0; j < nx; j++)
1373  {
1374  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1375  int jj = (j + xshift) % nx;
1376  double temp = FH(jj, ii)[_RE];
1377  FH(jj, ii)[_RE] = cos(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_RE] - sin(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_IM];
1378  FH(jj, ii)[_IM] = sin(sigmaf * (pow(x, 2) + pow(y, 2))) * temp + cos(sigmaf * (pow(x, 2) + pow(y, 2))) * FH(jj, ii)[_IM];
1379 
1380  }
1381  }
1382  fft2(FH, complexH, OPH_BACKWARD);
1383 
1384  return complexH;
1385 }
1386 
1387 double ophSig::sigGetParamAT()
1388 {
1389  //auto start_time = CUR_TIME;
1390 
1391  Real index = 0;
1392  if (is_CPU == true)
1393  {
1394  //std::cout << "Start Single Core CPU" << endl;
1395  index = sigGetParamAT_CPU();
1396 
1397  }
1398  else {
1399  //std::cout << "Start Multi Core GPU" << std::endl;
1400  index = sigGetParamAT_GPU();
1401 
1402  }
1403 
1404  //auto end_time = CUR_TIME;
1405 
1406  //auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1407 
1408  //LOG("Implement time : %.5lf sec\n", during_time);
1409  return index;
1410 }
1411 
1412 double ophSig::sigGetParamAT_CPU() {
1413 
1414  int i = 0, j = 0;
1415  Real max = 0; Real index = 0;
1416  Real_t NA_g = (Real_t)0.025;
1417  int nx = context_.pixel_number[_X];
1418  int ny = context_.pixel_number[_Y];
1419 
1420  OphComplexField Flr(nx, ny);
1421  OphComplexField Fli(nx, ny);
1422  OphComplexField Hsyn(nx, ny);
1423 
1424  OphComplexField Fo(nx, ny);
1425  OphComplexField Fon, yn, Ab_yn;
1426 
1427  OphRealField Ab_yn_half;
1428  OphRealField G(nx, ny);
1429  Real x = 1, y = 1;
1430  vector<Real> t, tn;
1431 
1432  for (i = 0; i < nx; i++)
1433  {
1434  for (j = 0; j < ny; j++)
1435  {
1436 
1437  x = (2 * M_PI*(i) / _cfgSig.height - M_PI*(nx - 1) / _cfgSig.height);
1438  y = (2 * M_PI*(j) / _cfgSig.width - M_PI*(ny - 1) / _cfgSig.width);
1439  G(i, j) = std::exp(-M_PI * pow((*context_.wave_length) / (2 * M_PI * NA_g), 2) * (pow(y, 2) + pow(x, 2)));
1440  Flr(i, j)[_RE] = (*ComplexH)(i, j)[_RE];
1441  Fli(i, j)[_RE] = (*ComplexH)(i, j)[_IM];
1442  Flr(i, j)[_IM] = 0;
1443  Fli(i, j)[_IM] = 0;
1444  }
1445  }
1446 
1447  fft2(Flr, Flr);
1448  fft2(Fli, Fli);
1449 
1450  int xshift = nx / 2;
1451  int yshift = ny / 2;
1452 
1453  for (i = 0; i < nx; i++)
1454  {
1455  int ii = (i + xshift) % nx;
1456  for (j = 0; j < ny; j++)
1457  {
1458  int jj = (j + yshift) % ny;
1459  Hsyn(i, j)[_RE] = Flr(i, j)[_RE] * G(i, j);
1460  Hsyn(i, j)[_IM] = Fli(i, j)[_RE] * G(i, j);
1461  /*Hsyn_copy1(i, j) = Hsyn(i, j);
1462  Hsyn_copy2(i, j) = Hsyn_copy1(i, j) * Hsyn(i, j);
1463  Hsyn_copy3(i, j) = pow(sqrt(Hsyn(i, j)[_RE] * Hsyn(i, j)[_RE] + Hsyn(i, j)[_IM] * Hsyn(i, j)[_IM]), 2) + pow(10, -300);
1464  Fo(ii, jj)[_RE] = Hsyn_copy2(i, j)[0] / Hsyn_copy3(i, j);
1465  Fo(ii, jj)[_IM] = Hsyn_copy2(i, j)[1] / Hsyn_copy3(i, j);*/
1466  Fo(ii, jj) = pow(Hsyn(i, j), 2) / (pow(abs(Hsyn(i, j)), 2) + pow(10, -300));
1467 
1468  }
1469  }
1470 
1471  t = linspace(0., 1., nx / 2 + 1);
1472  tn.resize(t.size());
1473  Fon.resize(1, t.size());
1474 
1475  int size = (int)tn.size();
1476  for (int i = 0; i < size; i++)
1477  {
1478  tn.at(i) = pow(t.at(i), 0.5);
1479  Fon(0, i)[_RE] = Fo(nx / 2 - 1, nx / 2 - 1 + i)[_RE];
1480  Fon(0, i)[_IM] = 0;
1481  }
1482  yn.resize(1, size);
1483  linInterp(t, Fon, tn, yn);
1484  fft1(yn, yn);
1485  Ab_yn.resize(yn.size[_X], yn.size[_Y]);
1486  absMat(yn, Ab_yn);
1487  Ab_yn_half.resize(1, nx / 4 + 1);
1488 
1489  for (int i = 0; i < nx / 4 + 1; i++)
1490  {
1491  Ab_yn_half(0, i) = Ab_yn(0, nx / 4 + i - 1)[_RE];
1492  if (i == 0) max = Ab_yn_half(0, 0);
1493  else
1494  {
1495  if (Ab_yn_half(0, i) > max)
1496  {
1497  max = Ab_yn_half(0, i);
1498  index = i;
1499  }
1500  }
1501  }
1502 
1503  index = -(((index + 1) - 120) / 10) / 140 + 0.1;
1504 
1505 
1506 
1507  return index;
1508 }
1509 
1510 double ophSig::sigGetParamSF_CPU(float zMax, float zMin, int sampN, float th)
1511 {
1512  int nx = context_.pixel_number[_X];
1513  int ny = context_.pixel_number[_Y];
1514 
1515  OphComplexField I(nx, ny);
1516  vector<Real> F;
1517  Real dz = (zMax - zMin) / sampN;
1518  Real f;
1519  Real_t z = 0;
1520  Real depth = 0;
1521  Real max = MIN_DOUBLE;
1522  int i, j, n = 0;
1523  Real ret1;
1524  Real ret2;
1525 
1526 
1527 
1528  for (n = 0; n < sampN + 1; n++)
1529  {
1530  OphComplexField field = *(ComplexH);
1531  z = ((n)* dz + zMin);
1532  f = 0;
1533  I = propagationHolo(field, z);
1534 
1535  for (i = 0; i < nx - 2; i++)
1536  {
1537  for (j = 0; j < ny - 2; j++)
1538  {
1539  ret1 = abs(I(i + 2, j)[_RE] - I(i, j)[_RE]);
1540  ret2 = abs(I(i, j + 2)[_RE] - I(i, j)[_RE]);
1541  if (ret1 >= th) { f += ret1 * ret1; }
1542  else if (ret2 >= th) { f += ret2 * ret2; }
1543  }
1544  }
1545 
1546  if (f > max) {
1547  max = f;
1548  depth = z;
1549  }
1550  }
1551 
1552 
1553  return depth;
1554 }
1555 
1556 bool ophSig::getComplexHFromPSDH(const char * fname0, const char * fname90, const char * fname180, const char * fname270)
1557 {
1558  auto start_time = CUR_TIME;
1559  string fname0str = fname0;
1560  string fname90str = fname90;
1561  string fname180str = fname180;
1562  string fname270str = fname270;
1563  int checktype = static_cast<int>(fname0str.rfind("."));
1564  OphRealField f0Mat[3], f90Mat[3], f180Mat[3], f270Mat[3];
1565 
1566  std::string f0type = fname0str.substr(checktype + 1, fname0str.size());
1567 
1568  uint16_t bitsperpixel;
1569  fileheader hf;
1570  bitmapinfoheader hInfo;
1571 
1572  if (f0type == "bmp")
1573  {
1574  FILE *f0, *f90, *f180, *f270;
1575  f0 = fopen(fname0str.c_str(), "rb");
1576  f90 = fopen(fname90str.c_str(), "rb");
1577  f180 = fopen(fname180str.c_str(), "rb");
1578  f270 = fopen(fname270str.c_str(), "rb");
1579  if (!f0)
1580  {
1581  LOG("bmp file open fail! (phase shift = 0)\n");
1582  return false;
1583  }
1584  if (!f90)
1585  {
1586  LOG("bmp file open fail! (phase shift = 90)\n");
1587  return false;
1588  }
1589  if (!f180)
1590  {
1591  LOG("bmp file open fail! (phase shift = 180)\n");
1592  return false;
1593  }
1594  if (!f270)
1595  {
1596  LOG("bmp file open fail! (phase shift = 270)\n");
1597  return false;
1598  }
1599  size_t nRead = fread(&hf, sizeof(fileheader), 1, f0);
1600  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f0);
1601 
1602  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
1603  if ((hInfo.height == 0) || (hInfo.width == 0))
1604  {
1605  LOG("bmp header is empty!\n");
1606  hInfo.height = context_.pixel_number[_X];
1607  hInfo.width = context_.pixel_number[_Y];
1608  if (hInfo.height == 0 || hInfo.width == 0)
1609  {
1610  LOG("check your parameter file!\n");
1611  return false;
1612  }
1613  }
1614  if ((context_.pixel_number[_Y] != (int)hInfo.height) || (context_.pixel_number[_X] != (int)hInfo.width)) {
1615  LOG("image size is different!\n");
1616  context_.pixel_number[_Y] = (int)hInfo.height;
1617  context_.pixel_number[_X] = (int)hInfo.width;
1618  LOG("changed parameter of size %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
1619  }
1620  bitsperpixel = hInfo.bitsperpixel;
1621  if (hInfo.bitsperpixel == 8)
1622  {
1623  _wavelength_num = 1;
1624  rgbquad palette[256];
1625  size_t nRead;
1626  nRead = fread(palette, sizeof(rgbquad), 256, f0);
1627  nRead = fread(palette, sizeof(rgbquad), 256, f90);
1628  nRead = fread(palette, sizeof(rgbquad), 256, f180);
1629  nRead = fread(palette, sizeof(rgbquad), 256, f270);
1630 
1631  f0Mat[0].resize(hInfo.height, hInfo.width);
1632  f90Mat[0].resize(hInfo.height, hInfo.width);
1633  f180Mat[0].resize(hInfo.height, hInfo.width);
1634  f270Mat[0].resize(hInfo.height, hInfo.width);
1635  ComplexH = new OphComplexField;
1636  ComplexH[0].resize(hInfo.height, hInfo.width);
1637  }
1638  else
1639  {
1640  _wavelength_num = 3;
1641  ComplexH = new OphComplexField[3];
1642  f0Mat[0].resize(hInfo.height, hInfo.width);
1643  f90Mat[0].resize(hInfo.height, hInfo.width);
1644  f180Mat[0].resize(hInfo.height, hInfo.width);
1645  f270Mat[0].resize(hInfo.height, hInfo.width);
1646  ComplexH[0].resize(hInfo.height, hInfo.width);
1647 
1648  f0Mat[1].resize(hInfo.height, hInfo.width);
1649  f90Mat[1].resize(hInfo.height, hInfo.width);
1650  f180Mat[1].resize(hInfo.height, hInfo.width);
1651  f270Mat[1].resize(hInfo.height, hInfo.width);
1652  ComplexH[1].resize(hInfo.height, hInfo.width);
1653 
1654  f0Mat[2].resize(hInfo.height, hInfo.width);
1655  f90Mat[2].resize(hInfo.height, hInfo.width);
1656  f180Mat[2].resize(hInfo.height, hInfo.width);
1657  f270Mat[2].resize(hInfo.height, hInfo.width);
1658  ComplexH[2].resize(hInfo.height, hInfo.width);
1659  }
1660 
1661  size_t size = hInfo.width * hInfo.height * (hInfo.bitsperpixel >> 3);
1662  uchar* f0data = (uchar*)malloc(sizeof(uchar) * size);
1663  uchar* f90data = (uchar*)malloc(sizeof(uchar) * size);
1664  uchar* f180data = (uchar*)malloc(sizeof(uchar) * size);
1665  uchar* f270data = (uchar*)malloc(sizeof(uchar) * size);
1666 
1667  nRead = fread(f0data, sizeof(uchar), size, f0);
1668  nRead = fread(f90data, sizeof(uchar), size, f90);
1669  nRead = fread(f180data, sizeof(uchar), size, f180);
1670  nRead = fread(f270data, sizeof(uchar), size, f270);
1671 
1672  fclose(f0);
1673  fclose(f90);
1674  fclose(f180);
1675  fclose(f270);
1676 
1677  uint16_t bpp = hInfo.bitsperpixel >> 3;
1678  uint32_t nLine = hInfo.width * bpp;
1679  for (int i = hInfo.height - 1; i >= 0; i--)
1680  {
1681  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
1682  {
1683  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1684  {
1685  int idx = hInfo.height - i - 1;
1686  int idx2 = i * nLine + bpp * j + z;
1687  f0Mat[z](idx, j) = (double)f0data[idx2];
1688  f90Mat[z](idx, j) = (double)f90data[idx2];
1689  f180Mat[z](idx, j) = (double)f180data[idx2];
1690  f270Mat[z](idx, j) = (double)f270data[idx2];
1691  }
1692  }
1693  }
1694  LOG("PSDH file load complete!\n");
1695 
1696  free(f0data);
1697  free(f90data);
1698  free(f180data);
1699  free(f270data);
1700 
1701  }
1702  else
1703  {
1704  LOG("wrong type (only BMP supported)\n");
1705  }
1706 
1707  // calculation complexH from 4 psdh and then normalize
1708  double normalizefactor = 1. / 256.;
1709  for (int z = 0; z < (hInfo.bitsperpixel >> 3); z++)
1710  {
1711  for (int i = 0; i < context_.pixel_number[_X]; i++)
1712  {
1713  for (int j = 0; j < context_.pixel_number[_Y]; j++)
1714  {
1715  ComplexH[z][j][i][_RE] = (f0Mat[z][j][i] - f180Mat[z][j][i])*normalizefactor;
1716  ComplexH[z][j][i][_IM] = (f90Mat[z][j][i] - f270Mat[z][j][i])*normalizefactor;
1717 
1718  }
1719  }
1720  }
1721  LOG("complex field obtained from 4 psdh\n");
1722 
1723  auto end_time = CUR_TIME;
1724 
1725  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1726 
1727  LOG("Implement time : %.5lf sec\n", during_time);
1728 
1729  return true;
1730 }
1731 
1732 bool ophSig::getComplexHFrom3ArbStepPSDH(
1733  const char* fname0, const char* fname1,
1734  const char* fname2, const char* fnameOI,
1735  const char* fnameRI, int nIter)
1736 {
1737  auto start_time = CUR_TIME;
1738  string fname0str = fname0;
1739  string fname1str = fname1;
1740  string fname2str = fname2;
1741  string fnameOIstr = fnameOI;
1742  string fnameRIstr = fnameRI;
1743  int checktype = static_cast<int>(fname0str.rfind("."));
1744  OphRealField f0Mat[3], f1Mat[3], f2Mat[3], fOIMat[3], fRIMat[3];
1745 
1746  std::string f0type = fname0str.substr(checktype + 1, fname0str.size());
1747 
1748  uint16_t bitsperpixel;
1749  fileheader hf;
1750  bitmapinfoheader hInfo;
1751 
1752  if (f0type == "bmp")
1753  {
1754  FILE *f0, *f1, *f2, *fOI, *fRI;
1755  f0 = fopen(fname0str.c_str(), "rb");
1756  f1 = fopen(fname1str.c_str(), "rb");
1757  f2 = fopen(fname2str.c_str(), "rb");
1758  fOI = fopen(fnameOIstr.c_str(), "rb");
1759  fRI = fopen(fnameRIstr.c_str(), "rb");
1760  if (!f0)
1761  {
1762  LOG("bmp file open fail! (first interference pattern)\n");
1763  return false;
1764  }
1765  if (!f1)
1766  {
1767  LOG("bmp file open fail! (second interference pattern)\n");
1768  return false;
1769  }
1770  if (!f2)
1771  {
1772  LOG("bmp file open fail! (third interference pattern)\n");
1773  return false;
1774  }
1775  if (!fOI)
1776  {
1777  LOG("bmp file open fail! (object wave intensity pattern)\n");
1778  return false;
1779  }
1780  if (!fRI)
1781  {
1782  LOG("bmp file open fail! (reference wave intensity pattern)\n");
1783  return false;
1784  }
1785 
1786  size_t nRead;
1787  nRead = fread(&hf, sizeof(fileheader), 1, f0);
1788  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f0);
1789  nRead = fread(&hf, sizeof(fileheader), 1, f1);
1790  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f1);
1791  nRead = fread(&hf, sizeof(fileheader), 1, f2);
1792  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f2);
1793  nRead = fread(&hf, sizeof(fileheader), 1, fOI);
1794  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, fOI);
1795  nRead = fread(&hf, sizeof(fileheader), 1, fRI);
1796  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, fRI);
1797 
1798  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
1799  if ((hInfo.height == 0) || (hInfo.width == 0))
1800  {
1801  LOG("bmp header is empty!\n");
1802  hInfo.height = context_.pixel_number[_X];
1803  hInfo.width = context_.pixel_number[_Y];
1804  if (hInfo.height == 0 || hInfo.width == 0)
1805  {
1806  LOG("check your parameter file!\n");
1807  return false;
1808  }
1809  }
1810  if ((context_.pixel_number[_Y] != (int)hInfo.height) || (context_.pixel_number[_X] != (int)hInfo.width)) {
1811  LOG("image size is different!\n");
1812  context_.pixel_number[_Y] = (int)hInfo.height;
1813  context_.pixel_number[_X] = (int)hInfo.width;
1814  LOG("changed parameter of size %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
1815  }
1816  bitsperpixel = hInfo.bitsperpixel;
1817  if (hInfo.bitsperpixel == 8)
1818  {
1819  _wavelength_num = 1;
1820  rgbquad palette[256];
1821 
1822  size_t nRead;
1823  nRead = fread(palette, sizeof(rgbquad), 256, f0);
1824  nRead = fread(palette, sizeof(rgbquad), 256, f1);
1825  nRead = fread(palette, sizeof(rgbquad), 256, f2);
1826  nRead = fread(palette, sizeof(rgbquad), 256, fOI);
1827  nRead = fread(palette, sizeof(rgbquad), 256, fRI);
1828 
1829  f0Mat[0].resize(hInfo.height, hInfo.width);
1830  f1Mat[0].resize(hInfo.height, hInfo.width);
1831  f2Mat[0].resize(hInfo.height, hInfo.width);
1832  fOIMat[0].resize(hInfo.height, hInfo.width);
1833  fRIMat[0].resize(hInfo.height, hInfo.width);
1834  ComplexH = new OphComplexField;
1835  ComplexH[0].resize(hInfo.height, hInfo.width);
1836  }
1837  else
1838  {
1839  _wavelength_num = 3;
1840  ComplexH = new OphComplexField[3];
1841  f0Mat[0].resize(hInfo.height, hInfo.width);
1842  f1Mat[0].resize(hInfo.height, hInfo.width);
1843  f2Mat[0].resize(hInfo.height, hInfo.width);
1844  fOIMat[0].resize(hInfo.height, hInfo.width);
1845  fRIMat[0].resize(hInfo.height, hInfo.width);
1846  ComplexH[0].resize(hInfo.height, hInfo.width);
1847 
1848  f0Mat[1].resize(hInfo.height, hInfo.width);
1849  f1Mat[1].resize(hInfo.height, hInfo.width);
1850  f2Mat[1].resize(hInfo.height, hInfo.width);
1851  fOIMat[1].resize(hInfo.height, hInfo.width);
1852  fRIMat[1].resize(hInfo.height, hInfo.width);
1853  ComplexH[1].resize(hInfo.height, hInfo.width);
1854 
1855  f0Mat[2].resize(hInfo.height, hInfo.width);
1856  f1Mat[2].resize(hInfo.height, hInfo.width);
1857  f2Mat[2].resize(hInfo.height, hInfo.width);
1858  fOIMat[2].resize(hInfo.height, hInfo.width);
1859  fRIMat[2].resize(hInfo.height, hInfo.width);
1860  ComplexH[2].resize(hInfo.height, hInfo.width);
1861  }
1862 
1863  size_t size = hInfo.width * hInfo.height * (hInfo.bitsperpixel >> 3);
1864 
1865  uchar* f0data = (uchar*)malloc(sizeof(uchar) * size);
1866  uchar* f1data = (uchar*)malloc(sizeof(uchar) * size);
1867  uchar* f2data = (uchar*)malloc(sizeof(uchar) * size);
1868  uchar* fOIdata = (uchar*)malloc(sizeof(uchar) * size);
1869  uchar* fRIdata = (uchar*)malloc(sizeof(uchar) * size);
1870 
1871  nRead = fread(f0data, sizeof(uchar), size, f0);
1872  nRead = fread(f1data, sizeof(uchar), size, f1);
1873  nRead = fread(f2data, sizeof(uchar), size, f2);
1874  nRead = fread(fOIdata, sizeof(uchar), size, fOI);
1875  nRead = fread(fRIdata, sizeof(uchar), size, fRI);
1876 
1877  fclose(f0);
1878  fclose(f1);
1879  fclose(f2);
1880  fclose(fOI);
1881  fclose(fRI);
1882 
1883  uint16_t bpp = hInfo.bitsperpixel >> 3;
1884 
1885  for (int i = hInfo.height - 1; i >= 0; i--)
1886  {
1887  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
1888  {
1889  for (int z = 0; z < (hInfo.bitsperpixel >> 3); z++)
1890  {
1891  int idx = hInfo.height - i - 1;
1892  size_t idx2 = i * size + bpp * j + z;
1893  f0Mat[z](idx, j) = (double)f0data[idx2];
1894  f1Mat[z](idx, j) = (double)f1data[idx2];
1895  f2Mat[z](idx, j) = (double)f2data[idx2];
1896  fOIMat[z](idx, j) = (double)fOIdata[idx2];
1897  fRIMat[z](idx, j) = (double)fRIdata[idx2];
1898  }
1899  }
1900  }
1901 
1902  LOG("PSDH_3ArbStep file load complete!\n");
1903 
1904  free(f0data);
1905  free(f1data);
1906  free(f2data);
1907  free(fOIdata);
1908  free(fRIdata);
1909 
1910  }
1911  else
1912  {
1913  LOG("wrong type (only BMP supported)\n");
1914  }
1915 
1916  // calculation complexH from 3 arbitrary step intereference patterns and the object wave intensity
1917  // prepare some variables
1918  double P[2] = { 0.0, }; // please see ref.
1919  double C[2] = { 2.0/M_PI, 2.0/M_PI };
1920  double alpha[2] = { 0.0, }; //phaseShift[j+1]-phaseShift[j]
1921  double ps[3] = { 0.0, }; // reference wave phase shift for each inteference pattern
1922  const int nX = context_.pixel_number[_X];
1923  const int nY = context_.pixel_number[_Y];
1924  const int nXY = nX * nY;
1925 
1926 
1927  // calculate difference between interference patterns
1928  OphRealField I01Mat, I02Mat, I12Mat, OAMat, RAMat;
1929  I01Mat.resize(nY, nX);
1930  I02Mat.resize(nY, nX);
1931  I12Mat.resize(nY, nX);
1932  OAMat.resize(nY, nX);
1933  RAMat.resize(nY, nX);
1934 
1935  double sin2m1h, sin2m0h, sin1m0h, sin0p1h, sin0p2h, cos0p1h, cos0p2h, sin1p2h, cos1p2h;
1936  double sinP, cosP;
1937  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1938  {
1939  // initialize
1940  P[0] = 0.0;
1941  P[1] = 0.0;
1942  C[0] = 2.0 / M_PI;
1943  C[1] = 2.0 / M_PI;
1944 
1945  // load current channel
1946  for (int i = 0; i < nX; i++)
1947  {
1948  for (int j = 0; j < nY; j++)
1949  {
1950  I01Mat[j][i] = (f0Mat[z][j][i] - f1Mat[z][j][i]) / 255.; // difference & normalize
1951  I02Mat[j][i] = (f0Mat[z][j][i] - f2Mat[z][j][i]) / 255.; // difference & normalize
1952  I12Mat[j][i] = (f1Mat[z][j][i] - f2Mat[z][j][i]) / 255.; // difference & normalize
1953  OAMat[j][i] = sqrt(fOIMat[z][j][i] / 255.); // normalize & then calculate amplitude from intensity
1954  RAMat[j][i] = sqrt(fRIMat[z][j][i] / 255.); // normalize & then calculate amplitude from intensity
1955  }
1956  }
1957 
1958  // calculate P
1959  for (int i = 0; i < nX; i++)
1960  {
1961  for (int j = 0; j < nY; j++)
1962  {
1963  P[0] += abs(I01Mat[j][i] / OAMat[j][i] / RAMat[j][i]);
1964  P[1] += abs(I12Mat[j][i] / OAMat[j][i] / RAMat[j][i]);
1965  }
1966  }
1967  P[0] = P[0] / (4.*((double) nXY));
1968  P[1] = P[1] / (4.*((double) nXY));
1969  LOG("P %f %f\n", P[0], P[1]);
1970 
1971  // iterative search
1972  for (int iter = 0; iter < nIter; iter++)
1973  {
1974  LOG("C %d %f %f\n", iter, C[0], C[1]);
1975  LOG("ps %d %f %f %f\n", iter, ps[0], ps[1], ps[2]);
1976 
1977  alpha[0] = 2.*asin(P[0] / C[0]);
1978  alpha[1] = 2.*asin(P[1] / C[1]);
1979 
1980  ps[0] = 0.0;
1981  ps[1] = ps[0] + alpha[0];
1982  ps[2] = ps[1] + alpha[1];
1983 
1984  sin2m1h = sin((ps[2] - ps[1]) / 2.);
1985  sin2m0h = sin((ps[2] - ps[0]) / 2.);
1986  sin1m0h = sin((ps[1] - ps[0]) / 2.);
1987  sin0p1h = sin((ps[0] + ps[1]) / 2.);
1988  sin0p2h = sin((ps[0] + ps[2]) / 2.);
1989  cos0p1h = cos((ps[0] + ps[1]) / 2.);
1990  cos0p2h = cos((ps[0] + ps[2]) / 2.);
1991  for (int i = 0; i < nX; i++)
1992  {
1993  for (int j = 0; j < nY; j++)
1994  {
1995  ComplexH[z][j][i][_RE] = (1. / (4.*RAMat[j][i]*sin2m1h))*((cos0p1h / sin2m0h)*I02Mat[j][i] - (cos0p2h / sin1m0h)*I01Mat[j][i]);
1996  ComplexH[z][j][i][_IM] = (1. / (4.*RAMat[j][i]*sin2m1h))*((sin0p1h / sin2m0h)*I02Mat[j][i] - (sin0p2h / sin1m0h)*I01Mat[j][i]);
1997  }
1998  }
1999 
2000  // update C
2001  C[0] = 0.0;
2002  C[1] = 0.0;
2003  sin1p2h = sin((ps[1] + ps[2]) / 2.);
2004  cos1p2h = cos((ps[1] + ps[2]) / 2.);
2005  for (int i = 0; i < nX; i++)
2006  {
2007  for (int j = 0; j < nY; j++)
2008  {
2009  sinP = ComplexH[z][j][i][_IM] / OAMat[j][i];
2010  cosP = ComplexH[z][j][i][_RE] / OAMat[j][i];
2011  C[0] += abs(sinP*cos0p1h - cosP*sin0p1h);
2012  C[1] += abs(sinP*cos1p2h - cosP*sin1p2h);
2013  }
2014  }
2015  LOG("C1 %d %f %f\n", iter, C[0], C[1]);
2016  C[0] = C[0] / ((double)nXY);
2017  C[1] = C[1] / ((double)nXY);
2018  LOG("C2 %d %f %f\n", iter, C[0], C[1]);
2019 
2021  for (int i = 0; i < nX; i++)
2022  {
2023  for (int j = 0; j < nY; j++)
2024  {
2025  ComplexH[z][j][i][_RE] = ComplexH[z][j][i][_RE] + 0.5;
2026  ComplexH[z][j][i][_IM] = ComplexH[z][j][i][_IM] + 0.5;
2027  }
2028  }
2029  }
2030  }
2031 
2032 
2033  LOG("complex field obtained from 3 interference patterns\n");
2034 
2035  auto end_time = CUR_TIME;
2036 
2037  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
2038 
2039  LOG("Implement time : %.5lf sec\n", during_time);
2040 
2041  return true;
2042 }
2043 
2044 
2045 void ophSig::ophFree(void) {
2046 
2047 }
OPH_BACKWARD
#define OPH_BACKWARD
Definition: define.h:67
ophSig::save
bool save(const char *real, const char *imag)
Save data as bmp or bin file.
Definition: ophSig.cpp:391
oph::abs
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:113
ophSig::cField2Buffer
void cField2Buffer(matrix< Complex< Real >> &src, Complex< Real > **dst, int nx, int ny)
Function for move data from matrix<Complex<Real>> to Complex<Real>
Definition: ophSig.cpp:56
bitmapinfoheader::dibheadersize
uint32_t dibheadersize
Definition: struct.h:57
Openholo::saveAsImg
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:135
Openholo::resCfg
ResolutionConfig resCfg
Definition: Openholo.h:487
ophSig::sigConvertOffaxis
bool sigConvertOffaxis(Real angleX, Real angleY)
Function for Convert complex hologram to off-axis hologram.
Definition: ophSig.cpp:996
hInfo
bitmapinfoheader hInfo
Definition: Openholo.cpp:423
ophSig::loadAsOhc
bool loadAsOhc(const char *fname)
Load data as ohc file.
Definition: ophSig.cpp:198
ophSig::_radius
Real_t _radius
Definition: ophSig.h:491
OphConfig::wave_length
Real * wave_length
Definition: Openholo.h:106
ophSig::cvtOffaxis_CPU
bool cvtOffaxis_CPU(Real angleX, Real angleY)
Definition: ophSig.cpp:1101
ophSig::ColorField2Buffer
void ColorField2Buffer(matrix< Complex< Real >> &src, Complex< Real > **dst, int nx, int ny)
Function for move Color data from matrix<Complex<Real>> to Complex<Real>
Definition: ophSig.cpp:71
ophSig::getComplexHFrom3ArbStepPSDH
bool getComplexHFrom3ArbStepPSDH(const char *f0, const char *f1, const char *f2, const char *fOI, const char *fRI, int nIter)
Extraction of complex field from 3 phase shifted interference patterns with arbitrary unknown shifts...
Definition: ophSig.cpp:1732
oph::uchar
unsigned char uchar
Definition: typedef.h:64
oph::matrix::mulElem
matrix< T > & mulElem(matrix< T > &p)
Definition: mat.h:286
ophSig::sigGetParamAT_GPU
double sigGetParamAT_GPU()
Extraction of distance parameter using axis transfomation by using GPU.
Definition: ophSig_GPU.cpp:456
Real
float Real
Definition: typedef.h:55
CUR_TIME
#define CUR_TIME
Definition: function.h:58
ophSig::fft1
void fft1(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Function for Fast Fourier transform 1D.
Definition: ophSig.cpp:113
ophSig::sigGetParamSF
double sigGetParamSF(float zMax, float zMin, int sampN, float th)
Extraction of distance parameter using sharpness functions.
Definition: ophSig.cpp:1080
ophSig::saveAsOhc
bool saveAsOhc(const char *fname)
Save data as ohc file.
Definition: ophSig.cpp:227
ophSig::sigConvertHPO_GPU
bool sigConvertHPO_GPU(Real depth, Real_t redRate)
Function for convert complex hologram to horizontal parallax only hologram by using GPU...
Definition: ophSig_GPU.cpp:119
ophSig::sigGetParamAT_CPU
double sigGetParamAT_CPU()
Extraction of distance parameter using axis transfomation by using CPU.
Definition: ophSig.cpp:1412
Openholo::checkExtension
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
ophSig::propagationHolo
bool propagationHolo(float depth)
Function for propagation hologram (class data)
Definition: ophSig.cpp:1054
oph::ivec2
structure for 2-dimensional integer vector and its arithmetic.
Definition: ivec.h:66
oph::matrix::resize
matrix< T > & resize(int x, int y)
Definition: mat.h:117
tinyxml2::XMLNode::FirstChild
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
ophSig::setMode
void setMode(bool is_CPU)
Function for select device.
Definition: ophSig.cpp:84
ophSig::fft2
void fft2(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Function for Fast Fourier transform 2D.
Definition: ophSig.cpp:149
ophSig::_foc
Real_t _foc[3]
Definition: ophSig.h:492
ophSig::getComplexHFromPSDH
bool getComplexHFromPSDH(const char *fname0, const char *fname90, const char *fname180, const char *fname270)
Extraction of complex field from 4 phase shifted interference patterns.
Definition: ophSig.cpp:1556
MIN_DOUBLE
#define MIN_DOUBLE
Definition: define.h:148
_Y
#define _Y
Definition: define.h:96
_IM
#define _IM
Definition: complex.h:58
ophSig::readConfig
bool readConfig(const char *fname)
Function for Read parameter.
Definition: ophSig.cpp:1265
Openholo::OHC_encoder
ImgEncoderOhc * OHC_encoder
OHC file format Variables for read and write.
Definition: Openholo.h:495
bitmapinfoheader::imagesize
uint32_t imagesize
Definition: struct.h:63
oph::ImgEncoderOhc::setFieldEncoding
void setFieldEncoding(const FldStore _fldStore, const FldCodeType _fldCodeType)
Definition: ImgCodecOhc.cpp:965
oph::ImgCodecOhc::getComplexFieldData
void getComplexFieldData(OphComplexField &cmplx_field, uint wavelen_idx)
Definition: ImgCodecOhc.h:74
ophSig::_cfgSig
ophSigConfig _cfgSig
Definition: ophSig.h:486
ophSig::sigConvertHPO
bool sigConvertHPO(Real depth, Real_t redRate)
Function for convert complex hologram to horizontal parallax only hologram.
Definition: ophSig.cpp:1014
Openholo::OHC_decoder
ImgDecoderOhc * OHC_decoder
Definition: Openholo.h:496
ophSig::sigConvertCAC_CPU
bool sigConvertCAC_CPU(double red, double green, double blue)
Function for Chromatic aberration compensation filter by using CPU .
Definition: ophSig.cpp:1191
ophSig::sigConvertCAC
bool sigConvertCAC(double red, double green, double blue)
Function for Chromatic aberration compensation filter.
Definition: ophSig.cpp:1035
_X
#define _X
Definition: define.h:92
ophSig::sigGetParamAT
double sigGetParamAT()
Extraction of distance parameter using axis transfomation.
Definition: ophSig.cpp:1387
ophSig::_wavelength_num
int _wavelength_num
Definition: ophSig.h:490
rgbquad
Definition: struct.h:69
true
return true
Definition: Openholo.cpp:434
ophSig::Color_propagationHolo_GPU
bool Color_propagationHolo_GPU(float depth)
Definition: ophSig_GPU.cpp:311
fclose
fclose(infile)
ophSigConfig::NA
Real_t NA
Definition: ophSig.h:73
oph::ImgEncoderOhc::addComplexFieldData
void addComplexFieldData(const OphComplexField &data)
Definition: ImgCodecOhc.cpp:1022
ophSig::ComplexH
OphComplexField * ComplexH
Definition: ophSig.h:487
ophSig::sigConvertCAC_GPU
bool sigConvertCAC_GPU(double red, double green, double blue)
Function for Chromatic aberration compensation filter by using GPU.
Definition: ophSig_GPU.cpp:181
oph::matrix::size
ivec2 size
Definition: mat.h:71
bitmapinfoheader::bitsperpixel
uint16_t bitsperpixel
Definition: struct.h:61
oph::ImgEncoderOhc::setNumOfWavlen
void setNumOfWavlen(const uint n_wavlens)
Definition: ImgCodecOhc.cpp:924
ophSig::is_CPU
bool is_CPU
Definition: ophSig.h:485
_RE
#define _RE
Definition: complex.h:55
ophSigConfig::z
Real_t z
Definition: ophSig.h:74
ophSigConfig::height
Real_t height
Definition: ophSig.h:72
oph::OphComplexField
oph::matrix< Complex< Real > > OphComplexField
Definition: mat.h:420
oph::matrix::mat
std::vector< T > * mat
Definition: mat.h:70
nRead
size_t nRead
Definition: Openholo.cpp:424
ophSig::wavelength_Set
void wavelength_Set(double wavelength)
Definition: ophSig.cpp:1247
ophSig::sigConvertHPO_CPU
bool sigConvertHPO_CPU(Real depth, Real_t redRate)
Function for convert complex hologram to horizontal parallax only hologram by using CPU...
Definition: ophSig.cpp:1145
Real_t
double Real_t
Definition: typedef.h:56
bitmapinfoheader::ypixelpermeter
uint32_t ypixelpermeter
Definition: struct.h:64
ophSig::Wavenumber_output
void Wavenumber_output(int &wavenumber)
Definition: ophSig.cpp:1260
ophSig::absMat
void absMat(matrix< Complex< T >> &src, matrix< T > &dst)
Function for extracts Complex absolute value.
Definition: ophSig.h:657
ophSig::load
bool load(const char *real, const char *imag)
Load bmp or bin file.
Definition: ophSig.cpp:257
ophSig::propagationHolo_CPU
bool propagationHolo_CPU(float depth)
Function for propagation hologram by using CPU.
Definition: ophSig.cpp:1311
oph::ImgCodecOhc::setFileName
bool setFileName(const std::string &_fname)
Definition: ImgCodecOhc.cpp:92
OphConfig::pixel_number
ivec2 pixel_number
Definition: Openholo.h:99
ophSig::propagationHolo_GPU
bool propagationHolo_GPU(float depth)
Function for propagation hologram by using GPU.
Definition: ophSig_GPU.cpp:255
tinyxml2::XMLDocument::LoadFile
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
fileheader::filesize
uint32_t filesize
Definition: struct.h:52
fileheader::fileoffset_to_pixelarray
uint32_t fileoffset_to_pixelarray
Definition: struct.h:54
bitmapinfoheader::width
uint32_t width
Definition: struct.h:58
ophSig::linInterp
void linInterp(vector< T > &X, matrix< Complex< T >> &src, vector< T > &Xq, matrix< Complex< T >> &dst)
Linear interpolation.
Definition: ophSig.cpp:90
ophSig::minOfMat
Real minOfMat(matrix< Real > &src)
Function for extracts minimum of matrix , where matrix is real number.
Definition: ophSig.h:817
ophSig::focal_length_Set
void focal_length_Set(double red, double green, double blue, double rad)
Definition: ophSig.cpp:1252
Openholo::loadAsImg
virtual uchar * loadAsImg(const char *fname)
Function for loading image files.
Definition: Openholo.cpp:321
bitmapinfoheader::height
uint32_t height
Definition: struct.h:59
bitmapinfoheader::planes
uint16_t planes
Definition: struct.h:60
ophSig::maxOfMat
Real maxOfMat(matrix< Real > &src)
Function for extracts maximum of matrix , where matrix is real number.
Definition: ophSig.h:783
oph::ImgEncoderOhc::setWavelength
void setWavelength(const Real _wavlen, const LenUnit _unit=LenUnit::m)
Definition: ImgCodecOhc.cpp:1012
oph::ImgEncoderOhc::setNumOfPixel
void setNumOfPixel(const uint _pxNumX, const uint _pxNumY)
Definition: ImgCodecOhc.cpp:878
ophSig::Data_output
void Data_output(uchar *data, int pos, int bitpixel)
Definition: ophSig.cpp:744
OPH_COMPRESSION
#define OPH_COMPRESSION
Definition: define.h:164
ophSig::Parameter_Set
void Parameter_Set(int nx, int ny, double width, double height, double NA)
Definition: ophSig.cpp:1238
ophSig::ophFree
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: ophSig.cpp:2045
oph::ImgCodecOhc::getOHCheader
void getOHCheader(oph::ohcHeader &_Header)
Definition: ImgCodecOhc.cpp:109
hf
fileheader hf
Definition: Openholo.cpp:422
ophSig::ophSig
ophSig(void)
Constructor.
Definition: ophSig.cpp:50
ophSig::sigGetParamSF_GPU
double sigGetParamSF_GPU(float zMax, float zMin, int sampN, float th)
Extraction of distance parameter using sharpness functions by using GPU.
Definition: ophSig_GPU.cpp:393
Openholo::context_
OphConfig context_
Definition: Openholo.h:486
ophSigConfig::width
Real_t width
Definition: ophSig.h:71
OPH_FORWARD
#define OPH_FORWARD
Definition: define.h:66
fileheader::signature
uint8_t signature[2]
Definition: struct.h:51
OPH_PLANES
#define OPH_PLANES
Definition: define.h:163
ophSig::cvtOffaxis_GPU
void cvtOffaxis_GPU(Real angleX, Real angleY)
Definition: ophSig_GPU.cpp:64
oph::ImgDecoderOhc::getWavelength
void getWavelength(std::vector< double_t > &wavlen_array)
Definition: ImgCodecOhc.cpp:289
tinyxml2::XMLNode::FirstChildElement
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
bitmapinfoheader::compression
uint32_t compression
Definition: struct.h:62
ophSig::linspace
vector< T > linspace(T first, T last, int len)
Generate linearly spaced vector.
Definition: ophSig.h:645
oph::uint
unsigned int uint
Definition: typedef.h:62
bitmapinfoheader::xpixelpermeter
uint32_t xpixelpermeter
Definition: struct.h:65
M_PI
#define M_PI
Definition: define.h:52
ophSig::sigGetParamSF_CPU
double sigGetParamSF_CPU(float zMax, float zMin, int sampN, float th)
Extraction of distance parameter using sharpness functions by using CPU.
Definition: ophSig.cpp:1510