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

Howard Y. H. Huang; L. Tian; Z. Zhang; Y. Liu; Z. Chen; G. Barbastathis

doi:10.1364/OE.20.014075

1. Introduction

Phase unwrapping is an essential step for quantitative measurement in a wide variety of interferometry based techniques such as speckle interferometry [1], electron holography [2], interferometric synthetic aperture radar (InSAR) [3], magnetic resonance imaging (MRI) [4], fringe projection [5], 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].

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–13] 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,15]. 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 _{_$L^{p}$}norm. This turns the problem of phase unwrapping into solving a discretized partial differential equation [6,16,17], in particular a Poisson equation when _{_$L^{2}$} 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], thus post-processing techniques [18,19] including the congruence process [20] have been applied to improve the algorithms. As a result, very good unwrapping results [21–23] have been obtained.

As the difficulty of unwrapping is caused by noise and phase residues in the wrapped phase map, another category of unwrapping algorithms tends to propose effective filters for noise filtering, and then use either path-following or path-independent algorithm to unwrap the denoised wrapped phase map [24,25]. In this paper, we propose a path-independent phase unwrapping algorithm employing the same idea of wrapped phase denoising. Our algorithm is inspired by compressive sensing, in particular the total-variation (TV) minimization principle. It should be emphasized that our proposed method doesn’t need to know the noise statistics. In the implementation, we use smoothness and differentiability of the unwrapped phase map as prior. These assumptions are always true for phase maps obtained from speckle interferograms representing surface deformation or strain. Unlike previous work [26] for direct wrapped phase map denoising, the proposed method first determines the phase gradient from the noisy wrapped phase map, TV minimization then exploits the fact that the differentiation process amplifies the noise content and smoothens the _{_{$2 π$}}jumps in the wrapped phase map. Thus, the noise can be successfully separated from the genuine phase gradient map; that is, TV effectively acts as a denoising process in the phase gradient domain.

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 $\dots$ 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 of_{_{$2 π$}}as:

ϕ_{w} = ϕ - 2 n π n \in Z, ϕ_{w} \in (- π, π] .

The task of phase unwrapping is thus to add a correct _{_{$2 n π$}} 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] 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]. 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] 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] 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.

Figure 1 shows a typical experimental wrapped phase map from speckle interferometry for dynamic deformation measurement [27,28]. In such applications for dynamic measurement, only one speckle pattern can be captured at each time interval since the object is continuously moving, thus the available information is less than that for static measurement where more than three phase-shifted speckle patterns can be captured at each interval for noise suppression. As a result, larger amount of noise will present in phase maps in dynamic measurement [27,28] and the corresponding phase unwrapping task will be more difficult. Figure 1(a) shows an experimental interference fringe representing the surface deflection. Figure 1(b) shows the initial noisy wrapped phase map obtained using a clustering method [27]. Normally such noisy wrapped phase map needs to go through a noise-filtering process to make it easier and time efficient for unwrapping. A common filtering method is to convert the wrapped phase map into fringes by cosine and sine functions, filter the fringes using low-pass filters (like convolution with 3x3 equal kernel) or median filters, and then reconstruct a wrapped phase map using the denoised fringes maps (this process is referred to as wrapped phase filtering in this paper). Figure 1(c) and 1(d) show the denoised wrapped phase map after applying the wrapped phase filter for 5 times and 50 times respectively. It can be observed that the wrapped phase filter can remove most of the noise; however, it fails to remove the phase residues within the “HN” area in Fig. 1(d) since the fringe maps of residues is smooth and both low-pass and median filters are ineffective. These residues have posed problems for some unwrapping algorithms. Figures 1(e) - 1(i) show the unwrapping results from a few popular path-following algorithms. It can be observed that the raster scanning method (Fig. 1(e)) severely propagates the errors to the right of the HN area, ruining the entire unwrapped phase map. Goldstein’s branch cut method [7], though it works well for most interferometric synthetic aperture radar (InSAR) wrapped phase maps, seems not very suitable for unwrapping phase maps where noise is clustered together in certain areas, hindering the generation of correct branch cuts. As a result, errors in the HN area dominate the unwrapped phase map (Fig. 1(f)). The quality guided [8] and mask cut [9] methods using quality maps from the maximum phase gradient have successfully unwrapped the noise-free areas first and confine the errors inside the HN area. However, it can be judged visually that the information of this area is incorrect (Fig. 1(g) and 1(h)). Flynn’s minimum discontinuity algorithm [10] is the only algorithm that can give a reasonable unwrapped phase map thanks to its capability to generate a result with minimum surface discontinuity (Fig. 1(i)).

Fig. 1 An interferometric fringe from shearography (a), the initial noisy wrapped phase map obtained (b), denoised wrapped phase by applying the wrapped phase filtering method 5 times (c) and 50 times (d), and the unwrapping results from raster scanning (e), Goldstein’s branch cut algorithm [result is displayed as logarithm because of the dominance of noise amplified by the algorithm at point indicated as “HN” in the picture] (f), Quality guided algorithm with maximum gradient (g), Mask cut algorithm with maximum gradient (h) and Flynn’s minimum discontinuity method (i).

Download Full Size | PDF

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 {ψ} = {\begin{cases} 2 \arctan [\frac{\sin (ψ)}{1 + \cos (ψ)}] \cos (ψ) \neq - 1 \\ π \cos (ψ) = - 1 \end{cases}

where _{_{$\arctan [\cdot]$}} is a normal arctangent function with range _{_{$(- π / 2, π / 2)$}}and the wrapping function _{_{$W {\cdot}$}} 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} = \frac{\partial ϕ}{\partial x} = W {\frac{\partial ϕ_{w}}{\partial x}}

ϕ_{y} = \frac{\partial ϕ}{\partial y} = W {\frac{\partial ϕ_{w}}{\partial y}} .

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 {\cdot}$}}, 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.

Figures 2(a) and 2(b) show the phase derivatives thus obtained from the wrapped phase map shown in Fig. 1(d). (A standard median filter is used to improve the image contrast). It can be clearly observed that the_{_{$2 π$}} phase jumps in the wrapped phase map have been removed and phase derivatives are smooth for most areas except the central part where noise is amplified. The next section describes how TV removes this amplified noise.

Fig. 2 The partial derivatives of Fig. 1(d) in x- (a) and y- directions (b), the TV denoised partial derivatives with regard to x- (c) and y- directions (d), the unwrapped phase map by integration the denoised phase derivatives (e) and its rewrapped phase map (f), the first residual wrapped phase (g), the second residual wrapped phase (h) and the final unwrapped phase map from the proposed algorithm (i).

Download Full Size | PDF

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] as we apply it here. Using TV, we apply a smoothness prior on the phase derivative _{_{$ϕ_{x}$}} and subsequently calculate a denoised estimate _{_{${\hat{ϕ}}_{x}$}} according to

{\hat{ϕ}}_{x} \underset{{\hat{ϕ}}_{x}}{= arg \min} {| | {\hat{ϕ}}_{x} - ϕ_{x} | |^{2} / (2 μ) + T V ({\hat{ϕ}}_{x})}

where _{_{$ϕ_{x}$}}is the input noisy phase derivative, _{_{$| | \cdot | |$}}is the Euclidean norm, _{_$μ$}is a regularization parameter, and _{_{$T V (\cdot)$}} is a functional that determines the total variation of its argument. If we define the TV functional as_{_{$T V (ψ) = \sum | | \nabla ψ | |$}}, i.e. the sum of the gradient magnitude over all pixels, it has been proved [29] that the solution to the minimization problem (5) is given by

{\hat{ϕ}}_{x} = ϕ_{x} - π_{μ κ} (ϕ_{x})

where _{_{$π_{μ κ} (ϕ_{x})$}}is the nonlinear projection of _{_{$ϕ_{x}$}}. Computing the nonlinear projection _{_{$π_{μ κ} (ϕ_{x})$}}is identical to solving the following problem

\underset{\vec{p}}{arg \min} {| | μ \nabla \cdot \vec{p} - ϕ_{x} | |^{2} : | \vec{p} |^{2} \leq 1}

and _{_{$π_{μ κ} (ϕ_{x})$}}can then be determined as

π_{μ κ} (ϕ_{x}) = μ \nabla \cdot \vec{p}

In [29], Chambolle proposed a semi-implicit gradient descent algorithm to determine _{_$\vec{p}$}as the convergent value of a vector sequence _{_{${\vec{p}}^{n}$}} of

{\vec{p}}^{n + 1} = \frac{{\vec{p}}^{n} + τ [\nabla (\nabla \cdot {\vec{p}}^{n} - ϕ_{x} / μ)]}{1 + τ | | \nabla (\nabla \cdot {\vec{p}}^{n} - ϕ_{x} / μ) | |}

where the iteration is initiated at _{_{${\vec{p}}^{1} = (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 _{_{$μ = 2$}}in our simulations and experiments. We also applied consistently a stopping criterion of _{_{$| | \nabla \cdot {\vec{p}}^{n + 1} - \nabla \cdot {\vec{p}}^{n} | |_{\infty} < 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_{_{${\hat{ϕ}}_{x}, {\hat{ϕ}}_{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 _{_{${\hat{ϕ}}_{x}$}}to obtain the unwrapped phase values for the top line while integrating from top to bottom using _{_{${\hat{ϕ}}_{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 _{_$\hat{ϕ}$} 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 {ϕ - \hat{ϕ}} = W {ϕ_{w} - \hat{ϕ}} .

If the difference between _{_$\hat{ϕ}$} and _{_$ϕ$}is within _{_{$2 π$}}for every pixel, then the residual wrapped phase map _{_{$ϕ_{r}$}}will only contain phase jumps from the phase residues. On the other hand, if the approximation _{_$\hat{ϕ}$}is underestimated or overestimated for more than_{_{$2 π$}}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 _{_$N$}is less than a preset criterion _{_{$N_{c}$}}(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 _{_{$ϕ_{r}$}}is required. On the other hand, if _{_$N$} exceeds _{_{$N_{c}$}}, 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 _{_$\hat{ϕ}$}, 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 _{_{${\hat{ϕ}}_{1}$}}and the 2nd, 3rd $\dots$ modes are calculated successively from the sequence of residual wrapped phase and assigned as _{_{${\hat{ϕ}}_{2}, {\hat{ϕ}}_{3} \dots$}}until no apparent _{_{$2 π$}} jumps remain in the final residual wrapped phase _{_{$ϕ_{r f}$}}. The final unwrapped phase can then be obtained by adding up all the modes and the final residual wrapped phase as

ϕ = \sum_{k = 1, 2 \dots} {\hat{ϕ}}_{k} + ϕ_{r f}

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$}<_{_{$N_{c}$}}is 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 _{_{$ϕ_{r f}$}}obtained by subtracting 2(g) to its approximated unwrapped phase. It can be observed and quantified that only phase jumps due to residues but no apparent_{_{$2 π$}} 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 for convenience of readers.

Fig. 3 Flowchart of the TV unwrapping algorithm.

Download Full Size | PDF

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 shows the ESPI phase which depicts the surface deflection of a test object [30]. As the experiments are carried out using two-step phase shifting method [30], 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.

Fig. 4 Evaluation of the TV unwrapping algorithm using experiment data from ESPI.

Download Full Size | PDF

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]. We generate two noisy wrapped phase maps. The first one, as shown in the first row of Fig. 5 , 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]. 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.

Fig. 5 Evaluation of the TV unwrapping algorithm using simulated data with height discontinuity.

Download Full Size | PDF

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 \times 201$}} pixels with the form

φ (x, y) = \frac{2 π}{λ} φ_{0} x e x p (- \frac{(x^{2} + y^{2})}{2 a^{2}})

where _{_{$φ_{0}$}}is the “height” information which controls the magnitude of phase, _{_{$x, y \in [- 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

I_{i} (x, y) = 1 + c o s [φ (x, y) + δ_{i}] + n o i s e (x, y)

where _{_{$i = 1, 2, 3, 4$}} are the image numbers and _{_{$δ_{i} = (i - 1) π / 2$}}are 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 _{_{$n o i s e (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 _{_{$n o i s e ~ 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 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.

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.

Download Full Size | PDF

The computation needed for unwrapping various noisy wrapped phases is also plotted in Fig. 7 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 . 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.

Fig. 7 Iteration times required to achieve a satisfying unwrapped phase V_s. 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.

Download Full Size | PDF

Fig. 8 Normalized mean squared error of the reconstructed phase V_s. 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.

Download Full Size | PDF

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 -11 , 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 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 _{_{$φ_{0}$}}approach 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.

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.

Download Full Size | PDF

Fig. 11 Reconstruction error of unwrapped phase V_s. 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.

Download Full Size | PDF

Fig. 10 Iteration times required to achieve a satisfying unwrapped phase V_s. 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.

Download Full Size | PDF

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 \times 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 \times 201 pixels$}}image 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]

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

Abstract

1. Introduction

2. Problem description and path-following solutions

3. Phase unwrapping using TV prior

3.1 Definition of the wrapping function

3.2 Determination of the phase gradient

3.3 TV minimization for denoising

3.4 Integration of phase derivative

3.5 Error metric

3.6 Final unwrapping result from iteration

3.7 Flowchart of the algorithm

4. Experiment evaluation

5. Simulation evaluation using wrapped phase with height discontinuity

6. Simulation evaluation of effect of noise and fringe density

7. Conclusion and future works

Acknowledgments

References and links

Cited By

Figures (11)

Equations (13)

Optics Express