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