Path independent demodulation method for single image interferograms with closed fringes within the function space C 2

In the last few years, works have been published about demodulating Single Fringe Pattern Images (SFPI) with closed fringes. The two best known methods are the regularized phase tracker (RPT), and the two-dimensional Hilbert Transform method (2D-HT). In both cases, the demodulation success depends strongly on the path followed to obtain the expected estimation. Therefore, both RPT and 2D-HT are path dependent methods. In this paper, we show a novel method to demodulate SFPI with closed fringes which follow arbitrary sequential paths. Through the work presented here, we introduce a new technique to demodulate SFPI with estimations within the function space C2; in other words, estimations where the phase curvature is continuous. The technique developed here, uses a frequency estimator which searches into a frequency discrete set. It uses a second order potential regularizer to force the demodulation system to look into the function space C2. The obtained estimator is a fast demodulator system for normalized SFPI with closed fringes. Some tests to demodulate SFPI with closed fringes using this technique following arbitrary paths are presented. The results are compared to those from RPT technique. Finally, an experimental normalized interferogram is demodulated with the herein suggested technique. © 2006 Optical Society of America OCIS codes: (120.2650) Fringe analysis; (120.5050) Phase measurement; (120.3940) Metrology. References and links 1. A. Mujeeb, V. U. Nayar and V. R. Ravindran, “Electronic Speckle Pattern Interferometry techniques for nondestructive evaluation: a review,” INSIGHT 45 275–281 (2006). 2. J. N. Butters and J. A. Leendertz, “A double exposure technique for speckle pattern interferometry,” J. Phys. E: Sci. Instrum. 4 277–279 (1971). 3. J. A. Quiroga and M. Servin, “Isotropic n-dimensional fringe pattern normalization,” Opt Commun. 224 221– 227 (2003). #74543 $15.00 USD Received 29 August 2006; revised 28 September 2006; accepted 29 September 2006 (C) 2006 OSA 16 October 2006 / Vol. 14, No. 21 / OPTICS EXPRESS 9687 4. J. A. Guerrero, J. L. Marroquin, and M. Rivera, “Adaptive monogenic filtering and normalization of ESPI fringe patterns,” Opt. Lett. 30 318–320 (2005). 5. M. Takeda, H. Ina, and S. Kobayashi, “Fourier transform methods of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72 156–160 (1982). 6. K. A. Stetson and W. R. Brohinsky, “Fringe-shifting technique for numerical analysis of time-average holograms of vibrating objects,” J. Opt. Soc. Am. A 5 1472 (1988). 7. J. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfel, A. D. White, and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surface and lenses,” Appl. Opt. 13 2693–2703 (1974). 8. M. Servin, J. L. Marroquin and F. J. Cuevas, “Fringe-follower regularized phase tracker for demodulation of closed-fringe interferograms,” J. Opt. Soc. Am. A 18 689–695 (2001). 9. K. G. Larkin, Donald J. Bone, and Michael A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. i. general background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18 1862-1870 (2001). 10. K. G. Larkin, “Natural demodulation of two-dimensional fringe patterns. II. Stationary phase analysis of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18 1871-1881 (2001). 11. K. Larkin, “Uniform estimation of orientation using local and nonlocal 2-D energy operators,” Opt. Express 13 8097-8121 (2005). 12. J. A. Quiroga, Manuel Servin and Francisco Cuevas, “Modulo 2 fringe orientation angle estimation by phase unwrapping with a regularized phase tracking algorithm,” J. Opt. Soc. Am. A 19 1524–1531 (2002). 13. J. Villa, I. De la Rosa and G. Miramontes, “Phase recovery from a single fringe pattern using an orientational vector-field-regularized estimator,” J. Opt. Soc. Am. A 22 2766–2773 (2005). 14. M. Servin and R. Rodriguez-Vera, “Two-dimensional phase locked loop demodulation of interferograms,” J. Mod. Opt. 4


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: I(x, y) = a(x, y) + b(x, y) cos[φ (x, y)], (1) 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: φ (x, y) = φ 0 (x, y) + ωx, 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 0 (η, ξ ) 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 φ (x, y), ω x and ω y at the site (x, y): where P(x, y, η, ξ ) = φ (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.
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 2 ; 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 2 , 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 φ ∈ C 2 .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.

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 2 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 φ (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 φ (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: and Then we define the discrete set Ω = {ω 1 , ω 2 } 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 φ (x) = 20.We also show the series S 1 and S 2 , which corresponds to the minimum values of the cost function (14).In graph (b) we see that exists only one frequency ω(x) ∈ Ω = {S 1 , S 2 } that minimizes Eq. ( 20).
Because we want that the estimated phase φ (x) belongs to the function space C 2 , 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 φ (•) 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 Ω = {ω 1 , ω 2 }.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 ω(x) with the estimator of Eq. ( 22) as follow: which is continuous in curvature since the frequency ω(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 ω 0 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.

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.4. For each subsequent site x + to the left and right, that isn't estimated, do the following: (a) Construct the frequency discrete set Ω = {ω 1 , ω 2 } using Eq. ( 17) and ( 18).
(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: Fig. 3.In graph 1) a signal is shown which phase is modeled as φ (x) = 1 2 ax 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.
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.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 φ with the FCT is φ ∈ C 2 .Let us demodulate the signal showed in 3(a) which its modulating phase is modeled as φ (x) = 1 2 ax 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 φ with the RPT is φ / ∈ C 2 while the FCT obtains a phase φ ∈ C 2 .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 φ ∈ C 2 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 φ ∈ C 2 .
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.

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.
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 φ (η, ξ ) = 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.Frame path 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.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.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.

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.

Discussion and Conclusions
With this paper we have demonstrated that estimating the phase into a function space C 2 , 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 2 , 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 φ as φ ∈ C 2 ).This is because phase profiles which are in C 2 , are typical phase profiles of single fringe pattern images with closed fringes.Furthermore, nowadays there is not one published work 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. 1 .
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.

Fig. 4 .
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.

Fig. 6 .
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.

Fig. 7 .
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.

Fig. 8 .
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.