1. Introduction

Applications of optical interferometry in areas such as digital contouring or strain analysis in mechanics are very important because they are non-invasive techniques and gives full field measurements [1]. For example, the electronic speckle pattern interferometry technique (ESPI), is widely used in this areas [2]. As the interesting information comes phase modulated in the generated image interferogram, it is necessary to apply a demodulation technique.

A Single Fringe Pattern Image (SFPI) is an image interferogram that is considered as a sinusoidal two-dimensional signal (without taking random fluctuations into account) given by the following function:

assuming that a(x,y) and b(x,y) are continuous smooth functions. (x,y) are integer values that represent a site on the image lattice L where the interferogram is recorded. However, in this paper we assume that we are working with a normalized fringe pattern as the following:

There are several fringe pattern normalization techniques. For example, Quiroga et. al. mention an isotropic fringe pattern normalization, while Guerrero et. al use monogenic filtering for the fringe pattern normalization [3, 4].

The demodulation process is based on recovering the encoded phase ϕ(x,y) given by the fringe pattern image. To do this, the only information we have is the intensity values I(x,y) of the fringe pattern. When it is possible to introduce a carrier frequency while the fringe pattern image is being recorded, the modulating phase is given as:

where ω is the introduced spatial carrier frequency in x direction. In this case, the modulating phase is recovered by the well known Fourier transform methods [5]. However, when non repetitive transitory phenomena are being analyzed, based on speckle interferometry systems, it is commonly impossible to introduce a carrier frequency, be it spatial or temporal [6, 7]. In these cases, there is no choice but to deal with a SFPI with closed fringes.

Two methods that have marked the way to demodulate SFPI with closed fringes are the Regularized phase tracker (RPT) by Servin et. al. [8] and the two-dimensional Hilbert Transform method (2D-HT) by Larkin et. al. [9]. 2D-HT defines the next vortex operator to obtain the interferogram’s quadrature in order to demodulate it:

where the modulating phase is obtained as follow:

The term β _2π=β _2π (x,y) in Eq. (4) is the modulo 2π fringe pattern’s orientation field and may be defined as follow:

Function exp[iψ(u,v)] is an spiral phase function given by the following expression:

As we can see in Eq. (4), we need to have the fringe orientation field β _2π(x,y) in order to use the vortex operator (4). Eq. (6) looks very easy to compute, but we obviously do not have access to the phase of the interferogram. Instead, we have access to the interferogram’s image and we have to estimate the orientation field β _2π(x,y) modulo 2π given the interferogram’s image. However, to estimate the fringe orientation field β _2π(x,y) from the image interferogram is almost as difficult as to estimate the phase using RPT method. Larkin et. al. suggested this quadrature operator in Ref. [9], but in their paper it is not presented how to obtain the fringe orientation field. However, in Ref. [10], an outline of the orientation estimation problem was given by Larkin and in [11] propose a way to obtain the orientation field.

Quiroga et. al. [12] and Villa et. al. [13] proposed separate sequential methods for the fringe orientation field estimation. For example, in Ref. [13] a linear orientation field estimator for β _2π(x,y) is proposed that works minimizing sequentially along the fringes the next cost functional with respect to n(x,y)=[n_x(x,y),n_y(x,y)] at each site (x,y):

where n ₀(η,ξ) is the vector field already estimated in neighborhood Γx,y around the site (x,y) and operator ∇(·) is the gradient operator. Once given the normalized vector field, the orientation field is obtained as:

On the other hand, RPT which has inspired all the fringe following sequential methods, obtains the interferogram’s phase by minimizing the following cost functional with respect to $\widehat{\varphi }$(x,y), ω_x and ω_y at the site (x,y):

where P(x,y,η,ξ)=$\widehat{\varphi }$(x,y)+ω_x(x-η)+ω_y(y-ξ) is a phase plane that fits the observed data in neighborhood Γ_x,y around the site (x,y), where I′(x,y) is the normalized observed data given by (2). In both cases (Eq. 8 and 10) the estimation process starts at a site (x,y) which is the initial seed that is sequentially propagated along a path on the image interferogram.

 

Fig. 1. Two estimated phases using the RPT from the same SFPI given. Panel labeled Wrong phase shows the obtained phase by following an arbitrary path, while panel labeled Expected phase is obtained using a fringe following path. The SFPI given is in panel labeled SFPI.

Download Full Size | PDF

The minimization process for Eq. (8) and (10) is risky and depends strongly on the initial seed given and on the path followed during the SFPI demodulation. These methods are called path dependent since from all the possible paths to follow in the two-dimensional interferogram’s space, only the path following the fringes yields the expected phase. To illustrate this, in Fig. 1, we see two global solutions given for the same SFPI following different paths with the RPT. In the panel labeled Wrong phase, we see a phase estimation following an arbitrary path, while in the panel labeled Expected phase, we see a phase estimation following the fringes, making the path dependence of the RPT evident.

In this article, we propose a novel idea, which constrains the phase estimation to look into the function space C ²; that is, estimated phases with a continuous phase curvature. With this restriction we remove the local ambiguities avoiding wrong solutions as that shown in Fig. 1 in the panel labeled Wrong phase. Searching into this function space C ², we obtain the expected modulated phase regardless of the path followed by the demodulator. Our technique consists in developing a local frequency estimator with second order regularizer potentials to obtain the expected phase $\widehat{\varphi }$∈C ². This frequency estimator searches the expected frequency into a discrete frequency set. In the next section, we will develop the local frequency estimator and we will show how to obtain the phase once given the estimated local frequency for each site in the image interferogram. Also, we are going to show with examples, that the technique developed is path independent and obtains the expected modulating phase.

2. Method

We are going to assume a one-dimensional space x. Then our normalized fringe pattern looks like:

To avoid ambiguities in the demodulation process, we limit the estimator’s degrees of freedom by restricting the demodulated phase ϕ withing the set C ² of twice differentiable functions. With this in mind, we start proposing the phase estimation by adjusting the local frequency. Then, we have the following relation:

where $\widehat{\varphi }$(x) is previously estimated, ω(x) is the local frequency between x and x ₊ and I′(x) is the normalized signal (Eq. 11). Site x ₊ is such that may take the following values:

depending on the demodulation direction taken. Thus, we propose to minimize the following cost function to estimate the frequency ω(x) that satisfies Eq. (12):

In Fig. 2(a) we show the graph for this cost function. As shown in Fig. 2(a), this cost function is ambiguous since it has multiple minimum values due to the cosine periodicity. But, all the minimum values of Eq. (14) can be generated with the following two series:

for n=0,1,2, …, where the wrapping operator W[·] wraps the phase $\widehat{\varphi }$(x) in the interval [-π,π]. These series are also showed in graph 2(a) to see its correspondence with the minimum values.

As the maximum frequency that a fringe pattern may have in a image interferogram is π radians per pixel, we are interested on the frequencies that are in the interval [-π,π]. These frequencies may be taken from the series of Eq. (15) and (16) as follow:

for ω(x)=(W[$\widehat{\varphi }$(x)] — arccos[I′(x ₊)])(x-x ₊) and

for ω(x)=(W[$\widehat{\varphi }$(x)]+arccos[I′(x ₊)])(x-x ₊). Then we define the discrete set Ω={ω ₁,ω ₂} as the domain of all possible frequencies that minimize (14) and are in the interval [-π,π].

 

Fig. 2. In this figure, (a) is the graph of the cost function (14) and (b) is the graph of the regularizer potential (20). The graph (a) was generated using the parameters I′(x)=0.66 and $\widehat{\varphi }$(x)=20. We also show the series S ₁ and S ₂, which corresponds to the minimum values of the cost function (14). In graph (b) we see that exists only one frequency ω(x)∈Ω={S ₁,S ₂} that minimizes Eq. (20).

Download Full Size | PDF

Because we want that the estimated phase $\widehat{\varphi }$(x) belongs to the function space C ², we take the priori assumption that the phase curvature in a site x is very smooth in its neighborhood. Then, we suggest the following potential regularizer to penalize the strong phase curvature variations between the neighbors x and x+1:

where, as we can see in this equation, the operator C _{$\widehat{\varphi }$}(·) estimates the backwards phase curvature in the given site. This, in frequency terms is reduced to:

if we take the demodulation direction to the right or

if we take the demodulation direction to the left. As we can see in Fig. 2(b), in the interval [-π,π], this potential regularizer removes the ambiguities given by cost function (14) since it has only a frequency in Ω that minimizes it. Therefore, we suggest the following estimator to obtain the frequency:

where the searched frequency must to be in the discrete set Ω={ω ₁,ω ₂}. V[ω(x)] is defined as in Eq. (20) or (21), depending on how the demodulation direction is taken. Finally we obtain the phase for the current site x ₊ given by the frequency $\widehat{\omega }$(x) with the estimator of Eq. (22) as follow:

which is continuous in curvature since the frequency $\widehat{\omega }$(x) obtained with estimator (22), gives the minimum possible phase curvature variation between the neighbors sites of x and x ₊.

Specialist readers may note a little similar Eq. (23) with the Phase Locked Loop (PLL) method applied to demodulate interferograms. In 1993, Servin et. al. showed the implementation of a PLL method to demodulate SFPI with carrier frequency [14]. Its iteration to update the phase in a site x+1 is given as follow:

where, as we can see, it must exist a carrier frequency ω ₀ in order to apply PLL method successfully. In literature, we can find second-order PLL implementations applied to fringe patterns as that reported in [15], but all PLL methods works only if a carrier frequency is given in the fringe pattern image. Therefore, there are strong differences between (24) and (23). First, the frequency for (23) is estimated using (22) in such a way that it preserves a continuous phase curvature and second, frequency estimator (22) is applicable to fringe pattern signals with and without carrier frequency. In contrast, PLL only assumes phase continuity.

2.1. Demodulation process.

The demodulation process, considering the previous development, may be summarized with the following steps:

1. Choose a site x as the initial seed to start the demodulation sequence.

2. Set the estimated phase in the initial seed as $\widehat{\varphi }$(x)=arccos[I′(x)].

3. Set the estimated phase in the initial seed neighborhood as $\widehat{\varphi }$(x+1)=arccos[I′(x+1)] and $\widehat{\varphi }$(x-1)=arccos[I′(x-1)].

4. For each subsequent site x ₊ to the left and right, that isn’t estimated, do the following:

(a) Construct the frequency discrete set Ω={ω ₁,ω ₂} using Eq. (17) and (18).

(b) Take the frequency $\widehat{\omega }$(x) as shown in Eq. (22).

(c) Set the estimated phase in the site x ₊ using Eq. (23).

In this way, we have a simple and fast demodulation process that we are going to call Frequency Curvature Tracker (FCT) because it estimates the frequency restricting its curvature to obtain the phase in each site x.

A more robust way to set the phase in the neighborhood of the initial seed in step 3, is the following: construct the frequency set for site x-1 and x+1 and choose the frequency combination that gives less phase curvature. In other words, we use the following criteria to choose the frequencies for sites x-1 and x+1:

The phase ϕ(x-1) and ϕ(x+1) is set according to Eq. (23), where Ω_l is the frequency discrete set to the left (site x-1) and Ω_r is the frequency discrete set to the right (site x+1). This is better than using arc cosine functions like in step 3 because the frequency around the initial seed with arc cosine functions, may abruptly change its sign.

 

Fig. 3. In graph 1) a signal is shown which phase is modeled as $\varphi \left(x\right)=\frac{1}{2}a{x}^2$. In graph 2) we show the obtained phase using the FCT and in graph 4) we show its curvature. In graph 3) we show the obtained phase using the RPT technique and in 5) we show its curvature. As we can see, the phase curvature presented in 4) is a continuous line while the phase curvature in 5) has one abrupt variation.

Download Full Size | PDF

Now, we are going to test the FCT and to campare it with the RPT applied to the onedimensional case to observe that the obtained phase $\widehat{\varphi }$ with the FCT is $\widehat{\varphi }$∈C ². Let us demodulate the signal showed in 3(a) which its modulating phase is modeled as $\varphi \left(x\right)=\frac{1}{2}a{x}^2$. In Fig. 3(b) we see the phase obtained with the FCT technique and Fig. 3(d) shows its phase curvature. In Fig. 3(c) we see the phase obtained with the RPT technique and Fig. 3(e) shows its phase curvature. As we can see, the phase curvature obtained with the FCT is a continuous line as expected. On the other hand, we can see that the phase curvature obtained with the RPT technique is a line with one abrupt change, therefore, the phase obtained $\widehat{\varphi }$ with the RPT is $\widehat{\varphi }$∉C ² while the FCT obtains a phase $\widehat{\varphi }$∈C ². Both phase profiles in Fig. 3(b) and 3(c) are compatible with the one-dimensional fringe pattern signal given in Fig. 3(a). However, estimated phase profile $\widehat{\varphi }$∈C ² shown in Fig. 3(b) is typical phase profile of a fringe pattern with closed fringes. As our final intention is to demodulate SFPI with closed fringes, we restrict phase solutions to those phase profiles which are $\widehat{\varphi }$∈C ².

Summing up, until here we have developed a demodulation method with a frequency estimator (Eq. 22) that searches in a discrete frequency set. It was tested for one-dimensional signals without carrier as is shown in Fig. 3. Now, we are going to apply the method in two-dimensional signals.

3. Two-dimensional case

In this case we are going to implement the frequency curvature tracker or FCT method, to demodulate square frames into the image interferogram lattice around a given site as a starting seed. To do this, we are going to use the FCT (see subsection 2.1) using a row by row scanning strategy to demodulate sequentially all sites in the frame. Each row and column is going to be treated as a one-dimensional signal.

 

Fig. 4. Steps for the row by row scanning strategy to demodulate a frame in the image interferogram lattice. Step 1 shows the initial seed with a circle, Steep 2 shows the column demodulation, Step 3 shows the row demodulation and Step 4 the row by row scanning to demodulate each row using the values of the column as an initial seed for each row.

Download Full Size | PDF

Let (η,ξ) be the initial seed in the frame, where η is the column number and ξ is the row number. Then we set the phase of the initial seed as $\widehat{\varphi }$(η,ξ)=arccos[I′(η,ξ)], where I′(·) is defined as in Eq. (2). For a graph illustration see Fig. 4(Step 1), where the initial seed is shown. The next step is to demodulate the column of the initial seed using the FCT described in subsection 2.1, taking the column as a one dimensional signal as shown in Fig. 4(Step 2); the same is done for the row of the initial seed as shown in Fig. 4(Step 3). Finally, we scan row by row and each row is demodulated with the FCT, using as an initial seed the value of the column as shown in Fig. 4(Step 4). Thus, we have a row by row demodulation process to demodulate a square frame in a image interferogram.

Our strategy to demodulate an image interferogram following an arbitrary path is by following the path with square frames as is shown in Fig. 5. Therefore, we explained before how to demodulate a frame. As we can see in Fig. 5, a frame in a site slope the frame of the neighbor site. Then, to maintain a spatial coherence between frames, only the seed of the first demodulated frame is set using the arc cosine function. The seed of the subsequent frame is already estimated in the previous frame.

4. Tests and results

Now let us to test FCT method with a simulated SFPI with closed fringes following an arbitrary path. We are going to compare this result with the result of RPT method following an arbitrary path. For this, we propose a free hand drawn path. In Fig. 6(a) we see the generated SFPI I ′(x,y) while in Fig. 6(b) the drawn path to follow. We can see the phase obtained using the FCT in Fig. 6(c) while in Fig. 6(d) the phase obtained with the RPT following the proposed path. As we can see, the RPT fails to obtain the expected modulated phase. On the other hand, we can see that the FCT obtains the expected modulated phase even if we are not following the fringes. In this sense, the FCT is more robust than the RPT since it obtains the expected phase without taking care on following fringes.

 

Fig. 5. In this figure we illustrate a path that is followed with square frames. Each site of the path, represted with a dark point, is the center of its frame and each frame is such that intersects the frame of the neighbor site and the neighbor site belong to this intersection.

Download Full Size | PDF

 

Fig. 6. (a) is the simulated image interferogram I′(x,y), (b) is the path to follow, (c) is the demodulated phase using the FCT in two dimensions and (d) is the demodulated phase using the RPT.

Download Full Size | PDF

Figure 7, shows two even more complicated fringe structures to demonstrate the FCT method’s ability to demodulate complicated fringes without following the fringes. The actual phase of the fringe patterns shown in Figs. 7(a) and 7(b), were generated using a linear combination of Gaussian functions.

Figure 8, shows an experimental SFPI generated with ESPI technique for contouring. As the FCT method assume that the modulating phase of the SFPI is continuous in curvature, we can not apply directly the FCT method to the image interferogram shown in Fig. 8(a). Firstly, we must use a normalization method able to normalize and filter the fringe pattern. In this case we have used the normalization method reported in Ref. [4] tuning properly the filters. Figure 8(b) shows the obtained normalized fringe pattern. The estimated phase is shown in Fig. 8(c) while Fig. 8(d) shows the estimated phase wrapped for illustration purposes.

4.1. Computational time cost report

We have tested the FCT method for various interferogram sizes and demodulated successfully. The computational time cost for the different interferogram image sizes is as follow: size of 1024×1024, 1.14 seconds; size of 512×512, 0.320 seconds; size of 256×256, 0.089 seconds; size of 128×128, 0.016 seconds. These tests were made on a personal computer with a 64bit processor at 2GHz and 1GB of memory RAM using double precision data. The FCT algorithm was programmed using a “C” compiler.

5. Discussion and Conclusions

With this paper we have demonstrated that estimating the phase into a function space C ², removes the local ambiguities given by the actual phase estimators as the RPT (see Fig. 2). To estimate the phase into the function space C ², we have developed an estimator based on the local frequency estimation within a discrete frequency set. The main advantage of this estimator is that we can to demodulate a SFPI with closed fringes regardless of the path to follow (see Fig. 6). Another advantage is that we do not need to use a minimization process like the steepest-descent algorithm (see Ref. [16]) for non linear systems to estimate the phase as the RPT method. As we can see in subsection 4.1, this method is a fast demodulation process. Finally, we can use a simple row by row scanning strategy like the one shown in Fig. 4 to demodulate the image interferogram. This is important because we do not need program complicated fringe following algorithms.

On the other hand, proposed method assumes that the modulating phase is continuous and that fringe pattern is normalized. Then, we can not apply this method directly to an image interferogram like that shown in Fig. 8. It is necessary to preprocess the SFPI to normalize and filter the fringes. This is a drawback of the FCT method as is, because it can not be applied to a noise contaminated experimental fringe pattern directly. However, this drawback does not diminish the importance of the suggested method and its ideas (estimating the phase $\widehat{\varphi }$ as $\widehat{\varphi }$∈C ²). This is because phase profiles which are in C ², are typical phase profiles of single fringe pattern images with closed fringes. Furthermore, nowadays there is not one publishedwork with a path independent technique to demodulate single image interferograms with closed fringes. At present, we are working with some ideas to apply these concepts within a modified robust FCT method to be applied directly to noise contaminated fringe patterns.

It should be clearly remarked, that as shown by the results of Figs. 6, 7 and 8, the method obtains the expected phase regardless of the path followed by the FCT. If the underlying modulating phase is continuous and at least one time differentiable, then the method is guaranteed to be path independent and able to recover the phase.

 

Fig. 7. A more complicated structure of fringes. (a) and (b) are simulated fringe patterns and (c) and (d) shows its estimated phase with the FCT using a row by row scanning strategy. The phase is shown wrapped for illustration purposes.

Download Full Size | PDF

 

Fig. 8. (a) is a real experimental SFPI, (b) is its normalized version, (c) is the obtained phase with the FCT using a row by row scanning strategy and (d) is the wrapped phase showed for illustration purposes.

Download Full Size | PDF