Openholo  v1.0
Open Source Digital Holographic Library
ophWRP.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 "ophwrp.h"
47 #include "sys.h"
48 #include "tinyxml2.h"
49 
51  : ophGen()
52  , scaledVertex(nullptr)
53 {
54  n_points = -1;
55  p_wrp_ = nullptr;
56  is_CPU = true;
57  is_ViewingWindow = false;
58 }
59 
61 {
62 }
63 
64 void ophWRP::setViewingWindow(bool is_ViewingWindow)
65 {
66  this->is_ViewingWindow = is_ViewingWindow;
67 }
68 
70 {
71 #ifdef CHECK_PROC_TIME
72  auto begin = CUR_TIME;
73 #endif
74  const uint pnX = context_.pixel_number[_X];
75  const uint pnY = context_.pixel_number[_Y];
76  const Real ppX = context_.pixel_pitch[_X];//wrp pitch
77  const Real ppY = context_.pixel_pitch[_Y];
78 
79  Real size = pnY * ppY * 0.8 / 2.0;
80 
82  if (scaledVertex) {
83  delete[] scaledVertex;
84  scaledVertex = nullptr;
85  }
86  scaledVertex = new Real[n_points * 3];
87  memset(scaledVertex, 0.0, sizeof(Real) * n_points * 3);
88 
89  if (is_ViewingWindow) {
90  transVW(scaledVertex, pc.vertex, n_points * 3);
91  }
92 
93  Real* x = new Real[n_points];
94  Real* y = new Real[n_points];
95  Real* z = new Real[n_points];
96 
97  int i;
98  int num_threads = 1;
99 #ifdef _OPENMP
100 #pragma omp parallel
101  {
102  num_threads = omp_get_num_threads(); // get number of Multi Threading
103 #pragma omp for private(i)
104 #endif
105  for (i = 0; i < n_points; i++) {
106  uint idx = 3 * i;
107  x[i] = is_ViewingWindow ? scaledVertex[idx + _X] : pc.vertex[idx + _X];
108  y[i] = is_ViewingWindow ? scaledVertex[idx + _Y] : pc.vertex[idx + _Y];
109  z[i] = is_ViewingWindow ? scaledVertex[idx + _Z] : pc.vertex[idx + _Z];
110  }
111 #ifdef _OPENMP
112  }
113 #endif
114  Real xmax = maxOfArr(x, n_points);
115  Real ymax = maxOfArr(y, n_points);
116  Real zmax = maxOfArr(z, n_points);
117 
118  int j;
119 #ifdef _OPENMP
120 #pragma omp parallel
121  {
122  num_threads = omp_get_num_threads(); // get number of Multi Threading
123 #pragma omp for private(j)
124 #endif
125  for (j = 0; j < n_points; ++j) { //Create Fringe Pattern
126  uint idx = 3 * j;
127  Real maxXY = (xmax > ymax) ? xmax : ymax;
128  Real *pSrc = is_ViewingWindow ? scaledVertex : pc.vertex;
129 
130  Real tmpX = pSrc[idx + _X] / maxXY * size;
131  Real tmpY = pSrc[idx + _Y] / maxXY * size;
132  Real tmpZ = pSrc[idx + _Z] / zmax * size;
133 
134  scaledVertex[idx + _X] = tmpX;
135  scaledVertex[idx + _Y] = tmpY;
136  scaledVertex[idx + _Z] = tmpZ;
137  z[j] = tmpZ;// pc.vertex[idx + _Z];
138  }
139 #ifdef _OPENMP
140  }
141 #endif
142  zmax_ = maxOfArr(z, n_points);
143 
144  delete[] x;
145  delete[] y;
146  delete[] z;
147 #ifdef CHECK_PROC_TIME
148  auto end = CUR_TIME;
149  LOG("\n%s : %lf(s) <%d threads>\n\n",
150  __FUNCTION__,
151  ((chrono::duration<Real>)(end - begin)).count(),
152  num_threads);
153 #endif
154 }
155 
156 int ophWRP::loadPointCloud(const char* pc_file)
157 {
158  n_points = ophGen::loadPointCloud(pc_file, &obj_);
159 
160  return n_points;
161 
162 }
163 
164 bool ophWRP::readConfig(const char* fname)
165 {
166  if (!ophGen::readConfig(fname))
167  return false;
168 
169  LOG("Reading....%s...", fname);
170 
171  auto start = CUR_TIME;
172 
173  using namespace tinyxml2;
174  /*XML parsing*/
175  tinyxml2::XMLDocument xml_doc;
176  XMLNode *xml_node;
177 
178  if (!checkExtension(fname, ".xml"))
179  {
180  LOG("file's extension is not 'xml'\n");
181  return false;
182  }
183  auto ret = xml_doc.LoadFile(fname);
184  if (ret != XML_SUCCESS)
185  {
186  LOG("Failed to load file \"%s\"\n", fname);
187  return false;
188  }
189 
190  xml_node = xml_doc.FirstChild();
191 
192  // about viewing window
193  auto next = xml_node->FirstChildElement("FieldLength");
194  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.fieldLength))
195  return false;
196 
197  // about point
198  next = xml_node->FirstChildElement("ScaleX");
199  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.scale[_X]))
200  return false;
201  next = xml_node->FirstChildElement("ScaleY");
202  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.scale[_Y]))
203  return false;
204  next = xml_node->FirstChildElement("ScaleZ");
205  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.scale[_Z]))
206  return false;
207  next = xml_node->FirstChildElement("Distance");
208  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.propagation_distance))
209  return false;
210  next = xml_node->FirstChildElement("LocationOfWRP");
211  if (!next || XML_SUCCESS != next->QueryDoubleText(&wrp_config_.wrp_location))
212  return false;
213  next = xml_node->FirstChildElement("NumOfWRP");
214  if (!next || XML_SUCCESS != next->QueryIntText(&wrp_config_.num_wrp))
215  return false;
216 
217  auto end = CUR_TIME;
218  auto during = ((chrono::duration<Real>)(end - start)).count();
219  LOG("%lf (s)..done\n", during);
220 
221  initialize();
222  return true;
223 }
224 
225 void ophWRP::addPixel2WRP(int x, int y, Complex<Real> temp)
226 {
227  const uint pnX = context_.pixel_number[_X];
228  const uint pnY = context_.pixel_number[_Y];
229 
230  if (x >= 0 && x < pnX && y >= 0 && y < pnY) {
231  uint adr = x + y * pnX;
232  if (adr == 0)
233  std::cout << ".0";
234  p_wrp_[adr] = p_wrp_[adr] + temp;
235  }
236 
237 }
238 
239 void ophWRP::addPixel2WRP(int x, int y, oph::Complex<Real> temp, oph::Complex<Real>* wrp)
240 {
241  long long int Nx = context_.pixel_number.v[0];
242  long long int Ny = context_.pixel_number.v[1];
243 
244  if (x >= 0 && x<Nx && y >= 0 && y< Ny) {
245  long long int adr = x + y*Nx;
246  wrp[adr] += temp[adr];
247  }
248 }
249 
250 oph::Complex<Real>* ophWRP::calSubWRP(double wrp_d, Complex<Real>* wrp, OphPointCloudData* pc)
251 {
252 
253  Real wave_num = context_.k; // wave_number
254  Real wave_len = context_.wave_length[0]; //wave_length
255 
256  int Nx = context_.pixel_number.v[0]; //slm_pixelNumberX
257  int Ny = context_.pixel_number.v[1]; //slm_pixelNumberY
258 
259  Real wpx = context_.pixel_pitch.v[0];//wrp pitch
260  Real wpy = context_.pixel_pitch.v[1];
261 
262 
263  int Nx_h = Nx >> 1;
264  int Ny_h = Ny >> 1;
265 
266  int num = n_points;
267 
268 
269 #ifdef _OPENMP
270  omp_set_num_threads(omp_get_num_threads());
271 #pragma omp parallel for
272 #endif
273 
274  for (int k = 0; k < num; k++) {
275 
276  uint idx = 3 * k;
277  uint color_idx = pc->n_colors * k;
278 
279  Real x = pc->vertex[idx + _X];
280  Real y = pc->vertex[idx + _Y];
281  Real z = pc->vertex[idx + _Z];
282  Real amplitude = pc->color[color_idx];
283 
284  float dz = wrp_d - z;
285  // float tw = (int)fabs(wave_len*dz / wpx / wpx / 2 + 0.5) * 2 - 1;
286  float tw = fabs(dz)*wave_len / wpx / wpx / 2;
287 
288  int w = (int)tw;
289 
290  int tx = (int)(x / wpx) + Nx_h;
291  int ty = (int)(y / wpy) + Ny_h;
292 
293  printf("num=%d, tx=%d, ty=%d, w=%d\n", k, tx, ty, w);
294 
295  for (int wy = -w; wy < w; wy++) {
296  for (int wx = -w; wx<w; wx++) {//WRP coordinate
297 
298  double dx = wx*wpx;
299  double dy = wy*wpy;
300  double dz = wrp_d - z;
301 
302  double sign = (dz>0.0) ? (1.0) : (-1.0);
303  double r = sign*sqrt(dx*dx + dy*dy + dz*dz);
304 
305  //double tmp_re,tmp_im;
306  Complex<Real> tmp;
307 
308  tmp._Val[_RE] = (amplitude*cosf(wave_num*r)*cosf(wave_num*wave_len*rand(0, 1))) / (r + 0.05);
309  tmp._Val[_IM] = (-amplitude*sinf(wave_num*r)*sinf(wave_num*wave_len*rand(0, 1))) / (r + 0.05);
310 
311  if (tx + wx >= 0 && tx + wx < Nx && ty + wy >= 0 && ty + wy < Ny)
312  {
313  addPixel2WRP(wx + tx, wy + ty, tmp, wrp);
314  }
315  }
316  }
317  }
318 
319  return wrp;
320 }
321 
322 void ophWRP::setMode(bool isCPU)
323 {
324  is_CPU = isCPU;
325 }
326 
328 {
329 #ifdef CHECK_PROC_TIME
330  auto begin = CUR_TIME;
331 #endif
332 
333  const uint pnX = context_.pixel_number[_X]; //slm_pixelNumberX
334  const uint pnY = context_.pixel_number[_Y]; //slm_pixelNumberY
335  const Real ppX = context_.pixel_pitch[_X]; //wrp pitch
336  const Real ppY = context_.pixel_pitch[_Y];
337  const uint pnXY = pnX * pnY;
338  const uint nChannel = context_.waveNum;
339  const Real distance = wrp_config_.propagation_distance;
340  int pnX_h = pnX >> 1;
341  int pnY_h = pnY >> 1;
342 
343  OphPointCloudData pc = obj_;
344  Real wrp_d = wrp_config_.wrp_location;
345 
346  // Memory Location for Result Image
347 
348  if (p_wrp_) {
349  delete[] p_wrp_;
350  p_wrp_ = nullptr;
351  }
352  p_wrp_ = new Complex<Real>[pnXY];
353  memset(p_wrp_, 0.0, sizeof(Complex<Real>) * pnXY);
354 
355  int num_threads = 1;
356 
357  for (uint ch = 0; ch < nChannel; ch++) {
358 
359  Real lambda = context_.wave_length[ch]; //wave_length
360  Real k = context_.k = 2 * M_PI / lambda;
361 
362  int i;
363 #ifdef _OPENMP
364 #pragma omp parallel
365  {
366  num_threads = omp_get_num_threads();
367 #pragma omp parallel for private(i)
368 #endif
369  for (i = 0; i < n_points; ++i) {
370  uint idx = 3 * i;
371  uint color_idx = pc.n_colors * i;
372 
373  Real x = scaledVertex[idx + _X];
374  Real y = scaledVertex[idx + _Y];
375  Real z = scaledVertex[idx + _Z];
376  Real amplitude = pc.color[color_idx];
377 
378  float dz = wrp_d - z;
379  float tw = (int)fabs(lambda * dz / ppX / ppX / 2 + 0.5) * 2 - 1;
380  // float tw = fabs(dz)*wave_len / ppX / ppX / 2;
381 
382  int w = (int)tw;
383 
384  int tx = (int)(x / ppX) + pnX_h;
385  int ty = (int)(y / ppY) + pnY_h;
386 
387  //cout << "num = " << k << ", tx = " << tx << ", ty = " << ty << ", w = " << w << endl;
388 
389  for (int wy = -w; wy < w; wy++) {
390  for (int wx = -w; wx < w; wx++) {//WRP coordinate
391 
392  double dx = wx * ppX;
393  double dy = wy * ppY;
394  double dz = wrp_d - z;
395 
396  double sign = (dz > 0.0) ? (1.0) : (-1.0);
397  double r = sign * sqrt(dx*dx + dy * dy + dz * dz);
398 
399  Complex<Real> tmp;
400  tmp[_RE] = (amplitude * cosf(k * r) * cosf(k * lambda * rand(0, 1))) / r;
401  tmp[_IM] = (-amplitude * sinf(k * r) * sinf(k * lambda * rand(0, 1))) / r;
402  if (tx + wx >= 0 && tx + wx < pnX && ty + wy >= 0 && ty + wy < pnY) {
403  int tmpX = wx + tx;
404  int tmpY = wy + ty;
405 
406  if (tmpX >= 0 && tmpX < pnX && tmpY >= 0 && tmpY < pnY) {
407  uint adr = tmpX + tmpY * pnX;
408  if (adr == 0)
409  std::cout << ".0";
410 
411 
412 #ifdef _OPENMP
413 #pragma omp atomic
414  p_wrp_[adr][_RE] += tmp[_RE];
415 #pragma omp atomic
416  p_wrp_[adr][_IM] += tmp[_IM];
417 #else
418  p_wrp_[adr] += tmp;
419 #endif
420  }
421  }
422  }
423  }
424  }
425 #ifdef _OPENMP
426  }
427 #endif
428  fresnelPropagation(p_wrp_, complex_H[ch], distance, ch);
429  memset(p_wrp_, 0.0, sizeof(Complex<Real>) * pnXY);
430  }
431  delete[] p_wrp_;
432  delete[] scaledVertex;
433  p_wrp_ = nullptr;
434  scaledVertex = nullptr;
435 
436 #ifdef CHECK_PROC_TIME
437  auto end = CUR_TIME;
438  LOG("\n%s : %lf(s) <%d threads>\n\n",
439  __FUNCTION__,
440  ((chrono::duration<Real>)(end - begin)).count(),
441  num_threads);
442 #endif
443  return 0.;
444 }
445 
447 {
448  resetBuffer();
449 
450  auto start_time = CUR_TIME;
451  LOG("1) Algorithm Method : WRP\n");
452  LOG("2) Generate Hologram with %s\n", is_CPU ?
453 #ifdef _OPENMP
454  "Multi Core CPU" :
455 #else
456  "Single Core CPU" :
457 #endif
458  "GPU");
459  LOG("3) Transform Viewing Window : %s\n", is_ViewingWindow ? "ON" : "OFF");
460 
461  auto begin = CUR_TIME;
462 
463  autoScaling();
464  is_CPU ? calculateWRPCPU() : calculateWRPGPU();
465 
466  auto end = CUR_TIME;
467  elapsedTime = ((std::chrono::duration<Real>)(end - begin)).count();
468 
469  LOG("Total Elapsed Time: %lf (s)\n", elapsedTime);
470 }
471 
472 Complex<Real>** ophWRP::calculateMWRP(void)
473 {
474  int wrp_num = wrp_config_.num_wrp;
475 
476  if (wrp_num < 1)
477  return nullptr;
478 
479  Complex<Real>** wrp_list = nullptr;
480 
481  Real k = context_.k; // wave_number
482  Real lambda = context_.wave_length[0]; //wave_length
483 
484  const uint pnX = context_.pixel_number[_X]; //slm_pixelNumberX
485  const uint pnY = context_.pixel_number[_Y]; //slm_pixelNumberY
486  const uint pnXY = pnX * pnY;
487  const Real ppX = context_.pixel_pitch[_X];//wrp pitch
488  const Real ppY = context_.pixel_pitch[_Y];
489 
490  Complex<Real>* wrp = new Complex<Real>[pnXY];
491 
492  OphPointCloudData pc = obj_;
493 
494  for (int i = 0; i<wrp_num; i++)
495  {
496 // wrp = calSubWRP(wrp_d, wrp, &pc);
497  wrp_list[i] = wrp;
498  }
499 
500  return wrp_list;
501 }
502 
503 void ophWRP::ophFree(void)
504 {
505  // delete[] obj_.vertex;
506  // delete[] obj_.color;
507 
508 }
509 
510 void ophWRP::transVW(Real* dst, Real* src, int size)
511 {
512  Real fieldLens = getFieldLens();
513  for (int i = 0; i < size; i++) {
514  dst[i] = (-fieldLens * src[i]) / (src[i] - fieldLens);
515  }
516 }
Real * vertex
Geometry of point clouds.
Definition: ophGen.h:433
Real * color
Color data of point clouds.
Definition: ophGen.h:435
Real k
Definition: Openholo.h:66
void setMode(bool is_CPU)
Set the value of a variable is_CPU(true or false)
Definition: ophWRP.cpp:322
Real * wave_length
Definition: Openholo.h:69
Real * scaledVertex
Definition: ophWRP.h:233
Real fieldLength
fieldLength variable for viewing window.
Definition: ophGen.h:513
OphPointCloudData obj_
Input Pointcloud Data.
Definition: ophWRP.h:232
float Real
Definition: typedef.h:55
void initialize(void)
Initialize variables for Hologram complex field, encoded data, normalized data.
Definition: ophGen.cpp:69
void setViewingWindow(bool is_ViewingWindow)
Set the value of a variable is_ViewingWindow(true or false)
Definition: ophWRP.cpp:64
#define CUR_TIME
Definition: function.h:58
int num_wrp
Number of wavefront recording plane(WRP)
Definition: ophGen.h:517
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:80
ophWRP(void)
Constructor.
Definition: ophWRP.cpp:50
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
virtual bool readConfig(const char *fname)
load to configuration file.
Definition: ophWRP.cpp:164
Real elapsedTime
Elapsed time of generate hologram.
Definition: ophGen.h:246
double calculateWRPGPU(void)
Definition: ophWRP_GPU.cpp:47
#define _Y
Definition: define.h:84
#define _IM
Definition: complex.h:57
vec3 scale
Scaling factor of coordinate of point cloud.
Definition: ophGen.h:515
#define _X
Definition: define.h:80
virtual int loadPointCloud(const char *pc_file)
load to point cloud data.
Definition: ophWRP.cpp:156
Real wrp_location
Location distance of WRP.
Definition: ophGen.h:519
oph::ivec2 pixel_number
Definition: Openholo.h:63
OphWRPConfig wrp_config_
structure variable for WRP hologram configuration
Definition: ophWRP.h:234
#define _RE
Definition: complex.h:54
void fresnelPropagation(OphConfig context, Complex< Real > *in, Complex< Real > *out, Real distance)
Fresnel propagation.
Definition: ophGen.cpp:768
virtual ~ophWRP(void)
Destructor.
Definition: ophWRP.cpp:60
uint waveNum
Definition: Openholo.h:68
int n_colors
Number of color channel.
Definition: ophGen.h:431
bool readConfig(const char *fname)
load to configuration file.
Definition: ophGen.cpp:150
double calculateWRPCPU(void)
Definition: ophWRP.cpp:327
Complex< Real > * p_wrp_
wrp buffer - complex type
Definition: ophWRP.h:230
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2157
Real rand(const Real min, const Real max, oph::ulong _SEED_VALUE=0)
Get random Real value from min to max.
Definition: function.h:287
int n_points
numbers of points
Definition: ophWRP.h:228
void autoScaling()
Definition: ophWRP.cpp:69
Data for Point Cloud.
Definition: ophGen.h:427
Real propagation_distance
Distance of Hologram plane.
Definition: ophGen.h:521
Real maxOfArr(const std::vector< Real > &arr)
Definition: function.h:86
void resetBuffer()
reset buffer
Definition: ophGen.cpp:363
OphConfig context_
Definition: Openholo.h:297
Complex< Real > ** complex_H
Definition: Openholo.h:298
#define _Z
Definition: define.h:88
int loadPointCloud(const char *pc_file, OphPointCloudData *pc_data_)
load to point cloud data.
Definition: ophGen.cpp:133
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
const Real & getFieldLens()
Definition: ophWRP.h:143
unsigned int uint
Definition: typedef.h:62
#define M_PI
Definition: define.h:52
Definition: ophGen.h:68
void generateHologram(void)
Generate a hologram, main funtion.
Definition: ophWRP.cpp:446
oph::vec2 pixel_pitch
Definition: Openholo.h:64
Complex< Real > ** calculateMWRP(void)
Generate multiple wavefront recording planes, main funtion.
Definition: ophWRP.cpp:472