## Abstract

An algorithm for dynamic phase retrieval in temporal speckle pattern interferometry using least squares method and windowed Fourier filtering is proposed. The least squares method is used to evaluate the phase change between two speckle patterns provided that the phase of either one speckle pattern has been estimated. The windowed Fourier filtering is used to eliminate the noise in the phase change. Based on these two techniques, the proposed algorithm determines the phase of the initial speckle pattern by phase shifting method at first, then the phase of the rest speckle patterns are retrieved by sequentially evaluating the phase changes between every two consecutive speckle patterns. The algorithm solves the problem of speckle decorrelation by refreshing the reference image frame by frame, and also avoids the problem of error accumulation during the reference image refreshing process by the windowed Fourier filtering. Two experimental results are presented to demonstrate the effectiveness and robustness of the proposed algorithm.

©2011 Optical Society of America

## 1. Introduction

Speckle pattern interferometry is a whole field, noncontact and high sensitivity interferometric technique appropriate for studies of objects with optically rough surfaces [1]. The method is based on the interference of a scattered light from the objects with a reference beam, which results in a random intensity field called a speckle pattern. The changes of the object state, for instance, deformation and displacement, are then measured by detecting the phase changes of the speckle pattern [2].

The most straightforward approach for dynamic phase retrieval from speckle patterns is to process each object state separately. A number of well-established phase retrieval methods for static or semi-static phenomena can be used [3]. Phase shifting method [4] is the most frequently used technique due to its high accuracy and easy implementation. This method requires at least three speckle patterns for a single state and is not suitable for dynamic processes. To cope with this problem, high speed temporal phase shifting method has been proposed and successfully implemented [5]. Spatial phase shifting method [6,7] is a good alternative which simultaneously acquires several phase shifting images for a single state. Both high speed temporal phase shifting and spatial phase shifting methods require complicated and usually more expensive optical setups. Fourier transform method [8] is able to retrieve the phase distribution from a single fringe pattern. This method requires a carrier to separate the spectrum of the fringe pattern in the frequency domain, which limits its measurement range. The phase tracking techniques [9,10] attempt to retrieve the phase distribution from a single fringe pattern even without a carrier. These techniques have no special requirements on the spectrum of the fringe pattern, but they often result in some nonlinear optimization problems and lead to complex and time-consuming computations.

In the past few years, temporal speckle pattern interferometry (TSPI) has been rapidly developed with emphasis on the relationship between object states. In TSPI, a sequence of speckle patterns is recorded throughout the entire measurement history of the object under study. Then the optical phase distributions coded in the temporal intensity are retrieved. The methods proposed for TSPI can be roughly categorized into two approaches. The first approach is transform-based [11–17]. Fourier transform [11,12], Hilbert transform [13,14], wavelet transform [15,16], and windowed Fourier transform [17] have been applied to the time axis. All of them suffer from the requirement of a temporal carrier which limits the measurement range. The second approach is reference-based [18–22]. The phase of the initial speckle pattern is first determined by the phase shifting method as a reference. The phase of the rest speckle patterns are evaluated through estimating the phase changes between the current speckle pattern and the reference by a least squares method [18–20] or a clustering algorithm [21,22]. This approach does not need a temporal carrier and therefore has a potential to be used more widely. The disadvantages of the current methods for this approach are that, during the measurement process, the background and the modulation of intensity should be estimated by the scanning phase technique [18,19], or these quantities are assumed to be constants [20–22]. Besides, the methods do not emphasize the problem of speckle decorrelation which plays an important role in speckle interferometry methods [1].

In this paper, an algorithm for phase retrieval in TSPI using least squares method and windowed Fourier filtering is proposed. The algorithm adopts the reference-based approach [18–22], i.e., determining the phase of the initial speckle pattern in advance and then evaluating the phase of the rest speckle patterns by estimating the phase changes. There are several innovations in the algorithm. First, the advanced iterative phase shifting method [23] is used to precisely estimate the phase of the initial speckle pattern without knowing the phase shift amounts. Second, a new least squares method is presented, which can simultaneously estimate the background, the modulation and the phase change. Third, the windowed Fourier filtering [24–26] is used to eliminate the noise in the phase change, which avoids the problem of error accumulation. Finally, the strategy of using refreshing reference image frame by frame solves the problem of speckle decorrelation. Based on these innovations, the proposed algorithm can successfully retrieve the phase distributions from the speckle patterns with the characteristics of (i) no requirement of the knowledge of phase shifts in the initial phase shifting stage; (ii) no requirement of the constant condition or pre-estimation of background intensity and modulation; (iii) no problems of speckle decorrelation; and (iv) no error accumulation during the solving process. Two experimental results are presented to demonstrate the effectiveness and robustness of the proposed algorithm.

## 2. Theory

In TSPI, a sequence of speckle patterns is recorded throughout the entire measurement process and can be expressed as

where $f\left(x,{t}_{i}\right)$ denotes the intensity distributions of the speckle pattern acquired at time ${t}_{i}$; $x=\left({x}_{1},{x}_{2}\right)$ is the pixel coordinate; $a\left(x,{t}_{i}\right)$ is the background; $b\left(x,{t}_{i}\right)$ is the modulation, $\phi \left(x,{t}_{i}\right)$ is the phase to be measured; $n\left(x,{t}_{i}\right)$ denotes the noise contained in the acquired speckle pattern which generally caused by the inherent speckle noise appeared in any techniques based on speckle pattern interference. *N* is the frame number of the recorded speckle patterns. Our task is to retrieve the phase term $\phi \left(x,{t}_{i}\right)$ from the speckle patterns $f\left(x,{t}_{i}\right)$. It can be achieved by the proposed algorithm which includes the following major components.

#### 2.1 Advanced iterative algorithm for determination of the phase of the initial speckle pattern

The initial state in most optical measurement can be seen as static or semi-static, so it is acceptable to use the phase shifting method to determine the phase distribution of the initial speckle pattern, $\phi \left(x,{t}_{0}\right)$. All existing phase shifting algorithms can be used [3]. The advanced iterative algorithm (AIA) proposed by Wang *et al.* [23] is selected, which can evaluate the phase through a few random phase shifting speckle patterns. To apply the AIA, four random phase shifting speckle patterns are acquired before the dynamic process starts. The speckle patterns can be expressed as

where ${f}_{k}\left(x,{t}_{0}\right)$ is the *k*-th phase shifting speckle pattern, ${\delta}_{k}$ is the corresponding random phase shift. The AIA solves the phase $\phi \left(x,{t}_{0}\right)$ by performing an iterative procedure. The details of the algorithm can be found in Ref. [23].

#### 2.2 Least squares method for the estimation of phase changes

The least squares method presented here simultaneously estimates the background, the modulation and the phase change between two speckle patterns provided that the phase of either one speckle pattern is estimated. Assuming the phase $\phi \left(x,{t}_{\alpha}\right)$ at time ${t}_{\alpha}$ has already been estimated, a speckle pattern $f\left(x,{t}_{\beta}\right)$ recorded at time ${t}_{\beta}$ can be rewritten as

where $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$ denotes the phase change between the two speckle patterns defined by

Note that the noise term in the speckle pattern $f\left(x,{t}_{\beta}\right)$ is omitted in Eq. (3) and the following equations for the sake of simplicity. Next a pixel at $\zeta =\left({\zeta}_{1},{\zeta}_{2}\right)$ is considered and its neighborhood is defined as${N}_{\zeta}=\left\{x|D\left(x,\zeta \right)\le \epsilon \right\}$, where $D\left(x,\zeta \right)=\mathrm{max}\left(\left|{x}_{1}-{\zeta}_{1}\right|,\left|{x}_{2}-{\zeta}_{2}\right|\right)$ and *ε* is a small positive number. In this neighborhood, since the speckle pattern $f\left(x,{t}_{\beta}\right)$ is spatially continuous, the background $a\left(x,{t}_{\beta}\right)$, the modulation $b\left(x,{t}_{\beta}\right)$ and the phase change $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$ can be approximated as constants,

With Eq. (5), $f\left(x,{t}_{\beta}\right)$ for $x\in {N}_{\zeta}$ can be approximated as follows:

where

It can be seen that there are three unknowns $a\left(\zeta ,{t}_{\beta}\right)$, $c\left(\zeta ,{t}_{\beta}\right)$, and $d\left(\zeta ,{t}_{\beta}\right)$ and *M* equations in Eq. (6), where *M* is the number of the pixels in ${N}_{\zeta}$. If $M\ge 3$, the unknowns can be solved in the least-squares sense. The residual error *r* accumulated from all the pixels in ${N}_{\zeta}$ can be written as

The least-squares criterion requires that

which yields

where

From Eqs. (10)–(13), the unknowns $a\left(\zeta ,{t}_{\beta}\right)$, $c\left(\zeta ,{t}_{\beta}\right)$, and $d\left(\zeta ,{t}_{\beta}\right)$ can be solved, and the phase change and the modulation can be subsequently determined by

The same procedure can be performed on every pixel of the speckle pattern. Therefore, the phase change of the entire speckle pattern, $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$, can be evaluated.

#### 2.3 Windowed Fourier filtering for elimination of noise in the phase changes

The phase change $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$ obtained above inevitably contains noise. In order to avoid the problem of error accumulation which will be discussed later in Section 2.4, it is necessary to filter the $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$ for noise elimination. The windowed Fourier filtering is used due to its effectiveness and accuracy, i.e., if the phase is locally quadratic, the windowed Fourier filtering does not distort the intrinsic signal [26]. According to [25], the phase change $\Delta {\phi}_{\beta ,\alpha}\left(x\right)$ is first converted into an exponential phase field $\Phi \left(x\right)$,

where $j=\sqrt{-1}$; next $\Phi \left(x\right)$ is transformed into the windowed Fourier space by performing a windowed Fourier transformation and the windowed Fourier spectrum is modified by a thresholding processing; last the filtered result $\overline{\Phi}\left(x\right)$ is obtained by performing an inverse windowed Fourier transformation on the modified spectrum, and the filtered phase change $\overline{\Delta {\phi}_{\beta ,\alpha}}\left(x\right)$ can be obtained as

More details about the windowed Fourier filtering can be found in Ref. [25].

#### 2.4 Proposed algorithm

As mentioned above, the phase change between two speckle patterns can be evaluated provided that the phase of either one speckle pattern is known. This conclusion can be reached only when the hypothesis expressed in Eq. (5) is satisfied. To make it possible, spatially, $\epsilon =1$ is used so that only the nearest 9 pixels are involved for the least squares method, which best satisfies the assumption of constant background and modulation in Eq. (5); temporally, ${t}_{\beta}-{t}_{\alpha}=1$ is used so that only the nearest frame is used for reference, which best satisfies the assumption of constant phase change in Eq. (5). The speckle decorrealtion [27] is automatically avoided since only the nearest frame is used for reference, i.e., the reference image is refreshed frame by frame. However, another problem of error accumulation [28,29] appears accompanying with the refreshing of the reference image due to the inevitable noise in the phase changes. To deal with this problem, the noise contained in the phase changes are eliminated by the windowed Fourier filtering as described in Section 2.3. The proposed integrated algorithm is detailed as follows:

- Step 1: set $n=0$; evaluate the phase $\phi \left(x,{t}_{0}\right)$ of the initial speckle pattern $f\left(x,{t}_{0}\right)$ by the AIA method;
- Step 2: $n=n+1$; estimate the phase change $\Delta {\phi}_{n,n-1}\left(x\right)$ between the
*n*-th speckle pattern and the (*n−*1)-th speckle pattern by the least squares method; - Step 3: obtain the filtered phase $\overline{\Delta {\phi}_{n,n-1}}\left(x\right)$ by performing the windowed Fourier filtering on $\Delta {\phi}_{n,n-1}\left(x\right)$ and retrieve the phase $\phi \left(x,{t}_{n}\right)$ of the speckle pattern $f\left(x,{t}_{n}\right)$ by calculating $\phi \left(x,{t}_{n}\right)=\phi \left(x,{t}_{n-1}\right)+\overline{\Delta {\phi}_{n,n-1}}\left(x\right)$;
- Step 4: repeat Step 2 and Step 3 until $n=N-1$.

In TSPI, the phase change between the current state and the initial state is often most interesting since it directly reflects the physical qualities under measurement, for instance, the displacement. The phase change between any speckle pattern, for example the *k*-th speckle pattern, and the initial speckle pattern can be calculated by

The phase result $\overline{\Delta {\phi}_{k,0}}\left(x\right)$ is automatically unwrapped because all the phase changes $\overline{\Delta {\phi}_{n,n-1}}\left(x\right)$ are usually small for an appropriate image acquiring frame rate.

## 3. Results and discussion

To verify the effectiveness and robustness of the proposed algorithm, two speckle pattern shearing interferometry experiments are performed. In the first experiment, the deformation of a fully clamped circular plate subjected to a central force is investigated. Before the force is applied, four random phase shifting speckle patterns are recorded. Then the force is applied to the plate, including several loading and unloading cycles. A sequence of speckle patterns is acquired at the frame rate of 30fps during the force acting process. The frame number is 117. The size of each speckle pattern is 366 × 371 pixels. Figures 1(a)-(c) show the speckle fringe patterns at three instances produced by subtracting the initial speckle pattern from them, respectively. The fringe patterns represent the contour of the derivative of the out-of-plane displacement of the plate. It can be seen that the quality of the fringe patterns is degraded by the speckle noise and the speckle decorrelation effect. It is difficult to retrieve the phase from the fringe patterns by conventional phase retrieval techniques, especially for the fringe pattern shown in Fig. 1(c). To employ the proposed algorithm, the phase of the initial speckle pattern is first determined from the four random phase shifting speckle patterns by the AIA, then the phase of other speckle patterns are retrieved by repeating Steps 2 and 3 described in Section 2.4 until all the speckle patterns are processed. Figures 1(d)-(f) show the obtained phase changes distributions calculated by Eq. (18) corresponding to the speckle fringe patterns shown in Figs. 1(a)-(c), respectively, which are seen to be satisfactory. The parameters of the windowed Fourier filtering used in this paper are ${\sigma}_{x}={\sigma}_{y}=20$, ${\xi}_{l}={\eta}_{l}=-0.2$, ${\xi}_{0}={\eta}_{0}=0.1$, ${\xi}_{h}={\eta}_{h}=0.2$, $thr=15$, respectively. The detailed meaning of these parameters can be found in [25]. Note that all phase results obtained in this paper are automatically unwrapped and are rewrapped for displaying purpose. A multimedia file linked to Fig. 1 (Media 1) gives the animation of the fringe patterns (left) and the phase results obtained by the proposed algorithm (right). The animation file is compressed due to the file size limitation. The original version of the animation can be found in [30]. To further understand the proposed algorithm, some comparative results are given. Figures 1(g)-(i) show the phase results obtained by taking the initial speckle pattern as a fixed reference image to estimate the phase of other speckle patterns. Speckle decorrelation is clearly observed in Fig. 1(i) and the fringe structures are partly destroyed. This shows the necessity of a refreshing reference speckle pattern. Figures 1(j)-(l) show the phase results obtained by refreshing the reference speckle pattern frame by frame but do not perform the windowed Fourier filtering. It can be seen that the real phase information is almost totally lost in the obtained results. Comparing with Figs. 1(d-f), this result reveals the effectiveness of the windowed Fourier filtering. All calculations carried out in this paper are performed on a personal Pentium Dual E8400 computer with 3.0GHz main frequency by Matlab programming. The total time for the phase retrieval in the above experiment is 379s corresponding to about 3s for each speckle pattern.

Another experiment is about nondestructive testing. A honeycomb aluminum panel is subjected to a thermal loading and shearography is used for disbond localization and evaluation. Four phase shifting speckle patterns are acquired before the thermal loading and then 171 frames of speckle patterns are sequentially recorded at the frame rate of 30fps accompanying with the thermal loading process. The size of each speckle pattern is 551 × 551 pixels. Figures 2(a)-(c) show the shearographic fringe patterns at three different instances, and Figs. 2(d)-(f) give the corresponding phase distributions. It can be seen that the disbond defects in the panel are distinctly shown in the phase maps. A multimedia file linked to Fig. 2 (Media 2) displays the animation of the fringe patterns (left) and the obtained phase results (right). The original version of the animation can also be found in [30]. The total time for the phase retrieval in this experiment is 1332s corresponding to about 8s for each speckle pattern.

## 4. Conclusions

An algorithm for phase retrieval in temporal speckle pattern interferometry using least squares method and windowed Fourier filtering is proposed. The algorithm can estimate the phase from a single speckle pattern provided the phase of the initial speckle pattern can be determined beforehand. A least squares method is presented to evaluate the phase change between two speckle patterns and the windowed Fourier filtering is employed to eliminate the noise contained in the phase change. Based on these two methods and another strategy of refreshing the reference image frame by frame, the phase distributions of the speckle patterns can be successfully retrieved without the influences of the speckle decorrelation and the error accumulation. Two experimental results are presented to demonstrate the effectiveness and robustness of the proposed algorithm. The computational efficiency of the proposed algorithm can be further enhanced with the technique of GPU [31] and we expect to achieve a real-time dynamic phase retrieval algorithm in the future studies.

## Acknowledgments

This work was partially supported by the Singapore MOE Academic Research Fund Tier 1 (RG11/10) and the National Natural Science Foundation of China (NSFC) (No. 10902066). The authors thank Prof. Dongsheng Zhang for his experimental support.

## References and links

**1. **R. S. Sirohi, Speckle Metrology (Marcel Dekker, New York, 1993).

**2. **R. Jones and C. Wykes, *Holographic and Speckle Interferometry* (Cambridge University Press, London, 1983).

**3. **D. Malacara, M. Servin, and Z. Malacara, Interferogram analysis for Optical Testing, 2th ed. (Marcel Dekker, 2003).

**4. **K. Creath, “Phase-shifting speckle interferometry,” Appl. Opt. **24**(18), 3053–3058 (1985). [CrossRef] [PubMed]

**5. **J. M. Huntley, G. H. Kaufmann, and D. Kerr, “Phase-shifted dynamic speckle pattern interferometry at 1 kHz,” Appl. Opt. **38**(31), 6556–6563 (1999). [CrossRef] [PubMed]

**6. **A. J. Haasteren and H. J. Frankena, “Real-time displacement measurement using a multicamera phase-stepping speckle interferometer,” Appl. Opt. **33**(19), 4137–4142 (1994). [CrossRef] [PubMed]

**7. **J. Millerd, N. Brock, J. Hayes, M. North-Morries, B. Kimbrough, and J. Wyant, “Pixelated phase-mask dynamic interferometers,” in *Fringe 2005. The 5th International Workshop on Automatic Processing of Fringe Patterns*, W. Osten, ed. (Springer, 2006), Session 5, pp. 640–647.

**8. **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**(1), 156–160 (1982). [CrossRef]

**9. **M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. **36**(19), 4540–4548 (1997). [CrossRef] [PubMed]

**10. **H. Wang and Q. Kemao, “Frequency guided methods for demodulation of a single fringe pattern,” Opt. Express **17**(17), 15118–15127 (2009). [CrossRef] [PubMed]

**11. **C. Joenathan, B. Franze, P. Haible, and H. J. Tiziani, “Speckle interferometry with temporal phase evaluation for measuring large-object deformation,” Appl. Opt. **37**(13), 2608–2614 (1998). [CrossRef] [PubMed]

**12. **G. H. Kaufmann and G. E. Galizzi, “Phase measurement in temporal speckle pattern interferometry: comparison between the phase-shifting and the Fourier transform methods,” Appl. Opt. **41**(34), 7254–7263 (2002). [CrossRef] [PubMed]

**13. **V. D. Madjarova, H. Kadono, and S. Toyooka, “Dynamic electronic speckle pattern interferometry (DESPI) phase analyses with temporal Hilbert transform,” Opt. Express **11**(6), 617–623 (2003). [CrossRef] [PubMed]

**14. **F. A. Marengo Rodriguez, A. Federico, and G. H. Kaufmann, “Hilbert transform analysis of a time series of speckle interferograms with a temporal carrier,” Appl. Opt. **47**(9), 1310–1316 (2008). [CrossRef] [PubMed]

**15. **X. Colonna de Lega, “Processing of non-stationary interference patterns: adapted phase shifting algorithms and wavelet analysis. Application to dynamic deformation measurements by holographic and speckle interferometry,” Ph.D. dissertation 1666 (Swiss Federal Institute of Technology, 1997).

**16. **A. Federico and G. H. Kaufmann, “Robust phase recovery in temporal speckle pattern interferometry using a 3D directional wavelet transform,” Opt. Lett. **34**(15), 2336–2338 (2009). [CrossRef] [PubMed]

**17. **Y. Fu, R. M. Groves, G. Pedrini, and W. Osten, “Kinematic and deformation parameter measurement by spatiotemporal analysis of an interferogram sequence,” Appl. Opt. **46**(36), 8645–8655 (2007). [CrossRef] [PubMed]

**18. **T. E. Carlsson and A. Wei, “Phase evaluation of speckle patterns during continuous deformation by use of phase-shifting speckle interferometry,” Appl. Opt. **39**(16), 2628–2637 (2000). [CrossRef] [PubMed]

**19. **W. An and T. E. Carlsson, “Speckle interferometry for measurement of continuous deformations,” Opt. Lasers Eng. **40**(5-6), 529–541 (2003). [CrossRef]

**20. **L. Bruno and A. Poggialini, “Phase shifting speckle interferometry for dynamic phenomena,” Opt. Express **16**(7), 4665–4670 (2008). [CrossRef] [PubMed]

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

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

**23. **Z. Wang and B. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. **29**(14), 1671–1673 (2004). [CrossRef] [PubMed]

**24. **Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. **43**(13), 2695–2702 (2004). [CrossRef] [PubMed]

**25. **Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: Principles, applications and implementations,” Opt. Lasers Eng. **45**(2), 304–317 (2007). [CrossRef]

**26. **Q. Kemao, H. X. Wang, and W. J. Gao, “Windowed Fourier transform for fringe pattern analysis: theoretical analyses,” Appl. Opt. **47**(29), 5408–5419 (2008). [CrossRef] [PubMed]

**27. **C. Joenathan, P. Haible, and H. J. Tiziani, “Speckle interferometry with temporal phase evaluation: influence of decorrelation, speckle size, and nonlinearity of the camera,” Appl. Opt. **38**(7), 1169–1178 (1999). [CrossRef] [PubMed]

**28. **A. Svanbro, J. M. Huntley, and A. Davila, “Optimal re-referencing rate for in-plane dynamic speckle interferometry,” Appl. Opt. **42**(2), 251–258 (2003). [CrossRef] [PubMed]

**29. **A. Davila, J. M. Huntley, G. H. Kaufmann, and D. Kerr, “High-speed dynamic speckle interferometry: phase errors due to intensity, velocity, and speckle decorrelation,” Appl. Opt. **44**(19), 3954–3962 (2005). [CrossRef] [PubMed]

**30. **Q. Kemao, “Windowed Fourier transform (WFT) for fringe pattern analysis,” http://www3.ntu.edu.sg/home/mkmqian/.

**31. **W. J. Gao, N. T. T. Huyen, H. S. Loi, and Q. Kemao, “Real-time 2D parallel windowed Fourier transform for fringe pattern analysis using Graphics Processing Unit,” Opt. Express **17**(25), 23147–23152 (2009). [CrossRef] [PubMed]