Robust Two-Step phase estimation using the Simplified Lissajous Ellipse Fitting method with Gabor Filter Banks preprocessing

TL;DR

This paper introduces a robust two-step phase estimation method that combines Gabor Filter Banks preprocessing with a simplified ellipse fitting algorithm, effectively handling noise and background variations in interferogram analysis.

Contribution

The paper proposes a novel simplified ellipse fitting approach integrated with Gabor Filter Banks for improved phase estimation in interferometry, reducing computational complexity and enhancing robustness.

Findings

Effective noise filtering and background normalization with Gabor Filter Banks.

Accurate phase step estimation using least squares and robust Leclerc's potential methods.

Validated performance on synthetic and experimental interferogram data.

Abstract

We present the Simplified Lissajous Ellipse Fitting (SLEF) method for the calculation of the random phase step and the phase distribution from two phase-shifted interferograms. We consider interferograms with spatial and temporal dependency of background intensities, amplitude modulations and noise. Given these problems, the use of the Gabor Filters Bank (GFB) allows us to filter--out the noise, normalize the amplitude and eliminate the background. The normalized patterns permit to implement the SLEF algorithm, which is based on reducing the number of estimated coefficients of the ellipse equation, from five terms to only two. Our method consists of three stages. First, we preprocess the interferograms with GFB methodology in order to normalize the fringe patterns. Second, we calculate the phase step by using the proposed SLEF technique and third, we estimate the phase distribution using a two--steps formula. For the calculation of the phase step, we present two alternatives: the use of the Least Squares (LS) method to approximate the values of the coefficients and, in order to improve the LS estimation, a robust estimation based on the Leclerc's potential. The SLEF method's performance is evaluated through synthetic and experimental data to demonstrate its feasibility.

Equations

$$I_k(p)=a_k(p)+b_k(p)\cos [\varphi (p)+{\delta }_k]+{\eta }_k(p),$$

$$I_{add}=I_1+I_2=2a+2b\cos \left(\varphi +\frac{\delta }{2}\right)\cos \left(\frac{\delta }{2}\right)$$

$$I_{sub}=I_1-I_2=2b\sin \left(\varphi +\frac{\delta }{2}\right)\sin \left(\frac{\delta }{2}\right),$$

$${\left(\frac{I_{add}-x_0}{{\alpha }_x}\right)}^2+{\left(\frac{I_{sub}-y_0}{{\alpha }_y}\right)}^2=1$$

$${\theta }_1{x}^2+{\theta }_2{y}^2+{\theta }_{3}x+{\theta }_{4}y+{\theta }_5=0$$

$$\delta =2\arctan \left(\sqrt{\frac{{\theta }_1}{{\theta }_2}}\right).$$

$$\varphi =\arctan \left(\frac{I_{sub}}{I_{add}+\frac{{\theta }_3}{2{\theta }_1}}\sqrt{\frac{{\theta }_2}{{\theta }_1}}\right)-\frac{\delta }{2}.$$

$$GF\{I\}(x,\omega )=I(x)\otimes [e^{-i\omega x}G(x,\sigma )]$$

$${\tilde{I}}_k(x)=GF\{I\}(x,{\omega }_k)=I(x)\otimes [e^{-i{\omega }_{k}x}G(x,{\sigma }_k)],$$

$${\tilde{I}}_k(x)=m_k(x)e^{-i{\psi }_k(x)}$$

$$k^{\ast }(x)=\underset{k}{\mathrm{argmax}}{m}_k(x).$$

$$m(x)=m_{k^{\ast }(x)}(x)$$

$$\psi (x)={\psi }_{k^{\ast }(x)}(x).$$

$$\hat{I}(x)=\cos [\psi (x)].$$

$${\hat{I}}_1(p)$$

$${\hat{I}}_2(p)$$

$${\alpha }_x$$

$${\alpha }_y$$

$${\left(\frac{{\hat{I}}_{add}}{{\alpha }_x}\right)}^2+{\left(\frac{{\hat{I}}_{sub}}{{\alpha }_y}\right)}^2=1,$$

$${\theta }_1{{\hat{I}}_{add}}^2+{\theta }_2{{\hat{I}}_{sub}}^2-1=0.$$

$$x={\hat{I}}_{add}-⟨{\hat{I}}_{add}⟩$$

$$y={\hat{I}}_{sub}-⟨{\hat{I}}_{sub}⟩$$

$${\theta }_1{x}^2+{\theta }_2{y}^2-1=\epsilon (x,y).$$

$$\underset{{\theta }_1,{\theta }_2}{\mathrm{argmin}}\frac{1}{2}\sum _{p\in \mathcal{L}}({\theta }_{1}x(p{)}^2+{\theta }_{2}y(p{)}^2-1{)}^2.$$

$$\underset{T}{\mathrm{argmin}}\frac{1}{2}||X^{T}T-\text{1}|{|}_2^2$$

$$T$$

$$X^T$$

$$T=(X{X}^T{)}^{-1}X\text{1}.$$

$$\phi(p)=\arctan\Bigg{[}\frac{\hat{I}_{1}(p)\cos(\delta)-\hat{I}_{2}(p)}{\hat{I}_{1}(p)\sin(\delta)}\Bigg{]}.$$

$$\underset{{\theta }_1,{\theta }_2}{\mathrm{argmin}}\sum _{p\in \mathcal{L}}\rho ({\theta }_1{x}_i(p{)}^2+{\theta }_2{y}_i(p{)}^2-1;\kappa )$$

$$\rho (z;\kappa )=1-\frac{1}{\kappa }\exp (-\kappa z^2)$$

$$T=(XW{X}^T{)}^{-1}XW\text{1}$$

$$W=\mathrm{diag}[w_1(x_1,y_1),w_2(x_2,y_2),\dots ,w_N(x_N,y_N)]$$

Full text

††journal: Optics Communications

1 Introduction

Phase shifting interferometry is widely used to obtain the phase distribution in interferometric measurements [1, 2, 3]. Even though the calculation can be performed in a single shot [4], the use of several phase shifted interferograms proved to make the measurement more robust to environmental variations [5]. Nowadays, the tendency has been to reduce the number of steps in order to measure dynamic events [6, 7, 8, 9, 10, 11, 12].

One of the main challenges in interferometry is the variation of the parameters of the intensity map. The spatial and temporal dependency of the background intensity, the amplitude modulation and noise are common in non-aligned arrangements [1]. These issues also apply to one-shot interferometry, where the use of optical components such as diffractive devices, polarizers or pixelated masks disturb the captured interferograms [6, 7, 8, 9].

Mathematically, the intensity model of these variable phase–shifted interferograms is given by

[TABLE]

where $k\in 1,2$ is the interferogram index, $p$ is pixel’s coordinates in the regular lattice $\mathcal{L}$, $a$ is the background component, $b$ is the fringe’s amplitude function, $\phi$ is the phase to be recovered, $\delta_{k}$ is the random phase step and $\eta_{k}$ is additive noise. For the case of two–step algorithms, we can assume that $\delta_{1}=0$ and $\delta_{2}=\delta$. The $a$ and $b$ dependencies on $k$ do not allow one to use the known algorithms of phase extraction since they assume temporally constant the background intensity and the amplitude modulation. In such cases, a preprocess is needed in order to normalize the patterns and compute the phase.

In this paper we propose a novel, simplified and more robust version of the Lissajous Ellipse Fitting (LEF) method proposed in [13] for estimating the arbitrary phase step between two phase–shifted fringe patterns with variable parameters. As mentioned before, due to the tendency of reducing the number of steps to estimate the phase distribution, the two–step algorithms have attracted considerably the attention in the past few years; for these reason, several algorithms have been proposed such as Refs. [10, 11, 12, 14, 15, 13, 16, 17, 18, 19, 20, 21, 22]. Particularly, various algorithms based on the LEF method have been designed to improve the estimation of the phase step and the phase distribution in two–step intertferometry, such as: iterative processes based on the least square technique [23], the use of the Gram-Schmidt orthonormalization to transforms the ellipse into a circle [24], the application of a Hilbert-Huang Transform (HHT) pre–filtering with the LEF algorithm [25] or the computation of the Euclidean distance from the points to the ellipse [26]; just to mention some of the novel techniques.

We named our proposed method as Simplified Lissajous Ellipse Fitting (SLEF). Our algorithm reduces the number of estimated coefficients of the ellipse equation, from five terms to only two. Consequently, it is improved the accuracy on the estimation of the relevant parameters by reducing the overfitting of the ellipse to residual noise. Our method consists of three stages: First, we preprocess the fringe patterns using a Gabor Filter Bank (GFB) in order to remove the background variation, normalize the amplitude modulation and filter–out noise [16]. Second, we calculate the phase step through the two term expression of the ellipse equation by using the Least Square (LS) method; alternatively, by minimizing a cost function based on the Leclerc’s potential to define a Robust Estimator (RE) of the coefficients. Third, we calculate the phase distribution using the two–steps algorithm reported in Ref. [27].

We will demonstrate that the use of a GFB as a filtering preprocess does not only improves the robustness of the phase extraction, but also, as it will be demonstrated in Section 3, it simplifies the computation of the ellipse coefficients and consequently the phase step estimation. We remark that we can replace the GFB based preprocessing with other normalization techniques that provide the elimination of the background component, normalize the amplitude modulation and filter–out the noise; for example the Windowed Fourier Transform [28], the HHT [29] or isotropic normalization [30] among others.

Our main contributions are:

A method that uses only 2 parameters that produce equivalent results as the more complex method that uses 5 parameters. 2. 2.

A robust estimation which allows us to overcome large residuals of the pre–filtering process. 3. 3.

We demonstrate by numerical experiments that the LEF methods are more accurate with normalized patterns with GFB than the HHT.

2 Related Methods

In this section, we will present a brief review of the methods to be used in our proposal (which is describe in Section 3): The Lissajous Ellipse Fitting and the Gabor Filters Bank.

2.1 Brief review of the Lissajous Ellipse Fitting (LEF) method

The Lissajous Ellipse Fitting (LEF) method consists on using the Lissajous figure to detect the phase step between two interferograms and estimate the phase [13, 25, 23, 24, 26].

Two phase-shifted interferograms can be represented as the Lissajous figure by plotting their pixel–wise corresponding intensities, see Figure 1. The relation between the major and the minor axes is the result of the phase–shift [13]: if $\delta=\pi/2$, the ellipse would become a circle.

For the case of two-step interferometry, one can consider that the background intensity and the amplitude term are spatially constant and timely invariant, , $a_{1}(p)=a_{2}(p)=a$ and $b_{1}(p)=b_{2}(p)=b$, and that the noise $\eta_{k}(p)$ is filtered–out. Then, by performing the addition and subtraction of these interferograms, one obtains:

[TABLE]

[TABLE]

where the spatial dependency of $\phi$ is omitted in order to simplify the notation.

By solving equations (2) and (3) for $\cos(\phi+\delta/2)$ and $\sin(\phi+\delta/2)$ respectively, and considering that $\cos^{2}(z)+\sin^{2}(z)=1$, one obtains the expression of an ellipse represented as:

[TABLE]

where $x_{0}=2a,y_{0}=0,\alpha_{x}=2b\cos(\delta/2)$ and $\alpha_{y}=2b\sin(\delta/2)$.

Then, one can rewrite (4) in the conical equation of the ellipse:

[TABLE]

where $\theta_{1}=\frac{1}{\alpha_{x}^{2}},\theta_{2}=\frac{1}{\alpha_{y}^{2}},\theta_{3}=-\frac{2x_{0}}{\alpha_{x}^{2}},\theta_{4}=-\frac{2y_{0}}{\alpha_{y}^{2}}$ and $\theta_{5}=\frac{x_{0}^{2}}{\alpha_{x}^{2}}+\frac{y_{0}^{2}}{\alpha_{y}^{2}}-1$. Thus, by solving the coefficients (vector $\theta$) by the least square method (as proposed in [13, 25, 23, 24, 26]), the phase step is computed as

[TABLE]

Hence, the phase distribution is calculated with

[TABLE]

In this work, we present a simplified extension of the LEF method and demonstrate its reliability with complex fringe pattern sets.

2.2 Brief review of Gabor Filters Bank (GFB)

As described by [31, 32, 33, 16, 34, 35], a Gabor Filter ($GF$) is a complex band–pass filter created from the modulation of a complex sinusoidal function with a Gaussian filter ($G$). The complex response from this filter is modeled as:

[TABLE]

where $I$ is the image to be filtered, in this case the fringe patterns, $x$ is the pixel’s index, $\omega$ is the tuned frequency of the filter, $\sigma$ is the Gaussian filter width (window size) and $\otimes$ denotes the convolution of the functions. In terms of frequency, the window size $\sigma$ of the filter represents the width of the band–pass filter centered at the $\omega$ frequency.

A GFB is a set of GFs defined by a set of frequencies $\{\omega_{k}\}_{k=1,2,\ldots}$ and windows sizes $\{\sigma_{k}\}_{k=1,2,\ldots}$. One of the filtered images of the bank would be given as

[TABLE]

which also can be expressed as

[TABLE]

where $m_{k}$ is the magnitude and $\psi_{k}$ is the phase of the response of the image to the $k^{th}$ filter. Then, the filter with the maximum response is estimated at each $x$ index as

[TABLE]

Hence, the local magnitude and phase of the filtered pattern would be

[TABLE]

Finally, the normalized and filtered fringe pattern is given by

[TABLE]

In fact, one could say that trough this method the phase of the fringe pattern is recovered, which is the goal, but the main issue is presented with closed fringes where the sign ambiguity is cannot be solved with only one fringe pattern.

In our synthetic experiments we use images of $256\times 256$ pixels. We use a GFB with ten orientations ($\theta_{k}=k*\pi/10$ for $k=0,2,\dots,9$), four frequencies corresponding to the periods of pixels $\tau=[7,10,15,25]$, a Gaussian window with standard deviation equal the half of the filter period and a window size equal the double of the period. Fig. 2 depicts the real component of the GFs corresponding to the first orientation.

3 Simplified Lissajous Ellipse Fitting (SLEF) Method

3.1 Least Squares based SLEF

Herein we introduce our extension to the LEF algorithm for estimating the actual phase step. We named our variant as SLEF. For this purpose, we will consider the intensity model of the interferograms expressed in (1).

Considering the implementation of the GFBs (explained in section 2.2), we eliminate the background variation, normalize the amplitude modulation and remove the noise (presented in (1)). The ideally normalized two-step interferograms are:

[TABLE]

where by effect of the normalization we have $a_{1}(p)=a_{2}(p)=0$ because of the background elimination, $b_{1}(p)=b_{2}(p)=1$ because of the amplitude normalization and $\eta_{1}(p)=\eta_{2}(p)=0$ because of the filtering process (similar to the one presented in (14)). For the two–step algorithm, $\delta_{1}=0$ and $\delta_{2}=\delta$.

Since we assume removed the background illumination variations, the center of the ellipse is at the origin because $x_{0}=2a=0$ and the eccentricity terms are given by

[TABLE]

Thus, we can simplify the (4) of the ellipse for the Lissajous pattern as

[TABLE]

which corresponds to the equation of the ellipse centered at the origin. Hence, the ellipse’s conical expression is given by

[TABLE]

Note that in this ideal case, (20) is fulfilled for all the pixels. In case that the pre–filtering process does not guarantee an effective normalization of the patterns, the scatter set of points $(\hat{I}_{add},\hat{I}_{sub})$ lay in an uncentered ellipse. Thus, we can center the points with

[TABLE]

where $\langle\cdot\rangle$ denotes the operator that computes the mean. In Figure 4(d) we present a sample of the effect of the filtering process on noisy and non-normalized patterns. The Lissajous ellipse, the line in red, can be considered as the mean of the observed values. Moreover, the spread points around the ideal ellipse are associated to the residual noise ($\varepsilon_{1}$ and $\varepsilon_{2}$) of the pre–filtering process. For this reason, we model such variations with an $\varepsilon$ term in (20):

[TABLE]

Our SLEF method is based in this simplified equation with only two free parameters, $\theta_{1}$ and $\theta_{2}$ and the use of a GFB for normalizing and pre–filtering the FPs. Then, because the presence of the residual $\varepsilon$, we choose to solve the overdetermined system (23) with the Least Square (LS) method

[TABLE]

This expression can be rewritten as

[TABLE]

where, in order to simplify the notation, we define $x_{i}=x(p_{i})$ and $y_{i}=y(p_{i})$. Thus

[TABLE]

and $N=\sharp\mathcal{L}$, 1 is a vector of $N$ entries equal $1$. Thus, the solution to (25) is given by

[TABLE]

Once we have solved system (25) for $\theta_{1}$ and $\theta_{2}$, we are in condition to calculate the phase step $\delta$ with (6). Finally, we compute the phase distribution with the formula for two–step phase shifting reported in Ref. [27]:

[TABLE]

3.2 Robust SLEF

In the previous subsection we estimate the parameters $\theta_{1}$ and $\theta_{2}$ with the LS method. This corresponds to assume a Gaussian distribution for the residual $\varepsilon$. By the examination of Figure 4(d), we noted that such residual is, in fact, non-Gaussian. Therefore, here we propose a robust procedure to improve the estimation of $\theta_{1}$ and $\theta_{2}$. Such robust estimator relies on the fact that the residual distribution has heavy tails [16, 36]. In general, the robust procedure can be formulated as the optimization problem:

[TABLE]

where $\rho$ is a robust potential and $\kappa$ is a positive parameter that controls the outlier rejection sensitivity. In this paper we use the Leclerc’s potential [37]:

[TABLE]

and we set $\kappa=0.1$.

According to [36, 38, 37], the optimization can be obtained by the iteration of the solution of a weighted linear system; i.e.,

[TABLE]

where $W$ is the diagonal matrix of weights

[TABLE]

with

[TABLE]

The solution $T^{\ast}$ is obtained by iterating (31) and (33). The initial conditions for the weight matrix is set $W=\mathop{\mathrm{diag}}[\text{1}]$. We observe that the system converged after just $3$ iterations.

4 Experiments and results

In order to evaluate the proposed algorithms performance, we use ten sets of synthetic patterns of $512\times 512$ with different phase steps $(\delta=[\pi/10,\pi/6,\pi/4,\pi/3,\pi/2])$ and five different Gaussian noise levels $(\sigma=[0.0,0.25,0.5,0.75,1.0])$ . These patterns present spatial and temporal dependency of background intensities, amplitude modulations and noise. In Figure 3 we present the used patterns. Figures 3(a) to 3(j) are the normalized noiseless patterns. It can be observed that the patterns high and low frequencies with circular fringes. Our main interest are these kind of patterns (circular ones) due to the sign error induced in the phase estimation, reason why we are using a two steps algorithm.

Figures 3(k) to 3(t) present the different noise levels applied to each synthetic fringe pattern as well as the background variations and the amplitude modulations. For illustrative purposes we present a different noise level applied to a different pattern, nevertheless, all the noise levels as well as the background variations and amplitude modulations were applied to all the samples.

In order to compare our proposal, we performed the phase step estimation by pre–filtering the synthetic patterns with the HHT as proposed in [25] (which we call LEF–HHT) by using the implementation the Enhaced Fast Empirical Decomposition (EFEMD) proposed in [29, 39].

Figure 4(a) depicts the Lissajous Pattern (LP) for noiseless and normalized fringe patterns with a phase shift of $\delta=\pi/3$; in this case, the LP is centered at the origin and it is perfectly well–marked. In order to remove the rotation of the ellipse, as seen in Figure 1, the LP is computed by using the addition and the subtraction of the fringe patterns [13]. These ideal figures correspond to the pattern presented in Figure 3(j). In Figure 4(b) we present the same pattern with spatial and temporal dependency of the background and amplitude modulation as well as Gaussian noise of $\sigma=0.5$, it is clear that a LP is not appreciated from the original data given the disturbances previously mentioned.

4.1 Pre–filtering process

The first step of our method consists on the preprocessing of the fringe patterns. For comparison purposes, in Figure 4(c) we show the filtered patterns using the HHT with its respective LP. It can be seen that the background term is retrieved but some level of noise remains. In this example, we applied a decomposition of seven modes with a sifting equal to $.01$. On the other hand, the LP pattern presents a dispersed cloud of points, nevertheless, the shape and the eccentricity of the ellipse is still noticeable. Finally, in Figure 4(d) we present the filtered pattern using the GFB and its respective LP. Even though the noise filtering is better as well as the normalization, there are several residuals as stablished in (23).

From obtained LP of the GFB filtered patterns, it is important to emphasize the following:

The LP center is at the origin. 2. 2.

The approximation of the points is close to the ideal ellipse. 3. 3.

The spread points are due to the residuals of the preprocess.

4.2 Phase step calculation

The second step is the estimation of the phase step $\delta$ using (6). To evaluate the feasibility of our proposed methods, we calculated the phase step estimation comparing them with the preprocessed LEF algorithm proposed in [25]. The first comparison consists on calculating the phase step at different Gaussian noise levels. For this analysis, we used ten different sets of images (See Figure 3) with five different Gaussian noise levels $(\sigma=[0.0,0.25,0.5,0.75,1.0])$, with spatial–temporal variations in their background intensity as well as the amplitude modulation. The phase step between the patterns was set to $\pi/3$. Each pattern was preprocessed using the HHT and the calculation of the phase step was using the same (6) using the computed parameters of each method: LEF–HHT, SLEF–LS and SLEF–RE .

Figure 5 depicts the computed Mean Absolute Error (MAE) of the analyzed sets of patterns pre–filtered with the HHT method. The total mean error of the algorithms are: $MAE_{LEF-HHT}=0.1021rad$, $MAE_{SLEF-LS}=0.1018rad$ and $MAE_{SLEF-RE}=0.0991rad$. It can be appreciated that the algorithms present similar behavior at high noise levels while the SLEF algorithms reduce the variance at low noise levels. The MAE in noise levels superior to $\sigma=0.5$ is around $0.1rad$ which could generate the presence of harmonics in the recovered phase.

Now, we pre–filtered the same synthetic patterns with the GFB to demonstrate the feasibility of this technique. The tested algorithms are the 5–term LEF which we will call LEF–GFB and the prosed SLEF–LS and SLEF–RE.

In Figure 6 we present the MAE of the analyzed sets of pre–filtered patterns. The total mean error of the algorithms are: $MAE_{LEF-GFB}=0.0869rad$, $MAE_{SLEF-LS}=0.0871rad$ and $MAE_{SLEF-RE}=0.0204rad$. It can be appreciated that the well-known LEF algorithm and the proposed SLEF–LS algorithm have practically the same behavior, meaning that if the patterns are well normalized, these algorithms are equivalent; with the SLEF–LS’s advantage of having simplified solution. Moreover, SLEF–RE is a more accurate (less error) and precise (less error variance) phase step estimator. Regardless of the noise presented, the MAE is smaller than $0.025rad$ in all levels.

The results presented in Figure 6 prove that the GFB pre–filtering is more robust to high noise levels, which improve the stability of the algorithms. On lower noise levels the performance is similar on both pre–filtering process. The SLEF–RE algorithm has the best performance when using GFB as normalizing method.

To evaluate the estimation of our method, we also performed a test consisting on the calculation of the phase step with different phase steps $\delta$. For this analysis, we used the same images with a fixed Gaussian noise of $\sigma=0.5$. The phase steps between the patterns were $\delta=[\pi/10,\pi/6,\pi/4,\pi/3,\pi/2]$. As before, each pattern was preprocessed using the GFB and the calculation of the phase step was using the same (6) with the computed parameters of each method: SLEF–LS, SLEF–RE and LEF–GFB.

Figure 7 shows the MAE resulted of the analyzed patterns. It can be observed that the LEF–GFB and SLEF–LS algorithms present the same behavior, proving their equivalency, while SLEF–RE have smaller variance and error. It is important to notice that the error increases for small phase steps, and it decreases as the step approaches to $\delta=\pi/2$. The error for phase steps in the interval $\delta\in(\pi/2,\pi)$ is the same as the presented in Figure 7, since the error tends to increase as the step gets closer to $\pi$.

We include Figure 8 to demonstrate the behavior of the LEF–HHT algorithm and our GFB based proposals. As mentioned before, the LEF–HHT performs better than the SLEF–LS in low noise level. The SLEF–RE demonstrated to be de most accurate the most accurate in all noise levels.

4.3 Phase map estimation

Finally, the third stage is the calculation of the phase using Eq. (7) for the LEF algorithm and Eq. (28) for the SLEF algorithms. To evaluate the phase error, we used the pattern presented in Figure 4. In this case, the phase shift is $\delta=\pi/3$ and the Gaussian noise presents a $\sigma=0.5$. The purpose of this experiment is to prove that the piston term induced in (7) can be avoided by using (28) instead. The fringe patterns were normalized with GFBs and since we already demonstrated that SLEF–LS and LEF are equivalent given such condition, we assume that the phase estimation using (28) would result the same with 2–terms or 5–terms.

In Figure 9 we present the calculation of one of the synthetic patterns shown in Figure 4(b). Figures 9(a) to 9(e) present the estimated phases. Figures 9(g) to 9(j) correspond to the wrapped error (which corresponds to the harmonics generated) between of the calculated phases with respect to the ideal phase . Figures 9(k) to 9(o) present the same wrapped error maps with contour lines for visualization purposes.

The first column corresponds to the results obtained using the LEF–HHT as proposed in [25]. The second column corresponds to the LEF–HHT algorithm but using (28) to estimate the phase map. It can be seen that the amount of errors is drastically reduced due to the use of a different equation for the phase estimation. In the third column we present the results of using the SLEF–HHT algorithm which is the use of the SLEF–RE algorithm but using the HHT as normalization method. Finally, the fourth and fifth columns correspond to the SLEF–LS and SLEF– RE methods with GFB pre–filtering.

For this particular experiment, the error of the estimation of the step using the LEF–HHT algorithm was $\delta_{error}=0.1194rad$ while the error for the SLEF–LS algorithm was $\delta_{error}=0.0567rad$ and the SLEF–RE was $\delta_{error}=0.0219rad$. The SLEF–RE method with HHT pre–filtering presented an error of $\delta_{error}=0.1092rad$.

The MAEs of the error surfaces are $MAE_{LEF-HHT}=1.1712$, $MAE_{LEF-HHT-2}=0.4722rad$, $MAE_{SLEF-HHT}=0.4846rad$, $MAE_{SLEF-LS}=0.37rad$ and $MAE_{SLEF-RE}=0.3569rad$. It is important to note that the high amount of harmonics in Figure 9(h) is due to used phase extraction equation which includes a piston term, in this case Eq. (7). If the phase was recovered by (28), the MAE would be reduced as presented in Figure 9(l).

4.4 Experimental results

In order to test the capabilities of our algorithm we implemented a Polarizing Cyclic Path Interferometer (PCPI) as the one proposed in [40] in its radial mode. Such system was implemented with a diode Laser with power of $400mW$ operating at $\lambda=532nm$. In the arrangement, a polarizer filter is placed at an angle of $45^{\circ}$ at the entrance of the interferometer, so the incoming beam will have perpendicular and a parallel components of the same intensity. Once the light has gone through the PCPI, we placed a quarter wave plate which generates crossed circular polarization states. At the output of the PCPI, an interference pattern is observed when placing an auxiliary linear polarizer, according the next equation

[TABLE]

which is in fact the same as (1), but in this case, $\Psi_{k}$ corresponds to the angle of the linear polarizer. This rotation produces the phase shift between the two interferograms as $\delta_{k}=2\Psi_{k}$.

For this experiment, the auxiliar linear polarizer was mounted on a graduated rotational mount. The two interferograms were recorded with a resolution of $410\times 410$ and a shift angle of $30^{\circ}$ between the captures, which means a phase step of $\delta=\pi/3$. Figures 10(a) and 10(b) present such patterns.

Figures 10(d) and 10(e) show the pre–filtered patterns using the HHT. Such result was obtained using the EFEMD algorithm with its automatic selection of modes. On the other hand, Figures 10(g) and 10(h) illustrate the normalized patterns using GFB. For these results applied a Winner–Take–All (WTA) strategy for the GFB corresponding to the periods equal to $[20,35,45,55]$ and eight orientations. Since the main problem of the GFB is the detection of low frequencies, we use

[TABLE]

where $\alpha=Mag/\max(Mag)$, $Mag$ is the magnitude map of the response of the GFB and $\tilde{I}_{k}$ is low–pass filtered version of the original image.

To prove the demonstrate the feasibility of our proposal, we present the phase obtained with the LEF–HHT, SLEF–LS and the SLEF–RE algorithms.

Figure 11(a) presents the estimated phase map using the LEF–HHT algorithm, were the pre–filtering was made with the HHT (Figures 10(d) and 10(e)). The estimation of the phase step by using the 5–term equation and the LS method, which gave a result of $\delta=1.2186$, a relative error of $16.4\%$. The calculation of the phase map was made with (7).

Figure 11(b) presents the estimated phase map using the SLEF–LS algorithm, were the pre–filtering was made with the GFB (Figures 10(g) and 10(h)). The estimation of the phase step by using the 2–term equation and the LS method (which under these conditions have the same behavior as the 5–term), which gave a result of $\delta=1.1673$, a relative error of $11.4\%$. The calculation of the phase map was made with (28).

Figure 11(c) presents the estimated phase map using the SLEF–RE algorithm, were the pre–filtering was made with the GFB (Figures 10(g) and 10(h)). The estimation of the phase step by using the 2–term equation and the robust estimator, which gave a result of $\delta=1.115$, a relative error of $6.4\%$. The calculation of the phase map was made with (28).

Figure 11 presents the estimated phase map using the LEF-HHT, SLEF–LS and SLEF–RE. For purposes of a qualitative comparison, we include the phase estimated with the well-known 4 steps algorithm, shown in Figure 11(d). The fringe patterns of this result were only filtered with a simple Gaussian filter of $\sigma=2$.

From a visual inspection, one can note that SLEF–LS and SLEF–RE algorithms present practically the same results. Meanwhile, LEF–HHT produces a result with a non smooth gray–scale transition; typical of a detunning. It is important to note the capabilities of the GFB of filtering out the noise, given that the obtained phase is clearly comparable to the 4–steps algorithm (which is more robust to noise than a two–steps method).

5 Discussions

As mentioned in section 2.1, two phase-shifted interferograms can be represented as the Lissajous figure by plotting their pixel–wise corresponding intensities. Nevertheless, this algorithm requires at least that the background intensity as well as the amplitude modulation to be temporary constants [13, 25, 23, 24, 26]. If this condition is not presented, it is required a preparation process for the interferograms in order to accomplish a constant modulation of the fringes. Preprocessing techniques such as Windowed Fourier Transform [28], GFB [16], the Hilbert–Huang Transform (HHT) [17] or isotropic normalization [30] allow to eliminate such variations. Our work presents the advantage of using normalized fringe patterns since it allows us to simplify the solution of the ellipse equation from a 5–terms equation to a 2–terms form. Liu et. al [25] present a similar approach even though, they still solve a 5–terms equation. We compare our approach with theirs in order to prove the equivalency of our 2–terms SLEF–LS algorithm with the 5–term pre–filtered technique. On the other hand, the SLE–RE algorithm present more robustness to the residuals of the pre–filtering process, giving more stability and accuracy in the phase step calculation process.

6 Conclusions

We introduced a simplified model of the Lissajous Ellipse Fitting to calculate the phase step and phase distribution of two randomly shifted interferograms. We focused on solving the problem of obtaining the phase of interferograms with spatial–temporal dependencies on their background intensities, amplitude modulations and noise. The main advantages of use of the GFB is that the phase estimation is robust to the mentioned issues since it filters–out the noise, normalizes the amplitude and eliminates the background. Given the normalized patterns, the ellipse equation can be simplified to a two–unknowns system instead of a five–unknowns system, we named this the SLEF algorithm. Our method consists of three stages: the preprocess the fringe patterns using a GFB, the estimation of the phase step through the estimation of the coefficients of the ellipse’s equation and the calculation the phase distribution. We remark that we can replace the GFB based preprocessing with other normalization techniques that provide the elimination of the background component, normalize the amplitude modulation and filter–out the noise. As presented, the estimation of the two terms of the ellipse equation can be done by the well–known LS method, which results are practically the same as the 5–terms. Also, we introduced a novel implementation of a robust estimator such as the Leclerc’s potential in order to improve the accuracy of the phase step estimation, this is mainly caused by residuals of the filtering process. The experimental results of the calculation of our algorithms to 100 pairs of images (10 different patterns with 5 different levels of noise and 5 different phase steps) prove that SLEF–LS is equivalent to the pre–filtered 5–term LEF algorithm and the robustness to residuals of our SLEF–RE algorithm, which improves significantly the accuracy of the phase step estimation.

Acknowledgements

VHFM thanks Consejo Nacional de Ciencia y Tecnología (Conacyt) for the provided postdoctoral grant. This research was supported in part by Conacyt, Mexico (Grant A1-S-43858) and the NVIDIA Academic program.