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