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