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