Openholo  v4.2
Open Source Digital Holographic Library
ophTriMesh.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 "ophTriMesh.h"
47 #include "tinyxml2.h"
48 #include "PLYparser.h"
49 #include "sys.h"
50 
51 #define _X1 0
52 #define _Y1 1
53 #define _Z1 2
54 #define _X2 3
55 #define _Y2 4
56 #define _Z2 5
57 #define _X3 6
58 #define _Y3 7
59 #define _Z3 8
60 
62  : ophGen()
63  , scaledMeshData(nullptr)
64  , no(nullptr)
65  , na(nullptr)
66  , nv(nullptr)
67  , angularSpectrum(nullptr)
68  , refAS(nullptr)
69  , rearAS(nullptr)
70  , convol(nullptr)
71  , phaseTerm(nullptr)
72  , is_ViewingWindow(false)
73  , meshData(nullptr)
74  , tempFreqTermX(nullptr)
75  , tempFreqTermY(nullptr)
76  , streamTriMesh(nullptr)
77  , angularSpectrum_GPU(nullptr)
78  , ffttemp(nullptr)
79 {
80  LOG("*** TRI MESH : BUILD DATE: %s %s ***\n\n", __DATE__, __TIME__);
81 }
82 
83 void ophTri::setViewingWindow(bool is_ViewingWindow)
84 {
85  this->is_ViewingWindow = is_ViewingWindow;
86 }
87 
88 bool ophTri::loadMeshData(const char* fileName, const char* ext)
89 {
90  auto begin = CUR_TIME;
91 
92  if (meshData != nullptr)
93  {
94  delete meshData;
95  }
96 
97  meshData = new OphMeshData;
98 
99  if (!strcmp(ext, "ply")) {
100  PLYparser meshPLY;
101  if (meshPLY.loadPLY(fileName, meshData->n_faces, &meshData->faces))
102  cout << "Mesh Data Load Finished.." << endl;
103  else
104  {
105  LOG("<FAILED> Loading ply file.\n");
106  return false;
107  }
108  }
109  else {
110  LOG("<FAILED> Wrong file ext.\n");
111  return false;
112  }
113 
114  triMeshArray = meshData->faces;
115 
116  LOG("%s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
117  return true;
118 }
119 
120 bool ophTri::readConfig(const char* fname)
121 {
122  if (!ophGen::readConfig(fname))
123  return false;
124 
125  bool bRet = true;
126  using namespace tinyxml2;
127  /*XML parsing*/
128  tinyxml2::XMLDocument xml_doc;
129  XMLNode *xml_node;
130 
131  if (!checkExtension(fname, ".xml"))
132  {
133  LOG("<FAILED> Wrong file ext.\n");
134  return false;
135  }
136  if (xml_doc.LoadFile(fname) != XML_SUCCESS)
137  {
138  LOG("<FAILED> Loading file.\n");
139  return false;
140  }
141  xml_node = xml_doc.FirstChild();
142 
143  char szNodeName[32] = { 0, };
144  sprintf(szNodeName, "ScaleX");
145  // about object
146  auto next = xml_node->FirstChildElement(szNodeName);
147  if (!next || XML_SUCCESS != next->QueryDoubleText(&objSize[_X]))
148  {
149  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
150  bRet = false;
151  }
152  sprintf(szNodeName, "ScaleY");
153  next = xml_node->FirstChildElement(szNodeName);
154  if (!next || XML_SUCCESS != next->QueryDoubleText(&objSize[_Y]))
155  {
156  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
157  bRet = false;
158  }
159  sprintf(szNodeName, "ScaleZ");
160  next = xml_node->FirstChildElement(szNodeName);
161  if (!next || XML_SUCCESS != next->QueryDoubleText(&objSize[_Z]))
162  {
163  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
164  bRet = false;
165  }
166 
167  sprintf(szNodeName, "LampDirectionX");
168  next = xml_node->FirstChildElement(szNodeName);
169  if (!next || XML_SUCCESS != next->QueryDoubleText(&illumination[_X]))
170  {
171  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
172  bRet = false;
173  }
174  sprintf(szNodeName, "LampDirectionY");
175  next = xml_node->FirstChildElement(szNodeName);
176  if (!next || XML_SUCCESS != next->QueryDoubleText(&illumination[_Y]))
177  {
178  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
179  bRet = false;
180  }
181  sprintf(szNodeName, "LampDirectionZ");
182  next = xml_node->FirstChildElement(szNodeName);
183  if (!next || XML_SUCCESS != next->QueryDoubleText(&illumination[_Z]))
184  {
185  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
186  bRet = false;
187  }
188 
189  sprintf(szNodeName, "Random_Phase");
190  // about extra functions
191  next = xml_node->FirstChildElement(szNodeName);
192  if (!next || XML_SUCCESS != next->QueryBoolText(&randPhase))
193  {
194  LOG("<FAILED> Not found node : \'%s\' (Boolean) \n", szNodeName);
195  bRet = false;
196  }
197 
198  sprintf(szNodeName, "Occlusion");
199  next = xml_node->FirstChildElement(szNodeName);
200  if (!next || XML_SUCCESS != next->QueryBoolText(&occlusion))
201  {
202  LOG("<FAILED> Not found node : \'%s\' (Boolean) \n", szNodeName);
203  bRet = false;
204  }
205  sprintf(szNodeName, "Texture");
206  next = xml_node->FirstChildElement(szNodeName);
207  if (!next || XML_SUCCESS != next->QueryBoolText(&textureMapping))
208  {
209  LOG("<FAILED> Not found node : \'%s\' (Boolean) \n", szNodeName);
210  bRet = false;
211  }
212  if (textureMapping == true)
213  {
214  sprintf(szNodeName, "TextureSizeX");
215  next = xml_node->FirstChildElement(szNodeName);
216  if (!next || XML_SUCCESS != next->QueryIntText(&texture.dim[_X]))
217  {
218  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
219  bRet = false;
220  }
221  sprintf(szNodeName, "TextureSizeY");
222  next = xml_node->FirstChildElement(szNodeName);
223  if (!next || XML_SUCCESS != next->QueryIntText(&texture.dim[_Y]))
224  {
225  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
226  bRet = false;
227  }
228  sprintf(szNodeName, "TexturePeriod");
229  next = xml_node->FirstChildElement(szNodeName);
230  if (!next || XML_SUCCESS != next->QueryDoubleText(&texture.period))
231  {
232  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
233  bRet = false;
234  }
235  }
236  initialize();
237 
238  LOG("**************************************************\n");
239  LOG(" Read Config (Tri Mesh) \n");
240  LOG("1) Illumination Direction : %.5lf / %.5lf / %.5lf\n", illumination[_X], illumination[_Y], illumination[_Z]);
241  LOG("2) Object Scale : %.5lf / %.5lf / %.5lf\n", objSize[_X], objSize[_Y], objSize[_Z]);
242  LOG("**************************************************\n");
243  return bRet;
244 }
245 
246 void ophTri::loadTexturePattern(const char* fileName, const char* ext)
247 {
249 
250  if (tempFreqTermX != nullptr) delete[] tempFreqTermX;
251  if (tempFreqTermY != nullptr) delete[] tempFreqTermY;
252  if (texture.pattern != nullptr) delete[] texture.pattern;
253  if (textFFT != nullptr) delete[] textFFT;
254 
255  // not used.
256  //uchar* image;
257  //image = loadAsImg(fileName);
258 
259  int bytesperpixel;
260  int size[2] = { 0,0 };
261  getImgSize(texture.dim[_X], texture.dim[_Y], bytesperpixel, fileName);
262 
263  texture.freq = 1 / texture.period;
264 
265  texture.pattern = new Complex<Real>[texture.dim[_X] * texture.dim[_Y]];
266  textFFT = new Complex<Real>[texture.dim[_X] * texture.dim[_Y]];
267  fft2(texture.dim, texture.pattern, OPH_FORWARD, OPH_ESTIMATE);
268  fftExecute(texture.pattern);
269  fft2(texture.pattern, textFFT, texture.dim[_X], texture.dim[_Y], OPH_FORWARD);
270 
271 
272  tempFreqTermX = new Real[N];
273  tempFreqTermY = new Real[N];
274 }
275 
276 void ophTri::initializeAS()
277 {
278  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
279  const int N = meshData->n_faces;
280 
281  if (scaledMeshData != nullptr) {
282  delete[] scaledMeshData;
283  }
284 
285  scaledMeshData = new Face[N];
286  memset(scaledMeshData, 0, sizeof(Face) * N);
287 
288  if (angularSpectrum != nullptr) {
289  delete[] angularSpectrum;
290  }
291  angularSpectrum = new Complex<Real>[pnXY];
292  memset(angularSpectrum, 0, sizeof(Complex<Real>) * pnXY);
293 
294  if (rearAS != nullptr) {
295  delete[] rearAS;
296  }
297  rearAS = new Complex<Real>[pnXY];
298  memset(rearAS, 0, sizeof(Complex<Real>) * pnXY);
299 
300  if (refAS != nullptr) {
301  delete[] refAS;
302  }
303  refAS = new Complex<Real>[pnXY];
304  memset(refAS, 0, sizeof(Complex<Real>) * pnXY);
305 
306  if (phaseTerm != nullptr) {
307  delete[] phaseTerm;
308  }
309  phaseTerm = new Complex<Real>[pnXY];
310  memset(phaseTerm, 0, sizeof(Complex<Real>) * pnXY);
311 
312 
313  if (convol != nullptr) {
314  delete[] convol;
315  }
316  convol = new Complex<Real>[pnXY];
317  memset(convol, 0, sizeof(Complex<Real>) * pnXY);
318 
319 
320  if (no != nullptr) {
321  delete[] no;
322  }
323  no = new vec3[N];
324  memset(no, 0, sizeof(vec3) * N);
325 
326 
327  if (na != nullptr) {
328  delete[] na;
329  }
330  na = new vec3[N];
331  memset(na, 0, sizeof(vec3) * N);
332 
333 
334  if (nv != nullptr) {
335  delete[] nv;
336  }
337  nv = new vec3[N * 3];
338  memset(nv, 0, sizeof(vec3) * N * 3);
339 }
340 
341 void ophTri::objSort(bool isAscending)
342 {
343  auto begin = CUR_TIME;
344  int N = meshData->n_faces;
345  Real* centerZ = new Real[N];
346 
347 #ifdef _OPENMP
348 #pragma omp parallel for
349 #endif
350  for (int i = 0; i < N; i++) {
351  centerZ[i] = (
352  scaledMeshData[i].vertices[_FIRST].point.pos[_Z] +
353  scaledMeshData[i].vertices[_SECOND].point.pos[_Z] +
354  scaledMeshData[i].vertices[_THIRD].point.pos[_Z]) / 3;
355  }
356 
357  bool bSwap = false;
358 
359  Face tmp;
360  if (isAscending)
361  {
362  while (true)
363  {
364  for (int i = 0; i < N - 1; i++)
365  {
366  bSwap = false;
367  int j = i + 1;
368  if (centerZ[i] > centerZ[j]) {
369  Real tmpZ = centerZ[i];
370  centerZ[i] = centerZ[j];
371  centerZ[j] = tmpZ;
372 
373  tmp = scaledMeshData[i];
374  scaledMeshData[i] = scaledMeshData[j];
375  scaledMeshData[j] = tmp;
376  bSwap = true;
377  }
378  }
379  if (!bSwap)
380  break;
381  }
382  }
383  else
384  {
385  while (true)
386  {
387  for (int i = 0; i < N - 1; i++)
388  {
389  bSwap = false;
390  int j = i + 1;
391  if (centerZ[i] < centerZ[j]) {
392  Real tmpZ = centerZ[i];
393  centerZ[i] = centerZ[j];
394  centerZ[j] = tmpZ;
395 
396  tmp = scaledMeshData[i];
397  scaledMeshData[i] = scaledMeshData[j];
398  scaledMeshData[j] = tmp;
399  bSwap = true;
400  }
401  }
402  if (!bSwap)
403  break;
404  }
405  }
406  delete[] centerZ;
407  LOG("<END> %s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
408 
409 }
410 
411 vec3 vecCross(const vec3& a, const vec3& b)
412 {
413  vec3 c;
414 
415  c(0) = a(0 + 1) * b(0 + 2) - a(0 + 2) * b(0 + 1);
416 
417  c(1) = a(1 + 1) * b(1 + 2) - a(1 + 2) * b(1 + 1);
418 
419  c(2) = a(2 + 1) * b(2 + 2) - a(2 + 2) * b(2 + 1);
420 
421 
422  return c;
423 }
424 
425 void ophTri::triTimeMultiplexing(char* dirName, uint ENCODE_METHOD, Real cenFx, Real cenFy, Real rangeFx, Real rangeFy, Real stepFx, Real stepFy)
426 {
427  TM = true;
428  char strFxFy[30];
430 
431  int nFx = floor(rangeFx / stepFx);
432  int nFy = floor(rangeFy / stepFy);
433  Real tFx, tFy, tFz;
434  for (int iFy = 0; iFy <= nFy; iFy++) {
435  for (int iFx = 0; iFx <= nFx; iFx++) {
436 
437  tFx = cenFx - rangeFx / 2 + iFx*stepFx;
438  tFy = cenFy - rangeFy / 2 + iFy*stepFy;
439  tFz = sqrt(1.0 / context_.wave_length[0] / context_.wave_length[0] - tFx*tFx - tFy*tFy);
440 
441  carrierWave[_X] = tFx*context_.wave_length[0];
442  carrierWave[_Y] = tFy*context_.wave_length[0];
443  carrierWave[_Z] = tFz*context_.wave_length[0];
444 
446 
449  normalize();
450  sprintf(strFxFy, "%s/holo_%d,%d.bmp", dirName, (int)tFx, (int)tFy);
451  save(strFxFy, 8, nullptr, m_vecEncodeSize[_X], m_vecEncodeSize[_Y]);
452 
455 
458  normalize();
459  sprintf(strFxFy, "%s/AS_%d,%d.bmp", dirName, (int)tFx, (int)tFy);
460  save(strFxFy, 8, nullptr, m_vecEncodeSize[_X], m_vecEncodeSize[_Y]);
461  }
462  }
463 }
464 
466 {
468  LOG("<FAILED> WRONG SHADING_FLAG\n");
469  return 0.0;
470  }
471 
472  resetBuffer();
473 
474  auto begin = CUR_TIME;
475  LOG("**************************************************\n");
476  LOG(" Generate Hologram \n");
477  LOG("1) Algorithm Method : Tri Mesh\n");
478  LOG("2) Generate Hologram with %s\n", m_mode & MODE_GPU ?
479  "GPU" :
480 #ifdef _OPENMP
481  "Multi Core CPU"
482 #else
483  "Single Core CPU"
484 #endif
485  );
486  LOG("3) Random Phase Use : %s\n", GetRandomPhase() ? "Y" : "N");
487  LOG("4) Shading Flag : %s\n", SHADING_FLAG == SHADING_FLAG::SHADING_FLAT ? "Flat" : "Continuous");
488  LOG("5) Number of Mesh : %llu\n", meshData->n_faces);
489  LOG("**************************************************\n");
490 
491  if (m_mode & MODE_GPU)
492  {
493  initialize_GPU();
494  prepareMeshData();
495  objSort(false);
496  generateAS_GPU(SHADING_FLAG);
497  }
498  else
499  {
500  initializeAS();
501  prepareMeshData();
502  objSort(false);
503  generateAS(SHADING_FLAG);
504  }
505 
506  /*
507 
508  if (!(m_mode & MODE_GPU)) {
509  fft2(context_.pixel_number, angularSpectrum, OPH_BACKWARD, OPH_ESTIMATE);
510  fft2(angularSpectrum, *(complex_H), context_.pixel_number[_X], context_.pixel_number[_Y], OPH_BACKWARD);
511  //fft2(context_.pixel_number, *(complex_H), OPH_FORWARD, OPH_ESTIMATE);
512  //fft2(*(complex_H), *(complex_H), context_.pixel_number[_X], context_.pixel_number[_Y], OPH_FORWARD);
513  //fftExecute((*complex_H));
514  //complex_H[0]= angularSpectrum;
515 
516  fftFree();
517  }
518  */
519  //fresnelPropagation(*(complex_H), *(complex_H), objShift[_Z]);
520  Real elapsed_time = ELAPSED_TIME(begin, CUR_TIME);
521  LOG("Total Elapsed Time: %.5lf (s)\n", elapsed_time);
522  m_nProgress = 0;
523  return elapsed_time;
524 }
525 
526 bool ophTri::generateAS(uint SHADING_FLAG)
527 {
528  auto begin = CUR_TIME;
529 
530  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
531  int N = meshData->n_faces;
532 
533  Point* freq = new Point[pnXY];
534  Point* fl = new Point[pnXY];
535  Real *freqTermX = new Real[pnXY];
536  Real *freqTermY = new Real[pnXY];
537 
538  findNormals(SHADING_FLAG);
539 
540  for (uint ch = 0; ch < context_.waveNum; ch++)
541  {
542  calGlobalFrequency(freq, context_.wave_length[ch]);
543 
544  for (int j = 0; j < N; j++)
545  {
546  Face mesh = scaledMeshData[j];
547 
548  geometric geom;// = { 0, };
549  //memcpy((void *)&mesh, &scaledMeshData[9 * j], sizeof(Real) * 9);
550 
551  if (!checkValidity(no[j])) // don't care
552  continue;
553  if (!findGeometricalRelations(mesh, no[j], geom))
554  continue;
555  if (!calFrequencyTerm(freq, fl, freqTermX, freqTermY, geom, context_.wave_length[ch]))
556  continue;
557 
558  switch (SHADING_FLAG)
559  {
560  case SHADING_FLAT:
561  refAS_Flat(no[j], freq, mesh, freqTermX, freqTermY, geom, context_.wave_length[ch]);
562  //refAS_Flat(no[j], freq, mesh, freqTermX, freqTermY, geom);
563  break;
564  case SHADING_CONTINUOUS:
565  refAS_Continuous(j, freqTermX, freqTermY);
566  break;
567  default:
568  LOG("<FAILED> WRONG SHADING_FLAG\n");
569  return false;
570  }
571  if (!refToGlobal(complex_H[ch], freq, fl, geom))
572  continue;
573 
574  m_nProgress = ((Real)(ch * N + j) / (Real)(N * context_.waveNum)) * 100;
575 
576  //m_nProgress = (int)((Real)j * 100 / ((Real)N));
577  }
578  }
579 #if 1
580  delete[] freq;
581  delete[] fl;
582  delete[] scaledMeshData;
583  delete[] freqTermX;
584  delete[] freqTermY;
585  delete[] refAS;
586  delete[] phaseTerm;
587  delete[] convol;
588  delete[] no;
589  delete[] na;
590  delete[] nv;
591  scaledMeshData = nullptr;
592  refAS = nullptr;
593  phaseTerm = nullptr;
594  convol = nullptr;
595 #else
596  for (int i = 0; i < 3; i++) {
597  delete[] freq[i];
598  delete[] fl[i];
599  }
600 
601  delete[] freq, fl, scaledMeshData, freqTermX, freqTermY, refAS, phaseTerm, convol;
602  scaledMeshData = nullptr;
603  refAS = nullptr;
604  phaseTerm = nullptr;
605  convol = nullptr;
606 #endif
607 
608  LOG("<END> %s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
609 
610  return true;
611 }
612 
613 void ophTri::calGlobalFrequency(Point* frequency, Real lambda)
614 {
615  auto begin = CUR_TIME;
616  const int pnX = context_.pixel_number[_X];
617  const int pnY = context_.pixel_number[_Y];
618  const Real ppX = context_.pixel_pitch[_X];
619  const Real ppY = context_.pixel_pitch[_Y];
620 
621  Real dfx = 1 / (ppX * pnX);
622  Real dfy = 1 / (ppY * pnY);
623 
624  int startX = pnX >> 1;
625  int startY = pnY >> 1;
626 
627  Real dfl = 1 / lambda;
628  Real sqdfl = dfl * dfl;
629 
630 #ifdef _OPENMP
631 #pragma omp parallel for firstprivate(pnX, startX, dfy, dfx, sqdfl)
632 #endif
633  for (int i = startY; i > -startY; i--) {
634  Real y = i * dfy;
635  Real yy = y * y;
636 
637  int base = (startY - i) * pnX; // for parallel
638 
639  for (int j = -startX; j < startX; j++) {
640  int idx = base + (j + startX); // for parallel
641  Real x = j * dfx;
642  Real xx = x * x;
643  frequency[idx].pos[_X] = x;
644  frequency[idx].pos[_Y] = y;
645  frequency[idx].pos[_Z] = sqrt(sqdfl - xx - yy);
646  }
647  }
648  LOG("<END> %s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
649 }
650 
651 bool ophTri::findNormals(uint SHADING_FLAG)
652 {
653  auto begin = CUR_TIME;
654 
655  int N = meshData->n_faces;
656 
657  Real normNo = 0.0;
658 
659 #ifdef _OPENMP
660 #pragma omp parallel for reduction(+:normNo)
661 #endif
662  for (int i = 0; i < N; i++) {
663  no[i] = vecCross(
664  {
665  scaledMeshData[i].vertices[_FIRST].point.pos[_X] - scaledMeshData[i].vertices[_SECOND].point.pos[_X],
666  scaledMeshData[i].vertices[_FIRST].point.pos[_Y] - scaledMeshData[i].vertices[_SECOND].point.pos[_Y],
667  scaledMeshData[i].vertices[_FIRST].point.pos[_Z] - scaledMeshData[i].vertices[_SECOND].point.pos[_Z]
668  },
669  {
670  scaledMeshData[i].vertices[_THIRD].point.pos[_X] - scaledMeshData[i].vertices[_SECOND].point.pos[_X],
671  scaledMeshData[i].vertices[_THIRD].point.pos[_Y] - scaledMeshData[i].vertices[_SECOND].point.pos[_Y],
672  scaledMeshData[i].vertices[_THIRD].point.pos[_Z] - scaledMeshData[i].vertices[_SECOND].point.pos[_Z]
673  }
674  );
675  Real tmp = norm(no[i]);
676  normNo += tmp * tmp;
677  }
678  normNo = sqrt(normNo);
679 
680 #ifdef _OPENMP
681 #pragma omp parallel for firstprivate(normNo)
682 #endif
683  for (int i = 0; i < N; i++) {
684  na[i] = no[i] / normNo;
685  }
686 
688  vec3* vertices = new vec3[N * 3];
689  vec3 zeros(0, 0, 0);
690  uint N3 = N * 3;
691  for (uint i = 0; i < N3; i++)
692  {
693  int idxFace = i / 3;
694  int idxVertex = i % 3;
695  std::memcpy(&vertices[i], &scaledMeshData[idxFace].vertices[idxVertex].point, sizeof(Point));
696  }
697 
698  //for (uint idx = 0; idx < N * 3; idx++) {
699  // memcpy(&vertices[idx], &scaledMeshData[idx * 3], sizeof(vec3));
700  //}
701  for (uint idx1 = 0; idx1 < N3; idx1++) {
702  if (vertices[idx1] == zeros)
703  continue;
704  vec3 sum = na[idx1 / 3];
705  uint count = 1;
706  uint* idxes = new uint[N3];
707  idxes[0] = idx1;
708  for (uint idx2 = idx1 + 1; idx2 < N3; idx2++) {
709  if (vertices[idx2] == zeros)
710  continue;
711  if ((vertices[idx1][0] == vertices[idx2][0])
712  && (vertices[idx1][1] == vertices[idx2][1])
713  && (vertices[idx1][2] == vertices[idx2][2])) {
714 
715  sum += na[idx2 / 3];
716  vertices[idx2] = zeros;
717  idxes[count++] = idx2;
718  }
719  }
720  vertices[idx1] = zeros;
721 
722  sum = sum / count;
723  sum = sum / norm(sum);
724  for (uint i = 0; i < count; i++)
725  nv[idxes[i]] = sum;
726 
727  delete[] idxes;
728  }
729 
730  delete[] vertices;
731  }
732 
733  LOG("<END> %s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
734 
735  return true;
736 }
737 
738 bool ophTri::checkValidity(vec3 no)
739 {
740  if (no[_Z] < 0 || (no[_X] == 0 && no[_Y] == 0 && no[_Z] == 0))
741  return false;
742  if (no[_Z] >= 0)
743  return true;
744 
745  return false;
746 }
747 
748 bool ophTri::findGeometricalRelations(Face mesh, vec3 no, geometric& geom)
749 {
750  vec3 n = no / norm(no);
751  Face mesh_local;
752 
753  Real th, ph;
754  if (n[_X] == 0 && n[_Z] == 0)
755  th = 0;
756  else
757  th = atan(n[_X] / n[_Z]);
758 
759  Real temp = n[_Y] / sqrt(n[_X] * n[_X] + n[_Z] * n[_Z]);
760  ph = atan(temp);
761 
762  Real costh = cos(th);
763  Real cosph = cos(ph);
764  Real sinth = sin(th);
765  Real sinph = sin(ph);
766 
767  geom.glRot[0] = costh; geom.glRot[1] = 0; geom.glRot[2] = -sinth;
768  geom.glRot[3] = -sinph * sinth; geom.glRot[4] = cosph; geom.glRot[5] = -sinph * costh;
769  geom.glRot[6] = cosph * sinth; geom.glRot[7] = sinph; geom.glRot[8] = cosph * costh;
770 
771  // get distance 1st pt(x, y, z) between (0, 0, 0)
772  for (int i = 0; i < 3; i++)
773  {
774  int idx = i * 3;
775  geom.glShift[i] = -(
776  geom.glRot[idx + 0] * mesh.vertices[_FIRST].point.pos[_X] +
777  geom.glRot[idx + 1] * mesh.vertices[_FIRST].point.pos[_Y] +
778  geom.glRot[idx + 2] * mesh.vertices[_FIRST].point.pos[_Z]);
779  }
780 
781  mesh_local.vertices[_FIRST].point.pos[_X] = (geom.glRot[0] * mesh.vertices[_FIRST].point.pos[_X] + geom.glRot[1] * mesh.vertices[_FIRST].point.pos[_Y] + geom.glRot[2] * mesh.vertices[_FIRST].point.pos[_Z]) + geom.glShift[_X];
782  mesh_local.vertices[_FIRST].point.pos[_Y] = (geom.glRot[3] * mesh.vertices[_FIRST].point.pos[_X] + geom.glRot[4] * mesh.vertices[_FIRST].point.pos[_Y] + geom.glRot[5] * mesh.vertices[_FIRST].point.pos[_Z]) + geom.glShift[_Y];
783  mesh_local.vertices[_FIRST].point.pos[_Z] = (geom.glRot[6] * mesh.vertices[_FIRST].point.pos[_X] + geom.glRot[7] * mesh.vertices[_FIRST].point.pos[_Y] + geom.glRot[8] * mesh.vertices[_FIRST].point.pos[_Z]) + geom.glShift[_Z];
784 
785  mesh_local.vertices[_SECOND].point.pos[_X] = (geom.glRot[0] * mesh.vertices[_SECOND].point.pos[_X] + geom.glRot[1] * mesh.vertices[_SECOND].point.pos[_Y] + geom.glRot[2] * mesh.vertices[_SECOND].point.pos[_Z]) + geom.glShift[_X];
786  mesh_local.vertices[_SECOND].point.pos[_Y] = (geom.glRot[3] * mesh.vertices[_SECOND].point.pos[_X] + geom.glRot[4] * mesh.vertices[_SECOND].point.pos[_Y] + geom.glRot[5] * mesh.vertices[_SECOND].point.pos[_Z]) + geom.glShift[_Y];
787  mesh_local.vertices[_SECOND].point.pos[_Z] = (geom.glRot[6] * mesh.vertices[_SECOND].point.pos[_X] + geom.glRot[7] * mesh.vertices[_SECOND].point.pos[_Y] + geom.glRot[8] * mesh.vertices[_SECOND].point.pos[_Z]) + geom.glShift[_Z];
788 
789  mesh_local.vertices[_THIRD].point.pos[_X] = (geom.glRot[0] * mesh.vertices[_THIRD].point.pos[_X] + geom.glRot[1] * mesh.vertices[_THIRD].point.pos[_Y] + geom.glRot[2] * mesh.vertices[_THIRD].point.pos[_Z]) + geom.glShift[_X];
790  mesh_local.vertices[_THIRD].point.pos[_Y] = (geom.glRot[3] * mesh.vertices[_THIRD].point.pos[_X] + geom.glRot[4] * mesh.vertices[_THIRD].point.pos[_Y] + geom.glRot[5] * mesh.vertices[_THIRD].point.pos[_Z]) + geom.glShift[_Y];
791  mesh_local.vertices[_THIRD].point.pos[_Z] = (geom.glRot[6] * mesh.vertices[_THIRD].point.pos[_X] + geom.glRot[7] * mesh.vertices[_THIRD].point.pos[_Y] + geom.glRot[8] * mesh.vertices[_THIRD].point.pos[_Z]) + geom.glShift[_Z];
792 
793  Real mul1 = mesh_local.vertices[_THIRD].point.pos[_X] * mesh_local.vertices[_SECOND].point.pos[_Y];
794  Real mul2 = mesh_local.vertices[_THIRD].point.pos[_Y] * mesh_local.vertices[_SECOND].point.pos[_X];
795 
796  if (mul1 == mul2)
797  return false;
798 
799  Real refTri[9] = { 0,0,0,1,1,0,1,0,0 };
800 
801  geom.loRot[0] = (refTri[_X3] * mesh_local.vertices[_SECOND].point.pos[_Y] - refTri[_X2] * mesh_local.vertices[_THIRD].point.pos[_Y]) / (mul1 - mul2);
802  geom.loRot[1] = (refTri[_X3] * mesh_local.vertices[_SECOND].point.pos[_X] - refTri[_X2] * mesh_local.vertices[_THIRD].point.pos[_X]) / (-mul1 + mul2);
803  geom.loRot[2] = (refTri[_Y3] * mesh_local.vertices[_SECOND].point.pos[_Y] - refTri[_Y2] * mesh_local.vertices[_THIRD].point.pos[_Y]) / (mul1 - mul2);
804  geom.loRot[3] = (refTri[_Y3] * mesh_local.vertices[_SECOND].point.pos[_X] - refTri[_Y2] * mesh_local.vertices[_THIRD].point.pos[_X]) / (-mul1 + mul2);
805 
806  if ((geom.loRot[0] * geom.loRot[3] - geom.loRot[1] * geom.loRot[2]) == 0)
807  return false;
808  return true;
809 
810 }
811 
812 bool ophTri::calFrequencyTerm(Point* frequency, Point* fl, Real *freqTermX, Real *freqTermY, geometric& geom, Real lambda)
813 {
814  // p.s. only 1 channel
815  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
816 
817  Real w = 1 / lambda;
818  Real ww = w * w;
819 
820  Real det = geom.loRot[0] * geom.loRot[3] - geom.loRot[1] * geom.loRot[2];
821  Real divDet = 1 / det;
822 
823  Real invLoRot[4];
824  invLoRot[0] = divDet * geom.loRot[3];
825  invLoRot[1] = -divDet * geom.loRot[2];
826  invLoRot[2] = -divDet * geom.loRot[1];
827  invLoRot[3] = divDet * geom.loRot[0];
828 
829  Real carrierWaveLoc[3];
830  Real glRot[9];
831  memcpy(glRot, geom.glRot, sizeof(glRot));
832  memcpy(carrierWaveLoc, carrierWave, sizeof(carrierWaveLoc));
833 
834 #ifdef _OPENMP
835 #pragma omp parallel for firstprivate(w, ww, glRot, carrierWaveLoc, invLoRot)
836 #endif
837  for (int i = 0; i < pnXY; i++)
838  {
839  fl[i].pos[_X] = glRot[0] * frequency[i].pos[_X] + glRot[1] * frequency[i].pos[_Y] + glRot[2] * frequency[i].pos[_Z];
840  fl[i].pos[_Y] = glRot[3] * frequency[i].pos[_X] + glRot[4] * frequency[i].pos[_Y] + glRot[5] * frequency[i].pos[_Z];
841  fl[i].pos[_Z] = sqrt(ww - fl[i].pos[_X] * fl[i].pos[_X] - fl[i].pos[_Y] * fl[i].pos[_Y]);
842 
843  Real flxShifted = fl[i].pos[_X] - w * (glRot[0] * carrierWaveLoc[_X] + glRot[1] * carrierWaveLoc[_Y] + glRot[2] * carrierWaveLoc[_Z]);
844  Real flyShifted = fl[i].pos[_Y] - w * (glRot[3] * carrierWaveLoc[_X] + glRot[4] * carrierWaveLoc[_Y] + glRot[5] * carrierWaveLoc[_Z]);
845 
846  freqTermX[i] = invLoRot[0] * flxShifted + invLoRot[1] * flyShifted;
847  freqTermY[i] = invLoRot[2] * flxShifted + invLoRot[3] * flyShifted;
848  }
849  return true;
850 }
851 
852 bool ophTri::calFrequencyTerm(Real** frequency, Real** fl, Real *freqTermX, Real *freqTermY, geometric& geom)
853 {
854  // p.s. only 1 channel
855  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
856 
857  Real waveLength = context_.wave_length[0];
858  Real w = 1 / waveLength;
859  Real ww = w * w;
860 
861  Real det = geom.loRot[0] * geom.loRot[3] - geom.loRot[1] * geom.loRot[2];
862  Real divDet = 1 / det;
863 
864  Real invLoRot[4];
865  invLoRot[0] = divDet * geom.loRot[3];
866  invLoRot[1] = -divDet * geom.loRot[2];
867  invLoRot[2] = -divDet * geom.loRot[1];
868  invLoRot[3] = divDet * geom.loRot[0];
869 
870  Real carrierWaveLoc[3];
871  Real glRot[9];
872  memcpy(glRot, geom.glRot, sizeof(glRot));
873  memcpy(carrierWaveLoc, carrierWave, sizeof(carrierWaveLoc));
874 
875 #ifdef _OPENMP
876 #pragma omp parallel for firstprivate(w, ww, glRot, carrierWaveLoc, invLoRot)
877 #endif
878  for (int i = 0; i < pnXY; i++) {
879  Real flxShifted;
880  Real flyShifted;
881 
882  fl[_X][i] = glRot[0] * frequency[_X][i] + glRot[1] * frequency[_Y][i] + glRot[2] * frequency[_Z][i];
883  fl[_Y][i] = glRot[3] * frequency[_X][i] + glRot[4] * frequency[_Y][i] + glRot[5] * frequency[_Z][i];
884  fl[_Z][i] = sqrt(ww - fl[_X][i] * fl[_X][i] - fl[_Y][i] * fl[_Y][i]);
885  flxShifted = fl[_X][i] - w * (glRot[0] * carrierWaveLoc[_X] + glRot[1] * carrierWaveLoc[_Y] + glRot[2] * carrierWaveLoc[_Z]);
886  flyShifted = fl[_Y][i] - w * (glRot[3] * carrierWaveLoc[_X] + glRot[4] * carrierWaveLoc[_Y] + glRot[5] * carrierWaveLoc[_Z]);
887 
888  freqTermX[i] = invLoRot[0] * flxShifted + invLoRot[1] * flyShifted;
889  freqTermY[i] = invLoRot[2] * flxShifted + invLoRot[3] * flyShifted;
890  }
891  return true;
892 }
893 
894 void ophTri::refAS_Flat(vec3 no, Point* frequency, Face mesh, Real* freqTermX, Real* freqTermY, geometric& geom, Real lambda)
895 {
898  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
899 
900  vec3 n = no / norm(no);
901 
902  Complex<Real> shadingFactor;
903  Real PI2 = M_PI * 2;
904  Real sqPI2 = PI2 * PI2;
905 
906  Complex<Real> term1(0, 0);
907  term1[_IM] = -PI2 / lambda * (
908  carrierWave[_X] * (geom.glRot[0] * geom.glShift[_X] + geom.glRot[3] * geom.glShift[_Y] + geom.glRot[6] * geom.glShift[_Z])
909  + carrierWave[_Y] * (geom.glRot[1] * geom.glShift[_X] + geom.glRot[4] * geom.glShift[_Y] + geom.glRot[7] * geom.glShift[_Z])
910  + carrierWave[_Z] * (geom.glRot[2] * geom.glShift[_X] + geom.glRot[5] * geom.glShift[_Y] + geom.glRot[8] * geom.glShift[_Z])
911  );
912 
913  if (illumination[_X] == 0 && illumination[_Y] == 0 && illumination[_Z] == 0) {
914  shadingFactor = exp(term1);
915  }
916  else {
917  vec3 normIllu = illumination / norm(illumination);
918  shadingFactor = (2 * (n[_X] * normIllu[_X] + n[_Y] * normIllu[_Y] + n[_Z] * normIllu[_Z]) + 0.3) * exp(term1);
919  if (shadingFactor[_RE] * shadingFactor[_RE] + shadingFactor[_IM] * shadingFactor[_IM] < 0)
920  shadingFactor = 0;
921  }
922  Real dfx = 1 / ssX;
923  Real dfy = 1 / ssY;
924 
925  if (occlusion) {
926  Complex<Real> term1(0, 0);
927  Real dfxy = dfx * dfy;
928 
929 #ifdef _OPENMP
930 #pragma omp parallel for firstprivate(PI2, dfxy, mesh, term1)
931 #endif
932  for (int i = 0; i < pnXY; i++) {
933  term1[_IM] = PI2 * (
934  frequency[i].pos[_X] * mesh.vertices[_FIRST].point.pos[_X] +
935  frequency[i].pos[_Y] * mesh.vertices[_FIRST].point.pos[_Y] +
936  frequency[i].pos[_Z] * mesh.vertices[_FIRST].point.pos[_Z]);
937  rearAS[i] = angularSpectrum[i] * exp(term1) * dfxy;
938  }
939 
940  refASInner_flat(freqTermX, freqTermY); // refAS main function including texture mapping
941 
942  if (randPhase) {
943 #ifdef _OPENMP
944 #pragma omp parallel for firstprivate(PI2, shadingFactor)
945 #endif
946  for (int i = 0; i < pnXY; i++) {
947  Complex<Real> phase(0, 0);
948  phase[_IM] = PI2 * rand(0.0, 1.0, i);
949  convol[i] = shadingFactor * exp(phase) - rearAS[i];
950  }
951  conv_fft2_scale(refAS, convol, refAS, context_.pixel_number);
952  }
953  else {
954  conv_fft2_scale(rearAS, refAS, convol, context_.pixel_number);
955 #ifdef _OPENMP
956 #pragma omp parallel for firstprivate(shadingFactor)
957 #endif
958  for (int i = 0; i < pnXY; i++) {
959  refAS[i] = refAS[i] * shadingFactor - convol[i];
960  }
961  }
962  }
963  else {
964  refASInner_flat(freqTermX, freqTermY); // refAS main function including texture mapping
965 
966  if (randPhase == true) {
967  Complex<Real> phase(0, 0);
968 #ifdef _OPENMP
969 #pragma omp parallel for firstprivate(PI2, shadingFactor)
970 #endif
971  for (int i = 0; i < pnXY; i++) {
972  Real randVal = rand(0.0, 1.0, i);
973  phase[_IM] = PI2 * randVal;
974  phaseTerm[i] = shadingFactor * exp(phase);
975  }
976  conv_fft2_scale(refAS, phaseTerm, refAS, context_.pixel_number);
977  }
978  else {
979 #ifdef _OPENMP
980 #pragma omp parallel for firstprivate(shadingFactor)
981 #endif
982  for (int i = 0; i < pnXY; i++) {
983  refAS[i] *= shadingFactor;
984  }
985  }
986  }
987 }
988 
989 void ophTri::refASInner_flat(Real* freqTermX, Real* freqTermY)
990 {
991  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
992  Real PI2 = M_PI * 2.0;
993 
994 //#ifndef _OPENMP
995 //#pragma omp parallel for private(i)
996 //#endif
997  for (int i = 0; i < pnXY; i++) {
998  if (textureMapping == true) {
999  refAS[i] = 0;
1000  for (int idxFy = -texture.dim[_Y] / 2; idxFy < texture.dim[_Y] / 2; idxFy++) {
1001  for (int idxFx = -texture.dim[_X] / 2; idxFx < texture.dim[_X] / 2; idxFy++) {
1002  textFreqX = idxFx * texture.freq;
1003  textFreqY = idxFy * texture.freq;
1004 
1005  tempFreqTermX[i] = freqTermX[i] - textFreqX;
1006  tempFreqTermY[i] = freqTermY[i] - textFreqY;
1007 
1008  if (tempFreqTermX[i] == -tempFreqTermY[i] && tempFreqTermY[i] != 0.0) {
1009  refTerm1[_IM] = PI2 * tempFreqTermY[i];
1010  refTerm2[_IM] = 1.0;
1011  refTemp = ((Complex<Real>)1.0 - exp(refTerm1)) / (4.0 * M_PI*M_PI*tempFreqTermY[i] * tempFreqTermY[i]) + refTerm2 / (PI2*tempFreqTermY[i]);
1012  refAS[i] = refAS[i] + textFFT[idxFx + texture.dim[_X] / 2 + (idxFy + texture.dim[_Y] / 2)*texture.dim[_X]] * refTemp;
1013 
1014  }
1015  else if (tempFreqTermX[i] == tempFreqTermY[i] && tempFreqTermX[i] == 0.0) {
1016  refTemp = (Real)(1.0 / 2.0);
1017  refAS[i] = refAS[i] + textFFT[idxFx + texture.dim[_X] / 2 + (idxFy + texture.dim[_Y] / 2)*texture.dim[_X]] * refTemp;
1018  }
1019  else if (tempFreqTermX[i] != 0.0 && tempFreqTermY[i] == 0.0) {
1020  refTerm1[_IM] = -PI2 * tempFreqTermX[i];
1021  refTerm2[_IM] = 1.0;
1022  refTemp = (exp(refTerm1) - (Complex<Real>)1.0) / (PI2*tempFreqTermX[i] * PI2*tempFreqTermX[i]) + (refTerm2 * exp(refTerm1)) / (PI2*tempFreqTermX[i]);
1023  refAS[i] = refAS[i] + textFFT[idxFx + texture.dim[_X] / 2 + (idxFy + texture.dim[_Y] / 2)*texture.dim[_X]] * refTemp;
1024  }
1025  else if (tempFreqTermX[i] == 0.0 && tempFreqTermY[i] != 0.0) {
1026  refTerm1[_IM] = PI2 * tempFreqTermY[i];
1027  refTerm2[_IM] = 1.0;
1028  refTemp = ((Complex<Real>)1.0 - exp(refTerm1)) / (4.0 * M_PI*M_PI*tempFreqTermY[i] * tempFreqTermY[i]) - refTerm2 / (PI2*tempFreqTermY[i]);
1029  refAS[i] = refAS[i] + textFFT[idxFx + texture.dim[_X] / 2 + (idxFy + texture.dim[_Y] / 2)*texture.dim[_X]] * refTemp;
1030  }
1031  else {
1032  refTerm1[_IM] = -PI2 * tempFreqTermX[i];
1033  refTerm2[_IM] = -PI2 * (tempFreqTermX[i] + tempFreqTermY[i]);
1034  refTemp = (exp(refTerm1) - (Complex<Real>)1.0) / (4.0 * M_PI*M_PI*tempFreqTermX[i] * tempFreqTermY[i]) + ((Complex<Real>)1.0 - exp(refTerm2)) / (4.0 * M_PI*M_PI*tempFreqTermY[i] * (tempFreqTermX[i] + tempFreqTermY[i]));
1035  refAS[i] = refAS[i] + textFFT[idxFx + texture.dim[_X] / 2 + (idxFy + texture.dim[_Y] / 2)*texture.dim[_X]] * refTemp;
1036  }
1037  }
1038  }
1039  }
1040  else {
1041  if (freqTermX[i] == -freqTermY[i] && freqTermY[i] != 0.0) {
1042  refTerm1[_IM] = PI2 * freqTermY[i];
1043  refTerm2[_IM] = 1.0;
1044  refAS[i] = ((Complex<Real>)1.0 - exp(refTerm1)) / (4.0 * M_PI*M_PI*freqTermY[i] * freqTermY[i]) + refTerm2 / (PI2*freqTermY[i]);
1045  }
1046  else if (freqTermX[i] == freqTermY[i] && freqTermX[i] == 0.0) {
1047  refAS[i] = (Real)(1.0 / 2.0);
1048  }
1049  else if (freqTermX[i] != 0 && freqTermY[i] == 0.0) {
1050  refTerm1[_IM] = -PI2 * freqTermX[i];
1051  refTerm2[_IM] = 1.0;
1052  refAS[i] = (exp(refTerm1) - (Complex<Real>)1.0) / (PI2 * freqTermX[i] * PI2 * freqTermX[i]) + (refTerm2 * exp(refTerm1)) / (PI2*freqTermX[i]);
1053  }
1054  else if (freqTermX[i] == 0 && freqTermY[i] != 0.0) {
1055  refTerm1[_IM] = PI2 * freqTermY[i];
1056  refTerm2[_IM] = 1.0;
1057  refAS[i] = ((Complex<Real>)1.0 - exp(refTerm1)) / (PI2 * PI2 * freqTermY[i] * freqTermY[i]) - refTerm2 / (PI2*freqTermY[i]);
1058  }
1059  else {
1060  refTerm1[_IM] = -PI2 * freqTermX[i];
1061  refTerm2[_IM] = -PI2 * (freqTermX[i] + freqTermY[i]);
1062  refAS[i] = (exp(refTerm1) - (Complex<Real>)1.0) / (PI2 * PI2 * freqTermX[i] * freqTermY[i]) + ((Complex<Real>)1.0 - exp(refTerm2)) / (4.0 * M_PI*M_PI*freqTermY[i] * (freqTermX[i] + freqTermY[i]));
1063  }
1064  }
1065  }
1066 }
1067 
1068 bool ophTri::refAS_Continuous(uint n, Real* freqTermX, Real* freqTermY)
1069 {
1070  const int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
1071 
1072  av = vec3(0.0, 0.0, 0.0);
1073  av[0] = nv[3 * n + 0][0] * illumination[0] + nv[3 * n + 0][1] * illumination[1] + nv[3 * n + 0][2] * illumination[2] + 0.1;
1074  av[2] = nv[3 * n + 1][0] * illumination[0] + nv[3 * n + 1][1] * illumination[1] + nv[3 * n + 1][2] * illumination[2] + 0.1;
1075  av[1] = nv[3 * n + 2][0] * illumination[0] + nv[3 * n + 2][1] * illumination[1] + nv[3 * n + 2][2] * illumination[2] + 0.1;
1076 
1077  Complex<Real> refTerm1(0, 0);
1078  Complex<Real> refTerm2(0, 0);
1079  Complex<Real> refTerm3(0, 0);
1080  Complex<Real> D1, D2, D3;
1081 
1082  for (int i = 0; i < pnXY; i++) {
1083  if (freqTermX[i] == 0.0 && freqTermY[i] == 0.0) {
1084  D1(1.0 / 3.0, 0);
1085  D2(1.0 / 5.0, 0);
1086  D3(1.0 / 2.0, 0);
1087  }
1088 
1089  else if (freqTermX[i] == 0.0 && freqTermY[i] != 0.0) {
1090  refTerm1[_IM] = -2 * M_PI*freqTermY[i];
1091  refTerm2[_IM] = 1;
1092 
1093  D1 = (refTerm1 - 1.0)*refTerm1.exp() / (8.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * freqTermY[i])
1094  - refTerm1 / (4.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * freqTermY[i]);
1095  D2 = -(M_PI*freqTermY[i] + refTerm2) / (4.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * freqTermY[i])*exp(refTerm1)
1096  + refTerm1 / (8.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * freqTermY[i]);
1097  D3 = exp(refTerm1) / (2.0 * M_PI*freqTermY[i]) + (1.0 - refTerm2) / (2.0 * M_PI*freqTermY[i]);
1098  }
1099  else if (freqTermX[i] != 0.0 && freqTermY[i] == 0.0) {
1100  refTerm1[_IM] = 4.0 * M_PI*M_PI*freqTermX[i] * freqTermX[i];
1101  refTerm2[_IM] = 1.0;
1102  refTerm3[_IM] = 2.0 * M_PI*freqTermX[i];
1103 
1104  D1 = (refTerm1 + 4.0 * M_PI*freqTermX[i] - 2.0 * refTerm2) / (8.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * freqTermY[i])*exp(-refTerm3)
1105  + refTerm2 / (4.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermX[i]);
1106  D2 = 1.0 / 2.0 * D1;
1107  D3 = ((refTerm3 + 1.0)*exp(-refTerm3) - 1.0) / (4.0 * M_PI*M_PI*freqTermX[i] * freqTermX[i]);
1108  }
1109  else if (freqTermX[i] == -freqTermY[i]) {
1110  refTerm1[_IM] = 1.0;
1111  refTerm2[_IM] = 2.0 * M_PI*freqTermX[i];
1112  refTerm3[_IM] = 2.0 * M_PI*M_PI*freqTermX[i] * freqTermX[i];
1113 
1114  D1 = (-2.0 * M_PI*freqTermX[i] + refTerm1) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermX[i])*exp(-refTerm2)
1115  - (refTerm3 + refTerm1) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermX[i]);
1116  D2 = (-refTerm1) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermX[i])*exp(-refTerm2)
1117  + (-refTerm3 + refTerm1 + 2.0 * M_PI*freqTermX[i]) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermX[i]);
1118  D3 = (-refTerm1) / (4.0 * M_PI*M_PI*freqTermX[i] * freqTermX[i])*exp(-refTerm2)
1119  + (-refTerm2 + 1.0) / (4.0 * M_PI*M_PI*freqTermX[i] * freqTermX[i]);
1120  }
1121  else {
1122  refTerm1[_IM] = -2.0 * M_PI*(freqTermX[i] + freqTermY[i]);
1123  refTerm2[_IM] = 1.0;
1124  refTerm3[_IM] = -2.0 * M_PI*freqTermX[i];
1125 
1126  D1 = exp(refTerm1)*(refTerm2 - 2.0 * M_PI*(freqTermX[i] + freqTermY[i])) / (8 * M_PI*M_PI*M_PI*freqTermY[i] * (freqTermX[i] + freqTermY[i])*(freqTermX[i] + freqTermY[i]))
1127  + exp(refTerm3)*(2.0 * M_PI*freqTermX[i] - refTerm2) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * freqTermY[i])
1128  + ((2.0 * freqTermX[i] + freqTermY[i])*refTerm2) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermX[i] * (freqTermX[i] + freqTermY[i])*(freqTermX[i] + freqTermY[i]));
1129  D2 = exp(refTerm1)*(refTerm2*(freqTermX[i] + 2.0 * freqTermY[i]) - 2.0 * M_PI*freqTermY[i] * (freqTermX[i] + freqTermY[i])) / (8.0 * M_PI*M_PI*M_PI*freqTermY[i] * freqTermY[i] * (freqTermX[i] + freqTermY[i])*(freqTermX[i] + freqTermY[i]))
1130  + exp(refTerm3)*(-refTerm2) / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * freqTermY[i] * freqTermY[i])
1131  + refTerm2 / (8.0 * M_PI*M_PI*M_PI*freqTermX[i] * (freqTermX[i] + freqTermY[i])* (freqTermX[i] + freqTermY[i]));
1132  D3 = -exp(refTerm1) / (4.0 * M_PI*M_PI*freqTermY[i] * (freqTermX[i] + freqTermY[i]))
1133  + exp(refTerm3) / (4.0 * M_PI*M_PI*freqTermX[i] * freqTermY[i])
1134  - 1.0 / (4.0 * M_PI*M_PI*freqTermX[i] * (freqTermX[i] + freqTermY[i]));
1135  }
1136  refAS[i] = (av[1] - av[0])*D1 + (av[2] - av[1])*D2 + av[0] * D3;
1137  }
1138  if (randPhase == true) {
1139  Complex<Real> phase(0, 0);
1140  Real PI2 = M_PI * 2.0;
1141  for (int i = 0; i < pnXY; i++) {
1142  Real randVal = rand(0.0, 1.0, i);
1143  phase[_IM] = PI2 * randVal;
1144  phaseTerm[i] = exp(phase);
1145  }
1146 
1147  conv_fft2_scale(refAS, phaseTerm, convol, context_.pixel_number);
1148  }
1149 
1150  return true;
1151 }
1152 
1153 bool ophTri::refToGlobal(Complex<Real> *dst, Point* frequency, Point* fl, geometric& geom)
1154 {
1155  const long long int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
1156 
1157  Real PI2 = M_PI * 2;
1158  Real det = geom.loRot[0] * geom.loRot[3] - geom.loRot[1] * geom.loRot[2];
1159 
1160  if (det == 0)
1161  return false;
1162  if (det < 0)
1163  det = -det;
1164 
1165 
1166  geometric g;
1167  memcpy(&g, &geom, sizeof(geometric));
1168 
1169 #ifdef _OPENMP
1170 #pragma omp parallel for firstprivate(det, g)
1171 #endif
1172  for (int i = 0; i < pnXY; i++)
1173  {
1174  Complex<Real> term1(0, 0);
1175  Complex<Real> term2(0, 0);
1176 
1177  if (frequency[i].pos[_Z] == 0)
1178  term2 = 0;
1179  else
1180  {
1181  term1[_IM] = PI2 * (fl[i].pos[_X] * g.glShift[_X] + fl[i].pos[_Y] * g.glShift[_Y] + fl[i].pos[_Z] * g.glShift[_Z]);
1182  term2 = refAS[i] / det * fl[i].pos[_Z] / frequency[i].pos[_Z] * exp(term1);
1183  }
1184  if (abs(term2) > MIN_DOUBLE)
1185  ; // do nothing
1186  else
1187  term2 = 0;
1188 
1189  dst[i] += term2;
1190 
1191  //complex_H[0][i] += term2;
1192  }
1193 
1194  return true;
1195 }
1196 
1197 bool ophTri::refToGlobal(Real** frequency, Real** fl, geometric& geom)
1198 {
1199  const int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
1200 
1201  Real PI2 = M_PI * 2;
1202  Real det = geom.loRot[0] * geom.loRot[3] - geom.loRot[1] * geom.loRot[2];
1203 
1204  if (det == 0)
1205  return false;
1206  if (det < 0)
1207  det = -det;
1208 
1209  geometric g;
1210  memcpy(&g, &geom, sizeof(geometric));
1211  Complex<Real> term1(0, 0);
1212  Complex<Real> term2(0, 0);
1213 
1214 #ifdef _OPENMP
1215 #pragma omp parallel for firstprivate(det, g, term1, term2)
1216 #endif
1217  for (int i = 0; i < pnXY; i++) {
1218  if (frequency[_Z][i] == 0)
1219  term2 = 0;
1220  else {
1221  term1[_IM] = PI2 * (fl[_X][i] * g.glShift[_X] + fl[_Y][i] * g.glShift[_Y] + fl[_Z][i] * g.glShift[_Z]);
1222  term2 = refAS[i] / det * fl[_Z][i] / frequency[_Z][i] * exp(term1);// *phaseTerm[i];
1223  }
1224  if (abs(term2) > MIN_DOUBLE) {}
1225  else { term2 = 0; }
1226  //angularSpectrum[i] += term2;
1227  complex_H[0][i] += term2;
1228  }
1229 
1230  return true;
1231 }
1232 
1233 void ophTri::reconTest(const char* fname)
1234 {
1235  const int pnXY = context_.pixel_number[_X] * context_.pixel_number[_Y];
1236 
1237  Complex<Real>* recon = new Complex<Real>[pnXY];
1238  fresnelPropagation((*complex_H), recon, context_.shift[_Z], 0);
1239  ophGen::encoding(ENCODE_AMPLITUDE, recon, nullptr);
1240 
1241  normalize();
1242  save(fname, 8, nullptr, m_vecEncodeSize[_X], m_vecEncodeSize[_Y]);
1243 
1244  delete[] recon;
1245 }
1246 // correct the output scale of the ophGen::conv_fft2
1247 void ophTri::conv_fft2_scale(Complex<Real>* src1, Complex<Real>* src2, Complex<Real>* dst, ivec2 size)
1248 {
1249  const int N = size[_X] * size[_Y];
1250  if (N <= 0) return;
1251 
1252 
1253  Complex<Real>* src1FT = new Complex<Real>[N];
1254  Complex<Real>* src2FT = new Complex<Real>[N];
1255  Complex<Real>* dstFT = new Complex<Real>[N];
1256 
1257 
1258  fft2(src1, src1FT, size[_X], size[_Y], OPH_FORWARD, (bool)OPH_ESTIMATE);
1259 
1260  fft2(src2, src2FT, size[_X], size[_Y], OPH_FORWARD, (bool)OPH_ESTIMATE);
1261 
1262  Real scale = (Real)N * (Real)N;
1263  for (int i = 0; i < N; i++)
1264  dstFT[i] = src1FT[i] * src2FT[i] * scale;
1265 
1266  fft2(dstFT, dst, size[_X], size[_Y], OPH_BACKWARD, (bool)OPH_ESTIMATE);
1267 
1268  delete[] src1FT;
1269  delete[] src2FT;
1270  delete[] dstFT;
1271 }
1272 
1273 void ophTri::prepareMeshData()
1274 {
1275  auto begin = CUR_TIME;
1276  const int N = meshData->n_faces;
1277  const int N3 = N * 3;
1278 
1279 
1280  Real x_max, x_min, y_max, y_min, z_max, z_min;
1281  x_max = y_max = z_max = MIN_DOUBLE;
1282  x_min = y_min = z_min = MAX_DOUBLE;
1283 
1284  for (int i = 0; i < N3; i++)
1285  {
1286  int idxFace = i / 3;
1287  int idxPoint = i % 3;
1288 
1289  if (x_max < triMeshArray[idxFace].vertices[idxPoint].point.pos[_X]) x_max = triMeshArray[idxFace].vertices[idxPoint].point.pos[_X];
1290  if (x_min > triMeshArray[idxFace].vertices[idxPoint].point.pos[_X]) x_min = triMeshArray[idxFace].vertices[idxPoint].point.pos[_X];
1291  if (y_max < triMeshArray[idxFace].vertices[idxPoint].point.pos[_Y]) y_max = triMeshArray[idxFace].vertices[idxPoint].point.pos[_Y];
1292  if (y_min > triMeshArray[idxFace].vertices[idxPoint].point.pos[_Y]) y_min = triMeshArray[idxFace].vertices[idxPoint].point.pos[_Y];
1293  if (z_max < triMeshArray[idxFace].vertices[idxPoint].point.pos[_Z]) z_max = triMeshArray[idxFace].vertices[idxPoint].point.pos[_Z];
1294  if (z_min > triMeshArray[idxFace].vertices[idxPoint].point.pos[_Z]) z_min = triMeshArray[idxFace].vertices[idxPoint].point.pos[_Z];
1295  }
1296 
1297  Real x_cen = (x_max + x_min) / 2;
1298  Real y_cen = (y_max + y_min) / 2;
1299  Real z_cen = (z_max + z_min) / 2;
1300  vec3 cen(x_cen, y_cen, z_cen);
1301 
1302  Real x_del = x_max - x_min;
1303  Real y_del = y_max - y_min;
1304  Real z_del = z_max - z_min;
1305 
1306  Real del = maxOfArr({ x_del, y_del, z_del });
1307 
1308  vec3 shift = getContext().shift;
1309  vec3 locObjSize = objSize;
1310 #ifdef _OPENMP
1311 #pragma omp parallel for firstprivate(cen, del, locObjSize, shift)
1312 #endif
1313  for (int i = 0; i < N3; i++)
1314  {
1315  int idxFace = i / 3;
1316  int idxPoint = i % 3;
1317  scaledMeshData[idxFace].vertices[idxPoint].point.pos[_X] = (triMeshArray[idxFace].vertices[idxPoint].point.pos[_X] - cen[_X]) / del * locObjSize[_X] + shift[_X];
1318  scaledMeshData[idxFace].vertices[idxPoint].point.pos[_Y] = (triMeshArray[idxFace].vertices[idxPoint].point.pos[_Y] - cen[_Y]) / del * locObjSize[_Y] + shift[_Y];
1319  scaledMeshData[idxFace].vertices[idxPoint].point.pos[_Z] = (triMeshArray[idxFace].vertices[idxPoint].point.pos[_Z] - cen[_Z]) / del * locObjSize[_Z] + shift[_Z];
1320  }
1321 
1322  LOG("<END> %s : %.5lf (sec)\n", __FUNCTION__, ELAPSED_TIME(begin, CUR_TIME));
1323 }
1324 
1325 
1326 void ophTri::encoding(unsigned int ENCODE_FLAG)
1327 {
1328  const uint pnX = context_.pixel_number[_X];
1329  const uint pnY = context_.pixel_number[_Y];
1330  const uint nChannel = context_.waveNum;
1331  Complex<Real>** dst = new Complex<Real>*[nChannel];
1332  for (uint ch = 0; ch < nChannel; ch++) {
1333  dst[ch] = new Complex<Real>[pnX * pnY];
1334  //fft2(context_.pixel_number, nullptr, OPH_BACKWARD);
1335  fft2(complex_H[ch], dst[ch], pnX, pnY, OPH_BACKWARD, true);
1336  ophGen::encoding(ENCODE_FLAG, dst[ch], m_lpEncoded[ch]);
1337  }
1338 
1339  for (uint ch = 0; ch < nChannel; ch++)
1340  delete[] dst[ch];
1341  delete[] dst;
1342 }
bool loadPLY(const std::string &fileName, ulonglong &n_points, Vertex **vertices)
Definition: PLYparser.cpp:142
ENCODE_FLAG
Definition: ophGen.h:84
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:113
OphConfig & getContext(void)
Function for getting the current context.
Definition: Openholo.h:229
#define M_PI
Definition: define.h:52
Point point
Definition: struct.h:103
#define OPH_FORWARD
Definition: define.h:66
w
Definition: Openholo.cpp:429
bytesperpixel
Definition: Openholo.cpp:431
vec3 shift
Definition: Openholo.h:102
#define MIN_DOUBLE
Definition: define.h:148
ulonglong n_faces
The number of faces in object.
Definition: ophGen.h:646
void setViewingWindow(bool is_ViewingWindow)
Set the value of a variable is_ViewingWindow(true or false)
Definition: ophTriMesh.cpp:83
SHADING_FLAG
Mesh object data scaling and shifting.
Definition: ophTriMesh.h:186
#define MODE_GPU
Definition: define.h:156
void initialize(void)
Initialize variables for Hologram complex field, encoded data, normalized data.
Definition: ophGen.cpp:145
Real glShift[3]
Definition: ophTriMesh.h:64
Real norm(const vec2 &a)
Definition: vec.h:417
Definition: struct.h:86
#define MAX_DOUBLE
Definition: define.h:140
vec2 ss
Definition: Openholo.h:104
float Real
Definition: typedef.h:55
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
Real loRot[4]
Definition: ophTriMesh.h:65
structure for 2-dimensional integer vector and its arithmetic.
Definition: ivec.h:66
Vertex vertices[3]
Definition: struct.h:117
#define _IM
Definition: complex.h:58
void loadTexturePattern(const char *fileName, const char *ext)
Definition: ophTriMesh.cpp:246
bool GetRandomPhase()
Function for getting the random phase.
Definition: ophGen.h:537
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
#define _Y2
Definition: ophTriMesh.cpp:55
Real generateHologram(uint SHADING_FLAG)
Hologram generation.
Definition: ophTriMesh.cpp:465
void reconTest(const char *fname)
vec2 pixel_pitch
Definition: Openholo.h:101
#define _SECOND
Definition: define.h:84
void normalize(void)
Normalization function to save as image file after hologram creation.
Definition: ophGen.cpp:677
void setEncodeMethod(unsigned int ENCODE_FLAG)
Definition: ophGen.h:262
void triTimeMultiplexing(char *dirName, uint ENCODE_METHOD, Real cenFx, Real cenFy, Real rangeFx, Real rangeFy, Real stepFx, Real stepFy)
Definition: ophTriMesh.cpp:425
bool save(const char *fname, uint8_t bitsperpixel=8, uchar *src=nullptr, uint px=0, uint py=0)
Function for saving image files.
ivec2 m_vecEncodeSize
Encoded hologram size, varied from encoding type.
Definition: ophGen.h:330
#define _X
Definition: define.h:92
unsigned int m_mode
Definition: ophGen.h:350
bool getImgSize(int &w, int &h, int &bytesperpixel, const char *fname)
Function for getting the image size.
Definition: Openholo.cpp:402
if(infile==nullptr)
Definition: Openholo.cpp:419
Complex< Real > * pattern
Definition: ophTriMesh.h:73
geometrical relations
Definition: ophTriMesh.h:62
Real ** m_lpEncoded
buffer to encoded.
Definition: ophGen.h:338
#define _X2
Definition: ophTriMesh.cpp:54
#define _RE
Definition: complex.h:55
void fftExecute(Complex< Real > *out, bool bReverse=false)
Execution functions to be called after fft1, fft2, and fft3.
Definition: Openholo.cpp:623
void fft2(ivec2 n, Complex< Real > *in, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Functions for performing fftw 2-dimension operations inside Openholo.
Definition: Openholo.cpp:559
#define _THIRD
Definition: define.h:88
bool loadMeshData(const char *fileName, const char *ext)
Mesh data load.
Definition: ophTriMesh.cpp:88
void fresnelPropagation(OphConfig context, Complex< Real > *in, Complex< Real > *out, Real distance)
Fresnel propagation.
Real glRot[9]
Definition: ophTriMesh.h:63
ivec2 dim
Definition: ophTriMesh.h:74
uint waveNum
Definition: Openholo.h:105
Data for triangular mesh.
Definition: ophGen.h:644
structure for 3-dimensional Real type vector and its arithmetic.
Definition: vec.h:466
Real period
Definition: ophTriMesh.h:75
#define _X3
Definition: ophTriMesh.cpp:57
Face * faces
Definition: ophGen.h:647
Complex< T > & exp()
Definition: complex.h:395
ivec2 pixel_number
Definition: Openholo.h:99
bool readConfig(const char *fname)
load to configuration file.
Definition: ophGen.cpp:220
#define _FIRST
Definition: define.h:80
Definition: struct.h:115
#define _Y
Definition: define.h:96
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
Real rand(const Real min, const Real max, oph::ulong _SEED_VALUE=0)
Get random Real value from min to max.
Definition: function.h:294
Complex< Real > ** complex_H
Definition: Openholo.h:489
#define OPH_ESTIMATE
Definition: define.h:76
Complex< Real > * AS
Binary Encoding - Error diffusion.
Definition: ophGen.h:426
bool readConfig(const char *fname)
Triangular mesh basc CGH configuration file load.
Definition: ophTriMesh.cpp:120
Real maxOfArr(const std::vector< Real > &arr)
Definition: function.h:87
void encoding()
Definition: ophGen.cpp:985
void resetBuffer()
reset buffer
Definition: ophGen.cpp:824
#define _Z
Definition: define.h:100
#define ELAPSED_TIME(x, y)
Definition: function.h:59
OphConfig context_
Definition: Openholo.h:485
#define OPH_BACKWARD
Definition: define.h:67
Real freq
Definition: ophTriMesh.h:76
int ENCODE_METHOD
Encoding method flag.
Definition: ophGen.h:332
Real sum(const vec2 &a)
Definition: vec.h:401
Real * wave_length
Definition: Openholo.h:106
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
unsigned int uint
Definition: typedef.h:62
vec3 vecCross(const vec3 &a, const vec3 &b)
Definition: ophTriMesh.cpp:411
#define CUR_TIME
Definition: function.h:58
Definition: ophGen.h:76
Real pos[3]
Definition: struct.h:87
ophTri(void)
Constructor.
Definition: ophTriMesh.cpp:61
#define _Y3
Definition: ophTriMesh.cpp:58