## Abstract

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

## 1. 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.

## 2. 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]:

*ϕ̂*(

*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–12]. The frequency and impulse response of this linear low-pass filter are:

*τ*∈ (0,2) [10,11]. To transform Eq. (1) as a phase unwrapping system, let us rewrite it as:

*τ*, and the previous estimation

*ϕ̂*(

*n*– 1) in the phase difference, Eq. (3) looks very similar to the following basic line unwrapping algorithm [2]:

*ϕ*(

*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:

## 3. 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:

*ϕ̂*(

*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.

#### 3.1. 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:

*ϕ̂*(

*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. 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

*L*

^{1}norm is ‖

*S*‖

_{1}= 1, ‖

*S*‖

_{1}= 4 and ‖

*S*‖

_{1}= 8, respectively.

#### 3.2. 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:

*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*).

## 4. 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
$\left({z}_{x}^{-1},{z}_{y}^{-1}\right)\in {\overline{U}}^{2}$, where *Ū*^{2} is the unit bidisc defined as

*z*-transfer function of the linear system in Eq. (7) (which includes Eq. (8) and Eq. (10)) we obtain the following

*Ū*

^{2}is

*τ*‖

*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 $\frac{1}{4}>\tau >0$.

## 5. 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*and

^{iu}*z*=

_{y}*e*, being

^{iv}*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:

*z*-transform Eq. (12), using

*τ*= 0.13.

## 6. 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–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 14*n* 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.

## 7. Conclusions and commentaries

We have shown a robust to noise, two-dimensional simultaneous phase unwrapping and low-pass 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.

## References and links

**1. **K. Itoh, “Analysis of the phase unwrapping algorithm.” Appl. Opt. **21**, 2470 (1982). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef]

**5. **J. L. Marroquin and M. Rivera, “Quadratic regularization functionals for phase unwrapping,” J. Opt. Soc. Am. A **12**, 2393–2400 (1995). [CrossRef]

**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). [CrossRef]

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

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

**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). [CrossRef] [PubMed]

**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). [CrossRef]