1. Introduction

Phase shifting (PS) is a popular interferometric analysis technique that consists of to acquire a set of fringe patterns. In this work we proposed an algorithm for computing, simultaneously, a filtered phase map and the phase shifts from a noisy fringe patterns set with irregular and unknown phase steps. First, we introduce our notation and the fringe set model: Let I = {I ₁,I ₂,… ,I_K } a set of K phase shifted fringe patterns, then the k^th fringe pattern, I_k , is modelled by

where r = [x,y]^T denotes a pixel position in the image lattice L, a_k is the background illumination component, b_k is the fringe contrast, f is the unknown phase, Δ = {Δ₁, Δ₂,…, Δ_K} is the phase shift vector and η_k represents additive independent noise. We assume that a normalization process is performed on the fringe pattern so that the illumination components are for filtering out: thus we have that a_k ≈ 0 and b_k is estimated. Therefore the normalized fringe pattern, Î, is approximated by

where it is preferable to left b in the model because we do not want to perform division by zero b values (or noise). However, without lost of generality, in the rest of the paper we assume b ≈ 1. and drop it from the rest of the derivation.

$\widehat{\eta }$ is a remainder residual of the background illumination and the additive noise with low frequency bandwidth and unknown distribution. Such a kind of normalization can be achieved by means of a procedure as the reported in [1].

Recently there has been a thorough research on methods for recovering phase from a single closed fringe pattern (see [2] and references therein). However PS methods are preferably used in stable acquisition conditions (i.e. the temporal dependency of the illumination components, a and b, is eliminated) and the phase shifts, Δ, can be introduced with high accuracy. In such a case the wrapped phase f̂, can be recovered with very simple algorithms; where f̂ = W(f) ≝ f + 2πn for an integer n such that f̂ ∈(-π,π] (where W denotes the wrapping operator). For instance, for a low level noise, f̂ can be computed from 4 fringe patterns (by assuming Δ_k = kπ/2) with:

In order to compute the real phase, it is needed to unwrap the calculated wrapped phase. Here we denote by f = W ^-1(f̂) the unwrapped phase computed with the unwrapping operator W ^-1. However the phase unwrapping process is an ill posed problem because a wrapped phase may correspond to multiple unwrapped phases, i.e. there exist many f̂ ≠ f such that f̂ = W(f̂). For this reason, the unwrapping operator W ^-1 is, in general, implemented as the minimization of a regularized cost function, see for instance [4].

In this paper we address the problem of computing the phase maps from noisy, arbitrary shifted fringe patterns. Before to introduce our approach, we present in section 2 a brief review of standard procedures that consists of estimating the phase shifts follows by the use of a generalized phase shifting method (based on least square procedures) for computing the phase. In section 3, we proposed a robust method for the simultaneous computation of the phase map and the phase shifts. Our method if formulated as the minimization of a half-quadratic cost function that reduces the outliers contribution to the final results. A set of experiments, designed to demonstrate our method performances, are presented in section 4. Our conclusion are finally given in section 5.

2. Out of step phase shifting

In this section we presents a brief review of methods for out of step phase shifting, i.e. methods foe computing the phase from fringe patterns with arbitrary phase shifts. In this case, the phase can be computed with the Bruning-Grievenkamp algorithm by using arbitrary (but known) phase shifts [5][6] (section 14.8.2 in [7]):

where

see also [8]. If phase shifts are unknown then one can use a method for estimating them in a previous stage. Here we denote by δ the computed phase shifts vector. The development of methods for computing arbitrary phase shifts has been a prolific research [9]–[22]. An instance of such methods is the Brug’s method. This method computes the relative phase difference, ${\delta }_{\mathit{\text{IJ}}}^b$ , between a pair of fringe patterns I and J by means of the arccos of their correlation [12], i.e.

where .* denotes the componentwise product of vectors and ⟨∙⟩ denotes the mean intensity value of the whole interferogram. More recently Cai et al. proposed:

for computing the relative phase shift between the interferograms I and J [18]; where b is the illumination contrast component (assumed equal for both interferograms). An alternative method was reported by Larkin in [17]. Therein, Larkin noted that the shifts estimated by his method (but the same is valid for (6), (7) or a similar method) should be refined by using the phase computed with the generalized phase step method (4). That leads one to an iterative procedure that is iterated until an accurate result is computed. This iterative procedure has the disadvantages of computing from scratch at each iteration, both: The phase and the phase shifts. Moreover the problem of adjusting, simultaneously, phase and shifts is ill-posed. To understood this claim, one can consider the case of 3 fringe patterns with unknown phase shifts, in such a case one have more unknowns (illumination components, phase and shifts) than data (that are, besides, corrupted with noise).

Last ill-posed problem could be solved by a minimum last squares procedure that implies to increase the number of computed fringe patterns with the subsequent increment in the complexity of the experimental procedure; or by using regularization techniques that implies to incorporate prior knowledge about the solution.

The estimation by least squares is equivalent to the maximum likelihood estimation when additive Gaussian residuals are assumed, $\widehat{\eta }$ [or η if one assumes model (1)]. As it is well known, the accuracy of least–squares based methods is reduced by small signal to noise rations (snr) and by non-gaussian residuals. A source for deviations from Gaussian residuals is a non-sinusoidal pure fringe profile. A significant deviation of the residual from the Gaussianity produces phase artifacts as fringe harmonics [23]–[26]. Therefore, as it is noted in [7], the best performance is achieved with uniformly spaced, into the interval (-π, π], phase shifts.

Other authors that have reported methods for the joint estimation of phase and shifts are Okada et al. [27] and Marroquin et al. [28]. In particular, in [27] three least-squares fitting procedures are proposed for computing an initial phase, the phase map and the phase steps, respectively. In Ref. [28] there was proposed a method for computing unknown phase steps, therein is shown that the steps are easily computed, with a closed formula, if the quadrature fringe pattern set is known. Thus, the method in [28] proposed a complexus nonlinear optimization procedure for, simultaneously, computing: the quadrature fringe pattern set, the local frequency and the corresponding phase shifts.

3. Adaptive out of step phase shifting

In this section we present a method that unifies in a single procedure the two main ingredients in out of step phase shifting: The computation of the phase shifts and the phase field. Our method is based on the minimization of a regularized cost function and guarantees convergence.

In particular, we present a half-quadratic cost function that uses as data the fringe pattern and, guesses for the phase shifts and the phase map. Then by assuming that such initial guesses are close enough to the correct values, we formulate as unknowns the corrections for the phase shifts and phase field. Half-quadratic regularization is a framework for formulating robust (in the real statistical sense of robust estimators) cost function that can be minimized by iterating solutions of linear systems and/or closed formulas. In such a formulation an outlier rejection mechanism is incorporated in order to weight the contribution of the data to the final estimation [2][4][29]–[33].

First, we consider that a good guess δ = {δ ₁,δ ₂,…,δ_K } of the unknown phase steps, Δ is available. Therefore, we have

where α = {α ₁,α ₂,…,α_K } are the unknown phase step residuals. If we expect small residual phases, then a coarse estimation of the real phase can be computed by means of the minimization of (4) and an unwrapping method [4]. Thus, as in the phase shift case, we have the follows relationship between the real phase, f, and its coarse estimation, ψ:

where ϕ is an unknown residual phase field that corrects artifacts produced by noise and the residual steps, α. Thus by substituting (8) and (9) in the normalized model (2), we obtain

Remember that we are assuming b_k ≈ 1 for simplifying the presentation, but the fringe contrast term can and should be used in the procedure.

The task here is to compute such a residual phase, ϕ_r , and therefore to compute the real phase, f. For such a purpose we need to estimate the phase steps residuals, α.

Recently was reported an efficient method for computing the residual phase from a coarse one for a single fringe pattern: Algorithm 1 in [2]. Such method transforms an, originally, no-linear optimization problem in a sequence of quadratic optimization problems. Here we extend such a formulation for computing the real phase, ψ + ϕ, and the true phase steps, δ + α, from a phase shifted pattern sets. Following [2], we assume that |ϕ_r + α_k | is relatively small such that the first order Taylor series can be used to define the residual error:

Note that this residual error can also be obtained from (10) by using the trigonometric identity cos(x + y) = cosxcosy - sinxsiny [as in the Grievenkamp’s residual error, Eq. (5)] and then by assuming y small enough such that cosy ≈ 1 and siny ≈ y; thus cos(x + y) ≈ cosx-ysinx. Precisely, the advantage of formulating our problem as the computation of the small unknown corrections (instead of the whole value) is that the residual error can be approximated by the linear Taylor expansion E_kr . It is precautious to expect that there could be some places (pixels) where, actually, the residual error were large. Large residual error can be the result of a corrupted initial phase (for instance, phase corrupted with fringe harmonics), an unprecise initial phase shifts or data corruption (for instance, occlusions in the illumination). Then we incorporate a mechanism for weighting the contribution of the data in the estimation process. Such outliers rejection mechanism is implemented by a half-quadratic potential controlled by a parameter that avoids an over–detection or under-detection of atypical data.

As was discussed in section 2, if no additional data are used for solving the ill posed problem of, then prior knowledge should be. We assume that unknown real phase f is smooth in the sense of have small second derivatives. Thus we penalize the correction phase ϕ to produces a smooth f. Therefore, we propose to compute the phase correction field, ϕ, the phase step corrections, α, and an outliers detection field, ω, by an alternating quadratic minimization of the cost function:

where ⟨q, r, s⟩ is a triple of pixels (in the lattice, L) in horizontal, vertical or diagonal position such that (q,r) and (r,s) are first neighbors, see [2][32]. Last term is named the thin plate potential. That is related with the potential energy of a thin plate physical model with straight rest position, see [32] and references therein.

The cost function (12), that extends the previously proposed in [2], for dealing with a set of fringe pattern as data. Moreover, we have included two terms (weighted by the parameter γ) that enforce small values for the step correction vector, α, and the residual phase field, ϕ. Such terms improved significantly the convergence and stability of the minimization process. The details of the phase refinement procedure are formalized in Algorithm 1. It is important to note in Algorithm 1 that once a residual (ϕ or α) is computed, then the corresponding base variable (ψ or δ) is updated. Such a strategy reduces iteratively the value of the unknown residual and the fitness of the first order Taylor approach, and consequently the algorithm performance.

Algorithm 1 Adaptive out of step phase shifts.

Let g = {g ₁,g ₂,… ,g_K } a fringe pattern set with expected phase steps equal to δ and ψ an initial coarse phase.

Now we discuss how to compute an effective initial guess for the Algorithm 1. The initial coarse phase, ψ, can be computed with standard algorithms by assuming correct phase steps, i.e. by neglecting the residual steps, α. Such a wrapped phase is corrupted with artifacts introduced by the residual steps. Then it recommended to use a robust algorithm for unwrapping the coarse wrapped phase. In particular, we use the half-quadratic convex unwrapping algorithm reported in [4]. The resultant unwrapped phase, ψ, may have a constant residual step, δ_dc , that can be approximated by a reduced search, i.e:

where D = {d_i = 2πi/N}, for i = 0,1,2,…,N - 1, is a N steps set (we use N = 20 in our experiments). Alternatively, the method in Eqs. (6) or (7) can be used for the same purpose. Taking into account that the accuracy of these methods are reduced if the fringe patterns are corrupted by additive independent noise, see experiments in Fig. 4. Nevertheless how the δ_dc is estimated, we initialize the coarse phase by assigning: ψ ← ψ + δ_dc . On the other hand, the steps δ_k ’s can be initialized with the ideal values or, to reduce the risk of large α residuals, can be estimated (as δ_dc ).

 

Fig. 1. First row: Synthetic fringe pattern set. Second row: Reconstruction using the computed phase map and the computed phase steps.

Download Full Size | PDF

 

Fig. 2. Coarse solution computed with Eq. (3) by assuming ideal steps (π/2). From left to right: Wrapped phase, unwrapped phase computed with the convex algorithm in [4] and its cosine (reconstructed fringe pattern).

Download Full Size | PDF

In the following we present details of the partial minimizations (steps 5,7 and 8) in Algorithm 1. First, we note that, for ϕ = 0 and ω fixed, cost function (12) can be written as:

where ψ̂_kr ≝ ψ_r + δ_k and the potential Q(∙) contains the α–independent terms. Equating to zero the partial gradient with respect to (w.r.t.) α and solving for α_k , we obtain a closed formula for computing the α’s optimum coefficients:

In a similar way, for ϕ = 0 and α = 0, we obtain a closed formula for computing the ω field:

Note that the weight is close to one for the sites (pixels) where the model fits well the data and is close to zero for a large error fit. Finally, ϕ is obtained by solving the linear system that results of equaling to zero the partial gradient w.r.t. ϕ of U(ϕ, α = 0, ω), keeping α = 0 and keeping fixed ω. In particular, we use a Gauss-Seidel scheme similar to the one proposed in the Appendix A of [2].

 

Fig. 3. From left to right: Rewrapped phase (for illustration purposes), unwrapped phase and its cosine: Results computed with the proposed method (first row) and by smoothing, with a thin plate regularization filter [32], the coarsely computed phase in Fig. 2

Download Full Size | PDF

4. Experiments

First experiment demonstrates the method performance in synthetic noisy test data. First row in Fig. 1 shows the noisy fringe pattern set generated from a synthetic phase that has low and high frequencies: a slight tilt with sharp Gaussian peaks. Then a coarse phase map in Fig. 2 is computed with (3) by assuming regular phase steps equal to π/2 [Fig 4(a)]. The real random phase shifts are plotted in Fig. 4(b). As one can note in the Fig. 2 phase artifacts, correlated with the fringe pattern, corrupt the computed phase. Such a coarse phase is used as initial guess for the proposed method. The results are shown in the first row of Fig. 3. The proposed method recovers the wide bandwidth phase by smoothing spurious artifacts and preserving real high frequencies. Moreover the method recovers effectively the phase steps [see the phaser plot in Fig. 4(c)]. We performed experiments for different size (cardinality) of the fringe pattern set, K. As it is expected, the method performance is improved as K grows given that cost function (12) effectively incorporates redundant information in the fringe pattern set. On the other hand, a simple low–pass filtering computed by the minimization a thin plate potential,

despite redundant information. Thus, the spurious artifacts are not eliminated and real high frequencies are over–smoothed, see second row in Fig. 3. Fig. 4. shows obtained results from an electronic speckle pattern interferometry (ESPI) set. The fringe pattern corresponds to a steel plate under mechanical stress.

 

Fig. 4. Phaser plots: (a) Ideal phase shifts (with steps equal to π/2), (b) real phase shifts, (c) phase shifts computed with the proposed method and (d) phase shift computed with the method in [18].

Download Full Size | PDF

5. Conclusion

We have presented a phase shifting robust method for irregular unknown phase steps. The method is based on a half-quadratic regularized phase refinement strategy. The method takes advantage of the redundant information in the phase shifted fringe pattern set for smoothing out artifacts produced by miss–calibrated phase steps and for preserving real high frequencies. The method, implemented as successive quadratic minimizations of a nonlinear cost function, guarantees convergence to a local minimum and is computationally efficient. The method performance was demonstrated by experiments.

The authors were partially supported by CONACYT, Mexico: M. Rivera (grants 40722 and 46270) and R. Bizuet (Scholarship).

 

Fig. 5. Real data experiment: (a) An original ESPI fringe pattern, (b) computed wrapped phase with a standard four steps (assuming phase steps equal to π/2), (c) refined phase computed with the proposed algorithm (rewrapped for illustration purposes) and (d) computed phase shifts.

Download Full Size | PDF