Openholo  v4.2
Open Source Digital Holographic Library
ophCascadedPropagation.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 "ophCascadedPropagation.h"
47 #include "sys.h"
48 #include "tinyxml2.h"
49 #include <string>
50 #include <cfloat>
51 
52 
53 
54 ophCascadedPropagation::ophCascadedPropagation()
55  : ready_to_propagate(false),
56  hologram_path(L"")
57 {
58 }
59 
60 ophCascadedPropagation::ophCascadedPropagation(const wchar_t* configfilepath)
61  : ready_to_propagate(false),
62  hologram_path(L"")
63 {
64  if (readConfig(configfilepath) && allocateMem())
65  {
66  string hologram_path_str;
67  hologram_path_str.assign(hologram_path.begin(), hologram_path.end());
68  ready_to_propagate = (sourcetype_ == SourceType::IMG && loadInputImg(hologram_path_str)) || (sourcetype_ == SourceType::OHC && loadAsOhc(hologram_path_str.c_str()));
69  }
70 }
71 
73 {
74 }
75 
77 {
78  deallocateMem();
79 }
80 
82 {
83  if (!isReadyToPropagate())
84  {
85  PRINT_ERROR("module not initialized");
86  return false;
87  }
88 
89  if (!propagateSlmToPupil())
90  {
91  PRINT_ERROR("failed to propagate to pupil plane");
92  return false;
93  }
94 
96  {
97  PRINT_ERROR("failed to propagate to retina plane");
98  return false;
99  }
100 
101  return true;
102 }
103 
104 bool ophCascadedPropagation::save(const wchar_t* pathname, uint8_t bitsperpixel)
105 {
106  wstring bufw(pathname);
107  string bufs;
108  bufs.assign(bufw.begin(), bufw.end());
109  oph::uchar* src = getIntensityfields(getRetinaWavefieldAll());
110  return saveAsImg(bufs.c_str(), bitsperpixel, src, getResX(), getResY());
111 }
112 
113 bool ophCascadedPropagation::saveAsOhc(const char * fname)
114 {
115  oph::uint nx = getResX();
116  oph::uint ny = getResY();
117  for (oph::uint i = 0; i < getNumColors(); i++)
118  memcpy(complex_H[i], wavefield_retina[i], nx * ny * sizeof(Complex<Real>));
119 
120  if (!context_.wave_length)
122  for (oph::uint i = 0; i < getNumColors(); i++)
123  context_.wave_length[i] = config_.wavelengths[i];
124  context_.pixel_pitch[0] = config_.dx;
125  context_.pixel_pitch[1] = config_.dy;
126  context_.pixel_number[0] = config_.nx;
127  context_.pixel_number[1] = config_.ny;
128 
129  return Openholo::saveAsOhc(fname);
130 }
131 
132 bool ophCascadedPropagation::loadAsOhc(const char * fname)
133 {
134  if (!Openholo::loadAsOhc(fname))
135  return -1;
136 
137  oph::uint nx = getResX();
138  oph::uint ny = getResY();
139  config_.num_colors = OHC_decoder->getNumOfWavlen();
140  for (oph::uint i = 0; i < getNumColors(); i++)
141  memcpy(wavefield_SLM[i], complex_H[i], nx * ny * sizeof(Complex<Real>));
142 
143  return 1;
144 }
145 
146 bool ophCascadedPropagation::allocateMem()
147 {
148  wavefield_SLM.resize(getNumColors());
149  wavefield_pupil.resize(getNumColors());
150  wavefield_retina.resize(getNumColors());
151  oph::uint nx = getResX();
152  oph::uint ny = getResY();
153  complex_H = new Complex<Real>*[getNumColors()];
154  for (oph::uint i = 0; i < getNumColors(); i++)
155  {
156  wavefield_SLM[i] = new oph::Complex<Real>[nx * ny];
157  wavefield_pupil[i] = new oph::Complex<Real>[nx * ny];
158  wavefield_retina[i] = new oph::Complex<Real>[nx * ny];
159  complex_H[i] = new oph::Complex<Real>[nx * ny];
160  }
161 
162  return true;
163 }
164 
165 void ophCascadedPropagation::deallocateMem()
166 {
167  for (auto e : wavefield_SLM)
168  delete[] e;
169  wavefield_SLM.clear();
170 
171  for (auto e : wavefield_pupil)
172  delete[] e;
173  wavefield_pupil.clear();
174 
175  for (auto e : wavefield_retina)
176  delete[] e;
177  wavefield_retina.clear();
178 
179  for (oph::uint i = 0; i < getNumColors(); i++)
180  delete[] complex_H[i];
181 }
182 
183 // read in hologram data
184 bool ophCascadedPropagation::loadInputImg(string hologram_path_str)
185 {
186  if (!checkExtension(hologram_path_str.c_str(), ".bmp"))
187  {
188  PRINT_ERROR("input file format not supported");
189  return false;
190  }
191  oph::uint nx = getResX();
192  oph::uint ny = getResY();
193  oph::uchar* data = new oph::uchar[nx * ny * getNumColors()];
194  if (!loadAsImgUpSideDown(hologram_path_str.c_str(), data)) // loadAsImg() keeps to fail
195  {
196  PRINT_ERROR("input file not found");
197  delete[] data;
198  return false;
199  }
200 
201  // copy data to wavefield
202  try {
203  oph::uint numColors = getNumColors();
204  for (oph::uint row = 0; row < ny; row++)
205  {
206  for (oph::uint col = 0; col < nx; col++)
207  {
208  for (oph::uint color = 0; color < numColors; color++)
209  {
210  // BGR to RGB & upside-down
211  wavefield_SLM[numColors - 1 - color][(ny - 1 - row) * nx+ col] = oph::Complex<Real>((Real)data[(row * nx + col) * numColors + color], 0);
212  }
213  }
214  }
215  }
216 
217  catch (...) {
218  PRINT_ERROR("failed to generate wavefield from bmp");
219  delete[] data;
220  return false;
221  }
222 
223  delete[] data;
224  return true;
225 }
226 
227 oph::uchar* ophCascadedPropagation::getIntensityfields(vector<oph::Complex<Real>*> waveFields)
228 {
229  oph::uint numColors = getNumColors();
230  if (numColors != 1 && numColors != 3)
231  {
232  PRINT_ERROR("invalid number of color channels");
233  return nullptr;
234  }
235 
236  oph::uint nx = getResX();
237  oph::uint ny = getResY();
238  oph::uchar* intensityField = new oph::uchar[nx * ny * numColors];
239  for (oph::uint color = 0; color < numColors; color++)
240  {
241  Real* intensityFieldUnnormalized = new Real[nx * ny];
242 
243  // find minmax
244  Real maxIntensity = 0.0;
245  Real minIntensity = REAL_IS_DOUBLE ? DBL_MAX : FLT_MAX;
246  for (oph::uint row = 0; row < ny; row++)
247  {
248  for (oph::uint col = 0; col < nx; col++)
249  {
250  intensityFieldUnnormalized[row * nx + col] = waveFields[color][row * nx + col].mag2();
251  maxIntensity = max(maxIntensity, intensityFieldUnnormalized[row * nx + col]);
252  minIntensity = min(minIntensity, intensityFieldUnnormalized[row * nx + col]);
253  }
254  }
255 
256  maxIntensity *= getNor(); // IS IT REALLY NEEDED?
257  if (maxIntensity <= minIntensity)
258  {
259  for (oph::uint row = 0; row < ny; row++)
260  {
261  for (oph::uint col = 0; col < nx; col++)
262  {
263  intensityField[(row * nx + col) * numColors + (numColors - 1 - color)] = 0; // flip RGB order
264  }
265  }
266  }
267  else
268  {
269  for (oph::uint row = 0; row < ny; row++)
270  {
271  for (oph::uint col = 0; col < nx; col++)
272  {
273  Real normalizedVal = (intensityFieldUnnormalized[row * nx + col] - minIntensity) / (maxIntensity - minIntensity);
274  normalizedVal = min(1.0, normalizedVal);
275 
276  // rotate 180 & RGB flip
277  intensityField[((ny - 1 - row) * nx + (nx - 1 - col)) * numColors + (numColors - 1 - color)] = (oph::uchar)(normalizedVal * 255);
278  }
279  }
280  }
281  delete[] intensityFieldUnnormalized;
282  }
283 
284  return intensityField;
285 }
286 
287 bool ophCascadedPropagation::readConfig(const wchar_t* fname)
288 {
289  /*XML parsing*/
290  tinyxml2::XMLDocument xml_doc;
291  tinyxml2::XMLNode *xml_node;
292  wstring fnamew(fname);
293  string fnames;
294  fnames.assign(fnamew.begin(), fnamew.end());
295 
296  if (!checkExtension(fnames.c_str(), ".xml"))
297  {
298  LOG("file's extension is not 'xml'\n");
299  return false;
300  }
301  auto ret = xml_doc.LoadFile(fnames.c_str());
302  if (ret != tinyxml2::XML_SUCCESS)
303  {
304  LOG("Failed to load file \"%s\"\n", fnames.c_str());
305  return false;
306  }
307 
308  xml_node = xml_doc.FirstChild();
309  auto next = xml_node->FirstChildElement("SourceType");
310  if (!next || !(next->GetText()))
311  return false;
312  string sourceTypeStr = (xml_node->FirstChildElement("SourceType"))->GetText();
313  if (sourceTypeStr == string("IMG"))
314  sourcetype_ = SourceType::IMG;
315  else if (sourceTypeStr == string("OHC"))
316  sourcetype_ = SourceType::OHC;
317 
318  next = xml_node->FirstChildElement("NumColors");
319  if (!next || tinyxml2::XML_SUCCESS != next->QueryUnsignedText(&config_.num_colors))
320  return false;
321  if (getNumColors() == 0 || getNumColors() > 3)
322  return false;
323 
324  if (config_.num_colors >= 1)
325  {
326  next = xml_node->FirstChildElement("WavelengthR");
327  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.wavelengths[0]))
328  return false;
329  }
330  if (config_.num_colors >= 2)
331  {
332  next = xml_node->FirstChildElement("WavelengthG");
333  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.wavelengths[1]))
334  return false;
335  }
336  if (config_.num_colors == 3)
337  {
338  next = xml_node->FirstChildElement("WavelengthB");
339  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.wavelengths[2]))
340  return false;
341  }
342 
343  next = xml_node->FirstChildElement("PixelPitchHor");
344  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.dx))
345  return false;
346  next = xml_node->FirstChildElement("PixelPitchVer");
347  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.dy))
348  return false;
349  if (getPixelPitchX() != getPixelPitchY())
350  {
351  PRINT_ERROR("current implementation assumes pixel pitches are same for X and Y axes");
352  return false;
353  }
354 
355  next = xml_node->FirstChildElement("ResolutionHor");
356  if (!next || tinyxml2::XML_SUCCESS != next->QueryUnsignedText(&config_.nx))
357  return false;
358  next = xml_node->FirstChildElement("ResolutionVer");
359  if (!next || tinyxml2::XML_SUCCESS != next->QueryUnsignedText(&config_.ny))
360  return false;
361 
362  if (!context_.wave_length)
364  for (oph::uint i = 0; i < getNumColors(); i++)
365  context_.wave_length[i] = config_.wavelengths[i];
366  context_.pixel_pitch[0] = config_.dx;
367  context_.pixel_pitch[1] = config_.dy;
368  context_.pixel_number[0] = config_.nx;
369  context_.pixel_number[1] = config_.ny;
370 
373  for (oph::uint i = 0; i < getNumColors(); i++)
375 
376  next = xml_node->FirstChildElement("FieldLensFocalLength");
377  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.field_lens_focal_length))
378  return false;
379  next = xml_node->FirstChildElement("DistReconstructionPlaneToPupil");
380  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.dist_reconstruction_plane_to_pupil))
381  return false;
382  next = xml_node->FirstChildElement("DistPupilToRetina");
383  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.dist_pupil_to_retina))
384  return false;
385  next = xml_node->FirstChildElement("PupilDiameter");
386  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.pupil_diameter))
387  return false;
388  next = xml_node->FirstChildElement("Nor");
389  if (!next || tinyxml2::XML_SUCCESS != next->QueryDoubleText(&config_.nor))
390  return false;
391 
392  next = xml_node->FirstChildElement("HologramPath");
393  if (!next || !(next->GetText()))
394  return false;
395  string holopaths = (xml_node->FirstChildElement("HologramPath"))->GetText();
396  hologram_path.assign(holopaths.begin(), holopaths.end());
397 
398  return true;
399 }
400 
402 {
403  auto start_time = CUR_TIME;
404  oph::uint numColors = getNumColors();
405  oph::uint nx = getResX();
406  oph::uint ny = getResY();
407  oph::Complex<Real>* buf = new oph::Complex<Real>[nx * ny];
408  for (oph::uint color = 0; color < numColors; color++)
409  {
410  fft2(getSlmWavefield(color), buf, nx, ny, OPH_FORWARD, false);
411 
412  Real k = 2 * M_PI / getWavelengths()[color];
414  Real dx1 = vw / (Real)nx;
415  Real dy1 = vw / (Real)ny;
416  for (oph::uint row = 0; row < ny; row++)
417  {
418  Real Y1 = ((Real)row - ((Real)ny - 1) * 0.5f) * dy1;
419  for (oph::uint col = 0; col < nx; col++)
420  {
421  Real X1 = ((Real)col - ((Real)nx - 1) * 0.5f) * dx1;
422 
423  // 1st propagation
424  oph::Complex<Real> t1 = oph::Complex<Real>(0, k / 2 / getFieldLensFocalLength() * (X1 * X1 + Y1 * Y1)).exp();
426  buf[row * nx + col] = t1 / t2 * buf[row * nx + col];
427 
428  // applying aperture: need some optimization later
429  if ((sqrt(X1 * X1 + Y1 * Y1) >= getPupilDiameter() / 2) || (row >= ny / 2 - 1))
430  buf[row * nx + col] = 0;
431 
433  oph::Complex<Real> t3 = oph::Complex<Real>(0, -k / 2 / f_eye * (X1 * X1 + Y1 * Y1)).exp();
434  buf[row * nx + col] *= t3;
435  }
436  }
437 
438  memcpy(getPupilWavefield(color), buf, sizeof(oph::Complex<Real>) * nx * ny);
439  }
440 
441  auto end_time = CUR_TIME;
442 
443  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
444 
445  LOG("SLM to Pupil propagation - Implement time : %.5lf sec\n", during_time);
446 
447  delete[] buf;
448  return true;
449 }
450 
452 {
453  auto start_time = CUR_TIME;
454  oph::uint numColors = getNumColors();
455  oph::uint nx = getResX();
456  oph::uint ny = getResY();
457  oph::Complex<Real>* buf = new oph::Complex<Real>[nx * ny];
458  for (oph::uint color = 0; color < numColors; color++)
459  {
460  memcpy(buf, getPupilWavefield(color), sizeof(oph::Complex<Real>) * nx * ny);
461 
462  Real k = 2 * M_PI / getWavelengths()[color];
464  Real dx1 = vw / (Real)nx;
465  Real dy1 = vw / (Real)ny;
466  for (oph::uint row = 0; row < ny; row++)
467  {
468  Real Y1 = ((Real)row - ((Real)ny - 1) * 0.5f) * dy1;
469  for (oph::uint col = 0; col < nx; col++)
470  {
471  Real X1 = ((Real)col - ((Real)nx - 1) * 0.5f) * dx1;
472 
473  // 2nd propagation
474  oph::Complex<Real> t1 = oph::Complex<Real>(0, k / 2 / getDistPupilToRetina() * (X1 * X1 + Y1 * Y1)).exp();
475  buf[row * nx + col] *= t1;
476  }
477  }
478 
479  fft2(buf, getRetinaWavefield(color), nx, ny, OPH_FORWARD, false);
480  }
481 
482  auto end_time = CUR_TIME;
483 
484  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
485 
486  LOG("Pupil to Retina propagation - Implement time : %.5lf sec\n", during_time);
487 
488  delete[] buf;
489  return true;
490 }
491 
493 {
494  if (wavefield_SLM.size() <= (size_t)id)
495  return nullptr;
496  return wavefield_SLM[id];
497 }
498 
500 {
501  if (wavefield_pupil.size() <= (size_t)id)
502  return nullptr;
503  return wavefield_pupil[id];
504 }
505 
507 {
508  if (wavefield_retina.size() <= (size_t)id)
509  return nullptr;
510  return wavefield_retina[id];
511 }
512 
513 vector<oph::Complex<Real>*> ophCascadedPropagation::getRetinaWavefieldAll()
514 {
515  return wavefield_retina;
516 }
517 
518 
#define M_PI
Definition: define.h:52
#define OPH_FORWARD
Definition: define.h:66
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:135
oph::uint getResX()
Returns horizontal resolution.
bool save(const wchar_t *pathname, uint8_t bitsperpixel)
Save wavefield at retina plane as Windows Bitmap file.
#define PRINT_ERROR(errorMsg)
Real getPupilDiameter()
Returns diameter of pupil in meter.
Real getDistPupilToRetina()
Returns distance from pupil plane to retina plane in meter.
void setPixelNumberOHC(const ivec2 pixel_number)
getter/setter for OHC file read and write
Definition: Openholo.h:501
#define REAL_IS_DOUBLE
Definition: typedef.h:49
unsigned char uchar
Definition: typedef.h:64
oph::Complex< Real > * getRetinaWavefield(oph::uint id)
Return monochromatic wavefield at retina plane.
Real getFieldLensFocalLength()
Returns focal length of field lens in meter.
bool propagateSlmToPupil()
Calculates 1st propagation (from SLM plane to pupil plane)
void setPixelPitchOHC(const vec2 pixel_pitch)
Definition: Openholo.h:504
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
float Real
Definition: typedef.h:55
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
def k(wvl)
Definition: Depthmap.py:16
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
vector< oph::Complex< Real > * > getRetinaWavefieldAll()
Return all wavefields at retina plane.
vec2 pixel_pitch
Definition: Openholo.h:101
virtual bool saveAsOhc(const char *fname)
Function to write OHC file.
oph::uint getResY()
Returns vertical resolution.
void addWaveLengthOHC(const Real wavelength)
Definition: Openholo.h:533
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
bool propagatePupilToRetina()
Calculates 2nd propagation (from pupil plane to retina plane)
ImgDecoderOhc * OHC_decoder
Definition: Openholo.h:495
oph::vec3 getWavelengths()
Returns wavelengths in meter.
bool isReadyToPropagate()
Returns if all data are prepared.
bool propagate()
Do cascaded propagation.
Complex< T > & exp()
Definition: complex.h:395
Real getPixelPitchY()
Returns vertical pixel pitch in meter.
ivec2 pixel_number
Definition: Openholo.h:99
oph::Complex< Real > * getSlmWavefield(oph::uint id)
Return monochromatic wavefield at SLM plane.
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: Openholo.cpp:280
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
virtual void ophFree(void)
Pure virtual function for override in child classes.
virtual bool saveAsOhc(const char *fname)
Function to write OHC file
Definition: Openholo.cpp:252
Complex< Real > ** complex_H
Definition: Openholo.h:489
oph::Complex< Real > * getPupilWavefield(oph::uint id)
Return monochromatic wavefield at pupil plane.
Real getPixelPitchX()
Returns horizontal pixel pitch in meter.
OphConfig context_
Definition: Openholo.h:485
Real * wave_length
Definition: Openholo.h:106
Real getDistObjectToPupil()
Returns distance from reconstruction plane to pupil plane in meter.
oph::uint getNumColors()
Returns number of colors.
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
unsigned int uint
Definition: typedef.h:62
#define CUR_TIME
Definition: function.h:58
Real getNor()
Returns Nor, which affects the range of output intensity.