Openholo  v4.2
Open Source Digital Holographic Library
ophRec.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 "ophRec.h"
47 #include "sys.h"
48 #include "function.h"
49 #include "tinyxml2.h"
50 #include "ImgControl.h"
51 #include <algorithm>
52 //#include "fftw3.h"
53 #include <omp.h>
54 
56  : Openholo()
57  , m_oldSimStep(0)
58  , m_nOldChannel(0)
59  , m_idx(0)
60 {
61 }
62 
64 {
65 }
66 #define IMAGE_VAL(x,y,c,w,n,mem) (mem[x*n + y*w*n + c])
67 vec3 image_sample(float xx, float yy, int c, size_t w, size_t h, double* in);
68 void circshift(Real* in, Real* out, int shift_x, int shift_y, int nx, int ny);
69 void circshift(Complex<Real>* in, Complex<Real>* out, int shift_x, int shift_y, int nx, int ny);
70 void ScaleBilnear(double* src, double* dst, int w, int h, int neww, int newh, double multiplyval = 1.0);
71 void reArrangeChannel(std::vector<double*>& src, double* dst, int pnx, int pny, int chnum);
72 void rotateCCW180(double* src, double* dst, int pnx, int pny, double mulival = 1.0);
73 
74 bool ophRec::readConfig(const char* fname)
75 {
76  bool bRet = true;
77  using namespace tinyxml2;
78  /*XML parsing*/
79  tinyxml2::XMLDocument xml_doc;
80  XMLNode *xml_node;
81 
82  if (!checkExtension(fname, ".xml"))
83  {
84  LOG("<FAILED> Wrong file ext.\n");
85  return false;
86  }
87 
88  auto ret = xml_doc.LoadFile(fname);
89  if (ret != XML_SUCCESS)
90  {
91  LOG("<FAILED> Loading file (%d)\n", ret);
92  return false;
93  }
94 
95  xml_node = xml_doc.FirstChild();
96 
97  int nWave = 1;
98  char szNodeName[32] = { 0, };
99  sprintf(szNodeName, "SLM_WaveNum");
100  auto next = xml_node->FirstChildElement(szNodeName); // OffsetInDepth
101  if (!next || XML_SUCCESS != next->QueryIntText(&nWave))
102  {
103  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
104  bRet = false;
105  }
106 
107  context_.waveNum = nWave;
109  context_.wave_length = new Real[nWave];
110 
111  for (int i = 1; i <= nWave; i++) {
112  sprintf(szNodeName, "SLM_WaveLength_%d", i);
113  next = xml_node->FirstChildElement(szNodeName);
114  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.wave_length[i - 1]))
115  {
116  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
117  bRet = false;
118  }
119  }
120 
121  sprintf(szNodeName, "SLM_PixelNumX");
122  next = xml_node->FirstChildElement(szNodeName);
123  if (!next || XML_SUCCESS != next->QueryIntText(&context_.pixel_number[_X]))
124  {
125  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
126  bRet = false;
127  }
128 
129  sprintf(szNodeName, "SLM_PixelNumY");
130  next = xml_node->FirstChildElement(szNodeName);
131  if (!next || XML_SUCCESS != next->QueryIntText(&context_.pixel_number[_Y]))
132  {
133  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
134  bRet = false;
135  }
136 
137  sprintf(szNodeName, "SLM_PixelPitchX");
138  next = xml_node->FirstChildElement("SLM_PixelPitchX");
139  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.pixel_pitch[_X]))
140  {
141  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
142  bRet = false;
143  }
144 
145  sprintf(szNodeName, "SLM_PixelPitchY");
146  next = xml_node->FirstChildElement("SLM_PixelPitchY");
147  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.pixel_pitch[_Y]))
148  {
149  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
150  bRet = false;
151  }
152 
153  // option
154  next = xml_node->FirstChildElement("IMG_Rotation");
155  if (!next || XML_SUCCESS != next->QueryBoolText(&imgCfg.rotate))
156  imgCfg.rotate = false;
157  next = xml_node->FirstChildElement("IMG_Merge");
158  if (!next || XML_SUCCESS != next->QueryBoolText(&imgCfg.merge))
159  imgCfg.merge = false;
160  next = xml_node->FirstChildElement("IMG_Flip");
161  if (!next || XML_SUCCESS != next->QueryIntText(&imgCfg.flip))
162  imgCfg.flip = 0;
163  next = xml_node->FirstChildElement("EyeLength");
164  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeLength))
165  bRet = false;
166  next = xml_node->FirstChildElement("EyePupilDiameter");
167  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyePupilDiaMeter))
168  bRet = false;
169  next = xml_node->FirstChildElement("EyeBoxSizeScaleFactor");
170  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeBoxSizeScale))
171  bRet = false;
172  next = xml_node->FirstChildElement("EyeBoxSizeX");
173  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeBoxSize[_X]))
174  bRet = false;
175  next = xml_node->FirstChildElement("EyeBoxSizeY");
176  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeBoxSize[_Y]))
177  bRet = false;
178  next = xml_node->FirstChildElement("EyeBoxUnit");
179  if (!next || XML_SUCCESS != next->QueryIntText(&rec_config.EyeBoxUnit))
180  bRet = false;
181  next = xml_node->FirstChildElement("EyeCenterX");
182  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeCenter[_X]))
183  bRet = false;
184  next = xml_node->FirstChildElement("EyeCenterY");
185  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeCenter[_Y]))
186  bRet = false;
187  next = xml_node->FirstChildElement("EyeCenterZ");
188  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeCenter[_Z]))
189  bRet = false;
190  next = xml_node->FirstChildElement("EyeFocusDistance");
191  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.EyeFocusDistance))
192  bRet = false;
193  next = xml_node->FirstChildElement("ResultSizeScale");
194  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.ResultSizeScale))
195  bRet = false;
196  next = xml_node->FirstChildElement("SimulationTo");
197  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.SimulationTo))
198  bRet = false;
199  next = xml_node->FirstChildElement("SimulationFrom");
200  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.SimulationFrom))
201  bRet = false;
202  next = xml_node->FirstChildElement("SimulationStep");
203  if (!next || XML_SUCCESS != next->QueryIntText(&rec_config.SimulationStep))
204  bRet = false;
205  next = xml_node->FirstChildElement("SimulationMode");
206  if (!next || XML_SUCCESS != next->QueryIntText(&rec_config.SimulationMode))
207  bRet = false;
208  next = xml_node->FirstChildElement("RatioAtRetina");
209  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.RatioAtRetina))
210  bRet = false;
211  next = xml_node->FirstChildElement("RatioAtPupil");
212  if (!next || XML_SUCCESS != next->QueryDoubleText(&rec_config.RatioAtPupil))
213  bRet = false;
214  next = xml_node->FirstChildElement("CreatePupilFieldImg");
215  if (!next || XML_SUCCESS != next->QueryBoolText(&rec_config.CreatePupilFieldImg))
216  bRet = false;
217  next = xml_node->FirstChildElement("CenteringRetinalImg");
218  if (!next || XML_SUCCESS != next->QueryBoolText(&rec_config.CenteringRetinaImg))
219  bRet = false;
220 
223 
224  rec_config.EyeBoxSize = rec_config.EyeBoxSizeScale * rec_config.EyePupilDiaMeter * rec_config.EyeBoxSize;
225  Initialize();
226 
227  LOG("**************************************************\n");
228  LOG(" Read Config (Common) \n");
229  LOG("1) SLM Number of Waves : %d\n", context_.waveNum);
230  for (uint i = 0; i < context_.waveNum; i++)
231  LOG(" 1-%d) SLM Wave length : %e\n", i + 1, context_.wave_length[i]);
232  LOG("2) SLM Resolution : %d x %d\n", context_.pixel_number[_X], context_.pixel_number[_Y]);
233  LOG("3) SLM Pixel Pitch : %e x %e\n", context_.pixel_pitch[_X], context_.pixel_pitch[_Y]);
234  LOG("4) Image Rotate : %s\n", imgCfg.rotate ? "Y" : "N");
235  LOG("5) Image Flip : %s\n", (imgCfg.flip == FLIP::NONE) ? "NONE" :
236  (imgCfg.flip == FLIP::VERTICAL) ? "VERTICAL" :
237  (imgCfg.flip == FLIP::HORIZONTAL) ? "HORIZONTAL" : "BOTH");
238  LOG("6) Image Merge : %s\n", imgCfg.merge ? "Y" : "N");
239  LOG("**************************************************\n");
240 
241  return bRet;
242 }
243 
244 bool ophRec::readImage(const char* path)
245 {
246  const int pnX = context_.pixel_number[_X];
247  const int pnY = context_.pixel_number[_Y];
248  const int N = pnX * pnY;
249 
250  int w, h, bytesperpixel;
251  bool ret = getImgSize(w, h, bytesperpixel, path);
252  uchar* imgload = new uchar[w * h * bytesperpixel];
253 
254  int nLine = ((w * bytesperpixel) + 3) & ~3;
255  int nSize = nLine * h;
256 
257  ret = loadAsImgUpSideDown(path, imgload);
258  if (!ret) {
259  LOG("Failed::Image Load: %s\n", path);
260  return false;
261  }
262  LOG("Succeed::Image Load: %s\n", path);
263 
264  uchar *tmp = new uchar[nSize];
265  memset(tmp, 0, sizeof(char) * nSize);
266 
267  if (w != pnX || h != pnY)
268  {
269  for (int ch = 0; ch < bytesperpixel; ch++)
270  imgScaleBilinear(imgload, tmp, w, h, pnX, pnY, ch);
271  }
272  else
273  memcpy(tmp, imgload, sizeof(char) * nSize);
274 
275  delete[] imgload;
276 
277  for (int ch = 0; ch < bytesperpixel; ch++)
278  {
279  int idx = bytesperpixel - 1 - ch;
280  for (int i = 0; i < N; i++)
281  {
282  int src = i * bytesperpixel;
283  complex_H[ch][i][_RE] = Real(tmp[src + idx]);
284  complex_H[ch][i][_IM] = 0.0;
285  }
286  }
287  delete[] tmp;
288  return true;
289 }
290 
291 bool ophRec::readImagePNA(const char *phase, const char *amplitude)
292 {
293  const int pnX = context_.pixel_number[_X];
294  const int pnY = context_.pixel_number[_Y];
295  const int N = pnX * pnY;
296 
297  int w, h, bytesperpixel;
298  bool ret = getImgSize(w, h, bytesperpixel, phase);
299  uchar* imgload = new uchar[w * h * bytesperpixel];
300 
301  int nLine = ((w * bytesperpixel) + 3) & ~3;
302  int nSize = nLine * h;
303 
304  ret = loadAsImgUpSideDown(phase, imgload);
305  if (!ret) {
306  LOG("Failed::Image Load: %s\n", phase);
307  return false;
308  }
309  LOG("Succeed::Image Load: %s\n", phase);
310 
311  uchar *phaseTmp = new uchar[nSize];
312  memset(phaseTmp, 0, sizeof(char) * nSize);
313 
314  if (w != pnX || h != pnY)
315  {
316  for (int ch = 0; ch < bytesperpixel; ch++)
317  imgScaleBilinear(imgload, phaseTmp, w, h, pnX, pnY, ch);
318  }
319  else
320  memcpy(phaseTmp, imgload, sizeof(char) * nSize);
321 
322  delete[] imgload;
323 
324 
325  ret = getImgSize(w, h, bytesperpixel, amplitude);
326  imgload = new uchar[w * h * bytesperpixel];
327 
328  nLine = ((w * bytesperpixel) + 3) & ~3;
329  nSize = nLine * h;
330 
331  ret = loadAsImgUpSideDown(amplitude, imgload);
332  if (!ret) {
333  LOG("Failed::Image Load: %s\n", amplitude);
334  return false;
335  }
336  LOG("Succeed::Image Load: %s\n", amplitude);
337 
338  uchar *ampTmp = new uchar[nSize];
339  memset(ampTmp, 0, sizeof(char) * nSize);
340 
341  if (w != pnX || h != pnY)
342  {
343  for (int ch = 0; ch < bytesperpixel; ch++)
344  imgScaleBilinear(imgload, ampTmp, w, h, pnX, pnY, ch);
345  }
346  else
347  memcpy(ampTmp, imgload, sizeof(char) * nSize);
348 
349  delete[] imgload;
350 
351  Real PI2 = M_PI * 2;
352 
353  for (int ch = 0; ch < bytesperpixel; ch++)
354  {
355  int idx = bytesperpixel - 1 - ch;
356  for (int i = 0; i < N; i++)
357  {
358  int src = i * bytesperpixel;
359  Real p = Real(phaseTmp[src + idx]);
360  Real a = Real(ampTmp[src + idx]);
361  p = p / 255.0 * PI2 - M_PI; // -pi ~ pi
362  a = a / 255.0; // 0 ~ 1
363  Complex<Real> tmp(0, p);
364  tmp.exp();
365  complex_H[ch][i] = a * tmp;
366 
367  }
368  }
369 
370  delete[] phaseTmp;
371  delete[] ampTmp;
372  return true;
373 }
374 
375 bool ophRec::readImageRNI(const char *real, const char *imag)
376 {
377  const int pnX = context_.pixel_number[_X];
378  const int pnY = context_.pixel_number[_Y];
379  const int N = pnX * pnY;
380 
381  int w, h, bytesperpixel;
382  bool ret = getImgSize(w, h, bytesperpixel, real);
383  uchar* imgload = new uchar[w * h * bytesperpixel];
384 
385  int nLine = ((w * bytesperpixel) + 3) & ~3;
386  int nSize = nLine * h;
387 
388  ret = loadAsImgUpSideDown(real, imgload);
389  if (!ret) {
390  LOG("Failed::Image Load: %s\n", real);
391  return false;
392  }
393  LOG("Succeed::Image Load: %s\n", real);
394 
395  uchar *realTmp = new uchar[nSize];
396  memset(realTmp, 0, sizeof(char) * nSize);
397 
398  if (w != pnX || h != pnY)
399  {
400  for (int ch = 0; ch < bytesperpixel; ch++)
401  imgScaleBilinear(imgload, realTmp, w, h, pnX, pnY, ch);
402  }
403  else
404  memcpy(realTmp, imgload, sizeof(char) * nSize);
405 
406  delete[] imgload;
407 
408 
409  ret = getImgSize(w, h, bytesperpixel, imag);
410  imgload = new uchar[w * h * bytesperpixel];
411  nLine = ((w * bytesperpixel) + 3) & ~3;
412  nSize = nLine * h;
413 
414  ret = loadAsImgUpSideDown(imag, imgload);
415  if (!ret) {
416  LOG("Failed::Image Load: %s\n", imag);
417  return false;
418  }
419  LOG("Succeed::Image Load: %s\n", imag);
420 
421  uchar *imagTmp = new uchar[nSize];
422  memset(imagTmp, 0, sizeof(char) * nSize);
423 
424  if (w != pnX || h != pnY)
425  {
426  for (int ch = 0; ch < bytesperpixel; ch++)
427  imgScaleBilinear(imgload, imagTmp, w, h, pnX, pnY, ch);
428  }
429  else
430  memcpy(imagTmp, imgload, sizeof(char) * nSize);
431 
432  delete[] imgload;
433 
434  for (int ch = 0; ch < bytesperpixel; ch++)
435  {
436  int idx = bytesperpixel - 1 - ch;
437  for (int i = 0; i < N; i++)
438  {
439  int src = i * bytesperpixel;
440  complex_H[ch][i][_RE] = (Real)realTmp[src + idx] / 255.0;
441  complex_H[ch][i][_IM] = (Real)imagTmp[src + idx] / 255.0;
442  }
443  }
444  delete[] realTmp;
445  delete[] imagTmp;
446  return true;
447 }
448 
450 {
451  LOG("Propagation to observer plane\n");
452 
453  const int channel = context_.waveNum;
454  const int pnX = context_.pixel_number[_X];
455  const int pnY = context_.pixel_number[_Y];
456  field_set_.resize(channel);
457  pp_set_.resize(channel);
458  fftw_complex *in = nullptr;
459  fftw_complex *out = nullptr;
460  fftw_plan fftw_plan_fwd = fftw_plan_dft_2d(pnY, pnX, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
461 
462  for (int ch = 0; ch < channel; ch++)
463  {
465  }
466 
467  fftw_destroy_plan(fftw_plan_fwd);
468 }
469 
471 {
472  LOG("Propagation to observer plane\n");
473 
474  const int nChannel = context_.waveNum;
475  const int pnX = context_.pixel_number[_X];
476  const int pnY = context_.pixel_number[_Y];
477  const int N = pnX * pnY;
478  const Real ppX = context_.pixel_pitch[_X];
479  const Real ppY = context_.pixel_pitch[_Y];
480 
481  field_set_.resize(nChannel);
482  pn_set_.resize(nChannel);
483  pp_set_.resize(nChannel);
484 
485  fftw_complex *in = nullptr, *out = nullptr;
486  fftw_plan fft_plan_fwd = fftw_plan_dft_2d(pnY, pnX, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
487 
488  double prop_z = rec_config.EyeCenter[_Z];
489  double f_field = rec_config.EyeCenter[_Z];
490 
491 
492  Real ssX = pnX * ppX;
493  Real ssY = pnY * ppY;
494 
495 
496  for (int ctr = 0; ctr < nChannel; ctr++)
497  {
498  LOG("Color Number: %d\n", ctr + 1);
499 
500  Complex<Real>* u = complex_H[ctr];
501  field_set_[ctr] = new Complex<Real>[N];
502  memset(field_set_[ctr], 0.0, sizeof(Complex<Real>) * N);
503 
504  Real lambda = context_.wave_length[ctr];
505 
506  Real k = 2 * M_PI / lambda;
507  Real kk = k / (prop_z * 2);
508  Real kpropz = k * prop_z;
509  Real lambdapropz = lambda * prop_z;
510  Real ss_res_x = fabs(lambdapropz / ppX);
511  Real ss_res_y = fabs(lambdapropz / ppY);
512  Real hss_res_x = ss_res_x / 2.0;
513  Real hss_res_y = ss_res_y / 2.0;
514  Real pp_res_x = ss_res_x / Real(pnX);
515  Real pp_res_y = ss_res_y / Real(pnY);
516  Real absppX = fabs(lambdapropz / (4 * ppX));
517  Real absppY = fabs(lambdapropz / (4 * ppY));
518 
519  int loopi;
520 #pragma omp parallel for private(loopi)
521  for (loopi = 0; loopi < N; loopi++)
522  {
523  int x = loopi % pnX;
524  int y = loopi / pnX;
525 
526  double xx_src = (-ssX / 2.0) + (ppX * double(x));
527  double yy_src = (ssY / 2.0) - ppY - (ppY * double(y));
528 
529  if (f_field != prop_z)
530  {
531  double effective_f = f_field * prop_z / (f_field - prop_z);
532  double sval = (xx_src*xx_src) + (yy_src*yy_src);
533  sval *= M_PI / lambda / effective_f;
534  Complex<Real> kernel(0, sval);
535  kernel.exp();
536  //exponent_complex(&kernel);
537 
538  double anti_aliasing_mask = fabs(xx_src) < fabs(lambda*effective_f / (4 * ppX)) ? 1 : 0;
539  anti_aliasing_mask *= fabs(yy_src) < fabs(lambda*effective_f / (4 * ppY)) ? 1 : 0;
540 
541  field_set_[ctr][x + y * pnX] = u[x + y * pnX] * kernel * anti_aliasing_mask;
542 
543 
544  }
545  else {
546 
547  field_set_[ctr][x + y * pnX] = u[x + y * pnX];
548 
549  }
550  }
551  //delete[] fringe_[ctr];
552  //free(fringe_[ctr]);
553 
554  fft2(field_set_[ctr], field_set_[ctr], pnX, pnY, FFTW_FORWARD, false);
555 
556 
557 
558 #pragma omp parallel for private(loopi)
559  for (loopi = 0; loopi < N; loopi++)
560  {
561  int x = loopi % pnX;
562  int y = loopi / pnX;
563 
564  double xx_res = (-ss_res_x / 2.0) + (pp_res_x * double(x));
565  double yy_res = (ss_res_y / 2.0) - pp_res_y - (pp_res_y * double(y));
566 
567  Complex<Real> tmp1(0, k*prop_z);
568  tmp1.exp();
569  //exponent_complex(&tmp1);
570 
571  Complex<Real> tmp2(0, lambda*prop_z);
572 
573  double ssval = (xx_res*xx_res) + (yy_res*yy_res);
574  ssval *= k / (2 * prop_z);
575  Complex<Real> tmp3(0, ssval);
576  //exponent_complex(&tmp3);
577  tmp3.exp();
578 
579  Complex<Real> coeff = tmp1 / tmp2 * tmp3;
580 
581  field_set_[ctr][x + y * pnX] *= coeff;
582 
583  }
584  pp_set_[ctr][0] = pp_res_x;
585  pp_set_[ctr][1] = pp_res_y;
586 
587  }
588 
589  fftw_destroy_plan(fft_plan_fwd);
590 }
591 
593 {
594  const int pnX = context_.pixel_number[_X];
595  const int pnY = context_.pixel_number[_Y];
596  const int N = pnX * pnY;
597  const int nWave = context_.waveNum;
598 
599  const Real ppX = context_.pixel_pitch[_X];
600  const Real ppY = context_.pixel_pitch[_Y];
601 
602  const Real simFrom = rec_config.SimulationFrom;
603  const Real simTo = rec_config.SimulationTo;
604  const int simStep = rec_config.SimulationStep;
605  const Real simGap = (simStep > 1) ? (simTo - simFrom) / (simStep - 1) : 0;
606 
607  const Real tx = 1 / ppX;
608  const Real ty = 1 / ppY;
609  const Real dx = tx / pnX;
610  const Real dy = ty / pnY;
611 
612  const Real htx = tx / 2;
613  const Real hty = ty / 2;
614  const Real hdx = dx / 2;
615  const Real hdy = dy / 2;
616  const Real baseX = -htx + hdx;
617  const Real baseY = -hty + hdy;
618 
619  Complex<Real>* tmp = new Complex<Real>[N];
620  Complex<Real>* src = nullptr;
621  Complex<Real>* dst = new Complex<Real>[N];
622  Complex<Real>** kernels = new Complex<Real>*[nWave];// [N];
623  for (int i = 0; i < nWave; i++)
624  {
625  kernels[i] = new Complex<Real>[N];
626  }
627 
628  LOG("%s : Get Spatial Kernel\n", __FUNCTION__);
629  auto begin = CUR_TIME;
630  for (int ch = 0; ch < nWave; ch++)
631  {
632  const Real lambda = context_.wave_length[ch];
633  const Real k = 2 * M_PI / lambda;
634  src = complex_H[ch];
635 
636  // Get Spatial Kernel
637  int i;
638 #ifdef _OPENMP
639 #pragma omp parallel for private(i) firstprivate(pnX, lambda, dx, dy, baseX, baseY)
640 #endif
641  for (i = 0; i < N; i++)
642  {
643  int x = i % pnX;
644  int y = i / pnX;
645 
646  Real curX = baseX + (x * dx);
647  Real curY = baseY + (y * dy);
648  Real xx = curX * lambda;
649  Real yy = curY * lambda;
650  Real powxx = xx * xx;
651  Real powyy = yy * yy;
652 
653  kernels[ch][i][_RE] = 0;
654  kernels[ch][i][_IM] = sqrt(1 - powxx - powyy);
655  }
656  }
657  LOG(" => %lf(s)\n", ELAPSED_TIME(begin, CUR_TIME));
658  LOG("%s : Simultation\n", __FUNCTION__);
659  begin = CUR_TIME;
660 
661  for (int step = 0; step < simStep; step++)
662  {
663  Real min = MAX_DOUBLE, max = MIN_DOUBLE;
664  for (int ch = 0; ch < nWave; ch++)
665  {
666  const Real lambda = context_.wave_length[ch];
667  const Real k = 2 * M_PI / lambda;
668  src = complex_H[ch];
669  fft2(src, tmp, pnX, pnY, FFTW_FORWARD, false);
670  Real z = simFrom + (step * simGap);
671  Real kz = k * z;
672 
673  Real* encode = new Real[N];
674  uchar* normal = new uchar[N];
675 
676  //m_vecEncoded[step] = new Real[N];
677  //m_vecNormalized[step] = new uchar[N];
678 
679  for (int i = 0; i < N; i++)
680  {
681  Complex<Real> kernel = kernels[ch][i];
682  kernel[_IM] *= kz;
683  kernel.exp();
684  tmp[i] *= kernel;
685  }
686 
687  fft2(tmp, dst, pnX, pnY, FFTW_BACKWARD, true);
688 
689  for (int i = 0; i < N; i++)
690  {
691  encode[i] = dst[i].mag();
692 
693  if (min > encode[i])
694  min = encode[i];
695  if (max < encode[i])
696  max = encode[i];
697  }
698 
699  m_vecEncoded.push_back(encode);
700  m_vecNormalized.push_back(normal);
701  }
702 
703  LOG("step: %d => max: %e / min: %e\n", step, max, min);
704  if (nWave == 3)
705  {
706  for (int ch = 0; ch < nWave; ch++)
707  {
708  int idx = step * nWave + ch;
709  normalize(m_vecEncoded[idx], m_vecNormalized[idx], pnX, pnY, max, min);
710  }
711  }
712  else
713  normalize(m_vecEncoded[step], m_vecNormalized[step], pnX, pnY);
714  }
715 
716  m_oldSimStep = simStep;
717  fftFree();
718  LOG(" => %lf(s)\n", ELAPSED_TIME(begin, CUR_TIME));
719 
720  for (int i = 0; i < nWave; i++)
721  delete[] kernels[i];
722  delete[] kernels;
723  delete[] tmp;
724  delete[] dst;
725 
726 
727 }
728 
730 {
731  LOG("Color Number: %d\n", chnum + 1);
732 
733  const int pnX = context_.pixel_number[_X];
734  const int pnY = context_.pixel_number[_Y];
735  const int N = pnX * pnY;
736  Real ppX = context_.pixel_pitch[_X];
737  Real ppY = context_.pixel_pitch[_Y];
738  Real ppX4 = ppX * 4;
739  Real ppY4 = ppY * 4;
740  Real lambda = context_.wave_length[chnum];
741  Real ssX = pnX * ppX;
742  Real ssY = pnY * ppY;
743  Real hssX = ssX / 2.0;
744  Real hssY = ssY / 2.0;
745  Real prop_z = rec_config.EyeCenter[_Z];
746  vec2 cxy = vec2(0);
747 
748  Real k = 2 * M_PI / lambda;
749  Real kk = k / (prop_z * 2);
750  Real kpropz = k * prop_z;
751  Real lambdapropz = lambda * prop_z;
752  Real ss_res_x = fabs(lambdapropz / ppX);
753  Real ss_res_y = fabs(lambdapropz / ppY);
754  Real hss_res_x = ss_res_x / 2.0;
755  Real hss_res_y = ss_res_y / 2.0;
756  Real pp_res_x = ss_res_x / Real(pnX);
757  Real pp_res_y = ss_res_y / Real(pnY);
758 
759  Complex<Real>* tmp = complex_H[chnum];
760  field_set_[chnum] = new Complex<Real>[N];
761  memset(field_set_[chnum], 0.0, sizeof(Complex<Real>) * N);
762 
763  //Real xx_src, yy_src, xx_res, yy_res;
764  //int x, y;
765 
766  Real absppX = fabs(lambdapropz / (ppX4));
767  Real absppY = fabs(lambdapropz / (ppY4));
768 
769  int i;
770 #pragma omp parallel for private(i)
771  for (i = 0; i < N; i++)
772  {
773  int x = i % pnX;
774  int y = i / pnX;
775 
776  Real xx_src = -hssX + (ppX * Real(x));
777  Real yy_src = hssY - ppY - (ppY * Real(y));
778  Real xxx = xx_src - cxy[_X];
779  Real yyy = yy_src - cxy[_Y];
780 
781  Real sval = xxx * xxx + yyy * yyy;
782  sval *= kk;
783  Complex<Real> kernel(0, sval);
784  kernel.exp();
785  //exponent_complex(&kernel);
786 
787  Real anti_aliasing_mask = fabs(xxx) < absppX ? 1 : 0;
788  anti_aliasing_mask *= fabs(yyy) < absppY ? 1 : 0;
789 
790  field_set_[chnum][x + y * pnX] = tmp[x + y * pnX] * kernel * anti_aliasing_mask;
791 
792  }
793  //free(fringe_[chnum]);
794 
795  //fftw_complex *in = nullptr, *out = nullptr;
796  fft2(field_set_[chnum], field_set_[chnum], pnX, pnY, FFTW_FORWARD, false);
797 
798 #pragma omp parallel for private(i)
799  for (i = 0; i < N; i++)
800  {
801  int x = i % pnX;
802  int y = i / pnX;
803 
804  Real xx_res = -hss_res_x + (pp_res_x * Real(x));
805  Real yy_res = hss_res_y - pp_res_y - (pp_res_y * Real(y));
806  Real xxx = xx_res - cxy[_X];
807  Real yyy = yy_res - cxy[_Y];
808 
809  Complex<Real> tmp1(0, kpropz);
810  tmp1.exp();
811  //exponent_complex(&tmp1);
812 
813  Complex<Real> tmp2(0, lambdapropz);
814 
815  Real ssval = xxx * xxx + yyy * yyy;
816  ssval *= kk;
817  Complex<Real> tmp3(0, ssval);
818  tmp3.exp();
819  //exponent_complex(&tmp3);
820 
821  Complex<Real> coeff = tmp1 / tmp2 * tmp3;
822 
823  field_set_[chnum][x + y * pnX] *= coeff;
824 
825  }
826 
827  //pn_set_[chnum] = PIXEL_NUMBER;
828  pp_set_[chnum][_X] = pp_res_x;
829  pp_set_[chnum][_Y] = pp_res_y;
830 }
831 
833 {
834  LOG("Simulation start\n");
835  FILE *fp;
836  const int pnX = context_.pixel_number[_X];
837  const int pnY = context_.pixel_number[_Y];
838  const Real ppX = context_.pixel_pitch[_X];
839  const Real ppY = context_.pixel_pitch[_Y];
840  const int nChannel = context_.waveNum;
841  const int simStep = rec_config.SimulationStep;
842  const Real simFrom = rec_config.SimulationFrom;
843  const Real simTo = rec_config.SimulationTo;
844  const int simMode = rec_config.SimulationMode;
845  vec2 boxSize = rec_config.EyeBoxSize;
846  vec3 eyeCenter = rec_config.EyeCenter;
847  Real eyeLen = rec_config.EyeLength;
848  Real eyeFocusDistance = rec_config.EyeFocusDistance;
849  bool bCreatePupilImg = rec_config.CreatePupilFieldImg;
850  bool bCenteringRetinaImg = rec_config.CenteringRetinaImg;
851  Real resultSizeScale = rec_config.ResultSizeScale;
852  Real eyePupil = rec_config.EyePupilDiaMeter;
853  bool bSimPos[3];
854  Real ratioAtRetina = rec_config.RatioAtRetina;
855 
856  memcpy(&bSimPos[0], rec_config.SimulationPos, sizeof(bSimPos));
857 
858  vec3 var_step;
859  std::vector<vec3> var_vals;
860  var_vals.resize(simStep);
861 
862  if (simStep > 1)
863  {
864  var_step = (simTo - simFrom) / (simStep - 1);
865  for (int i = 0; i < simStep; i++)
866  var_vals[i] = simFrom + var_step * i;
867 
868  }
869  else
870  var_vals[0] = (simTo + simFrom) / 2.0;
871 
872  Real lambda, k;
873  int pn_e_x, pn_e_y;
874  Real pp_e_x, pp_e_y, ss_e_x, ss_e_y;
875  int pn_p_x, pn_p_y;
876  Real pp_p_x, pp_p_y, ss_p_x, ss_p_y;
877  vec2 eye_box_range_x, eye_box_range_y;
878  ivec2 eye_box_range_idx, eye_box_range_idy;
879  ivec2 eye_shift_by_pn;
880 
881  field_ret_set_.resize(nChannel);
882  pn_ret_set_.resize(nChannel);
883  pp_ret_set_.resize(nChannel);
884  ss_ret_set_.resize(nChannel);
885  recon_set.resize(simStep * nChannel);
886  img_set.resize(simStep * nChannel);
887  img_size.resize(simStep * nChannel);
888  focus_recon_set.resize(simStep);
889  focus_img_set.resize(simStep);
890  focus_img_size.resize(simStep);
891 
892  std::string varname2;
893 
894  for (int vtr = 0; vtr < simStep; vtr++)
895  {
896  if (simMode == 0) {
897  eyeFocusDistance = var_vals[vtr][_X];
898  }
899  else {
900  if (bSimPos[_X]) eyeCenter[_X] = var_vals[vtr][_X];
901  if (bSimPos[_Y]) eyeCenter[_Y] = var_vals[vtr][_X];
902  if (bSimPos[_Z]) eyeCenter[_Z] = var_vals[vtr][_X];
903  }
904 
905  for (int ctr = 0; ctr < nChannel; ctr++)
906  {
907  lambda = context_.wave_length[ctr];
908  k = 2 * M_PI / lambda;
909 
910  pn_e_x = pnX;
911  pn_e_y = pnY;
912  pp_e_x = pp_set_[ctr][_X];
913  pp_e_y = pp_set_[ctr][_Y];
914  ss_e_x = Real(pn_e_x) * pp_e_x;
915  ss_e_y = Real(pn_e_y) * pp_e_y;
916 
917  eye_shift_by_pn[0] = round(eyeCenter[_X] / pp_e_x);
918  eye_shift_by_pn[1] = round(eyeCenter[_Y] / pp_e_y);
919 
920  Complex<Real>* hh_e_shift = new Complex<Real>[pn_e_x * pn_e_y];
921  memset(hh_e_shift, 0.0, sizeof(Complex<Real>) * pn_e_x * pn_e_y);
922  circshift(field_set_[ctr], hh_e_shift, -eye_shift_by_pn[0], eye_shift_by_pn[1], pn_e_x, pn_e_y);
923 
924 
925  if (rec_config.EyeBoxUnit == 0)
926  {
927  eye_box_range_x[0] = -boxSize[_X] / 2.0;
928  eye_box_range_x[1] = boxSize[_X] / 2.0;
929  eye_box_range_y[0] = -boxSize[_Y] / 2.0;
930  eye_box_range_y[1] = boxSize[_Y] / 2.0;
931  eye_box_range_idx[0] = floor((eye_box_range_x[0] + ss_e_x / 2.0) / pp_e_x);
932  eye_box_range_idx[1] = floor((eye_box_range_x[1] + ss_e_x / 2.0) / pp_e_x);
933  eye_box_range_idy[0] = pn_e_y - floor((eye_box_range_y[1] + ss_e_y / 2.0) / pp_e_y);
934  eye_box_range_idy[1] = pn_e_y - floor((eye_box_range_y[0] + ss_e_y / 2.0) / pp_e_y);
935 
936  }
937  else {
938 
939  int temp = floor((pn_e_x - boxSize[_X]) / 2.0);
940  eye_box_range_idx[0] = temp;
941  eye_box_range_idx[1] = temp + boxSize[_X] - 1;
942  temp = floor((pn_e_y - boxSize[_Y]) / 2.0);
943  eye_box_range_idy[0] = temp;
944  eye_box_range_idy[1] = temp + boxSize[_Y] - 1;
945  }
946 
947  pn_p_x = eye_box_range_idx[1] - eye_box_range_idx[0] + 1;
948  pn_p_y = eye_box_range_idy[1] - eye_box_range_idy[0] + 1;
949  pp_p_x = pp_e_x;
950  pp_p_y = pp_e_y;
951  ss_p_x = pn_p_x * pp_p_x;
952  ss_p_y = pn_p_y * pp_p_y;
953 
954  int N = pn_p_x * pn_p_y;
955  int N2 = pn_e_x * pn_e_y;
956  Complex<Real>* hh_p = new Complex<Real>[N];
957  memset(hh_p, 0.0, sizeof(Complex<Real>) * N);
958 
959  int cropidx = 0;
960 
961  for (int p = 0; p < N2; p++)
962  {
963  int x = p % pn_e_x;
964  int y = p / pn_e_x;
965 
966  if (x >= eye_box_range_idx[0] - 1 && x <= eye_box_range_idx[1] - 1 && y >= eye_box_range_idy[0] - 1 && y <= eye_box_range_idy[1] - 1)
967  {
968  int xx = cropidx % pn_p_x;
969  int yy = cropidx / pn_p_x;
970  hh_p[yy*pn_p_x + xx] = hh_e_shift[p];
971  cropidx++;
972  }
973  }
974  delete[] hh_e_shift;
975 
976  Real f_eye = eyeLen * (eyeCenter[_Z] - eyeFocusDistance) / (eyeLen + (eyeCenter[_Z] - eyeFocusDistance));
977  Real effective_f = f_eye * eyeLen / (f_eye - eyeLen);
978 
979  Complex<Real>* hh_e_ = new Complex<Real>[N];
980  memset(hh_e_, 0.0, sizeof(Complex<Real>) * N);
981 
982  int loopp;
983  //#pragma omp parallel for private(loopp)
984  for (loopp = 0; loopp < N; loopp++)
985  {
986  int x = loopp % pn_p_x;
987  int y = loopp / pn_p_x;
988 
989  Real XE = -ss_p_x / 2.0 + (pp_p_x *x);
990  Real YE = ss_p_y / 2.0 - pp_p_y - (pp_p_y * y);
991 
992  Real sval = (XE*XE) + (YE*YE);
993  sval *= M_PI / lambda / effective_f;
994  Complex<Real> eye_propagation_kernel(0, sval);
995  eye_propagation_kernel.exp();
996 
997  Real eye_lens_anti_aliasing_mask = fabs(XE) < fabs(lambda*effective_f / (4 * pp_e_x)) ? 1.0 : 0.0;
998  eye_lens_anti_aliasing_mask *= fabs(YE) < fabs(lambda*effective_f / (4 * pp_e_y)) ? 1.0 : 0.0;
999 
1000  Real eye_pupil_mask = sqrt(XE*XE + YE * YE) < (eyePupil / 2.0) ? 1.0 : 0.0;
1001 
1002  hh_e_[x + y * pn_p_x] = hh_p[x + y * pn_p_x] * eye_propagation_kernel * eye_lens_anti_aliasing_mask * eye_pupil_mask;
1003  }
1004 
1005  delete[] hh_p;
1006 
1007  fft2(hh_e_, hh_e_, pn_p_x, pn_p_y, FFTW_FORWARD, false);
1008 
1009  Real pp_ret_x, pp_ret_y;
1010  int pn_ret_x, pn_ret_y;
1011  vec2 ret_size_xy;
1012 
1013  pp_ret_x = lambda * eyeLen / ss_p_x;
1014  pp_ret_y = lambda * eyeLen / ss_p_y;
1015  pn_ret_x = pn_p_x;
1016  pn_ret_y = pn_p_y;
1017  ret_size_xy[0] = pp_ret_x * pn_ret_x;
1018  ret_size_xy[1] = pp_ret_y * pn_ret_y;
1019 
1020  field_ret_set_[ctr] = new Real[pn_p_x * pn_p_y];
1021  memset(field_ret_set_[ctr], 0.0, sizeof(Real)*pn_p_x*pn_p_y);
1022 
1023  //#pragma omp parallel for private(loopp)
1024  for (loopp = 0; loopp < pn_p_x*pn_p_y; loopp++)
1025  {
1026  int x = loopp % pn_p_x;
1027  int y = loopp / pn_p_x;
1028 
1029  Real XR = ret_size_xy[0] / 2.0 + (pp_ret_x * x);
1030  Real YR = ret_size_xy[1] / 2.0 - pp_ret_y - (pp_ret_y * y);
1031 
1032  Real sval = (XR*XR) + (YR*YR);
1033  sval *= k / (2 * eyeLen);
1034  Complex<Real> val1(0, sval);
1035  val1.exp();
1036 
1037  Complex<Real> val2(0, k * eyeLen);
1038  val2.exp();
1039  Complex<Real> val3(0, lambda * eyeLen);
1040 
1041  field_ret_set_[ctr][x + pn_p_x * y] = (hh_e_[x + pn_p_x * y] * (val1 * val2 / val3)).mag();
1042 
1043  }
1044  delete[] hh_e_;
1045 
1046  pp_ret_set_[ctr] = vec2(pp_ret_x, pp_ret_y);
1047  pn_ret_set_[ctr] = ivec2(pn_ret_x, pn_ret_y);
1048  ss_ret_set_[ctr] = pp_ret_set_[ctr] * pn_ret_set_[ctr];
1049 
1050  if (bCreatePupilImg)
1051  {
1052  Real pp_min = (pp_e_x > pp_e_y) ? pp_e_y : pp_e_x;
1053  Real ss_max = (ss_e_x > ss_e_y) ? ss_e_x : ss_e_y;
1054  Real pn_tar = ceil(ss_max / pp_min);
1055  pp_min = ss_max / pn_tar;
1056  Real pn_e_tar_x = round(ss_e_x / pp_min);
1057  Real pn_e_tar_y = round(ss_e_y / pp_min);
1058 
1059  Real resize_scale_x = pn_e_tar_x * resultSizeScale;
1060  Real resize_scale_y = pn_e_tar_y * resultSizeScale;
1061 
1062  int N = resize_scale_x * resize_scale_y;
1063 
1064  recon_set[vtr * nChannel + ctr] = new Real[N];
1065  img_set[vtr * nChannel + ctr] = new uchar[N];
1066 
1067  memset(recon_set[vtr * nChannel + ctr], 0.0, sizeof(Real) * N);
1068  memset(img_set[vtr * nChannel + ctr], 0, sizeof(uchar) * N);
1069 
1070  GetPupilFieldImage(field_set_[ctr], recon_set[vtr * nChannel + ctr]
1071  , pn_e_x, pn_e_y, pp_e_x, pp_e_y, resize_scale_x, resize_scale_y);
1072 
1073  std::string fname = std::string("./").append("").append("/FIELD/");
1074  fname = fname.append("FIELD_COLOR_").append(std::to_string(ctr + 1)).append("_").append("").append("_SAT_").append("").append("_").append(std::to_string(vtr + 1)).append(".bmp");
1075  normalize(recon_set[vtr * nChannel + ctr], img_set[vtr * nChannel + ctr], (int)resize_scale_x, (int)resize_scale_y);
1076  img_size[vtr * nChannel + ctr][_X] = (int)resize_scale_x;
1077  img_size[vtr * nChannel + ctr][_Y] = (int)resize_scale_y;
1078  }
1079 
1080  } // end of ctr
1081 
1082  Real pnx_max = 0.0, pny_max = 0.0;
1083  for (int i = 0; i < nChannel; i++)
1084  {
1085  pnx_max = (pn_ret_set_[i][0] > pnx_max ? pn_ret_set_[i][0] : pnx_max);
1086  pny_max = (pn_ret_set_[i][1] > pny_max ? pn_ret_set_[i][1] : pny_max);
1087  }
1088 
1089  Real retinal_image_shift_x = eyeCenter[_X] * eyeLen / eyeCenter[_Z];
1090  Real retinal_image_shift_y = eyeCenter[_Y] * eyeLen / eyeCenter[_Z];
1091 
1092  res_set_.resize(nChannel);
1093  res_set_norm_255_.resize(nChannel);
1094 
1095  int loopi;
1096  //#pragma omp parallel for private(loopi)
1097  for (loopi = 0; loopi < nChannel; loopi++)
1098  {
1099  Real lambda = context_.wave_length[loopi];
1100  Real* hh_ret_ = new Real[pn_ret_set_[loopi][0] * pn_ret_set_[loopi][1]];
1101  memset(hh_ret_, 0.0, sizeof(Real)*pn_ret_set_[loopi][0] * pn_ret_set_[loopi][1]);
1102 
1103  if (bCenteringRetinaImg)
1104  {
1105  Real retinal_image_shift_by_pn_x = round(retinal_image_shift_x / pp_ret_set_[loopi][0]);
1106  Real retinal_image_shift_by_pn_y = round(retinal_image_shift_y / pp_ret_set_[loopi][1]);
1107 
1108  circshift(field_ret_set_[loopi], hh_ret_, -retinal_image_shift_by_pn_x, retinal_image_shift_by_pn_y, pn_ret_set_[loopi][0], pn_ret_set_[loopi][1]);
1109 
1110  }
1111  else
1112  hh_ret_ = field_ret_set_[loopi];
1113 
1114  delete[] field_ret_set_[loopi];
1115 
1116  int size = (int)(pnx_max * pny_max);
1117  res_set_[loopi] = new Real[size];
1118  memset(res_set_[loopi], 0.0, sizeof(Real) * size);
1119  ScaleBilnear(hh_ret_, res_set_[loopi], pn_ret_set_[loopi][0], pn_ret_set_[loopi][1], pnx_max, pny_max, lambda * lambda);
1120 
1121  }
1122 
1123  Real maxvalue = res_set_[0][0];
1124  Real minvalue = res_set_[0][0];
1125  for (int i = 0; i < nChannel; i++)
1126  {
1127  for (int j = 0; j < pnx_max*pny_max; j++)
1128  {
1129  maxvalue = res_set_[i][j] > maxvalue ? res_set_[i][j] : maxvalue;
1130  minvalue = res_set_[i][j] < minvalue ? res_set_[i][j] : minvalue;
1131  }
1132  }
1133 
1134  for (int j = 0; j < pnx_max*pny_max; j++)
1135  {
1136  res_set_[0][j] = (res_set_[0][j] - minvalue) / (maxvalue - minvalue) * 255;
1137  }
1138 
1139  //#pragma omp parallel for private(loopi)
1140  for (loopi = 0; loopi < nChannel; loopi++)
1141  {
1142  for (int j = 0; j < pnx_max*pny_max; j++)
1143  {
1144  if (res_set_[loopi][j] > ratioAtRetina * maxvalue)
1145  res_set_[loopi][j] = ratioAtRetina * maxvalue;
1146 
1147  res_set_[loopi][j] = (res_set_[loopi][j] - minvalue) / (ratioAtRetina*maxvalue - minvalue);
1148 
1149  }
1150  }
1151 
1152  Real ret_size_x = pnx_max * resultSizeScale;
1153  Real ret_size_y = pny_max * resultSizeScale;
1154 
1155  int size = (int)(ret_size_x * ret_size_y);
1156 
1157  //#pragma omp parallel for private(loopi)
1158  for (loopi = 0; loopi < nChannel; loopi++)
1159  {
1160  Real* res_set_norm = new Real[size];
1161  memset(res_set_norm, 0.0, sizeof(Real) * size);
1162 
1163  ScaleBilnear(res_set_[loopi], res_set_norm, pnx_max, pny_max, ret_size_x, ret_size_y);
1164 
1165  res_set_norm_255_[loopi] = new Real[size];
1166  memset(res_set_norm_255_[loopi], 0.0, sizeof(Real) * size);
1167 
1168  rotateCCW180(res_set_norm, res_set_norm_255_[loopi], ret_size_x, ret_size_y, 255.0);
1169 
1170  delete[] res_set_norm;
1171  }
1172 
1173  for (int i = 0; i < nChannel; i++)
1174  delete[] res_set_[i];
1175 
1176  int N = ret_size_x * ret_size_y;
1177 
1178  focus_recon_set[vtr] = new Real[N * nChannel];
1179  focus_img_set[vtr] = new uchar[N * nChannel];
1180  memset(focus_recon_set[vtr], 0.0, sizeof(Real) * N * nChannel);
1181  memset(focus_img_set[vtr], 0, sizeof(uchar) * N * nChannel);
1182 
1183  if (nChannel == 1)
1184  memcpy(focus_recon_set[vtr], res_set_norm_255_[0], sizeof(Real)*ret_size_x*ret_size_y);
1185  else
1186  reArrangeChannel(res_set_norm_255_, focus_recon_set[vtr], ret_size_x, ret_size_y, nChannel);
1187 
1188  std::string fname = std::string("./").append("").append("/").append("").append("_SAT_").append("").append("_").append(varname2).append("_").append(std::to_string(vtr + 1)).append(".bmp");
1189 
1190  normalize(focus_recon_set[vtr], focus_img_set[vtr], ret_size_x, ret_size_y);
1191 
1192  focus_img_size[vtr][_X] = (int)ret_size_x;
1193  focus_img_size[vtr][_Y] = (int)ret_size_y;
1194 
1195  m_vecEncoded.push_back(focus_recon_set[vtr]);
1196  m_vecNormalized.push_back(focus_img_set[vtr]);
1197 
1198 
1199  for (int i = 0; i < nChannel; i++)
1200  delete[] res_set_norm_255_[i];
1201  } // end of vtr
1202 }
1203 
1204 bool ophRec::save(const char * fname, uint8_t bitsperpixel, uchar* src, uint px, uint py)
1205 {
1206  bool bOK = false;
1207 
1208  if (fname == nullptr) return bOK;
1209 
1210  uchar* source = src;
1211  bool bAlloc = false;
1212  const uint nChannel = context_.waveNum;
1213 
1214  ivec2 p(px, py);
1215  if (px == 0 && py == 0)
1217 
1218 
1219  std::string file = fname;
1220  std::replace(file.begin(), file.end(), '\\', '/');
1221 
1222  // split path
1223  std::vector<std::string> components;
1224  std::stringstream ss(file);
1225  std::string item;
1226  char token = '/';
1227 
1228  while (std::getline(ss, item, token)) {
1229  components.push_back(item);
1230  }
1231 
1232  std::string dir;
1233 
1234  for (size_t i = 0; i < components.size() - 1; i++)
1235  {
1236  dir += components[i];
1237  dir += "/";
1238  }
1239 
1240  std::string filename = components[components.size() - 1];
1241 
1242  // find extension
1243  bool hasExt;
1244  size_t ext_pos = file.rfind(".");
1245  hasExt = (ext_pos == string::npos) ? false : true;
1246 
1247  if (!hasExt)
1248  filename.append(".bmp");
1249 
1250  std::string fullpath = dir + filename;
1251 
1252  if (nChannel == 1) {
1253  saveAsImg(fullpath.c_str(), bitsperpixel, src, p[_X], p[_Y]);
1254  }
1255  else if (nChannel == 3) {
1256  if (imgCfg.merge) {
1257  uint nSize = (((p[_X] * bitsperpixel >> 3) + 3) & ~3) * p[_Y];
1258  uchar *source = new uchar[nSize];
1259  bAlloc = true;
1260  for (uint i = 0; i < nChannel; i++) {
1261  mergeColor(i, p[_X], p[_Y], src, source);
1262  }
1263  saveAsImg(fullpath.c_str(), bitsperpixel, source, p[_X], p[_Y]);
1264  if (bAlloc) delete[] source;
1265  }
1266  else {
1267  char path[FILENAME_MAX] = { 0, };
1268  sprintf(path, "%s%s", dir.c_str(), filename.c_str());
1269  saveAsImg(path, bitsperpixel / nChannel, src, p[_X], p[_Y]);
1270  }
1271  }
1272  else return false;
1273 
1274  return bOK;
1275 }
1276 
1277 void ophRec::SaveImage(const char *path, const char *ext)
1278 {
1279  int nSimStep = rec_config.SimulationStep;
1280  int nChannel = context_.waveNum;
1281  bool bCreatePupilImg = rec_config.CreatePupilFieldImg;
1282  Real simTo = rec_config.SimulationTo;
1283  Real simFrom = rec_config.SimulationFrom;
1284  Real step = 0.0;
1285 
1286  int pnX;
1287  int pnY;
1288 
1289  if (m_idx == 0)
1290  {
1291  pnX = context_.pixel_number[_X];
1292  pnY = context_.pixel_number[_Y];
1293  }
1294  else
1295  {
1296  pnX = focus_img_size[0][_X];
1297  pnY = focus_img_size[0][_Y];
1298  }
1299  if (nSimStep > 1)
1300  {
1301  step = (simTo - simFrom) / (nSimStep - 1);
1302  }
1303 
1304  char tmpPath[FILENAME_MAX] = { 0, };
1305  bool bMultiStep = nSimStep > 1 ? true : false;
1306  uint nSize = (((pnX * nChannel) + 3) & ~3) * pnY;
1307  uchar *tmp = new uchar[nSize];
1308 
1309  for (int i = 0; i < nSimStep; i++)
1310  {
1311  sprintf(tmpPath, "%s\\FOCUS_%.4f.%s", path, bMultiStep ?
1312  simFrom + (step * i) : (simFrom + simTo) / 2, ext);
1313 
1314  if (nChannel == 3)
1315  {
1316  if (imgCfg.merge)
1317  {
1318  for (int ch = 0; ch < nChannel; ch++)
1319  {
1320  mergeColor(ch, pnX, pnY, m_vecNormalized[i * nChannel + ch], tmp);
1321  }
1322  saveAsImg(tmpPath, 8 * nChannel, tmp, pnX, pnY); // save RGB
1323  }
1324  else
1325  {
1326  for (int ch = 0; ch < nChannel; ch++)
1327  {
1328  memset(tmp, 0, sizeof(uchar) * nSize);
1329  sprintf(tmpPath, "%s\\FOCUS_%.4f (%d).%s", path, bMultiStep ? simFrom + step * i : (simFrom + simTo) / 2, ch, ext);
1330  mergeColor(ch, pnX, pnY, m_vecNormalized[i * nChannel + ch], tmp);
1331  saveAsImg(tmpPath, 8 * nChannel, tmp, pnX, pnY); // save RGB
1332  }
1333  }
1334  }
1335  else
1336  {
1337  memcpy(tmp, m_vecNormalized[i], sizeof(uchar) * nSize);
1338  saveAsImg(tmpPath, 8 * nChannel, tmp, pnX, pnY); // save Grayscale
1339  }
1340  }
1341  delete[] tmp;
1342 }
1343 
1344 void ophRec::GetPupilFieldImage(Complex<Real>* src, Real* dst, int pnx, int pny, Real ppx, Real ppy, Real scaleX, Real scaleY)
1345 {
1346  const int N = pnx * pny;
1347  Real ratioAtPupil = rec_config.RatioAtPupil;
1348  Real eyePupil = rec_config.EyePupilDiaMeter;
1349  vec3 eyeCenter = rec_config.EyeCenter;
1350 
1351 
1352  Real* dimg = new Real[N];
1353  memset(dimg, 0.0, sizeof(Real) * N);
1354 
1355  Real maxvalue = src[0].mag();
1356  for (int k = 0; k < N; k++)
1357  {
1358  Real val = src[k].mag();
1359  maxvalue = val > maxvalue ? val : maxvalue;
1360  dimg[k] = val;
1361  }
1362 
1363  Real SAT_VAL = maxvalue * ratioAtPupil;
1364  Real hss_x = (pnx * ppx) / 2.0;
1365  Real hss_y = (pny * ppy) / 2.0;
1366  Real halfPupil = eyePupil / 2.0;
1367 
1368  Real maxv = dimg[0], minv = dimg[0];
1369  for (int k = 0; k < N; k++)
1370  {
1371  if (dimg[k] > SAT_VAL)
1372  dimg[k] = SAT_VAL;
1373 
1374  maxv = (dimg[k] > maxv) ? dimg[k] : maxv;
1375  minv = (dimg[k] < minv) ? dimg[k] : minv;
1376  }
1377 
1378  for (int k = 0; k < N; k++)
1379  {
1380  int x = k % pnx;
1381  int y = k / pnx;
1382 
1383  Real xx = -hss_x + (ppx * x);
1384  Real yy = -hss_y + (pny - 1) * ppy - (ppy * y);
1385  Real xxx = xx - eyeCenter[_X];
1386  Real yyy = yy - eyeCenter[_Y];
1387  bool eye_pupil_mask = sqrt(xxx * xxx + yyy * yyy) < halfPupil ? 1.0 : 0.0;
1388 
1389  Real val = dimg[k];
1390  Real field_out = (val - minv) / (maxv - minv) * 255.0;
1391  Real field_in = (val - minv) / (maxv - minv) * 127.0 + 128.0;
1392  dimg[k] = field_in * eye_pupil_mask + field_out * (1 - eye_pupil_mask);
1393  //dimg[k] = field_in * eye_pupil_mask;
1394 
1395  }
1396  ScaleBilnear(dimg, dst, pnx, pny, scaleX, scaleY);
1397 
1398  delete[] dimg;
1399 }
1400 
1401 void ophRec::getVarname(int vtr, vec3& var_vals, std::string& varname2)
1402 {
1403  bool bSimPos[3];
1404  memcpy(&bSimPos[0], rec_config.SimulationPos, sizeof(bSimPos));
1405 
1406  std::string varname;
1407  varname.clear();
1408  varname2.clear();
1409  if (rec_config.SimulationMode == 0) {
1410 
1411  varname.append("FOCUS");
1412 
1413  LOG("Step # %d %s = %.8f \n", vtr + 1, varname.c_str(), var_vals[0]);
1414  char temp[100];
1415  sprintf(temp, "%.5f", var_vals[0]);
1416  varname2.clear();
1417  varname2.append(varname).append("_").append(temp);
1418  varname2.replace(varname2.find('.'), 1, "_");
1419 
1420  }
1421  else {
1422 
1423  varname.append("POSITION");
1424 
1425  if (bSimPos[_X]) varname.append("_X");
1426  if (bSimPos[_Y]) varname.append("_Y");
1427  if (bSimPos[_Z]) varname.append("_Z");
1428 
1429  LOG("Step # %d %s = ", vtr + 1, varname.c_str());
1430  varname2.clear();
1431  varname2.append(varname);
1432 
1433  if (bSimPos[_X]) {
1434  LOG("%.8f ", var_vals[0]);
1435  char temp[100];
1436  sprintf(temp, "%.5f", var_vals[0]);
1437  varname2.append("_").append(temp);
1438  varname2.replace(varname2.find('.'), 1, "_");
1439  }
1440  if (bSimPos[_Y]) {
1441  LOG("%.8f ", var_vals[1]);
1442  char temp[100];
1443  sprintf(temp, "%.5f", var_vals[1]);
1444  varname2.append("_").append(temp);
1445  varname2.replace(varname2.find('.'), 1, "_");
1446  }
1447  if (bSimPos[_Z]) {
1448  LOG("%.8f ", var_vals[2]);
1449  char temp[100];
1450  sprintf(temp, "%.5f", var_vals[2]);
1451  varname2.append("_").append(temp);
1452  varname2.replace(varname2.find('.'), 1, "_");
1453  }
1454  LOG("\n");
1455  }
1456 }
1457 
1459 {
1460  for (vector<Real*>::iterator it = m_vecEncoded.begin(); it != m_vecEncoded.end(); it++)
1461  {
1462  delete[](*it);
1463  }
1464 
1465  for (vector<uchar*>::iterator it = m_vecNormalized.begin(); it != m_vecNormalized.end(); it++)
1466  {
1467  delete[](*it);
1468  }
1469  m_vecEncoded.clear();
1470  m_vecNormalized.clear();
1471 }
1472 
1474 {
1475 
1476  Clear();
1477  auto begin = CUR_TIME;
1478 
1479  LOG("1) Algorithm Method : Angular Spectrum\n");
1480  LOG("2) Reconstruct Image with %s\n", m_mode & MODE_GPU ?
1481  "GPU" :
1482 #ifdef _OPENMP
1483  "Multi Core CPU"
1484 #else
1485  "Single Core CPU"
1486 #endif
1487  );
1488 
1489  m_mode & MODE_GPU ? ASM_Propagation_GPU() : ASM_Propagation();
1490 
1491  LOG("Total Elapsed Time: %lf (s)\n", ELAPSED_TIME(begin, CUR_TIME));
1492  return true;
1493 }
1494 
1496 {
1497  const int nChannel = context_.waveNum;
1498  const int N = context_.pixel_number[_X] * context_.pixel_number[_Y];
1499 
1500  // Memory Location for Result Image
1501  if (complex_H != nullptr) {
1502  for (int i = 0; i < m_nOldChannel; i++) {
1503  if (complex_H[i] != nullptr) {
1504  delete[] complex_H[i];
1505  complex_H[i] = nullptr;
1506  }
1507  }
1508  delete[] complex_H;
1509  complex_H = nullptr;
1510  }
1511  complex_H = new Complex<Real>*[nChannel];
1512  for (int i = 0; i < nChannel; i++) {
1513  complex_H[i] = new Complex<Real>[N];
1514  memset(complex_H[i], 0, sizeof(Complex<Real>) * N);
1515  }
1516  m_nOldChannel = nChannel;
1517 }
1518 
1520 {
1522 
1523  for (vector<Real*>::iterator it = m_vecEncoded.begin(); it != m_vecEncoded.end(); it++)
1524  {
1525  delete (*it);
1526  }
1527  m_vecEncoded.clear();
1528 
1529  for (vector<uchar*>::iterator it = m_vecNormalized.begin(); it != m_vecNormalized.end(); it++)
1530  {
1531  delete (*it);
1532  }
1533  m_vecNormalized.clear();
1534 }
1535 
1536 void circshift(Real* in, Real* out, int shift_x, int shift_y, int nx, int ny)
1537 {
1538  int ti, tj;
1539  for (int i = 0; i < nx; i++)
1540  {
1541  for (int j = 0; j < ny; j++)
1542  {
1543  ti = (i + shift_x) % nx;
1544  if (ti < 0)
1545  ti = ti + nx;
1546  tj = (j + shift_y) % ny;
1547  if (tj < 0)
1548  tj = tj + ny;
1549 
1550  out[ti + tj * nx] = in[i + j * nx];
1551  }
1552  }
1553 }
1554 
1555 void circshift(Complex<Real>* in, Complex<Real>* out, int shift_x, int shift_y, int nx, int ny)
1556 {
1557  int ti, tj;
1558  for (int i = 0; i < nx; i++)
1559  {
1560  for (int j = 0; j < ny; j++)
1561  {
1562  ti = (i + shift_x) % nx;
1563  if (ti < 0)
1564  ti = ti + nx;
1565  tj = (j + shift_y) % ny;
1566  if (tj < 0)
1567  tj = tj + ny;
1568 
1569  out[ti + tj * nx] = in[i + j * nx];
1570  }
1571  }
1572 }
1573 
1574 vec3 image_sample(float xx, float yy, int c, size_t w, size_t h, double* in)
1575 {
1576  int x1 = (int)floor(xx);
1577  int x2 = (int)ceil(xx);
1578  int y1 = (int)floor(yy);
1579  int y2 = (int)ceil(yy);
1580 
1581  if (x1 < 0 || x1 >= (int)w || x2 < 0 || x2 >= (int)w) return vec3(0);
1582  if (y1 < 0 || y1 >= (int)h || y2 < 0 || y2 >= (int)h) return vec3(0);
1583 
1584  vec3 ret(0);
1585  double v1, v2, v3, v4, tvx, tvy;
1586 
1587  tvx = xx - floor(xx);
1588  tvy = yy - floor(yy);
1589  int tc = min(c, 3);
1590  for (int i = 0; i < tc; i++) {
1591  v1 = IMAGE_VAL(x1, y1, i, w, c, in);
1592  v2 = IMAGE_VAL(x2, y1, i, w, c, in);
1593  v3 = IMAGE_VAL(x1, y2, i, w, c, in);
1594  v4 = IMAGE_VAL(x2, y2, i, w, c, in);
1595  v1 = (1.0 - tvx)*v1 + tvx * v2;
1596  v3 = (1.0 - tvx)*v3 + tvx * v4;
1597  v1 = (1.0 - tvy)*v1 + tvy * v3;
1598 
1599  ret[i] = v1;
1600  }
1601 
1602  return ret;
1603 }
1604 
1605 void ScaleBilnear(double* src, double* dst, int w, int h, int neww, int newh, double multiplyval)
1606 {
1607  for (int y = 0; y < newh; y++)
1608  {
1609  for (int x = 0; x < neww; x++)
1610  {
1611  float gx = (x / (float)neww) * (w - 1);
1612  float gy = (y / (float)newh) * (h - 1);
1613 
1614  vec3 ret = multiplyval * image_sample(gx, gy, 1, w, h, src);
1615  dst[x + y * neww] = ret[0];
1616 
1617  /*
1618  int gxi = (int)gx;
1619  int gyi = (int)gy;
1620 
1621  double a00 = src[gxi + 0 + gyi * w];
1622  double a01 = src[gxi + 1 + gyi * w];
1623  double a10 = src[gxi + 0 + (gyi + 1)*w];
1624  double a11 = src[gxi + 1 + (gyi + 1)*w];
1625 
1626  float dx = gx - gxi;
1627  float dy = gy - gyi;
1628 
1629  dst[x + y * neww] = multiplyval * (a00 * (1 - dx)*(1 - dy) + a01 * dx*(1 - dy) + a10 * (1 - dx)*dy + a11 * dx*dy);
1630  */
1631 
1632  }
1633  }
1634 }
1635 
1636 
1637 
1638 void rotateCCW180(double* src, double* dst, int pnx, int pny, double mulival)
1639 {
1640  for (int i = 0; i < pnx*pny; i++)
1641  {
1642  int x = i % pnx;
1643  int y = i / pnx;
1644 
1645  int newx = pnx - x - 1;
1646  int newy = pny - y - 1;
1647 
1648  dst[newx + newy * pnx] = mulival * src[x + y * pnx];
1649 
1650  }
1651 
1652 }
1653 
1654 void reArrangeChannel(std::vector<double*>& src, double* dst, int pnx, int pny, int chnum)
1655 {
1656  for (int k = 0; k < pnx*pny; k++)
1657  {
1658  for (int c = 0; c < chnum; c++)
1659  {
1660  if (c == 0)
1661  dst[k*chnum + 2] = src[c][k];
1662  else if (c == 2)
1663  dst[k*chnum + 0] = src[c][k];
1664  else
1665  dst[k*chnum + c] = src[c][k];
1666  }
1667 
1668  }
1669 }
Real SimulationTo
Definition: ophRec.h:68
Real RatioAtRetina
Definition: ophRec.h:72
Abstract class.
Definition: Openholo.h:145
#define M_PI
Definition: define.h:52
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:135
void reArrangeChannel(std::vector< double *> &src, double *dst, int pnx, int pny, int chnum)
Definition: ophRec.cpp:1654
w
Definition: Openholo.cpp:429
int SimulationStep
Definition: ophRec.h:70
void fftFree(void)
Resource release method.
Definition: Openholo.cpp:677
Real SimulationFrom
Definition: ophRec.h:69
bytesperpixel
Definition: Openholo.cpp:431
void GetPupilFieldFromVWHologram()
Definition: ophRec.cpp:470
#define MIN_DOUBLE
Definition: define.h:148
void Clear()
Definition: ophRec.cpp:1458
unsigned char uchar
Definition: typedef.h:64
HORIZONTAL
Definition: ImgControl.h:25
void SaveImage(const char *path, const char *ext="bmp")
Definition: ophRec.cpp:1277
void rotateCCW180(double *src, double *dst, int pnx, int pny, double mulival=1.0)
Definition: ophRec.cpp:1638
#define MODE_GPU
Definition: define.h:156
void Initialize()
Definition: ophRec.cpp:1495
NONE
Definition: ImgControl.h:25
vec3 EyeCenter
Definition: ophRec.h:65
#define MAX_DOUBLE
Definition: define.h:140
vec2 ss
Definition: Openholo.h:104
return true
Definition: Openholo.cpp:434
float Real
Definition: typedef.h:55
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
void getVarname(int vtr, vec3 &var_vals, std::string &varname2)
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: ophRec.cpp:1519
#define _IM
Definition: complex.h:58
def k(wvl)
Definition: Depthmap.py:16
h
Definition: Openholo.cpp:430
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
bool readImage(const char *path)
Definition: ophRec.cpp:244
VERTICAL
Definition: ImgControl.h:25
virtual ~ophRec(void)
Destructor.
Definition: ophRec.cpp:63
#define FFTW_BACKWARD
Definition: fftw3.h:380
bool rotate
Definition: Openholo.h:124
vec2 EyeBoxSize
Definition: ophRec.h:63
ImageConfig imgCfg
Definition: Openholo.h:487
bool CreatePupilFieldImg
Definition: ophRec.h:74
vec2 pixel_pitch
Definition: Openholo.h:101
void GetPupilFieldFromHologram()
Definition: ophRec.cpp:449
Real RatioAtPupil
Definition: ophRec.h:73
void ASM_Propagation()
Definition: ophRec.cpp:592
#define _X
Definition: define.h:92
void circshift(Real *in, Real *out, int shift_x, int shift_y, int nx, int ny)
Definition: ophRec.cpp:1536
bool ReconstructImage()
Definition: ophRec.cpp:1473
Real EyeFocusDistance
Definition: ophRec.h:66
void imgScaleBilinear(uchar *src, uchar *dst, int w, int h, int neww, int newh, int channels=1)
Function for change image size.
Definition: Openholo.cpp:437
vec3 image_sample(float xx, float yy, int c, size_t w, size_t h, double *in)
Definition: ophRec.cpp:1574
bool getImgSize(int &w, int &h, int &bytesperpixel, const char *fname)
Function for getting the image size.
Definition: Openholo.cpp:402
#define _RE
Definition: complex.h:55
bool readImagePNA(const char *phase, const char *amplitude)
Definition: ophRec.cpp:291
Real EyePupilDiaMeter
Definition: ophRec.h:61
bool readImageRNI(const char *real, const char *imaginary)
Definition: ophRec.cpp:375
void fft2(ivec2 n, Complex< Real > *in, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Functions for performing fftw 2-dimension operations inside Openholo.
Definition: Openholo.cpp:559
#define FFTW_FORWARD
Definition: fftw3.h:379
ophRec(void)
Constructor.
Definition: ophRec.cpp:55
bool SimulationPos[3]
Definition: ophRec.h:77
#define IMAGE_VAL(x, y, c, w, n, mem)
Definition: ophRec.cpp:66
int SimulationMode
Definition: ophRec.h:71
uint waveNum
Definition: Openholo.h:105
void GetPupilFieldImage(Complex< Real > *src, double *dst, int pnx, int pny, double ppx, double ppy, double scaleX, double scaleY)
ivec2 pixel_number
Definition: Openholo.h:99
Real EyeLength
Definition: ophRec.h:60
bool save(const char *fname, uint8_t bitsperpixel, uchar *src, uint px, uint py)
#define _Y
Definition: define.h:96
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
bool CenteringRetinaImg
Definition: ophRec.h:75
void ASM_Propagation_GPU()
Definition: ophRec_GPU.cpp:49
bool readConfig(const char *fname)
Definition: ophRec.cpp:74
Complex< Real > ** complex_H
Definition: Openholo.h:489
void ScaleBilnear(double *src, double *dst, int w, int h, int neww, int newh, double multiplyval=1.0)
Definition: ophRec.cpp:1605
#define _Z
Definition: define.h:100
int EyeBoxUnit
Definition: ophRec.h:64
#define ELAPSED_TIME(x, y)
Definition: function.h:59
void Perform_Simulation()
Definition: ophRec.cpp:832
OphConfig context_
Definition: Openholo.h:485
void Propagation_Fresnel_FFT(int chnum)
Definition: ophRec.cpp:729
Real ResultSizeScale
Definition: ophRec.h:67
Real * wave_length
Definition: Openholo.h:106
bool loadAsImgUpSideDown(const char *fname, uchar *dst)
Function for loading image files | Output image data upside down.
Definition: Openholo.cpp:353
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
bool mergeColor(int idx, int width, int height, uchar *src, uchar *dst)
Function for generate RGB image from each grayscale image.
Definition: Openholo.cpp:103
unsigned int uint
Definition: typedef.h:62
#define CUR_TIME
Definition: function.h:58
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: Openholo.cpp:805
void normalize(T *src, uchar *dst, int x, int y)
bool merge
Definition: Openholo.h:125
#define FFTW_ESTIMATE
Definition: fftw3.h:392