Openholo  v1.0
Open Source Digital Holographic Library
ophGen.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 "ophGen.h"
47 #include <windows.h>
48 #include "sys.h"
49 #include "function.h"
50 #include <cuda_runtime.h>
51 #include <cufft.h>
52 
53 #include "tinyxml2.h"
54 #include "PLYparser.h"
55 
57  : Openholo()
58  , holo_encoded(nullptr)
59  , holo_normalized(nullptr)
60  , nOldChannel(0)
61  , elapsedTime(0.0)
62 {
63 }
64 
66 {
67 }
68 
70 {
71  LOG("%s...\n", __FUNCTION__);
72  // Output Image Size
73  const uint pnX = context_.pixel_number[_X];
74  const uint pnY = context_.pixel_number[_Y];
75  const uint pnXY = pnX * pnY;
76  const int nChannel = context_.waveNum;
77 
78  // Memory Location for Result Image
79  if (complex_H != nullptr) {
80  for (uint i = 0; i < nOldChannel; i++) {
81  if (complex_H[i] != nullptr) {
82  delete[] complex_H[i];
83  complex_H[i] = nullptr;
84  }
85  }
86  delete[] complex_H;
87  complex_H = nullptr;
88  }
89 
90  complex_H = new Complex<Real>*[nChannel];
91  for (uint i = 0; i < nChannel; i++) {
92  complex_H[i] = new Complex<Real>[pnXY];
93  memset(complex_H[i], 0, sizeof(Complex<Real>) * pnXY);
94  }
95 
96  if (holo_encoded != nullptr) {
97  for (uint i = 0; i < nOldChannel; i++) {
98  if (holo_encoded[i] != nullptr) {
99  delete[] holo_encoded[i];
100  holo_encoded[i] = nullptr;
101  }
102  }
103  delete[] holo_encoded;
104  holo_encoded = nullptr;
105  }
106  holo_encoded = new Real*[nChannel];
107  for (uint i = 0; i < nChannel; i++) {
108  holo_encoded[i] = new Real[pnXY];
109  memset(holo_encoded[i], 0, sizeof(Real) * pnXY);
110  }
111 
112  if (holo_normalized != nullptr) {
113  for (uint i = 0; i < nOldChannel; i++) {
114  if (holo_normalized[i] != nullptr) {
115  delete[] holo_normalized[i];
116  holo_normalized[i] = nullptr;
117  }
118  }
119  delete[] holo_normalized;
120  holo_normalized = nullptr;
121  }
122  holo_normalized = new uchar*[nChannel];
123  for (uint i = 0; i < nChannel; i++) {
124  holo_normalized[i] = new uchar[pnXY];
125  memset(holo_normalized[i], 0, sizeof(uchar) * pnXY);
126  }
127 
128  nOldChannel = nChannel;
129  encode_size[_X] = pnX;
130  encode_size[_Y] = pnY;
131 }
132 
133 int ophGen::loadPointCloud(const char* pc_file, OphPointCloudData *pc_data_)
134 {
135  LOG("Reading....%s...\n", pc_file);
136 
137  auto start = CUR_TIME;
138 
139  PLYparser plyIO;
140  if (!plyIO.loadPLY(pc_file, pc_data_->n_points, pc_data_->n_colors, &pc_data_->vertex, &pc_data_->color, &pc_data_->phase, pc_data_->isPhaseParse))
141  return -1;
142 
143  auto end = CUR_TIME;
144  auto during = ((std::chrono::duration<Real>)(end - start)).count();
145 
146  LOG("%.5lfsec...done\n", during);
147  return pc_data_->n_points;
148 }
149 
150 bool ophGen::readConfig(const char* fname)
151 {
152  LOG("Reading....%s...", fname);
153  using namespace tinyxml2;
154 
155  auto start = CUR_TIME;
156  tinyxml2::XMLDocument xml_doc;
157  XMLNode *xml_node = nullptr;
158 
159  if (!checkExtension(fname, ".xml"))
160  {
161  LOG("file's extension is not 'xml'\n");
162  return false;
163  }
164  auto ret = xml_doc.LoadFile(fname);
165  if (ret != XML_SUCCESS)
166  {
167  LOG("Failed to load file \"%s\"\n", fname);
168  return false;
169  }
170  xml_node = xml_doc.FirstChild();
171 
172  int nWave = 1;
173  auto next = xml_node->FirstChildElement("SLM_WaveNum"); // OffsetInDepth
174  if (!next || XML_SUCCESS != next->QueryIntText(&nWave))
175  return false;
176 
177  context_.waveNum = nWave;
179  context_.wave_length = new Real[nWave];
180 
181  char szNodeName[32] = { 0, };
182  for (int i = 1; i <= nWave; i++) {
183  wsprintfA(szNodeName, "SLM_WaveLength_%d", i);
184  next = xml_node->FirstChildElement(szNodeName);
185  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.wave_length[i - 1]))
186  return false;
187  }
188  next = xml_node->FirstChildElement("SLM_PixelNumX");
189  if (!next || XML_SUCCESS != next->QueryIntText(&context_.pixel_number[_X]))
190  return false;
191  next = xml_node->FirstChildElement("SLM_PixelNumY");
192  if (!next || XML_SUCCESS != next->QueryIntText(&context_.pixel_number[_Y]))
193  return false;
194  next = xml_node->FirstChildElement("SLM_PixelPitchX");
195  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.pixel_pitch[_X]))
196  return false;
197  next = xml_node->FirstChildElement("SLM_PixelPitchY");
198  if (!next || XML_SUCCESS != next->QueryDoubleText(&context_.pixel_pitch[_Y]))
199  return false;
200 
203 
206 
208  for (int i = 0; i < nWave; i++)
210 
211  auto end = CUR_TIME;
212  auto during = ((chrono::duration<Real>)(end - start)).count();
213  LOG("%lf (s)...done\n", during);
214 
215  return true;
216 }
217 
218 void ophGen::propagationAngularSpectrum(int ch, Complex<Real>* input_u, Real propagation_dist, Real k, Real lambda)
219 {
220  const int pnX = context_.pixel_number[_X];
221  const int pnY = context_.pixel_number[_Y];
222  const Real ppX = context_.pixel_pitch[_X];
223  const Real ppY = context_.pixel_pitch[_Y];
224  const Real ssX = context_.ss[_X] = pnX * ppX;
225  const Real ssY = context_.ss[_Y] = pnY * ppY;
226 
227  for (int i = 0; i < pnX * pnY; i++) {
228  Real x = i % pnX;
229  Real y = i / pnX;
230 
231  Real fxx = (-1.0 / (2.0*ppX)) + (1.0 / ssX) * x;
232  Real fyy = (1.0 / (2.0*ppY)) - (1.0 / ssY) - (1.0 / ssY) * y;
233 
234  Real sval = sqrt(1 - (lambda*fxx)*(lambda*fxx) - (lambda*fyy)*(lambda*fyy));
235  sval *= k * propagation_dist;
236  Complex<Real> kernel(0, sval);
237  kernel.exp();
238 
239  int prop_mask = ((fxx * fxx + fyy * fyy) < (k * k)) ? 1 : 0;
240 
241  Complex<Real> u_frequency;
242  if (prop_mask == 1)
243  u_frequency = kernel * input_u[i];
244 
245 #pragma omp atomic
246  complex_H[ch][i][_RE] += u_frequency[_RE];
247 #pragma omp atomic
248  complex_H[ch][i][_IM] += u_frequency[_IM];
249  }
250 }
251 
253 {
254  for (uint ch = 0; ch < context_.waveNum; ch++)
256 }
257 
258 bool ophGen::save(const char * fname, uint8_t bitsperpixel, uchar* src, uint px, uint py)
259 {
260  if (fname == nullptr) return false;
261 
262  uchar* source = src;
263  const uint nChannel = context_.waveNum;
264  ivec2 p(px, py);
265  if (px == 0 && py == 0)
267 
268  char path[_MAX_PATH] = { 0, };
269  char drive[_MAX_DRIVE] = { 0, };
270  char dir[_MAX_DIR] = { 0, };
271  char file[_MAX_FNAME] = { 0, };
272  char ext[_MAX_EXT] = { 0, };
273  _splitpath_s(fname, drive, dir, file, ext);
274 
275  sprintf_s(path, "%s", fname);
276 
277  for (uint ch = 0; ch < nChannel; ch++) {
278  if (src == nullptr)
279  source = holo_normalized[ch];
280  if (nChannel > 1) {
281  sprintf_s(path, "%s%s%s_%d%s", drive, dir, file, ch + 1, ext);
282  }
283 
284  if (checkExtension(path, ".bmp")) // when the extension is bmp
285  Openholo::saveAsImg(path, bitsperpixel, source, p[_X], p[_Y]);
286  else if (
287  checkExtension(path, ".jpg") ||
288  checkExtension(path, ".gif") ||
289  checkExtension(path, ".png")) {
290  Openholo::saveAsImg(path, bitsperpixel, source, p[_X], p[_Y]);
291  }
292  else { // when extension is not .ohf, .bmp - force bmp
293  char buf[256];
294  memset(buf, 0x00, sizeof(char) * 256);
295  sprintf_s(buf, "%s.bmp", path);
296 
297  Openholo::saveAsImg(buf, bitsperpixel, source, p[_X], p[_Y]);
298  }
299  }
300  return true;
301 }
302 
303 bool ophGen::save(const char * fname, uint8_t bitsperpixel, uint px, uint py, uint fnum, uchar* args ...)
304 {
305  std::string file = fname;
306  std::string name;
307  std::string ext;
308 
309  size_t ex = file.rfind(".");
310  if (ex == -1) ex = file.length();
311 
312  name = file.substr(0, ex);
313  ext = file.substr(ex, file.length() - 1);
314 
315  va_list ap;
316  __crt_va_start(ap, args);
317 
318  for (uint i = 0; i < fnum; i++) {
319  name.append(std::to_string(i)).append(ext);
320  if (i == 0) {
321  save(name.c_str(), bitsperpixel, args, px, py);
322  continue;
323  }
324  uchar* data = __crt_va_arg(ap, uchar*);
325  save(name.c_str(), bitsperpixel, data, px, py);
326  }
327 
328  __crt_va_end(ap);
329 
330  return true;
331 }
332 
333 void* ophGen::load(const char * fname)
334 {
335  if (checkExtension(fname, ".bmp")) {
336  return Openholo::loadAsImg(fname);
337  }
338  else { // when extension is not .bmp
339  return nullptr;
340  }
341 
342  return nullptr;
343 }
344 
345 bool ophGen::loadAsOhc(const char * fname)
346 {
347  if (!Openholo::loadAsOhc(fname)) return false;
348 
349  const uint nChannel = context_.waveNum;
351 
352  holo_encoded = new Real*[nChannel];
353  holo_normalized = new uchar*[nChannel];
354  for (uint ch = 0; ch < nChannel; ch++) {
355  holo_encoded[ch] = new Real[pnXY];
356  memset(holo_encoded, 0, sizeof(Real) * pnXY);
357  holo_normalized[ch] = new uchar[pnXY];
358  memset(holo_normalized, 0, sizeof(uchar) * pnXY);
359  }
360  return true;
361 }
362 
364 {
365  const uint pnX = context_.pixel_number[_X];
366  const uint pnY = context_.pixel_number[_Y];
367  const uint pnXY = pnX * pnY;
368 
369  for (uint ch = 0; ch < context_.waveNum; ch++) {
370  if(complex_H[ch])
371  memset(complex_H[ch], 0., sizeof(Complex<Real>) * pnXY);
372  if (holo_encoded[ch])
373  memset(holo_encoded[ch], 0., sizeof(Real) * encode_size[_X] * encode_size[_Y]);
374  if (holo_normalized[ch])
375  memset(holo_normalized[ch], 0, sizeof(uchar) * encode_size[_X] * encode_size[_Y]);
376  }
377 
378 }
379 
380 #define for_i(itr, oper) for(int i=0; i<itr; i++){ oper }
381 
382 
383 void ophGen::encoding(unsigned int ENCODE_FLAG, Complex<Real>* holo, bool bShift)
384 {
385  bool bChangeSize = false;
386  const uint pnX = context_.pixel_number[_X];
387  const uint pnY = context_.pixel_number[_Y];
388  const uint nChannel = context_.waveNum;
389  const uint pnXY = pnX * pnY;
390  Complex<Real>* dst = nullptr;
391  if (bShift) {
392  dst = new Complex<Real>[pnXY];
393  }
394 
395  for (int ch = 0; ch < nChannel; ch++) {
396  if (bShift) {
397  memset(dst, 0.0, sizeof(Complex<Real>) * pnXY);
399  fftwShift(complex_H[ch], dst, pnX, pnY, OPH_BACKWARD);
400  holo = complex_H[ch];
401  }
402  else
403  holo = complex_H[ch];
404 
406  encode_size[_X] = pnX * 3;
407  encode_size[_Y] = pnY;
408  bChangeSize = true;
409  }
410  else if (ENCODE_FLAG == ENCODE_TWOPHASE) {
411  encode_size[_X] = pnX * 2;
412  encode_size[_Y] = pnY;
413  bChangeSize = true;
414  }
415  else {
416  encode_size[_X] = pnX;
417  encode_size[_Y] = pnY;
418  }
419 
420  /* initialize */
421  if (bChangeSize) {
422  if (holo_encoded != nullptr) delete[] holo_encoded;
423  if (holo_normalized != nullptr) delete[] holo_normalized;
424 
425  holo_encoded[ch] = new Real[encode_size[_X] * encode_size[_Y]];
427  }
428 
429  memset(holo_encoded[ch], 0, sizeof(Real) * encode_size[_X] * encode_size[_Y]);
430  memset(holo_normalized[ch], 0, sizeof(uchar) * encode_size[_X] * encode_size[_Y]);
431 
432  switch (ENCODE_FLAG)
433  {
434  case ENCODE_SIMPLENI:
435  LOG("Simple Numerical Interference Encoding..");
436  numericalInterference((holo), holo_encoded[ch], pnXY);
437  LOG("Done.\n.");
438  break;
439  case ENCODE_REAL:
440  LOG("Real Part Encoding..");
441  realPart<Real>((holo), holo_encoded[ch], pnXY);
442  LOG("Done.\n.");
443  break;
444  case ENCODE_BURCKHARDT:
445  LOG("Burckhardt Encoding..");
446  burckhardt((holo), holo_encoded[ch], pnXY);
447  LOG("Done.\n.");
448  break;
449  case ENCODE_TWOPHASE:
450  LOG("Two Phase Encoding..");
451  twoPhaseEncoding((holo), holo_encoded[ch], pnXY);
452  LOG("Done.\n.");
453  break;
454  case ENCODE_PHASE:
455  LOG("Phase Encoding..");
456  getPhase((holo), holo_encoded[ch], pnXY);
457  LOG("Done.\n.");
458  break;
459  case ENCODE_AMPLITUDE:
460  LOG("Amplitude Encoding..");
461  getAmplitude((holo), holo_encoded[ch], pnXY);
462  LOG("Done.\n.");
463  break;
464  case ENCODE_SSB:
465  case ENCODE_OFFSSB:
466  LOG("error: PUT PASSBAND\n");
467  cin.get();
468  return;
470  LOG("Symmetrization Encoding..");
471  encodeSymmetrization((holo), holo_encoded[ch], ivec2(0, 1));
472  LOG("Done.\n.");
473  break;
474  default:
475  LOG("error: WRONG ENCODE_FLAG\n");
476  cin.get();
477  return;
478  }
479  }
480  if(dst)
481  delete[] dst;
482 }
483 
484 void ophGen::encoding(unsigned int ENCODE_FLAG, unsigned int passband, Complex<Real>* holo)
485 {
486  holo == nullptr ? holo = *complex_H : holo;
487 
488  const uint pnX = encode_size[_X] = context_.pixel_number[_X];
489  const uint pnY = encode_size[_Y] = context_.pixel_number[_Y];
490  const uint pnXY = pnX * pnY;
491  const uint nChannel = context_.waveNum;
492 
493  for (uint ch = 0; ch < nChannel; ch++) {
494  /* initialize */
495  int encode_size = pnXY;
496  if (holo_encoded[ch] != nullptr) delete[] holo_encoded[ch];
497  holo_encoded[ch] = new Real[encode_size];
498  memset(holo_encoded[ch], 0, sizeof(Real) * encode_size);
499 
500  if (holo_normalized[ch] != nullptr) delete[] holo_normalized[ch];
501  holo_normalized[ch] = new uchar[encode_size];
502  memset(holo_normalized[ch], 0, sizeof(uchar) * encode_size);
503 
504  switch (ENCODE_FLAG)
505  {
506  case ENCODE_SSB:
507  LOG("Single Side Band Encoding..");
508  singleSideBand((holo), holo_encoded[ch], context_.pixel_number, passband);
509  LOG("Done.");
510  break;
511  case ENCODE_OFFSSB:
512  LOG("Off-axis Single Side Band Encoding..");
514  singleSideBand((holo), holo_encoded[ch], context_.pixel_number, passband);
515  LOG("Done.\n");
516  break;
517  default:
518  LOG("error: WRONG ENCODE_FLAG\n");
519  cin.get();
520  return;
521  }
522  }
523 }
524 
526 {
527  const uint pnX = encode_size[_X] = context_.pixel_number[_X];
528  const uint pnY = encode_size[_Y] = context_.pixel_number[_Y];
529  const uint pnXY = pnX * pnY;
530  const uint nChannel = context_.waveNum;
531 
533  else if (ENCODE_METHOD == ENCODE_TWOPHASE) encode_size[_X] *= 2;
534 
535  for (uint ch = 0; ch < nChannel; ch++) {
536  /* initialize */
537  if (holo_encoded[ch] != nullptr) delete[] holo_encoded[ch];
538  holo_encoded[ch] = new Real[encode_size[_X] * encode_size[_Y]];
539  memset(holo_encoded[ch], 0, sizeof(Real) * encode_size[_X] * encode_size[_Y]);
540 
541  if (holo_normalized[ch] != nullptr) delete[] holo_normalized[ch];
543  memset(holo_normalized[ch], 0, sizeof(uchar) * encode_size[_X] * encode_size[_Y]);
544 
545 
546  switch (ENCODE_METHOD)
547  {
548  case ENCODE_SIMPLENI:
549  cout << "Simple Numerical Interference Encoding.." << endl;
551  break;
552  case ENCODE_REAL:
553  cout << "Real Part Encoding.." << endl;
554  realPart<Real>(complex_H[ch], holo_encoded[ch], pnXY);
555  break;
556  case ENCODE_BURCKHARDT:
557  cout << "Burckhardt Encoding.." << endl;
558  burckhardt(complex_H[ch], holo_encoded[ch], pnXY);
559  break;
560  case ENCODE_TWOPHASE:
561  cout << "Two Phase Encoding.." << endl;
562  twoPhaseEncoding(complex_H[ch], holo_encoded[ch], pnXY);
563  break;
564  case ENCODE_PHASE:
565  cout << "Phase Encoding.." << endl;
566  getPhase(complex_H[ch], holo_encoded[ch], pnXY);
567  break;
568  case ENCODE_AMPLITUDE:
569  cout << "Amplitude Encoding.." << endl;
570  getAmplitude(complex_H[ch], holo_encoded[ch], pnXY);
571  break;
572  case ENCODE_SSB:
573  cout << "Single Side Band Encoding.." << endl;
575  break;
576  case ENCODE_OFFSSB:
577  cout << "Off-axis Single Side Band Encoding.." << endl;
580  break;
582  cout << "Symmetrization Encoding.." << endl;
583  encodeSymmetrization(complex_H[ch], holo_encoded[ch], ivec2(0, 1));
584  break;
585  default:
586  cout << "error: WRONG ENCODE_FLAG" << endl;
587  cin.get();
588  return;
589  }
590  }
591 }
592 
593 void ophGen::numericalInterference(oph::Complex<Real>* holo, Real* encoded, const int size)
594 {
595  Real* temp1 = new Real[size];
596  absCplxArr<Real>(holo, temp1, size);
597 
598  Real* ref = new Real;
599  *ref = maxOfArr(temp1, size);
600 
601  oph::Complex<Real>* temp2 = new oph::Complex<Real>[size];
602  for_i(size,
603  temp2[i] = holo[i] + *ref;
604  );
605 
606  Real* temp3 = new Real[size];
607  absCplxArr<Real>(temp2, temp3, size);
608 
609  for_i(size,
610  encoded[i] = temp3[i] * temp3[i];
611  );
612 
613  delete[] temp1;
614  delete[] temp2;
615  delete[] temp3;
616  delete ref;
617 }
618 
619 void ophGen::twoPhaseEncoding(oph::Complex<Real>* holo, Real* encoded, const int size)
620 {
621  Complex<Real>* normCplx = new Complex<Real>[size];
622  oph::normalize<Real>(holo, normCplx, size);
623 
624  Real* amp = new Real[size];
625  oph::getAmplitude(normCplx, encoded, size);
626 
627  Real* pha = new Real[size];
628  oph::getPhase(normCplx, pha, size);
629 
630  for_i(size, *(pha + i) += M_PI;);
631 
632  Real* delPhase = new Real[size];
633  for_i(size, *(delPhase + i) = acos(*(amp + i)););
634 
635  for_i(size,
636  *(encoded + i * 2) = *(pha + i) + *(delPhase + i);
637  *(encoded + i * 2 + 1) = *(pha + i) - *(delPhase + i);
638  );
639 
640  delete[] normCplx;
641  delete[] amp;
642  delete[] pha;
643  delete[] delPhase;
644 }
645 
646 void ophGen::burckhardt(oph::Complex<Real>* holo, Real* encoded, const int size)
647 {
648  Complex<Real>* norm = new Complex<Real>[size];
649  oph::normalize(holo, norm, size);
650 
651  Real* phase = new Real[size];
652  oph::getPhase(norm, phase, size);
653 
654  Real* ampl = new Real[size];
655  oph::getAmplitude(norm, ampl, size);
656 
657  Real* A1 = new Real[size];
658  memsetArr<Real>(A1, 0, 0, size - 1);
659  Real* A2 = new Real[size];
660  memsetArr<Real>(A2, 0, 0, size - 1);
661  Real* A3 = new Real[size];
662  memsetArr<Real>(A3, 0, 0, size - 1);
663 
664  for_i(size,
665  if (*(phase + i) >= 0 && *(phase + i) < (2 * M_PI / 3))
666  {
667  *(A1 + i) = *(ampl + i)*(cos(*(phase + i)) + sin(*(phase + i)) / sqrt(3));
668  *(A2 + i) = 2 * sin(*(phase + i)) / sqrt(3);
669  }
670  else if (*(phase + i) >= (2 * M_PI / 3) && *(phase + i) < (4 * M_PI / 3))
671  {
672  *(A2 + i) = *(ampl + i)*(cos(*(phase + i) - (2 * M_PI / 3)) + sin(*(phase + i) - (2 * M_PI / 3)) / sqrt(3));
673  *(A3 + i) = 2 * sin(*(phase + i) - (2 * M_PI / 3)) / sqrt(3);
674  }
675  else if (*(phase + i) >= (4 * M_PI / 3) && *(phase + i) < (2 * M_PI))
676  {
677  *(A3 + i) = *(ampl + i)*(cos(*(phase + i) - (4 * M_PI / 3)) + sin(*(phase + i) - (4 * M_PI / 3)) / sqrt(3));
678  *(A1 + i) = 2 * sin(*(phase + i) - (4 * M_PI / 3)) / sqrt(3);
679  }
680  );
681 
682  for_i(size,
683  *(encoded + (3 * i)) = *(A1 + i);
684  *(encoded + (3 * i + 1)) = *(A2 + i);
685  *(encoded + (3 * i + 2)) = *(A3 + i);
686  );
687 
688 }
689 
690 void ophGen::singleSideBand(oph::Complex<Real>* holo, Real* encoded, const ivec2 holosize, int SSB_PASSBAND)
691 {
692  int size = holosize[_X] * holosize[_Y];
693 
694  oph::Complex<Real>* AS = new oph::Complex<Real>[size];
695  fft2(holosize, holo, OPH_FORWARD, OPH_ESTIMATE);
696  fftwShift(holo, AS, holosize[_X], holosize[_Y], OPH_FORWARD, false);
697  //fftExecute(temp);
698 
699  switch (SSB_PASSBAND)
700  {
701  case SSB_LEFT:
702  for (int i = 0; i < holosize[_Y]; i++)
703  {
704  for (int j = holosize[_X] / 2; j < holosize[_X]; j++)
705  {
706  AS[i*holosize[_X] + j] = 0;
707  }
708  }
709  break;
710  case SSB_RIGHT:
711  for (int i = 0; i < holosize[_Y]; i++)
712  {
713  for (int j = 0; j < holosize[_X] / 2; j++)
714  {
715  AS[i*holosize[_X] + j] = 0;
716  }
717  }
718  break;
719  case SSB_TOP:
720  for (int i = size / 2; i < size; i++)
721  {
722  AS[i] = 0;
723  }
724  break;
725 
726  case SSB_BOTTOM:
727  for (int i = 0; i < size / 2; i++)
728  {
729  AS[i] = 0;
730  }
731  break;
732  }
733 
734  oph::Complex<Real>* filtered = new oph::Complex<Real>[size];
735  fft2(holosize, AS, OPH_BACKWARD, OPH_ESTIMATE);
736  fftwShift(AS, filtered, holosize[_X], holosize[_Y], OPH_BACKWARD, false);
737 
738  //fftExecute(filtered);
739 
740 
741  Real* realFiltered = new Real[size];
742  oph::realPart<Real>(filtered, realFiltered, size);
743 
744  oph::normalize(realFiltered, encoded, size);
745 
746  delete[] AS, filtered , realFiltered;
747 }
748 
749 
750 void ophGen::freqShift(oph::Complex<Real>* src, Complex<Real>* dst, const ivec2 holosize, int shift_x, int shift_y)
751 {
752  int size = holosize[_X] * holosize[_Y];
753 
754  Complex<Real>* AS = new oph::Complex<Real>[size];
755  fft2(holosize, src, OPH_FORWARD, OPH_ESTIMATE);
756  fftwShift(src, AS, holosize[_X], holosize[_Y], OPH_FORWARD);
757  //fftExecute(AS);
758 
759  Complex<Real>* shifted = new oph::Complex<Real>[size];
760  circShift<Complex<Real>>(AS, shifted, shift_x, shift_y, holosize.v[_X], holosize.v[_Y]);
761 
762  fft2(holosize, shifted, OPH_BACKWARD, OPH_ESTIMATE);
763  fftwShift(shifted, dst, holosize[_X], holosize[_Y], OPH_BACKWARD);
764  //fftExecute(dst);
765 }
766 
767 
768 void ophGen::fresnelPropagation(OphConfig context, Complex<Real>* in, Complex<Real>* out, Real distance)
769 {
770  const int pnX = context.pixel_number[_X];
771  const int pnY = context.pixel_number[_Y];
772  const uint pnXY = pnX * pnY;
773 
774  Complex<Real>* in2x = new Complex<Real>[pnXY * 4];
775  Complex<Real> zero(0, 0);
776  memsetArr<Complex<Real>>(in2x, zero, 0, pnXY * 4 - 1);
777 
778  uint idxIn = 0;
779 
780  for (int idxNy = pnY / 2; idxNy < pnY + (pnY / 2); idxNy++) {
781  for (int idxNx = pnX / 2; idxNx < pnX + (pnX / 2); idxNx++) {
782  in2x[idxNy * pnX * 2 + idxNx] = in[idxIn];
783  idxIn++;
784  }
785  }
786 
787  Complex<Real>* temp1 = new Complex<Real>[pnXY * 4];
788 
789  fft2({ pnX * 2, pnY * 2 }, in2x, OPH_FORWARD, OPH_ESTIMATE);
790  fftwShift(in2x, temp1, pnX, pnY, OPH_FORWARD);
791  //fftExecute(temp1);
792 
793  Real* fx = new Real[pnXY * 4];
794  Real* fy = new Real[pnXY * 4];
795 
796  uint i = 0;
797  for (int idxFy = -pnY; idxFy < pnY; idxFy++) {
798  for (int idxFx = -pnX; idxFx < pnX; idxFx++) {
799  fx[i] = idxFx / (2 * pnX * context.pixel_pitch[_X]);
800  fy[i] = idxFy / (2 * pnY * context.pixel_pitch[_Y]);
801  i++;
802  }
803  }
804 
805  Complex<Real>* prop = new Complex<Real>[pnXY * 4];
806  memsetArr<Complex<Real>>(prop, zero, 0, pnXY * 4 - 1);
807 
808  Real sqrtPart;
809 
810  Complex<Real>* temp2 = new Complex<Real>[pnXY * 4];
811 
812  for (int i = 0; i < pnXY * 4; i++) {
813  sqrtPart = sqrt(1 / (context.wave_length[0] * context.wave_length[0]) - fx[i] * fx[i] - fy[i] * fy[i]);
814  prop[i][_IM] = 2 * M_PI * distance;
815  prop[i][_IM] *= sqrtPart;
816  temp2[i] = temp1[i] * exp(prop[i]);
817  }
818 
819  Complex<Real>* temp3 = new Complex<Real>[pnXY * 4];
820  fft2({ pnX * 2, pnY * 2 }, temp2, OPH_BACKWARD, OPH_ESTIMATE);
821  fftwShift(temp2, temp3, pnX * 2, pnY * 2, OPH_BACKWARD);
822  //fftExecute(temp3);
823 
824  uint idxOut = 0;
825 
826  for (int idxNy = pnY / 2; idxNy < pnY + (pnY / 2); idxNy++) {
827  for (int idxNx = pnX / 2; idxNx < pnX + (pnX / 2); idxNx++) {
828  out[idxOut] = temp3[idxNy * pnX * 2 + idxNx];
829  idxOut++;
830  }
831  }
832 
833  delete[] in2x;
834  delete[] temp1;
835  delete[] fx;
836  delete[] fy;
837  delete[] prop;
838  delete[] temp2;
839  delete[] temp3;
840 }
841 
842 void ophGen::fresnelPropagation(Complex<Real>* in, Complex<Real>* out, Real distance, uint channel)
843 {
844 #ifdef CHECK_PROC_TIME
845  auto begin = CUR_TIME;
846 #endif
847  const int pnX = context_.pixel_number[_X];
848  const int pnY = context_.pixel_number[_Y];
849  const Real ppX = context_.pixel_pitch[_X];
850  const Real ppY = context_.pixel_pitch[_Y];
851  const uint pnXY = pnX * pnY;
852  const Real lambda = context_.wave_length[channel];
853 
854  Complex<Real>* in2x = new Complex<Real>[pnXY * 4];
855  Complex<Real> zero(0, 0);
856  memsetArr<Complex<Real>>(in2x, zero, 0, pnXY * 4 - 1);
857 
858  uint idxIn = 0;
859  int idxnY = pnY / 2;
860 
861 #ifdef _OPENMP
862 #pragma omp parallel
863  {
864 #pragma omp parallel for private(idxnY) reduction(+:idxIn)
865 #endif
866  for (idxnY = pnY / 2; idxnY < pnY + (pnY / 2); idxnY++) {
867  for (int idxnX = pnX / 2; idxnX < pnX + (pnX / 2); idxnX++) {
868  in2x[idxnY * pnX * 2 + idxnX] = in[idxIn++];
869  }
870  }
871 #ifdef _OPENMP
872  }
873 #endif
874 
875  Complex<Real>* temp1 = new Complex<Real>[pnXY * 4];
876 
877  fft2({ pnX * 2, pnY * 2 }, in2x, OPH_FORWARD, OPH_ESTIMATE);
878  fftwShift(in2x, temp1, pnX*2, pnY*2, OPH_FORWARD, false);
879 
880  Real* fx = new Real[pnXY * 4];
881  Real* fy = new Real[pnXY * 4];
882 
883  uint i = 0;
884  for (int idxFy = -pnY; idxFy < pnY; idxFy++) {
885  for (int idxFx = -pnX; idxFx < pnX; idxFx++) {
886  fx[i] = idxFx / (2 * pnX * ppX);
887  fy[i] = idxFy / (2 * pnY * ppY);
888  i++;
889  }
890  }
891 
892  Complex<Real>* prop = new Complex<Real>[pnXY * 4];
893  memsetArr<Complex<Real>>(prop, zero, 0, pnXY * 4 - 1);
894 
895  Real sqrtPart;
896 
897  Complex<Real>* temp2 = new Complex<Real>[pnXY * 4];
898 
899  for (int i = 0; i < pnXY * 4; i++) {
900  sqrtPart = sqrt(1 / (lambda * lambda) - fx[i] * fx[i] - fy[i] * fy[i]);
901  prop[i][_IM] = 2 * M_PI * distance;
902  prop[i][_IM] *= sqrtPart;
903  temp2[i] = temp1[i] * exp(prop[i]);
904  }
905 
906  Complex<Real>* temp3 = new Complex<Real>[pnXY * 4];
907  fft2({ pnX * 2, pnY * 2 }, temp2, OPH_BACKWARD, OPH_ESTIMATE);
908  fftwShift(temp2, temp3, pnX * 2, pnY * 2, OPH_BACKWARD, false);
909 
910  uint idxOut = 0;
911  // 540 ~ 1620
912  // 960 ~ 2880
913  // 540 * 1920 * 2 + 960
914  for (int idxnY = pnY / 2; idxnY < pnY + (pnY / 2); idxnY++) {
915  for (int idxnX = pnX / 2; idxnX < pnX + (pnX / 2); idxnX++) {
916  out[idxOut++] = temp3[idxnY * pnX * 2 + idxnX];
917  }
918  }
919  //delete[] in;
920  delete[] in2x;
921  delete[] temp1;
922  delete[] fx;
923  delete[] fy;
924  delete[] prop;
925  delete[] temp2;
926  delete[] temp3;
927 #ifdef CHECK_PROC_TIME
928  auto end = CUR_TIME;
929  LOG("\n%s : %lf(s)\n\n",
930  __FUNCTION__,
931  ((chrono::duration<Real>)(end - begin)).count()
932  );
933 #endif
934 }
935 
936 void ophGen::waveCarry(Real carryingAngleX, Real carryingAngleY, Real distance)
937 {
938  const uint pnX = context_.pixel_number[_X];
939  const uint pnY = context_.pixel_number[_Y];
940  const uint pnXY = pnX * pnY;
941  const uint nChannel = context_.waveNum;
942  Real dfx = 1 / context_.pixel_pitch[_X] / pnX;
943  Real dfy = 1 / context_.pixel_pitch[_Y] / pnY;
944  Real* fx = new Real[pnXY];
945  Real* fy = new Real[pnXY];
946  Real* fz = new Real[pnXY];
947  uint i = 0;
948  for (uint ch = 0; ch < nChannel; ch++) {
949  Real lambda = context_.wave_length[ch];
950 
951  for (int idxFy = pnY / 2; idxFy > -pnY / 2; idxFy--) {
952  for (int idxFx = -pnX / 2; idxFx < pnX / 2; idxFx++) {
953  fx[i] = idxFx * dfx;
954  fy[i] = idxFy * dfy;
955  fz[i] = sqrt((1 / lambda)*(1 / lambda) - fx[i] * fx[i] - fy[i] * fy[i]);
956 
957  i++;
958  }
959  }
960 
961  Complex<Real>* carrier = new Complex<Real>[pnXY];
962 
963  for (int i = 0; i < pnXY; i++) {
964  carrier[i][_RE] = 0;
965  carrier[i][_IM] = distance * tan(carryingAngleX)*fx[i] + distance * tan(carryingAngleY)*fy[i];
966  complex_H[ch][i] = complex_H[ch][i] * exp(carrier[i]);
967  }
968 
969  delete[] carrier;
970  }
971  delete[] fx;
972  delete[] fy;
973  delete[] fz;
974 }
975 
976 void ophGen::encodeSideBand(bool bCPU, ivec2 sig_location)
977 {
978  if (complex_H == nullptr) {
979  LOG("Not found diffracted data.");
980  return;
981  }
982 
983  const uint pnX = context_.pixel_number[_X];
984  const uint pnY = context_.pixel_number[_Y];
985 
986  int cropx1, cropx2, cropx, cropy1, cropy2, cropy;
987  if (sig_location[1] == 0) { //Left or right half
988  cropy1 = 1;
989  cropy2 = pnY;
990  }
991  else {
992  cropy = (int)floor(((Real)pnY) / 2);
993  cropy1 = cropy - (int)floor(((Real)cropy) / 2);
994  cropy2 = cropy1 + cropy - 1;
995  }
996 
997  if (sig_location[0] == 0) { // Upper or lower half
998  cropx1 = 1;
999  cropx2 = pnX;
1000  }
1001  else {
1002  cropx = (int)floor(((Real)pnX) / 2);
1003  cropx1 = cropx - (int)floor(((Real)cropx) / 2);
1004  cropx2 = cropx1 + cropx - 1;
1005  }
1006 
1007  cropx1 -= 1;
1008  cropx2 -= 1;
1009  cropy1 -= 1;
1010  cropy2 -= 1;
1011 
1012  if (bCPU)
1013  encodeSideBand_CPU(cropx1, cropx2, cropy1, cropy2, sig_location);
1014  else
1015  encodeSideBand_GPU(cropx1, cropx2, cropy1, cropy2, sig_location);
1016 }
1017 
1018 void ophGen::encodeSideBand_CPU(int cropx1, int cropx2, int cropy1, int cropy2, oph::ivec2 sig_location)
1019 {
1020  const uint pnX = context_.pixel_number[_X];
1021  const uint pnY = context_.pixel_number[_Y];
1022  const uint pnXY = pnX * pnY;
1023  const uint nChannel = context_.waveNum;
1024 
1025  Complex<Real>* h_crop = new Complex<Real>[pnXY];
1026 
1027  for (uint ch = 0; ch < nChannel; ch++) {
1028 
1029  memset(h_crop, 0.0, sizeof(Complex<Real>) * pnXY);
1030 
1031  int p = 0;
1032 #pragma omp parallel for private(p)
1033  for (p = 0; p < pnXY; p++)
1034  {
1035  int x = p % pnX;
1036  int y = p / pnX;
1037  if (x >= cropx1 && x <= cropx2 && y >= cropy1 && y <= cropy2)
1038  h_crop[p] = complex_H[ch][p];
1039  }
1040 
1041  Complex<Real> *in = nullptr;
1042 
1043  fft2(ivec2(pnX, pnY), in, OPH_BACKWARD);
1044  fftwShift(h_crop, h_crop, pnX, pnY, OPH_BACKWARD, true);
1045 
1046  memset(holo_encoded[ch], 0.0, sizeof(Real) * pnXY);
1047  int i = 0;
1048 #pragma omp parallel for private(i)
1049  for (i = 0; i < pnXY; i++) {
1050  Complex<Real> shift_phase(1, 0);
1051  getShiftPhaseValue(shift_phase, i, sig_location);
1052 
1053  holo_encoded[ch][i] = (h_crop[i] * shift_phase).real();
1054  }
1055  }
1056  delete[] h_crop;
1057 }
1058 
1059 extern "C"
1060 {
1073  void cudaFFT(CUstream_st* stream, int nx, int ny, cufftDoubleComplex* in_filed, cufftDoubleComplex* output_field, int direction, bool bNormailized = false);
1074 
1089  void cudaCropFringe(CUstream_st* stream, int nx, int ny, cufftDoubleComplex* in_field, cufftDoubleComplex* out_field, int cropx1, int cropx2, int cropy1, int cropy2);
1090 
1108  void cudaGetFringe(CUstream_st* stream, int pnx, int pny, cufftDoubleComplex* in_field, cufftDoubleComplex* out_field, int sig_locationx, int sig_locationy,
1109  Real ssx, Real ssy, Real ppx, Real ppy, Real PI);
1110 }
1111 
1112 void ophGen::encodeSymmetrization(Complex<Real>* holo, Real* encoded, const ivec2 sig_loc)
1113 {
1114  const int pnX = context_.pixel_number[_X];
1115  const int pnY = context_.pixel_number[_Y];
1116  const int pnXY = pnX * pnY;
1117 
1118  int cropx1, cropx2, cropx, cropy1, cropy2, cropy;
1119  if (sig_loc[1] == 0) //Left or right half
1120  {
1121  cropy1 = 1;
1122  cropy2 = pnY;
1123 
1124  }
1125  else {
1126 
1127  cropy = floor(pnY / 2);
1128  cropy1 = cropy - floor(cropy / 2);
1129  cropy2 = cropy1 + cropy - 1;
1130  }
1131 
1132  if (sig_loc[0] == 0) // Upper or lower half
1133  {
1134  cropx1 = 1;
1135  cropx2 = pnX;
1136 
1137  }
1138  else {
1139 
1140  cropx = floor(pnX / 2);
1141  cropx1 = cropx - floor(cropx / 2);
1142  cropx2 = cropx1 + cropx - 1;
1143  }
1144 
1145  cropx1 -= 1;
1146  cropx2 -= 1;
1147  cropy1 -= 1;
1148  cropy2 -= 1;
1149 
1150  Complex<Real>* h_crop = new Complex<Real >[pnXY];
1151  memset(h_crop, 0.0, sizeof(Complex<Real>) * pnXY);
1152  int i;
1153 #ifdef _OPENMP
1154 #pragma omp parallel for private(i)
1155 #endif
1156  for (i = 0; i < pnXY; i++) {
1157  int x = i % pnX;
1158  int y = i / pnX;
1159  if (x >= cropx1 && x <= cropx2 && y >= cropy1 && y <= cropy2)
1160  h_crop[i] = holo[i];
1161  }
1162  fftw_complex *in = nullptr, *out = nullptr;
1163  fftw_plan plan = fftw_plan_dft_2d(pnX, pnY, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
1164  fftwShift(h_crop, h_crop, pnX, pnY, -1, true);
1165  fftw_destroy_plan(plan);
1166  fftw_cleanup();
1167 
1168 #ifdef _OPENMP
1169 #pragma omp parallel for private(i)
1170 #endif
1171  for (i = 0; i < pnXY; i++) {
1172  Complex<Real> shift_phase(1, 0);
1173  getShiftPhaseValue(shift_phase, i, sig_loc);
1174  encoded[i] = (h_crop[i] * shift_phase)._Val[_RE];
1175  }
1176 
1177  delete[] h_crop;
1178 }
1179 
1180 void ophGen::encodeSideBand_GPU(int cropx1, int cropx2, int cropy1, int cropy2, oph::ivec2 sig_location)
1181 {
1182  const int pnX = context_.pixel_number[_X];
1183  const int pnY = context_.pixel_number[_Y];
1184  const int pnXY = pnX * pnY;
1185  const Real ppX = context_.pixel_pitch[_X];
1186  const Real ppY = context_.pixel_pitch[_Y];
1187  const Real ssX = context_.ss[_X] = pnX * ppX;
1188  const Real ssY = context_.ss[_Y] = pnY * ppY;
1189  const uint nChannel = context_.waveNum;
1190 
1191  cufftDoubleComplex *k_temp_d_, *u_complex_gpu_;
1192  cudaStream_t stream_;
1193  cudaStreamCreate(&stream_);
1194 
1195  cudaMalloc((void**)&u_complex_gpu_, sizeof(cufftDoubleComplex) * pnXY);
1196  cudaMalloc((void**)&k_temp_d_, sizeof(cufftDoubleComplex) * pnXY);
1197 
1198  for (uint ch = 0; ch < nChannel; ch++) {
1199  cudaMemcpy(u_complex_gpu_, complex_H[ch], sizeof(cufftDoubleComplex) * pnXY, cudaMemcpyHostToDevice);
1200 
1201  cudaMemsetAsync(k_temp_d_, 0, sizeof(cufftDoubleComplex) * pnXY, stream_);
1202  cudaCropFringe(stream_, pnX, pnY, u_complex_gpu_, k_temp_d_, cropx1, cropx2, cropy1, cropy2);
1203 
1204  cudaMemsetAsync(u_complex_gpu_, 0, sizeof(cufftDoubleComplex) * pnXY, stream_);
1205  cudaFFT(stream_, pnX, pnY, k_temp_d_, u_complex_gpu_, 1, true);
1206 
1207  cudaMemsetAsync(k_temp_d_, 0, sizeof(cufftDoubleComplex) * pnXY, stream_);
1208  cudaGetFringe(stream_, pnX, pnY, u_complex_gpu_, k_temp_d_, sig_location[0], sig_location[1], ssX, ssY, ppX, ppY, M_PI);
1209 
1210  cufftDoubleComplex* sample_fd = (cufftDoubleComplex*)malloc(sizeof(cufftDoubleComplex) * pnXY);
1211  memset(sample_fd, 0.0, sizeof(cufftDoubleComplex) * pnXY);
1212 
1213  cudaMemcpyAsync(sample_fd, k_temp_d_, sizeof(cufftDoubleComplex) * pnXY, cudaMemcpyDeviceToHost), stream_;
1214  memset(holo_encoded[ch], 0.0, sizeof(Real) * pnXY);
1215 
1216  for (int i = 0; i < pnX * pnY; i++)
1217  holo_encoded[ch][i] = sample_fd[i].x;
1218 
1219  cudaFree(sample_fd);
1220  }
1221  cudaStreamDestroy(stream_);
1222 }
1223 
1224 void ophGen::getShiftPhaseValue(oph::Complex<Real>& shift_phase_val, int idx, oph::ivec2 sig_location)
1225 {
1226  const int pnX = context_.pixel_number[_X];
1227  const int pnY = context_.pixel_number[_Y];
1228  const Real ppX = context_.pixel_pitch[_X];
1229  const Real ppY = context_.pixel_pitch[_Y];
1230  const Real ssX = context_.ss[_X] = pnX * ppX;
1231  const Real ssY = context_.ss[_Y] = pnY * ppY;
1232 
1233  if (sig_location[1] != 0)
1234  {
1235  int r = idx / pnX;
1236  int c = idx % pnX;
1237  Real yy = (ssY / 2.0) - (ppY)*r - ppY;
1238 
1239  oph::Complex<Real> val;
1240  if (sig_location[1] == 1)
1241  val[_IM] = 2 * M_PI * (yy / (4 * ppY));
1242  else
1243  val[_IM] = 2 * M_PI * (-yy / (4 * ppY));
1244 
1245  val.exp();
1246  shift_phase_val *= val;
1247  }
1248 
1249  if (sig_location[0] != 0)
1250  {
1251  int r = idx / pnX;
1252  int c = idx % pnX;
1253  Real xx = (-ssX / 2.0) - (ppX)*c - ppX;
1254 
1255  oph::Complex<Real> val;
1256  if (sig_location[0] == -1)
1257  val[_IM] = 2 * M_PI * (-xx / (4 * ppX));
1258  else
1259  val[_IM] = 2 * M_PI * (xx / (4 * ppX));
1260 
1261  val.exp();
1262  shift_phase_val *= val;
1263  }
1264 }
1265 
1266 void ophGen::getRandPhaseValue(Complex<Real>& rand_phase_val, bool rand_phase)
1267 {
1268  if (rand_phase)
1269  {
1270  rand_phase_val[_RE] = 0.0;
1271  Real min, max;
1272 #if REAL_IS_DOUBLE & true
1273  min = 0.0;
1274  max = 1.0;
1275 #else
1276  min = 0.f;
1277  max = 1.f;
1278 #endif
1279  rand_phase_val[_IM] = 2 * M_PI * oph::rand(min, max);
1280  rand_phase_val.exp();
1281 
1282  }
1283  else {
1284  rand_phase_val[_RE] = 1.0;
1285  rand_phase_val[_IM] = 0.0;
1286  }
1287 }
1288 
1289 void ophGen::setResolution(ivec2 resolution)
1290 {
1291  // 기존 해상도와 다르면 버퍼를 다시 생성.
1292  if (context_.pixel_number != resolution) {
1293  setPixelNumber(resolution);
1294  Openholo::setPixelNumberOHC(resolution);
1295  initialize();
1296  }
1297 }
1298 
1300 {
1301  if (holo_encoded) delete[] holo_encoded;
1302  if (holo_normalized) delete[] holo_normalized;
1303 }
#define OPH_BACKWARD
Definition: define.h:67
ENCODE_FLAG
Definition: ophGen.h:72
Real * vertex
Geometry of point clouds.
Definition: ophGen.h:433
Abstract class.
Definition: Openholo.h:83
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:90
Real * color
Color data of point clouds.
Definition: ophGen.h:435
void cudaCropFringe(CUstream_st *stream, int nx, int ny, cufftDoubleComplex *in_field, cufftDoubleComplex *out_field, int cropx1, int cropx2, int cropy1, int cropy2)
Crop input data according to x, y coordinates on GPU.
void setResolution(ivec2 resolution)
Function for setting buffer size.
Definition: ophGen.cpp:1289
void encodeSideBand_CPU(int cropx1, int cropx2, int cropy1, int cropy2, ivec2 sig_location)
Encode the CGH according to a signal location parameter on the CPU.
Definition: ophGen.cpp:1018
void encodeSymmetrization(Complex< Real > *holo, Real *encoded, const ivec2 sig_loc)
Encoding method.
Definition: ophGen.cpp:1112
Real * wave_length
Definition: Openholo.h:69
ulonglong n_points
Number of points.
Definition: ophGen.h:429
void setPixelNumberOHC(const ivec2 pixel_number)
getter/setter for OHC file read and write
Definition: Openholo.h:311
SSB_PASSBAND
Definition: ophGen.h:227
unsigned char uchar
Definition: typedef.h:64
cufftDoubleComplex * u_complex_gpu_
float Real
Definition: typedef.h:55
void initialize(void)
Initialize variables for Hologram complex field, encoded data, normalized data.
Definition: ophGen.cpp:69
void twoPhaseEncoding(Complex< Real > *holo, Real *encoded, const int size)
Definition: ophGen.cpp:619
#define CUR_TIME
Definition: function.h:58
Real norm(const vec2 &a)
Definition: vec.h:417
void setPixelPitchOHC(const vec2 pixel_pitch)
Definition: Openholo.h:314
vec2 ss
Definition: Openholo.h:67
void cudaFFT(CUstream_st *stream, int nx, int ny, cufftDoubleComplex *in_filed, cufftDoubleComplex *output_field, int direction, bool bNormailized=false)
Convert data from the spatial domain to the frequency domain using 2D FFT on GPU. ...
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:80
#define OPH_ESTIMATE
Definition: define.h:76
void numericalInterference(Complex< Real > *holo, Real *encoded, const int size)
Encoding method.
Definition: ophGen.cpp:593
void getPhase(oph::Complex< Real > *src, Real *dst, const int &size)
Definition: function.h:316
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
void getRandPhaseValue(Complex< Real > &rand_phase_val, bool rand_phase)
Assign random phase value if RANDOM_PHASE == 1.
Definition: ophGen.cpp:1266
void normalize(const Complex< T > *src, Complex< T > *dst, const int &size)
Normalize all elements of Complex<T>* src from 0 to 1.
Definition: function.h:170
#define PI
Definition: ophACPAS.h:9
#define _Y
Definition: define.h:84
#define _IM
Definition: complex.h:57
void setWavelengthOHC(const Real wavelength, const LenUnit wavelength_unit)
Definition: Openholo.h:317
ImgEncoderOhc * OHC_encoder
OHC file format Variables for read and write.
Definition: Openholo.h:304
void encodeSideBand(bool bCPU, ivec2 sig_location)
Encode the CGH according to a signal location parameter.
Definition: ophGen.cpp:976
void normalize(void)
Normalization function to save as image file after hologram creation.
Definition: ophGen.cpp:252
bool save(const char *fname, uint8_t bitsperpixel=8, uchar *src=nullptr, uint px=0, uint py=0)
Function for saving image files.
Definition: ophGen.cpp:258
#define _X
Definition: define.h:80
ivec2 encode_size
Encoded hologram size, varied from encoding type.
Definition: ophGen.h:240
void singleSideBand(Complex< Real > *holo, Real *encoded, const ivec2 holosize, int passband)
Encoding method.
Definition: ophGen.cpp:690
void encodeSideBand_GPU(int cropx1, int cropx2, int cropy1, int cropy2, ivec2 sig_location)
Encode the CGH according to a signal location parameter on the GPU.
Definition: ophGen.cpp:1180
void getShiftPhaseValue(Complex< Real > &shift_phase_val, int idx, ivec2 sig_location)
Calculate the shift phase value.
Definition: ophGen.cpp:1224
void * load(const char *fname)
Function for loading image files.
Definition: ophGen.cpp:333
Real ** holo_encoded
buffer to encoded.
Definition: ophGen.h:248
void burckhardt(Complex< Real > *holo, Real *encoded, const int size)
Definition: ophGen.cpp:646
oph::ivec2 pixel_number
Definition: Openholo.h:63
#define _RE
Definition: complex.h:54
void waveCarry(Real carryingAngleX, Real carryingAngleY, Real distance)
Wave carry.
Definition: ophGen.cpp:936
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:416
void freqShift(Complex< Real > *src, Complex< Real > *dst, const ivec2 holosize, int shift_x, int shift_y)
Frequency shift.
Definition: ophGen.cpp:750
void fresnelPropagation(OphConfig context, Complex< Real > *in, Complex< Real > *out, Real distance)
Fresnel propagation.
Definition: ophGen.cpp:768
void cudaGetFringe(CUstream_st *stream, int pnx, int pny, cufftDoubleComplex *in_field, cufftDoubleComplex *out_field, int sig_locationx, int sig_locationy, Real ssx, Real ssy, Real ppx, Real ppy, Real PI)
Encode the CGH according to a signal location parameter on the GPU.
uint waveNum
Definition: Openholo.h:68
int n_colors
Number of color channel.
Definition: ophGen.h:431
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: Openholo.cpp:249
bool readConfig(const char *fname)
load to configuration file.
Definition: ophGen.cpp:150
bool isPhaseParse
Selects wheter to parse the phase data.
Definition: ophGen.h:439
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2157
virtual uchar * loadAsImg(const char *fname)
Function for loading image files.
Definition: Openholo.cpp:199
cufftDoubleComplex * k_temp_d_
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
void getAmplitude(oph::Complex< Real > *src, Real *dst, const int &size)
Definition: function.h:324
Data for Point Cloud.
Definition: ophGen.h:427
uchar ** holo_normalized
buffer to normalized.
Definition: ophGen.h:250
Real maxOfArr(const std::vector< Real > &arr)
Definition: function.h:86
void encoding()
Definition: ophGen.cpp:525
void resetBuffer()
reset buffer
Definition: ophGen.cpp:363
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
void setPixelNumber(ivec2 n)
Function for setting the output resolution.
Definition: Openholo.h:171
#define FFTW_ESTIMATE
Definition: fftw3.h:392
ophGen(void)
Constructor.
Definition: ophGen.cpp:56
OphConfig context_
Definition: Openholo.h:297
#define FFTW_BACKWARD
Definition: fftw3.h:380
#define OPH_FORWARD
Definition: define.h:66
Complex< Real > ** complex_H
Definition: Openholo.h:298
virtual bool loadAsOhc(const char *fname)
Function to read OHC file.
Definition: ophGen.cpp:345
void propagationAngularSpectrum(int ch, Complex< Real > *input_u, Real propagation_dist, Real k, Real lambda)
Angular spectrum propagation method.
Definition: ophGen.cpp:218
cudaStream_t stream_
int ENCODE_METHOD
Encoding method flag.
Definition: ophGen.h:242
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
unsigned int uint
Definition: typedef.h:62
#define for_i(itr, oper)
Definition: ophGen.cpp:380
#define M_PI
Definition: define.h:52
Real * phase
Phase value of point clouds.
Definition: ophGen.h:437
virtual void ophFree(void)
Pure virtual function for override in child classes.
Definition: ophGen.cpp:1299
virtual ~ophGen(void)=0
Destructor.
Definition: ophGen.cpp:65
oph::vec2 pixel_pitch
Definition: Openholo.h:64