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