Openholo  v4.2
Open Source Digital Holographic Library
ophACPAS.cpp
Go to the documentation of this file.
1 #define OPH_DM_EXPORT
2 
3 #include "ophACPAS.h"
4 #include <fstream>
5 #include <iostream>
6 #include <string>
7 #include <iomanip>
8 #include "sys.h"
9 #include "tinyxml2.h"
10 #include "PLYparser.h"
11 
12 using namespace std;
13 
15  : ophGen()
16  , m_pHologram(nullptr)
17  , n_points(-1)
18 
19 {
20 }
21 
23 {
24  delete[] m_pHologram;
25 }
26 
27 
28 bool ophACPAS::readConfig(const char* fname)
29 {
30  if (!ophGen::readConfig(fname))
31  return false;
32 
33  bool bRet = true;
34  /*XML parsing*/
35 
36  using namespace tinyxml2;
37  tinyxml2::XMLDocument xml_doc;
38  XMLNode* xml_node = nullptr;
39 
40  if (!checkExtension(fname, ".xml"))
41  {
42  LOG("<FAILED> Wrong file ext.\n");
43  return false;
44  }
45  auto ret = xml_doc.LoadFile(fname);
46  if (ret != XML_SUCCESS)
47  {
48  LOG("<FAILED> Loading file.\n");
49  return false;
50  }
51 
52  xml_node = xml_doc.FirstChild();
53 
54  char szNodeName[32] = { 0, };
55  sprintf(szNodeName, "ScaleX");
56  // about point
57  auto next = xml_node->FirstChildElement(szNodeName);
58  if (!next || XML_SUCCESS != next->QueryDoubleText(&pc_config_.scale[_X]))
59  {
60  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
61  bRet = false;
62  }
63  sprintf(szNodeName, "ScaleY");
64  next = xml_node->FirstChildElement(szNodeName);
65  if (!next || XML_SUCCESS != next->QueryDoubleText(&pc_config_.scale[_Y]))
66  {
67  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
68  bRet = false;
69  }
70  sprintf(szNodeName, "ScaleZ");
71  next = xml_node->FirstChildElement(szNodeName);
72  if (!next || XML_SUCCESS != next->QueryDoubleText(&pc_config_.scale[_Z]))
73  {
74  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
75  bRet = false;
76  }
77  sprintf(szNodeName, "Distance");
78  next = xml_node->FirstChildElement(szNodeName);
79  if (!next || XML_SUCCESS != next->QueryDoubleText(&pc_config_.distance))
80  {
81  LOG("<FAILED> Not found node : \'%s\' (Double) \n", szNodeName);
82  bRet = false;
83  }
84 
85  // ACPAS Environment
86  sprintf(szNodeName, "CGH Width");
87  next = xml_node->FirstChildElement(szNodeName);
88  if (!next || XML_SUCCESS != next->QueryIntText(&env.CghWidth))
89  {
90  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
91  bRet = false;
92  }
93  sprintf(szNodeName, "CGH Height");
94  next = xml_node->FirstChildElement(szNodeName);
95  if (!next || XML_SUCCESS != next->QueryIntText(&env.CghHeight))
96  {
97  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
98  bRet = false;
99  }
100  sprintf(szNodeName, "CGH Scale");
101  next = xml_node->FirstChildElement(szNodeName);
102  if (!next || XML_SUCCESS != next->QueryFloatText(&env.CGHScale))
103  {
104  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
105  bRet = false;
106  }
107  sprintf(szNodeName, "Segmentation Size");
108  next = xml_node->FirstChildElement(szNodeName);
109  if (!next || XML_SUCCESS != next->QueryIntText(&env.SegmentationSize))
110  {
111  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
112  bRet = false;
113  }
114  sprintf(szNodeName, "FFT Segmentation Size");
115  next = xml_node->FirstChildElement(szNodeName);
116  if (!next || XML_SUCCESS != next->QueryIntText(&env.fftSegmentationSize))
117  {
118  LOG("<FAILED> Not found node : \'%s\' (Integer) \n", szNodeName);
119  bRet = false;
120  }
121  sprintf(szNodeName, "Red WaveLength");
122  next = xml_node->FirstChildElement(szNodeName);
123  if (!next || XML_SUCCESS != next->QueryFloatText(&env.rWaveLength))
124  {
125  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
126  bRet = false;
127  }
128 
129  sprintf(szNodeName, "Tilting angle on x axis");
130  next = xml_node->FirstChildElement(szNodeName);
131  if (!next || XML_SUCCESS != next->QueryFloatText(&env.ThetaX))
132  {
133  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
134  bRet = false;
135  }
136  sprintf(szNodeName, "Tilting angle on y axis");
137  next = xml_node->FirstChildElement(szNodeName);
138  if (!next || XML_SUCCESS != next->QueryFloatText(&env.ThetaY))
139  {
140  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
141  bRet = false;
142  }
143  sprintf(szNodeName, "Default depth");
144  next = xml_node->FirstChildElement(szNodeName);
145  if (!next || XML_SUCCESS != next->QueryFloatText(&env.DefaultDepth))
146  {
147  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
148  bRet = false;
149  }
150 
151  sprintf(szNodeName, "3D point interval on x axis");
152  next = xml_node->FirstChildElement(szNodeName);
153  if (!next || XML_SUCCESS != next->QueryFloatText(&env.xInterval))
154  {
155  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
156  bRet = false;
157  }
158  sprintf(szNodeName, "3D point interval on y axis");
159  next = xml_node->FirstChildElement(szNodeName);
160  if (!next || XML_SUCCESS != next->QueryFloatText(&env.yInterval))
161  {
162  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
163  bRet = false;
164  }
165  sprintf(szNodeName, "Hologram interval on xi axis");
166  next = xml_node->FirstChildElement(szNodeName);
167  if (!next || XML_SUCCESS != next->QueryFloatText(&env.xiInterval))
168  {
169  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
170  bRet = false;
171  }
172  sprintf(szNodeName, "Hologram interval on eta axis");
173  next = xml_node->FirstChildElement(szNodeName);
174  if (!next || XML_SUCCESS != next->QueryFloatText(&env.etaInterval))
175  {
176  LOG("<FAILED> Not found node : \'%s\' (Float) \n", szNodeName);
177  bRet = false;
178  }
179 
180  initialize();
181 
182  context_.k = (2 * M_PI) / context_.wave_length[0];
185 
189 
190  LOG("**************************************************\n");
191  LOG(" Read Config (Accurate Phase-Added Stereogram) \n");
192  LOG("1) Focal Length : %.5lf\n", pc_config_.distance);
193  LOG("2) Object Scale : %.5lf / %.5lf / %.5lf\n", pc_config_.scale[_X], pc_config_.scale[_Y], pc_config_.scale[_Z]);
194  LOG("**************************************************\n");
195 
196  return bRet;
197 }
198 
199 int ophACPAS::loadPointCloud(const char* pc_file)
200 {
201  n_points = ophGen::loadPointCloud(pc_file, &pc_data_);
202  return n_points;
203 }
204 
205 int ophACPAS::save(const char * fname, uint8_t bitsperpixel, uchar* src, uint px, uint py)
206 {
207  if (fname == nullptr) return -1;
208 
209  uchar* source = src;
210  ivec2 p(px, py);
211  const uint pnX = context_.pixel_number[_X];
212  const uint pnY = context_.pixel_number[_Y];
213  const uint nChannel = context_.waveNum;
214 
215  for (uint ch = 0; ch < nChannel; ch++) {
216  if (src == nullptr)
217  source = m_lpNormalized[ch];
218  if (px == 0 && py == 0)
219  p = ivec2(pnX, pnY);
220 
221  if (checkExtension(fname, ".bmp")) // when the extension is bmp
222  return Openholo::saveAsImg(fname, bitsperpixel, source, p[_X], p[_Y]) ? 1 : -1;
223  else { // when extension is not .ohf, .bmp - force bmp
224  char buf[256];
225  memset(buf, 0x00, sizeof(char) * 256);
226  sprintf(buf, "%s.bmp", fname);
227 
228  return Openholo::saveAsImg(buf, bitsperpixel, source, p[_X], p[_Y]) ? 1 : -1;
229  }
230  }
231  return -1;
232 }
233 
234 
235 void ophACPAS::DataInit()
236 {
237  if (m_pHologram == nullptr)
239  memset(m_pHologram, 0x00, sizeof(double) * getContext().pixel_number[_X] * getContext().pixel_number[_Y]);
240 
241  for (int i = 0; i < NUMTBL; i++) {
242  float theta = (float)M2_PI * (float)(i + i - 1) / (float)(2 * NUMTBL);
243  m_COStbl[i] = (float)cos(theta);
244  m_SINtbl[i] = (float)sin(theta);
245  }
246 }
247 
248 int ophACPAS::ACPASCalcuation(unsigned char * cghfringe)
249 {
250  double Max = -1E9, Min = 1E9;
251  int cghwidth = getContext().pixel_number[_X];
252  int cghheight = getContext().pixel_number[_Y];
253 
254  DataInit();
255  ACPAS();
256 
257  for (int i = 0; i < cghheight; i++) {
258  for (int j = 0; j < cghwidth; j++) {
259  if (Max < m_pHologram[i*cghwidth + j]) Max = m_pHologram[i*cghwidth + j];
260  if (Min > m_pHologram[i*cghwidth + j]) Min = m_pHologram[i*cghwidth + j];
261  }
262  }
263 
264  for (int i = 0; i < cghheight; i++) {
265  for (int j = 0; j < cghwidth; j++) {
266  double temp = 1.0*(((m_pHologram[i*cghwidth + j] - Min) / (Max - Min))*255. + 0.5);
267  if (temp >= 255.0) cghfringe[i*cghwidth + j] = 255;
268  else cghfringe[i*cghwidth + j] = (unsigned char)(temp);
269  }
270  }
271 
272  return 0;
273 }
274 
275 void ophACPAS::ACPAS()
276 {
277  float xiInterval = env.xiInterval;
278  float etaInterval = env.etaInterval;
279  float rLamda = env.rWaveLength;
280  float rWaveNum = env.rWaveNumber;
281  float thetaX = env.ThetaX;
282  float thetaY = env.ThetaY;
283 
284  int hsegSize = env.SegmentationSize >> 1;
285  int dsegSize = env.SegmentationSize * env.SegmentationSize;
286  int segNumx = (int)(env.CghWidth / env.SegmentationSize);
287  int segNumy = (int)(env.CghHeight/ env.SegmentationSize);
288  int hsegNumx = segNumx >> 1;
289  int hsegNumy = segNumy >> 1;
290 
291  int FFThsegSize = env.fftSegmentationSize >> 1;
292  int FFTdsegSize = env.fftSegmentationSize * env.fftSegmentationSize;
293 
294  float *SFrequency_cx = new float[segNumx]; // Center phase
295  float *SFrequency_cy = new float[segNumy]; // Center phase
296  int *PickPoint_cx = new int[segNumx];
297  int *PickPoint_cy = new int[segNumy];
298  int *Coefficient_cx = new int[segNumx];
299  int *Coefficient_cy = new int[segNumy];
300  float *dPhaseSFy = new float[segNumx];
301  float *dPhaseSFx = new float[segNumy];
302  float *xc = new float[segNumx];
303  float *yc = new float[segNumy];
304  float sf_base = 1.0 / (env.xiInterval * env.fftSegmentationSize);
305  float theta;
306  float theta_s, theta_c;
307  int dtheta_s, dtheta_c;
308  int idx_c, idx_s;
309  double *Compensation_cx = new double[segNumx];
310  double *Compensation_cy = new double[segNumy];
311 
312  //CString mm;
313 
314  fftw_complex *in, *out;
315  fftw_plan plan;
316 
317  double **inRe = new double *[segNumy * segNumx];
318  double **inIm = new double *[segNumy * segNumx];
319  for (int i = 0; i < segNumy; i++) {
320  for (int j = 0; j < segNumx; j++) {
321  inRe[i * segNumx + j] = new double[FFTdsegSize];
322  inIm[i * segNumx + j] = new double[FFTdsegSize];
323  memset(inRe[i*segNumx + j], 0x00, sizeof(double) * FFTdsegSize);
324  memset(inIm[i*segNumx + j], 0x00, sizeof(double) * FFTdsegSize);
325  }
326  }
327 
328  in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTdsegSize);
329  out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * FFTdsegSize);
330  memset(in, 0x00, sizeof(fftw_complex) * FFTdsegSize);
331  memset(m_pHologram, 0x00, sizeof(double) * env.CghWidth * env.CghHeight);
332 
333  for (int i = 0; i < segNumy; i++)
334  yc[i] = ((i - hsegNumy) * env.SegmentationSize + hsegSize) * env.etaInterval;
335  for (int i = 0; i < segNumx; i++)
336  xc[i] = (((i- hsegNumx) * env.SegmentationSize) + hsegSize) * env.xiInterval;
337 
338  for (int i = 0; i < pc_data_.n_points; i++)
339  {
340  float x = pc_data_.vertices[i].point.pos[_X] * pc_config_.scale[_X];
341  float y = pc_data_.vertices[i].point.pos[_Y] * pc_config_.scale[_Y];
342  float z = pc_data_.vertices[i].point.pos[_Z] * pc_config_.scale[_Z] - pc_config_.distance;
343  float amplitude = pc_data_.vertices[i].color.color[_R]; // 1 channel
344  float phase = pc_data_.vertices[i].phase;
345 
346  for (int j = 0; j < segNumy; j++)
347  {
348  float theta_cy = (yc[j] - y) - z;
349  SFrequency_cy[j] = (theta_cy + thetaY) / env.rWaveLength;
350  PickPoint_cy[j] = (SFrequency_cy[j] >= 0 ) ? (int)(SFrequency_cy[j] / sf_base + 0.5) : (int)(SFrequency_cy[j] / sf_base - 0.5);
351  Coefficient_cy[j] = (abs(PickPoint_cy[j]) < FFThsegSize) ? ((env.fftSegmentationSize - PickPoint_cy[j]) % env.fftSegmentationSize) : 0;
352  Compensation_cy[j] = (float)(2 * PI * ((yc[j] - y) * SFrequency_cy[j] + PickPoint_cy[j] * sf_base * FFThsegSize * xiInterval));
353  }
354  for (int j = 0; j < segNumx; j++)
355  {
356  float theta_cx = (xc[j] - x) / z;
357  SFrequency_cx[j] = (theta_cx + thetaX) / env.rWaveLength;
358  PickPoint_cx[j] = (SFrequency_cx[j] >= 0) ? (int)(SFrequency_cx[j] / sf_base + 0.5) : (int)(SFrequency_cx[j] / sf_base - 0.5);
359  Coefficient_cx[j] = (abs(PickPoint_cx[j]) < FFThsegSize) ? ((env.fftSegmentationSize - PickPoint_cx[j]) % env.fftSegmentationSize) : 0;
360  Compensation_cx[j] = (float)(2 * PI * ((xc[j] - x) * SFrequency_cx[j] + PickPoint_cx[j] * sf_base * FFThsegSize * etaInterval));
361  }
362  for (int j = 0; j < segNumy; j++) {
363  for (int k = 0; k < segNumx; k++) {
364  int idx = j* segNumx + k;
365  int idx2 = Coefficient_cy[j] * env.fftSegmentationSize + Coefficient_cx[k];
366  float R = sqrt((xc[k] - x) * (xc[k] - x) + (yc[j] - y) * (yc[j] - y) + z * z);
367  theta = rWaveNum * R + phase + Compensation_cy[j] + Compensation_cx[k];
368  //+ dPhaseSFy[j] + dPhaseSFx[k];
369  theta_c = theta;
370  theta_s = theta + PI;
371  dtheta_c = ((int)(theta_c * NUMTBL / M2_PI));
372  dtheta_s = ((int)(theta_s * NUMTBL / M2_PI));
373  idx_c = (dtheta_c) & (NUMTBL2);
374  idx_s = (dtheta_s) & (NUMTBL2);
375  inRe[idx][idx2] += (double)(amplitude * m_COStbl[idx_c]);
376  inIm[idx][idx2] += (double)(amplitude * m_SINtbl[idx_s]);
377  }
378  }
379  }
380 
381  plan = fftw_plan_dft_2d(env.SegmentationSize, env.SegmentationSize, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
382 
383  for (int j = 0; j < segNumy; j++) {
384  for (int k = 0; k < segNumx; k++) {
385  int idx = j * segNumx + k;
386  memset(in, 0x00, sizeof(fftw_complex) * FFTdsegSize);
387  for (int l = 0; l < env.fftSegmentationSize; l++) {
388  for (int m = 0; m < env.fftSegmentationSize; m++) {
389  int idx2 = l * env.fftSegmentationSize + m;
390  in[l * env.fftSegmentationSize + m][0] = inRe[idx][idx2];
391  in[l * env.fftSegmentationSize + m][1] = inIm[idx][idx2];
392  }
393  }
394  fftw_execute(plan);
395  for (int l = 0; l < env.SegmentationSize; l++) {
396  for (int m = 0; m < env.SegmentationSize; m++) {
397  m_pHologram[(j * env.SegmentationSize + l) * env.CghWidth + (k * env.SegmentationSize + m)] +=
398  out[(l + FFThsegSize - hsegSize) * env.fftSegmentationSize + (m + FFThsegSize - hsegSize)][0];// - out[l * SEGSIZE + m][1];
399  }
400  }
401  }
402  }
403 
404  fftw_destroy_plan(plan);
405  fftw_free(in);
406  fftw_free(out);
407  delete[] dPhaseSFy;
408  delete[] dPhaseSFx;
409  delete[] SFrequency_cx;
410  delete[] SFrequency_cy;
411  delete[] PickPoint_cx;
412  delete[] PickPoint_cy;
413  delete[] Coefficient_cx;
414  delete[] Coefficient_cy;
415  delete[] xc;
416  delete[] yc;
417  for (int i = 0; i < segNumy; i++) {
418  for (int j = 0; j < segNumx; j++) {
419  delete[] inRe[i * segNumx + j];
420  delete[] inIm[i * segNumx + j];
421  }
422  }
423  delete[] inRe;
424  delete[] inIm;
425 }
Color color
Definition: struct.h:104
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
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:135
Real k
Definition: Openholo.h:103
#define NUMTBL2
Definition: ophACPAS.h:16
#define _R
Definition: define.h:108
Real distance
Offset value of point cloud.
Definition: ophGen.h:559
int ACPASCalcuation(unsigned char *cghfringe)
Definition: ophACPAS.cpp:248
bool readConfig(const char *fname)
Definition: ophACPAS.cpp:28
int loadPointCloud(const char *pc_file)
Definition: ophACPAS.cpp:199
ulonglong n_points
Number of points.
Definition: ophGen.h:584
void setPixelNumberOHC(const ivec2 pixel_number)
getter/setter for OHC file read and write
Definition: Openholo.h:501
unsigned char uchar
Definition: typedef.h:64
#define NUMTBL
Definition: ophACPAS.h:15
float m_COStbl[NUMTBL]
Definition: ophACPAS.h:72
void initialize(void)
Initialize variables for Hologram complex field, encoded data, normalized data.
Definition: ophGen.cpp:145
float m_SINtbl[NUMTBL]
Definition: ophACPAS.h:73
void setPixelPitchOHC(const vec2 pixel_pitch)
Definition: Openholo.h:504
vec2 ss
Definition: Openholo.h:104
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
#define M2_PI
Definition: ophACPAS.h:10
structure for 2-dimensional integer vector and its arithmetic.
Definition: ivec.h:66
def k(wvl)
Definition: Depthmap.py:16
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
#define FFTW_BACKWARD
Definition: fftw3.h:380
Real phase
Definition: struct.h:105
ophACPAS()
Definition: ophACPAS.cpp:14
void setWavelengthOHC(const Real wavelength, const LenUnit wavelength_unit)
Definition: Openholo.h:507
vec2 pixel_pitch
Definition: Openholo.h:101
void ACPAS()
Definition: ophACPAS.cpp:275
vec3 scale
Scaling factor of coordinate of point cloud.
Definition: ophGen.h:557
#define _X
Definition: define.h:92
int save(const char *fname, uint8_t bitsperpixel, uchar *src, uint px, uint py)
Definition: ophACPAS.cpp:205
double * m_pHologram
Definition: ophACPAS.h:70
virtual ~ophACPAS()
Definition: ophACPAS.cpp:22
void DataInit()
Definition: ophACPAS.cpp:235
uint waveNum
Definition: Openholo.h:105
ivec2 pixel_number
Definition: Openholo.h:99
bool readConfig(const char *fname)
load to configuration file.
Definition: ophGen.cpp:220
#define _Y
Definition: define.h:96
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
Vertex * vertices
Data of point clouds.
Definition: ophGen.h:588
#define PI
Definition: ophACPAS.h:9
uchar ** m_lpNormalized
buffer to normalized.
Definition: ophGen.h:340
#define _Z
Definition: define.h:100
int fftSegmentationSize
Definition: ophACPAS.h:26
OphConfig context_
Definition: Openholo.h:485
Real * wave_length
Definition: Openholo.h:106
int loadPointCloud(const char *pc_file, OphPointCloudData *pc_data_)
load to point cloud data.
Definition: ophGen.cpp:206
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
unsigned int uint
Definition: typedef.h:62
Definition: ophGen.h:76
#define FFTW_ESTIMATE
Definition: fftw3.h:392
Real pos[3]
Definition: struct.h:87
Real color[3]
Definition: struct.h:95