Openholo  v4.2
Open Source Digital Holographic Library
ophIFTA.cpp
Go to the documentation of this file.
1 #include "ophIFTA.h"
2 #include "sys.h"
3 #include <omp.h>
4 #include <stdlib.h>
5 #include "include.h"
6 #include "tinyxml2.h"
7 
9 {
10  imgRGB = nullptr;
11  imgDepth = nullptr;
12  imgOutput = nullptr;
13  width = 0;
14  height = 0;
15  bytesperpixel = 0;
16  nearDepth = 0.0;
17  farDepth = 0.0;
18  nDepth = 0;
19  nIteration = 0;
20 }
21 
22 
24 {
25  if (imgOutput) delete[] imgOutput;
26  if (imgRGB) delete[] imgRGB;
27  if (imgDepth) delete[] imgDepth;
28 }
29 
30 using namespace oph;
32 {
33  if ((!imgRGB && m_config.num_of_depth == 1) || (!imgRGB && !imgDepth))
34  return 0.0;
35 
36  srand(time(NULL));
37 
38  auto begin = CUR_TIME;
39 
40  const uint nWave = context_.waveNum;
41  const int pnX = context_.pixel_number[_X];
42  const int pnY = context_.pixel_number[_Y];
43  const Real ppX = context_.pixel_pitch[_X];
44  const Real ppY = context_.pixel_pitch[_Y];
45  const long long int pnXY = pnX * pnY;
46  const Real ssX = context_.ss[_X];
47  const Real ssY = context_.ss[_Y];
48  const Real nDepth = m_config.num_of_depth;
49  Real nearDepth = m_config.near_depthmap;
50  Real farDepth = m_config.far_depthmap;
51  int nIteration = m_config.num_of_iteration;
52  int num_thread = 1;
53 
54  Real *waveRatio = new Real[nWave];
55  for (uint i = 0; i < nWave; i++) {
56  waveRatio[i] = context_.wave_length[nWave - 1] / context_.wave_length[i];
57  }
58 
59  // constants
60  Real d = (nDepth == 1) ? 0.0 : (farDepth - nearDepth) / (nDepth - 1);
61  Real z = farDepth;
62  Real pi2 = 2 * M_PI;
63  uchar m = getMax(imgDepth, pnX, pnY);
64  Real *depth_quant = new Real[pnXY];
65  memset(depth_quant, 0, sizeof(Real) * pnXY);
66 
67  for (int depth = nDepth; depth > 0; depth--) {
68 #ifdef _OPENMP
69 #pragma omp parallel for
70 #endif
71  for (int i = 0; i < pnXY; i++) {
72  depth_quant[i] += (imgDepth[i] > (depth*m / nDepth)) ? 1 : 0;
73  }
74  }
75  if (imgOutput) {
76  delete[] imgOutput;
77  imgOutput = nullptr;
78  }
79  imgOutput = new uchar[pnXY * bytesperpixel];
80  for (int depth = 0; depth < nDepth; depth++) {
81  z = farDepth - (d*depth);
82 
83  for (uint ch = 0; ch < nWave; ch++) {
84  uchar *img = new uchar[pnXY];
85  separateColor(ch, pnX, pnY, imgRGB, img);
86 #ifdef _OPENMP
87 #pragma omp parallel for
88 #endif
89  for (int i = 0; i < pnXY; i++) {
90  Real val = (depth_quant[i] == depth) ? 1.0 : 0.0;
91  img[i] = val * (Real)img[i];
92  }
93  Real lambda = context_.wave_length[ch];
94  Real k = 2 * M_PI / lambda;
95  Real hssX = lambda * z / ppX;
96  Real hssY = lambda * z / ppY;
97  Real hppX = hssX / pnX;
98  Real hppY = hssY / pnY;
99 
100  Real hStartX = -hssX / 2;
101  Real hStartY = -hssY / 2;
102  Real startX = -ssX / 2;
103  Real startY = -ssY / 2;
104 
105  Complex<Real> c1(0.0, lambda * z);
106  Complex<Real> c2(0.0, -k / (2 * lambda));
107  Complex<Real> c3(0.0, -k / (2 * z));
108  uchar *img_tmp = nullptr;
109  uchar *imgScaled = nullptr;
110 
111  // blue¿¡¼­´Â ¸®»çÀÌÁî ÇÒ ÇÊ¿ä°¡ ¾ø´Ù.
112  if (ch < 2)
113  {
114  int scaleX = round(pnX * 4 * waveRatio[ch]);
115  int scaleY = round(pnY * 4 * waveRatio[ch]);
116 
117  int ww = pnX * 4;
118  int hh = pnY * 4;
119  img_tmp = new uchar[ww * hh];
120  imgScaled = new uchar[scaleX * scaleY];
121  imgScaleBilinear(img, imgScaled, pnX, pnY, scaleX, scaleY);
122  delete[] img;
123  Real ppY2 = 2 * ppY;
124  Real ppX2 = 2 * ppX;
125 
126  memset(img_tmp, 0, sizeof(uchar) * ww * hh);
127 
128  int h1 = round((hh - scaleY) / 2);
129  int w1 = round((ww - scaleX) / 2);
130 
131  // À̹ÌÁö¸¦ Áß¾ÓÀ¸·Î Á¶Á¤
132 #ifdef _OPENMP
133 #pragma omp parallel for
134 #endif
135  for (int y = 0; y < scaleY; y++) {
136  for (int x = 0; x < scaleX; x++) {
137  img_tmp[(y + h1)*ww + x + w1] = imgScaled[y*scaleX + x];
138  }
139  }
140  img = new uchar[pnXY];
141  imgScaleBilinear(img_tmp, img, ww, hh, pnX, pnY);
142  }
143  Real *target = new Real[pnXY];
144  Complex<Real> *result1 = new Complex<Real>[pnXY];
145  Complex<Real> *result2 = new Complex<Real>[pnXY];
146  Complex<Real> *kernel = new Complex<Real>[pnXY];
147  Complex<Real> *kernel2 = new Complex<Real>[pnXY];
148 #ifdef _OPENMP
149 #pragma omp parallel for
150 #endif
151  for (int y = 0; y < pnY; y++) {
152  int offset = y * pnX;
153 
154  Real ty = hStartY + (y * hppY);
155  Real yy = startY + (y * ppY);
156 
157  for (int x = 0; x < pnX; x++) {
158  target[offset + x] = (Real)img[offset + x];
159  Real ran;
160 #ifdef ANOTHER_RAND
161  ran = rand(0.0, 1.0);
162 #else
163  ran = ((Real)rand() / RAND_MAX) * 1.0;
164 #endif
165  Complex<Real> tmp, c4;
166  if (ran < 1.0) {
167  tmp(0.0, ran * 2 * M_PI);
168  c4 = tmp.exp();
169  }
170  else
171  c4(1.0, 0.0);
172 
173  Real tx = hStartX + (x * hppX);
174  Real txy = tx * tx + ty * ty;
175  Real xx = startX + (x * ppX);
176  Real xy = xx * xx + yy * yy;
177 
178  Complex<Real> t = (c2 * txy).exp();
179  kernel[offset + x] = c1 * t;
180  kernel2[offset + x] = (c3 * xy).exp();
181  result1[offset + x] = target[offset + x] * c4;
182  result1[offset + x] *= kernel[offset + x];
183  }
184  }
185  Complex<Real> *tmp = new Complex<Real>[pnXY];
186  Complex<Real> *tmp2 = new Complex<Real>[pnXY];
187  memset(tmp, 0, sizeof(Complex<Real>) * pnXY);
188  memset(tmp2, 0, sizeof(Complex<Real>) * pnXY);
189  fftShift(pnX, pnY, result1, tmp);
190  fft2(ivec2(pnX, pnY), tmp, OPH_FORWARD, OPH_ESTIMATE);
191  fftExecute(tmp, true);
192  memset(result1, 0, sizeof(Complex<Real>) * pnXY);
193 
194  if (nIteration == 0) {
195 #ifdef _OPENMP
196 #pragma omp parallel for
197 #endif
198  for (int j = 0; j < pnXY; j++) {
199  result2[j] = tmp[j] / kernel2[j];
200  complex_H[ch][j] += result2[j];
201  }
202  }
203  else {
204 #ifdef _OPENMP
205 #pragma omp parallel for
206 #endif
207  for (int y = 0; y < pnXY; y++) {
208  result2[y] = tmp[y] / kernel2[y];
209  }
210 
211  memset(tmp, 0, sizeof(Complex<Real>) * pnXY);
212 
213  for (int i = 0; i < nIteration; i++) {
214  for (int j = 0; j < pnXY; j++) {
215  result2[j] = tmp[j] / kernel2[j];
216  tmp[j][_RE] = 0.0;
217  tmp[j][_IM] = result2[j].angle();
218  result1[j] = tmp[j].exp() * kernel2[j];
219  }
220 
221  // FFT
222  memset(tmp, 0, sizeof(Complex<Real>) * pnXY);
223  fft2(ivec2(pnX, pnY), result1, OPH_FORWARD, OPH_ESTIMATE);
224  fftExecute(tmp);
225  fftShift(pnX, pnY, tmp, result1);
226 
227  for (int j = 0; j < pnXY; j++) {
228  result1[j] /= kernel[j];
229  Complex<Real> aa(0.0, result1[j].angle());
230  aa = aa.exp();
231  result2[j] = (target[j] / 255.0) * aa;
232  result2[j] *= kernel[j];
233  }
234  fftShift(pnX, pnY, result2, tmp);
235  fft2(ivec2(pnX, pnY), tmp, OPH_FORWARD, OPH_ESTIMATE);
236  fftExecute(tmp2, true);
237 
238  for (int j = 0; j < pnXY; j++) {
239  result2[j] = tmp2[j] / kernel2[j];
240  }
241  LOG("Iteration (%d / %d)\n", i + 1, nIteration);
242  }
243  memset(img, 0, pnXY);
244  for (int j = 0; j < pnXY; j++) {
245  complex_H[ch][j] += result2[j];
246  }
247  }
248  m_nProgress = (int)((Real)(depth*nWave + ch) * 100 / ((Real)nDepth * nWave));
249 
250  delete[] kernel;
251  delete[] kernel2;
252  delete[] tmp;
253  delete[] tmp2;
254  delete[] result2;
255  delete[] result1;
256  delete[] target;
257  delete[] img_tmp;
258  delete[] imgScaled;
259  delete[] img;
260 
261  LOG("Color Channel (%d / %d) %lf(s)\n", ch + 1, nWave, ELAPSED_TIME(begin, CUR_TIME));
262  }
263  LOG("Depth Level (%d / %d) %lf(s)\n", depth + 1, (int)nDepth, ELAPSED_TIME(begin, CUR_TIME));
264  }
265  delete[] depth_quant;
266  delete[] waveRatio;
267  auto end = CUR_TIME;
268  LOG("\nTotal Elapsed Time : %lf(s)\n\n", ELAPSED_TIME(begin, end));
269  return ELAPSED_TIME(begin, end);
270 }
271 
272 bool ophIFTA::normalize()
273 {
274  const uint nWave = context_.waveNum;
275  const uint pnX = context_.pixel_number[_X];
276  const uint pnY = context_.pixel_number[_Y];
277  const long long int pnXY = pnX * pnY;
278  const Real pi2 = M_PI * 2;
279 
280  for (uint ch = 0; ch < nWave; ch++) {
281  for (int i = 0; i < pnXY; i++) {
282  m_lpNormalized[ch][i] = m_lpEncoded[ch][i] / pi2 * 255;
283  }
284  }
285  return true;
286 }
287 
288 void ophIFTA::encoding(unsigned int ENCODE_FLAG, unsigned int SSB_PASSBAND)
289 {
290  auto begin = CUR_TIME;
291  const uint pnX = context_.pixel_number[_X];
292  const uint pnY = context_.pixel_number[_Y];
293  const uint nChannel = context_.waveNum;
294  Complex<Real>* dst = new Complex<Real>[pnX * pnY];
295 
296  for (uint ch = 0; ch < nChannel; ch++) {
297  fft2(context_.pixel_number, complex_H[ch], OPH_BACKWARD);
298  fft2(complex_H[ch], dst, pnX, pnY, OPH_BACKWARD);
299 
300  if (ENCODE_FLAG == ophGen::ENCODE_SSB) {
301  ivec2 location;
302  switch (SSB_PASSBAND) {
303  case SSB_TOP:
304  location = ivec2(0, 1);
305  break;
306  case SSB_BOTTOM:
307  location = ivec2(0, -1);
308  break;
309  case SSB_LEFT:
310  location = ivec2(-1, 0);
311  break;
312  case SSB_RIGHT:
313  location = ivec2(1, 0);
314  break;
315  }
316 
317  encodeSideBand(true, location);
318  }
319  else ophGen::encoding(ENCODE_FLAG, SSB_PASSBAND, dst);
320  }
321  delete[] dst;
322  auto end = CUR_TIME;
323  LOG("Elapsed Time: %lf(s)\n", ELAPSED_TIME(begin, end));
324 }
325 
326 bool ophIFTA::readConfig(const char* fname)
327 {
328  if (!ophGen::readConfig(fname))
329  return false;
330 
331  using namespace tinyxml2;
332  /*XML parsing*/
333  tinyxml2::XMLDocument xml_doc;
334  XMLNode *xml_node;
335 
336  if (!checkExtension(fname, ".xml"))
337  {
338  LOG("file's extension is not 'xml'\n");
339  return false;
340  }
341  if (xml_doc.LoadFile(fname) != XML_SUCCESS)
342  {
343  LOG("Failed to load file \"%s\"\n", fname);
344  return false;
345  }
346 
347  xml_node = xml_doc.FirstChild();
348 
349  // about viewing window
350  auto next = xml_node->FirstChildElement("DepthLevel");
351  if (!next || XML_SUCCESS != next->QueryIntText(&m_config.num_of_depth))
352  m_config.num_of_depth = 1;
353  next = xml_node->FirstChildElement("NearOfDepth");
354  if (!next || XML_SUCCESS != next->QueryDoubleText(&m_config.near_depthmap))
355  return false;
356  next = xml_node->FirstChildElement("FarOfDepth");
357  if (!next || XML_SUCCESS != next->QueryDoubleText(&m_config.far_depthmap))
358  return false;
359  next = xml_node->FirstChildElement("NumOfIteration");
360  if (!next || XML_SUCCESS != next->QueryIntText(&m_config.num_of_iteration))
361  m_config.num_of_iteration = 1;
362 
363  initialize();
364  return true;
365 }
366 
367 
368 bool ophIFTA::readImage(const char* fname, bool bRGB)
369 {
370  bool ret = getImgSize(width, height, bytesperpixel, fname);
371 
372  if (ret) {
373  const uint pnX = context_.pixel_number[_X];
374  const uint pnY = context_.pixel_number[_Y];
375  uchar *imgIFTA = nullptr;
376 
377  if (imgIFTA != nullptr) {
378  delete[] imgIFTA;
379  imgIFTA = nullptr;
380  }
381  uchar *imgTmp = loadAsImg(fname);
382  imgIFTA = new uchar[pnX * pnY * bytesperpixel];
383  memset(imgIFTA, 0, pnX * pnY * bytesperpixel);
384  imgScaleBilinear(imgTmp, imgIFTA, width, height, pnX, pnY, bytesperpixel);
385 
386  delete[] imgTmp;
387 
388  if (bRGB) imgRGB = imgIFTA;
389  else imgDepth = imgIFTA;
390  }
391  return ret;
392 }
393 
394 uchar ophIFTA::getMax(uchar *src, int width, int height)
395 {
396  const long long int size = width * height;
397  uchar max = 0;
398 
399  for (int i = 0; i < size; i++) {
400  if (*(src + i) > max) max = *(src + i);
401  }
402  return max;
403 }
ophIFTA()
Definition: ophIFTA.cpp:8
#define M_PI
Definition: define.h:52
#define OPH_FORWARD
Definition: define.h:66
virtual ~ophIFTA()
Definition: ophIFTA.cpp:23
T angle(void)
Definition: complex.h:387
bytesperpixel
Definition: Openholo.cpp:431
unsigned char uchar
Definition: typedef.h:64
uchar getMax(uchar *src, int width, int height)
Definition: ophIFTA.cpp:394
float Real
Definition: typedef.h:55
structure for 2-dimensional integer vector and its arithmetic.
Definition: ivec.h:66
#define _IM
Definition: complex.h:58
def k(wvl)
Definition: Depthmap.py:16
bool readConfig(const char *fname)
Definition: ophIFTA.cpp:326
const XMLNode * FirstChild() const
Get the first child node, or null if none exists.
Definition: tinyxml2.h:761
#define _X
Definition: define.h:92
#define _RE
Definition: complex.h:55
bool readImage(const char *fname, bool bRGB)
Definition: ophIFTA.cpp:368
bool normalize()
Definition: ophIFTA.cpp:272
Complex< T > & exp()
Definition: complex.h:395
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
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 angle(const std::vector< Complex< T >> &src, std::vector< T > &dst)
Definition: function.h:159
#define OPH_ESTIMATE
Definition: define.h:76
void encoding()
Definition: ophGen.cpp:985
#define ELAPSED_TIME(x, y)
Definition: function.h:59
#define OPH_BACKWARD
Definition: define.h:67
Definition: Bitmap.h:49
Real generateHologram()
Definition: ophIFTA.cpp:31
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
unsigned int uint
Definition: typedef.h:62
#define CUR_TIME
Definition: function.h:58