Regularized least squares phase sampling interferometry

In phase sampling interferometry, existing temporal analysis methods are sensitive to border effects and cannot deal with missing data. In this work we propose a quadrature filter that allows a reliable dynamic phase measurement for every sample, even in the cases involving few samples or missing data. The method is based on the use of a regularized least squares cost function that enforces the quadrature character of the filter. A comparison with existing techniques shows the effectiveness of the proposed method. ©2011 Optical Society of America OCIS codes: (120.3180) Interferometry; (120.3940) Metrology; (120.5050) Phase measurement. References and links 1. D. Malacara, M. Servín, and Z. Malacara, Interferogram Analysis For Optical Testing, 2 ed. (CRC Press, 2005. 2. M. Takeda and H. Yamamoto, “Fourier-transform speckle profilometry: three-dimensional shape measurements of diffuse objects with large height steps and/or spatially isolated surfaces,” Appl. Opt. 33(34), 7829–7837 (1994). 3. 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). 4. P. Haible, M. P. Kothiyal, and H. J. Tiziani, “Heterodyne temporal speckle-pattern interferometry,” Appl. Opt. 39(1), 114–117 (2000). 5. P. D. Ruiz, J. M. Huntley, and G. H. Kaufmann, “Adaptive phase-shifting algorithm for temporal phase evaluation,” J. Opt. Soc. Am. A 20(2), 325–332 (2003). 6. J. L. Marroquin, J. E. Figueroa, and M. Servin, “Robust quadrature filters,” J. Opt. Soc. Am. A 14(4), 779–791 (1997). 7. J. C. Estrada, M. Servin, and J. A. Quiroga, “Easy and straightforward construction of wideband phase-shifting algorithms for interferometry,” Opt. Lett. 34(4), 413–415 (2009). 8. M. Servin, J. C. Estrada, and J. A. Quiroga, “The general theory of phase shifting algorithms,” Opt. Express 17(24), 21867–21881 (2009). 9. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20(5), 925–934 (2003). 10. A. Papoulis, Signal analysis, McGraw-Hill (1977) 11. J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. 4, 350–352 (1984). 12. http://www.mathworks.com/help/techdoc/ref/peaks.html Introduction Phase shifting interferometry (PSI) is an experimental technique for phase measurement based on the introduction of a linear phase shift between a set of interferograms. In the temporal case, if we describe the interferogram as a rectangular 2D matrix, the measured intensity at every pixel is given by a set of discrete temporal samples       0 cos 1... , g t b m t t t N       (1) #137146 $15.00 USD Received 25 Oct 2010; revised 22 Dec 2010; accepted 23 Dec 2010; published 2 Mar 2011 (C) 2011 OSA 14 March 2011 / Vol. 19, No. 6 / OPTICS EXPRESS 5002 where b is the background or DC term, m the modulation or AC term, 0  the carrier frequency, which we will assume known, and   t  is the phase at every sample. The goal of this work is the development of a demodulation technique for the reliable measurement of the temporal phase   t  for all sampling points, even in cases involving small numbers of samples or missing data. When the number of samples is small, typically 4 11 N   , we speak of PSI methods [1], but they require a constant phase and do not tolerate missing data. On the other hand, when the number of samples is large, 2 3 10 10 N   , we speak of temporal analysis methods [2–6] which have many problems with the interferogram borders and missing data, or are quadrature filters in an asymptotic way as in [6]. Another possibility for analyzing the temporal signal [Eq. (1)] is the use of a running PSI method tuned at the carrier frequency [3–7]. For example if we use a tree step PSI method we could demodulate the phase locally for each of three consecutive samples. Although this method will deal well with borders, missing data cannot be handled and can even impede the use of this strategy. For example, in the case of an 8 N  samples signal, if samples 3 and 6 are missing, we cannot use a three step PSI for any set of three consecutive samples. In experimental methods missing data appear in the case of a saturated signal and also in heterodyne temporal speckle-pattern interferometry when temporal decorrelation appears. Also all the temporal results presented in this paper can be directly applied to the spatial case where missing data and discontinuities due to occlusions or shadows are very common in projected fringe profilometry. In the literature, all these techniques are presented as different methods, however from the point of view of signal processing they are all quadrature filters (QF). For this reason, we are going to introduce briefly the basics of the quadrature phase detection in the next section and the problems associated with finite length interferograms. Quadrature phase detection of finite interferograms Any linear PSI method can be described as a linear filter, defined by its impulse response   h t or its spectral response   H  [8]. They are related by the discrete time Fourier transform (DTFT) and its inverse given by       it t H DTFT h h t e    


Introduction
Phase shifting interferometry (PSI) is an experimental technique for phase measurement based on the introduction of a linear phase shift between a set of interferograms.In the temporal case, if we describe the interferogram as a rectangular 2D matrix, the measured intensity at every pixel is given by a set of discrete temporal samples for all sampling points, even in cases involving small numbers of samples or missing data.
When the number of samples is small, typically 4 11 N  , we speak of PSI methods [1], but they require a constant phase and do not tolerate missing data.On the other hand, when the number of samples is large, 23 10 10 N , we speak of temporal analysis methods [2-6]  which have many problems with the interferogram borders and missing data, or are quadrature filters in an asymptotic way as in [6].
Another possibility for analyzing the temporal signal [Eq.(1)] is the use of a running PSI method tuned at the carrier frequency [3][4][5][6][7].For example if we use a tree step PSI method we could demodulate the phase locally for each of three consecutive samples.Although this method will deal well with borders, missing data cannot be handled and can even impede the use of this strategy.For example, in the case of an 8 N  samples signal, if samples 3 and 6 are missing, we cannot use a three step PSI for any set of three consecutive samples.In experimental methods missing data appear in the case of a saturated signal and also in heterodyne temporal speckle-pattern interferometry when temporal decorrelation appears.Also all the temporal results presented in this paper can be directly applied to the spatial case where missing data and discontinuities due to occlusions or shadows are very common in projected fringe profilometry.
In the literature, all these techniques are presented as different methods, however from the point of view of signal processing they are all quadrature filters (QF).For this reason, we are going to introduce briefly the basics of the quadrature phase detection in the next section and the problems associated with finite length interferograms.

Quadrature phase detection of finite interferograms
Any linear PSI method can be described as a linear filter, defined by its impulse response   ht or its spectral response

 
H  [8].They are related by the discrete time Fourier transform (DTFT) and its inverse given by where   will be two lobes centred at   and    respectively and confined in the intervals will generate an analytic signal from which it is possible to obtain the modulating phase [8].
Note that in Eq. (3), which, expressed in the temporal domain will be where denotes the convolution product.In Eq. ( 5) we have applied that if the modulating phase is locally monochromatic, the application of a linear filter can be approximated as [10].Equation (5) states that as long as the quadrature condition holds, the phase at every temporal sample can be recovered.However we must remember that we are always dealing with finite sequences, thus the signal is actually the product of an infinite sequence by a window   wt : , where the subscript w stands for windowed.Therefore, we are not measuring In consequence, the signal spectrum is always convolved with the temporal window spectrum.This situation is alleviated in two extreme cases.In the first case, and . In this case the approximation of Eq. ( 5) becomes an exact result given by and we can speak of the monochromatic signal condition.This is the case of all linear PSI methods in which the phase must be constant within the temporal filter extension and the phase can be obtained from   0 qt .In the second extreme case, the signal duration is very long and and we recover an analytical signal like Eq. (5) from which it is possible to measure the phase for all samples.In this situation, we speak of the monochromatic window condition, and this is the typical case analyzed by the temporal methods.
To our knowledge, there is a gap in the literature when we are not in one of these two extreme cases.The method we propose to fill this gap is what we have denominated Regularized Least Squares Quadrature Filter.This method allows the phase recovery for cases that cover from the PSI methods with constant phase and reduced number of samples to the temporal phase measurement with a large number of samples.The proposed Quadrature Filter can also deal with missing data and border effects very efficiently.
Before presenting the Regularized Least Squares Quadrature Filter, we want to discuss a last caveat regarding the demodulation of a temporal phase using quadrature filters.The problem is that for any finite interferogram, it is not possible to demodulate the phase for all points, regardless of the number of samples.Formally speaking, there is a limit in the allowed phase variation that can fit in a given number of samples without aliasing.In classical .However, in our case we are dealing with the general case in which the window spectrum is smearing the signal spectrum, and therefore even for a small bandwidth interferogram we can have a wide spectrum and the classical criterion must be slightly modified.For example, for a rectangular window The DTFT is given by If we use the first zero as the window spectrum extension, the bandwidth of the rectangular window is 2 , in first approximation its bandwidth can be estimated as the sum of which for a big number of samples corresponds to the classical criterion.If we define the ratio it can be used to determine if the temporal phase can be recovered for all samples.Moreover, in our experience, we must impose 2 P  for a reliable phase demodulation.

Regularized least squares phase shifting interferometry
In this work, the objective is to obtain from the temporal signal [Eq.( 1)] an analytical signal given by and from it, calculate the phase  at every sample.
As we mentioned, in a general case in which the signal or the window are not monochromatic, the signal spectrum is smeared.Mathematically that means that we have an ill-posed problem and that we must introduce extra information.Regularized Bayesian estimation theory offers a good framework to deal with this problem [6].Using this idea, it is possible to find a least squares solution to   ft introducing some global constraints like smoothness in the recovered data or continuity in the first derivative.Following this framework we propose to calculate   ft from the samples   gt by minimizing a functional with the form R f, f g is the data term, λ is the regularization parameter and the potential   Vf introduces the apriori knowledge about the solution.Marroquin et al. [6] were the first to explore the idea of using a regularized least squares solution.However, the obtained PSI filters were not truly QF especially for small number of samples and small regularization parameters.Mathematically, the PSI filters presented in [6] were Quadrature Filter only asymptotically for big values of the regularization parameter and number of samples.
In regularization methods, the typical election for the data term is a quadratic cost function where the asterisk indicates complex conjugate and

 
Mt is a binary mask with a valid measurement at sample t.However this data term needs the measurements   gt to be DC filtered, for this reason, in this work we are going to use a DC suppressing data term using first differences. Here,

 
Mt indicates a valid measurement at samples t and 1 t  .In Bayesian estimation theory, two popular potentials are the membrane potential and the thinplate potential The membrane potential enforces continuity in the recovered signal, while the thinplate enforces continuity in the first difference.
The introduction of the conjugate * f in the data term Eq. ( 15) is the main difference with the work of Marroquin [6] and it has important consequences, especially for a small number of samples and low values of the carrier.Qualitatively, we can see that if we introduce the model Eq.(12) in the data term Eqs. ( 14) or (15), assuming a perfect measurement, we have in both cases   * ,0 R f f  , indicating that the minimization of Eq. (13) using Eqs.( 14) or (15) is a true QF.
We can find the solution to the minimization problem Eq. ( 13) setting and treating f and f  as independent variables.In the membrane case, we get the next two difference equations, where where     Finally, solving Eq. ( 21) for

 
F  we obtain the infinite-time spectral response of the membrane solution As can be check using condition (3),

 
M H  represents a DC suppressing QF tuned at 0  , for any combination of λ and 0  .This is what we call Regularized Least Squares Quadrature filter (RLSQF).
In the case of using the thinplate potential, the difference equations are the same as Eq. ( 18) but exchanging The infinite spectral response of the LSRQF filter using the thin-plate potential,

 
ft  in the data term Eq. ( 15) is responsible for this quadrature behaviour, and it is the big difference with the robust quadrature filter (RQF) solution of Marroquin et al [6].For example, using the notation of this work, the infinite spectral response of a DC resistant RQF using a thin-plate potential is From ( 26) we can see that  breaking the quadrature conditions [8].In Fig. 1 we illustrate graphically this effect plotting the infinite frequency response of two DC suppressing RQF and RLSQF filters using both a thin plate potential with tuning frequency 0 2 /   and regularization parameter 1   .As can be seen, the RQF deviates from the quadrature condition for 2 /

 
. The main consequence is the generation of an error in the recovered phase with spatial frequency twice the fringe period [8].The infinite spectral response permits the complete characterization of the RLSQF for sequences that are signal or window monochromatic.However, for intermediate cases with a small number of samples, discontinuities or missing data, it is necessary to work in the direct space and obtain the impulse response for each particular case.This can be done reordering the difference Eqs.(18) and setting proper boundary conditions.In our case we are going to assume the next boundary conditions: . In the case of the DC suppressing data term Eq. ( 15), the mask at 1 t  and tN  is set to zero, M N  If we reorder the difference Eqs.(18), and move to the right side the data we get a linear system given by 2, AL  fg (28 where A is a 22 NN  banded matrix, (29) Solving the linear system, we find the 22 NN  matrix , from which we must extract only the information regarding the A very interesting fact about matrix T Q is that each of its rows can be interpreted as the position varying impulse response for the sample k.Using MATLAB notation we can define , and from it obtain the analytic signal for temporal sample k as The time-varying character of the k h impulse responses responds to the distance to the border and possible missed data.For each impose response k h we can compute a frequency response , and all of them represent true QF filters.For long sequences with no missing data, the frequency response for the central samples will be close to the infinite impulse response.However, near the borders or missing data there can be dramatic changes in the frequency response shape.In Fig. 2, we show two sample-dependent frequency responses for a 15 N  samples DC suppressing RLSQF filter using a thin-plate potential with tuning frequency 0 0.6

 
, regularization parameter 1 Mt for all samples.In particular, we show the frequency response of the samples 1 k  and 8 k  , together with the infinite impulse response given by Eq. ( 26).As can be seen in Fig. 2   It is worth noting that following [6], the presented linear method can be extended to the non-linear case with applications in magnitude equalization or piecewise-smooth QF for discontinuous phase distributions.
We are going to finish this section showing the relation of the proposed method with the classical least squares phase determination [11].If we rewrite Eq. ( 12) as Hence, the proposed method can be interpreted as the regularized extension of the classical least squares phase demodulation in which the 2 a and 3 a coefficients are not constant.In our approach, the regularization term imposes the apriori knowledge about 2 a and 3 a .Although we could use 2 a and 3 a and obtain same results, the complex notation Eq. ( 12) offer more compact results and makes it much easier to complete the frequency analysis and visualization of the quadrature properties of the proposed filter.

Experimental results
The first example is a 15 N  samples signal given by where   gt is given by Eq. ( 1), the modulating phase is nt is a zero-mean normally distributed random variable with standard deviation 5.This is a temporal PSI example with few samples difficult to handle by the Fourier Transform (FT) and RQF methods.As the phase is not constant, this case cannot be handled by standard N samples PSI techniques.Figure 3 shows the results of the demodulation of   n gt using three methods, RQF, FT and the proposed RLSQF.In all the examples of this section, the FT is applied in the transformed space using a DC suppressing QF filter tuned at the interferogram carrier 0  given by   The RQF and the RLSQF are implemented in the direct-space, using the time-varying impulse responses obtained from the rows of matrix T Q (see ref [6] for particular details about the RQF), and the phase is calculated using Eq.(31).In all examples of this section we used a second order thin plate potential with tuning frequency given by the interferogram carrier and regularization parameter 1   .Figure 3a shows the spectra of the interferogram together with the frequency responses of the three demodulation filters.In the case of the RQF and the RLSQF, we are plotting the frequency response for the central sample 8 k  .As can be seen, the signal spectrum is spread due to the small window size.Figure 3b shows the demodulated phase   t  plotted together with the actual phase given by Eq. (34).Although the filters bandwidth match the signal spectrum, the results of the FT and the RQF are incorrect.In the case of the FT, although the quadrature condition is fulfilled, the non-periodic behaviour of   n gt generates the smoothing at the borders.In the case of the RQF, the quadrature condition for 0.75

 
is not accomplished, generating ripples in the recovered phase.In the case of the RLSQF, the boundaries are handled by the time-varying impulse  31)) and the quadrature condition is for all samples, making possible a good demodulation of   t  .In this case, the filling factor is 19 P  (see Eq. ( 11)) indicating a good sampling.Fig. 3. Temporal PSI example for few samples.A) the signal spectrum (dashed line) is plotted together with the frequency response of the FT (blue), RQF (green) and RLSQF (red) filters.In the case of the RQF and LSRQF method, we are plotting the response for the sample k=8.B) Demodulated phase for the FT (blue), RQF (green) and LSRQF (red) methods compared with the actual phase (dashed line).
In the next example, we are going to demonstrate the performance of the RLSQF in dealing with low carrier frequency and missing values.In this case the signal has 60 N  samples and is given by Although the number samples is relatively high, the presence of a missing data and the low carrier makes it very difficult to handle this problem using traditional temporal demodulation techniques.Figure 4 shows the demodulation results obtained using the FT, RQF and RLSQF techniques.Figure 4a shows a plot of the interferogram   n gt with the central part masked.Figure 4b shows the demodulated phase compared with the actual phase (36).As can be seen, not only is the recovered phase correct in the region with valid points, but also, in the masked area, the interpolation is smooth and consistent with the boundary conditions at both sides of it.In this case, the boundary conditions imposed by the processing mask make the FT method unreliable.On the other hand the RQF has a poor quadrature performance for 0.25

 
, generating quadrature errors clearly visible near the masked region.In this case, the filling factor is 5 P  (see Eq. ( 11)) indicating a good sampling.Finally, we are going to present a case in which the interpolation capability suggested in Fig. 4 can be used to demodulate a saturated signal.If we mask the non-linear areas, we can use the RLSQF to demodulate-interpolate the interferogram.Figure 5    0 n gt for those values greater than 255 and smaller than 1. Figure 5 shows a panel of two snapshots of the temporal analysis for 6 t  and 75 t  .Each panel row presents the interferogram, the mask marking saturated areas and the demodulated phase along time.As can be seen the quadrature and interpolating character of the RLSQF permits a good phase recovery even in the case of a saturation problem.For this example the minimum of the fill factor is   min 3 P  .

Conclusions
In this work, we have presented a new quadrature filter using a regularized least squares method.The proposed technique guaranties the quadrature conditions for any combination of carrier frequency, regularization parameter and boundary conditions imposed by the processing mask or missed data.The method is especially useful in temporal PSI when we cannot assume a constant phase or a very long processing window, and with possible missed data.Comparison with existing temporal phase measurement methods demonstrates the performance of the proposed technique.
The MATLAB code to run all the examples of this paper can be downloaded from: goo.gl/RhJq criterion for the allowable interferogram bandwidth becomes 2,

H
 represents a DC suppressing QF for any combination of tuning frequencies and regularization parameter.The presence of

Fig. 2 .
Fig. 2. Frequency responses for a N=15 samples DC suppressing RLSQF filter (see text for details).In this plot, we show the RLSQF frequency responses for samples k=1 (blue) and k=8 (green) together with the infinite frequency response (red).

Fig. 4 .
Fig. 4. Temporal PSI example for few samples and low carrier frequency and missing data.A) plot of the masked signal B) Demodulated phase for the FT (blue), RQF (green) and LSRQF (red) methods compared with the actual phase (dashed black line)

Fig. 5 .
Fig. 5. Temporal demodulation of a time-varying saturated interferogram.From the 150 samples we depict the results for t=6 and t=75 (see text for details).In each row, we depict the instantaneous interferogram, the processing mask and the demodulated phase.
is the background or DC term, m the modulation or AC term, 0 is the phase at every sample.The goal of this work is the development of a demodulation technique for the reliable measurement of the temporal phase   #137146 -$15.00USD Received 25 Oct 2010; revised 22 Dec 2010; accepted 23 Dec 2010; published 2 Mar 2011 (C) 2011 OSA 14 March 2011 / Vol. 19, No. 6 / OPTICS EXPRESS 5002 where b t  t  a big number of samples, the typical criterion is to require that the signal bandwidth must be smaller than the carrier.In a first approximation, the interferogram bandwidth (halfwidth) can be approximated as max (C) 2011 OSA 14 March 2011 / Vol. 19, No. 6 / OPTICS EXPRESS 5004 temporal analysis with g   where the modulating phase is [12]s the demodulation results for a 120×140 time-varying interferogram given by where peaks() is the well known MATLAB test function[12]-.In this example, the AC and DC terms are given by nt is a zero-mean normally distributed random variable with standard deviation 25.To simulate the saturation we have set #137146 -$15.00USD Received 25 Oct 2010; revised 22 Dec 2010; accepted 23 Dec 2010; published 2 Mar 2011 (C) 2011 OSA 14 March 2011 / Vol. 19, No. 6 / OPTICS EXPRESS 5012 and   n gt and