Noise robust linear dynamic system for phase unwrapping and smoothing

Phaseunwrapping techniques remove the modulus 2 π ambiguities of wrapped phase maps. The present work shows a first-order feedback system for phase unwrapping and smoothing. This system is a fast phase unwrapping system which also allows filtering some noise since in deed it is an Infinite Impulse Response (IIR) low-pass filter. In other words, our system is capable of low-pass filtering the wrapped phase as the unwrapping process proceeds. We demonstrate the temporal stability of this unwrapping feedback system, as well as its low-pass filtering capabilities. Our system even outperforms the most common and used unwrapping methods that we tested, such as the Flynn’s method, the Goldstain’s method, and the Ghiglia least-squares method (weighted or unweighted). The comparisons with these methods show that our system filters-out some noise while preserving the dynamic range of the phase-data. Its application areas may cover: optical metrology, synthetic aperture radar systems, magnetic resonance, and those imaging systems where information is obtained as a demodulated wrapped phase map. © 2011 Optical Society of America OCIS codes: (100.5088) Phase unwrapping; (100.2650) Fringe Analysis. References and links 1. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72, 156 (1982). 2. J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. 13, 2693 (1974). 3. Q. Kemao, “Two-dimensional windowed fourier transform for fringe pattern analysis: Principles, applications and implementations,” Optics And Lasers In Engineering 45, 304–317 (2007). 4. M. Servin, J. C. Estrada, and J. A. Quiroga, “The general theory of phase shifting algorithms,” Opt. Express 17, 21867–21881 (2009). 5. L. N. Mertz, “Speckle imaging, photon by photon,” Appl. Opt. 18, 611–614 (1979). 6. B. R. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,” J. Opt. Soc. Am. 69, 393–399 (1979). 7. K. A. Stetson, J. Wahid, and P. Gauthier “Noise-immune phase unwrapping by use of calculated wrap regions,” Appl. Opt.36, 4830–4838 (1997). 8. K. Itoh, “Analysis of the phase unwrapping algorithm,” Appl. Opt. 21, 2470–2470 (1982). 9. T. R. Judge and P. J. Bryanston-Cross, “A review of phase unwrapping techniques in fringe analysis,” Optics and Lasers in Engineering 21, 199 – 239 (1994). 10. D. C. Ghiglia and M. D. Pritt, Two-dimensional Phase Unwrapping; Theory, Algorithms, and Software (Wil yInterscience, 1998). #138076 $15.00 USD Received 11 Nov 2010; revised 15 Jan 2011; accepted 25 Jan 2011; published 3 Mar 2011 (C) 2011 OSA 14 March 2011 / Vol. 19, No. 6 / OPTICS EXPRESS 5126 11. D. C. Ghiglia, G. A. Mastin, and L. A. Romero, “Cellular-automata method for phase unwrapping,” J. Opt. Soc. Am. A 4, 267–280 (1987). 12. J. M. Huntley, “Noise-immune phase unwrapping algorithm,” Appl. Opt. 28, 3268–3270 (1989). 13. D. C. Ghiglia 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). 14. J. L. Marroquin and M. Rivera, “Quadratic regularization functionals for phase unwrapping,” J. Opt. Soc. Am. A 12, 2393–2400 (1995). 15. J. L. Marroquin, M. Tapia, R. Rodriguez-Vera, and M. Servin, “Parallel algorithms for phase unwrapping based on Markov random field models,” J. Opt. Soc. Am. A 12, 2578–2585 (1995). 16. K. M. Hung and T. Yamada, “Phase unwrapping by regions using least-squares approach,” Optical Engineering 37, 2965–2970 (1998). 17. V. V. Volkov and Y. Zhu, “Deterministic phase unwrapping in the presence of noise,” Opt. Lett. 28, 2156–2158 (2003). 18. 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). 19. J. G. Proakis and D. G. Manolakis, Digital Signal Processing. Principles, Algorithms, and Applications (PrenticeHall, October 5, 1995), 3rd ed.


Introduction
There are several applications where information of interest is phase modulated within the received signals and recovered as a wrapped phase.For example, in optical metrology the information of testing events is phase modulated within the interferograms by the use of optical interferometers [3].Other examples where information is phase modulated occur in synthetic aperture radars, magnetic resonance images, acoustic imaging, and X-ray crystallography.In all these cases the phase information is generally obtained modulus 2π radians, and this is called the wrapped demodulated phase.The wrapping process occurs because the searched phase information is recovered containing phase jumps (modulo 2π) when the phase information goes below −π radians or above π radians.In optical interferometry most techniques give the recovered phase modulus 2π [1][2][3][4].Under noise-free conditions, the phase unwrapping procedure may be performed by sequential line integration of the phase gradient [5,6,[8][9][10].However, the presence of high noise introduces spurious phase jumps when using these line integration techniques.Therefore more robust implementations of phase unwrapping methods that tolerates the presence of noise have been proposed [7,[11][12][13][14][15][16][17][18].However, most of them are numerically-heavy iterative methods, and/or complex ad-hoc implementations are required.Hence, a much faster, and robust phase unwrapping system is necessary.
In this paper we present a first-order dynamic system for phase unwrapping, which is very simple, stable, and fast.This is a Infinite Impulse Response (IIR) first-order dynamic system.In the following sections, we show the analysis and design of our phase unwrapping system using well known techniques from Digital Signal Processing (DSP).We will make an analysis of its impulse response to show its low-pass filtering properties in the frequency domain.After this, some tests and comparisons against well established techniques will be shown, such as the Flynn's method, the Goldstein's method, and the Ghiglia's leas-squares method.All of them are completely exposed in the Ghiglia's book [10].Finally, we draw some commentaries and conclusions.

The phase unwrapping system: analysis, and design
We start with the one-dimensional (1D) view of our dynamic phase unwrapping system.In a continuous domain the unwrapping process can be stated as a line integration procedure [8,10]: Fig. 1.Block diagram of the first order phase unwrapping system (4).The block system z −1 delays its input one sample.
where φ (t) is the output (unwrapped phase) obtained from a previously unwrapped phase φ (t 0 ), plus the integration of φ ′ w (t) = W dφ w (t) dt ; being φ w (t) the wrapped phase that we want to unwrap.The function W [•] is the modulus 2π operator; the main use of this operator in the phase derivative (or differences in a discreet domain) is to wrap or remove the outliers of the difference operator generated by the 2π phase jumps of the wrapped phase.This 1D line integration may be also stated as the solution of the following continuous system which is a first order differential equation that may be solved numerically in the discrete domain using the following recursive difference system: where φ (n) is the output unwrapped phase and φ (n − 1) its unwrapped previous value.However, in 2D we know from previous works that with noisy data the simple line integration (3) fails most of cases.This is because the phase difference operation acts as a high-pass filter, augmenting even more the levels of noise [19].We introduce in Eq. (3) two significant changes as explained below, obtaining our proposed difference system as φ The two introduced changes are: the wrapped data φ w (n − 1) at the previous site n − 1 by the unwrapped data φ (n − 1) at the same site n − 1, and the parameter τ.As will be shown below, the parameter τ is a regularization parameter that filters out most of the noise.In addition to this, we see that the wrapping operator W [•] takes the difference between a noisy data and a noisereduced estimation.Thus, the probability of generating spurious phase jumps in the unwrapped phase is much lower.A block diagram of our system (Eq.4) is depicted in Fig. 1.Now, considering that the digital modulating phase generates a well sampled sinusoidal signal, it means that the absolute value of the phase differences (local frequency) is less than π radians since this is the Nyquist limit, or in other words, the maximum frequency rate that can be sampled without aliasing.Then, just for analysis purposes we can write the following: and the linear system of Eq. ( 4) can be rewritten as To solve this difference equation we can use the z-transform, and obtain that the response of our dynamic phase unwrapping (4) for any arbitrary input φ w (n) is: ´ µ ´ µ Fig. 2. Impulse response, and frequency response of the dynamic phase unwrapping system (4), respectively.
where φ (−1) is the unwrapped phase initial value.This is a causal system, and from its response we have that the system's impulse response is where u(n) is the unitary step function.Now, to demonstrate the stability of this system we must prove that its impulse response is absolutely summable (see Ref. [19]).Thus, we can proof that the following is true: since this is a geometric series, proving in this way that the system's impulse response is absolutely summable, and thus we can say that the system (4) is an stable system.To see the frequency response of the system (4) we can take directly the Fourier transform of its impulse response and show that its frequency response is being F {•} the Fourier transform operator.In Fig. 2, panel (a) shows the graph of the system's impulse response, an panel (b) shows the graph of its frequency response.Note from this figure, that the systems impulse response is attenuated by τ, but its frequency response varies its wideband with τ and it has always a maximum of one in magnitude.Thus, the frequency response of the system behaves like a low-pass filter and its bandwidth is controlled by the parameter τ.

Two-dimensional extension
In the previous analysis we have designed a 1D phase unwrapping system (Eq.( 5)) and we have shown that this system is an stable system which low-pass filter the wrapped input.Now, we extend our 1D system to a two-dimensional (2D) unwrapping system.In this case we have a 2D neighborhood whose sites are in different directions.Trying to follow the same ideas from the 1D view, let us denote the vector r = (x, y) as the 2D site being unwrapped; being x and y are integers.Now let us denote with m(r) = {r 0 , r 1 , . . ., r N−1 } to the previously unwrapped sites around r.In Fig. 3, we show a graphic representation of the spatial support that we take around the site r.In that figure, Ω is a neighborhood around the site r.There, we show with dark points the set m(r) ⊂ Ω of previously unwrapped sites, and with a low gray, the sites that Fig. 3. Particular 2D graphic representation of neighborhood Ω around the site r that is going to be unwrapped.
are not unwrapped yet.Thus, following this scheme, a 2D extension for the 1D dynamic system (4) is Where |m(r)| denote the cardinality of the set m(r) ⊂ Ω.For example, counting the dark points of Fig. 3, we see that for that particular case the dimension of m(r) is |m(r)| = 4. Thus, the 2D dynamic phase unwrapping system (11) can be seen as the mean of the unwrapped phase obtained with the 1D system (4) from all possible directions around r.The unwrapping system (11) must process each site sequentially starting from an initial given site or point.At the starting point, we use the wrapped phase of that point as the initial unwrapped phase.To visit sequentially all sites of the 2D phase domain, one can use any 2D scanning strategies.However, when the phase domain has holes, or the phase domain is not a regular form, we recommend the standard flood-fill scanning used to color sequentially connected regions.

Tests and results
The following test compares our unwrapping dynamic system (11) with the Flynn's method, the Goldstein's method, and the least-squares method.These methods are well described in the Ghiglia's book of phase unwrapping [10].In practice these are some of the most used methods.The Fig. 4 shows the results.In Fig. 4.A it is shown the quadratic simulated phase that we numerically generated, Fig. 4.B shows the simulated phase plus noise with a normal distribution with mean zero, an variance of 1.3 radians.Fig. 4.C shows the unwrapped phase surface obtained with the herein proposed system (the dynamic phase unwrapping system), Fig. 4.D shows the unwrapped surface with the Flynn's method, Fig. 4.E with the Goldstein's method, and Fig. 4.F the Ghiglia's least-squares method.In all these cases, the input is the noisy phase shown in Fig. 4.B.For illustration purposes, and for comparing the output of each method with its input, we have projected bellow each phase surface the image of the modulus 2π wrapped phase corresponding to each unwrapped phase surface.The scanning strategy used by our phase unwrapping system was the simple row by row scanning, and we set τ = 0.2.In those results, we see that the Flynn's method obtains better results than the Goldstein's and the least-squares method.Also, we can see that our phase unwrapping system obtains a better result than the Flynn's method, and we see that it reduces some levels of noise.Counting the rings of the iso-phase contours from the projected wrapped phase under each unwrapped phase surface, we can see that only the unwrapped phase of the least-squares method has reduced the dynamic range.This is a drawback of the least-squares method that is already known in this area, and it is important to remark that our phase unwrapping system preserves the dynamic range and reduces the levels of noise.In these tests, the computational time spent for our dynamic phase unwrapping system was of 0.064 seconds, for the Flynn's method was of 0.858 seconds, for the Goldstein method was of 0.068, and for the least-squares method was of 0.105 seconds (by means of direct transforms [10]); with a phase map of 256 × 256 pixels.All these test were made in a computer with a 32-bits CPU and 2GB of memory.This times shows that our phase unwrapping system is as fast as other tested methods and even reduces some levels of noise.Finally, in Fig. 5 we show a result from an experimentally obtained wrapped phase image.We show in Fig. 5.A the experimental wrapped phase image, in Fig. 5.B the unwrapped phase surface obtained with our dynamic phase unwrapping system, and in Fig. 5.B with the Flynn's method.We compare only with the Flynn's method since from the previous test it was the one with best results than the others (the Goldstein, and least-squares method).To obtain the result shown in Fig. 5.A, we pass two times our phase unwrapping system; the first time using τ = 0.6, and the second time using τ = 0.2 with the output of the first pass as input.We can see here that our phase unwrapping system is able to reduce considerably the levels of noise without practically affect the dynamic range of the unwrapped phase.This because the projected wrapped phase of Fig. 5.B has the same iso-phase contours as the input wrapped phase used (Fig. 5.A).

Conclusions
As it was shown, our phase unwrapping system filters directly the wrapped phase while the unwrapping process proceeds.The bandwidth of our phase unwrapping system is controlled with a parameter 0 < τ < 1.The spatial support of the 2D extension allow us to unwrap the phase following arbitrary scanning strategies, however, for the results shown here we used a simple row by row scanning strategy.Nevertheless, we implemented the flood-fill scanning strategy for the phase unwrapping processing obtaining the same results.As reported in the last paragraph of the section 3, our phase unwrapping system takes around 0.065 seconds of computational time to unwrap an filter the wrapped phase, whereas the other phase unwrapping methods that we programmed took approximately the same time but these did not removed any level of noise.Hence, we consider that our phase unwrapping system is a noise tolerant fast phase unwrapping system.
Important remarks: As it was shown in section 2, the phase unwrapping system developed in this paper is a linear IIR causal system [19].This particular property reduces the noise levels while the phase unwrapping process proceeds.However, recalling from digital system's theory, causal IIR systems are systems with past memory [19].This memory is the initial value condition that results from the solution of the system's difference equation (see Eq. ( 7)).When this memory is initially zero, it is said that the system is relaxed.Hence, if we start the processing with a relaxed system the first processed values will underestimated by the influence of this memory-less system.In our case, we take the value of the input at the starting point as initial memory.Even thus, one can not avoid completely the initial effect that has a linear relaxed system with the first processed samples.Nevertheless, numerically we realized that starting the system with the value of the input at the starting point is better than using a completely relaxed system.The reader may experiment the meaning of this words making a simple 1D test and programming the system of Eq. 4. By passing two times our phase unwrapping system as described in section 3 for the last test, the effect of the system with the first processed samples is reduced since in the second pass the system is not relaxed, it will have the memory from the first pass.

Fig. 4 .
Fig. 4. Phase unwrapping of simulated data.A) is the simulated wrapped phase, B) the simulated wrapped phase plus noise with normal distribution, mean zero and variance of 1.3 radians, C) is the unwrapped surface with the herein proposed method, D) with the Flynn's method, E) with the Goldstein's method, and F) with the leas-squares method.To locate the dynamic range, bellow each unwrapped surface is projected the image corresponding to the modulus 2π wrapped phase of each unwrapped phase surface.

B
) Proposed method C) Fynn's method A) Experimental phase image

Fig. 5 .
Fig. 5. Test using an experimentally obtained phase map.The phase map dimensions are 512 × 512.A) shows the experimental phase map image, B) the unwrapped phase surface obtained with our dynamic phase unwrapping system, and C) with the Flynn's method.