Openholo  v1.0
Open Source Digital Holographic Library
ophSigPU.cpp
Go to the documentation of this file.
1 #include "ophSigPU.h"
2 
4 {
5 }
6 
7 bool ophSigPU::setPUparam(int maxBoxRadius)
8 {
9  MaxBoxRadius = maxBoxRadius;
10  return false;
11 }
12 
13 bool ophSigPU::loadPhaseOriginal(const char * fname, int bitpixel)
14 {
15  string fnamestr = fname;
16  int checktype = static_cast<int>(fnamestr.rfind("."));
17  matrix<Real> phaseMat;
18 
19  std::string filetype = fnamestr.substr(checktype + 1, fnamestr.size());
20 
21  if (filetype == "bmp")
22  {
23  FILE *fphase;
24  fileheader hf;
25  bitmapinfoheader hInfo;
26  fopen_s(&fphase, fnamestr.c_str(), "rb");
27  if (!fphase)
28  {
29  LOG("real bmp file open fail!\n");
30  return false;
31  }
32  fread(&hf, sizeof(fileheader), 1, fphase);
33  fread(&hInfo, sizeof(bitmapinfoheader), 1, fphase);
34 
35  if (hf.signature[0] != 'B' || hf.signature[1] != 'M') { LOG("Not BMP File!\n"); }
36  if ((hInfo.height == 0) || (hInfo.width == 0))
37  {
38  LOG("bmp header is empty!\n");
39  hInfo.height = _cfgSig.rows;
40  hInfo.width = _cfgSig.cols;
41  if (_cfgSig.rows == 0 || _cfgSig.cols == 0)
42  {
43  LOG("check your parameter file!\n");
44  return false;
45  }
46  }
47  if ((_cfgSig.rows != hInfo.height) || (_cfgSig.cols != hInfo.width)) {
48  LOG("image size is different!\n");
49  _cfgSig.rows = hInfo.height;
50  _cfgSig.cols = hInfo.width;
51  LOG("changed parameter of size %d x %d\n", _cfgSig.cols, _cfgSig.rows);
52  }
53  hInfo.bitsperpixel = bitpixel;
54  if (bitpixel == 8)
55  {
56  rgbquad palette[256];
57  fread(palette, sizeof(rgbquad), 256, fphase);
58 
59  phaseMat.resize(hInfo.height, hInfo.width);
60  }
61  else
62  {
63  LOG("currently only 8 bitpixel is supported.");
64  /*
65  phaseMat[0].resize(hInfo.height, hInfo.width);
66  phaseMat[1].resize(hInfo.height, hInfo.width);
67  phaseMat[2].resize(hInfo.height, hInfo.width); */
68  }
69 
70  uchar* phasedata = (uchar*)malloc(sizeof(uchar)*hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8));
71 
72  fread(phasedata, sizeof(uchar), hInfo.width*hInfo.height*(hInfo.bitsperpixel / 8), fphase);
73 
74  fclose(fphase);
75 
76  for (int i = hInfo.height - 1; i >= 0; i--)
77  {
78  for (int j = 0; j < static_cast<int>(hInfo.width); j++)
79  {
80  for (int z = 0; z < (hInfo.bitsperpixel / 8); z++)
81  {
82  phaseMat(hInfo.height - i - 1, j) = (double)phasedata[i*hInfo.width*(hInfo.bitsperpixel / 8) + (hInfo.bitsperpixel / 8)*j + z];
83  }
84  }
85  }
86  LOG("file load complete!\n");
87 
88  free(phasedata);
89  }
90  else if (filetype == "bin")
91  {
92  if (bitpixel == 8)
93  {
94 
95  ifstream fphase(fnamestr, ifstream::binary);
96  phaseMat.resize(_cfgSig.rows, _cfgSig.cols);
97  int total = _cfgSig.rows*_cfgSig.cols;
98  double *phasedata = new double[total];
99  int i = 0;
100  fphase.read(reinterpret_cast<char*>(phasedata), sizeof(double) * total);
101 
102  for (int col = 0; col < _cfgSig.cols; col++)
103  {
104  for (int row = 0; row < _cfgSig.rows; row++)
105  {
106  phaseMat(row, col) = phasedata[_cfgSig.rows*col + row];
107  }
108  }
109 
110  fphase.close();
111  delete[]phasedata;
112  }
113  else if (bitpixel == 24)
114  {
115  LOG("currently only 8 bitpixel is supported.");
116  /*
117  phaseMat[0].resize(_cfgSig.rows, _cfgSig.cols);
118  phaseMat[1].resize(_cfgSig.rows, _cfgSig.cols);
119  phaseMat[2].resize(_cfgSig.rows, _cfgSig.cols);
120 
121  int total = _cfgSig.rows*_cfgSig.cols;
122 
123 
124  string RGB_name[] = { "_B","_G","_R" };
125  double *phasedata = new double[total];
126 
127  for (int z = 0; z < (bitpixel / 8); z++)
128  {
129  ifstream fphase(strtok((char*)fnamestr.c_str(), ".") + RGB_name[z] + "bin", ifstream::binary);
130 
131  fphase.read(reinterpret_cast<char*>(phasedata), sizeof(double) * total);
132 
133  for (int col = 0; col < _cfgSig.cols; col++)
134  {
135  for (int row = 0; row < _cfgSig.rows; row++)
136  {
137  phaseMat[z](row, col) = phasedata[_cfgSig.rows*col + row];
138  }
139  }
140  fphase.close();
141  }
142  delete[] phasedata; */
143  }
144  }
145  else
146  {
147  LOG("wrong type\n");
148  }
149 
152  //nomalization
153  Nr = _cfgSig.rows;
154  Nc = _cfgSig.cols;
155  PhaseOriginal.resize(Nr, Nc);
156  PhaseUnwrapped.resize(Nr, Nc);
157  for (int i = 0; i < _cfgSig.rows; i++)
158  {
159  for (int j = 0; j < _cfgSig.cols; j++)
160  {
161  PhaseOriginal(i, j) = phaseMat(i, j) / 255.0*M_PI*2 - M_PI;
162  PhaseUnwrapped(i, j) = 0;
163  }
164  }
165 
166  LOG("data nomalization complete\n");
167 
168  return true;
169 }
170 
172 {
173  return false;
174 }
175 
176 bool ophSigPU::runPU(void)
177 {
178  auto start_time = CUR_TIME;
179 
180  matrix<Real> residue(Nr,Nc);
181  residue.zeros();
182  matrix<Real> branchcut(Nr, Nc);
183  branchcut.zeros();
184  phaseResidues(residue);
185  branchCuts(residue, branchcut);
186  floodFill(branchcut);
187 
188  auto end_time = CUR_TIME;
189 
190  auto during_time = ((std::chrono::duration<Real>)(end_time - start_time)).count();
191 
192  LOG("Implement time : %.5lf sec\n", during_time);
193 
194  return false;
195 }
196 
197 bool ophSigPU::savePhaseUnwrapped(const char * fname)
198 {
199  oph::uchar* phaseData;
200  phaseData = (oph::uchar*)malloc(sizeof(oph::uchar) * Nr * Nc);
201 
202  string fnamestr = fname;
203 
204  double maxPhase = 0;
205  double minPhase = 0;
206  for (int i = 0; i < Nr; i++)
207  {
208  for (int j = 0; j < Nc; j++)
209  {
210  if (PhaseUnwrapped(i, j) > maxPhase)
211  {
212  maxPhase = PhaseUnwrapped(i, j);
213  }
214  if (PhaseUnwrapped(i, j) < minPhase)
215  {
216  minPhase = PhaseUnwrapped(i, j);
217  }
218  }
219  }
220 
221 
222  for (int i = 0; i < Nr; i++)
223  {
224  for (int j = 0; j < Nc; j++)
225  {
226  phaseData[i*Nc + j] = (uchar)(((PhaseUnwrapped(Nr - i -1,j)-minPhase)/(maxPhase-minPhase))*255.0);
227  }
228  }
229  saveAsImg(fnamestr.c_str(), 8, phaseData, Nc, Nr);
230  return TRUE;
231 }
232 
233 
234 bool ophSigPU::readConfig(const char * fname)
235 {
236  return false;
237 }
238 
239 void ophSigPU::phaseResidues(matrix<Real>& outputResidue)
240 {
241  /*
242  This code was written with slight modification based on
243  https://kr.mathworks.com/matlabcentral/fileexchange/22504-2d-phase-unwrapping-algorithms
244  */
245 
246  matrix<Real> below(Nr, Nc);
247  matrix<Real> right(Nr, Nc);
248  matrix<Real> belowright(Nr, Nc);
249  below.zeros(); right.zeros(); belowright.zeros();
250 
251  for (int i = 0; i < Nr-1; i++)
252  {
253  for (int j = 0; j < Nc-1; j++)
254  {
255  below(i, j) = PhaseOriginal(i + 1, j);
256  right(i, j) = PhaseOriginal(i, j + 1);
257  belowright(i, j) = PhaseOriginal(i + 1, j + 1);
258  }
259  }
260  for (int i = 0; i < Nr-1; i++)
261  {
262  below(i, Nc-1) = PhaseOriginal(i + 1, Nc-1);
263  }
264  for (int j = 0; j < Nc - 1; j++)
265  {
266  right(Nr-1, j) = PhaseOriginal(Nr-1, j+1);
267  }
268 
269  double res1, res2, res3, res4;
270  double temp_residue;
271  for (int i = 0; i < Nr; i++)
272  {
273  for (int j = 0; j < Nc; j++)
274  {
275  res1 = mod2pi(PhaseOriginal(i, j) - below(i, j));
276  res2 = mod2pi(below(i, j) - belowright(i,j));
277  res3 = mod2pi(belowright(i, j) - right(i, j));
278  res4 = mod2pi(right(i, j) - PhaseOriginal(i, j));
279 
280  temp_residue = res1 + res2 + res3 + res4;
281  if (temp_residue >= 6.)
282  {
283  outputResidue(i, j) = 1;
284  }
285  else if (temp_residue <= -6.)
286  {
287  outputResidue(i, j) = -1;
288  }
289  else
290  {
291  outputResidue(i, j) = 0;
292  }
293  }
294  }
295 
296  for (int i = 0; i < Nr; i++)
297  {
298  outputResidue(i, 0) = 0;
299  outputResidue(i, Nc - 1) = 0;
300  }
301  for (int j = 0; j < Nc; j++)
302  {
303  outputResidue(0, j) = 0;
304  outputResidue(Nr - 1, j) = 0;
305  }
306 }
307 
308 void ophSigPU::branchCuts(matrix<Real>& inputResidue, matrix<Real>& outputBranchCuts)
309 {
310  /*
311  This code was written with slight modification based on
312  https://kr.mathworks.com/matlabcentral/fileexchange/22504-2d-phase-unwrapping-algorithms
313  */
314 
315  int clusterCounter = 1;
316  int satelliteResidue = 0;
317  matrix<int> residueBinary(Nr, Nc);
318  for (int i = 0; i < Nr; i++)
319  {
320  for (int j = 0; j < Nc; j++)
321  {
322  if (inputResidue(i, j) != 0)
323  {
324  residueBinary(i, j) = 1;
325  }
326  else
327  {
328  residueBinary(i, j) = 0;
329  }
330  }
331  }
332  matrix<Real> residueBalanced(Nr, Nc);
333  residueBalanced.zeros();
334  matrix<int> adjacentResidue(Nr, Nc);
335  adjacentResidue.zeros();
336  int missedResidue = 0;
337 
338  int adjacentResidueCount = 0;
339 
340  for (int i = 0; i < Nr; i++)
341  {
342  for (int j = 0; j < Nc; j++)
343  {
344  if (residueBinary(i, j) == 1)
345  {
346  int rActive = i;
347  int cActive = j;
348 
349  int radius = 1;
350  int countNearbyResidueFlag = 1;
351  clusterCounter = 1;
352  adjacentResidue.zeros();
353  int chargeCounter = inputResidue(rActive, cActive);
354  if (residueBalanced(rActive, cActive) != 1)
355  {
356  while (chargeCounter != 0)
357  {
358  for (int m = rActive - radius; m < rActive + radius + 1; m++)
359  {
360  for (int n = cActive - radius; n < cActive + radius + 1; n++)
361  {
362  if (((abs(m - rActive) == radius) | (abs(n - cActive) == radius)) & (chargeCounter != 0))
363  {
364  if ((m < 1) | (m >= Nr - 1) | (n < 1) | (m >= Nc - 1))
365  {
366  if (m >= Nr - 1) { m = Nr - 1; }
367  if (n >= Nc - 1) { n = Nc - 1; }
368  if (n < 1) { n = 0; }
369  if (m < 1) { m = 0; }
370  placeBranchCutsInternal(outputBranchCuts, rActive, cActive, m, n);
371  clusterCounter += 1;
372  chargeCounter = 0;
373  residueBalanced(rActive, cActive) = 1;
374  }
375  if (residueBinary(m, n))
376  {
377  if (countNearbyResidueFlag == 1) { adjacentResidue(m, n) = 1; }
378  placeBranchCutsInternal(outputBranchCuts, rActive, cActive, m, n);
379  clusterCounter += 1;
380  if (residueBalanced(m, n) == 0)
381  {
382  residueBalanced(m, n) = 1;
383  chargeCounter += inputResidue(m, n);
384  }
385  if (chargeCounter == 0) { residueBalanced(rActive, cActive) = 1; }
386  }
387 
388  }
389  }
390  }
391 
392  double sumAdjacentResidue = 0;
393  int adjacentSize = 0;
394  for (int ii = 0; ii < Nr; ii++)
395  {
396  for (int jj = 0; jj < Nc; jj++)
397  {
398  sumAdjacentResidue += adjacentResidue(ii, jj);
399  }
400  }
401  if (sumAdjacentResidue == 0)
402  {
403  radius += 1;
404  rActive = i;
405  cActive = j;
406  }
407  else
408  {
409  vector<int> rAdjacent, cAdjacent;
410  if (countNearbyResidueFlag == 1)
411  {
412  findNZ(adjacentResidue, rAdjacent, cAdjacent);
413  adjacentSize = rAdjacent.size();
414  rActive = rAdjacent[0];
415  cActive = cAdjacent[0];
416  adjacentResidueCount = 1;
417  residueBalanced(rActive, cActive) = 1;
418  countNearbyResidueFlag = 0;
419  }
420  else
421  {
422  adjacentResidueCount += 1;
423  if (adjacentResidueCount <= adjacentSize)
424  {
425  rActive = rAdjacent[adjacentResidueCount-1];
426  cActive = cAdjacent[adjacentResidueCount-1];
427  }
428  else
429  {
430  radius += 1;
431  rActive = i;
432  cActive = j;
433  adjacentResidue.zeros();
434  countNearbyResidueFlag = 1;
435  }
436  }
437  }
438  if (radius > MaxBoxRadius)
439  {
440  if (clusterCounter != 1)
441  {
442  missedResidue += 1;
443  }
444  else
445  {
446  satelliteResidue += 1;
447  }
448  chargeCounter = 0;
449  while (clusterCounter == 1)
450  {
451  rActive = i;
452  cActive = j;
453  for (int m = rActive - radius; m < rActive + radius + 1; m++)
454  {
455  for (int n = cActive - radius; n < cActive + radius + 1; n++)
456  {
457  if (((abs(m - rActive) == radius) | (abs(n - cActive) == radius) ))
458  {
459  if ((m < 1) | (m >= Nr - 1) | (n < 1) | (m >= Nc - 1))
460  {
461  if (m >= Nr - 1) { m = Nr - 1; }
462  if (n >= Nc - 1) { n = Nc - 1; }
463  if (n < 1) { n = 0; }
464  if (m < 1) { m = 0; }
465  placeBranchCutsInternal(outputBranchCuts, rActive, cActive, m, n);
466  clusterCounter += 1;
467  residueBalanced(rActive, cActive) = 1;
468  }
469  if (residueBinary(m, n))
470  {
471  placeBranchCutsInternal(outputBranchCuts, rActive, cActive, m, n);
472  clusterCounter += 1;
473  residueBalanced(rActive, cActive) = 1;
474  }
475 
476  }
477  }
478  }
479  radius += 1;
480  }
481  }
482  }
483  }
484  }
485  }
486  }
487 
488  LOG("Branch cut operation completed. \n");
489  LOG("Satellite residues accounted for = %d\n", satelliteResidue);
490 }
491 
492 void ophSigPU::placeBranchCutsInternal(matrix<Real>& branchCuts, int r1, int c1, int r2, int c2)
493 {
494  branchCuts(r1, c1) = 1;
495  branchCuts(r2, c2) = 1;
496  double rdist = abs(r2 - r1);
497  double cdist = abs(c2 - c1);
498  int rsign = (r2 > r1) ? 1 : -1;
499  int csign = (c2 > c1) ? 1 : -1;
500  int r = 0;
501  int c = 0;
502  if (rdist > cdist)
503  {
504  for (int i = 0; i <= rdist; i++)
505  {
506  r = r1 + i*rsign;
507  c = c1 + (int)(round(((double)(i))*((double)(csign))*cdist / rdist));
508  branchCuts(r, c) = 1;
509  }
510  }
511  else
512  {
513  for (int j = 0; j <= cdist; j++)
514  {
515  c = c1 + j*csign;
516  r = r1 + (int)(round(((double)(j))*((double)(rsign))*rdist / cdist));
517  branchCuts(r, c) = 1;
518  }
519  }
520 }
521 
522 void ophSigPU::floodFill(matrix<Real>& inputBranchCuts)
523 {
524  /*
525  This code was written with slight modification based on
526  https://kr.mathworks.com/matlabcentral/fileexchange/22504-2d-phase-unwrapping-algorithms
527  */
528 
529  matrix<int> adjoin(Nr, Nc);
530  adjoin.zeros();
531 
532  matrix<int> unwrappedBinary(Nr, Nc);
533  unwrappedBinary.zeros();
534 
535  // set ref phase
536  bool refPhaseFound = false;
537  for (int i = 1; (i < Nr - 1) & !refPhaseFound; i++)
538  {
539  for (int j = 1; (j < Nc - 1) & !refPhaseFound; j++)
540  {
541  if (inputBranchCuts(i, j) == 0)
542  {
543  adjoin(i - 1, j) = 1;
544  adjoin(i + 1, j) = 1;
545  adjoin(i, j - 1) = 1;
546  adjoin(i, j + 1) = 1;
547  unwrappedBinary(i, j) = 1;
548  PhaseUnwrapped(i, j) = PhaseOriginal(i, j);
549  refPhaseFound = true;
550  }
551  }
552  }
553 
554  // floodfill
555  int countLimit = 0;
556  int adjoinStuck = 0;
557  vector<int> rAdjoin;
558  vector<int> cAdjoin;
559  int rActive = 0;
560  int cActive = 0;
561  double phaseRef = 0;
562  while ((matrixPartialSum(adjoin, 1, 1, Nr - 2, Nc - 2) > 0) & (countLimit < 100))
563  {
564  //while (countLimit < 100)
565  //{
566  LOG("%d\n", matrixPartialSum(adjoin, 1, 1, Nr - 2, Nc - 2));
567  findNZ(adjoin, rAdjoin, cAdjoin);
568  if (adjoinStuck == rAdjoin.size())
569  {
570  countLimit += 1;
571  LOG("countLimit %d\n", countLimit);
572  }
573  else
574  {
575  countLimit = 0;
576  }
577  adjoinStuck = rAdjoin.size();
578  for (int i = 0; i < rAdjoin.size(); i++)
579  {
580  rActive = rAdjoin[i];
581  cActive = cAdjoin[i];
582  if ((rActive <= Nr - 2) & (rActive >= 1) & (cActive <= Nc - 2) & (cActive >= 1))
583  {
584  if ((inputBranchCuts(rActive + 1, cActive) == 0) & (unwrappedBinary(rActive + 1, cActive) == 1))
585  {
586  phaseRef = PhaseUnwrapped(rActive + 1, cActive);
587  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
588  unwrappedBinary(rActive, cActive) = 1;
589  adjoin(rActive, cActive) = 0;
590  if ((unwrappedBinary(rActive - 1, cActive) == 0) & (inputBranchCuts(rActive - 1, cActive) == 0))
591  {
592  adjoin(rActive - 1, cActive) = 1;
593  }
594  if ((unwrappedBinary(rActive, cActive-1) == 0) & (inputBranchCuts(rActive, cActive-1) == 0))
595  {
596  adjoin(rActive, cActive-1) = 1;
597  }
598  if ((unwrappedBinary(rActive, cActive+1) == 0) & (inputBranchCuts(rActive, cActive+1) == 0))
599  {
600  adjoin(rActive, cActive+1) = 1;
601  }
602  }
603  if ((inputBranchCuts(rActive - 1, cActive) == 0) & (unwrappedBinary(rActive - 1, cActive) == 1))
604  {
605  phaseRef = PhaseUnwrapped(rActive - 1, cActive);
606  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
607  unwrappedBinary(rActive, cActive) = 1;
608  adjoin(rActive, cActive) = 0;
609  if ((unwrappedBinary(rActive + 1, cActive) == 0) & (inputBranchCuts(rActive + 1, cActive) == 0))
610  {
611  adjoin(rActive + 1, cActive) = 1;
612  }
613  if ((unwrappedBinary(rActive, cActive - 1) == 0) & (inputBranchCuts(rActive, cActive - 1) == 0))
614  {
615  adjoin(rActive, cActive - 1) = 1;
616  }
617  if ((unwrappedBinary(rActive, cActive + 1) == 0) & (inputBranchCuts(rActive, cActive + 1) == 0))
618  {
619  adjoin(rActive, cActive + 1) = 1;
620  }
621  }
622  if ((inputBranchCuts(rActive, cActive +1) == 0) & (unwrappedBinary(rActive, cActive+1) == 1))
623  {
624  phaseRef = PhaseUnwrapped(rActive, cActive+1);
625  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
626  unwrappedBinary(rActive, cActive) = 1;
627  adjoin(rActive, cActive) = 0;
628  if ((unwrappedBinary(rActive + 1, cActive) == 0) & (inputBranchCuts(rActive + 1, cActive) == 0))
629  {
630  adjoin(rActive + 1, cActive) = 1;
631  }
632  if ((unwrappedBinary(rActive, cActive - 1) == 0) & (inputBranchCuts(rActive, cActive - 1) == 0))
633  {
634  adjoin(rActive, cActive - 1) = 1;
635  }
636  if ((unwrappedBinary(rActive -1 , cActive) == 0) & (inputBranchCuts(rActive-1, cActive) == 0))
637  {
638  adjoin(rActive-1, cActive) = 1;
639  }
640  }
641  if ((inputBranchCuts(rActive, cActive-1) == 0) & (unwrappedBinary(rActive, cActive-1) == 1))
642  {
643  phaseRef = PhaseUnwrapped(rActive, cActive-1);
644  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
645  unwrappedBinary(rActive, cActive) = 1;
646  adjoin(rActive, cActive) = 0;
647  if ((unwrappedBinary(rActive + 1, cActive) == 0) & (inputBranchCuts(rActive + 1, cActive) == 0))
648  {
649  adjoin(rActive + 1, cActive) = 1;
650  }
651  if ((unwrappedBinary(rActive -1, cActive) == 0) & (inputBranchCuts(rActive-1, cActive) == 0))
652  {
653  adjoin(rActive-1, cActive) = 1;
654  }
655  if ((unwrappedBinary(rActive, cActive + 1) == 0) & (inputBranchCuts(rActive, cActive + 1) == 0))
656  {
657  adjoin(rActive, cActive + 1) = 1;
658  }
659  }
660 
661  }
662 
663  }
664  //}
665  }
666 
667  adjoin.zeros();
668  for (int i = 1; i <= Nr - 2; i++)
669  {
670  for (int j = 1; j < Nc - 2; j++)
671  {
672  if ((inputBranchCuts(i, j) == 1) &
673  ((inputBranchCuts(i + 1, j) == 0) |
674  (inputBranchCuts(i - 1, j) == 0) |
675  (inputBranchCuts(i, j - 1) == 0) |
676  (inputBranchCuts(i, j + 1) == 0)))
677  {
678  adjoin(i, j) = 1;
679  }
680  }
681  }
682  findNZ(adjoin, rAdjoin, cAdjoin);
683  for (int i = 0; i < rAdjoin.size(); i++)
684  {
685  rActive = rAdjoin[i];
686  cActive = cAdjoin[i];
687 
688  if (unwrappedBinary(rActive + 1, cActive) == 1)
689  {
690  phaseRef = PhaseUnwrapped(rActive + 1, cActive);
691  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
692  unwrappedBinary(rActive, cActive) = 1;
693  adjoin(rActive, cActive) = 0;
694  }
695  if (unwrappedBinary(rActive - 1, cActive) == 1)
696  {
697  phaseRef = PhaseUnwrapped(rActive - 1, cActive);
698  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
699  unwrappedBinary(rActive, cActive) = 1;
700  adjoin(rActive, cActive) = 0;
701  }
702  if (unwrappedBinary(rActive, cActive+1) == 1)
703  {
704  phaseRef = PhaseUnwrapped(rActive, cActive+1);
705  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
706  unwrappedBinary(rActive, cActive) = 1;
707  adjoin(rActive, cActive) = 0;
708  }
709  if (unwrappedBinary(rActive, cActive-1) == 1)
710  {
711  phaseRef = PhaseUnwrapped(rActive, cActive-1);
712  PhaseUnwrapped(rActive, cActive) = unwrap(phaseRef, PhaseOriginal(rActive, cActive));
713  unwrappedBinary(rActive, cActive) = 1;
714  adjoin(rActive, cActive) = 0;
715  }
716  }
717  LOG("Floodfill completed\n");
718 }
719 
720 double ophSigPU::unwrap(double phaseRef, double phaseInput)
721 {
722  double diff = phaseInput - phaseRef;
723  double modval = mod2pi(diff);
724  return phaseRef + modval;
725 }
726 
727 double ophSigPU::mod2pi(double phase)
728 {
729  double temp;
730  temp = fmod(phase, 2 * M_PI);
731  if (temp > M_PI)
732  {
733  temp = temp - 2 * M_PI;
734  }
735  if (temp < -M_PI)
736  {
737  temp = temp + 2 * M_PI;
738  }
739  return temp;
740 }
741 
742 void ophSigPU::findNZ(matrix<int>& inputMatrix, vector<int>& row, vector<int>& col)
743 {
744  row.clear();
745  col.clear();
746  for (int i = 0; i < inputMatrix.size(_X); i++)
747  {
748  for (int j = 0; j < inputMatrix.size(_Y); j++)
749  {
750  if (inputMatrix(i,j) != 0)
751  {
752  row.push_back(i);
753  col.push_back(j);
754  }
755  }
756  }
757 }
758 
759 int ophSigPU::matrixPartialSum(matrix<int>& inputMatrix, int r1, int c1, int r2, int c2)
760 {
761  int outputsum = 0;
762  for (int i = r1; i <= r2; i++)
763  {
764  for (int j = c1; j <= c2; j++)
765  {
766  outputsum += inputMatrix(i, j);
767  }
768  }
769  return outputsum;
770 }
void abs(const oph::Complex< T > &src, oph::Complex< T > &dst)
Definition: function.h:112
matrix< Real > PhaseOriginal
Definition: ophSigPU.h:81
void floodFill(matrix< Real > &inputBranchCuts)
Definition: ophSigPU.cpp:522
virtual bool saveAsImg(const char *fname, uint8_t bitsperpixel, uchar *src, int width, int height)
Function for creating image files.
Definition: Openholo.cpp:90
double mod2pi(double phase)
Definition: ophSigPU.cpp:727
unsigned char uchar
Definition: typedef.h:64
void findNZ(matrix< int > &inputMatrix, vector< int > &row, vector< int > &col)
Definition: ophSigPU.cpp:742
int Nr
Definition: ophSigPU.h:79
int MaxBoxRadius
Definition: ophSigPU.h:78
#define CUR_TIME
Definition: function.h:58
ophSigPU(void)
Definition: ophSigPU.cpp:3
int Nc
Definition: ophSigPU.h:80
#define _Y
Definition: define.h:84
ophSigConfig _cfgSig
Definition: ophSig.h:108
void placeBranchCutsInternal(matrix< Real > &branchCuts, int r1, int c1, int r2, int c2)
Definition: ophSigPU.cpp:492
#define _X
Definition: define.h:80
double unwrap(double phaseRef, double phaseInput)
Definition: ophSigPU.cpp:720
matrix< Real > PhaseUnwrapped
Definition: ophSigPU.h:82
Definition: struct.h:69
int cols
Definition: ophSig.h:63
int rows
Definition: ophSig.h:64
void phaseResidues(matrix< Real > &outputResidue)
Definition: ophSigPU.cpp:239
uint16_t bitsperpixel
Definition: struct.h:61
bool savePhaseUnwrapped(const char *fname)
Save the unwrapped phase data to image file.
Definition: ophSigPU.cpp:197
bool runPU(void)
Run phase unwrapping algorithm.
Definition: ophSigPU.cpp:176
int matrixPartialSum(matrix< int > &inputMatrix, int r1, int c1, int r2, int c2)
Definition: ophSigPU.cpp:759
bool setPUparam(int maxBoxRadius)
Set parameters for Goldstein branchcut algorithm.
Definition: ophSigPU.cpp:7
void branchCuts(matrix< Real > &inputResidue, matrix< Real > &outputBranchCuts)
Definition: ophSigPU.cpp:308
uint32_t width
Definition: struct.h:58
uint32_t height
Definition: struct.h:59
bool readConfig(const char *fname)
Read configure file.
Definition: ophSigPU.cpp:234
uint8_t signature[2]
Definition: struct.h:51
bool loadPhaseOriginal(void)
Definition: ophSigPU.cpp:171
#define M_PI
Definition: define.h:52