Openholo  v4.2
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 
1303  return true;
1304 }
1305 
1306 
1307 
1308 bool ophSig::propagationHolo_CPU(float depth) {
1309  int i, j;
1310  Real x, y, sigmaf;
1311  int nx = context_.pixel_number[_X];
1312  int ny = context_.pixel_number[_Y];
1313 
1314  OphComplexField FH(nx, ny);
1315  //OphComplexField FFZP(nx, ny);
1316  int xshift = nx / 2;
1317  int yshift = ny / 2;
1318 
1319  for (int z = 0; z < _wavelength_num; z++) {
1320 
1321  sigmaf = (depth * context_.wave_length[z]) / (4 * M_PI);
1322 
1323  /*FH.resize(nx, ny);
1324  FFZP.resize(nx, ny);*/
1325 
1326  fft2(ComplexH[z], FH, OPH_FORWARD);
1327 
1328  for (i = 0; i < ny; i++)
1329  {
1330  int ii = (i + yshift) % ny;
1331  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1332 
1333  for (j = 0; j < nx; j++)
1334  {
1335  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1336  int jj = (j + xshift) % nx;
1337  double temp = FH(jj, ii)[_RE];
1338  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];
1339  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];
1340 
1341  }
1342  }
1343 
1344  fft2(FH, ComplexH[z], OPH_BACKWARD);
1345  }
1346  return true;
1347 }
1348 
1349 OphComplexField ophSig::propagationHolo(OphComplexField complexH, float depth) {
1350  int i, j;
1351  Real x, y, sigmaf;
1352  int nx = context_.pixel_number[_X];
1353  int ny = context_.pixel_number[_Y];
1354 
1355  OphComplexField FH(nx, ny);
1356  int xshift = nx / 2;
1357  int yshift = ny / 2;
1358 
1359 
1360  sigmaf = (depth * (*context_.wave_length)) / (4 * M_PI);
1361 
1362  fft2(complexH, FH);
1363 
1364  for (i = 0; i < ny; i++)
1365  {
1366  int ii = (i + yshift) % ny;
1367  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1368 
1369  for (j = 0; j < nx; j++)
1370  {
1371  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1372  int jj = (j + xshift) % nx;
1373  double temp = FH(jj, ii)[_RE];
1374  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];
1375  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];
1376 
1377  }
1378  }
1379  fft2(FH, complexH, OPH_BACKWARD);
1380 
1381  return complexH;
1382 }
1383 
1384 double ophSig::sigGetParamAT()
1385 {
1386  //auto start_time = CUR_TIME;
1387 
1388  Real index = 0;
1389  if (is_CPU == true)
1390  {
1391  //std::cout << "Start Single Core CPU" << endl;
1392  index = sigGetParamAT_CPU();
1393 
1394  }
1395  else {
1396  //std::cout << "Start Multi Core GPU" << std::endl;
1397  index = sigGetParamAT_GPU();
1398 
1399  }
1400 
1401  //auto end_time = CUR_TIME;
1402 
1403  //auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1404 
1405  //LOG("Implement time : %.5lf sec\n", during_time);
1406  return index;
1407 }
1408 
1409 double ophSig::sigGetParamAT_CPU() {
1410 
1411  int i = 0, j = 0;
1412  Real max = 0; Real index = 0;
1413  Real_t NA_g = (Real_t)0.025;
1414  int nx = context_.pixel_number[_X];
1415  int ny = context_.pixel_number[_Y];
1416 
1417  OphComplexField Flr(nx, ny);
1418  OphComplexField Fli(nx, ny);
1419  OphComplexField Hsyn(nx, ny);
1420 
1421  OphComplexField Fo(nx, ny);
1422  OphComplexField Fon, yn, Ab_yn;
1423 
1424  OphRealField Ab_yn_half;
1425  OphRealField G(nx, ny);
1426  Real x = 1, y = 1;
1427  vector<Real> t, tn;
1428 
1429  for (i = 0; i < nx; i++)
1430  {
1431  for (j = 0; j < ny; j++)
1432  {
1433 
1434  x = (2 * M_PI*(i) / _cfgSig.height - M_PI*(nx - 1) / _cfgSig.height);
1435  y = (2 * M_PI*(j) / _cfgSig.width - M_PI*(ny - 1) / _cfgSig.width);
1436  G(i, j) = std::exp(-M_PI * pow((*context_.wave_length) / (2 * M_PI * NA_g), 2) * (pow(y, 2) + pow(x, 2)));
1437  Flr(i, j)[_RE] = (*ComplexH)(i, j)[_RE];
1438  Fli(i, j)[_RE] = (*ComplexH)(i, j)[_IM];
1439  Flr(i, j)[_IM] = 0;
1440  Fli(i, j)[_IM] = 0;
1441  }
1442  }
1443 
1444  fft2(Flr, Flr);
1445  fft2(Fli, Fli);
1446 
1447  int xshift = nx / 2;
1448  int yshift = ny / 2;
1449 
1450  for (i = 0; i < nx; i++)
1451  {
1452  int ii = (i + xshift) % nx;
1453  for (j = 0; j < ny; j++)
1454  {
1455  int jj = (j + yshift) % ny;
1456  Hsyn(i, j)[_RE] = Flr(i, j)[_RE] * G(i, j);
1457  Hsyn(i, j)[_IM] = Fli(i, j)[_RE] * G(i, j);
1458  /*Hsyn_copy1(i, j) = Hsyn(i, j);
1459  Hsyn_copy2(i, j) = Hsyn_copy1(i, j) * Hsyn(i, j);
1460  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);
1461  Fo(ii, jj)[_RE] = Hsyn_copy2(i, j)[0] / Hsyn_copy3(i, j);
1462  Fo(ii, jj)[_IM] = Hsyn_copy2(i, j)[1] / Hsyn_copy3(i, j);*/
1463  Fo(ii, jj) = pow(Hsyn(i, j), 2) / (pow(abs(Hsyn(i, j)), 2) + pow(10, -300));
1464 
1465  }
1466  }
1467 
1468  t = linspace(0., 1., nx / 2 + 1);
1469  tn.resize(t.size());
1470  Fon.resize(1, t.size());
1471 
1472  int size = (int)tn.size();
1473  for (int i = 0; i < size; i++)
1474  {
1475  tn.at(i) = pow(t.at(i), 0.5);
1476  Fon(0, i)[_RE] = Fo(nx / 2 - 1, nx / 2 - 1 + i)[_RE];
1477  Fon(0, i)[_IM] = 0;
1478  }
1479  yn.resize(1, size);
1480  linInterp(t, Fon, tn, yn);
1481  fft1(yn, yn);
1482  Ab_yn.resize(yn.size[_X], yn.size[_Y]);
1483  absMat(yn, Ab_yn);
1484  Ab_yn_half.resize(1, nx / 4 + 1);
1485 
1486  for (int i = 0; i < nx / 4 + 1; i++)
1487  {
1488  Ab_yn_half(0, i) = Ab_yn(0, nx / 4 + i - 1)[_RE];
1489  if (i == 0) max = Ab_yn_half(0, 0);
1490  else
1491  {
1492  if (Ab_yn_half(0, i) > max)
1493  {
1494  max = Ab_yn_half(0, i);
1495  index = i;
1496  }
1497  }
1498  }
1499 
1500  index = -(((index + 1) - 120) / 10) / 140 + 0.1;
1501 
1502 
1503 
1504  return index;
1505 }
1506 
1507 double ophSig::sigGetParamSF_CPU(float zMax, float zMin, int sampN, float th)
1508 {
1509  int nx = context_.pixel_number[_X];
1510  int ny = context_.pixel_number[_Y];
1511 
1512  OphComplexField I(nx, ny);
1513  vector<Real> F;
1514  Real dz = (zMax - zMin) / sampN;
1515  Real f;
1516  Real_t z = 0;
1517  Real depth = 0;
1518  Real max = MIN_DOUBLE;
1519  int i, j, n = 0;
1520  Real ret1;
1521  Real ret2;
1522 
1523 
1524 
1525  for (n = 0; n < sampN + 1; n++)
1526  {
1527  OphComplexField field = *(ComplexH);
1528  z = ((n)* dz + zMin);
1529  f = 0;
1530  I = propagationHolo(field, z);
1531 
1532  for (i = 0; i < nx - 2; i++)
1533  {
1534  for (j = 0; j < ny - 2; j++)
1535  {
1536  ret1 = abs(I(i + 2, j)[_RE] - I(i, j)[_RE]);
1537  ret2 = abs(I(i, j + 2)[_RE] - I(i, j)[_RE]);
1538  if (ret1 >= th) { f += ret1 * ret1; }
1539  else if (ret2 >= th) { f += ret2 * ret2; }
1540  }
1541  }
1542 
1543  if (f > max) {
1544  max = f;
1545  depth = z;
1546  }
1547  }
1548 
1549 
1550  return depth;
1551 }
1552 
1553 bool ophSig::getComplexHFromPSDH(const char * fname0, const char * fname90, const char * fname180, const char * fname270)
1554 {
1555  auto start_time = CUR_TIME;
1556  string fname0str = fname0;
1557  string fname90str = fname90;
1558  string fname180str = fname180;
1559  string fname270str = fname270;
1560  int checktype = static_cast<int>(fname0str.rfind("."));
1561  OphRealField f0Mat[3], f90Mat[3], f180Mat[3], f270Mat[3];
1562 
1563  std::string f0type = fname0str.substr(checktype + 1, fname0str.size());
1564 
1565  uint16_t bitsperpixel;
1566  fileheader hf;
1567  bitmapinfoheader hInfo;
1568 
1569  if (f0type == "bmp")
1570  {
1571  FILE *f0, *f90, *f180, *f270;
1572  f0 = fopen(fname0str.c_str(), "rb");
1573  f90 = fopen(fname90str.c_str(), "rb");
1574  f180 = fopen(fname180str.c_str(), "rb");
1575  f270 = fopen(fname270str.c_str(), "rb");
1576  if (!f0)
1577  {
1578  LOG("bmp file open fail! (phase shift = 0)\n");
1579  return false;
1580  }
1581  if (!f90)
1582  {
1583  LOG("bmp file open fail! (phase shift = 90)\n");
1584  return false;
1585  }
1586  if (!f180)
1587  {
1588  LOG("bmp file open fail! (phase shift = 180)\n");
1589  return false;
1590  }
1591  if (!f270)
1592  {
1593  LOG("bmp file open fail! (phase shift = 270)\n");
1594  return false;
1595  }
1596  size_t nRead = fread(&hf, sizeof(fileheader), 1, f0);
1597  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f0);
1598 
1599  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
1600  if ((hInfo.height == 0) || (hInfo.width == 0))
1601  {
1602  LOG("bmp header is empty!\n");
1603  hInfo.height = context_.pixel_number[_X];
1604  hInfo.width = context_.pixel_number[_Y];
1605  if (hInfo.height == 0 || hInfo.width == 0)
1606  {
1607  LOG("check your parameter file!\n");
1608  return false;
1609  }
1610  }
1611  if ((context_.pixel_number[_Y] != (int)hInfo.height) || (context_.pixel_number[_X] != (int)hInfo.width)) {
1612  LOG("image size is different!\n");
1613  context_.pixel_number[_Y] = (int)hInfo.height;
1614  context_.pixel_number[_X] = (int)hInfo.width;
1615  LOG("changed parameter of size %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
1616  }
1617  bitsperpixel = hInfo.bitsperpixel;
1618  if (hInfo.bitsperpixel == 8)
1619  {
1620  _wavelength_num = 1;
1621  rgbquad palette[256];
1622  size_t nRead;
1623  nRead = fread(palette, sizeof(rgbquad), 256, f0);
1624  nRead = fread(palette, sizeof(rgbquad), 256, f90);
1625  nRead = fread(palette, sizeof(rgbquad), 256, f180);
1626  nRead = fread(palette, sizeof(rgbquad), 256, f270);
1627 
1628  f0Mat[0].resize(hInfo.height, hInfo.width);
1629  f90Mat[0].resize(hInfo.height, hInfo.width);
1630  f180Mat[0].resize(hInfo.height, hInfo.width);
1631  f270Mat[0].resize(hInfo.height, hInfo.width);
1632  ComplexH = new OphComplexField;
1633  ComplexH[0].resize(hInfo.height, hInfo.width);
1634  }
1635  else
1636  {
1637  _wavelength_num = 3;
1638  ComplexH = new OphComplexField[3];
1639  f0Mat[0].resize(hInfo.height, hInfo.width);
1640  f90Mat[0].resize(hInfo.height, hInfo.width);
1641  f180Mat[0].resize(hInfo.height, hInfo.width);
1642  f270Mat[0].resize(hInfo.height, hInfo.width);
1643  ComplexH[0].resize(hInfo.height, hInfo.width);
1644 
1645  f0Mat[1].resize(hInfo.height, hInfo.width);
1646  f90Mat[1].resize(hInfo.height, hInfo.width);
1647  f180Mat[1].resize(hInfo.height, hInfo.width);
1648  f270Mat[1].resize(hInfo.height, hInfo.width);
1649  ComplexH[1].resize(hInfo.height, hInfo.width);
1650 
1651  f0Mat[2].resize(hInfo.height, hInfo.width);
1652  f90Mat[2].resize(hInfo.height, hInfo.width);
1653  f180Mat[2].resize(hInfo.height, hInfo.width);
1654  f270Mat[2].resize(hInfo.height, hInfo.width);
1655  ComplexH[2].resize(hInfo.height, hInfo.width);
1656  }
1657 
1658  size_t size = hInfo.width * hInfo.height * (hInfo.bitsperpixel >> 3);
1659  uchar* f0data = (uchar*)malloc(sizeof(uchar) * size);
1660  uchar* f90data = (uchar*)malloc(sizeof(uchar) * size);
1661  uchar* f180data = (uchar*)malloc(sizeof(uchar) * size);
1662  uchar* f270data = (uchar*)malloc(sizeof(uchar) * size);
1663 
1664  nRead = fread(f0data, sizeof(uchar), size, f0);
1665  nRead = fread(f90data, sizeof(uchar), size, f90);
1666  nRead = fread(f180data, sizeof(uchar), size, f180);
1667  nRead = fread(f270data, sizeof(uchar), size, f270);
1668 
1669  fclose(f0);
1670  fclose(f90);
1671  fclose(f180);
1672  fclose(f270);
1673 
1674  uint16_t bpp = hInfo.bitsperpixel >> 3;
1675  uint32_t nLine = hInfo.width * bpp;
1676  for (int i = hInfo.height - 1; i >= 0; i--)
1677  {
1678  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
1679  {
1680  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1681  {
1682  int idx = hInfo.height - i - 1;
1683  int idx2 = i * nLine + bpp * j + z;
1684  f0Mat[z](idx, j) = (double)f0data[idx2];
1685  f90Mat[z](idx, j) = (double)f90data[idx2];
1686  f180Mat[z](idx, j) = (double)f180data[idx2];
1687  f270Mat[z](idx, j) = (double)f270data[idx2];
1688  }
1689  }
1690  }
1691  LOG("PSDH file load complete!\n");
1692 
1693  free(f0data);
1694  free(f90data);
1695  free(f180data);
1696  free(f270data);
1697 
1698  }
1699  else
1700  {
1701  LOG("wrong type (only BMP supported)\n");
1702  }
1703 
1704  // calculation complexH from 4 psdh and then normalize
1705  double normalizefactor = 1. / 256.;
1706  for (int z = 0; z < (hInfo.bitsperpixel >> 3); z++)
1707  {
1708  for (int i = 0; i < context_.pixel_number[_X]; i++)
1709  {
1710  for (int j = 0; j < context_.pixel_number[_Y]; j++)
1711  {
1712  ComplexH[z][j][i][_RE] = (f0Mat[z][j][i] - f180Mat[z][j][i])*normalizefactor;
1713  ComplexH[z][j][i][_IM] = (f90Mat[z][j][i] - f270Mat[z][j][i])*normalizefactor;
1714 
1715  }
1716  }
1717  }
1718  LOG("complex field obtained from 4 psdh\n");
1719 
1720  auto end_time = CUR_TIME;
1721 
1722  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1723 
1724  LOG("Implement time : %.5lf sec\n", during_time);
1725 
1726  return true;
1727 }
1728 
1729 bool ophSig::getComplexHFrom3ArbStepPSDH(
1730  const char* fname0, const char* fname1,
1731  const char* fname2, const char* fnameOI,
1732  const char* fnameRI, int nIter)
1733 {
1734  auto start_time = CUR_TIME;
1735  string fname0str = fname0;
1736  string fname1str = fname1;
1737  string fname2str = fname2;
1738  string fnameOIstr = fnameOI;
1739  string fnameRIstr = fnameRI;
1740  int checktype = static_cast<int>(fname0str.rfind("."));
1741  OphRealField f0Mat[3], f1Mat[3], f2Mat[3], fOIMat[3], fRIMat[3];
1742 
1743  std::string f0type = fname0str.substr(checktype + 1, fname0str.size());
1744 
1745  uint16_t bitsperpixel;
1746  fileheader hf;
1747  bitmapinfoheader hInfo;
1748 
1749  if (f0type == "bmp")
1750  {
1751  FILE *f0, *f1, *f2, *fOI, *fRI;
1752  f0 = fopen(fname0str.c_str(), "rb");
1753  f1 = fopen(fname1str.c_str(), "rb");
1754  f2 = fopen(fname2str.c_str(), "rb");
1755  fOI = fopen(fnameOIstr.c_str(), "rb");
1756  fRI = fopen(fnameRIstr.c_str(), "rb");
1757  if (!f0)
1758  {
1759  LOG("bmp file open fail! (first interference pattern)\n");
1760  return false;
1761  }
1762  if (!f1)
1763  {
1764  LOG("bmp file open fail! (second interference pattern)\n");
1765  return false;
1766  }
1767  if (!f2)
1768  {
1769  LOG("bmp file open fail! (third interference pattern)\n");
1770  return false;
1771  }
1772  if (!fOI)
1773  {
1774  LOG("bmp file open fail! (object wave intensity pattern)\n");
1775  return false;
1776  }
1777  if (!fRI)
1778  {
1779  LOG("bmp file open fail! (reference wave intensity pattern)\n");
1780  return false;
1781  }
1782 
1783  size_t nRead;
1784  nRead = fread(&hf, sizeof(fileheader), 1, f0);
1785  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f0);
1786  nRead = fread(&hf, sizeof(fileheader), 1, f1);
1787  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f1);
1788  nRead = fread(&hf, sizeof(fileheader), 1, f2);
1789  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, f2);
1790  nRead = fread(&hf, sizeof(fileheader), 1, fOI);
1791  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, fOI);
1792  nRead = fread(&hf, sizeof(fileheader), 1, fRI);
1793  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, fRI);
1794 
1795  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
1796  if ((hInfo.height == 0) || (hInfo.width == 0))
1797  {
1798  LOG("bmp header is empty!\n");
1799  hInfo.height = context_.pixel_number[_X];
1800  hInfo.width = context_.pixel_number[_Y];
1801  if (hInfo.height == 0 || hInfo.width == 0)
1802  {
1803  LOG("check your parameter file!\n");
1804  return false;
1805  }
1806  }
1807  if ((context_.pixel_number[_Y] != (int)hInfo.height) || (context_.pixel_number[_X] != (int)hInfo.width)) {
1808  LOG("image size is different!\n");
1809  context_.pixel_number[_Y] = (int)hInfo.height;
1810  context_.pixel_number[_X] = (int)hInfo.width;
1811  LOG("changed parameter of size %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
1812  }
1813  bitsperpixel = hInfo.bitsperpixel;
1814  if (hInfo.bitsperpixel == 8)
1815  {
1816  _wavelength_num = 1;
1817  rgbquad palette[256];
1818 
1819  size_t nRead;
1820  nRead = fread(palette, sizeof(rgbquad), 256, f0);
1821  nRead = fread(palette, sizeof(rgbquad), 256, f1);
1822  nRead = fread(palette, sizeof(rgbquad), 256, f2);
1823  nRead = fread(palette, sizeof(rgbquad), 256, fOI);
1824  nRead = fread(palette, sizeof(rgbquad), 256, fRI);
1825 
1826  f0Mat[0].resize(hInfo.height, hInfo.width);
1827  f1Mat[0].resize(hInfo.height, hInfo.width);
1828  f2Mat[0].resize(hInfo.height, hInfo.width);
1829  fOIMat[0].resize(hInfo.height, hInfo.width);
1830  fRIMat[0].resize(hInfo.height, hInfo.width);
1831  ComplexH = new OphComplexField;
1832  ComplexH[0].resize(hInfo.height, hInfo.width);
1833  }
1834  else
1835  {
1836  _wavelength_num = 3;
1837  ComplexH = new OphComplexField[3];
1838  f0Mat[0].resize(hInfo.height, hInfo.width);
1839  f1Mat[0].resize(hInfo.height, hInfo.width);
1840  f2Mat[0].resize(hInfo.height, hInfo.width);
1841  fOIMat[0].resize(hInfo.height, hInfo.width);
1842  fRIMat[0].resize(hInfo.height, hInfo.width);
1843  ComplexH[0].resize(hInfo.height, hInfo.width);
1844 
1845  f0Mat[1].resize(hInfo.height, hInfo.width);
1846  f1Mat[1].resize(hInfo.height, hInfo.width);
1847  f2Mat[1].resize(hInfo.height, hInfo.width);
1848  fOIMat[1].resize(hInfo.height, hInfo.width);
1849  fRIMat[1].resize(hInfo.height, hInfo.width);
1850  ComplexH[1].resize(hInfo.height, hInfo.width);
1851 
1852  f0Mat[2].resize(hInfo.height, hInfo.width);
1853  f1Mat[2].resize(hInfo.height, hInfo.width);
1854  f2Mat[2].resize(hInfo.height, hInfo.width);
1855  fOIMat[2].resize(hInfo.height, hInfo.width);
1856  fRIMat[2].resize(hInfo.height, hInfo.width);
1857  ComplexH[2].resize(hInfo.height, hInfo.width);
1858  }
1859 
1860  size_t size = hInfo.width * hInfo.height * (hInfo.bitsperpixel >> 3);
1861 
1862  uchar* f0data = (uchar*)malloc(sizeof(uchar) * size);
1863  uchar* f1data = (uchar*)malloc(sizeof(uchar) * size);
1864  uchar* f2data = (uchar*)malloc(sizeof(uchar) * size);
1865  uchar* fOIdata = (uchar*)malloc(sizeof(uchar) * size);
1866  uchar* fRIdata = (uchar*)malloc(sizeof(uchar) * size);
1867 
1868  nRead = fread(f0data, sizeof(uchar), size, f0);
1869  nRead = fread(f1data, sizeof(uchar), size, f1);
1870  nRead = fread(f2data, sizeof(uchar), size, f2);
1871  nRead = fread(fOIdata, sizeof(uchar), size, fOI);
1872  nRead = fread(fRIdata, sizeof(uchar), size, fRI);
1873 
1874  fclose(f0);
1875  fclose(f1);
1876  fclose(f2);
1877  fclose(fOI);
1878  fclose(fRI);
1879 
1880  uint16_t bpp = hInfo.bitsperpixel >> 3;
1881 
1882  for (int i = hInfo.height - 1; i >= 0; i--)
1883  {
1884  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
1885  {
1886  for (int z = 0; z < (hInfo.bitsperpixel >> 3); z++)
1887  {
1888  int idx = hInfo.height - i - 1;
1889  size_t idx2 = i * size + bpp * j + z;
1890  f0Mat[z](idx, j) = (double)f0data[idx2];
1891  f1Mat[z](idx, j) = (double)f1data[idx2];
1892  f2Mat[z](idx, j) = (double)f2data[idx2];
1893  fOIMat[z](idx, j) = (double)fOIdata[idx2];
1894  fRIMat[z](idx, j) = (double)fRIdata[idx2];
1895  }
1896  }
1897  }
1898 
1899  LOG("PSDH_3ArbStep file load complete!\n");
1900 
1901  free(f0data);
1902  free(f1data);
1903  free(f2data);
1904  free(fOIdata);
1905  free(fRIdata);
1906 
1907  }
1908  else
1909  {
1910  LOG("wrong type (only BMP supported)\n");
1911  }
1912 
1913  // calculation complexH from 3 arbitrary step intereference patterns and the object wave intensity
1914  // prepare some variables
1915  double P[2] = { 0.0, }; // please see ref.
1916  double C[2] = { 2.0/M_PI, 2.0/M_PI };
1917  double alpha[2] = { 0.0, }; //phaseShift[j+1]-phaseShift[j]
1918  double ps[3] = { 0.0, }; // reference wave phase shift for each inteference pattern
1919  const int nX = context_.pixel_number[_X];
1920  const int nY = context_.pixel_number[_Y];
1921  const int nXY = nX * nY;
1922 
1923 
1924  // calculate difference between interference patterns
1925  OphRealField I01Mat, I02Mat, I12Mat, OAMat, RAMat;
1926  I01Mat.resize(nY, nX);
1927  I02Mat.resize(nY, nX);
1928  I12Mat.resize(nY, nX);
1929  OAMat.resize(nY, nX);
1930  RAMat.resize(nY, nX);
1931 
1932  double sin2m1h, sin2m0h, sin1m0h, sin0p1h, sin0p2h, cos0p1h, cos0p2h, sin1p2h, cos1p2h;
1933  double sinP, cosP;
1934  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1935  {
1936  // initialize
1937  P[0] = 0.0;
1938  P[1] = 0.0;
1939  C[0] = 2.0 / M_PI;
1940  C[1] = 2.0 / M_PI;
1941 
1942  // load current channel
1943  for (int i = 0; i < nX; i++)
1944  {
1945  for (int j = 0; j < nY; j++)
1946  {
1947  I01Mat[j][i] = (f0Mat[z][j][i] - f1Mat[z][j][i]) / 255.; // difference & normalize
1948  I02Mat[j][i] = (f0Mat[z][j][i] - f2Mat[z][j][i]) / 255.; // difference & normalize
1949  I12Mat[j][i] = (f1Mat[z][j][i] - f2Mat[z][j][i]) / 255.; // difference & normalize
1950  OAMat[j][i] = sqrt(fOIMat[z][j][i] / 255.); // normalize & then calculate amplitude from intensity
1951  RAMat[j][i] = sqrt(fRIMat[z][j][i] / 255.); // normalize & then calculate amplitude from intensity
1952  }
1953  }
1954 
1955  // calculate P
1956  for (int i = 0; i < nX; i++)
1957  {
1958  for (int j = 0; j < nY; j++)
1959  {
1960  P[0] += abs(I01Mat[j][i] / OAMat[j][i] / RAMat[j][i]);
1961  P[1] += abs(I12Mat[j][i] / OAMat[j][i] / RAMat[j][i]);
1962  }
1963  }
1964  P[0] = P[0] / (4.*((double) nXY));
1965  P[1] = P[1] / (4.*((double) nXY));
1966  LOG("P %f %f\n", P[0], P[1]);
1967 
1968  // iterative search
1969  for (int iter = 0; iter < nIter; iter++)
1970  {
1971  LOG("C %d %f %f\n", iter, C[0], C[1]);
1972  LOG("ps %d %f %f %f\n", iter, ps[0], ps[1], ps[2]);
1973 
1974  alpha[0] = 2.*asin(P[0] / C[0]);
1975  alpha[1] = 2.*asin(P[1] / C[1]);
1976 
1977  ps[0] = 0.0;
1978  ps[1] = ps[0] + alpha[0];
1979  ps[2] = ps[1] + alpha[1];
1980 
1981  sin2m1h = sin((ps[2] - ps[1]) / 2.);
1982  sin2m0h = sin((ps[2] - ps[0]) / 2.);
1983  sin1m0h = sin((ps[1] - ps[0]) / 2.);
1984  sin0p1h = sin((ps[0] + ps[1]) / 2.);
1985  sin0p2h = sin((ps[0] + ps[2]) / 2.);
1986  cos0p1h = cos((ps[0] + ps[1]) / 2.);
1987  cos0p2h = cos((ps[0] + ps[2]) / 2.);
1988  for (int i = 0; i < nX; i++)
1989  {
1990  for (int j = 0; j < nY; j++)
1991  {
1992  ComplexH[z][j][i][_RE] = (1. / (4.*RAMat[j][i]*sin2m1h))*((cos0p1h / sin2m0h)*I02Mat[j][i] - (cos0p2h / sin1m0h)*I01Mat[j][i]);
1993  ComplexH[z][j][i][_IM] = (1. / (4.*RAMat[j][i]*sin2m1h))*((sin0p1h / sin2m0h)*I02Mat[j][i] - (sin0p2h / sin1m0h)*I01Mat[j][i]);
1994  }
1995  }
1996 
1997  // update C
1998  C[0] = 0.0;
1999  C[1] = 0.0;
2000  sin1p2h = sin((ps[1] + ps[2]) / 2.);
2001  cos1p2h = cos((ps[1] + ps[2]) / 2.);
2002  for (int i = 0; i < nX; i++)
2003  {
2004  for (int j = 0; j < nY; j++)
2005  {
2006  sinP = ComplexH[z][j][i][_IM] / OAMat[j][i];
2007  cosP = ComplexH[z][j][i][_RE] / OAMat[j][i];
2008  C[0] += abs(sinP*cos0p1h - cosP*sin0p1h);
2009  C[1] += abs(sinP*cos1p2h - cosP*sin1p2h);
2010  }
2011  }
2012  LOG("C1 %d %f %f\n", iter, C[0], C[1]);
2013  C[0] = C[0] / ((double)nXY);
2014  C[1] = C[1] / ((double)nXY);
2015  LOG("C2 %d %f %f\n", iter, C[0], C[1]);
2016 
2018  for (int i = 0; i < nX; i++)
2019  {
2020  for (int j = 0; j < nY; j++)
2021  {
2022  ComplexH[z][j][i][_RE] = ComplexH[z][j][i][_RE] + 0.5;
2023  ComplexH[z][j][i][_IM] = ComplexH[z][j][i][_IM] + 0.5;
2024  }
2025  }
2026  }
2027  }
2028 
2029 
2030  LOG("complex field obtained from 3 interference patterns\n");
2031 
2032  auto end_time = CUR_TIME;
2033 
2034  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
2035 
2036  LOG("Implement time : %.5lf sec\n", during_time);
2037 
2038  return true;
2039 }
2040 
2041 
2042 void ophSig::ophFree(void) {
2043 
2044 }
ophSig::save
bool save(const char *real, const char *imag)
Save data as bmp or bin file.
Definition: ophSig.cpp:391
oph::OphComplexField
oph::matrix< Complex< Real > > OphComplexField
Definition: mat.h:420
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
M_PI
#define M_PI
Definition: define.h:52
OPH_FORWARD
#define OPH_FORWARD
Definition: define.h:66
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
nRead
size_t nRead
Definition: Openholo.cpp:424
ophSig::sigConvertOffaxis
bool sigConvertOffaxis(Real angleX, Real angleY)
Function for Convert complex hologram to off-axis hologram.
Definition: ophSig.cpp:996
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
ophSig::cvtOffaxis_CPU
bool cvtOffaxis_CPU(Real angleX, Real angleY)
Definition: ophSig.cpp:1101
MIN_DOUBLE
#define MIN_DOUBLE
Definition: define.h:148
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:1729
oph::uchar
unsigned char uchar
Definition: typedef.h:64
ophSig::sigGetParamAT_GPU
double sigGetParamAT_GPU()
Extraction of distance parameter using axis transfomation by using GPU.
Definition: ophSig_GPU.cpp:456
ophSig::_foc
Real_t _foc[3]
Definition: ophSig.h:492
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
true
return true
Definition: Openholo.cpp:434
Real
float Real
Definition: typedef.h:55
ophSig::sigGetParamAT_CPU
double sigGetParamAT_CPU()
Extraction of distance parameter using axis transfomation by using CPU.
Definition: ophSig.cpp:1409
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
_IM
#define _IM
Definition: complex.h:58
Real_t
double Real_t
Definition: typedef.h:56
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::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:1553
fclose
fclose(infile)
ophSig::readConfig
bool readConfig(const char *fname)
Function for Read parameter.
Definition: ophSig.cpp:1265
bitmapinfoheader::imagesize
uint32_t imagesize
Definition: struct.h:63
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
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
ophSig::sigGetParamAT
double sigGetParamAT()
Extraction of distance parameter using axis transfomation.
Definition: ophSig.cpp:1384
ophSig::_wavelength_num
int _wavelength_num
Definition: ophSig.h:490
_X
#define _X
Definition: define.h:92
fileheader::signature
uint8_t signature[2]
Definition: struct.h:51
rgbquad
Definition: struct.h:69
ophSig::Color_propagationHolo_GPU
bool Color_propagationHolo_GPU(float depth)
Definition: ophSig_GPU.cpp:311
ophSigConfig::NA
Real_t NA
Definition: ophSig.h:73
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
bitmapinfoheader::bitsperpixel
uint16_t bitsperpixel
Definition: struct.h:61
OPH_PLANES
#define OPH_PLANES
Definition: define.h:163
ophSig::ComplexH
OphComplexField * ComplexH
Definition: ophSig.h:487
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::OphRealField
oph::matrix< Real > OphRealField
Definition: mat.h:418
Openholo::OHC_decoder
ImgDecoderOhc * OHC_decoder
Definition: Openholo.h:495
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
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:1308
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
_Y
#define _Y
Definition: define.h:96
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
OPH_COMPRESSION
#define OPH_COMPRESSION
Definition: define.h:164
hInfo
bitmapinfoheader hInfo
Definition: Openholo.cpp:423
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
ophSig::Data_output
void Data_output(uchar *data, int pos, int bitpixel)
Definition: ophSig.cpp:744
ophSig::Parameter_Set
void Parameter_Set(int nx, int ny, double width, double height, double NA)
Definition: ophSig.cpp:1238
hf
fileheader hf
Definition: Openholo.cpp:422
ophSig::ophFree
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: ophSig.cpp:2042
Openholo::OHC_encoder
ImgEncoderOhc * OHC_encoder
OHC file format Variables for read and write.
Definition: Openholo.h:494
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:485
ophSigConfig::width
Real_t width
Definition: ophSig.h:71
OPH_BACKWARD
#define OPH_BACKWARD
Definition: define.h:67
ophSig::cvtOffaxis_GPU
void cvtOffaxis_GPU(Real angleX, Real angleY)
Definition: ophSig_GPU.cpp:64
OphConfig::wave_length
Real * wave_length
Definition: Openholo.h:106
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
CUR_TIME
#define CUR_TIME
Definition: function.h:58
bitmapinfoheader::xpixelpermeter
uint32_t xpixelpermeter
Definition: struct.h:65
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:1507