Fast two-dimensional simultaneous phase unwrapping and low-pass filtering

Here, we present a fast algorithm for two-dimensional (2D) phase unwrapping which behaves as a recursive linear filter. This linear behavior allows us to easily find its frequency response and stability conditions. Previously, we published a robust to noise recursive 2D phase unwrapping system with smoothing capabilities. But our previous approach was rather heuristic in the sense that not general 2D theory was given. Here an improved and better understood version of our previous 2D recursive phase unwrapper is presented. In addition, a full characterization of it is shown in terms of its frequency response and stability. The objective here is to extend our previous unwrapping algorithm and give a more solid theoretical foundation to it. © 2012 Optical Society of America OCIS codes: (100.5088) Phase unwrapping; (100.2650) Fringe Analysis. References and links 1. K. Itoh, “Analysis of the phase unwrapping algorithm.” Appl. Opt. 21, 2470 (1982). 2. D. C. Ghihlia and M. D. Pritt, Two-dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley, 1998). 3. T. Judge, “A review of phase unwrapping techniques in fringe analysis,” Opt. Lasers Eng. 21, 199–239 (1994). 4. D. C. Ghihlia and L. A. Romero, “Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods,” J. Opt. Soc. Am. A 11 107–117 (1994). 5. J. L. Marroquin and M. Rivera, “Quadratic regularization functionals for phase unwrapping,” J. Opt. Soc. Am. A 12, 2393–2400 (1995). 6. M. Servin, F. J. Cuevas, D. Malacara, J. L. Marroquin, and R. Rodriguez-Vera, “Phase unwrapping through demodulation by use of the regularized phase-tracking technique,” Appl. Opt. 38, 1934–1941 (1999). 7. R. Goldstein, H. Zebker, and C. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988). 8. T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14, 2692–2701 (1997). 9. J. C. Estrada, M. Servin, and J. A. Quiroga, “Noise robust linear dynamic system for phase unwrapping and smoothing,” Opt. Express 19, 5126–5133 (2011). 10. B. Jähne, Digital Image Processing (Springer, 2005). 11. J. G. Proakis and D. G. Manolakis, Digital Signal Processing. Principles, Algorothims, and Applications, 3rd ed. (Prentice-Hall, October 5, 1995). 12. W.-S. Lu and A. Antoniou, Two-Dimensional Digital Filters (Marcel Dekker, Inc., 1992). 13. R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction theory,” Trans. ASME, Ser. D 83, 95–107 (1961). #157939 $15.00 USD Received 9 Nov 2011; revised 2 Dec 2011; accepted 4 Jan 2012; published 20 Jan 2012 (C) 2012 OSA 30 January 2012 / Vol. 20, No. 3 / OPTICS EXPRESS 2556


Introduction
Phase unwrapping is a process that removes the modulus 2π discontinuities of wrapped phase maps.Itoh provided a simple analysis of the one-dimensional phase unwrapping algorithm and showed that the true phase can be recovered by line integration of its wrapping differences [1].Most sequential unwrapping algorithms search the 2π discontinuities of the wrapped phase map in the spatial domain [2].The simplest unwrapping system integrates the wrapped phase differences on a line by line basis.However, in noisy conditions, erroneous phase unwrapped values may propagate outside from high noise regions corrupting the rest of the unwrapped phase [3].On the other hand, non-sequential unwrapping methods, such as the least squares method and its regularized version, behaves better on regions with high noise [4,5].However, the main drawback of least squares methods is that the unwrapped phase is obtained with a reduced dynamic range [6].Branch cut methods are the most used phase unwrapping systems [4,7,8].In a nutshell, branch cut methods is a two step process: firstly the wrapped phase inconsistencies are marked, and secondly it unwraps the phase avoiding those marked inconsistencies.Nevertheless, the problem of instability under high noise remains.That is because high noise make the marking of these phase inconsistencies very difficult.
Previously, we presented a robust to noise two-dimensional (2D) recursive phase unwrapping system and shown that it may also smooths the unwrapped phase [9].However, an important point that was not shown there, was its 2D frequency response, and stability.This omission was because we did not realize at that time that our system [9] may be perfectly studied and analyzed using standard 2D linear system theory [12].Here, we continue that work showing a improved version of our previous 2D recursive unwrapping system.Our presentation is as follows: in the next sections we show how a modified 1D first order recursive linear filter can be transformed into a 1D phase unwrapper.Afterwards, this 1D design is generalized to a 2D recursive unwrapping system considering it as a linear 2D recursive filter.We follow presenting its stability analysis and frequency response in 2D.Finally, we show some unwrapping results using experimental data and we end by drawing some conclusions.

One-dimensional recursive filters for phase unwrapping
From our previous [9], let us consider the simplest 1D general form of a first order recursive linear filter [10,11]: where φ (n) is the filters output at the current n-site, φ (n−1) is the previous output (or estimate), and τ is a parameter that controls the cut-off frequency of the filter [10][11][12].The frequency and impulse response of this linear low-pass filter are: This linear system is stable if its impulse response is absolutely summable.Thus, the system in Eq. ( 1) is stable only for τ ∈ (0, 2) [10,11].To transform Eq. ( 1) as a phase unwrapping system, let us rewrite it as: Apart from τ, and the previous estimation φ (n − 1) in the phase difference, Eq. (3) looks very similar to the following basic line unwrapping algorithm [2]: where φ (n) is the input (wrapped phase) and φ (n) the output (unwrapped phase).This important observation was firstly given in Ref. [9].The operator W [•] is the wrapping modulus 2π operator.As the input is wrapped, the phase difference has a principal value plus a residue multiple of 2π [1].The wrapping operator removes the residue and the unwrapped phase is obtained by integrating the principal values [1].Similarly, when dealing with a wrapped phase φ (n) in Eq. ( 3), it has a principal value plus a residue multiple of 2π.Then, using the wrapping operator in Eq. ( 3) we can obtain the unwrapping phase as the low-pass filtering of its principal values.Thus, the recursive first order unwrapper system is: In brief, the transition from a recursive linear low-pass filter and a low-pass filtered phase unwrapper is based on the following property of the unwrapping operator: Thus, for synthesis and analysis purposes, the wrapping operator is not taken into account and the unwrapping system is treated as a linear Infinite Impulse Response (IIR) recursive filter.

Two-dimensional recursive filter construction for phase unwrapping
Previously [9], we used sets to describe our recursive phase unwrapping system.For convenience, here we change the notation for an easier to read.Thus, our recursive system in Eq. ( 5) may be seeing as a sum of two terms, where the first term will be a predictor and the second a corrector.This Predictor and corrector concepts are actually very used in estimation theory, for example, in Kalman processing systems [13].The phase predictor is an estimation based on the previously unwrapped values, while the corrector adds input data to correct the current estimation.Thus, our recursive system in 2D has the following general form: where φ (x, y) is the 2D unwrapped pixel at (x, y), φp (x, y) is the 2D phase predictor, and φc (x, y) is the 2D corrector.Parameter τ controls the bandwidth of the system.

The predictor
A recursive filter takes information only from its previously estimated values as predictor.Therefore, for our predictor φp (x, y), we propose the mean of the 3 × 3 neighborhood of unwrapped-only pixels around pixel (x, y), marked by the indicators s(m, n) as follows: This is a convolution of φ (x, y) with an adapting kernel S whose elements are s(m, n).The elements s(m, n) indicates with ones the unwrapped pixels and with zeros the wrapped ones.Finally, S 1 is the L 1 norm of S; i.e. the sum of its elements.The kernel S is adapted for each visited pixel (x, y) being unwrapped, as illustrated in Fig.

The corrector
The corrector must include previously unwrapped phase and new wrapped input to correct our estimated prediction.Besides, the corrector must remove the residues modulus 2π from phase differences as the 1D system in Eq. ( 5) does.Taking all these criteria into account, our 2D corrector is defined as: Where kernel S is the complement of S whose elements are s(m, n), in such a way that S 1 + S 1 = 9.Then, s(m, n) is the complement of s(m, n).

System stability
The recursive filter shown in Eq. ( 7) is applied by setting an stable value for the parameter τ, and visiting each pixel following a predefined scanning strategy.For any bounded input, an stable recursive filter must obtain a bounded output.For this analysis we consider the system's corrector in Eq. ( 10) without the wrapping operator, as shown in Eq. (6).In 2D, one uses the stability Shanks stability theorem which says that a 2D recursive filter is stable if the denominator of its z-transform is non zero for (z −1 x , z −1 y ) ∈ U 2 , where U 2 is the unit bidisc defined as Taking the z-transfer function of the linear system in Eq. ( 7) (which includes Eq. ( 8) and Eq. ( 10)) we obtain the following From here, we can demonstrate that the denominator into the unit bidisc Then, the denominator of Eq. ( 12) is not zero only if |1 − τ S 1 | < 1, giving the following stability condition for τ: To guarantee that our 2D unwrapping system is stable in all neighborhood cases found in a sequential scanning, we choose our worst configuration case which is S 1 = 8 (see Fig. 1(a)).Then, in this worst case scenario the stability range is 1 4 > τ > 0.

Frequency response
In Eq. ( 12) we have shown the 2D z-transform of our recursive phase unwrapping system.Using the z-transform, the frequency response is obtained by substituting z x = e iu and z y = e iv , being u and v the spatial frequencies along the x and y axis.Figure 1 shows three possible neighborhood cases found in a sequentially scanning strategy.As example, we are going to obtain the frequency response of the most probable case, which is the one shown in panel 1(b).For this case, its 2D frequency response is: τ(e iu + e i(u+v) + e iv + e i(u−v) + 1) For illustration purposes, Figs.1(d), 1(e) and 1(f), shows the 2D graphic of the power spectrum of the frequency response corresponding to each presented case.To obtain those graphics, first obtain its corresponding frequency response from its z-transform Eq. ( 12), using τ = 0.13.

Test and results
We tested our recursive phase unwrapping system (Eq.( 7)) with a wrapped phase obtained from the experimental interferogram of 512×512 pixels shown in Fig. 2(a).To show the robust phase unwrapping capabilities of our system (Eq.( 7)), we use a simple row by row scanning strategy.In contrast, branch-cut methods need to use complex 2D scanning to avoid noise generated phase inconsistencies [7,8].We compare our results with the one of Flynn's phase unwrapping method, which is a quite popular sequential branch cut unwrapping method and it is more stable than the Goldstein's method (see Ref. [2,[7][8][9]).The result of this test is shown in Fig. 2. In Fig. 2(b), we show the unwrapped phase using our recursive phase unwrapping system, and Fig. 2(c) shows the unwrapped phase using the Flynn's algorithm.We show the unwrapped phase surface in three dimensions and the re-wrapped phase projected under this surface for comparison purposes with the original data.Here we pass our phase unwrapping system using τ = 0.013.As we can see in Fig. 2(c), our proposed 2D recursive phase unwrapping system obtains a cleaner unwrapped phase, whereas the Flynn's unwrapped phase is fare more noised and damaged.The computation load to unwrap a single pixel using our recursive system requires S 1 + 9 arithmetic operations (without taking into account the operations needed by the wrapping modulus 2π operator).Therefore, for S 1 = 1 (Fig. 1(c)) it requires 10 operations and for S 1 = 8 (Fig. 1(a)) it requires 17 operations.However our most probable case (Fig. 1(b)) requires 14 operations.Thus, taking the most probable case as the average, to unwrap a phase map of n = M × N pixels one needs an average of 14n operations.In this way, compared with the Flynn's method (and the least-squares methods), our recursive phase unwrapping system is faster than the Flynn's method, because the Flynn's method requires far more operations to unwrap a 2D phase map, as you can see in Ref. [8].The computational time taken for our recursive phase unwrapping system was 0.01 seconds, while the computational time taken by the Flynn's algorithm was 4.48 seconds.Both algorithms were programmed in C -language and compiled in a 64bit CPU.Fig. 2. In (a) we have the experimental wrapped phase used as input for the phase unwrapping systems, in (b) shows the recovery phase using the Flynn's phase unwrapping method (c) shows the recovery unwrapped phase using our proposed phase unwrapping recursive system.

Conclusions and commentaries
We have shown a robust to noise, two-dimensional simultaneous phase unwrapping and lowpass filtering system.It means that the presented recursive system unwraps and smoothes the phase map in the same process.Actually, this is a 2D Infinite Impulse Response (IIR) linear system.This is a very fast phase unwrapping system; its computational time for a 512 × 512 phase map was 0.01 seconds The presented phase unwrapping system is able to recover the unwrapped phase following a simple row by row scanning strategy with no previous labeling of phase inconsistencies as needed with the branch cut methods.
This work extends our previous method [9] by showing the 2D stability and frequency analysis of a new recursive phase unwrapping system.This frequency and stability analysis were not given before; this is the first time that a 2D phase unwrapping system is synthesized and gauged using frequency and stability criteria from linear systems theory.The 2D recursive unwrapping system shown here, is a substantial improvement of our previous 2D recursive system [9], because in our new case we conceptually decompose the 2D recursive filter [9] into a sum of a predictor and a corrector.This decomposition substantially improves the understanding of the inner workings of our new recursive phase unwrapping estimator.Finally our new 2D recursive phase unwrapping system demonstrates how such a system may be regarded, analyzed and gauged as a linear Infinite Impulse Response (IIR) filter.
Because some phase differences in Eq. ( 10) are taken at a distances of 2 pixels, our system is theoretically limited to spatial frequencies lower than π/2 radians.However, in practice our recursive phase unwrapping system has not this limit due to the recursive inertia of the previously unwrapped pixels already processed.

1 .
The Panels (a), (b) and (c) shows three possible neighborhood configurations found in a sequential scanning around (x, y), Fig.1(a) presents a single previously unwrapped pixel, Fig.1(b) presents 4 unwrapped pixels and Fig.1(c) presents 8.For these cases, the kernel is adapted as

Fig. 1 .
Fig. 1.In (a), (b) and (c), it is shown three different neighborhoods that we can find following a sequential scanning strategy.In (d), (e) and (f), we show their 2D power spectrums of its frequency responses for τ = 0.13.The pixel (x, y) visited is at the center, and the power spectrums are shown between the range (−π, π) in both frequency directions.