Statistical phase-shifting step estimation algorithm ... - Semantic Scholar
Statistical phase-shifting step estimation algorithm based on the continuous wavelet transform for high-resolution interferometry metrology Bicheng Chen and Cemal Basaran* Electronic Packaging Laboratory, State University of New York at Buffalo, Buffalo, New York 14260, USA *Corresponding author: [email protected] Received 1 November 2010; accepted 13 December 2010; posted 21 December 2010 (Doc. ID 137421); published 31 January 2011
We propose a statistical phase-shifting estimation algorithm for temporal phase-shifting interferometry (PSI) based on the continuous wavelet transform (CWT). The proposed algorithm explores spatial information redundancy in the intraframe interferogram dataset using the phase recovery property on the power ridge of the CWT. Despite the errors introduced by the noise of the interferogram, the statistical part of the algorithm is utilized to give a sound estimation of the phase-shifting step. It also introduces the usage of directional statistics as the statistical model, which was validated, so as to offer a better estimation compared with other statistical models. The algorithm is implemented in computer codes, and the validations of the algorithm were performed on numerical simulated signals and actual phase-shifted moiré interferograms. The major advantage of the proposed algorithm is that it imposes weaker conditions on the presumptions in the temporal PSI, which, under most circumstances, requires uniform and precalibrated phase-shifting steps. Compared with other existing deterministic estimation algorithms, the proposed algorithm estimates the phase-shifting step statistically. The proposed algorithm allows the temporal PSI to operate under dynamic loading conditions and arbitrary phase steps and also without precalibration of the phase shifter. The proposed method can serve as a benchmark method for comparing the accuracy of the different phase-step estimation methods. © 2011 Optical Society of America OCIS codes: 050.5080, 000.5490, 100.5070, 100.7410, 120.2880, 120.3180.
1. Introduction
Phase-shifting interferometry (PSI) is an important technique in high-precision optical metrology. Temporal PSI was first introduced by Bruning et al. [1] to overcome the λ=2 sensitivity limits of the fringe counting methods in the interferometry for the application in the surface quality measurement of the optical components. In the Bruning et al. paper, error sources are enumerated, including electrical noise in the detectors, mechanical drifts, atmospheric turbulence and the noise in the laser source. It is also emphasized that the error can be significant when the sensitivity reaches to λ=100. The arctangent function is used in Bruning’s paper to reconstruct the 0003-6935/11/040586-08$15.00/0 © 2011 Optical Society of America 586
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
phase information. The least-squares estimation (LSE) method was first introduced by Morgan [2] to minimize the intensity errors nonlinearly between the imaged intensity map and the estimated one on the parameter of the measured phase. The commonly used version of LSE was established by Greivenkamp [3]. Greivenkamp uses the same interference equation but minimizes the intensity errors according to three parameters (background, cosine, and sine function of the phase shifter). Greivenkamp’s LSE is linear and removes several restrictions when designing the phase steps. Lai and Yatagai apply additional interference setups to detect the shifted phase, while using Fourier transform methods [4]. Larkin added preprocessing and postprocessing to improve the robustness of Lai and Yatagai’s Fourier-transform-based method [5]. Farrell and Player once again pointed out that
the error in the phase steps can cause a distortion of the results, and they suggested an interpixel ellipse fitting estimation method [6]. The phase step can be unknown and unequal. Wei et al. gave a leastsquares algorithm based on the interpixel relationship by means of a simple intensity transform between the sampled interferograms [7]. Iterative LSE was proposed by Kong and Kim, and considers both the phase and phase shift as parameters [8]. The iterative method and its application was also documented by Rivera et al. [9]. Based on the initially induced phase step (e.g., from the linearity table of piezoelectric actuator), Greivenkamp’s LSE is first used to estimate the wavefront phases, and a modified Greivenkamp linear LSE is used to estimate the phase step afterward. The iterative process is therefore built on these two LSE models to obtain the final wavefront phases upon the convergence condition. The idea is generalized and extended by Wang and Han [10]. Another phase-shift extraction method, based on the assumed object wave, is proposed by Cai et al. [11,12]. The Cai et al. method uses the mean of the differences between two interferograms to calculate the phase difference. An iterative process is available using the estimated phase steps to reconstruct the phase distribution iteratively. Xu et al. proposed a noniterative version based on the same principle [13]. The methods reviewed above can be divided into two categories: temporal methods and temporalspatial methods. The temporal methods rely on an established phase step to give pixelwise estimations on the phase distribution. The temporal-spatial methods use the spatial data redundancy to first estimate the phase step and then to reconstruct the phase distribution from the temporal samples [5]. The proposed algorithm in this paper belongs to the temporal-spatial method. The flow chart of the proposed algorithm is shown in Fig. 1. Compared with the existing algorithms, the proposed algorithm utilizes the continuous wavelet transform (CWT) to estimate the wavefront phases first, which is operated using spatial data correlation within an interferogram frame. Later, the phase steps are estimated statistically upon the phase difference distributions between temporal sampled interferogram frames. The parameter estimation part of the algorithm was demonstrated here using the maximum likelihood estimation (MLE). However, any statistical parameter estimation method or nonparametric estimation can be used instead. The proposed method allows an arbitrary phase step but does not involve any iterative process, which differs from most of the existing temporal-spatial methods. Furthermore, the proposed algorithm is a framework for building a robust phase-shifting algorithm, which can allow the PSI to operate under a dynamic loading condition. According to our literature survey, this is the first time that the CWT is reported to be used to help in the estimation of the phase-shifting steps.
Fig. 1. Flow chart of the proposed phase-shifting step estimation algorithm.
2. Interferometry Phase Least-Squares Estimation
The structure and notation of the interference equation are introduced in this section. They are critical elements to get an understanding of the proposed algorithm. This section also lays out some essential 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
587
concepts for the proposed algorithm, especially on the formulation and parameterization of the LSE used in the field of phase-shift estimations. In common circumstances, the temporal PSI samples multiple intensity values on the same spatial point at different times and with shifted phase steps. For the detector using a charge coupled device (CCD), each sampling is an accumulation process of the electrical charges that are proportional to the intensity of the incident photon flux integrated over the exposure time. The reconstruction of the temporal phaseshifted interferograms using the least-squares formula is widely used [3,14]. This serves in the last step in the proposed algorithm. The interference equation for the temporal PSI can be written as [14] I k ði; jÞ ¼ aði; jÞ þ bði; jÞ cos½ϕði; jÞ þ Δϕk ;
ð1Þ
where k is the number of phase shifting k ¼ 0, 1, …, K − 1 (K is the total shifting number) and ði; jÞ are the pixel coordinates of the image collected by the CCD. I k ði; jÞ, aði; jÞ, bði; jÞ denote the intensity value, background intensity, and the fringe visibility. Equation (1) can be written into a linear system by using a trigonometric identity: I k ði; jÞ ¼ aði; jÞ þ bði; jÞ cosðΔϕk Þ cos½ϕði; jÞ − bði; jÞ sinðΔϕk Þ sin½ϕði; jÞ:
ð2Þ
The LSEs of the wavefront phases are performed on each individual spatial sampling point ði; jÞ. For a phase step larger than 3, the linear equation system is overdetermined and the mean square error (MSE) is defined on the difference between the sampled intensity values and estimated intensity values, as shown in Eq. (3). The least-squares principle is then applied, as shown in Eq. (4): MSE ¼ Σ½I e ði; jÞ − I k ði; jÞ2 ; i;j
∂I k ði; kÞ I k ði; jÞ ¼ 0; i;j ∂cn
−2Σ
n ¼ 1; 2; 3; …;
ð3Þ ð4Þ
where 2
3
aði; jÞ cði; jÞ ¼ 4 bði; jÞ cos ϕði; jÞ 5; bði; jÞ sin ϕði; jÞ
ð5Þ
where I e ði; jÞ is the estimated intensity. The minimization is carried out on the parameters cði; jÞ in Eq. (5). The LSE formulation here follows Greivenkamp’s choice of parameters. As shown in Eq. (5), this kind of LSE is effective for combining the temporal information to isolate the phase information ϕði; jÞ from the background noises aði; jÞ and bði; jÞ. Therefore, it is used in the proposed algorithm for the interferogram reconstruction (see Fig. 1). However, as shown in 588
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
Eq. (4), we assume that the MSE is only affected by the background variation aði; jÞ, the fringe visibility bði; jÞ, and the phase term ϕði; jÞ. Under this formulation, the phase-shift step Δϕk is assumed to be accurate and should not contribute to the error. In real experimental conditions, the piezoelectric actuator calibration error [15–17] and the fast drifting interference patterns under dynamic loading conditions can make the static phase-shift step assumption invalid. Because the phase-shift step is assumed accurate in traditional PSI, the inference of the phase-shift step is critical to the accuracy of the measurement, which depends on the way of introducing the phase shifting itself. Phase shifting in the interferometry can be introduced by the translation and the tilting of the optical components (i.e., mirror, wedge, and grating) [18]. For example, in moiré interferometry with phase shifting [19], a nanoscale resolution piezoelectric motion actuator and controller are used here to motorize the movement of the beam splitting grating to introduce phase shifting. Under this configuration, the estimation of the phase steps can be calculated from the relation between the shifted phase and the translation step of the beam splitting grating [19–21], and it is given as δ Δϕ ¼ 4π ; ð6Þ Λ where Δϕ is the shifted phase; Λ is the beam splitting period in nanometer; δ is the translation step of the beam splitter grating. For the estimation of the phase step to be unbiased, the phase shifting should only be introduced by the piezoelectric device and there should be no errors from the miscalibration of the actuator and its feedback controller, which is difficult (and sometimes impossible) to satisfy in a real-world testing situation, in which environmental vibration and thermal noise can be large enough to introduce a large bias into Eq. (6). The proposed algorithm is invented to solve this problem effectively by considering the shifted phase as a variable, which can be exacted from the measurement datasets. 3. Wavefront Phase Estimation Using the Continuous Wavelet Transform
In the proposed algorithm, the CWT is integrated into the algorithm for wavefront phase estimation (see Fig. 1). The CWT is implemented in a scanning manner (i.e., 1D CWT) to demonstrate the working principle of the proposed phase-shifting estimation framework. Nonetheless, 2D CWT can be integrated in the framework to replace the 1D CWT used here. Watkins et al. have shown that there are several properties with CWT suitable for extracting phases from the sinusoidal modulated signal, and the literature has reported several successful applications using it for processing the interference fringe patterns [20,22–26]. Compared with the Fouriertransform-based phase extraction methods, the wavelet-based methods give higher resolution when extracting localized phase information [20]. The
definition of CWT is to convolute the signal with a scaling wavelet family: Z W f ða; bÞ ¼
þ∞
−∞
f ðtÞψ a;b ðtÞdt;
ð7Þ
where f ðtÞ is the input signal, ψ a;b ðtÞ is the wavelet function, and a and b are the wavelet scaling and translational parameters, respectively. W f ða; bÞ is the wavelet transformation as a function of a; b. f is the complex conjugate of function f . The function ψ a;b is scaled and translated from a mother function ψ, and the relationship between them is 1 x−b ψ a;b ðxÞ ¼ pffiffiffi ψ : a a
ð8Þ
The complex Morlet wavelet is chosen here because of its sinusoidal feature weighted by the Gaussian distribution function. The mother wavelet function can be written as 2 1 t ; ψðtÞ ¼ pffiffiffiffiffiffi expð−jω0 tÞ exp − 2 2π
ð9Þ
where ω0 ¼ 5:336 following [27] in order to get the normalized result. Although the complex Morlet wavelet is not admissible in the strict sense, it can be considered admissible in a loosely satisfying sense. Starting with results from the phase shifting [refer to the interference equation Eq. (1)], we set the input signal as f ðtÞ ¼ cosðϕðtÞÞ. It has been shown that the point of maximum magnitude (called the “ridge”) of the CWT W f ða; bÞ has the following properties, max
that is, the phase at position b can be computed by choosing the largest response of wavelet banks [25,28]: arg W f ðamax ; bÞ ¼ ϕðbÞ; ð10Þ where amax is the scaling at W f ða; bÞ . Equamax tion (10) shows that the ridge of the maximum power of the CWT can retrieve the phase information. Phase estimation from the CWT suffers the errors from the inappropriate choices of the scaling factor a, the finite length of the signal input, and the spatial discontinuities of the input signal. In practice, the implementation of different ridge detection schemata [29,30] can effectively reduce the errors introduced by inappropriate choices of the scaling factor a. It is assumed that the errors happening here follow the statistical distribution. The statistical analysis in the next section validates this.
the deterministic nature of the phase step. This is described in the interference intensity equation [Eq. (1)]. This is also the basic assumption for PSI [11,12]. The phase step is assumed to be constant over the space at the specific sampling moment, but the measured phase step is prone to include errors. Hence, the phase step can be treated as a statistical variable, which is expressed in Eq. (11): Δϕk ði; jÞ ¼ δk þ εk ði; jÞ;
ð11Þ
where Δϕk ði; jÞ is the phase difference between two interferograms calculated by the CWT, δk is the true value of the phase-shifting step, and εk ði; jÞ is the error for each estimation calculated from the phase differences using CWT. After aggregating the phase difference results, the probability distribution of the measured phase-shift step can be established. The error term in Eq. (11) is assumed to follow certain statistical distribution. As shown in the flow chart (see Fig. 1), both parametric estimation and nonparametric estimation can be used to estimate δk . In parametric estimation, a parameterized population distribution is required to fit the sample probability distribution. By contrast, nonparametric estimation makes no assumptions on the population distribution. In this paper, the parametric estimation, combined with the MLE, is chosen as the statistical estimation method to demonstrate the feasibility and effectiveness of the proposed algorithm. However, any other statistical estimation can be used in the proposed algorithm, as shown in Fig. 1. The MLE is generally consistent, which means that a larger sample size will generally help increase the accuracy of the estimation [31]. The principle of MLE is to find the estimator vector θ to maximize the likelihood function Lðθ; ωÞ, which is related to the probability function pðω; θÞ around the observation vector ω [32–34]. Considering that the observations are independent to one another (i.e., pðω; θÞ ¼ Πpn ðωn ; θÞ, where n is the number of obn servations), the MLE is as shown in Eq. (12) (the log-likelihood function is used): X ∇θ lnðpn ðωn ; θÞÞ ¼ 0:
ð12Þ
n
The variance of the MLE estimator can be judged by the asymptotic covariance matrix, which can be derived from the Fisher information. The Fisher information is the variance of the Fisher score sθ ¼ ∇θ lnðpn ðωn ; θÞÞ. Under regular conditions, the Fisher information matrix can be simplified as shown in Eq. (13) and the asymptotic covariance matrix is equal to the inverse of the Fisher information matrix [31]:
4. Statistical Analysis of the Phase Step Distribution
The phase information itself offers no statistical properties, but the phase difference at the same sampling location over the entire interferogram should follow the same statistical distribution because of
FðθÞ ¼ Eð½∇θ lnðpðω; θÞÞ½∇θ lnðpðω; θÞÞT jθÞ;
ð13Þ
where Eð…jθÞ denotes the expectation over a variable with respect to the probability function pðω; θÞ. 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
589
Table 1.
Normal distribution Laplace distribution Von Mises distribution a
Probability Distribution Candidates for the Parametric Parameter Estimation
Probability Density Function ðϕ−μÞ2 1 ffi f ðϕÞ ¼ pffiffiffiffiffiffiffi exp − 2 2σ 2πσ 2 jϕ−μj 1 f ðϕÞ ¼ 2b exp − b f ðϕÞ ¼ 2πI10 ðκÞ exp½κðϕ − μÞ
MLE Parameter 1 P μ ^ ¼ ni¼1 ϕi
MLE Parameter 2 P σ^2 ¼ n1 ni¼1 ðϕi − μ ^Þ2
μ ^ ¼ medianðϕi Þ Pn sin ϕi μ ^ ¼ arctan Pni¼1
^ ¼ 1 Pn jϕi − μ b ^j i¼1 n Pn cosðϕi −^ μÞ a I1 ðκÞ i¼1 n I0 ðκÞ ¼
i¼1
cos αϕ
Where I α is the modified Bessel function of the first kind.
For parametric parameter estimation, it is necessary to make an assumption on the population distribution. There are several distributions that are good candidates for the phase step, and they are listed in Table 1. Although all three candidates belong to the family of exponential distributions, the von Mises distribution (or circular normal distribution) considers the data as circular data, which is more closely related to the data structure of the phase information [35,36]. 5. Validation Using Simulated and Experimental Data
In this section, the proposed algorithm is validated using both simulated and experimental data. The proposed algorithm is implemented in MATLAB code. Unlike the experimental data, whose true phase-shift step is unknown, the advantage of using the numerically simulated data is that they provide true values to be compared with. In the validation process of the numerical data, the simulated input signal was a single- frequency sinusoidal signal with a signal length of L ¼ 1008 (which was the same length as the detector of the lab setup in the experimental data below). The data are synthetically signal generated with Eq. (14): s½n ¼
1 n 1 þ cos 2π · 10 · þ δk ; 2 L
ð14Þ
where s½n is the synthetic signal; n is the spatial coordinate; δk ¼ 0, 2=5π, 4=5π, 6=5π are the phase-
Fig. 2. (Color online) Input signal at δ1 ¼ 0 and δ3 ¼ 4=5π. 590
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
shifting steps; and L is the length of the signal. The plot of the signals is as shown in Fig. 2. The CWT was performed on the input signals here at δ1 ¼ 0 and δ3 ¼ 4=5π. The differences of the phase were taken both in angular scalar values and in vector values. For the normal distribution and Laplace distribution, the scalar angular values were used. For the von Mises distribution, which is based on circular statistics, the vector values were used. The wavefront phases from CWT were plotted in Fig. 3. The distribution fittings, using a different distribution on the phase-shifting step, were plotted in Fig. 4. The parameters were estimated by the MLE using Eq. (12), and the asymptotic covariance matrix was calculated using Eq. (13). The estimated parameters and variances are listed in Table 2. The results from the numerical testing showed that the choice of the population distribution was important in order to obtain a good result. The normal distribution assumption on the phase difference distribution was not a fair assumption on the simulated signal. The von Mises distribution gave a better estimation. The Laplace distribution gave the exact phase-shift steps with very small estimation variance (refer to Table 2). The reason for this is that CWT operation introduces only a very small amount of the error during operation. The error from the CWT was explained earlier. In contrast, the real data in the next section showed that the error term εk ði; jÞ in Eq. (11) followed different error patterns because
Fig. 3. (Color online) CWT estimation of the phase difference for δ3 ¼ 4=5π compared with the true value.
included 101,808 individual sampling points (see Fig. 5). From our perspective, the size of the samples is large enough for a statistically significant inference. The algorithm was performed on the experimental data in the same manner as it was for the simulated data. The probability distributions of the phase-shift step are plotted in Fig. 6, and the estimated parameters and variances calculated from the Fisher information are listed in Table 3. 6. Discussion
Fig. 4. (Color online) Probability distribution of the CWT estimation of the phase-shift step in the numerical simulation for δ3 ¼ 4=5π.
of the optical defects and the variation of the spatial frequency in the interference pattern, which was not presented in the numerical simulation shown above. The experimental data used here were taken from a real moiré interferometry experiment. For details of the moiré interferometry and the experiment configuration used to obtain the data, please refer to our previous paper [19]. On the one hand, the phaseshifting step was introduced by the translation of the piezoelectric device by 100 nm at each temporal sampling, which could be inferred from Eq. (6); on the other hand, the phase drift was also introduced by the loading condition of the experiment (the thermal expansion of the experiment stage, and air flow), which was undeterminable before the actual experiments. The specimen was a copper sheet (50 mm × 25 mm with 1 mm thickness) under current stressing (26 A) in the longitudinal direction. Because this paper addresses the phase-shifting step estimation, the details of the experiment will not be included. The phase step estimated from Eq. (6) was 1:5080 rad. The CCD sensor was a Pulnix TM1040 which has one million valid pixel points. The area of interest was the image row 300 to 400, which Table 2.
Distribution Normal Laplace Von Mises
The results show that the phase-shift step is a statistical variable instead of a deterministic one. In the numerical case, the proposed algorithm gave the exact value of the phase-shift value. In the experimental case, all three estimations from different underlying population distribution assumptions gave different numbers, and they deviated from the original estimation, which was inferred from a control feedback circuit of the piezoelectric actuator using Eq. (6). The variance was due to the error term in Eq. (11), which can be introduced either by any optical defects within the instruments or by a digital sampling process. Based upon the comparison between three preferred distribution candidates, both the Laplace and von Mises distributions have significant values in estimation and accuracy. The validation only implemented MLE and Fisher information covariance calculation. However, most of the parametric and nonparametric parameter estimation methods can be used to replace the statistical estimation block. The decoupling of different components is one of the major advantages of the proposed algorithm. The MLE calculated the results based on the likelihood function built on the sample frequency. A larger sample size, usually from a digital detector, decreased the variance of the estimated parameter.
Parameter Estimate for the Numerical Simulation
Parameter
Estimate
Standard Error
Real
μ ^ σ^ μ ^ ^ b
2.524 0.340 2.5133 0.1523 2.5204 9.7344
0.011 0.008 0.0002 0.0048 0.0104 0.4208
2.5133 N/A 2.5133 N/A 2.5133 N/A
μ ^ ^κ
Fig. 5. Moiré interferometry data at (a) δ1 and (b) δ3 . 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
591
7. Conclusions
Fig. 6. (Color online) Probability of the CWT estimation of the phase-shifting step in the moiré interferometry experimental data for δ3.
In contrast, the deterministic method cannot explore this feature from the PSI data source. The significance of the study includes building the phase step estimation on the statistical ground and, also, being able to establish the phase-shift step variable directly from the interferograms. From the same perspective, the existing algorithm can also be interpreted as building a statistical model on the intensity value. For example, the LSE on the MSEs of the intensity values, using phase steps as the minimization parameters [8,10], can be treated equally as an LSE on the phase steps, considering that the noise is included in the intensity term and the noise follows a normal distribution [34]. Because the intensity values were determined by multiple variables [refer to Eq. (1)], it will be difficult to establish a large number of samples with the same statistical structures of the errors. The estimation from the proposed algorithm gives a simple structure of the estimated parameter [refer to Eq. (11)]. This allows the researchers and engineers to design or utilize a more specific and precise statistical estimation procedure. Table 3. Parameter Estimate for the Moiré Interferometry Experimental Data
Distribution Parameter μ ^ σ^ μ ^ ^ b
Normal Laplace Von Mises
592
μ ^ ^κ
Estimate
Standard Error
Original Estimation
3.1432 0.6869 3.1001 0.4573 3.1032 3.0955
0.0022 0.0015 0.0016 0.0014 0.0020 0.0121
3.016 N/A 3.016 N/A 3.016 N/A
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
In this paper, a statistical phase-shifting step estimation algorithm, based on the CWT for high-resolution PSI was proposed, and a statistical analysis on the phase-shifting step in temporal PSI, based on CWT, was included. The proposed algorithm is an open framework, treating the phase-shifting step as a regular statistical variable. Inside this framework, the statistical estimation component can be easily replaced by other advanced parametric or nonparametric methods. Every point in the interferogram can be included in the samples to build up the statistical population. In the algorithm, the wavefront phases are exacted based on the maximum power ridge property of the CWT. The error could by introduced largely by the disruption in the data in the intensity map caused by optical defects on the light path. This was validated by utilizing the algorithm both in the numerical simulated data and experimental data. In the numerically simulated data, the exact phase-shifting step was obtained up to four significant figures. In the experimental study, the experimental data from the phase-shifting moiré interferometry were put though the algorithm. This gave a larger variance compared with the numerical one. However, a large sampling size (around 100,000 samples) for each interferogram enables a statistically significant inference based on MLE and the Fisher information based on asymptotic variance calculation. The validation showed the proposed method was enabled to estimate the phase steps accurately. The proposed methods can not only be used as a stand-alone method to estimate the phase, but also to provide a new statistical tool for temporal PSI. This project has been sponsored by the United States Navy Office of Naval Research (ONR) Advanced Electrical Power System under the direction of Terry Ericsen. References 1. 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–2703 (1974). 2. C. J. Morgan, “Least-squares estimation in phasemeasurement interferometry,” Opt. Lett. 7, 368–370 (1982). 3. J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. 23, 350–352 (1984). 4. G. Lai and T. Yatagai, “Generalized phase-shifting interferometry,” J. Opt. Soc. Am. A 8, 822–827 (1991). 5. K. Larkin, “A self-calibrating phase-shifting algorithm based on the natural demodulation of two-dimensional fringe patterns,” Opt. Express 9, 236–253 (2001). 6. C. T. Farrell and M. A. Player, “Phase-step insensitive algorithms for phase-shifting interferometry,” Meas. Sci. Technol. 5, 648–654 (1994). 7. C. Wei, M. Chen, and Z. Wang, “General phase-stepping algorithm with automatic calibration of phase steps,” Opt. Eng. 38, 1357–1360 (1999). 8. I.-B. Kong and S.-W. Kim, “General algorithm of phaseshifting interferometry by iterative least-squares fitting,” Opt. Eng. 34, 183–188 (1995).
9. M. Rivera, R. Bizuet, A. Martinez, and J. A. Rayas, “Halfquadratic cost function for computing arbitrary phase shifts and phase: adaptive out of step phase shifting,” Opt. Express 14, 3204–3213 (2006). 10. Z. Wang and B. Han, “Advanced iterative algorithm for randomly phase-shifted interferograms with intra- and interframe intensity variations,” Opt. Lasers Eng. 45, 274–280 (2007). 11. L. Z. Cai, Q. Liu, and X. L. Yang, “Phase-shift extraction and wave-front reconstruction in phase-shifting interferometry with arbitrary phase steps,” Opt. Lett. 28, 1808–1810 (2003). 12. L. Z. Cai, Q. Liu, and X. L. Yang, “Generalized phase-shifting interferometry with arbitrary unknown phase steps for diffraction objects,” Opt. Lett. 29, 183–185 (2004). 13. X. F. Xu, L. Z. Cai, X. F. Meng, G. Y. Dong, and X. X. Shen, “Fast blind extraction of arbitrary unknown phase shifts by an iterative tangent approach in generalized phase-shifting interferometry,” Opt. Lett. 31, 1966–1968 (2006). 14. H. Guo, Z. Zhao, and M. Chen, “Efficient iterative algorithm for phase-shifting interferometry,” Opt. Lasers Eng. 45, 281–292 (2007). 15. Y.-Y. Cheng and J. C. Wyant, “Phase shifter calibration in phaseshifting interferometry,” Appl. Opt. 24, 3049–3052 (1985). 16. C. Ai and J. C. Wyant, “Effect of piezoelectric transducer nonlinearity on phase shift interferometry,” Appl. Opt. 26, 1112–1116 (1987). 17. B. Gutmann and H. Weber, “Phase-shifter calibration and error detection in phase-shifting applications: a new method,” Appl. Opt. 37, 7624–7631 (1998). 18. D. Malacara, S. Manuel, and M. Zacarias, Interferogram Analysis for Optical Testing, 2nd ed., Optical Engineering and Science (Taylor & Francis, 2005). 19. B. Chen and C. Basaran, “Automatic full strain field Moire interferometry measurement with nano-scale resolution,” Exp. Mech. 48, 665–673 (2008). 20. H. Liu, “Phase reconstruction of phase shifted moire interferograms using continuous wavelet transforms,” Ph.D. dissertation (State University of New York at Buffalo, 2003). 21. C. Han and B. Han, “Phase-shifting in achromatic moiré interferometry system,” Opt. Express 15, 9970–9976 (2007).
22. L. R. Watkins, S. M. Tan, and T. H. Barnes, “Interferometer profile extraction using continuous wavelet transform,” Electron. Lett. 33, 2116–2117 (1997). 23. L. R. Watkins, S. M. Tan, and T. H. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24, 905–907 (1999). 24. H. Liu, A. N. Cartwright, and C. Basaran, “Moire interferogram phase extraction: a ridge detection algorithm for continuous wavelet transforms,” Appl. Opt. 43, 850–857 (2004). 25. L. R. Watkins, “Phase recovery from fringe patterns using the continuous wavelet transform,” Opt. Lasers Eng. 45, 298–303 (2007). 26. C. A. Sciammarella and T. Kim, “Frequency modulation interpretation of fringes and computation of strains,” Exp. Mech. 45, 393–403 (2005). 27. M. Vetterli and J. Kovacevic, Wavelets and Subband Coding (Prentice-Hall, 1995). 28. M. Afifi, A. Fassi-Fihri, M. Marjane, K. Nassim, M. Sidki, and S. Rachafi, “Paul wavelet-based algorithm for optical phase distribution evaluation,” Opt. Commun. 211, 47–51 (2002). 29. N. Delprat, B. Escudie, P. Guillemain, R. Kronland-Martinet, P. Tchamitchian, and B. Torresani, “Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies,” IEEE Trans. Inf. Theory 38, 644–664 (1992). 30. R. A. Carmona, W. L. Hwang, and B. Torresani, “Characterization of signals by the ridges of their wavelet transforms,” IEEE Trans. Signal Process. 45, 2586–2590 (1997). 31. A. van den Bos, Parameter Estimation for Scientists and Engineers (Wiley-Interscience, 2007). 32. A. W. F. Edwards, Likelihood: An Account of the Statistical Concept of Likelihood and Its Application to Scientific Inference (Cambridge University, 1972). 33. L. R. Shenton and K. O. Bowman, Maximum Likelihood Estimation in Small Samples, Griffin’s Statistical Monographs & Courses (Macmillan, 1977). 34. S. M. Kay, Fundamentals of Statistical Signal Processing, Prentice Hall Signal Processing Series (Prentice-Hall, 1993). 35. E. Batschelet, Circular Statistics in Biology, Mathematics in Biology (Academic, 1981). 36. S. R. S. A. Jammalamadaka, Topics in Circular Statistics, Series on Multivariate Analysis (World Scientific, 2001).
1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
593
We propose a statistical phase-shifting estimation algorithm for temporal phase-shifting interferometry (PSI) based on the continuous wavelet transform (CWT). The proposed algorithm explores spatial information redundancy in the intraframe interferogram dataset using the phase recovery property on the power ridge of the CWT. Despite the errors introduced by the noise of the interferogram, the statistical part of the algorithm is utilized to give a sound estimation of the phase-shifting step. It also introduces the usage of directional statistics as the statistical model, which was validated, so as to offer a better estimation compared with other statistical models. The algorithm is implemented in computer codes, and the validations of the algorithm were performed on numerical simulated signals and actual phase-shifted moiré interferograms. The major advantage of the proposed algorithm is that it imposes weaker conditions on the presumptions in the temporal PSI, which, under most circumstances, requires uniform and precalibrated phase-shifting steps. Compared with other existing deterministic estimation algorithms, the proposed algorithm estimates the phase-shifting step statistically. The proposed algorithm allows the temporal PSI to operate under dynamic loading conditions and arbitrary phase steps and also without precalibration of the phase shifter. The proposed method can serve as a benchmark method for comparing the accuracy of the different phase-step estimation methods. © 2011 Optical Society of America OCIS codes: 050.5080, 000.5490, 100.5070, 100.7410, 120.2880, 120.3180.
1. Introduction
Phase-shifting interferometry (PSI) is an important technique in high-precision optical metrology. Temporal PSI was first introduced by Bruning et al. [1] to overcome the λ=2 sensitivity limits of the fringe counting methods in the interferometry for the application in the surface quality measurement of the optical components. In the Bruning et al. paper, error sources are enumerated, including electrical noise in the detectors, mechanical drifts, atmospheric turbulence and the noise in the laser source. It is also emphasized that the error can be significant when the sensitivity reaches to λ=100. The arctangent function is used in Bruning’s paper to reconstruct the 0003-6935/11/040586-08$15.00/0 © 2011 Optical Society of America 586
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
phase information. The least-squares estimation (LSE) method was first introduced by Morgan [2] to minimize the intensity errors nonlinearly between the imaged intensity map and the estimated one on the parameter of the measured phase. The commonly used version of LSE was established by Greivenkamp [3]. Greivenkamp uses the same interference equation but minimizes the intensity errors according to three parameters (background, cosine, and sine function of the phase shifter). Greivenkamp’s LSE is linear and removes several restrictions when designing the phase steps. Lai and Yatagai apply additional interference setups to detect the shifted phase, while using Fourier transform methods [4]. Larkin added preprocessing and postprocessing to improve the robustness of Lai and Yatagai’s Fourier-transform-based method [5]. Farrell and Player once again pointed out that
the error in the phase steps can cause a distortion of the results, and they suggested an interpixel ellipse fitting estimation method [6]. The phase step can be unknown and unequal. Wei et al. gave a leastsquares algorithm based on the interpixel relationship by means of a simple intensity transform between the sampled interferograms [7]. Iterative LSE was proposed by Kong and Kim, and considers both the phase and phase shift as parameters [8]. The iterative method and its application was also documented by Rivera et al. [9]. Based on the initially induced phase step (e.g., from the linearity table of piezoelectric actuator), Greivenkamp’s LSE is first used to estimate the wavefront phases, and a modified Greivenkamp linear LSE is used to estimate the phase step afterward. The iterative process is therefore built on these two LSE models to obtain the final wavefront phases upon the convergence condition. The idea is generalized and extended by Wang and Han [10]. Another phase-shift extraction method, based on the assumed object wave, is proposed by Cai et al. [11,12]. The Cai et al. method uses the mean of the differences between two interferograms to calculate the phase difference. An iterative process is available using the estimated phase steps to reconstruct the phase distribution iteratively. Xu et al. proposed a noniterative version based on the same principle [13]. The methods reviewed above can be divided into two categories: temporal methods and temporalspatial methods. The temporal methods rely on an established phase step to give pixelwise estimations on the phase distribution. The temporal-spatial methods use the spatial data redundancy to first estimate the phase step and then to reconstruct the phase distribution from the temporal samples [5]. The proposed algorithm in this paper belongs to the temporal-spatial method. The flow chart of the proposed algorithm is shown in Fig. 1. Compared with the existing algorithms, the proposed algorithm utilizes the continuous wavelet transform (CWT) to estimate the wavefront phases first, which is operated using spatial data correlation within an interferogram frame. Later, the phase steps are estimated statistically upon the phase difference distributions between temporal sampled interferogram frames. The parameter estimation part of the algorithm was demonstrated here using the maximum likelihood estimation (MLE). However, any statistical parameter estimation method or nonparametric estimation can be used instead. The proposed method allows an arbitrary phase step but does not involve any iterative process, which differs from most of the existing temporal-spatial methods. Furthermore, the proposed algorithm is a framework for building a robust phase-shifting algorithm, which can allow the PSI to operate under a dynamic loading condition. According to our literature survey, this is the first time that the CWT is reported to be used to help in the estimation of the phase-shifting steps.
Fig. 1. Flow chart of the proposed phase-shifting step estimation algorithm.
2. Interferometry Phase Least-Squares Estimation
The structure and notation of the interference equation are introduced in this section. They are critical elements to get an understanding of the proposed algorithm. This section also lays out some essential 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
587
concepts for the proposed algorithm, especially on the formulation and parameterization of the LSE used in the field of phase-shift estimations. In common circumstances, the temporal PSI samples multiple intensity values on the same spatial point at different times and with shifted phase steps. For the detector using a charge coupled device (CCD), each sampling is an accumulation process of the electrical charges that are proportional to the intensity of the incident photon flux integrated over the exposure time. The reconstruction of the temporal phaseshifted interferograms using the least-squares formula is widely used [3,14]. This serves in the last step in the proposed algorithm. The interference equation for the temporal PSI can be written as [14] I k ði; jÞ ¼ aði; jÞ þ bði; jÞ cos½ϕði; jÞ þ Δϕk ;
ð1Þ
where k is the number of phase shifting k ¼ 0, 1, …, K − 1 (K is the total shifting number) and ði; jÞ are the pixel coordinates of the image collected by the CCD. I k ði; jÞ, aði; jÞ, bði; jÞ denote the intensity value, background intensity, and the fringe visibility. Equation (1) can be written into a linear system by using a trigonometric identity: I k ði; jÞ ¼ aði; jÞ þ bði; jÞ cosðΔϕk Þ cos½ϕði; jÞ − bði; jÞ sinðΔϕk Þ sin½ϕði; jÞ:
ð2Þ
The LSEs of the wavefront phases are performed on each individual spatial sampling point ði; jÞ. For a phase step larger than 3, the linear equation system is overdetermined and the mean square error (MSE) is defined on the difference between the sampled intensity values and estimated intensity values, as shown in Eq. (3). The least-squares principle is then applied, as shown in Eq. (4): MSE ¼ Σ½I e ði; jÞ − I k ði; jÞ2 ; i;j
∂I k ði; kÞ I k ði; jÞ ¼ 0; i;j ∂cn
−2Σ
n ¼ 1; 2; 3; …;
ð3Þ ð4Þ
where 2
3
aði; jÞ cði; jÞ ¼ 4 bði; jÞ cos ϕði; jÞ 5; bði; jÞ sin ϕði; jÞ
ð5Þ
where I e ði; jÞ is the estimated intensity. The minimization is carried out on the parameters cði; jÞ in Eq. (5). The LSE formulation here follows Greivenkamp’s choice of parameters. As shown in Eq. (5), this kind of LSE is effective for combining the temporal information to isolate the phase information ϕði; jÞ from the background noises aði; jÞ and bði; jÞ. Therefore, it is used in the proposed algorithm for the interferogram reconstruction (see Fig. 1). However, as shown in 588
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
Eq. (4), we assume that the MSE is only affected by the background variation aði; jÞ, the fringe visibility bði; jÞ, and the phase term ϕði; jÞ. Under this formulation, the phase-shift step Δϕk is assumed to be accurate and should not contribute to the error. In real experimental conditions, the piezoelectric actuator calibration error [15–17] and the fast drifting interference patterns under dynamic loading conditions can make the static phase-shift step assumption invalid. Because the phase-shift step is assumed accurate in traditional PSI, the inference of the phase-shift step is critical to the accuracy of the measurement, which depends on the way of introducing the phase shifting itself. Phase shifting in the interferometry can be introduced by the translation and the tilting of the optical components (i.e., mirror, wedge, and grating) [18]. For example, in moiré interferometry with phase shifting [19], a nanoscale resolution piezoelectric motion actuator and controller are used here to motorize the movement of the beam splitting grating to introduce phase shifting. Under this configuration, the estimation of the phase steps can be calculated from the relation between the shifted phase and the translation step of the beam splitting grating [19–21], and it is given as δ Δϕ ¼ 4π ; ð6Þ Λ where Δϕ is the shifted phase; Λ is the beam splitting period in nanometer; δ is the translation step of the beam splitter grating. For the estimation of the phase step to be unbiased, the phase shifting should only be introduced by the piezoelectric device and there should be no errors from the miscalibration of the actuator and its feedback controller, which is difficult (and sometimes impossible) to satisfy in a real-world testing situation, in which environmental vibration and thermal noise can be large enough to introduce a large bias into Eq. (6). The proposed algorithm is invented to solve this problem effectively by considering the shifted phase as a variable, which can be exacted from the measurement datasets. 3. Wavefront Phase Estimation Using the Continuous Wavelet Transform
In the proposed algorithm, the CWT is integrated into the algorithm for wavefront phase estimation (see Fig. 1). The CWT is implemented in a scanning manner (i.e., 1D CWT) to demonstrate the working principle of the proposed phase-shifting estimation framework. Nonetheless, 2D CWT can be integrated in the framework to replace the 1D CWT used here. Watkins et al. have shown that there are several properties with CWT suitable for extracting phases from the sinusoidal modulated signal, and the literature has reported several successful applications using it for processing the interference fringe patterns [20,22–26]. Compared with the Fouriertransform-based phase extraction methods, the wavelet-based methods give higher resolution when extracting localized phase information [20]. The
definition of CWT is to convolute the signal with a scaling wavelet family: Z W f ða; bÞ ¼
þ∞
−∞
f ðtÞψ a;b ðtÞdt;
ð7Þ
where f ðtÞ is the input signal, ψ a;b ðtÞ is the wavelet function, and a and b are the wavelet scaling and translational parameters, respectively. W f ða; bÞ is the wavelet transformation as a function of a; b. f is the complex conjugate of function f . The function ψ a;b is scaled and translated from a mother function ψ, and the relationship between them is 1 x−b ψ a;b ðxÞ ¼ pffiffiffi ψ : a a
ð8Þ
The complex Morlet wavelet is chosen here because of its sinusoidal feature weighted by the Gaussian distribution function. The mother wavelet function can be written as 2 1 t ; ψðtÞ ¼ pffiffiffiffiffiffi expð−jω0 tÞ exp − 2 2π
ð9Þ
where ω0 ¼ 5:336 following [27] in order to get the normalized result. Although the complex Morlet wavelet is not admissible in the strict sense, it can be considered admissible in a loosely satisfying sense. Starting with results from the phase shifting [refer to the interference equation Eq. (1)], we set the input signal as f ðtÞ ¼ cosðϕðtÞÞ. It has been shown that the point of maximum magnitude (called the “ridge”) of the CWT W f ða; bÞ has the following properties, max
that is, the phase at position b can be computed by choosing the largest response of wavelet banks [25,28]: arg W f ðamax ; bÞ ¼ ϕðbÞ; ð10Þ where amax is the scaling at W f ða; bÞ . Equamax tion (10) shows that the ridge of the maximum power of the CWT can retrieve the phase information. Phase estimation from the CWT suffers the errors from the inappropriate choices of the scaling factor a, the finite length of the signal input, and the spatial discontinuities of the input signal. In practice, the implementation of different ridge detection schemata [29,30] can effectively reduce the errors introduced by inappropriate choices of the scaling factor a. It is assumed that the errors happening here follow the statistical distribution. The statistical analysis in the next section validates this.
the deterministic nature of the phase step. This is described in the interference intensity equation [Eq. (1)]. This is also the basic assumption for PSI [11,12]. The phase step is assumed to be constant over the space at the specific sampling moment, but the measured phase step is prone to include errors. Hence, the phase step can be treated as a statistical variable, which is expressed in Eq. (11): Δϕk ði; jÞ ¼ δk þ εk ði; jÞ;
ð11Þ
where Δϕk ði; jÞ is the phase difference between two interferograms calculated by the CWT, δk is the true value of the phase-shifting step, and εk ði; jÞ is the error for each estimation calculated from the phase differences using CWT. After aggregating the phase difference results, the probability distribution of the measured phase-shift step can be established. The error term in Eq. (11) is assumed to follow certain statistical distribution. As shown in the flow chart (see Fig. 1), both parametric estimation and nonparametric estimation can be used to estimate δk . In parametric estimation, a parameterized population distribution is required to fit the sample probability distribution. By contrast, nonparametric estimation makes no assumptions on the population distribution. In this paper, the parametric estimation, combined with the MLE, is chosen as the statistical estimation method to demonstrate the feasibility and effectiveness of the proposed algorithm. However, any other statistical estimation can be used in the proposed algorithm, as shown in Fig. 1. The MLE is generally consistent, which means that a larger sample size will generally help increase the accuracy of the estimation [31]. The principle of MLE is to find the estimator vector θ to maximize the likelihood function Lðθ; ωÞ, which is related to the probability function pðω; θÞ around the observation vector ω [32–34]. Considering that the observations are independent to one another (i.e., pðω; θÞ ¼ Πpn ðωn ; θÞ, where n is the number of obn servations), the MLE is as shown in Eq. (12) (the log-likelihood function is used): X ∇θ lnðpn ðωn ; θÞÞ ¼ 0:
ð12Þ
n
The variance of the MLE estimator can be judged by the asymptotic covariance matrix, which can be derived from the Fisher information. The Fisher information is the variance of the Fisher score sθ ¼ ∇θ lnðpn ðωn ; θÞÞ. Under regular conditions, the Fisher information matrix can be simplified as shown in Eq. (13) and the asymptotic covariance matrix is equal to the inverse of the Fisher information matrix [31]:
4. Statistical Analysis of the Phase Step Distribution
The phase information itself offers no statistical properties, but the phase difference at the same sampling location over the entire interferogram should follow the same statistical distribution because of
FðθÞ ¼ Eð½∇θ lnðpðω; θÞÞ½∇θ lnðpðω; θÞÞT jθÞ;
ð13Þ
where Eð…jθÞ denotes the expectation over a variable with respect to the probability function pðω; θÞ. 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
589
Table 1.
Normal distribution Laplace distribution Von Mises distribution a
Probability Distribution Candidates for the Parametric Parameter Estimation
Probability Density Function ðϕ−μÞ2 1 ffi f ðϕÞ ¼ pffiffiffiffiffiffiffi exp − 2 2σ 2πσ 2 jϕ−μj 1 f ðϕÞ ¼ 2b exp − b f ðϕÞ ¼ 2πI10 ðκÞ exp½κðϕ − μÞ
MLE Parameter 1 P μ ^ ¼ ni¼1 ϕi
MLE Parameter 2 P σ^2 ¼ n1 ni¼1 ðϕi − μ ^Þ2
μ ^ ¼ medianðϕi Þ Pn sin ϕi μ ^ ¼ arctan Pni¼1
^ ¼ 1 Pn jϕi − μ b ^j i¼1 n Pn cosðϕi −^ μÞ a I1 ðκÞ i¼1 n I0 ðκÞ ¼
i¼1
cos αϕ
Where I α is the modified Bessel function of the first kind.
For parametric parameter estimation, it is necessary to make an assumption on the population distribution. There are several distributions that are good candidates for the phase step, and they are listed in Table 1. Although all three candidates belong to the family of exponential distributions, the von Mises distribution (or circular normal distribution) considers the data as circular data, which is more closely related to the data structure of the phase information [35,36]. 5. Validation Using Simulated and Experimental Data
In this section, the proposed algorithm is validated using both simulated and experimental data. The proposed algorithm is implemented in MATLAB code. Unlike the experimental data, whose true phase-shift step is unknown, the advantage of using the numerically simulated data is that they provide true values to be compared with. In the validation process of the numerical data, the simulated input signal was a single- frequency sinusoidal signal with a signal length of L ¼ 1008 (which was the same length as the detector of the lab setup in the experimental data below). The data are synthetically signal generated with Eq. (14): s½n ¼
1 n 1 þ cos 2π · 10 · þ δk ; 2 L
ð14Þ
where s½n is the synthetic signal; n is the spatial coordinate; δk ¼ 0, 2=5π, 4=5π, 6=5π are the phase-
Fig. 2. (Color online) Input signal at δ1 ¼ 0 and δ3 ¼ 4=5π. 590
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
shifting steps; and L is the length of the signal. The plot of the signals is as shown in Fig. 2. The CWT was performed on the input signals here at δ1 ¼ 0 and δ3 ¼ 4=5π. The differences of the phase were taken both in angular scalar values and in vector values. For the normal distribution and Laplace distribution, the scalar angular values were used. For the von Mises distribution, which is based on circular statistics, the vector values were used. The wavefront phases from CWT were plotted in Fig. 3. The distribution fittings, using a different distribution on the phase-shifting step, were plotted in Fig. 4. The parameters were estimated by the MLE using Eq. (12), and the asymptotic covariance matrix was calculated using Eq. (13). The estimated parameters and variances are listed in Table 2. The results from the numerical testing showed that the choice of the population distribution was important in order to obtain a good result. The normal distribution assumption on the phase difference distribution was not a fair assumption on the simulated signal. The von Mises distribution gave a better estimation. The Laplace distribution gave the exact phase-shift steps with very small estimation variance (refer to Table 2). The reason for this is that CWT operation introduces only a very small amount of the error during operation. The error from the CWT was explained earlier. In contrast, the real data in the next section showed that the error term εk ði; jÞ in Eq. (11) followed different error patterns because
Fig. 3. (Color online) CWT estimation of the phase difference for δ3 ¼ 4=5π compared with the true value.
included 101,808 individual sampling points (see Fig. 5). From our perspective, the size of the samples is large enough for a statistically significant inference. The algorithm was performed on the experimental data in the same manner as it was for the simulated data. The probability distributions of the phase-shift step are plotted in Fig. 6, and the estimated parameters and variances calculated from the Fisher information are listed in Table 3. 6. Discussion
Fig. 4. (Color online) Probability distribution of the CWT estimation of the phase-shift step in the numerical simulation for δ3 ¼ 4=5π.
of the optical defects and the variation of the spatial frequency in the interference pattern, which was not presented in the numerical simulation shown above. The experimental data used here were taken from a real moiré interferometry experiment. For details of the moiré interferometry and the experiment configuration used to obtain the data, please refer to our previous paper [19]. On the one hand, the phaseshifting step was introduced by the translation of the piezoelectric device by 100 nm at each temporal sampling, which could be inferred from Eq. (6); on the other hand, the phase drift was also introduced by the loading condition of the experiment (the thermal expansion of the experiment stage, and air flow), which was undeterminable before the actual experiments. The specimen was a copper sheet (50 mm × 25 mm with 1 mm thickness) under current stressing (26 A) in the longitudinal direction. Because this paper addresses the phase-shifting step estimation, the details of the experiment will not be included. The phase step estimated from Eq. (6) was 1:5080 rad. The CCD sensor was a Pulnix TM1040 which has one million valid pixel points. The area of interest was the image row 300 to 400, which Table 2.
Distribution Normal Laplace Von Mises
The results show that the phase-shift step is a statistical variable instead of a deterministic one. In the numerical case, the proposed algorithm gave the exact value of the phase-shift value. In the experimental case, all three estimations from different underlying population distribution assumptions gave different numbers, and they deviated from the original estimation, which was inferred from a control feedback circuit of the piezoelectric actuator using Eq. (6). The variance was due to the error term in Eq. (11), which can be introduced either by any optical defects within the instruments or by a digital sampling process. Based upon the comparison between three preferred distribution candidates, both the Laplace and von Mises distributions have significant values in estimation and accuracy. The validation only implemented MLE and Fisher information covariance calculation. However, most of the parametric and nonparametric parameter estimation methods can be used to replace the statistical estimation block. The decoupling of different components is one of the major advantages of the proposed algorithm. The MLE calculated the results based on the likelihood function built on the sample frequency. A larger sample size, usually from a digital detector, decreased the variance of the estimated parameter.
Parameter Estimate for the Numerical Simulation
Parameter
Estimate
Standard Error
Real
μ ^ σ^ μ ^ ^ b
2.524 0.340 2.5133 0.1523 2.5204 9.7344
0.011 0.008 0.0002 0.0048 0.0104 0.4208
2.5133 N/A 2.5133 N/A 2.5133 N/A
μ ^ ^κ
Fig. 5. Moiré interferometry data at (a) δ1 and (b) δ3 . 1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
591
7. Conclusions
Fig. 6. (Color online) Probability of the CWT estimation of the phase-shifting step in the moiré interferometry experimental data for δ3.
In contrast, the deterministic method cannot explore this feature from the PSI data source. The significance of the study includes building the phase step estimation on the statistical ground and, also, being able to establish the phase-shift step variable directly from the interferograms. From the same perspective, the existing algorithm can also be interpreted as building a statistical model on the intensity value. For example, the LSE on the MSEs of the intensity values, using phase steps as the minimization parameters [8,10], can be treated equally as an LSE on the phase steps, considering that the noise is included in the intensity term and the noise follows a normal distribution [34]. Because the intensity values were determined by multiple variables [refer to Eq. (1)], it will be difficult to establish a large number of samples with the same statistical structures of the errors. The estimation from the proposed algorithm gives a simple structure of the estimated parameter [refer to Eq. (11)]. This allows the researchers and engineers to design or utilize a more specific and precise statistical estimation procedure. Table 3. Parameter Estimate for the Moiré Interferometry Experimental Data
Distribution Parameter μ ^ σ^ μ ^ ^ b
Normal Laplace Von Mises
592
μ ^ ^κ
Estimate
Standard Error
Original Estimation
3.1432 0.6869 3.1001 0.4573 3.1032 3.0955
0.0022 0.0015 0.0016 0.0014 0.0020 0.0121
3.016 N/A 3.016 N/A 3.016 N/A
APPLIED OPTICS / Vol. 50, No. 4 / 1 February 2011
In this paper, a statistical phase-shifting step estimation algorithm, based on the CWT for high-resolution PSI was proposed, and a statistical analysis on the phase-shifting step in temporal PSI, based on CWT, was included. The proposed algorithm is an open framework, treating the phase-shifting step as a regular statistical variable. Inside this framework, the statistical estimation component can be easily replaced by other advanced parametric or nonparametric methods. Every point in the interferogram can be included in the samples to build up the statistical population. In the algorithm, the wavefront phases are exacted based on the maximum power ridge property of the CWT. The error could by introduced largely by the disruption in the data in the intensity map caused by optical defects on the light path. This was validated by utilizing the algorithm both in the numerical simulated data and experimental data. In the numerically simulated data, the exact phase-shifting step was obtained up to four significant figures. In the experimental study, the experimental data from the phase-shifting moiré interferometry were put though the algorithm. This gave a larger variance compared with the numerical one. However, a large sampling size (around 100,000 samples) for each interferogram enables a statistically significant inference based on MLE and the Fisher information based on asymptotic variance calculation. The validation showed the proposed method was enabled to estimate the phase steps accurately. The proposed methods can not only be used as a stand-alone method to estimate the phase, but also to provide a new statistical tool for temporal PSI. This project has been sponsored by the United States Navy Office of Naval Research (ONR) Advanced Electrical Power System under the direction of Terry Ericsen. References 1. 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–2703 (1974). 2. C. J. Morgan, “Least-squares estimation in phasemeasurement interferometry,” Opt. Lett. 7, 368–370 (1982). 3. J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. 23, 350–352 (1984). 4. G. Lai and T. Yatagai, “Generalized phase-shifting interferometry,” J. Opt. Soc. Am. A 8, 822–827 (1991). 5. K. Larkin, “A self-calibrating phase-shifting algorithm based on the natural demodulation of two-dimensional fringe patterns,” Opt. Express 9, 236–253 (2001). 6. C. T. Farrell and M. A. Player, “Phase-step insensitive algorithms for phase-shifting interferometry,” Meas. Sci. Technol. 5, 648–654 (1994). 7. C. Wei, M. Chen, and Z. Wang, “General phase-stepping algorithm with automatic calibration of phase steps,” Opt. Eng. 38, 1357–1360 (1999). 8. I.-B. Kong and S.-W. Kim, “General algorithm of phaseshifting interferometry by iterative least-squares fitting,” Opt. Eng. 34, 183–188 (1995).
9. M. Rivera, R. Bizuet, A. Martinez, and J. A. Rayas, “Halfquadratic cost function for computing arbitrary phase shifts and phase: adaptive out of step phase shifting,” Opt. Express 14, 3204–3213 (2006). 10. Z. Wang and B. Han, “Advanced iterative algorithm for randomly phase-shifted interferograms with intra- and interframe intensity variations,” Opt. Lasers Eng. 45, 274–280 (2007). 11. L. Z. Cai, Q. Liu, and X. L. Yang, “Phase-shift extraction and wave-front reconstruction in phase-shifting interferometry with arbitrary phase steps,” Opt. Lett. 28, 1808–1810 (2003). 12. L. Z. Cai, Q. Liu, and X. L. Yang, “Generalized phase-shifting interferometry with arbitrary unknown phase steps for diffraction objects,” Opt. Lett. 29, 183–185 (2004). 13. X. F. Xu, L. Z. Cai, X. F. Meng, G. Y. Dong, and X. X. Shen, “Fast blind extraction of arbitrary unknown phase shifts by an iterative tangent approach in generalized phase-shifting interferometry,” Opt. Lett. 31, 1966–1968 (2006). 14. H. Guo, Z. Zhao, and M. Chen, “Efficient iterative algorithm for phase-shifting interferometry,” Opt. Lasers Eng. 45, 281–292 (2007). 15. Y.-Y. Cheng and J. C. Wyant, “Phase shifter calibration in phaseshifting interferometry,” Appl. Opt. 24, 3049–3052 (1985). 16. C. Ai and J. C. Wyant, “Effect of piezoelectric transducer nonlinearity on phase shift interferometry,” Appl. Opt. 26, 1112–1116 (1987). 17. B. Gutmann and H. Weber, “Phase-shifter calibration and error detection in phase-shifting applications: a new method,” Appl. Opt. 37, 7624–7631 (1998). 18. D. Malacara, S. Manuel, and M. Zacarias, Interferogram Analysis for Optical Testing, 2nd ed., Optical Engineering and Science (Taylor & Francis, 2005). 19. B. Chen and C. Basaran, “Automatic full strain field Moire interferometry measurement with nano-scale resolution,” Exp. Mech. 48, 665–673 (2008). 20. H. Liu, “Phase reconstruction of phase shifted moire interferograms using continuous wavelet transforms,” Ph.D. dissertation (State University of New York at Buffalo, 2003). 21. C. Han and B. Han, “Phase-shifting in achromatic moiré interferometry system,” Opt. Express 15, 9970–9976 (2007).
22. L. R. Watkins, S. M. Tan, and T. H. Barnes, “Interferometer profile extraction using continuous wavelet transform,” Electron. Lett. 33, 2116–2117 (1997). 23. L. R. Watkins, S. M. Tan, and T. H. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24, 905–907 (1999). 24. H. Liu, A. N. Cartwright, and C. Basaran, “Moire interferogram phase extraction: a ridge detection algorithm for continuous wavelet transforms,” Appl. Opt. 43, 850–857 (2004). 25. L. R. Watkins, “Phase recovery from fringe patterns using the continuous wavelet transform,” Opt. Lasers Eng. 45, 298–303 (2007). 26. C. A. Sciammarella and T. Kim, “Frequency modulation interpretation of fringes and computation of strains,” Exp. Mech. 45, 393–403 (2005). 27. M. Vetterli and J. Kovacevic, Wavelets and Subband Coding (Prentice-Hall, 1995). 28. M. Afifi, A. Fassi-Fihri, M. Marjane, K. Nassim, M. Sidki, and S. Rachafi, “Paul wavelet-based algorithm for optical phase distribution evaluation,” Opt. Commun. 211, 47–51 (2002). 29. N. Delprat, B. Escudie, P. Guillemain, R. Kronland-Martinet, P. Tchamitchian, and B. Torresani, “Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies,” IEEE Trans. Inf. Theory 38, 644–664 (1992). 30. R. A. Carmona, W. L. Hwang, and B. Torresani, “Characterization of signals by the ridges of their wavelet transforms,” IEEE Trans. Signal Process. 45, 2586–2590 (1997). 31. A. van den Bos, Parameter Estimation for Scientists and Engineers (Wiley-Interscience, 2007). 32. A. W. F. Edwards, Likelihood: An Account of the Statistical Concept of Likelihood and Its Application to Scientific Inference (Cambridge University, 1972). 33. L. R. Shenton and K. O. Bowman, Maximum Likelihood Estimation in Small Samples, Griffin’s Statistical Monographs & Courses (Macmillan, 1977). 34. S. M. Kay, Fundamentals of Statistical Signal Processing, Prentice Hall Signal Processing Series (Prentice-Hall, 1993). 35. E. Batschelet, Circular Statistics in Biology, Mathematics in Biology (Academic, 1981). 36. S. R. S. A. Jammalamadaka, Topics in Circular Statistics, Series on Multivariate Analysis (World Scientific, 2001).
1 February 2011 / Vol. 50, No. 4 / APPLIED OPTICS
593
