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

Rob ust T wo-Step phase estimation using the Simplified Lissajous Ellipse Fitting method with Gabor Filters Bank preprocessing V ´ ıctor H. Flores a,b, ∗ , Mariano Riv era a a Centr o de In vestigaci ´ on en Matem ´ aticas AC, 36023, Guanajuato, Gto., M ´ exico. b Departamento de Ingenier ´ ıa Rob ´ otica, Universidad P olit ´ ecnica del Bicentenario, 36283, Silao, Gto., M ´ exico. Abstract W e present the Simplified Lissajous Ellipse Fitting (SLEF) method for the calculation of the random phase step and the phase distribution from tw o phase-shifted interferograms. W e consider interferograms with spatial and temporal dependency of back- ground intensities, amplitude modulations and noise. Giv en 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 coe ffi cients of the ellipse equation, from fi ve 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 distrib ution using a two–steps formula. For the calculation of the phase step, we present two alternati ves: the use of the Least Squares (LS) method to approximate the v alues of the coe ffi cients and, in order to improv e the LS estimation, a robust estimation based on the Leclerc’ s potential. The SLEF method’ s performance is ev aluated through synthetic and experimental data to demonstrate its feasibility . K e ywor ds: T w o-Step, Phase Shifting, Lissajous Ellipse Fitting, Gabor Filters Bank 2010 MSC: 78-04, 78-05 1. Introduction Phase shifting interferometry is widely used to obtain the phase distribution in interferometric measurements [1 – 3]. Even though the calculation can be performed in a single shot [4], the use of sev eral phase shifted interferograms prov ed to make the mea- surement more robust to en vironmental v ariations [5]. No wa- days, the tendency has been to reduce the number of steps in order to measure dynamic ev ents [6 – 12]. One of the main challenges in interferometry is the v ariation of the parameters of the intensity map. The spatial and temporal dependency of the background intensity , the amplitude modu- lation 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 di ff ractiv e de vices, polarizers or pixelated masks disturb the captured interferograms [6 – 9]. Mathematically , the intensity model of these variable phase– shifted interferograms is giv en by I k ( p ) = a k ( p ) + b k ( p ) cos[ φ ( p ) + δ k ] + η k ( p ) , (1) where k ∈ 1 , 2 is the interferogram index, p is pixel’ s coordi- nates in the regular lattice L , a is the background component, b is the fringe’ s amplitude function, φ is the phase to be recov- ered, δ k is the random phase step and η k is additi ve noise. For ∗ Corresponding author Email addr ess: victor.flores@cimat.mx (V ´ ıctor H. Flores) the case of two–step algorithms, we can assume that δ 1 = 0 and δ 2 = δ . 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 ro- bust v ersion of the Lissajous Ellipse Fitting (LEF) method pro- posed 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 algo- rithms have attracted considerably the attention in the past few years; for these reason, sev eral algorithms ha ve been proposed such as Refs. [10 – 22]. Particularly , various algorithms based on the LEF method have been designed to improv e the estima- tion of the phase step and the phase distribution in two–step intertferometry , such as: iterativ e processes based on the least square technique [23], the use of the Gram-Schmidt orthonor- malization to transforms the ellipse into a circle [24], the appli- cation of a Hilbert-Huang Transform (HHT) pre–filtering with the LEF algorithm [25] or the computation of the Euclidean dis- tance from the points to the ellipse [26]; just to mention some of the nov el techniques. W e named our proposed method as Simplified Lissajous El- lipse Fitting (SLEF). Our algorithm reduces the number of es- timated coe ffi cients of the ellipse equation, from fi ve terms to only two. Consequently , it is improv ed the accuracy on the es- Pr eprint submitted to Optics Communications J anuary 29, 2020 timation of the relev ant 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 ex- pression of the ellipse equation by using the Least Square (LS) method; alternativ ely , by minimizing a cost function based on the Leclerc’ s potential to define a Robust Estimator (RE) of the coe ffi cients. Third, we calculate the phase distrib ution using the two–steps algorithm reported in Ref. [27]. W e 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 coe ffi cients and con- sequently the phase step estimation. W e remark that we can replace the GFB based preprocessing with other normalization techniques that pro vide the elimination of the background com- ponent, normalize the amplitude modulation and filter –out the noise; for example the W indowed F ourier T ransform [28], the HHT [29] or isotropic normalization [30] among others. Our main contributions are: 1. A method that uses only 2 parameters that produce equi v- alent results as the more complex method that uses 5 pa- rameters. 2. A robust estimation which allo ws us to ov ercome large residuals of the pre–filtering process. 3. W e 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 revie w 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 r evie w of the Lissajous Ellipse F itting (LEF) method The Lissajous Ellipse Fitting (LEF) method consists on using the Lissajous figure to detect the phase step between two inter- ferograms and estimate the phase [13, 23 – 26]. T wo phase-shifted interferograms can be represented as the Lissajous figure by plotting their pixel–wise corresponding in- tensities, see Figure 1. The relation between the major and the minor axes is the result of the phase–shift [13]: if δ = π/ 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 spa- tially constant and timely inv ariant, , a 1 ( p ) = a 2 ( p ) = a and b 1 ( p ) = b 2 ( p ) = b , and that the noise η k ( p ) is filtered–out. Then, by performing the addition and subtraction of these in- terferograms, one obtains: I add = I 1 + I 2 = 2 a + 2 b cos  φ + δ 2  cos  δ 2  (2) Figure 1: Resulting Lissajous ellipse of mapping pixel-wise the interferograms corresponding intensities of two ideal interferograms. I sub = I 1 − I 2 = 2 b sin  φ + δ 2  sin  δ 2  , (3) where the spatial dependency of φ is omitted in order to sim- plify the notation. By solving equations (2) and (3) for cos( φ + δ/ 2) and sin( φ + δ/ 2) respectiv ely , and considering that cos 2 ( z ) + sin 2 ( z ) = 1, one obtains the expression of an ellipse represented as: I add − x 0 α x ! 2 + I sub − y 0 α y ! 2 = 1 (4) where x 0 = 2 a , y 0 = 0 , α x = 2 b cos( δ/ 2) and α y = 2 b sin( δ/ 2). Then, one can rewrite (4) in the conical equation of the el- lipse: θ 1 x 2 + θ 2 y 2 + θ 3 x + θ 4 y + θ 5 = 0 (5) where θ 1 = 1 α 2 x , θ 2 = 1 α 2 y , θ 3 = − 2 x 0 α 2 x , θ 4 = − 2 y 0 α 2 y and θ 5 = x 2 0 α 2 x + y 2 0 α 2 y − 1. Thus, by solving the coe ffi cients (v ector θ ) by the least square method (as proposed in [13, 23 – 26]), the phase step is computed as δ = 2 arctan       r θ 1 θ 2       . (6) Hence, the phase distribution is calculated with φ = arctan         I sub I add + θ 3 2 θ 1 r θ 2 θ 1         − δ 2 . (7) In this work, we present a simplified extension of the LEF method and demonstrate its reliability with complex fringe pat- tern sets. 2.2. Brief re view of Gabor F ilters Bank (GFB) As described by [16, 31–35], a Gabor Filter ( G F ) is a com- plex band–pass filter created from the modulation of a comple x sinusoidal function with a Gaussian filter ( G ). The complex response from this filter is modeled as: G F { I } ( x , ω ) = I ( x ) ⊗ [ e − i ω x G ( x , σ )] (8) where I is the image to be filtered, in this case the fringe pat- terns, x is the pixel’ s index, ω is the tuned frequency of the filter , σ is the Gaussian filter width (window size) and ⊗ denotes the con volution of the functions. In terms of frequency , the window 2 size σ of the filter represents the width of the band–pass filter centered at the ω frequency . A GFB is a set of GFs defined by a set of frequencies { ω k } k = 1 , 2 ,... and windo ws sizes { σ k } k = 1 , 2 ,... . One of the filtered images of the bank would be gi ven as ˜ I k ( x ) = G F { I } ( x , ω k ) = I ( x ) ⊗ [ e − i ω k x G ( x , σ k )] , (9) which also can be expressed as ˜ I k ( x ) = m k ( x ) e − i ψ k ( x ) (10) where m k is the magnitude and ψ 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 k ∗ ( x ) = argmax k m k ( x ) . (11) Hence, the local magnitude and phase of the filtered pattern would be m ( x ) = m k ∗ ( x ) ( x ) (12) ψ ( x ) = ψ k ∗ ( x ) ( x ) . (13) Finally , the normalized and filtered fringe pattern is gi ven by ˆ I ( x ) = cos[ ψ ( x )] . (14) In fact, one could say that trough this method the phase of the fringe pattern is recov ered, 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 × 256 pixels. W e use a GFB with ten orientations ( θ k = k ∗ π/ 10 for k = 0 , 2 , . . . , 9), four frequencies corresponding to the periods of pixels τ = [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 compo- nent of the GFs corresponding to the first orientation. Figure 2: Real component of the GFs corresponding to the first orientations and periods equal to [7 , 10 , 15 , 25] pixels, respecti vely . 3. Simplified Lissajous Ellipse Fitting (SLEF) Method 3.1. Least Squares based SLEF Herein we introduce our extension to the LEF algorithm for es- timating the actual phase step. W e 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 (e xplained in section 2.2), we eliminate the background variation, normalize the amplitude modulation and remov e the noise (presented in (1)). The ideally normalized two-step interferograms are: ˆ I 1 ( p ) = cos[ φ ( p )] (15) ˆ I 2 ( p ) = cos[ φ ( p ) + δ ] , (16) where by e ff ect of the normalization we hav e 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 η 1 ( p ) = η 2 ( p ) = 0 because of the filtering process (similar to the one presented in (14)). For the two–step algorithm, δ 1 = 0 and δ 2 = δ . Since we assume removed the background illumination v ari- ations, the center of the ellipse is at the origin because x 0 = 2 a = 0 and the eccentricity terms are gi ven by α x = 2 cos  δ 2  (17) α y = 2 sin  δ 2  . (18) Thus, we can simplify the (4) of the ellipse for the Lissajous pattern as ˆ I add α x ! 2 + ˆ I sub α y ! 2 = 1 , (19) which corresponds to the equation of the ellipse centered at the origin. Hence, the ellipse’ s conical expression is giv en by θ 1 ˆ I 2 add + θ 2 ˆ I 2 sub − 1 = 0 . (20) Note that in this ideal case, (20) is fulfilled for all the pix- els. In case that the pre–filtering process does not guarantee an e ff ectiv e normalization of the patterns, the scatter set of points ( ˆ I add , ˆ I sub ) lay in an uncentered ellipse. Thus, we can center the points with x = ˆ I add − h ˆ I add i (21) y = ˆ I sub − h ˆ I sub i (22) where h·i denotes the operator that computes the mean. In Fig- ure 4d we present a sample of the e ff ect 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. Moreo ver , the spread points around the ideal ellipse are associated to the residual noise ( ε 1 and ε 2 ) of the pre–filtering process. For this reason, we model such variations with an ε term in (20): θ 1 x 2 + θ 2 y 2 − 1 = ε ( x , y ) . (23) Our SLEF method is based in this simplified equation with only two free parameters, θ 1 and θ 2 and the use of a GFB for 3 normalizing and pre–filtering the FPs. Then, because the pres- ence of the residual ε , we choose to solve the overdetermined system (23) with the Least Square (LS) method arg min θ 1 ,θ 2 1 2 X p ∈L ( θ 1 x ( p ) 2 + θ 2 y ( p ) 2 − 1) 2 . (24) This expression can be re written as arg min T 1 2 || X T T − 1 || 2 2 (25) where, in order to simplify the notation, we define x i = x ( p i ) and y i = y ( p i ). Thus T de f = [ θ 1 , θ 2 ] T (26a) X T de f =                  x 2 1 y 2 1 x 2 2 y 2 2 . . . . . . x 2 N y 2 N                  (26b) and N = ] L , 1 is a v ector of N entries equal 1. Thus, the solution to (25) is giv en by T = ( X X T ) − 1 X 1 . (27) Once we hav e solved system (25) for θ 1 and θ 2 , we are in condition to calculate the phase step δ with (6). Finally , we compute the phase distrib ution with the formula for two–step phase shifting reported in Ref. [27]: φ ( p ) = arctan " ˆ I 1 ( p ) cos( δ ) − ˆ I 2 ( p ) ˆ I 1 ( p ) sin( δ ) # . (28) 3.2. Robust SLEF In the pre vious subsection we estimate the parameters θ 1 and θ 2 with the LS method. This corresponds to assume a Gaussian distribution for the residual ε . By the examination of Figure 4d, we noted that such residual is, in fact, non-Gaussian. Therefore, here we propose a robust procedure to improve the estimation of θ 1 and θ 2 . Such robust estimator relies on the fact that the residual distrib ution has hea vy tails [16, 36]. In general, the ro- bust procedure can be formulated as the optimization problem: arg min θ 1 ,θ 2 X p ∈L ρ ( θ 1 x i ( p ) 2 + θ 2 y i ( p ) 2 − 1; κ ) (29) where ρ is a robust potential and κ is a positi ve parameter that controls the outlier rejection sensitivity . In this paper we use the Leclerc’ s potential [37]: ρ ( z ; κ ) = 1 − 1 κ exp( − κ z 2 ) (30) and we set κ = 0 . 1. According to [36 – 38], the optimization can be obtained by the iteration of the solution of a weighted linear system; i.e. , T = ( X W X T ) − 1 X W 1 (31) where W is the diagonal matrix of weights W = diag[ w 1 ( x 1 , y 1 ) , w 2 ( x 2 , y 2 ) , . . . , w N ( x N , y N )] (32) with w i ( x i , y i ) = e xp( − 2 k [ θ 1 x 2 i + θ 2 y 2 i − 1] 2 ) . (33) The solution T ∗ is obtained by iterating (31) and (33). The initial conditions for the weight matrix is set W = diag[ 1 ]. W e observe that the system con ver ged after just 3 iterations. 4. Experiments and results In order to e v aluate the proposed algorithms performance, we use ten sets of synthetic patterns of 512 × 512 with di ff erent phase steps ( δ = [ π/ 10 , π/ 6 , π/ 4 , π/ 3 , π/ 2]) and fiv e di ff erent Gaussian noise le vels ( σ = [0 . 0 , 0 . 25 , 0 . 5 , 0 . 75 , 1 . 0]) . These patterns present spatial and temporal dependenc y of background intensities, amplitude modulations and noise. In Figure 3 we present the used patterns. Figures 3a to 3j 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 in- duced in the phase estimation, reason why we are using a two steps algorithm. Figures 3k to 3t present the di ff erent noise levels applied to each synthetic fringe pattern as well as the background v aria- tions and the amplitude modulations. For illustrati ve purposes we present a di ff erent noise le vel applied to a di ff erent pattern, nev ertheless, all the noise lev els as well as the background v ari- ations and amplitude modulations were applied to all the sam- ples. 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 4a depicts the Lissajous P attern (LP) for noiseless and normalized fringe patterns with a phase shift of δ = π/ 3; in this case, the LP is centered at the origin and it is perfectly well–marked. In order to remov e 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 3j. In Figure 4b we present the same pattern with spatial and temporal depen- dency of the background and amplitude modulation as well as Gaussian noise of σ = 0 . 5, it is clear that a LP is not appre- ciated from the original data gi ven the disturbances previously mentioned. 4.1. Pre–filtering pr ocess The first step of our method consists on the preprocessing of the fringe patterns. For comparison purposes, in Figure 4c we sho w the filtered patterns using the HHT with its respec- tiv e LP . It can be seen that the background term is retrie ved b ut 4 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) σ = 0 (l) σ = 0 . 25 (m) σ = 0 . 5 (n) σ = 0 . 75 (o) σ = 1 . 0 (p) σ = 0 (q) σ = 0 . 25 (r) σ = 0 . 5 (s) σ = 0 . 75 (t) σ = 1 . 0 Figure 3: Synthetic patterns. some level of noise remains. In this example, we applied a de- composition of seven modes with a sifting equal to . 01. On the other hand, the LP pattern presents a dispersed cloud of points, nev ertheless, the shape and the eccentricity of the ellipse is still noticeable. Finally , in Figure 4d we present the filtered pattern using the GFB and its respecti ve LP . Even though the noise fil- tering 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 impor- tant to emphasize the following: 1. The LP center is at the origin. 2. The approximation of the points is close to the ideal el- lipse. 3. The spread points are due to the residuals of the prepro- cess. 4.2. Phase step calculation The second step is the estimation of the phase step δ using (6). T o e valuate the feasibility of our proposed methods, we calculated the phase step estimation comparing them with the preprocessed LEF algorithm proposed in [25]. The first com- parison consists on calculating the phase step at di ff erent Gaus- sian noise le vels. For this analysis, we used ten di ff erent sets of images (See Figure 3) with fi ve di ff erent Gaussian noise le v- els ( σ = [0 . 0 , 0 . 25 , 0 . 5 , 0 . 75 , 1 . 0]), with spatial–temporal v ari- ations in their background intensity as well as the amplitude modulation. The phase step between the patterns was set to π/ 3. Each pattern w as 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 . (a) Ideal patterns (b) Experimental patterns (c) Filtered patterns using HHT (d) Filtered patterns using GFB Figure 4: Lissajous figures for di ff erent interferogram cases. The phase step between the patterns is δ = π/ 3. 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: M A E LE F − H H T = 0 . 1021 rad , M A E S L E F − LS = 0 . 1018 rad and M AE S L E F − RE = 0 . 0991 rad . It can be appreciated that the algorithms present similar behavior at high noise lev els while the SLEF algorithms reduce the variance at low noise le vels. The MAE in noise lev- els superior to σ = 0 . 5 is around 0 . 1 r ad which could generate the presence of harmonics in the recov ered 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: M A E LE F − G F B = 0 . 0869 rad , M A E S L E F − LS = 0 . 0871 rad and M A E S L E F − RE = 0 . 0204 rad . It can be appreciated that the well-known LEF algorithm and the proposed SLEF–LS algo- rithm have practically the same behavior , meaning that if the patterns are well normalized, these algorithms are equi v alent; 5 Figure 5: Mean Absolute Error of the phase step calculation using HHT prepro- cess for the LEF–HHT , SLEF–LS and SLEF–RE algorithms at di ff erent noise lev els. Figure 6: Mean Absolute Error of the phase step calculation using GFB prepro- cess for the LEF–GFB, SLEF–LS and SLEF–RE algorithms at di ff erent noise lev els. with the SLEF–LS’ s advantage of having simplified solution. Moreov er , SLEF–RE is a more accurate (less error) and pre- cise (less error variance) phase step estimator . Regardless of the noise presented, the MAE is smaller than 0 . 025 rad in all lev els. The results presented in Figure 6 prov e that the GFB pre– filtering is more robust to high noise levels, which improve the stability of the algorithms. On lower noise lev els the perfor- mance is similar on both pre–filtering process. The SLEF–RE algorithm has the best performance when using GFB as nor- malizing method. T o e v aluate the estimation of our method, we also performed a test consisting on the calculation of the phase step with dif- ferent phase steps δ . For this analysis, we used the same im- ages with a fix ed Gaussian noise of σ = 0 . 5. The phase steps between the patterns were δ = [ π/ 10 , π/ 6 , π/ 4 , π/ 3 , π/ 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 , pro ving their equi v alency , while SLEF– RE ha ve 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 δ = π/ 2. The error for phase steps in the Figure 7: Mean Absolute Error of the phase step calculation for the LEF–GFB SLEF–LS and SLEF–RE algorithms at di ff erent phase steps. interval δ ∈ ( π/ 2 , π ) is the same as the presented in Figure 7, since the error tends to increase as the step gets closer to π . Figure 8: Comparison of the 5–term LEF algorithm with HHT pre–filtering with the proposed algorithms SLEF–LS and SLEF–RE. W e include Figure 8 to demonstrate the behavior of the LEF–HHT algorithm and our GFB based proposals. As men- tioned 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 lev els. 4.3. Phase map estimation Finally , the third stage is the calculation of the phase us- ing Eq. (7) for the LEF a lgorithm and Eq. (28) for the SLEF algorithms. T o ev aluate the phase error , we used the pattern presented in Figure 4. In this case, the phase shift is δ = π/ 3 and the Gaussian noise presents a σ = 0 . 5. The purpose of this experiment is to prov e 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 equiv alent given such condition, we as- sume 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 4b. Figures 9a to 9e present the esti- mated phases. Figures 9g to 9j correspond to the wrapped error (which corresponds to the harmonics generated) between of the calculated phases with respect to the ideal phase . Figures 9k to 9o present the same wrapped error maps with contour lines for visualization purposes. 6 (a) Phase LEF-HHT (b) Phase LEF-HHT -2 (c) Phase SLEF-HHT (d) Phase SLEF-LS (e) Phase SLEF-RE (f) Error LEF-HHT (g) Error LEF-HHT -2 (h) Error SLEF-HHT (i) Error SLEF-LS (j) Error SLEF-RE (k) Error LEF-HHT (contour) (l) Error LEF-HHT -2 (contour) (m) Error SLEF-HHT (contour) (n) Error SLEF-LS (contour) (o) Error SLEF-RE (contour) Figure 9: Phase estimation using LEF-GFB, SLEF-LS and SLEF-RE algorithms. (a), (b), (c) and (d) present the estimated phases. (e), (f), (g) and (h) are the wrapped error maps. (i), (j), (k) and (l) are the 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 cor- responds to the LEF–HHT algorithm b ut using (28) to estimate the phase map. It can be seen that the amount of errors is drasti- cally reduced due to the use of a di ff erent 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 al- gorithm 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 e xperiment, the error of the estimation of the step using the LEF–HHT algorithm was δ err or = 0 . 1194 r ad while the error for the SLEF–LS algorithm was δ err or = 0 . 0567 rad and the SLEF–RE was δ err or = 0 . 0219 r ad . The SLEF–RE method with HHT pre–filtering presented an error of δ err or = 0 . 1092 rad . The MAEs of the error surfaces are M A E LE F − H H T = 1 . 1712, M A E LE F − H H T − 2 = 0 . 4722 r ad , M AE S L E F − H H T = 0 . 4846 r ad , M A E S L E F − LS = 0 . 37 r ad and M A E S L E F − RE = 0 . 3569 r ad . It is important to note that the high amount of harmonics in Figure 9h 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 9l. 4.4. Experimental results In order to test the capabilities of our algorithm we imple- mented a Polarizing Cyclic Path Interferometer (PCPI) as the one proposed in [40] in its radial mode. Such system was im- plemented with a diode Laser with power of 400 mW operating at λ = 532 nm . In the arrangement, a polarizer filter is placed at an angle of 45 ◦ at the entrance of the interferometer , so the incoming beam will have perpendicular and a parallel compo- nents of the same intensity . Once the light has gone through the PCPI, we placed a quarter wa ve plate which generates crossed circular polarization states. At the output of the PCPI, an in- terference pattern is observed when placing an auxiliary linear polarizer , according the next equation I k ( p ) = a k ( p ) + b k ( p ) cos[ φ ( p ) + 2 Ψ k ] + η k ( p ) , (34) which is in fact the same as (1), b ut in this case, Ψ k corresponds to the angle of the linear polarizer . This rotation produces the phase shift between the two interferograms as δ k = 2 Ψ k . 7 (a) I 1 , δ = 0 (b) I 2 , δ = π/ 3 (c) Lissajous Pattern (d) ˆ I 1 , HHT pre–filtering (e) ˆ I 2 , HHT pre–filtering (f) Lissajous Pattern (g) ˆ I 1 , GFB pre–filtering (h) ˆ I 2 , GFB pre–filtering (i) Lissajous Pattern Figure 10: Experimental results from a PCSI and filtered patterns. (a) and (b) are the experimental patterns, (c) and (d) are the patterns filtered with the HHT , (e) and (f) are the patterns filtered with GFB For this experiment, the auxiliar linear polarizer was mounted on a graduated rotational mount. The two interferograms were recorded with a resolution of 410 × 410 and a shift angle of 30 ◦ between the captures, which means a phase step of δ = π/ 3. Figures 10a and 10b present such patterns. Figures 10d and 10e show the pre–filtered patterns using the HHT . Such result w as obtained using the EFEMD algorithm with its automatic selection of modes. On the other hand, Fig- ures 10g and 10h illustrate the normalized patterns using GFB. For these results applied a Winner –T ake–All (WT A) strate gy 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 ˆ I k ← α ˆ I k + (1 − α ) ˜ I k , (35) where α = M ag / max( M ag ), M ag is the magnitude map of the response of the GFB and ˜ I k is low–pass filtered version of the original image. T o pro ve the demonstrate the feasibility of our proposal, we present the phase obtained with the LEF–HHT , SLEF–LS and the SLEF–RE algorithms. Figure 11a presents the estimated phase map using the LEF– HHT algorithm, were the pre–filtering w as made with the HHT (Figures 10d and 10e). The estimation of the phase step by us- ing the 5–term equation and the LS method, which ga ve a result (a) φ LE F − H H T (b) φ S L E F − LS (c) φ S L E F − R E (d) φ 4 − ste p Figure 11: Phase maps estimated with (a) LEF–HHT , (b) SLEF–LS and (c) SLEF–RE (d) Four Steps of δ = 1 . 2186, a relativ e error of 16 . 4%. The calculation of the phase map was made with (7). Figure 11b presents the estimated phase map using the SLEF– LS algorithm, were the pre–filtering was made with the GFB (Figures 10g and 10h). The estimation of the phase step by us- ing the 2–term equation and the LS method (which under these conditions have the same beha vior as the 5–term), which gav e a result of δ = 1 . 1673, a relative error of 11 . 4%. The calculation of the phase map was made with (28). Figure 11c presents the estimated phase map using the SLEF– RE algorithm, were the pre–filtering was made with the GFB (Figures 10g and 10h). The estimation of the phase step by us- ing the 2–term equation and the robust estimator , which gave a result of δ = 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 qualitati ve comparison, we include the phase estimated with the well-known 4 steps algorithm, shown in Figure 11d. The fringe patterns of this result were only filtered with a simple Gaussian filter of σ = 2. From a visual inspection, one can note that SLEF–LS and SLEF–RE algorithms present practically the same results. Mean- while, 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 algo- rithm (which is more robust to noise than a two–steps method). 8 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. Ne vertheless, this algorithm re- quires at least that the background intensity as well as the ampli- tude modulation to be temporary constants [13, 23 – 26]. If this condition is not presented, it is required a preparation process for the interferograms in order to accomplish a constant mod- ulation of the fringes. Preprocessing techniques such as Win- dowed Fourier Transform [28], GFB [16], the Hilbert–Huang T ransform (HHT) [17] or isotropic normalization [30] allow to eliminate such variations. Our work presents the adv antage 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. W e compare our approach with theirs in order to prove the equiv alency of our 2–terms SLEF–LS algorithm with the 5–term pre–filtered tech- nique. 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 pro- cess. 6. Conclusions W e introduced a simplified model of the Lissajous Ellipse Fit- ting to calculate the phase step and phase distribution of two randomly shifted interferograms. W e focused on solving the problem of obtaining the phase of interferograms with spatial– temporal dependencies on their background intensities, ampli- tude modulations and noise. The main advantages of use of the GFB is that the phase estimation is robust to the mentioned is- sues since it filters–out the noise, normalizes the amplitude and eliminates the background. Giv en the normalized patterns, the ellipse equation can be simplified to a two–unkno wns system instead of a five–unkno wns 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 coe ffi cients of the ellipse’ s equation and the calculation the phase distribution. W e re- mark that we can replace the GFB based preprocessing with other normalization techniques that provide the elimination of the background component, normalize the amplitude modula- tion 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 im- prov e the accuracy of the phase step estimation, this is mainly caused by residuals of the filtering process. The e xperimental results of the calculation of our algorithms to 100 pairs of im- ages (10 di ff erent patterns with 5 di ff erent levels of noise and 5 di ff erent phase steps) prove that SLEF–LS is equiv alent to the pre–filtered 5–term LEF algorithm and the robustness to resid- uals of our SLEF–RE algorithm, which improves significantly the accuracy of the phase step estimation. Acknowledgements VHFM thanks Consejo Nacional de Ciencia y T ecnolog ´ ıa (Cona- cyt) for the provided postdoctoral grant. This research was sup- ported in part by Conacyt, Mexico (Grant A1-S-43858) and the NVIDIA Academic program. References References [1] D. Malacara, Optical shop testing, John Wile y & Sons, 2007. [2] M. Serv ´ ın, J. A. Quiroga, M. Padilla, Fringe pattern analysis for opti- cal metrology: theory , algorithms, and applications, John W iley & Sons, 2014. [3] K. Creath, T emporal phase measurement methods, Interferogram analysis 96. [4] M. T akeda, H. Ina, S. Kobayashi, Fourier-transform method of fringe- pattern analysis for computer -based topography and interferometry , JosA 72 (1) (1982) 156–160. [5] P . Carr ´ e, Installation et utilisation du comparateur photo ´ electrique et in- terf ´ erentiel du Bureau International des Poids et Mesures, Metrologia 2 (1) (1966) 13. [6] G. Rodriguez-Zurita, N.-I. T oto-Arellano, C. Meneses-Fabian, J. F . V ´ azquez-Castillo, One-shot phase-shifting interferometry: fiv e, sev en, and nine interferograms, Optics letters 33 (23) (2008) 2788–2790. [7] N.-I. T oto-Arellano, V . H. Flores-Mu ˜ noz, B. Lopez-Ortiz, Dynamic phase imaging of microscopic measurements using parallel interferograms gen- erated from a cyclic shear interferometer, Optics express 22 (17) (2014) 20185–20192. [8] V . Flores-Mu ˜ noz, N. T oto-Arellano, B. Lopez-Ortiz, A. M. Garc ´ ıa, G. Rodr ´ ıguez-Zurita, Measurement of red blood cell characteristic us- ing parallel phase shifting interferometry , Optik-International Journal for Light and Electron Optics 126 (24) (2015) 5307–5309. [9] B. T . Kimbrough, Pixelated mask spatial carrier phase shifting interfer- ometry algorithms and associated errors, Applied Optics 45 (19) (2006) 4554–4562. [10] X. Meng, L. Cai, X. Xu, X. Y ang, X. Shen, G. Dong, Y . W ang, T wo- step phase-shifting interferometry and its application in image encryption, Optics letters 31 (10) (2006) 1414–1416. [11] J. V argas, J. A. Quiroga, C. Sorzano, J. Estrada, J. Carazo, T wo-step inter - ferometry by a regularized optical flow algorithm, Optics letters 36 (17) (2011) 3485–3487. [12] T . M. Kreis, W . P . Jueptner, Fourier transform ev aluation of interference patterns: demodulation and sign ambiguity , in: Laser Interferometry IV : Computer-Aided Interferometry , V ol. 1553, International Society for Op- tics and Photonics, 1992, pp. 263–274. [13] C. Farrell, M. Player , Phase step measurement and v ariable step algo- rithms in phase-shifting interferometry , Measurement Science and T ech- nology 3 (10) (1992) 953. [14] J. V argas, J. A. Quiroga, C. Sorzano, J. Estrada, J. Carazo, T wo-step de- modulation based on the Gram–Schmidt orthonormalization method, Op- tics letters 37 (3) (2012) 443–445. [15] J. Deng, H. W ang, F . Zhang, D. Zhang, L. Zhong, X. Lu, T wo-step phase demodulation algorithm based on the extreme value of interference, Op- tics letters 37 (22) (2012) 4669–4671. [16] M. Rivera, O. Dalmau, A. Gonzalez, F . Hernandez-Lopez, T wo-step fringe pattern analysis with a Gabor filter bank, Optics and Lasers in En- gineering 85 (2016) 29–37. [17] M. T rusiak, K. Patorski, T wo-shot fringe pattern phase-amplitude demod- ulation using Gram-Schmidt orthonormalization with Hilbert-Huang pre- filtering, Optics express 23 (4) (2015) 4672–4690. [18] X. Meng, L. Cai, Y . W ang, X. Y ang, X. Xu, G. Dong, X. Shen, X. Cheng, W avefront reconstruction by two-step generalized phase-shifting interfer- ometry , Optics Communications 281 (23) (2008) 5701–5705. [19] M. Wielgus, Z. Sunderland, K. Patorski, T wo-frame tilt-shift error esti- mation and phase demodulation algorithm, Optics letters 40 (15) (2015) 3460–3463. 9 [20] R. K ulkarni, P . Rastogi, T wo-step phase demodulation algorithm based on quadratic phase parameter estimation using state space analysis, Optics and Lasers in Engineering 110 (2018) 41–46. [21] F . Liu, Y . Wu, F . Wu, Phase shifting interferometry from two normalized interferograms with random tilt phase-shift, Optics express 23 (15) (2015) 19932–19946. [22] D. Saide, M. T rusiak, K. Patorski, Ev aluation of adaptively enhanced two- shot fringe pattern phase and amplitude demodulation methods, Applied optics 56 (19) (2017) 5489–5500. [23] Y . Zhang, X. Tian, R. Liang, Random two-step phase shifting interfer- ometry based on Lissajous ellipse fitting and least squares technologies, Optics express 26 (12) (2018) 15059–15071. [24] Y . Zhang, X. Tian, R. Liang, T wo-step random phase retriev al approach based on Gram-Schmidt orthonormalization and lissajous ellipse fitting method, Optics express 27 (3) (2019) 2575–2588. [25] F . Liu, J. W ang, Y . W u, F . W u, M. T rusiak, K. Patorski, Y . W an, Q. Chen, X. Hou, Simultaneous extraction of phase and phase shift from two in- terferograms using Lissajous figure and ellipse fitting technology with Hilbert–Huang prefiltering, Journal of Optics 18 (10) (2016) 105604. [26] C. Meneses-Fabian, F . A. Lara-Cortes, Phase retriev al by Euclidean dis- tance in self-calibrating generalized phase-shifting interferometry of three steps, Optics express 23 (10) (2015) 13589–13604. [27] L. Muravsky , O. Ostash, T . V oronyak, I. Andreiko, et al., T wo-frame phase-shifting interferometry for retriev al of smooth surface and its dis- placements, Optics and Lasers in Engineering 49 (3) (2011) 305–312. [28] Q. Kemao, T wo-dimensional windowed Fourier transform for fringe pat- tern analysis: principles, applications and implementations, Optics and Lasers in Engineering 45 (2) (2007) 304–317. [29] M. Trusiak, M. W ielgus, K. Patorski, Advanced processing of optical fringe patterns by automated selectiv e reconstruction and enhanced fast empirical mode decomposition, Optics and Lasers in Engineering 52 (2014) 230–240. [30] J. A. Quiroga, M. Servin, Isotropic n-dimensional fringe pattern normal- ization, Optics communications 224 (4-6) (2003) 221–227. [31] J. G. Daugman, Uncertainty relation for resolution in space, spatial fre- quency , and orientation optimized by two-dimensional visual cortical fil- ters, JOSA A 2 (7) (1985) 1160–1169. [32] J. G. Daugman, Complete discrete 2-d gabor transforms by neural net- works for image analysis and compression, IEEE Transactions on acous- tics, speech, and signal processing 36 (7) (1988) 1169–1179. [33] J. G. Daugman, High confidence visual recognition of persons by a test of statistical independence, IEEE transactions on pattern analysis and ma- chine intelligence 15 (11) (1993) 1148–1161. [34] M. Rivera, Robust fringe pattern analysis method for transient phenom- ena, Optics and Lasers in Engineering 108 (2018) 19–27. [35] W . Jun, A. Asundi, Strain contouring with Gabor filters: filter bank de- sign, Applied optics 41 (34) (2002) 7229–7236. [36] P . J. Huber, Robust estimation of a location parameter , in: Breakthroughs in statistics, Springer , 1992, pp. 492–518. [37] M. J. Black, A. Rangarajan, On the unification of line processes, outlier rejection, and robust statistics with applications in early vision, Interna- tional Journal of Computer V ision 19 (1) (1996) 57–91. [38] P . Charbonnier, L. Blanc-F ´ eraud, G. Aubert, M. Barlaud, Deterministic edge-preserving regularization in computed imaging, IEEE Transactions on image processing 6 (2) (1997) 298–311. [39] S. M. Bhuiyan, R. R. Adhami, J. F . Khan, Fast and adapti ve bidimensional empirical mode decomposition using order-statistics filter based en velope estimation, EURASIP Journal on Advances in Signal Processing 2008 (1) (2008) 728356. [40] N.-I. T oto-Arellano, 4d measurements of biological and synthetic struc- tures using a dynamic interferometer , Journal of Modern optics 64 (sup4) (2017) S20–S29. 10