Openholo  v1.0
Open Source Digital Holographic Library
ophSig.cpp
Go to the documentation of this file.
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install, copy or use the software.
7 //
8 //
9 // License Agreement
10 // For Open Source Digital Holographic Library
11 //
12 // Openholo library is free software;
13 // you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
14 //
15 // Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
16 // E-mail : contact.openholo@gmail.com
17 // Web : http://www.openholo.org
18 //
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
21 //
22 // 1. Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
24 //
25 // 2. Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other materials provided with the distribution.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the copyright holder or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 // This software contains opensource software released under GNU Generic Public License,
41 // NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
42 // Check whether software you use contains licensed software.
43 //
44 //M*/
45 
46 #include "ophSig.h"
47 #include "include.h"
48 
49 
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);
210  int wavelength_num = OHC_decoder->getNumOfWavlen();
211 
213 
215 
217 
218  for (int i = 0; i < _wavelength_num; i++)
219  {
220  context_.wave_length[i] = wavelengthArray[(_wavelength_num - 1) - i];
221 
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  OHC_encoder->getOHCheader(header);
237 
239 
240  OHC_encoder->setFieldEncoding(FldStore::Directly, FldCodeType::RI);
241 
243 
244  for (int i = _wavelength_num - 1; i >= 0; i--)
245  {
246  //int wl = context_.wave_length[i] * 1000000000;
247  OHC_encoder->setWavelength(context_.wave_length[i], LenUnit::nm);
248 
250  }
251 
252  if (!OHC_encoder->save()) return false;
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] = { "","","" };
263 
264  if (_wavelength_num > 1) {
265  RGB_name[0] = "_B";
266  RGB_name[1] = "_G";
267  RGB_name[2] = "_R";
268  }
269 
270  int checktype = static_cast<int>(realname.rfind("."));
271 
272  OphRealField* realMat = new OphRealField[_wavelength_num];
273  OphRealField* imagMat = new OphRealField[_wavelength_num];
274 
275  std::string realtype = realname.substr(checktype + 1, realname.size());
276  std::string imgtype = imagname.substr(checktype + 1, realname.size());
277 
279 
280  if (realtype != imgtype) {
281  LOG("failed : The data type between real and imaginary is different!\n");
282  return false;
283  }
284  if (realtype == "bmp")
285  {
286  realname = real;
287  imagname = imag;
288 
289  oph::uchar* realdata = loadAsImg(realname.c_str());
290  oph::uchar* imagdata = loadAsImg(imagname.c_str());
291 
292  if (realdata == 0 && imagdata == 0) {
293  cout << "failed : hologram data load was failed." << endl;
294  return false;
295  }
296 
297  for (int z = 0; z < _wavelength_num; z++)
298  {
299  realMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
300  imagMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
301  for (int i = context_.pixel_number[_X] - 1; i >= 0; i--)
302  {
303  for (int j = 0; j < context_.pixel_number[_Y]; j++)
304  {
305  realMat[z](context_.pixel_number[_X] - i - 1, j) = (double)realdata[(i * context_.pixel_number[_Y] + j)*_wavelength_num + z];
306  imagMat[z](context_.pixel_number[_X] - i - 1, j) = (double)imagdata[(i * context_.pixel_number[_Y] + j)*_wavelength_num + z];
307  }
308  }
309  }
310  delete[] realdata;
311  delete[] imagdata;
312  }
313  else if (realtype == "bin")
314  {
315  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
316  double *imagdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
317 
318  for (int z = 0; z < _wavelength_num; z++)
319  {
320  realname = real;
321  imagname = imag;
322 
323  realname.insert(checktype, RGB_name[z]);
324  imagname.insert(checktype, RGB_name[z]);
325 
326  ifstream freal(realname.c_str(), ifstream::binary);
327  ifstream fimag(imagname.c_str(), ifstream::binary);
328 
329  freal.read(reinterpret_cast<char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
330  fimag.read(reinterpret_cast<char*>(imagdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
331 
332  realMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
333  imagMat[z].resize(context_.pixel_number[_X], context_.pixel_number[_Y]);
334 
335  for (int i = 0; i < context_.pixel_number[_X]; i++)
336  {
337  for (int j = 0; j < context_.pixel_number[_Y]; j++)
338  {
339  realMat[z](i, j) = realdata[i + j * context_.pixel_number[_X]];
340  imagMat[z](i, j) = imagdata[i + j * context_.pixel_number[_X]];
341  }
342  }
343  freal.close();
344  fimag.close();
345  }
346  delete[] realdata;
347  delete[] imagdata;
348  }
349  else
350  {
351  LOG("Error: wrong type\n");
352  }
353  //nomalization
354  double realout, imagout;
355  for (int z = 0; z < _wavelength_num; z++)
356  {
357  meanOfMat(realMat[z], realout); meanOfMat(imagMat[z], imagout);
358 
359  realMat[z] / realout; imagMat[z] / imagout;
360  absMat(realMat[z], realMat[z]);
361  absMat(imagMat[z], imagMat[z]);
362  realout = maxOfMat(realMat[z]); imagout = maxOfMat(imagMat[z]);
363  realMat[z] / realout; imagMat[z] / imagout;
364  realout = minOfMat(realMat[z]); imagout = minOfMat(imagMat[z]);
365  realMat[z] - realout; imagMat[z] - imagout;
366 
368 
369  for (int i = 0; i < context_.pixel_number[_X]; i++)
370  {
371  for (int j = 0; j < context_.pixel_number[_Y]; j++)
372  {
373  ComplexH[z](i, j)._Val[_RE] = realMat[z](i, j);
374  ComplexH[z](i, j)._Val[_IM] = imagMat[z](i, j);
375  }
376  }
377  }
378  LOG("Reading Openholo Complex Field File...%s, %s\n", realname.c_str(), imagname.c_str());
379 
380  return true;
381 }
382 
383 
384 
385 bool ophSig::save(const char *real, const char *imag)
386 {
387  string realname = real;
388  string imagname = imag;
389 
390  char* RGB_name[3] = { "","","" };
391 
392  if (_wavelength_num > 1) {
393  RGB_name[0] = "_B";
394  RGB_name[1] = "_G";
395  RGB_name[2] = "_R";
396  }
397 
398  int checktype = static_cast<int>(realname.rfind("."));
399  string type = realname.substr(checktype + 1, realname.size());
400  if (type == "bin")
401  {
402  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
403  double *imagdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
404 
405  for (int z = 0; z < _wavelength_num; z++)
406  {
407  realname = real;
408  imagname = imag;
409  realname.insert(checktype, RGB_name[z]);
410  imagname.insert(checktype, RGB_name[z]);
411  std::ofstream cos(realname.c_str(), std::ios::binary);
412  std::ofstream sin(imagname.c_str(), std::ios::binary);
413 
414  if (!cos.is_open()) {
415  LOG("real file not found.\n");
416  cos.close();
417  delete[] realdata;
418  delete[] imagdata;
419  return FALSE;
420  }
421 
422  if (!sin.is_open()) {
423  LOG("imag file not found.\n");
424  sin.close();
425  delete[] realdata;
426  delete[] imagdata;
427  return FALSE;
428  }
429 
430  for (int i = 0; i < context_.pixel_number[_X]; i++)
431  {
432  for (int j = 0; j < context_.pixel_number[_Y]; j++)
433  {
434  realdata[i + j * context_.pixel_number[_X]] = ComplexH[z](i, j)._Val[_RE];
435  imagdata[i + j * context_.pixel_number[_X]] = ComplexH[z](i, j)._Val[_IM];
436  }
437  }
438  cos.write(reinterpret_cast<const char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
439  sin.write(reinterpret_cast<const char*>(imagdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
440 
441  cos.close();
442  sin.close();
443  }
444  delete[]realdata;
445  delete[]imagdata;
446 
447  LOG("Writing Openholo Complex Field...%s, %s\n", realname.c_str(), imagname.c_str());
448  }
449  else if (type == "bmp")
450  {
451  oph::uchar* realdata;
452  oph::uchar* imagdata;
453  int _pixelbytesize = 0;
454  int _width = context_.pixel_number[_Y], _height = context_.pixel_number[_X];
455  int bitpixel = _wavelength_num * 8;
456 
457  if (bitpixel == 8)
458  {
459  _pixelbytesize = _height * _width;
460  }
461  else
462  {
463  _pixelbytesize = _height * _width * 3;
464  }
465  int _filesize = 0;
466 
467 
468  FILE *freal, *fimag;
469  fopen_s(&freal, realname.c_str(), "wb");
470  fopen_s(&fimag, imagname.c_str(), "wb");
471 
472  if ((freal == nullptr) || (fimag == nullptr))
473  {
474  LOG("file not found\n");
475  return FALSE;
476  }
477 
478  if (bitpixel == 8)
479  {
480  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height);
481  imagdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height);
482  _filesize = _pixelbytesize + sizeof(bitmap8bit);
483 
484  bitmap8bit *pbitmap = (bitmap8bit*)calloc(1, sizeof(bitmap8bit));
485  memset(pbitmap, 0x00, sizeof(bitmap8bit));
486 
487  pbitmap->fileheader.signature[0] = 'B';
488  pbitmap->fileheader.signature[1] = 'M';
489  pbitmap->fileheader.filesize = _filesize;
490  pbitmap->fileheader.fileoffset_to_pixelarray = sizeof(bitmap8bit);
491 
492  for (int i = 0; i < 256; i++) {
493  pbitmap->rgbquad[i].rgbBlue = i;
494  pbitmap->rgbquad[i].rgbGreen = i;
495  pbitmap->rgbquad[i].rgbRed = i;
496  }
497 
498 
500  for (int i = _height - 1; i >= 0; i--)
501  {
502  for (int j = 0; j < _width; j++)
503  {
504  if (ComplexH[0].mat[_height - i - 1][j]._Val[_RE] < 0)
505  {
506  ComplexH[0].mat[_height - i - 1][j]._Val[_RE] = 0;
507  }
508 
509  if (ComplexH[0].mat[_height - i - 1][j]._Val[_IM] < 0)
510  {
511  ComplexH[0].mat[_height - i - 1][j]._Val[_IM] = 0;
512  }
513  }
514  }
515 
516  double minVal, iminVal, maxVal, imaxVal;
517  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
518  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
519  if ((i == 0) && (j == 0))
520  {
521  minVal = ComplexH[0](i, j)._Val[_RE];
522  maxVal = ComplexH[0](i, j)._Val[_RE];
523  }
524  else {
525  if (ComplexH[0](i, j)._Val[_RE] < minVal)
526  {
527  minVal = ComplexH[0](i, j).real();
528  }
529  if (ComplexH[0](i, j)._Val[_RE] > maxVal)
530  {
531  maxVal = ComplexH[0](i, j).real();
532  }
533  }
534  if ((i == 0) && (j == 0)) {
535  iminVal = ComplexH[0](i, j)._Val[_IM];
536  imaxVal = ComplexH[0](i, j)._Val[_IM];
537  }
538  else {
539  if (ComplexH[0](i, j)._Val[_IM] < iminVal)
540  {
541  iminVal = ComplexH[0](i, j)._Val[_IM];
542  }
543  if (ComplexH[0](i, j)._Val[_IM] > imaxVal)
544  {
545  imaxVal = ComplexH[0](i, j)._Val[_IM];
546  }
547  }
548  }
549  }
550  for (int i = _height - 1; i >= 0; i--)
551  {
552  for (int j = 0; j < _width; j++)
553  {
554  realdata[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)._Val[_RE] - minVal) / (maxVal - minVal) * 255 + 0.5);
555  imagdata[i*_width + j] = (uchar)((ComplexH[0](_height - i - 1, j)._Val[_IM] - iminVal) / (imaxVal - iminVal) * 255 + 0.5);
556  }
557  }
558 
559  pbitmap->bitmapinfoheader.dibheadersize = sizeof(bitmapinfoheader);
560  pbitmap->bitmapinfoheader.width = _width;
561  pbitmap->bitmapinfoheader.height = _height;
562  pbitmap->bitmapinfoheader.planes = OPH_PLANES;
563  pbitmap->bitmapinfoheader.bitsperpixel = bitpixel;
564  pbitmap->bitmapinfoheader.compression = OPH_COMPRESSION;
565  pbitmap->bitmapinfoheader.imagesize = _pixelbytesize;
566  pbitmap->bitmapinfoheader.ypixelpermeter = 0;
567  pbitmap->bitmapinfoheader.xpixelpermeter = 0;
568  pbitmap->bitmapinfoheader.numcolorspallette = 256;
569 
570  fwrite(pbitmap, 1, sizeof(bitmap8bit), freal);
571  fwrite(realdata, 1, _pixelbytesize, freal);
572 
573  fwrite(pbitmap, 1, sizeof(bitmap8bit), fimag);
574  fwrite(imagdata, 1, _pixelbytesize, fimag);
575 
576  fclose(freal);
577  fclose(fimag);
578  free(pbitmap);
579  }
580  else
581  {
582  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height * bitpixel / 3);
583  imagdata = (oph::uchar*)malloc(sizeof(oph::uchar) * _width * _height * bitpixel / 3);
584  _filesize = _pixelbytesize + sizeof(fileheader) + sizeof(bitmapinfoheader);
585 
586  fileheader *hf = (fileheader*)calloc(1, sizeof(fileheader));
587  bitmapinfoheader *hInfo = (bitmapinfoheader*)calloc(1, sizeof(bitmapinfoheader));
588 
589  hf->signature[0] = 'B';
590  hf->signature[1] = 'M';
591  hf->filesize = _filesize;
592  hf->fileoffset_to_pixelarray = sizeof(fileheader) + sizeof(bitmapinfoheader);
593 
594  double minVal, iminVal, maxVal, imaxVal;
595  for (int z = 0; z < 3; z++)
596  {
597  for (int j = 0; j < ComplexH[0].size[_Y]; j++) {
598  for (int i = 0; i < ComplexH[0].size[_X]; i++) {
599  if ((i == 0) && (j == 0))
600  {
601  minVal = ComplexH[z](i, j)._Val[_RE];
602  maxVal = ComplexH[z](i, j)._Val[_RE];
603  }
604  else {
605  if (ComplexH[z](i, j)._Val[_RE] < minVal)
606  {
607  minVal = ComplexH[z](i, j)._Val[_RE];
608  }
609  if (ComplexH[z](i, j)._Val[_RE] > maxVal)
610  {
611  maxVal = ComplexH[z](i, j)._Val[_RE];
612  }
613  }
614  if ((i == 0) && (j == 0)) {
615  iminVal = ComplexH[z](i, j)._Val[_IM];
616  imaxVal = ComplexH[z](i, j)._Val[_IM];
617  }
618  else {
619  if (ComplexH[z](i, j)._Val[_IM] < iminVal)
620  {
621  iminVal = ComplexH[z](i, j)._Val[_IM];
622  }
623  if (ComplexH[z](i, j)._Val[_IM] > imaxVal)
624  {
625  imaxVal = ComplexH[z](i, j)._Val[_IM];
626  }
627  }
628  }
629  }
630 
631  for (int i = _height - 1; i >= 0; i--)
632  {
633  for (int j = 0; j < _width; j++)
634  {
635  realdata[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)._Val[_RE] - minVal) / (maxVal - minVal) * 255);
636  imagdata[3 * j + 3 * i * _width + z] = (uchar)((ComplexH[z](_height - i - 1, j)._Val[_IM] - iminVal) / (imaxVal - iminVal) * 255);
637 
638  }
639  }
640  }
641  hInfo->dibheadersize = sizeof(bitmapinfoheader);
642  hInfo->width = _width;
643  hInfo->height = _height;
644  hInfo->planes = OPH_PLANES;
645  hInfo->bitsperpixel = bitpixel;
646  hInfo->compression = OPH_COMPRESSION;
647  hInfo->imagesize = _pixelbytesize;
648  hInfo->ypixelpermeter = 0;
649  hInfo->xpixelpermeter = 0;
650 
651  fwrite(hf, 1, sizeof(fileheader), freal);
652  fwrite(hInfo, 1, sizeof(bitmapinfoheader), freal);
653  fwrite(realdata, 1, _pixelbytesize, freal);
654 
655  fwrite(hf, 1, sizeof(fileheader), fimag);
656  fwrite(hInfo, 1, sizeof(bitmapinfoheader), fimag);
657  fwrite(imagdata, 1, _pixelbytesize, fimag);
658 
659  fclose(freal);
660  fclose(fimag);
661  free(hf);
662  free(hInfo);
663  }
664 
665  free(realdata);
666  free(imagdata);
667  std::cout << "file save bmp complete\n" << endl;
668 
669  }
670  else {
671  LOG("failed : The Invalid data type! - %s\n", type);
672  }
673  return TRUE;
674 }
675 
676 bool ophSig::save(const char *fname)
677 {
678  string fullname = fname;
679 
680  char* RGB_name[3] = { "","","" };
681 
682  if (_wavelength_num > 1) {
683  RGB_name[0] = "_B";
684  RGB_name[1] = "_G";
685  RGB_name[2] = "_R";
686  }
687  int checktype = static_cast<int>(fullname.rfind("."));
688 
689  if (fullname.substr(checktype + 1, fullname.size()) == "bmp")
690  {
691  oph::uchar* realdata;
692  realdata = (oph::uchar*)malloc(sizeof(oph::uchar) * context_.pixel_number[_X] * context_.pixel_number[_Y] * _wavelength_num);
693 
694  double gamma = 0.5;
695  double maxIntensity = 0.0;
696  double realVal = 0.0;
697  double imagVal = 0.0;
698  double intensityVal = 0.0;
699 
700  for (int z = 0; z < _wavelength_num; z++)
701  {
702  for (int j = 0; j < context_.pixel_number[_Y]; j++) {
703  for (int i = 0; i < context_.pixel_number[_X]; i++) {
704  realVal = ComplexH[z](i, j)._Val[_RE];
705  imagVal = ComplexH[z](i, j)._Val[_RE];
706  intensityVal = realVal*realVal + imagVal*imagVal;
707  if (intensityVal > maxIntensity) {
708  maxIntensity = intensityVal;
709  }
710  }
711  }
712  for (int i = context_.pixel_number[_X] - 1; i >= 0; i--)
713  {
714  for (int j = 0; j < context_.pixel_number[_Y]; j++)
715  {
716  realVal = ComplexH[z](context_.pixel_number[_X] - i - 1, j)._Val[_RE];
717  imagVal = ComplexH[z](context_.pixel_number[_X] - i - 1, j)._Val[_IM];
718  intensityVal = realVal*realVal + imagVal*imagVal;
719  realdata[(i*context_.pixel_number[_Y] + j)* _wavelength_num + z] = (uchar)(pow(intensityVal / maxIntensity, gamma)*255.0);
720  }
721  }
722  //sprintf(str, "_%.2u", z);
723  //realname.insert(checktype, RGB_name[z]);
724  }
725  saveAsImg(fullname.c_str(), _wavelength_num * 8, realdata, context_.pixel_number[_X], context_.pixel_number[_Y]);
726 
727  delete[] realdata;
728  }
729  else if (fullname.substr(checktype + 1, fullname.size()) == "bin")
730  {
731  double *realdata = new double[context_.pixel_number[_X] * context_.pixel_number[_Y]];
732 
733  for (int z = 0; z < _wavelength_num; z++)
734  {
735  fullname = fname;
736  fullname.insert(checktype, RGB_name[z]);
737  std::ofstream cos(fullname.c_str(), std::ios::binary);
738 
739  if (!cos.is_open()) {
740  LOG("Error: file name not found.\n");
741  cos.close();
742  delete[] realdata;
743  return FALSE;
744  }
745 
746  for (int i = 0; i < context_.pixel_number[_X]; i++)
747  {
748  for (int j = 0; j < context_.pixel_number[_Y]; j++)
749  {
750  realdata[context_.pixel_number[_Y] * i + j] = ComplexH[z](i, j)._Val[_RE];
751  }
752  }
753  cos.write(reinterpret_cast<const char*>(realdata), sizeof(double) * context_.pixel_number[_X] * context_.pixel_number[_Y]);
754 
755  cos.close();
756  }
757  delete[] realdata;
758  }
759 
760  LOG("Writing Openholo Complex Field...%s\n", fullname.c_str());
761  return TRUE;
762 }
763 bool ophSig::sigConvertOffaxis(Real angleX, Real angleY) {
764  auto start_time = CUR_TIME;
765 
766  if (is_CPU == true)
767  {
768  std::cout << "Start Single Core CPU" << endl;
769  cvtOffaxis_CPU(angleX,angleY);
770  }
771  else {
772  std::cout << "Start Multi Core GPU" << std::endl;
773  cvtOffaxis_GPU(angleX, angleY);
774  }
775  auto end_time = CUR_TIME;
776  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
777  LOG("Implement time : %.5lf sec\n", during_time);
778  return true;
779 }
780 
781 bool ophSig::sigConvertHPO(Real depth, Real_t redRate) {
782  auto start_time = CUR_TIME;
783  if (is_CPU == true)
784  {
785  std::cout << "Start Single Core CPU" << endl;
786  sigConvertHPO_CPU(depth,redRate);
787 
788  }
789  else {
790  std::cout << "Start Multi Core GPU" << std::endl;
791 
792  sigConvertHPO_GPU(depth, redRate);
793 
794  }
795 
796  auto end_time = CUR_TIME;
797  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
798  LOG("Implement time : %.5lf sec\n", during_time);
799  return true;
800 }
801 
802 bool ophSig::sigConvertCAC(double red, double green, double blue){
803  auto start_time = CUR_TIME;
804  if (is_CPU == true)
805  {
806  std::cout << "Start Single Core CPU" << endl;
807  sigConvertCAC_CPU(red, green, blue);
808 
809  }
810  else {
811  std::cout << "Start Multi Core GPU" << std::endl;
812  sigConvertCAC_GPU(red, green, blue);
813 
814  }
815  auto end_time = CUR_TIME;
816  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
817  LOG("Implement time : %.5lf sec\n", during_time);
818  return true;
819 }
820 
821 bool ophSig::propagationHolo(float depth) {
822  auto start_time = CUR_TIME;
823  if (is_CPU == true)
824  {
825  std::cout << "Start Single Core CPU" << endl;
826  propagationHolo_CPU(depth);
827 
828  }
829  else {
830  std::cout << "Start Multi Core GPU" << std::endl;
831  propagationHolo_GPU(depth);
832 
833  }
834  auto end_time = CUR_TIME;
835  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
836  LOG("Implement time : %.5lf sec\n", during_time);
837  return true;
838 }
839 
840 double ophSig::sigGetParamSF(float zMax, float zMin, int sampN, float th) {
841  auto start_time = CUR_TIME;
842  double out = 0;
843  if (is_CPU == true)
844  {
845  std::cout << "Start Single Core CPU" << endl;
846  out = sigGetParamSF_CPU(zMax,zMin,sampN,th);
847 
848  }
849  else {
850  std::cout << "Start Multi Core GPU" << std::endl;
851  out = sigGetParamSF_GPU(zMax, zMin, sampN, th);
852 
853  }
854  auto end_time = CUR_TIME;
855  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
856  LOG("Implement time : %.5lf sec\n", during_time);
857  return out;
858 }
859 
860 
861 bool ophSig::cvtOffaxis_CPU(Real angleX, Real angleY) {
862 
863  int nx = context_.pixel_number[_X];
864  int ny = context_.pixel_number[_Y];
865  OphRealField H1(nx,ny);
866  Real x, y;
867  Complex<Real> F;
868  H1.resize(nx, ny);
869 
870  for (int i = 0; i < nx; i++)
871  {
872  y = (_cfgSig.height / (nx - 1)*i - _cfgSig.height / 2);
873  for (int j = 0; j < ny; j++)
874  {
875  x = (_cfgSig.width / (ny - 1)*j - _cfgSig.width / 2);
876 
877  //(*ComplexH)(i, j)._Val[_RE] = cos(((2 * M_PI) / *context_.wave_length)*(x*sin(angle[_X]) + y *sin(angle[_Y])));
878  //(*ComplexH)(i, j)._Val[_IM] = sin(((2 * M_PI) / *context_.wave_length)*(x*sin(angle[_X]) + y *sin(angle[_Y])));
879  F._Val[_RE] = cos(((2 * M_PI) / *context_.wave_length)*(x*sin(angleX) + y *sin(angleY)));
880  F._Val[_IM] = sin(((2 * M_PI) / *context_.wave_length)*(x*sin(angleX) + y *sin(angleY)));
881  H1(i, j) = ((*ComplexH)(i, j) * F)._Val[_RE];
882  }
883  }
884  double out = minOfMat(H1);
885  H1 - out;
886  out = maxOfMat(H1);
887  H1 / out;
888  //normalizeMat(H1, H1);
889 
890 
891  for (int i = 0; i < nx; i++)
892  {
893  for (int j = 0; j < ny; j++)
894  {
895  (*ComplexH)(i, j)._Val[_RE] = H1(i, j);
896  (*ComplexH)(i, j)._Val[_IM] = 0;
897  }
898  }
899 
900 
901 
902  return true;
903 }
904 
905 bool ophSig::sigConvertHPO_CPU(Real depth, Real_t redRate) {
906 
907 
908  int nx = context_.pixel_number[_X];
909  int ny = context_.pixel_number[_Y];
910 
911  Real wl = *context_.wave_length;
912  Real NA = _cfgSig.width/(2*depth);
913 
914  int xshift = nx / 2;
915  int yshift = ny / 2;
916 
917  Real y;
918 
919  Real_t NA_g = NA * redRate;
920 
921  OphComplexField F1(nx, ny);
922  OphComplexField FH(nx, ny);
923 
924 
925  Real Rephase = -(1 / (4 * M_PI)*pow((wl / NA_g), 2));
926  Real Imphase = ((1 / (4 * M_PI))*depth*wl);
927 
928  for (int i = 0; i < ny; i++)
929  {
930  int ii = (i + yshift) % ny;
931 
932  for (int j = 0; j < nx; j++)
933  {
934  y = (2 * M_PI * (j) / _cfgSig.height - M_PI * (nx - 1) / _cfgSig.height);
935  int jj = (j + xshift) % nx;
936  F1(jj, ii)._Val[_RE] = std::exp(Rephase*pow(y, 2))*cos(Imphase*pow(y, 2));
937  F1(jj, ii)._Val[_IM] = std::exp(Rephase*pow(y, 2))*sin(Imphase*pow(y, 2));
938  }
939  }
940  fft2((*ComplexH), FH, OPH_FORWARD);
941  F1.mulElem(FH);
942  fft2(F1, (*ComplexH), OPH_BACKWARD);
943 
944 
945 
946 
947  return true;
948 
949 }
950 
951 bool ophSig::sigConvertCAC_CPU(double red, double green, double blue) {
952 
953  Real x, y;
954  //OphComplexField exp, conj, FH_CAC;
955  int nx = context_.pixel_number[_X];
956  int ny = context_.pixel_number[_Y];
957 
958  if (_wavelength_num != 3) {
959  _wavelength_num = 3;
960  delete[] context_.wave_length;
962  }
963 
964  context_.wave_length[0] = blue;
965  context_.wave_length[1] = green;
966  context_.wave_length[2] = red;
967 
968  OphComplexField FFZP(nx, ny);
969  OphComplexField FH(nx, ny);
970 
971  for (int z = 0; z < _wavelength_num; z++)
972  {
973  double sigmaf = ((_foc[2] - _foc[z]) * context_.wave_length[z]) / (4 * M_PI);
974  int xshift = nx / 2;
975  int yshift = ny / 2;
976 
977  for (int i = 0; i < ny; i++)
978  {
979  int ii = (i + yshift) % ny;
980  y = (2 * M_PI * i) / _radius - (M_PI*(ny - 1)) / _radius;
981  for (int j = 0; j < nx; j++)
982  {
983  x = (2 * M_PI * j) / _radius - (M_PI*(nx - 1)) / _radius;
984 
985  int jj = (j + xshift) % nx;
986 
987  FFZP(jj, ii)._Val[_RE] = cos(sigmaf * (pow(x, 2) + pow(y, 2)));
988  FFZP(jj, ii)._Val[_IM] = -sin(sigmaf * (pow(x, 2) + pow(y, 2))); //conjugate 때문에 -붙음여
989  }
990  }
991  fft2(ComplexH[z], FH, OPH_FORWARD);
992  FH.mulElem(FFZP);
993  fft2(FH, ComplexH[z], OPH_BACKWARD);
994  }
995 
996  return true;
997 }
998 
999 bool ophSig::readConfig(const char* fname)
1000 {
1001  //LOG("Reading....%s...\n", fname);
1002 
1003  /*XML parsing*/
1004  tinyxml2::XMLDocument xml_doc;
1005  tinyxml2::XMLNode* xml_node;
1006 
1007  if (!checkExtension(fname, ".xml"))
1008  {
1009  LOG("file's extension is not 'xml'\n");
1010  return false;
1011  }
1012  auto ret = xml_doc.LoadFile(fname);
1013  if (ret != tinyxml2::XML_SUCCESS)
1014  {
1015  LOG("Failed to load file \"%s\"\n", fname);
1016  return false;
1017  }
1018 
1019  xml_node = xml_doc.FirstChild();
1020 
1021  (xml_node->FirstChildElement("pixel_number_x"))->QueryIntText(&context_.pixel_number[_X]);
1022  (xml_node->FirstChildElement("pixel_number_y"))->QueryIntText(&context_.pixel_number[_Y]);
1023  (xml_node->FirstChildElement("width"))->QueryFloatText(&_cfgSig.width);
1024  (xml_node->FirstChildElement("height"))->QueryFloatText(&_cfgSig.height);
1025  (xml_node->FirstChildElement("wavelength_num"))->QueryIntText(&_wavelength_num);
1026 
1028 
1029  (xml_node->FirstChildElement("wavelength"))->QueryDoubleText(context_.wave_length);
1030  (xml_node->FirstChildElement("NA"))->QueryFloatText(&_cfgSig.NA);
1031  (xml_node->FirstChildElement("z"))->QueryFloatText(&_cfgSig.z);
1032  (xml_node->FirstChildElement("radius_of_lens"))->QueryFloatText(&_radius);
1033  (xml_node->FirstChildElement("focal_length_R"))->QueryFloatText(&_foc[2]);
1034  (xml_node->FirstChildElement("focal_length_G"))->QueryFloatText(&_foc[1]);
1035  (xml_node->FirstChildElement("focal_length_B"))->QueryFloatText(&_foc[0]);
1036 
1037  return true;
1038 }
1039 
1040 
1041 
1042 bool ophSig::propagationHolo_CPU(float depth) {
1043  int i, j;
1044  Real x, y, sigmaf;
1045  int nx = context_.pixel_number[_X];
1046  int ny = context_.pixel_number[_Y];
1047 
1048  OphComplexField FH(nx, ny);
1049  //OphComplexField FFZP(nx, ny);
1050  int xshift = nx / 2;
1051  int yshift = ny / 2;
1052 
1053  for (int z = 0; z < _wavelength_num; z++) {
1054 
1055  sigmaf = (depth * context_.wave_length[z]) / (4 * M_PI);
1056 
1057  /*FH.resize(nx, ny);
1058  FFZP.resize(nx, ny);*/
1059 
1060  fft2(ComplexH[z], FH, OPH_FORWARD);
1061 
1062  for (i = 0; i < ny; i++)
1063  {
1064  int ii = (i + yshift) % ny;
1065  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1066 
1067  for (j = 0; j < nx; j++)
1068  {
1069  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1070  int jj = (j + xshift) % nx;
1071  double temp = FH(jj, ii)._Val[_RE];
1072  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];
1073  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];
1074 
1075  }
1076  }
1077 
1078  fft2(FH, ComplexH[z], OPH_BACKWARD);
1079  }
1080  return true;
1081 }
1082 
1084  int i, j;
1085  Real x, y, sigmaf;
1086  int nx = context_.pixel_number[_X];
1087  int ny = context_.pixel_number[_Y];
1088 
1089  OphComplexField FH(nx, ny);
1090  int xshift = nx / 2;
1091  int yshift = ny / 2;
1092 
1093 
1094  sigmaf = (depth * (*context_.wave_length)) / (4 * M_PI);
1095 
1096  fft2(complexH, FH);
1097 
1098  for (i = 0; i < ny; i++)
1099  {
1100  int ii = (i + yshift) % ny;
1101  y = (2 * M_PI * (i)) / _cfgSig.width - (M_PI*(ny - 1)) / (_cfgSig.width);
1102 
1103  for (j = 0; j < nx; j++)
1104  {
1105  x = (2 * M_PI * (j)) / _cfgSig.height - (M_PI*(nx - 1)) / (_cfgSig.height);
1106  int jj = (j + xshift) % nx;
1107  double temp = FH(jj, ii)._Val[_RE];
1108  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];
1109  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];
1110 
1111  }
1112  }
1113  fft2(FH, complexH, OPH_BACKWARD);
1114 
1115  return complexH;
1116 }
1117 
1119 {
1120  //auto start_time = CUR_TIME;
1121 
1122  Real index = 0;
1123  if (is_CPU == true)
1124  {
1125  //std::cout << "Start Single Core CPU" << endl;
1126  index = sigGetParamAT_CPU();
1127 
1128  }
1129  else {
1130  //std::cout << "Start Multi Core GPU" << std::endl;
1131  index = sigGetParamAT_GPU();
1132 
1133  }
1134 
1135  //auto end_time = CUR_TIME;
1136 
1137  //auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1138 
1139  //LOG("Implement time : %.5lf sec\n", during_time);
1140  return index;
1141 }
1142 
1144 
1145  int i = 0, j = 0;
1146  Real max = 0; Real index = 0;
1147  Real_t NA_g = (Real_t)0.025;
1148  int nx = context_.pixel_number[_X];
1149  int ny = context_.pixel_number[_Y];
1150 
1151  OphComplexField Flr(nx, ny);
1152  OphComplexField Fli(nx, ny);
1153  OphComplexField Hsyn(nx, ny);
1154 
1155  OphComplexField Fo(nx, ny);
1156  OphComplexField Fon, yn, Ab_yn;
1157 
1158  OphRealField Ab_yn_half;
1159  OphRealField G(nx, ny);
1160  Real x = 1, y = 1;
1161  vector<Real> t, tn;
1162 
1163  for (i = 0; i < nx; i++)
1164  {
1165  for (j = 0; j < ny; j++)
1166  {
1167 
1168  x = (2 * M_PI*(i) / _cfgSig.height - M_PI*(nx - 1) / _cfgSig.height);
1169  y = (2 * M_PI*(j) / _cfgSig.width - M_PI*(ny - 1) / _cfgSig.width);
1170  G(i, j) = std::exp(-M_PI * pow((*context_.wave_length) / (2 * M_PI * NA_g), 2) * (pow(y, 2) + pow(x, 2)));
1171  Flr(i, j)._Val[_RE] = (*ComplexH)(i, j)._Val[_RE];
1172  Fli(i, j)._Val[_RE] = (*ComplexH)(i, j)._Val[_IM];
1173  Flr(i, j)._Val[_IM] = 0;
1174  Fli(i, j)._Val[_IM] = 0;
1175  }
1176  }
1177 
1178  fft2(Flr, Flr);
1179  fft2(Fli, Fli);
1180 
1181  int xshift = nx / 2;
1182  int yshift = ny / 2;
1183 
1184  for (i = 0; i < nx; i++)
1185  {
1186  int ii = (i + xshift) % nx;
1187  for (j = 0; j < ny; j++)
1188  {
1189  int jj = (j + yshift) % ny;
1190  Hsyn(i, j)._Val[_RE] = Flr(i, j)._Val[_RE] * G(i, j);
1191  Hsyn(i, j)._Val[_IM] = Fli(i, j)._Val[_RE] * G(i, j);
1192  /*Hsyn_copy1(i, j) = Hsyn(i, j);
1193  Hsyn_copy2(i, j) = Hsyn_copy1(i, j) * Hsyn(i, j);
1194  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);
1195  Fo(ii, jj)._Val[_RE] = Hsyn_copy2(i, j)._Val[0] / Hsyn_copy3(i, j);
1196  Fo(ii, jj)._Val[_IM] = Hsyn_copy2(i, j)._Val[1] / Hsyn_copy3(i, j);*/
1197  Fo(ii, jj) = pow(Hsyn(i, j), 2) / (pow(abs(Hsyn(i, j)), 2) + pow(10, -300));
1198 
1199  }
1200  }
1201 
1202  t = linspace(0., 1., nx / 2 + 1);
1203  tn.resize(t.size());
1204  Fon.resize(1, t.size());
1205 
1206  for (int i = 0; i < tn.size(); i++)
1207  {
1208  tn.at(i) = pow(t.at(i), 0.5);
1209  Fon(0, i)._Val[_RE] = Fo(nx / 2 - 1, nx / 2 - 1 + i)._Val[_RE];
1210  Fon(0, i)._Val[_IM] = 0;
1211  }
1212  yn.resize(1, tn.size());
1213  linInterp(t, Fon, tn, yn);
1214  fft1(yn, yn);
1215  Ab_yn.resize(yn.size[_X], yn.size[_Y]);
1216  absMat(yn, Ab_yn);
1217  Ab_yn_half.resize(1, nx / 4 + 1);
1218 
1219  for (int i = 0; i < nx / 4 + 1; i++)
1220  {
1221  Ab_yn_half(0, i) = Ab_yn(0, nx / 4 + i - 1)._Val[_RE];
1222  if (i == 0) max = Ab_yn_half(0, 0);
1223  else
1224  {
1225  if (Ab_yn_half(0, i) > max)
1226  {
1227  max = Ab_yn_half(0, i);
1228  index = i;
1229  }
1230  }
1231  }
1232 
1233  index = -(((index + 1) - 120) / 10) / 140 + 0.1;
1234 
1235 
1236 
1237  return index;
1238 }
1239 
1240 double ophSig::sigGetParamSF_CPU(float zMax, float zMin, int sampN, float th) {
1241 
1242  int nx = context_.pixel_number[_X];
1243  int ny = context_.pixel_number[_Y];
1244 
1245  OphComplexField I(nx, ny);
1246  vector<Real> F;
1247  Real dz = (zMax - zMin) / sampN;
1248  Real f;
1249  Real_t z = 0;
1250  Real depth = 0;
1251  Real max = MIN_DOUBLE;
1252  int i, j, n = 0;
1253  Real ret1;
1254  Real ret2;
1255 
1256  for (n = 0; n < sampN + 1; n++)
1257  {
1258  z = ((n)* dz + zMin);
1259  f = 0;
1260  I = propagationHolo((*ComplexH), z);
1261 
1262  for (i = 0; i < nx - 2; i++)
1263  {
1264  for (j = 0; j < ny - 2; j++)
1265  {
1266  ret1 = abs(I(i + 2, j)._Val[_RE] - I(i, j)._Val[_RE]);
1267  ret2 = abs(I(i, j + 2)._Val[_RE] - I(i, j)._Val[_RE]);
1268  if (ret1 >= th) { f += ret1 * ret1; }
1269  else if (ret2 >= th) { f += ret2 * ret2; }
1270  }
1271  }
1272  //cout << (float)n / sampN * 100 << " %" << endl;
1273 
1274  if (f > max) {
1275  max = f;
1276  depth = z;
1277  }
1278  }
1279 
1280 
1281  return depth;
1282 }
1283 
1284 bool ophSig::getComplexHFromPSDH(const char * fname0, const char * fname90, const char * fname180, const char * fname270)
1285 {
1286  auto start_time = CUR_TIME;
1287  string fname0str = fname0;
1288  string fname90str = fname90;
1289  string fname180str = fname180;
1290  string fname270str = fname270;
1291  int checktype = static_cast<int>(fname0str.rfind("."));
1292  OphRealField f0Mat[3], f90Mat[3], f180Mat[3], f270Mat[3];
1293 
1294  std::string f0type = fname0str.substr(checktype + 1, fname0str.size());
1295 
1296  uint16_t bitsperpixel;
1297  fileheader hf;
1298  bitmapinfoheader hInfo;
1299 
1300  if (f0type == "bmp")
1301  {
1302  FILE *f0, *f90, *f180, *f270;
1303  fopen_s(&f0, fname0str.c_str(), "rb"); fopen_s(&f90, fname90str.c_str(), "rb");
1304  fopen_s(&f180, fname180str.c_str(), "rb"); fopen_s(&f270, fname270str.c_str(), "rb");
1305  if (!f0)
1306  {
1307  LOG("bmp file open fail! (phase shift = 0)\n");
1308  return false;
1309  }
1310  if (!f90)
1311  {
1312  LOG("bmp file open fail! (phase shift = 90)\n");
1313  return false;
1314  }
1315  if (!f180)
1316  {
1317  LOG("bmp file open fail! (phase shift = 180)\n");
1318  return false;
1319  }
1320  if (!f270)
1321  {
1322  LOG("bmp file open fail! (phase shift = 270)\n");
1323  return false;
1324  }
1325  fread(&hf, sizeof(fileheader), 1, f0);
1326  fread(&hInfo, sizeof(bitmapinfoheader), 1, f0);
1327 
1328  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
1329  if ((hInfo.height == 0) || (hInfo.width == 0))
1330  {
1331  LOG("bmp header is empty!\n");
1332  hInfo.height = context_.pixel_number[_X];
1333  hInfo.width = context_.pixel_number[_Y];
1334  if (hInfo.height == 0 || hInfo.width == 0)
1335  {
1336  LOG("check your parameter file!\n");
1337  return false;
1338  }
1339  }
1340  if ((context_.pixel_number[_Y] != hInfo.height) || (context_.pixel_number[_X] != hInfo.width)) {
1341  LOG("image size is different!\n");
1342  context_.pixel_number[_Y] = hInfo.height;
1343  context_.pixel_number[_X] = hInfo.width;
1344  LOG("changed parameter of size %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
1345  }
1346  bitsperpixel = hInfo.bitsperpixel;
1347  if (hInfo.bitsperpixel == 8)
1348  {
1349  _wavelength_num = 1;
1350  rgbquad palette[256];
1351  fread(palette, sizeof(rgbquad), 256, f0);
1352  fread(palette, sizeof(rgbquad), 256, f90);
1353  fread(palette, sizeof(rgbquad), 256, f180);
1354  fread(palette, sizeof(rgbquad), 256, f270);
1355 
1356  f0Mat[0].resize(hInfo.height, hInfo.width);
1357  f90Mat[0].resize(hInfo.height, hInfo.width);
1358  f180Mat[0].resize(hInfo.height, hInfo.width);
1359  f270Mat[0].resize(hInfo.height, hInfo.width);
1360  ComplexH = new OphComplexField;
1361  ComplexH[0].resize(hInfo.height, hInfo.width);
1362  }
1363  else
1364  {
1365  _wavelength_num = 3;
1366  ComplexH = new OphComplexField[3];
1367  f0Mat[0].resize(hInfo.height, hInfo.width);
1368  f90Mat[0].resize(hInfo.height, hInfo.width);
1369  f180Mat[0].resize(hInfo.height, hInfo.width);
1370  f270Mat[0].resize(hInfo.height, hInfo.width);
1371  ComplexH[0].resize(hInfo.height, hInfo.width);
1372 
1373  f0Mat[1].resize(hInfo.height, hInfo.width);
1374  f90Mat[1].resize(hInfo.height, hInfo.width);
1375  f180Mat[1].resize(hInfo.height, hInfo.width);
1376  f270Mat[1].resize(hInfo.height, hInfo.width);
1377  ComplexH[1].resize(hInfo.height, hInfo.width);
1378 
1379  f0Mat[2].resize(hInfo.height, hInfo.width);
1380  f90Mat[2].resize(hInfo.height, hInfo.width);
1381  f180Mat[2].resize(hInfo.height, hInfo.width);
1382  f270Mat[2].resize(hInfo.height, hInfo.width);
1383  ComplexH[2].resize(hInfo.height, hInfo.width);
1384  }
1385 
1386  uchar* f0data = (uchar*)malloc(sizeof(uchar)*hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8));
1387  uchar* f90data = (uchar*)malloc(sizeof(uchar)*hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8));
1388  uchar* f180data = (uchar*)malloc(sizeof(uchar)*hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8));
1389  uchar* f270data = (uchar*)malloc(sizeof(uchar)*hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8));
1390 
1391  fread(f0data, sizeof(uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), f0);
1392  fread(f90data, sizeof(uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), f90);
1393  fread(f180data, sizeof(uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), f180);
1394  fread(f270data, sizeof(uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), f270);
1395 
1396  fclose(f0);
1397  fclose(f90);
1398  fclose(f180);
1399  fclose(f270);
1400 
1401  for (int i = hInfo.height - 1; i >= 0; i--)
1402  {
1403  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
1404  {
1405  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1406  {
1407  f0Mat[z](hInfo.height - i - 1, j) = (double)f0data[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
1408  f90Mat[z](hInfo.height - i - 1, j) = (double)f90data[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
1409  f180Mat[z](hInfo.height - i - 1, j) = (double)f180data[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
1410  f270Mat[z](hInfo.height - i - 1, j) = (double)f270data[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
1411  }
1412  }
1413  }
1414  LOG("PSDH file load complete!\n");
1415 
1416  free(f0data);
1417  free(f90data);
1418  free(f180data);
1419  free(f270data);
1420 
1421  }
1422  else
1423  {
1424  LOG("wrong type (only BMP supported)\n");
1425  }
1426 
1427  // calculation complexH from 4 psdh and then normalize
1428  double normalizefactor = 1. / 256.;
1429  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
1430  {
1431  for (int i = 0; i < context_.pixel_number[_X]; i++)
1432  {
1433  for (int j = 0; j < context_.pixel_number[_Y]; j++)
1434  {
1435  ComplexH[z][j][i]._Val[_RE] = (f0Mat[z][j][i] - f180Mat[z][j][i])*normalizefactor;
1436  ComplexH[z][j][i]._Val[_IM] = (f90Mat[z][j][i] - f270Mat[z][j][i])*normalizefactor;
1437 
1438  }
1439  }
1440  }
1441  LOG("complex field obtained from 4 psdh\n");
1442 
1443  auto end_time = CUR_TIME;
1444 
1445  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
1446 
1447  LOG("Implement time : %.5lf sec\n", during_time);
1448 
1449  return true;
1450 }
1451 
1452 void ophSig::ophFree(void) {
1453 
1454 }
Real_t * _foc
Definition: ophSig.h:124
oph::matrix< Real > OphRealField
Definition: mat.h:419
#define OPH_BACKWARD
Definition: define.h:67
bool save(const char *real, const char *imag)
Save data as bmp or bin file.
Definition: ophSig.cpp:385
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:112
void cField2Buffer(matrix< Complex< Real >> &src, Complex< Real > **dst, int nx, int ny)
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:90
bool sigConvertOffaxis(Real angleX, Real angleY)
Function for Convert complex hologram to off-axis hologram.
Definition: ophSig.cpp:763
bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: ophSig.cpp:200
Real_t _radius
Definition: ophSig.h:123
Real * wave_length
Definition: Openholo.h:69
void meanOfMat(matrix< T > &src, T &dst)
Function for returns exponent e(x), where x is real number.
Definition: ophSig.h:291
bool cvtOffaxis_CPU(Real angleX, Real angleY)
Definition: ophSig.cpp:861
void ColorField2Buffer(matrix< Complex< Real >> &src, Complex< Real > **dst, int nx, int ny)
Definition: ophSig.cpp:73
unsigned char uchar
Definition: typedef.h:64
double sigGetParamAT_GPU()
Definition: ophSig_GPU.cpp:377
float Real
Definition: typedef.h:55
#define CUR_TIME
Definition: function.h:58
void fft1(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Definition: ophSig.cpp:115
double sigGetParamSF(float zMax, float zMin, int sampN, float th)
Extraction of distance parameter using sharpness functions.
Definition: ophSig.cpp:840
bool saveAsOhc(const char *fname)
Function to write OHC file
Definition: ophSig.cpp:229
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:80
bool propagationHolo(float depth)
Function for Chromatic aberration compensation filter.
Definition: ophSig.cpp:821
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
void setMode(bool is_CPU)
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
void getOHCheader(ohcHeader &_Header)
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:1284
#define MIN_DOUBLE
Definition: define.h:132
#define _Y
Definition: define.h:84
#define _IM
Definition: complex.h:57
bool readConfig(const char *fname)
Function for Read parameter.
Definition: ophSig.cpp:999
ImgEncoderOhc * OHC_encoder
OHC file format Variables for read and write.
Definition: Openholo.h:304
uint32_t imagesize
Definition: struct.h:63
void setFieldEncoding(const FldStore _fldStore, const FldCodeType _fldCodeType)
void getComplexFieldData(OphComplexField &cmplx_field, uint wavelen_idx)
Definition: ImgCodecOhc.h:81
ophSigConfig _cfgSig
Definition: ophSig.h:108
bool sigConvertHPO(Real depth, Real_t redRate)
Function for convert complex hologram to horizontal parallax only hologram.
Definition: ophSig.cpp:781
ImgDecoderOhc * OHC_decoder
Definition: Openholo.h:305
bool sigConvertCAC(double red, double green, double blue)
Function for Chromatic aberration compensation filter.
Definition: ophSig.cpp:802
#define _X
Definition: define.h:80
double sigGetParamAT()
Extraction of distance parameter using axis transfomation.
Definition: ophSig.cpp:1118
int _wavelength_num
Definition: ophSig.h:118
double sigGetParamAT_CPU()
Definition: ophSig.cpp:1143
Definition: struct.h:69
Real_t NA
Definition: ophSig.h:67
void addComplexFieldData(const OphComplexField &data)
OphComplexField * ComplexH
Definition: ophSig.h:109
uint16_t bitsperpixel
Definition: struct.h:61
void setNumOfWavlen(const uint n_wavlens)
bool sigConvertHPO_CPU(Real depth, Real_t redRate)
Definition: ophSig.cpp:905
oph::ivec2 pixel_number
Definition: Openholo.h:63
bool is_CPU
Definition: ophSig.h:107
bool sigConvertHPO_GPU(Real depth, Real_t redRate)
Definition: ophSig_GPU.cpp:118
#define _RE
Definition: complex.h:54
Real_t z
Definition: ophSig.h:68
Real_t height
Definition: ophSig.h:66
oph::matrix< Complex< Real > > OphComplexField
Definition: mat.h:421
bool sigConvertCAC_GPU(double red, double green, double blue)
Definition: ophSig_GPU.cpp:185
bool sigConvertCAC_CPU(double red, double green, double blue)
Definition: ophSig.cpp:951
double Real_t
Definition: typedef.h:56
uint32_t ypixelpermeter
Definition: struct.h:64
void absMat(matrix< Complex< T >> &src, matrix< T > &dst)
Function for Linear interpolation 1D.
Definition: ophSig.h:181
bool load(const char *real, const char *imag)
Load bmp or bin file.
Definition: ophSig.cpp:257
double sigGetParamSF_GPU(float zMax, float zMin, int sampN, float th)
Definition: ophSig_GPU.cpp:314
bool propagationHolo_CPU(float depth)
Definition: ophSig.cpp:1042
bool setFileName(const std::string &_fname)
Definition: ImgCodecOhc.cpp:90
bool propagationHolo_GPU(float depth)
Definition: ophSig_GPU.cpp:257
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2157
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)
Definition: ophSig.cpp:92
Real minOfMat(matrix< Real > &src)
Function for extracts min of matrix.
Definition: ophSig.h:336
virtual uchar * loadAsImg(const char *fname)
Function for loading image files.
Definition: Openholo.cpp:199
uint32_t height
Definition: struct.h:59
uint16_t planes
Definition: struct.h:60
Real maxOfMat(matrix< Real > &src)
Function for extracts max of matrix.
Definition: ophSig.h:307
void setWavelength(const Real _wavlen, const LenUnit _unit=LenUnit::m)
void setNumOfPixel(const uint _pxNumX, const uint _pxNumY)
#define OPH_COMPRESSION
Definition: define.h:146
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: ophSig.cpp:1452
ophSig(void)
Constructor.
Definition: ophSig.cpp:50
OphConfig context_
Definition: Openholo.h:297
Real_t width
Definition: ophSig.h:65
#define OPH_FORWARD
Definition: define.h:66
double sigGetParamSF_CPU(float zMax, float zMin, int sampN, float th)
Definition: ophSig.cpp:1240
uint8_t signature[2]
Definition: struct.h:51
#define OPH_PLANES
Definition: define.h:145
void cvtOffaxis_GPU(Real angleX, Real angleY)
Definition: ophSig_GPU.cpp:64
void getWavelength(std::vector< double_t > &wavlen_array)
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:160
unsigned int uint
Definition: typedef.h:62
uint32_t xpixelpermeter
Definition: struct.h:65
#define M_PI
Definition: define.h:52