Openholo  v4.2
Open Source Digital Holographic Library
ophSigCH.cpp
Go to the documentation of this file.
1 #include "ophSigCH.h"
2 
4  : Nz(0)
5  , MaxIter(0)
6  , Tau(0)
7  , TolA(0)
8  , TvIter(0)
9 {
10 }
11 
12 bool ophSigCH::saveNumRec(const char *fname)
13 {
14  double gamma = 0.5;
15 
16  oph::uchar* intensityData;
17  intensityData = (oph::uchar*)malloc(sizeof(oph::uchar) * _cfgSig.rows * _cfgSig.cols);
18  double maxIntensity = 0.0;
19  double realVal = 0.0;
20  double imagVal = 0.0;
21  double intensityVal = 0.0;
22  for (int k = 0; k < Nz; k++)
23  {
24  for (int i = 0; i < _cfgSig.rows; i++)
25  {
26  for (int j = 0; j < _cfgSig.cols; j++)
27  {
28  realVal = NumRecRealImag(i, j + k*_cfgSig.cols);
29  imagVal = NumRecRealImag(i + _cfgSig.rows, j + k*_cfgSig.cols);
30  intensityVal = realVal*realVal + imagVal*imagVal;
31  if (intensityVal > maxIntensity)
32  {
33  maxIntensity = intensityVal;
34  }
35  }
36  }
37  }
38 
39  string fnamestr = fname;
40  int checktype = static_cast<int>(fnamestr.rfind("."));
41  char str[10];
42 
43  for (int k = 0; k < Nz; k++)
44  {
45  fnamestr = fname;
46  for (int i = _cfgSig.rows - 1; i >= 0; i--)
47  {
48  for (int j = 0; j < _cfgSig.cols; j++)
49  {
50  realVal = NumRecRealImag(i, j + k*_cfgSig.cols);
51  imagVal = NumRecRealImag(i + _cfgSig.rows, j + k*_cfgSig.cols);
52  intensityVal = realVal*realVal + imagVal*imagVal;
53  intensityData[i*_cfgSig.cols + j] = (uchar)(pow(intensityVal / maxIntensity,gamma)*255.0);
54  }
55  }
56  sprintf(str, "_%.2d", k);
57  fnamestr.insert(checktype, str);
58  saveAsImg(fnamestr.c_str(), 8, intensityData, _cfgSig.cols, _cfgSig.rows);
59  }
60 
61  return true;
62 }
63 
64 bool ophSigCH::setCHparam(vector<Real> &z, int maxIter, double tau, double tolA, int tvIter) {
65  MaxIter = maxIter;
66  Tau = tau;
67  TolA = tolA;
68  TvIter = tvIter;
69  Z.resize(z.size());
70  Z = z;
71  Nz = Z.size();
72  return true;
73 }
74 
75 bool ophSigCH::runCH(int complexHidx)
76 {
77  auto start_time = CUR_TIME;
78 
79  matrix<Real> realimagplaneinput(_cfgSig.rows*2, _cfgSig.cols);
80  c2ri(ComplexH[complexHidx], realimagplaneinput);
81 
82  NumRecRealImag.resize(_cfgSig.rows * 2, _cfgSig.cols*Nz);
83  twist(realimagplaneinput, NumRecRealImag);
84 
85  auto end_time = CUR_TIME;
86 
87  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
88 
89  LOG("Implement time : %.5lf sec\n", during_time);
90 
91  return true;
92 }
93 
94 void ophSigCH::tvdenoise(matrix<Real>& input, double lam, int iters, matrix<Real>& output)
95 {
96  double dt = 0.25;
97  int nr = input.size(_X);
98  int nc = input.size(_Y);
99  matrix<Real> divp(nr, nc);
100  divp.zeros();
101 
102  matrix<Real> z(nr, nc);
103  matrix<Real> z1, z2;
104  z1.resize(nr, nc);
105  z2.resize(nr, nc);
106 
107  matrix<Real> p1, p2;
108  p1.resize(nr, nc);
109  p2.resize(nr, nc);
110 
111  matrix<Real> denom(nr, nc);
112 
113  for (int iter = 0; iter < iters; iter++)
114  {
115  for (int i = 0; i < nr; i++)
116  {
117  for (int j = 0; j < nc; j++)
118  {
119  z(i,j) = divp(i, j) - input(i, j)*lam;
120  }
121  }
122 
124  for (int i = 0; i < nr - 1; i++)
125  {
126  for (int j = 0; j < nc - 1; j++)
127  {
128  z1(i, j) = z(i, j + 1) - z(i, j);
129  z2(i, j) = z(i + 1, j) - z(i, j);
130  denom(i, j) = 1 + dt*sqrt(z1(i, j)*z1(i, j) + z2(i, j)*z2(i, j));
131  }
132  }
133  for (int i = 0; i < nr-1; i++)
134  {
135  z1(i, nc-1) = 0.0;
136  z2(i, nc-1) = z(i + 1, nc-1) - z(i, nc-1);
137  denom(i, nc-1) = 1 + dt*sqrt(z1(i, nc-1)*z1(i, nc-1) + z2(i, nc-1)*z2(i, nc-1));
138  }
139  for (int j = 0; j < nc - 1; j++)
140  {
141  z1(nr-1, j) = z(nr-1, j+1) - z(nr-1, j);
142  z2(nr-1, j) = 0.0;
143  denom(nr-1, j) = 1 + dt*sqrt(z1(nr-1, j)*z1(nr-1, j) + z2(nr-1, j)*z2(nr-1, j));
144  }
145  denom(nr-1, nc-1) = 1.0;
146  z1(nr - 1, nc - 1) = 0.0;
147  z2(nr - 1, nc - 1) = 0.0;
149 
150 
151  for (int i = 0; i < nr ; i++)
152  {
153  for (int j = 0; j < nc; j++)
154  {
155  p1(i, j) = (p1(i, j) + dt*z1(i, j)) / denom(i, j);
156  p2(i, j) = (p2(i, j) + dt*z2(i, j)) / denom(i, j);
157  }
158  }
159 
161  for (int i = 1; i < nr; i++)
162  {
163  for (int j = 1; j < nc; j++)
164  {
165  divp(i, j) = p1(i, j) - p1(i, j - 1) + p2(i, j) - p2(i - 1, j);
166  }
167  }
168  for (int i = 1; i < nr; i++)
169  {
170  divp(i, 0) = p2(i, 0) - p2(i-1, 0);
171  }
172  for (int j = 1; j < nc; j++)
173  {
174  divp(0, j) = p1(0, j) - p1(0, j - 1);
175  }
176  divp(0, 0) = 0.0;
178 
179  }
180 
181  for (int i = 0; i < input.size[_X]; i++)
182  {
183  for (int j = 0; j < input.size[_Y]; j++)
184  {
185  output(i, j) = input(i,j) - divp(i, j)/lam;
186  }
187  }
188 
189 }
190 
191 double ophSigCH::tvnorm(matrix<Real>& input)
192 {
193  double sqrtsum = 0.0;
194  int nr = input.size[_X];
195  int nc = input.size[_Y];
196 
197  for (int i = 0; i < nr-1; i++)
198  {
199  for (int j = 0; j < nc-1; j++)
200  {
201  sqrtsum += sqrt((input(i,j)-input(i+1,j))*(input(i,j)-input(i+1,j)) + (input(i,j)-input(i,j+1))*(input(i,j)-input(i,j+1)));
202  }
203  }
204  for (int i = 0; i < nr-1; i++)
205  {
206  sqrtsum += sqrt(pow(input(i, nc-1) - input(i + 1, nc-1), 2) + pow(input(i, nc-1) - input(i, 0), 2));
207  }
208  for (int j = 0; j < nc - 1; j++)
209  {
210  sqrtsum += sqrt(pow(input(nr-1, j) - input(0, j), 2) + pow(input(nr-1, j) - input(nr-1, j+1), 2));
211  }
212  sqrtsum += sqrt(pow(input(nr-1, nc-1) - input(0, nc-1), 2) + pow(input(nr-1, nc-1) - input(nr-1, 0), 2));
213 
214  return sqrtsum;
215 }
216 
217 void ophSigCH::c2ri(matrix<Complex<Real>>& complexinput, matrix<Real>& realimagoutput)
218 {
219  for (int i = 0; i < complexinput.size[_X]; i++)
220  {
221  for (int j = 0; j < complexinput.size[_Y]; j++)
222  {
223  realimagoutput(i, j) = complexinput(i, j)[_RE];
224  realimagoutput(i + complexinput.size[_X], j) = complexinput(i, j)[_IM];
225  }
226  }
227 
228 }
229 
230 void ophSigCH::ri2c(matrix<Real>& realimaginput, matrix<Complex<Real>>& complexoutput)
231 {
232  for (int i = 0; i < complexoutput.size[_X]; i++)
233  {
234  for (int j = 0; j < complexoutput.size[_Y]; j++)
235  {
236  complexoutput(i, j)[_RE] = realimaginput(i, j);
237  complexoutput(i, j)[_IM] = realimaginput(i + complexoutput.size[_X], j);
238  }
239  }
240 }
241 
242 void ophSigCH::volume2plane(matrix<Real>& realimagvolumeinput, vector<Real> z, matrix<Real>& realimagplaneoutput)
243 {
244  int nz = z.size();
245  int nr = realimagvolumeinput.size(_X) >> 1; // real imag
246  int nc = realimagvolumeinput.size(_Y) / nz;
247 
248  matrix<Complex<Real>> complexTemp(nr,nc);
249  matrix<Complex<Real>> complexAccum(nr, nc);
250  complexAccum.zeros();
251 
252  for (int k = 0; k < nz; k++)
253  {
254  int temp = k * nc;
255  for (int i = 0; i < nr; i++)
256  {
257  for (int j = 0; j < nc; j++)
258  {
259  complexTemp(i, j)[_RE] = realimagvolumeinput(i, j + temp);
260  complexTemp(i, j)[_IM] = realimagvolumeinput(i + nr, j + temp);
261  }
262  }
263  OphComplexField result = propagationHoloAS(complexTemp, static_cast<float>(z.at(k)));
264  complexAccum = complexAccum + result;
265  }
266  c2ri(complexAccum, realimagplaneoutput);
267 }
268 
269 
270 void ophSigCH::plane2volume(matrix<Real>& realimagplaneinput, vector<Real> z, matrix<Real>& realimagplaneoutput)
271 {
272  int nr = realimagplaneinput.size(_X)/2;
273  int nc = realimagplaneinput.size(_Y);
274  int nz = z.size();
275 
276  matrix<Complex<Real>> complexplaneinput(nr, nc);
277  for (int i = 0; i < nr; i++)
278  {
279  for (int j = 0; j < nc; j++)
280  {
281  complexplaneinput(i, j)[_RE] = realimagplaneinput(i, j);
282  complexplaneinput(i, j)[_IM] = realimagplaneinput(i + nr, j);
283  }
284  }
285 
286 
287  matrix<Complex<Real>> temp(nr, nc);
288  for (int k = 0; k < nz; k++)
289  {
290  temp = propagationHoloAS(complexplaneinput, static_cast<float>(-z.at(k)));
291  for (int i = 0; i < nr; i++)
292  {
293  for (int j = 0; j < nc; j++)
294  {
295  realimagplaneoutput(i, j + k*nc) = temp(i, j)[_RE];
296  realimagplaneoutput(i + nr, j + k*nc) = temp(i, j)[_IM];
297  }
298  }
299  }
300 
301 }
302 
303 void ophSigCH::convert3Dto2D(matrix<Complex<Real>>* complex3Dinput, int nz, matrix<Complex<Real>>& complex2Doutput)
304 {
305  int nr = complex3Dinput[0].size(_X);
306  int nc = complex3Dinput[0].size(_Y);
307 
308  for (int i = 0; i < nr; i++)
309  {
310  for (int j = 0; j < nc; j++)
311  {
312  for (int k = 0; k < nz; k++)
313  {
314  complex2Doutput(i, j + k*nc) = complex3Dinput[k](i, j);
315  }
316  }
317  }
318 }
319 
320 void ophSigCH::convert2Dto3D(matrix<Complex<Real>>& complex2Dinput, int nz, matrix<Complex<Real>>* complex3Doutput)
321 {
322  int nr = complex2Dinput.size(_X);
323  int nc = complex2Dinput.size(_Y)/nz;
324 
325  for (int k = 0; k < nz; k++)
326  {
327  for (int i = 0; i < nr; i++)
328  {
329  for (int j = 0; j < nc; j++)
330  {
331  complex3Doutput[k](i, j) = complex2Dinput(i, j + k*nc);
332  }
333  }
334  }
335 }
336 
337 void ophSigCH::twist(matrix<Real>& realimagplaneinput, matrix<Real>& realimagvolumeoutput)
338 {
339  //
340  int nrv = realimagvolumeoutput.size(_X);
341  int ncv = realimagvolumeoutput.size(_Y);
342  int nrp = realimagplaneinput.size(_X);
343  int ncp = realimagplaneinput.size(_Y);
344 
345 
346  // TWiST
347  double stopCriterion = 1;
348  double tolA = TolA;
349  int maxiter = MaxIter;
350  int miniter = MaxIter;
351  int init = 2;
352  bool enforceMonotone = 1;
353  bool compute_mse = 0;
354  bool plot_ISNR = 0;
355  bool verbose = 1;
356  double alpha = 0;
357  double beta = 0;
358  bool sparse = 1;
359  double tolD = 0.001;
360  double phi_l1 = 0;
361  double psi_ok = 0;
362  double lam1 = 1e-4;
363  double lamN = 1;
364  double tau = Tau;
365  int tv_iters = TvIter;
366 
367  bool for_ever = 1;
368  double max_svd = 1;
369 
370  // twist parameters
371  double rho0 = (1. - lam1 / lamN) / (1. + lam1 / lamN);
372  if (alpha == 0) {
373  alpha = 2. / (1 + sqrt(1 - pow(rho0, 2)));
374  }
375  if (beta == 0) {
376  beta = alpha*2. / (lam1 + lamN);
377  }
378 
379  double prev_f = 0.0;
380  double f = 0.0;
381 
382  double criterion = 0.0;
383 
384  // initialization
385  plane2volume(realimagplaneinput, Z, realimagvolumeoutput);
386 
387  // compute and sotre initial value of the objective function
388  matrix<Real> resid(nrp, ncp);
389  volume2plane(realimagvolumeoutput, Z, resid);
390  for (int i = 0; i < nrp; i++)
391  {
392  for (int j = 0; j < ncp; j++)
393  {
394  resid(i, j) = realimagplaneinput(i, j) - resid(i, j);
395  }
396  }
397  prev_f = 0.5*matrixEleSquareSum(resid) + tau*tvnorm(realimagvolumeoutput);
398 
399  //
400  int iter = 0;
401  bool cont_outer = 1;
402 
403  if (verbose)
404  {
405  LOG("initial objective = %10.6e\n", prev_f);
406  }
407 
408  int IST_iters = 0;
409  int TwIST_iters = 0;
410 
411  // initialize
412  matrix<Real> xm1, xm2;
413  xm1.resize(nrv, ncv);
414  xm2.resize(nrv, ncv);
415  xm1 = realimagvolumeoutput;
416  xm2 = realimagvolumeoutput;
417 
418  // TwIST iterations
419  matrix<Real> grad, temp_volume;
420  grad.resize(nrv, ncv);
421  temp_volume.resize(nrv, ncv);
422  while (cont_outer)
423  {
424  plane2volume(resid, Z, grad);
425  while (for_ever)
426  {
427  for (int i = 0; i < nrv; i++)
428  {
429  for (int j = 0; j < ncv; j++)
430  {
431  temp_volume(i, j) = xm1(i, j) + grad(i, j) / max_svd;
432  }
433  }
434  tvdenoise(temp_volume, 2.0 / (tau / max_svd), tv_iters, realimagvolumeoutput);
435 
436  if ((IST_iters >= 2) | (TwIST_iters != 0))
437  {
438  if (sparse)
439  {
440  for (int i = 0; i < nrv; i++)
441  {
442  for (int j = 0; j < ncv; j++)
443  {
444  if (realimagvolumeoutput(i, j) == 0)
445  {
446  xm1(i, j) = 0.0;
447  xm2(i, j) = 0.0;
448  }
449  }
450  }
451 
452  }
453  // two step iteration
454  for (int i = 0; i < nrv; i++)
455  {
456  for (int j = 0; j < ncv; j++)
457  {
458  xm2(i, j) = (alpha - beta)*xm1(i, j) + (1.0 - alpha)*xm2(i, j) + beta*realimagvolumeoutput(i, j);
459  }
460  }
461  // compute residual
462  volume2plane(xm2, Z, resid);
463 
464  for (int i = 0; i < nrp; i++)
465  {
466  for (int j = 0; j < ncp; j++)
467  {
468  resid(i, j) = realimagplaneinput(i, j) - resid(i, j);
469  }
470  }
471  f = 0.5*matrixEleSquareSum(resid) + tau*tvnorm(xm2);
472  if ((f > prev_f) & (enforceMonotone))
473  {
474  TwIST_iters = 0;
475  }
476  else
477  {
478  TwIST_iters += 1;
479  IST_iters = 0;
480  realimagvolumeoutput = xm2;
481  if (TwIST_iters % 10000 == 0)
482  {
483  max_svd = 0.9*max_svd;
484  }
485  break;
486  }
487  }
488  else
489  {
490  volume2plane(realimagvolumeoutput, Z, resid);
491  for (int i = 0; i < nrp; i++)
492  {
493  for (int j = 0; j < ncp; j++)
494  {
495  resid(i, j) = realimagplaneinput(i, j) - resid(i, j);
496  }
497  }
498  f = 0.5*matrixEleSquareSum(resid) + tau*tvnorm(realimagvolumeoutput);
499  if (f > prev_f)
500  {
501  max_svd = 2 * max_svd;
502  if (verbose)
503  {
504  LOG("Incrementing S = %2.2e\n", max_svd);
505  }
506  IST_iters = 0;
507  TwIST_iters = 0;
508  }
509  else
510  {
511  TwIST_iters += 1;
512  break;
513  }
514  }
515 
516  }
517  xm2 = xm1;
518  xm1 = realimagvolumeoutput;
519 
520  criterion = (abs(f - prev_f)) / prev_f;
521 
522  cont_outer = ((iter <= maxiter) & (criterion > tolA));
523  if (iter <= miniter)
524  {
525  cont_outer = 1;
526  }
527 
528  iter += 1;
529  prev_f = f;
530 
531  if (verbose)
532  {
533  LOG("Iteration=%4d, objective=%9.5e, criterion=%7.3e\n", iter, f, criterion / tolA);
534  }
535  }
536 
537  if (verbose)
538  {
539  double sum_abs_x=0.0;
540  for (int i = 0; i < nrv; i++)
541  {
542  for (int j = 0; j < ncv; j++)
543  {
544  sum_abs_x += abs(realimagvolumeoutput(i, j));
545  }
546  }
547  LOG("\n Finishied the main algorithm!\n");
548  LOG("||Ax-y||_2 = %10.3e\n", matrixEleSquareSum(resid));
549  LOG("||x||_1 = %10.3e\n", sum_abs_x);
550  LOG("Objective function = %10.3e\n", f);
551  }
552 }
553 
554 double ophSigCH::matrixEleSquareSum(matrix<Real>& input)
555 {
556  double output = 0.0;
557  for (int i = 0; i < input.size(_X); i++)
558  {
559  for (int j = 0; j < input.size(_Y); j++)
560  {
561  output += input(i, j)*input(i, j);
562  }
563  }
564  return output;
565 }
566 
567 bool ophSigCH::readConfig(const char* fname)
568 {
569  LOG("CH Reading....%s...\n", fname);
570 
571  /*XML parsing*/
572  tinyxml2::XMLDocument xml_doc;
573  tinyxml2::XMLNode *xml_node;
574 
575  if (!checkExtension(fname, ".xml"))
576  {
577  LOG("file's extension is not 'xml'\n");
578  return false;
579  }
580  auto ret = xml_doc.LoadFile(fname);
581  if (ret != tinyxml2::XML_SUCCESS)
582  {
583  LOG("Failed to load file \"%s\"\n", fname);
584  return false;
585  }
586 
587  xml_node = xml_doc.FirstChild();
588 
589  (xml_node->FirstChildElement("rows"))->QueryIntText(&_cfgSig.rows);
590  (xml_node->FirstChildElement("cols"))->QueryIntText(&_cfgSig.cols);
591  (xml_node->FirstChildElement("width"))->QueryFloatText(&_cfgSig.width);
592  (xml_node->FirstChildElement("height"))->QueryFloatText(&_cfgSig.height);
593  (xml_node->FirstChildElement("wavelength"))->QueryDoubleText(&_cfgSig.wavelength[0]);
594  (xml_node->FirstChildElement("nz"))->QueryIntText(&Nz);
595  (xml_node->FirstChildElement("maxiter"))->QueryIntText(&MaxIter);
596  (xml_node->FirstChildElement("tau"))->QueryDoubleText(&Tau);
597  (xml_node->FirstChildElement("tolA"))->QueryDoubleText(&TolA);
598  (xml_node->FirstChildElement("tv_iters"))->QueryIntText(&TvIter);
599 
600  double zmin, dz;
601  (xml_node->FirstChildElement("zmin"))->QueryDoubleText(&zmin);
602  (xml_node->FirstChildElement("dz"))->QueryDoubleText(&dz);
603 
604  Z.resize(Nz);
605  for (int i = 0; i < Nz; i++)
606  {
607  Z.at(i) = zmin + i*dz;
608  }
609 
610  return true;
611 }
612 
613 bool ophSigCH::loadCHtemp(const char * real, const char * imag, uint8_t bitpixel)
614 {
615  string realname = real;
616  string imagname = imag;
617  int checktype = static_cast<int>(realname.rfind("."));
618  matrix<Real> realMat[3], imagMat[3];
619 
620  std::string realtype = realname.substr(checktype + 1, realname.size());
621  std::string imgtype = imagname.substr(checktype + 1, realname.size());
622  size_t nRead;
623 
624  if (realtype != imgtype) {
625  LOG("failed : The data type between real and imaginary is different!\n");
626  return false;
627  }
628  if (realtype == "bmp")
629  {
630  FILE* freal = nullptr;
631  FILE* fimag = nullptr;
632  fileheader hf;
634  freal = fopen(realname.c_str(), "rb");
635  fimag = fopen(imagname.c_str(), "rb");
636  if (freal == nullptr)
637  {
638  LOG("real bmp file open fail!\n");
639  return false;
640  }
641  if (fimag == nullptr)
642  {
643  LOG("imaginary bmp file open fail!\n");
644  return false;
645  }
646  nRead = fread(&hf, sizeof(fileheader), 1, freal);
647  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, freal);
648  nRead = fread(&hf, sizeof(fileheader), 1, fimag);
649  nRead = fread(&hInfo, sizeof(bitmapinfoheader), 1, fimag);
650 
651  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
652  if ((hInfo.height == 0) || (hInfo.width == 0))
653  {
654  LOG("bmp header is empty!\n");
657  if (_cfgSig.rows == 0 || _cfgSig.cols == 0)
658  {
659  LOG("check your parameter file!\n");
660  return false;
661  }
662  }
663  if ((_cfgSig.rows != (int)hInfo.height) || (_cfgSig.cols != (int)hInfo.width)) {
664  LOG("image size is different!\n");
667  LOG("changed parameter of size %d x %d\n", _cfgSig.cols, _cfgSig.rows);
668  }
669  hInfo.bitsperpixel = bitpixel;
670  if (bitpixel == 8)
671  {
672  rgbquad palette[256];
674  nRead = fread(palette, sizeof(rgbquad), 256, freal);
675  nRead = fread(palette, sizeof(rgbquad), 256, fimag);
676 
677  realMat[0].resize(hInfo.height, hInfo.width);
678  imagMat[0].resize(hInfo.height, hInfo.width);
679  ComplexH[0].resize(hInfo.height, hInfo.width);
680  }
681  else
682  {
683  ComplexH = new OphComplexField[3];
684  realMat[0].resize(hInfo.height, hInfo.width);
685  imagMat[0].resize(hInfo.height, hInfo.width);
686  ComplexH[0].resize(hInfo.height, hInfo.width);
687 
688  realMat[1].resize(hInfo.height, hInfo.width);
689  imagMat[1].resize(hInfo.height, hInfo.width);
690  ComplexH[1].resize(hInfo.height, hInfo.width);
691 
692  realMat[2].resize(hInfo.height, hInfo.width);
693  imagMat[2].resize(hInfo.height, hInfo.width);
694  ComplexH[2].resize(hInfo.height, hInfo.width);
695  }
696  size_t size = hInfo.width * hInfo.height * (hInfo.bitsperpixel >> 3);
697  uchar* realdata = (uchar*)malloc(sizeof(uchar) * size);
698  uchar* imagdata = (uchar*)malloc(sizeof(uchar) * size);
699  nRead = fread(realdata, sizeof(uchar), size, freal);
700  nRead = fread(imagdata, sizeof(uchar), size, fimag);
701 
702  fclose(freal);
703  fclose(fimag);
704 
705  for (int i = hInfo.height - 1; i >= 0; i--)
706  {
707  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
708  {
709  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
710  {
711  realMat[z](hInfo.height - i - 1, j) = (double)realdata[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
712  imagMat[z](hInfo.height - i - 1, j) = (double)imagdata[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
713  }
714  }
715  }
716  LOG("file load complete!\n");
717 
718  free(realdata);
719  free(imagdata);
720  }
721  else if (realtype == "bin")
722  {
723  if (bitpixel == 8)
724  {
725 
726  ifstream freal(realname, ifstream::binary);
727  ifstream fimag(imagname, ifstream::binary);
728  realMat[0].resize(_cfgSig.rows, _cfgSig.cols); imagMat[0].resize(_cfgSig.rows, _cfgSig.cols);
729  ComplexH[0].resize(_cfgSig.rows, _cfgSig.cols);
730  int total = _cfgSig.rows*_cfgSig.cols;
731  double *realdata = new double[total];
732  double *imagdata = new double[total];
733  int i = 0;
734  freal.read(reinterpret_cast<char*>(realdata), sizeof(double) * total);
735  fimag.read(reinterpret_cast<char*>(imagdata), sizeof(double) * total);
736 
737  for (int col = 0; col < _cfgSig.cols; col++)
738  {
739  for (int row = 0; row < _cfgSig.rows; row++)
740  {
741  realMat[0](row, col) = realdata[_cfgSig.rows*col + row];
742  imagMat[0](row, col) = imagdata[_cfgSig.rows*col + row];
743  }
744  }
745 
746  freal.close();
747  fimag.close();
748  delete[]realdata;
749  delete[]imagdata;
750  }
751  else if (bitpixel == 24)
752  {
753  realMat[0].resize(_cfgSig.rows, _cfgSig.cols);
754  imagMat[0].resize(_cfgSig.rows, _cfgSig.cols);
755  ComplexH[0].resize(_cfgSig.rows, _cfgSig.cols);
756 
757  realMat[1].resize(_cfgSig.rows, _cfgSig.cols);
758  imagMat[1].resize(_cfgSig.rows, _cfgSig.cols);
759  ComplexH[1].resize(_cfgSig.rows, _cfgSig.cols);
760 
761  realMat[2].resize(_cfgSig.rows, _cfgSig.cols);
762  imagMat[2].resize(_cfgSig.rows, _cfgSig.cols);
763  ComplexH[2].resize(_cfgSig.rows, _cfgSig.cols);
764 
765  int total = _cfgSig.rows*_cfgSig.cols;
766 
767 
768  string RGB_name[] = { "_B","_G","_R" };
769  double *realdata = new double[total];
770  double *imagdata = new double[total];
771  char* context = nullptr;
772 
773  for (int z = 0; z < (bitpixel / 8); z++)
774  {
775  ifstream freal(strtok((char*)realname.c_str(), ".") + RGB_name[z] + "bin", ifstream::binary);
776  ifstream fimag(strtok((char*)imagname.c_str(), ".") + RGB_name[z] + "bin", ifstream::binary);
777 
778  freal.read(reinterpret_cast<char*>(realdata), sizeof(double) * total);
779  fimag.read(reinterpret_cast<char*>(imagdata), sizeof(double) * total);
780 
781  for (int col = 0; col < _cfgSig.cols; col++)
782  {
783  for (int row = 0; row < _cfgSig.rows; row++)
784  {
785  realMat[z](row, col) = realdata[_cfgSig.rows*col + row];
786  imagMat[z](row, col) = imagdata[_cfgSig.rows*col + row];
787  }
788  }
789  freal.close();
790  fimag.close();
791  }
792  delete[] realdata;
793  delete[] imagdata;
794  }
795  }
796  else
797  {
798  LOG("wrong type\n");
799  }
800 
803  //nomalization
804  //double realout, imagout;
805  for (int z = 0; z < (bitpixel) / 8; z++)
806  {
807  for (int i = 0; i < _cfgSig.rows; i++)
808  {
809  for (int j = 0; j < _cfgSig.cols; j++)
810  {
811  ComplexH[z](i, j)[_RE] = realMat[z](i, j)/255.0*2.0-1.0;
812  ComplexH[z](i, j)[_IM] = imagMat[z](i, j)/255.0*2.0-1.0;
813 
814  }
815  }
816  }
817  LOG("data nomalization complete\n");
818 
819  return true;
820 }
821 
822 matrix<Complex<Real>> ophSigCH::propagationHoloAS(matrix<Complex<Real>> complexH, float depth) {
823  int nr = _cfgSig.rows;
824  int nc = _cfgSig.cols;
825 
826  double dr = _cfgSig.height / _cfgSig.rows;
827  double dc = _cfgSig.width / _cfgSig.cols;
828 
829  int nr2 = 2 * nr;
830  int nc2 = 2 * nc;
831 
832  oph::matrix<oph::Complex<Real>> src2(nr2,nc2); // prepare complex matrix with 2x size (to prevent artifacts caused by circular convolution)
833  src2 * 0; // initialize to 0
834 
835  int iStart = nr / 2 - 1;
836  int jStart = nc / 2 - 1;
837  for (int i = 0; i < nr; i++)
838  {
839  for (int j = 0; j < nc; j++)
840  {
841  src2(i+iStart, j+jStart)[_RE] = complexH(i, j)[_RE];
842  src2(i+iStart, j+jStart)[_IM] = complexH(i, j)[_IM];
843 
844  }
845  }
846 
847  double dfr = 1.0 / (((double)nr2)*dr); // spatial frequency step in src2
848  double dfc = 1.0 / (((double)nc2)*dc);
849 
850  matrix<Complex<Real>> propKernel(nr2, nc2);
851  double fz = 0;
852  for (int i = 0; i < nr2; i++)
853  {
854  for (int j = 0; j < nc2; j++)
855  {
856  fz = sqrt(pow(1.0 / _cfgSig.wavelength[0], 2) - pow((i - nr2 / 2.0 + 1.0)*dfr, 2) - pow((j - nc2 / 2.0 + 1.0)*dfc, 2));
857  propKernel(i, j)[_RE] = cos(2 * M_PI*depth*fz);
858  propKernel(i, j)[_IM] = sin(2 * M_PI*depth*fz);
859  }
860  }
861 
862  matrix<Complex<Real>> src2temp(nr2, nc2);
863  matrix<Complex<Real>> FFZP(nr2, nc2);
864  matrix<Complex<Real>> FFZP2(nr2, nc2);
865  fftShift(src2, src2temp);
866  fft2(src2temp, FFZP, OPH_FORWARD);
867  fftShift(FFZP, FFZP2);
868  FFZP2.mulElem(propKernel);
869 
870  fftShift(FFZP2, FFZP);
871  fft2(FFZP, src2temp, OPH_BACKWARD);
872  fftShift(src2temp, src2);
873 
874  matrix<Complex<Real>> dst(nr, nc);
875  for (int i = 0; i < nr; i++)
876  {
877  for (int j = 0; j < nc; j++)
878  {
879  dst(i,j)[_RE] = src2(i + iStart , j + jStart)[_RE];
880  dst(i,j)[_IM] = src2(i + iStart , j + jStart)[_IM];
881  }
882  }
883  return dst;
884 }
oph::matrix< Complex< Real > > OphComplexField
Definition: mat.h:420
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:113
#define M_PI
Definition: define.h:52
#define OPH_FORWARD
Definition: define.h:66
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:135
size_t nRead
Definition: Openholo.cpp:424
vector< Real > Z
Definition: ophSigCH.h:82
int MaxIter
Definition: ophSigCH.h:77
unsigned char uchar
Definition: typedef.h:64
void plane2volume(matrix< Real > &realimagplaneinput, vector< Real > z, matrix< Real > &realimagplaneoutput)
Definition: ophSigCH.cpp:270
void ri2c(matrix< Real > &realimaginput, matrix< Complex< Real >> &complexoutput)
Definition: ophSigCH.cpp:230
void c2ri(matrix< Complex< Real >> &complexinput, matrix< Real > &realimagoutput)
Definition: ophSigCH.cpp:217
Real wavelength[3]
Definition: ophSig.h:76
double TolA
Definition: ophSigCH.h:79
bool runCH(int complexHidx)
Definition: ophSigCH.cpp:75
bool checkExtension(const char *fname, const char *ext)
Functions for extension checking.
Definition: Openholo.cpp:86
#define _IM
Definition: complex.h:58
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
ophSigCH(void)
Constructor.
Definition: ophSigCH.cpp:3
void fft2(matrix< Complex< T >> &src, matrix< Complex< T >> &dst, int sign=OPH_FORWARD, uint flag=OPH_ESTIMATE)
Function for Fast Fourier transform 2D.
Definition: ophSig.cpp:149
fclose(infile)
bool saveNumRec(const char *fname)
Definition: ophSigCH.cpp:12
ophSigConfig _cfgSig
Definition: ophSig.h:486
#define _X
Definition: define.h:92
matrix< Complex< Real > > propagationHoloAS(matrix< Complex< Real >> complexH, float depth)
Definition: ophSigCH.cpp:822
uint8_t signature[2]
Definition: struct.h:51
Definition: struct.h:69
int cols
Definition: ophSig.h:69
void convert2Dto3D(matrix< Complex< Real >> &complex2Dinput, int nz, matrix< Complex< Real >> *complex3Doutput)
Definition: ophSigCH.cpp:320
int rows
Definition: ophSig.h:70
double matrixEleSquareSum(matrix< Real > &input)
Definition: ophSigCH.cpp:554
uint16_t bitsperpixel
Definition: struct.h:61
OphComplexField * ComplexH
Definition: ophSig.h:487
void twist(matrix< Real > &realimagplaneinput, matrix< Real > &realimagvolumeoutput)
Definition: ophSigCH.cpp:337
bool readConfig(const char *fname)
Definition: ophSigCH.cpp:567
#define _RE
Definition: complex.h:55
int TvIter
Definition: ophSigCH.h:80
Real_t height
Definition: ophSig.h:72
void fftShift(matrix< Complex< Real >> &src, matrix< Complex< Real >> &dst)
Shift zero-frequency component to center of spectrum.
Definition: ophSig.h:849
void volume2plane(matrix< Real > &realimagvolumeinput, vector< Real > z, matrix< Real > &realimagplaneoutput)
Definition: ophSigCH.cpp:242
#define _Y
Definition: define.h:96
XMLError LoadFile(const char *filename)
Definition: tinyxml2.cpp:2150
bool setCHparam(vector< Real > &z, int maxIter, double tau, double tolA, int tvIter)
Definition: ophSigCH.cpp:64
uint32_t width
Definition: struct.h:58
bitmapinfoheader hInfo
Definition: Openholo.cpp:423
uint32_t height
Definition: struct.h:59
int Nz
Definition: ophSigCH.h:76
void tvdenoise(matrix< Real > &input, double lam, int iters, matrix< Real > &output)
Definition: ophSigCH.cpp:94
bool loadCHtemp(const char *real, const char *imag, uint8_t bitpixel)
Definition: ophSigCH.cpp:613
fileheader hf
Definition: Openholo.cpp:422
matrix< Real > NumRecRealImag
Definition: ophSigCH.h:81
Real_t width
Definition: ophSig.h:71
#define OPH_BACKWARD
Definition: define.h:67
const XMLElement * FirstChildElement(const char *name=0) const
Definition: tinyxml2.cpp:940
void convert3Dto2D(matrix< Complex< Real >> *complex3Dinput, int nz, matrix< Complex< Real >> &complex2Doutput)
Definition: ophSigCH.cpp:303
double tvnorm(matrix< Real > &input)
Definition: ophSigCH.cpp:191
#define CUR_TIME
Definition: function.h:58
double Tau
Definition: ophSigCH.h:78