OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 13 — Jun. 18, 2012
  • pp: 14075–14089
« Show journal navigation

Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising

Howard Y. H. Huang, L. Tian, Z. Zhang, Y. Liu, Z. Chen, and G. Barbastathis  »View Author Affiliations


Optics Express, Vol. 20, Issue 13, pp. 14075-14089 (2012)
http://dx.doi.org/10.1364/OE.20.014075


View Full Text Article

Acrobat PDF (2265 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Phase unwrapping is a challenging task for interferometry based techniques in the presence of noise. The majority of existing phase unwrapping techniques are path-following methods, which explicitly or implicitly define an intelligent path and integrate phase difference along the path to mitigate the effect of erroneous pixels. In this paper, a path-independent unwrapping method is proposed where the unwrapped phase gradient is determined from the wrapped phase and subsequently denoised by a TV minimization based method. Unlike the wrapped phase map where 2π phase jumps are present, the gradient of the unwrapped phase map is smooth and slowly-varying at noise-free areas. On the other hand, the noise is greatly amplified by the differentiation process, which makes it easier to separate from the smooth phase gradient. Thus an approximate unwrapped phase can be obtained by integrating the denoised phase gradient. The final unwrapped phase map is subsequently determined by adding the first few modes of the unwrapped phase. The proposed method is most suitable for unwrapping phase maps without abrupt phase changes. Its capability has been demonstrated both numerically and by experimental data from shearography and electronic speckle pattern interferometry (ESPI).

© 2012 OSA

1. Introduction

Phase unwrapping is an essential step for quantitative measurement in a wide variety of interferometry based techniques such as speckle interferometry [1

1. Y. H. Huang, Y. S. Liu, S. Y. Hung, C. G. Li, and F. Janabi-Sharifi, “Dynamic phase evaluation in sparse-sampled temporal speckle pattern sequence,” Opt. Lett. 36(4), 526–528 (2011). [CrossRef] [PubMed]

], electron holography [2

2. E. Volkl, L. F. Allard, and D. C. Joy, eds., Introduction to Electron holography (Plenum, 1999).

], interferometric synthetic aperture radar (InSAR) [3

3. B. Osmanoglu, T. H. Dixon, S. Wdowinski, and E. Cabral-Cano, “On the importance of path for phase unwrapping in synthetic aperture radar interferometry,” Appl. Opt. 50(19), 3205–3220 (2011). [CrossRef] [PubMed]

], magnetic resonance imaging (MRI) [4

4. J. Langley and Q. Zhao, “Unwrapping magnetic resonance phase maps with Chebyshev polynomials,” Magn. Reson. Imaging 27(9), 1293–1301 (2009). [CrossRef] [PubMed]

], fringe projection [5

5. C. J. Tay, C. Quan, T. Wu, and Y. H. Huang, “Integrated method for 3-D rigid-body displacement measurement using fringe projection,” Opt. Eng. 43(5), 1152–1159 (2004). [CrossRef]

], etc. In the noise-free case, the unwrapping process is as simple as conducting a raster scan and adding an integral multiple of 2πto each pixel in the wrapped phase map as necessary to confine the phase change of adjacent pixels to less thanπ. However, real experiment data always contain noise, which considerably inhibits the phase unwrapping task [6

6. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software, (John Wiley & Sons, 1998).

].

If the wrapped phase map contains noise, the unwrapping problem becomes ill-posed: different unwrapping paths result in different unwrapped phase maps. A well-chosen path tends to go first through pixels unaffected by noise and then tackle pixels with noise and phase residues. Numerous path-following algorithms have been developed during the last two decades [7

7. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. 23(4), 713–720 (1988). [CrossRef]

13

13. M. Gdeisat, M. Arevalillo-Herráez, D. Burton, and F. Lilley, “Three-dimensional phase unwrapping using the Hungarian algorithm,” Opt. Lett. 34(19), 2994–2996 (2009). [CrossRef] [PubMed]

] and widely applied in various fields. More details on typical path-following methods will be given in section 2.

The path-following algorithms mentioned above are effective. However, the choice for best path is not always obvious in applications, especially when the noise statistics are unknown or, worse, in the case of object-dependent noise. To mitigate this problem, path-independent algorithms have been proposed. An intuitive approach for path-independent unwrapping is to approximate the unwrapped phase map from among a set of smooth basis functions and minimize the difference between the wrapped phase map and the proposed unwrapped phase [14

14. S. S. Gorthi, G. Rajshekhar, and P. Rastogi, “Strain estimation in digital holographic interferometry using piecewise polynomial phase approximation based method,” Opt. Express 18(2), 560–565 (2010). [CrossRef] [PubMed]

,15

15. J. Bioucas-Dias, V. Katkovnik, J. Astola, and K. Egiazarian, “Absolute phase estimation: adaptive local denoising and global unwrapping,” Appl. Opt. 47(29), 5358–5369 (2008). [CrossRef] [PubMed]

]. Then, however, the success of the unwrapping method depends on the degree to which the unknown surface can be well fitted by the chosen set of basis functions.

Another, more systematic approach is to minimize the difference between the gradient of the wrapped phase and the gradient of the unwrapped phase according to a Lpnorm. This turns the problem of phase unwrapping into solving a discretized partial differential equation [6

6. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software, (John Wiley & Sons, 1998).

,16

16. D. C. Ghiglia and L. A. Romero, “Direct phase estimation from phase differences using fast elliptic partial differential equation solvers,” Opt. Lett. 14(20), 1107–1109 (1989). [CrossRef] [PubMed]

,17

17. M. A. Schofield and Y. M. Zhu, “Fast phase unwrapping algorithm for interferometric applications,” Opt. Lett. 28(14), 1194–1196 (2003). [CrossRef] [PubMed]

], in particular a Poisson equation when L2 norm is employed. Methods in this category seem to underestimate the magnitude of the resultant unwrapped phase due to the strong smoothing effect of Poisson solvers [6

6. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software, (John Wiley & Sons, 1998).

], thus post-processing techniques [18

18. R. Yamaki and A. Hirose, “Singularity-spreading phase unwrapping,” IEEE Trans. Geosci. Rem. Sens. 45(10), 3240–3251 (2007). [CrossRef]

,19

19. S. Tomioka, S. Heshmat, N. Miyamoto, and S. Nishiyama, “Phase unwrapping for noisy phase maps using rotational compensator with virtual singular points,” Appl. Opt. 49(25), 4735–4745 (2010). [CrossRef] [PubMed]

] including the congruence process [20

20. M. D. Pritt, “Congruence in least-squares phase unwrapping,” in Proceedings Vol. II: Remote Sensing - a Scientific Vision for Sustainable Development, IGARSS '97 - 1997 International Geoscience and Remote Sensing Symposium, 1997), pp.875–877.

] have been applied to improve the algorithms. As a result, very good unwrapping results [21

21. J. Strand and T. Taxt, “Performance evaluation of two-dimensional phase unwrapping algorithms,” Appl. Opt. 38(20), 4333–4344 (1999). [CrossRef] [PubMed]

23

23. J. Parkhurst, G. Price, P. Sharrock, and C. Moore, “Phase unwrapping algorithms for use in a true real-time optical body sensor system for use during radiotherapy,” Appl. Opt. 50(35), 6430–6439 (2011). [CrossRef] [PubMed]

] have been obtained.

The first approximation (1st mode) of the unwrapped phase map is obtained by a simple integration of the denoised phase gradient. To correct for additional errors due to the smoothing effect of TV minimization we subtract the wrapped phase map from the 1st mode unwrapped phase map to determine the residual wrapped phase map and iteratively apply the proposed method to retrieve the 2nd, 3rd modes of the unwrapped phase maps until no apparent 2π jumps are detected in the residual wrapped phase map. In our experiments, 1-2 iterations were typically sufficient to determine the final correct unwrapped phase map.

2. Problem description and path-following solutions

Phase extraction methods normally retrieve a wrapped phase ϕw from an extended arctangent function which confines the wrapped phase within(π,π]and makes it different from the genuine unwrapped phase ϕ representing the physical quantity by an integer multiple of2πas:
ϕw=ϕ2nπ         nZ, ϕw(π,π].
(1)
The task of phase unwrapping is thus to add a correct 2nπ to each pixel of the wrapped phase such that a continuous unwrapped phase map is reconstructed.

The presence of noise introduces phase residues in the wrapped phase map; that causes closed-loop line integrals to assume non-zero values, which can be interpreted as poles (singularities) appearing in the reconstructed function. Thus, the unwrapped phase values of phase residues cannot be uniquely determined and the error caused by the residues propagates to adjacent pixels. For this reason, path-following phase unwrapping algorithms tend to define branch-cuts to be avoided by the path integrals so as to first unwrap intact pixels while keeping residues until the last stage. Thus, they can mitigate error propagation.

Goldstein’s branch cut algorithm [7

7. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. 23(4), 713–720 (1988). [CrossRef]

] is such an approach. After creating branch-cuts for all the poles and connecting unbalanced poles to the nearest boundary, the unwrapping process then takes a path that does not cross the branch-cuts. This method works well for wrapped phase with sparsely distributed residues, however, when residues are clustered at certain areas, the branch-cuts may be incorrectly produced and the unwrapping result may be problematic. Another type of algorithms introduces an auxiliary quality map to guide the unwrapping path from high quality pixels to low quality ones [8

8. Y. Xu and C. Ai, ““Simple and effective phase unwrapping technique,” Interferometry IV: Techniques and Analysis,” Proc. SPIE 2003, 254–263 (1993). [CrossRef]

]. The quality maps can be generated from the wrapped phase map itself (for example by calculating the phase gradient and assigning high gradient pixels as low quality) or calculated from the original raw intensity data (for example the contrast or correlation coefficient of the initial data). Mask cut methods [9

9. C. Prati, M. Giani, and N. Leuratti, “SAR interferometry: a 2-D phase unwrapping technique based on phase and absolute values information,” in Proceedings of the 1990 International Geoscience and Remote Sensing Symposium (IEEE, 1990), pp. 2043–2046.

] try to combine the branch cut method with the quality-guided method for path generation. Contrary to the above mentioned explicit path generating method, Flynn [10

10. T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14(10), 2692–2702 (1997). [CrossRef]

] proposed a minimum discontinuity method which uses graph theory to implicitly detect existing edges of 2πphase jumps (“jump cuts”) and extend them to form close loops. By iteratively adding or subtracting 2πto one of the part separated by the loop, the phase jumps number is gradually minimized, and an unwrapped phase is finally generated. Flynn’s algorithm is most suitable for unwrapping smooth phase maps as the algorithm tends to minimize the discontinuity in the unwrapped result.

3. Phase unwrapping using TV prior

In this section, we describe our algorithm for path-independent phase unwrapping using total-variation (TV) as the denoising prior. Our motivation is similar to Flynn’s idea of minimizing discontinuities in the unwrapped phase map, except instead we minimize the total-variation of the phase derivative map using the TV solver. However, unlike Flynn’s algorithm which implicitly generates a series of jump cuts, our algorithm works pixel-wise without any path-dependent process. After obtaining the TV minimized phase gradient map, an approximate unwrapped phase map can be determined by a simple integration. Then we apply an iterative process to correct for possible filtering effects by the TV solver. The initial TV denoising as well as iterative correction algorithms are detailed step-by-step in the following sub-sections.

3.1 Definition of the wrapping function

A simple wrapping function which wraps a phase map into the range of (π,π] by adding an integer multiple of 2π to each pixel can be defined as:
W{ψ}={2arctan[sin(ψ)1+cos(ψ)]          cos(ψ)1             π                          cos(ψ)=1  
(2)
where arctan[] is a normal arctangent function with range (π/2,  π/2)and the wrapping function W{} has an extended range of (π, π]. This wrapping function enables comparison between the wrapped and unwrapped phase maps and will be frequently used in the following derivations.

3.2 Determination of the phase gradient

Instead of trying to determine the correct unwrapped phase, we try to determine the correct gradient of the unwrapped phase and then follow up with an integration step to recover the final unwrapped phase. Noting that the only difference between the wrapped phase and the unwrapped phase is the 2π jumps, the discretized unwrapped phase derivatives can then be obtained from the wrapped phase map as

ϕx=ϕx=W{ϕwx}
(3)
ϕy=ϕy=W{ϕwy}.
(4)

The wrapped phase derivatives obtained in this way have two merits. The first one is that the 2πjumps due to wrapping have been removed by the wrapping function W{}, making the derivative map admissible to noise-removal filters. The second merit is that the differentiation process amplifies the noise. Intuition usually suggests that such action is to be avoided; however, the core concept of our denoising is that the amplified noise will be easier to be picked up and removed by the TV solver. Thus, surprisingly, we can turn noise amplification to our advantage.

3.3 TV minimization for denoising

As we already mentioned, usually amplified noise is considered deadly because it makes the job of denoising filters even tougher. However, that is not true for total variation (TV) minimization [19

19. S. Tomioka, S. Heshmat, N. Miyamoto, and S. Nishiyama, “Phase unwrapping for noisy phase maps using rotational compensator with virtual singular points,” Appl. Opt. 49(25), 4735–4745 (2010). [CrossRef] [PubMed]

] as we apply it here. Using TV, we apply a smoothness prior on the phase derivative ϕx and subsequently calculate a denoised estimate ϕ^x according to
ϕ^x= arg min     ϕ^x {||ϕ^xϕx||2/(2μ)+TV(ϕ^x)}
(5)
where ϕxis the input noisy phase derivative, ||||is the Euclidean norm, μis a regularization parameter, and TV() is a functional that determines the total variation of its argument. If we define the TV functional asTV(ψ)=||ψ||, i.e. the sum of the gradient magnitude over all pixels, it has been proved [29

29. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis. 20(1/2), 89–97 (2004). [CrossRef]

] that the solution to the minimization problem (5) is given by
ϕ^x=ϕxπμκ(ϕx)
(6)
where πμκ(ϕx)is the nonlinear projection of ϕx. Computing the nonlinear projection πμκ(ϕx)is identical to solving the following problem
arg min     p {||μpϕx||2:  |p|21}
(7)
and πμκ(ϕx)can then be determined as

πμκ(ϕx)=μ p.
(8)

In [29

29. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis. 20(1/2), 89–97 (2004). [CrossRef]

], Chambolle proposed a semi-implicit gradient descent algorithm to determine pas the convergent value of a vector sequence pn of
pn+1=pn+τ[(pnϕx/μ)]1+τ||(pnϕx/μ)||
(9)
where the iteration is initiated at p1=(0, 0)and the parameter τis routinely set to τ=0.25.

In the TV minimization algorithm proposed above, the regularization parameter μand a stopping criterion for the iteration in Eq. (9) still need to be defined. We found that the quality of reconstruction is quite tolerant to any value ofμ within the range of [0.5, 5]so we routinely used μ=2in our simulations and experiments. We also applied consistently a stopping criterion of ||pn+1pn||<0.002. The results of applying TV denoising on the partial derivatives maps in Fig. 2(a) and 2(b) are shown in Fig. 2(c) and 2(d), respectively. It can be visually judged that the noise caused by phase residues has been totally removed at the cost of slight distortion in the phase derivative information; we will show how to remove this distortion finally in sections 3.5-3.6.

3.4 Integration of phase derivative

After obtaining the denoised phase derivative mapsϕ^x,ϕ^y, a simple integration gives a good approximation of the unwrapped phase map. However, the unwrapped phase values on the top and left lines of the image should be given as starting values. They can be obtained by setting the top-left corner unwrapped phase as the same as the wrapped phase, and integrating from left to right using ϕ^xto obtain the unwrapped phase values for the top line while integrating from top to bottom using ϕ^y to obtain the unwrapped phase values for the left line. The integration over the whole image can then be conducted based on the unwrapped boundary information. Figure 2(e) shows the unwrapped phase map ϕ^ thus obtained. To visually evaluate the fidelity of the approximated unwrapped phase map to the original wrapped phase map, we rewrap the result in Fig. 2(e) into Fig. 2(f). By comparing Fig. 1(d) and Fig. 2(f), it can be observed that the rewrapped phase map is similar to the original wrapped phase, but not identical, which means that TV denoising has slightly modify the true phase gradient, and further method for the error correction may be required.

3.5 Error metric

To quantify how well this 1st approximation (1st mode) obtained so far represents the correct unwrapped phase map, we define a residual wrapped phase map as

ϕr=W{ϕϕ^}=W{ϕwϕ^}.
(10)

If the difference between ϕ^ and ϕis within 2πfor every pixel, then the residual wrapped phase map ϕrwill only contain phase jumps from the phase residues. On the other hand, if the approximation ϕ^is underestimated or overestimated for more than2πfor some areas, apparent 2πphase jumps will present and form loops in the residual wrapped phaseϕr. Thus a simple algorithm can be designed to count the total number of 2πphase jumps and denote it as N. If Nis less than a preset criterion Nc(set to 0.1% of the total pixel number by default), we can conclude that the approximated unwrapped phase is within 2πfrom the correct unwrapped phase, and no further processing on the residual wrapped phase ϕris required. On the other hand, if N exceeds Nc, the residual wrapped phase still contains loops of 2πjumps and these 2πjumps should be further eliminated by adding more modes to the approximated unwrapped phase ϕ^, which will be described in the following section.

Figure 2(g) shows the residual wrapped phase of 2(e) where 2πjumps due to insufficient unwrapping can be clearly observed and quantified. This residual wrapped phase map will be further processed to extract its approximated unwrapped phase.

3.6 Final unwrapping result from iteration

When apparent 2πphase jumps are present in the residual wrapped phase map, the process described in sections 3.2-3.4 needs to be applied to this residual wrapped phase map again to extract an approximated unwrapped phase map for it. In such case, the 1st mode is assigned as ϕ^1and the 2nd, 3rd modes are calculated successively from the sequence of residual wrapped phase and assigned as ϕ^2,ϕ^3until no apparent 2π jumps remain in the final residual wrapped phase ϕrf. The final unwrapped phase can then be obtained by adding up all the modes and the final residual wrapped phase as

ϕ=k=1,2ϕ^k+ϕrf.
(11)

For a sufficiently smooth original unwrapped phase map, normally 1-2 iterations suffice to converge to a jump-free reconstruction. However, to avoid long lasting cycles when the wrapped phase is extremely noisy and N<Ncis never satisfied, a forced stopping criterion is implemented which stops the cycling when the total iteration number exceeds a certain number like 500. In our shearography example shown in Fig. 2, two iterations were required. Figure 2(h) shows the final residual wrapped phase ϕrfobtained by subtracting 2(g) to its approximated unwrapped phase. It can be observed and quantified that only phase jumps due to residues but no apparent2π phase jumps are present. Thus the final unwrapped phase map is obtained by Eq. (11) and shown in Fig. 2(i) where a satisfying result similar to Flynn’s minimum discontinuity method (Fig. 1(i)) has been obtained but the need to define any implicit unwrapping path has been eliminated.

3.7 Flowchart of the algorithm

A flowchart describing the whole algorithm is given in Fig. 3
Fig. 3 Flowchart of the TV unwrapping algorithm.
for convenience of readers.

4. Experiment evaluation

It should be noted that the examples given in Fig. 1 and 2 are real experiment data from shearography. To further evaluate the proposed TV unwrapping method, we use experiment data from electronic speckle pattern interferometry (ESPI) as well.

Figure 4
Fig. 4 Evaluation of the TV unwrapping algorithm using experiment data from ESPI.
shows the ESPI phase which depicts the surface deflection of a test object [30

30. Y. H. Huang, S. Y. Hung, F. Janabi-Sharifi, W. Wang, and Y. S. Liu, “Quantitative phase retrieval in dynamic laser speckle interferometry,” Opt. Lasers Eng. 50(4), 534–539 (2012). [CrossRef]

]. As the experiments are carried out using two-step phase shifting method [30

30. Y. H. Huang, S. Y. Hung, F. Janabi-Sharifi, W. Wang, and Y. S. Liu, “Quantitative phase retrieval in dynamic laser speckle interferometry,” Opt. Lasers Eng. 50(4), 534–539 (2012). [CrossRef]

], the initial wrapped phase (top left image) contains more noise than normal four step phase shifting method or temporal phase analysis method. The wrapped phase filter described in section 2 has been applied to the initial wrapped phase for 5 times (top middle image) and 20 times (top right image) to remove the noise. Note that in the top right image, phase residues can be identified at area with dense fringes and these residues will cause problem in unwrapping. However, thanks to the denoising capability of TV unwrapping algorithm, the unwrapped results for middle and right columns have been directly obtained and shown in Fig. 4. The unwrapping of the initial noisy wrapped phase shown in top left column is also trivial. One can always filter the noisy wrapped phase first and unwrap the denoised wrapped phase map (like middle and right column in Fig. 4), then use the resultant unwrapped phase as a reference to obtain the unwrapped phase map as shown in left column of Fig. 4.

5. Simulation evaluation using wrapped phase with height discontinuity

The proposed TV unwrapping algorithm has also been evaluated using simulated wrapped phase map contains height discontinuity according to the parameters described in section 4 in reference [24

24. J. F. Weng and Y. L. Lo, “Integration of robust filters and phase unwrapping algorithms for image reconstruction of objects containing height discontinuities,” Opt. Express 20(10), 10896–10920 (2012). [CrossRef] [PubMed]

]. We generate two noisy wrapped phase maps. The first one, as shown in the first row of Fig. 5
Fig. 5 Evaluation of the TV unwrapping algorithm using simulated data with height discontinuity.
, contains speckle noise only with a parameter 0.18 (i.e. imnoise (each specklegram, ‘speckle’, 0.18) in Matlab code). The second wrapped phase, as shown in the second row of Fig. 5, contains speckle noise of 0.08, residual noise of 0.35 and edge noise of 0.01 (i.e. imnoise (each specklegram, ‘speckle’, 0.08), imnoise (each specklegram, ‘salt & pepper’, 0.35) and imnoise(edges at each specklegram, ‘salt & pepper’, 0.01) in Matlab code), which is the same as the simulation in reference [24

24. J. F. Weng and Y. L. Lo, “Integration of robust filters and phase unwrapping algorithms for image reconstruction of objects containing height discontinuities,” Opt. Express 20(10), 10896–10920 (2012). [CrossRef] [PubMed]

]. In the evaluation, the proposed TV unwrapping algorithm can directly obtained a correct unwrapped phase for the case of speckle noise only shown in first row. To unwrap the second row, the initial noisy wrapped phase map should be filtered for several times by the wrapped phase filters described in section 2, and then unwrapped using the TV unwrapping algorithm. After obtaining an unwrapped phase of the slightly denoised wrapped phase, it can be used as a reference and the unwrapped phase of the initial noisy wrapped phase can be readily obtained as shown in second row of Fig. 5.

6. Simulation evaluation of effect of noise and fringe density

Simulation is always more controllable and convenient, we thus conduct numerical simulations to investigate the performance of proposed unwrapping method for different amount of noise and different fringe density. For the simulation results shown in this section, we generated “original” phase maps of 201×201 pixels with the form
φ(x,y)=2πλφ0xexp((x2+y2)2a2)
(12)
where φ0is the “height” information which controls the magnitude of phase, x,y[1, 1]is the normalized distance units, a = 0.3 is the phase map spread factor, and the wavelengthλis set to 1.

We then generate four phase-shifted fringe maps according to
Ii(x,y)=1+cos[φ(x,y)+δi]+noise(x,y)
(13)
where i=1,2,3,4 are the image numbers and δi=(i1)π/2are the phase-shift amounts for each image. Standard four-step phase shifting method is then used to recover the wrapped phase map from the four fringe maps with different statistics for noise(x,y). Before applying the TV unwrapping algorithm, the wrapped phase filtering process is applied three times to obtain a denoised wrapped phase.

In the first simulation, we add to the fringe maps zero-mean Gaussian noise~N(0,σ2) with varying standard deviation σ. We then unwrap the noisy wrapped phase map using the proposed TV unwrapping method. The results are shown in Fig. 6
Fig. 6 Typical fringe maps (top row) with Gaussian noise, wrapped phase map (middle row) and unwrapped phase map (bottom row) obtained by the proposed TV unwrapping method. The fringe density controlling parameter φ0 in Eq. (12) with is fixed at 20 for these simulations.
where the signal to noise ratio (SNR) is defined as the ratio of the average intensity in Eq. (13) to the standard deviation of noise σ. It can be seen that the algorithm is able to deal with very noisy wrapped phase map. However, the unwrapping for the case of SNR = 0.3 is incomplete due to the severe noise and the compulsory stopping criteria.

The computation needed for unwrapping various noisy wrapped phases is also plotted in Fig. 7
Fig. 7 Iteration times required to achieve a satisfying unwrapped phase Vs. Gaussian noise level using simulated fringe maps according to Eq. (12) with φ0 fixed at 20. The mean and standard deviation values are calculated from 100 random realizations.
where the simulation was repeated 100 times with different random realizations; the mean and standard deviation of the number of iterations are reported. As the ground true phase map is available in this simulation, we also calculate the normalized mean squared error for the reconstructed phase as shown in Fig. 8
Fig. 8 Normalized mean squared error of the reconstructed phase Vs. Gaussian noise level using simulated fringe maps according to Eq. (12) with φ0 fixed at 20. The mean and standard deviation values are calculated from 100 random realizations.
. It can be observed from Figs. 7 and 8 that for SNR larger than 1, very little computation is needed and the reconstructed phase is almost error free. For SNR ranging from 1 to 0.5, the computation needed is increased dramatically, however, the reconstructed phase error remains very small (<0.005) which means that the unwrapping is still correct. For SNR smaller than 0.5, both the computation needed and the phase error increase, indicating that the unwrapping algorithm can no longer render correct results.

It is straightforward to reason that if the magnitude of unwrapped phase is increased, then the fringes will become denser and detrimental effect of phase residues on unwrapping will be more apparent. Thus we also simulate noisy fringe maps with different fringe density to investigate this problem. In Figs. 9
Fig. 9 Typical fringe maps (top row) with Gaussian noise with SNR of 0.8, wrapped phase map (middle row) and unwrapped phase map (bottom row) obtained by the proposed TV unwrapping method. The fringe density controlling parameter φ0 in Eq. (12) is changed from 10 to 40 at step 10.
-11
Fig. 11 Reconstruction error of unwrapped phase Vs. phase magnitude using simulated fringe maps according to Eq. (12) with fixed Gaussian noise andφ0 ranging from 10 to 40. The mean and standard deviation values are calculated from 100 random realizations.
, uniform Gaussian noise with SNR of 0.8 is added to the fringe maps with varying phase magnitude. Figure 9 shows the typical unwrapping results, while Figs. 10
Fig. 10 Iteration times required to achieve a satisfying unwrapped phase Vs. phase magnitude using simulated fringe maps according to Eq. (12) with fixed Gaussian noise andφ0 ranging from 10 to 40. The mean and standard deviation values are calculated from 100 random realizations.
and 11 show the computations and reconstructed phase error respectively. Unlike the increase of noise which will affect the whole wrapped phase map, the increase of phase magnitude only affects the central region where the fringe density is the highest. Thus the unwrapping algorithm may be subject to underestimation or other kind of failure at this region. It is also noticed from Figs. 10 and 11 that the computation needed and the reconstructed phase error both increase dramatically when the phase magnitude φ0approach 25, this is partially due to the fact that the fringe density in the central area is approaching 4 pixels per fringe, which is merely two times of the sampling limit. At such high fringe density, the wrapped phase filtering process has begun to produce abundant phase residues at central area and cause both computation and reconstruction error to increase.

It is also observed from the simulations that the first mode is sufficient for extracting the correct unwrapped phase with most noise level and phase amplitude. Second mode is only necessary at extremely noisy cases. In terms of computation time, the whole unwrapping process for the 576×768 pixels image shown in Fig. 1(d) (maximum phase gradient is 0.5 rad/pixel) needs 38 iterations and costs 2.04 seconds in an i7 2.8G laptop, and the simulated 201×201 pixelsimage with Gaussian noise and φ0 = 20 as shown in 2nd column of Fig. 9 (maximum phase gradient is about 1.25 rad/pixel) costs only 0.12 seconds. The algorithm can still be optimized for more efficient computation, especially for the TV iteration part.

7. Conclusion and future works

In conclusion, we have demonstrated a novel use for TV minimization to simultaneously remove phase noise and preserve phase gradient. This feature of TV minimization has made it outperform other kind of denoising methods and render TV unwrapping an effective method, especially for continuous phase map without abrupt change. The first mode unwrapped phase obtained from TV has provided a good approximation to the correct unwrapped phase and made the residual wrapped phase contains much sparser 2πloops. Based on the residual wrapped phase, the unwrapping task is much easier for other kind of unwrapping methods. Thus the TV unwrapping method can also been combined with existing unwrapping algorithms to tackle more complex problem. Future work will include extending the proposed TV minimization scheme for unwrapping other complex wrapped phase such as that from InSAR, and investigation of alternative compressive sensing methods for phase unwrapping.

Acknowledgments

The authors would like to thank Dr. L. J. Chen for invaluable discussions and use of an unwrapping program to produce results for comparison. This research was funded by the National Research Foundation (NRF) of Singapore through the Singapore-MIT Alliance for Research and Technology (SMART) Center for Environmental Sensing and Modeling (CENSAM), and the Chevron Energy Technology Company. Yi Liu also gratefully acknowledges the Xerox MIT Fellowship.

References and links

1.

Y. H. Huang, Y. S. Liu, S. Y. Hung, C. G. Li, and F. Janabi-Sharifi, “Dynamic phase evaluation in sparse-sampled temporal speckle pattern sequence,” Opt. Lett. 36(4), 526–528 (2011). [CrossRef] [PubMed]

2.

E. Volkl, L. F. Allard, and D. C. Joy, eds., Introduction to Electron holography (Plenum, 1999).

3.

B. Osmanoglu, T. H. Dixon, S. Wdowinski, and E. Cabral-Cano, “On the importance of path for phase unwrapping in synthetic aperture radar interferometry,” Appl. Opt. 50(19), 3205–3220 (2011). [CrossRef] [PubMed]

4.

J. Langley and Q. Zhao, “Unwrapping magnetic resonance phase maps with Chebyshev polynomials,” Magn. Reson. Imaging 27(9), 1293–1301 (2009). [CrossRef] [PubMed]

5.

C. J. Tay, C. Quan, T. Wu, and Y. H. Huang, “Integrated method for 3-D rigid-body displacement measurement using fringe projection,” Opt. Eng. 43(5), 1152–1159 (2004). [CrossRef]

6.

D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software, (John Wiley & Sons, 1998).

7.

R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. 23(4), 713–720 (1988). [CrossRef]

8.

Y. Xu and C. Ai, ““Simple and effective phase unwrapping technique,” Interferometry IV: Techniques and Analysis,” Proc. SPIE 2003, 254–263 (1993). [CrossRef]

9.

C. Prati, M. Giani, and N. Leuratti, “SAR interferometry: a 2-D phase unwrapping technique based on phase and absolute values information,” in Proceedings of the 1990 International Geoscience and Remote Sensing Symposium (IEEE, 1990), pp. 2043–2046.

10.

T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14(10), 2692–2702 (1997). [CrossRef]

11.

I. Shalem and I. Yavneh, “A multilevel graph algorithm for two dimensional phase unwrapping,” Comput. Vis. Sci. 11(2), 89–100 (2008). [CrossRef]

12.

T. M. Venema and J. D. Schmidt, “Optical phase unwrapping in the presence of branch points,” Opt. Express 16(10), 6985–6998 (2008). [CrossRef] [PubMed]

13.

M. Gdeisat, M. Arevalillo-Herráez, D. Burton, and F. Lilley, “Three-dimensional phase unwrapping using the Hungarian algorithm,” Opt. Lett. 34(19), 2994–2996 (2009). [CrossRef] [PubMed]

14.

S. S. Gorthi, G. Rajshekhar, and P. Rastogi, “Strain estimation in digital holographic interferometry using piecewise polynomial phase approximation based method,” Opt. Express 18(2), 560–565 (2010). [CrossRef] [PubMed]

15.

J. Bioucas-Dias, V. Katkovnik, J. Astola, and K. Egiazarian, “Absolute phase estimation: adaptive local denoising and global unwrapping,” Appl. Opt. 47(29), 5358–5369 (2008). [CrossRef] [PubMed]

16.

D. C. Ghiglia and L. A. Romero, “Direct phase estimation from phase differences using fast elliptic partial differential equation solvers,” Opt. Lett. 14(20), 1107–1109 (1989). [CrossRef] [PubMed]

17.

M. A. Schofield and Y. M. Zhu, “Fast phase unwrapping algorithm for interferometric applications,” Opt. Lett. 28(14), 1194–1196 (2003). [CrossRef] [PubMed]

18.

R. Yamaki and A. Hirose, “Singularity-spreading phase unwrapping,” IEEE Trans. Geosci. Rem. Sens. 45(10), 3240–3251 (2007). [CrossRef]

19.

S. Tomioka, S. Heshmat, N. Miyamoto, and S. Nishiyama, “Phase unwrapping for noisy phase maps using rotational compensator with virtual singular points,” Appl. Opt. 49(25), 4735–4745 (2010). [CrossRef] [PubMed]

20.

M. D. Pritt, “Congruence in least-squares phase unwrapping,” in Proceedings Vol. II: Remote Sensing - a Scientific Vision for Sustainable Development, IGARSS '97 - 1997 International Geoscience and Remote Sensing Symposium, 1997), pp.875–877.

21.

J. Strand and T. Taxt, “Performance evaluation of two-dimensional phase unwrapping algorithms,” Appl. Opt. 38(20), 4333–4344 (1999). [CrossRef] [PubMed]

22.

E. Zappa and G. Busca, “Comparison of eight unwrapping algorithms applied to Fourier-transform profilometry,” Opt. Lasers Eng. 46(2), 106–116 (2008). [CrossRef]

23.

J. Parkhurst, G. Price, P. Sharrock, and C. Moore, “Phase unwrapping algorithms for use in a true real-time optical body sensor system for use during radiotherapy,” Appl. Opt. 50(35), 6430–6439 (2011). [CrossRef] [PubMed]

24.

J. F. Weng and Y. L. Lo, “Integration of robust filters and phase unwrapping algorithms for image reconstruction of objects containing height discontinuities,” Opt. Express 20(10), 10896–10920 (2012). [CrossRef] [PubMed]

25.

M. A. Navarro, J. C. Estrada, M. Servin, J. A. Quiroga, and J. Vargas, “Fast two-dimensional simultaneous phase unwrapping and low-pass filtering,” Opt. Express 20(3), 2556–2561 (2012). [CrossRef] [PubMed]

26.

Y. H. Huang, F. Janabi-Sharifi, Y. Liu, and Y. Y. Hung, “Dynamic phase measurement in shearography by clustering method and Fourier filtering,” Opt. Express 19(2), 606–615 (2011). [CrossRef] [PubMed]

27.

Y. H. Huang, S. P. Ng, L. Liu, Y. S. Chen, and M. Y. Y. Hung, “Shearographic phase retrieval using one single specklegram: a clustering approach,” Opt. Eng. 47(5), 054301 (2008). [CrossRef]

28.

Y. H. Huang, Y. S. Liu, S. Y. Hung, C. G. Li, and F. Janabi-Sharifi, “Dynamic phase evaluation in sparse-sampled temporal speckle pattern sequence,” Opt. Lett. 36(4), 526–528 (2011). [CrossRef] [PubMed]

29.

A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis. 20(1/2), 89–97 (2004). [CrossRef]

30.

Y. H. Huang, S. Y. Hung, F. Janabi-Sharifi, W. Wang, and Y. S. Liu, “Quantitative phase retrieval in dynamic laser speckle interferometry,” Opt. Lasers Eng. 50(4), 534–539 (2012). [CrossRef]

OCIS Codes
(090.2880) Holography : Holographic interferometry
(120.2650) Instrumentation, measurement, and metrology : Fringe analysis
(120.5050) Instrumentation, measurement, and metrology : Phase measurement
(280.6730) Remote sensing and sensors : Synthetic aperture radar
(100.5088) Image processing : Phase unwrapping
(120.6165) Instrumentation, measurement, and metrology : Speckle interferometry, metrology

ToC Category:
Instrumentation, Measurement, and Metrology

History
Original Manuscript: April 23, 2012
Revised Manuscript: May 23, 2012
Manuscript Accepted: June 1, 2012
Published: June 11, 2012

Citation
Howard Y. H. Huang, L. Tian, Z. Zhang, Y. Liu, Z. Chen, and G. Barbastathis, "Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising," Opt. Express 20, 14075-14089 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-13-14075


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. Y. H. Huang, Y. S. Liu, S. Y. Hung, C. G. Li, and F. Janabi-Sharifi, “Dynamic phase evaluation in sparse-sampled temporal speckle pattern sequence,” Opt. Lett.36(4), 526–528 (2011). [CrossRef] [PubMed]
  2. E. Volkl, L. F. Allard, and D. C. Joy, eds., Introduction to Electron holography (Plenum, 1999).
  3. B. Osmanoglu, T. H. Dixon, S. Wdowinski, and E. Cabral-Cano, “On the importance of path for phase unwrapping in synthetic aperture radar interferometry,” Appl. Opt.50(19), 3205–3220 (2011). [CrossRef] [PubMed]
  4. J. Langley and Q. Zhao, “Unwrapping magnetic resonance phase maps with Chebyshev polynomials,” Magn. Reson. Imaging27(9), 1293–1301 (2009). [CrossRef] [PubMed]
  5. C. J. Tay, C. Quan, T. Wu, and Y. H. Huang, “Integrated method for 3-D rigid-body displacement measurement using fringe projection,” Opt. Eng.43(5), 1152–1159 (2004). [CrossRef]
  6. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software, (John Wiley & Sons, 1998).
  7. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci.23(4), 713–720 (1988). [CrossRef]
  8. Y. Xu and C. Ai, ““Simple and effective phase unwrapping technique,” Interferometry IV: Techniques and Analysis,” Proc. SPIE2003, 254–263 (1993). [CrossRef]
  9. C. Prati, M. Giani, and N. Leuratti, “SAR interferometry: a 2-D phase unwrapping technique based on phase and absolute values information,” in Proceedings of the 1990 International Geoscience and Remote Sensing Symposium (IEEE, 1990), pp. 2043–2046.
  10. T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A14(10), 2692–2702 (1997). [CrossRef]
  11. I. Shalem and I. Yavneh, “A multilevel graph algorithm for two dimensional phase unwrapping,” Comput. Vis. Sci.11(2), 89–100 (2008). [CrossRef]
  12. T. M. Venema and J. D. Schmidt, “Optical phase unwrapping in the presence of branch points,” Opt. Express16(10), 6985–6998 (2008). [CrossRef] [PubMed]
  13. M. Gdeisat, M. Arevalillo-Herráez, D. Burton, and F. Lilley, “Three-dimensional phase unwrapping using the Hungarian algorithm,” Opt. Lett.34(19), 2994–2996 (2009). [CrossRef] [PubMed]
  14. S. S. Gorthi, G. Rajshekhar, and P. Rastogi, “Strain estimation in digital holographic interferometry using piecewise polynomial phase approximation based method,” Opt. Express18(2), 560–565 (2010). [CrossRef] [PubMed]
  15. J. Bioucas-Dias, V. Katkovnik, J. Astola, and K. Egiazarian, “Absolute phase estimation: adaptive local denoising and global unwrapping,” Appl. Opt.47(29), 5358–5369 (2008). [CrossRef] [PubMed]
  16. D. C. Ghiglia and L. A. Romero, “Direct phase estimation from phase differences using fast elliptic partial differential equation solvers,” Opt. Lett.14(20), 1107–1109 (1989). [CrossRef] [PubMed]
  17. M. A. Schofield and Y. M. Zhu, “Fast phase unwrapping algorithm for interferometric applications,” Opt. Lett.28(14), 1194–1196 (2003). [CrossRef] [PubMed]
  18. R. Yamaki and A. Hirose, “Singularity-spreading phase unwrapping,” IEEE Trans. Geosci. Rem. Sens.45(10), 3240–3251 (2007). [CrossRef]
  19. S. Tomioka, S. Heshmat, N. Miyamoto, and S. Nishiyama, “Phase unwrapping for noisy phase maps using rotational compensator with virtual singular points,” Appl. Opt.49(25), 4735–4745 (2010). [CrossRef] [PubMed]
  20. M. D. Pritt, “Congruence in least-squares phase unwrapping,” in Proceedings Vol. II: Remote Sensing - a Scientific Vision for Sustainable Development, IGARSS '97 - 1997 International Geoscience and Remote Sensing Symposium, 1997), pp.875–877.
  21. J. Strand and T. Taxt, “Performance evaluation of two-dimensional phase unwrapping algorithms,” Appl. Opt.38(20), 4333–4344 (1999). [CrossRef] [PubMed]
  22. E. Zappa and G. Busca, “Comparison of eight unwrapping algorithms applied to Fourier-transform profilometry,” Opt. Lasers Eng.46(2), 106–116 (2008). [CrossRef]
  23. J. Parkhurst, G. Price, P. Sharrock, and C. Moore, “Phase unwrapping algorithms for use in a true real-time optical body sensor system for use during radiotherapy,” Appl. Opt.50(35), 6430–6439 (2011). [CrossRef] [PubMed]
  24. J. F. Weng and Y. L. Lo, “Integration of robust filters and phase unwrapping algorithms for image reconstruction of objects containing height discontinuities,” Opt. Express20(10), 10896–10920 (2012). [CrossRef] [PubMed]
  25. M. A. Navarro, J. C. Estrada, M. Servin, J. A. Quiroga, and J. Vargas, “Fast two-dimensional simultaneous phase unwrapping and low-pass filtering,” Opt. Express20(3), 2556–2561 (2012). [CrossRef] [PubMed]
  26. Y. H. Huang, F. Janabi-Sharifi, Y. Liu, and Y. Y. Hung, “Dynamic phase measurement in shearography by clustering method and Fourier filtering,” Opt. Express19(2), 606–615 (2011). [CrossRef] [PubMed]
  27. Y. H. Huang, S. P. Ng, L. Liu, Y. S. Chen, and M. Y. Y. Hung, “Shearographic phase retrieval using one single specklegram: a clustering approach,” Opt. Eng.47(5), 054301 (2008). [CrossRef]
  28. Y. H. Huang, Y. S. Liu, S. Y. Hung, C. G. Li, and F. Janabi-Sharifi, “Dynamic phase evaluation in sparse-sampled temporal speckle pattern sequence,” Opt. Lett.36(4), 526–528 (2011). [CrossRef] [PubMed]
  29. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis.20(1/2), 89–97 (2004). [CrossRef]
  30. Y. H. Huang, S. Y. Hung, F. Janabi-Sharifi, W. Wang, and Y. S. Liu, “Quantitative phase retrieval in dynamic laser speckle interferometry,” Opt. Lasers Eng.50(4), 534–539 (2012). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited