Open Access
Add to BibSonomyAbstract
It has been demonstrated that the vectorial fringe-direction field is very important to demodulate fringe patterns without a dominant (or carrier) frequency. Unfortunately, the computation of this direction-filed is by far the most difficult task in the full interferogram phase-demodulation process. In this paper we present an algorithm to estimate this fringe-direction vector-field of a single n-dimensional fringe pattern. Despite that our theoretical results are valid at any dimension in the Euclidean space, we present some computer-simulated results in three dimensions because it is the most useful case in practical applications. As herein demonstrated, our method is based on linear matrix and vector analysis, this translates into a low computational cost.
©2010 Optical Society of America
1. Introduction
It is widely known in signal processing that the determination of the signal in quadrature is of main importance to extract the phase of a signal. For example, if we are dealing with a single n-dimensional fringe pattern which can be represented by
where x = (x 1, x 2, …,xn) is the coordinate vector in the region of valid data L, ϕ(x) the phase of interest, a(x) the background illumination, and b(x) the contrast. The last two terms are usually, for convenience, suppressed by means of a normalization procedure such that from now in the text we may consider f(x) ≈ cosϕ(x). The signal in quadrature fq(x) = sinϕ (x) is useful because the phase can be determined by means of the inverse tangent function. Processing a single fringe pattern without a dominant frequency the vortex operator [1] can be a solution to recover fq(x) in the two-dimensional case, however, for more dimensions this operator is not obviously established. Fortunately, for the general case of n-dimensions, the signal in quadrature can be determined by means of the general n-dimensional quadrature transform [2] which is is defined as
where
This vector contains the direction cosines that point out to the fringe direction, where e k are the standard vectors in ℝn. The key point using this transformation is the determination of n ϕ, however, the direct access to this vector field is not available. The first approximation to it can be by means of the fringe pattern gradient ∇f, defining the following vector field:
This relation indicates that the unit vector field n f is parallel to n ϕ but it changes the sign at every fringe contour. This simple difference between these vector fields implies a very difficult problem to solve and has been a widely studied topic in two-dimensions [3, 4, 5]. The relation between n f and n ϕ can also be established in the two-dimensional case defining the angles θ and α which represent the orientation and direction angles respectively, where
and W (2θ) = W (2α). The symbol W represents the wrapping operator. This relation indicates that α can be computed from θ by means of an unwrapping process [4]. Now, the vector field n ϕ is then defined as
For more than two dimensions, however, the relation between the angles of n ϕ and n f is not so obvious. An alternative solution is by determining the fringe pattern quadrature sign QS{f} = −sgn(sinϕ) which can be obtained with [6]
where
represents orientation angle subtended by the fringe pattern gradient projection on the (k,k+1) plane with the kth coordinate axis, and φk the direction angle which can obtained as in reference [4]. There are two main inconveniences using the methodology of reference [6] to determine the fringe direction: first, several unwrapping processes must be performed, and second, the unwrapping method itself is slow and complicated due to the algorithm to solve a nonlinear system in order to minimize a regularized cost-function. An improvement of the algorithm reported in [6] was proposed in reference [7], the cost function were simplified making some approximations and the time processing reduced, however, the optimization of the proposed cost function still requires to solve a nonlinear system.
2. The -dimensional regularized fringe direction-estimator
As mentioned before, n ϕ and n f are parallel but n f changes the sign with respect to n ϕ at every fringe contour. The key idea of the method presented here is based on our previously reported work [5]. For ℝn consider we compute from n f a set of n − 1 unit vectors d k, where k = 1,2…,n−1, and
such that n f and d k form an orthonormal basis for ℝn. The set of vectors d k can be obtained from n f in the following way: When calculating the null space of n f by means of its QR decomposition, Q will be formed by a set of orthonormal column vectors, that is Q=(a 1 a 2…a n) where a 1 and n f are parallel, so that set d k can be selected as d 1 = a 2, d 2 = a 3, …,d (n−1) = a n. By observing Figure (1), which is the case for ℝ3, we note that n ϕ ⊥ d k. Once d k is computed for all x ∊ L the idea is to compute a smooth vector field p(x) = (p 1, p 2, …, pn)T that points out to the same direction of n ϕ(x). The first restriction to construct our estimator is that d k(x) ⊥ p(x), or which is the same
On the other hand, to avoid abrupt sign changes we most restrict p(x) to be smooth. Taking into a count these restrictions, the strategy of our algorithm requires consider a subset Γ ∊ L, which contains already estimated sites around a given site x to be estimated. The vector field p(x) can be locally estimated by means of a scanning strategy and the following cost function
Fig. 1. Relation between vectors n ϕ and n f in a 3D fringe pattern. (A) They point out to the same direction, (B) or they have opposite sign, but they are parallel at every site in the fringe pattern. Note that d k and n f form an orthonormal set.
In this cost function x̃ represents the coordinates in the region Γ rounding x, p(x̃) the already estimated vectors in Γ, s(x̃) a Boolean function used to indicate if the site x̃ ∊ Γ has already been estimated, and µ a regularization parameter that controls the smoothness of the estimated vector field. To compute p in a given site x we set ∇U x(p) = 0, which represents a simple linear system of n equations that can be written in matrix form as
where
To estimate the full vector field p(x) in L we start by setting the field s(x) = 0 (x∀L), then the linear system (12) solved for every site in L. By observing our algorithm we note that in the first site to be estimated, Equation 12 represents an homogeneous system, so it is necessary to estimate p otherwise. In practice we only select p = n f. Once a site x is estimated the corresponding indicator s(x) is set equal to 1. As the estimation of p(x) in a given site requires already estimated values of neighbors, it is necessary a proper scanning strategy. One way to realize it is by following fringe contours, for this reason we use a quality map based scanning as the reported in [8] for phase unwrapping. For our purposes we use as the quality map the magnitude of fringe pattern gradient. Unlike previously reported methods for direction estimation [6, 7], from the computational point of view our method has the following advantages: once the region Γ has been defined, the processing time is fixed for every site in the fringe image and it works efficiently because a simple linear system Ap = b have to be solved by means of any direct method, being A of size n×n. This is not the case for methods in references [6, 7] that require to solve a non-linear system by means of iterative methods.
The following algorithm describes our method for fringe direction estimation.
- Compute n f and d k, and set s = 0 for every site in the fringe image. Define the size of Γ.
- Choose a site in the field n f for the first estimation, use p = n f for such a site and set s = 1.
- Compute p from Equation 12 in an adjacent neighbor of a previously computed site (for example, following the scanning strategy reported in [5]), and set s = 1.
- Repeat step (3) until all sites are computed.
3. Numerical experiments
In this section we present the results obtained applying our method for estimating direction-fields of three-dimensional fringe patterns. In the two following experiments we used a size of 7 × 7 × 3 for Γ (in the x,y and z directions respectively), and µ = 1. The first experiment was a simple 100 × 100 noisy simulated phase-shifted fringe pattern with N = 50 equally spaced phase-steeps ranging from 0 to 2π. For this experiment we used uniformly-distributed additive noise ranging from 0 to 1. The modulating phase was a semi-sphere which generates circular fringes. The sequence were generated according to I(x, y, z) = cos[ϕ(x,y)+κ(z)]. The function κ(z) defines the phase shift such that , where z = 0,1, …, N − 1. The phase ϕ(x,y) was calculated with
Figure 2 (A) shows some samples of the sequence where the z-axis indicates the phase-shift, while Figure 2 (B) shows the corresponding gray-level-codified direction-angles computed with the proposed method. Black represents −π rad and white π rad. The processing time in this experiment was about 88 seconds (using a 2.4 GHz Pentium D based computer), and the direction angles were computed using
In this experiment we carried out a quantitative evaluation of our fringe direction-estimator computing the normalized mean-square error (NMSE), which is defined as
where n ϕ is the theoretical fringe-direction vector-field. In this case the error were ε = 0.0055. It is important to remark that, as mentioned by Larkin [1], the interferogram demodulation does not require an accurate estimate of the fringe direction-field. The second was a simulated load-stepping photoelastic experiment using the model of a circular disc under diametral compression to evaluate the relative retardation [9, 10]. Figure 3 (A) shows some samples of a sequence increasing the load compression. In this case it was a 180 × 400 image size with 20 load-steeps. Figure 3 (B) shows the corresponding three-dimensional phase-map using our n-dimensional fringe direction-estimator and the quadrature transform [2]. In this experiment the computation of the fringe-direction required bout 233 seconds.
Fig. 2. (A) Sequence of phase-shifted interferograms where z-axis indicates the phase-shift. (B) Gray-level-codified direction-maps (Black represents −π rad and white π rad).
Fig. 3. (A) Simulated photoelastic fringe patterns by load-steeping. (B) Three-dimensional phase-map obtained using the n-dimensional fringe direction-estimator and the quadrature transform [2].
4. Conclusions
We have proposed a method to determine the fringe-direction vector-field of a single n-dimensional fringe pattern. Our proposed theoretical approach was validated presenting some simulated experiments of three-dimensional fringe patterns. As mentioned, the fringe direction estimation of a n-dimensional fringe pattern is by far the most difficult task in the process of phase demodulation, even more, for more than two-dimensions it can be a strong computational effort. Unlike already reported techniques, our proposal can be easily described and implemented by means of linear vector and matrix analysis, and can be understood naturally regardless of the problem’s dimension which allows the possibility of being applied in problems of future research. An additional attractive feature of our method is that in the demodulation of interferograms does not require a precise estimate of the fringe-direction vector-field, so it can be used in many real applications.
Acknowledgements
We acknowledge the support for the realization of this work to the Consejo Nacional de Ciencia y Tecnología (CONACYT), México, through the project CB-2008-01/102041, and the Ministerio de Ciencia e Innovación of Spain trough the project DPI2009-09023.
References and links
1. K. G. Larkin, D. J. Bone, and M. 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). [CrossRef]
2. M. Servín, J. A. Quiroga, and J. L. Marroquín, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20, 925–934 (2003). [CrossRef]
3. X. Zhou, J. P. Baird, and J. F. Arnold, “Fringe orientation estimation by use of a gaussian gradient-filter and neighboring-direction averaging,” Appl. Opt. 38,, 795–804 (1999). [CrossRef]
4. J. A. Quiroga, M. Servín, and F. 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). [CrossRef]
5. Jesús Villa, Ismael De la Rosa, Gerardo Miramontes, and Juan Antonio Quiroga, “Phase recovery from a single fringe pattern using an orientational vector field regularized estimator,” J. Opt. Soc. Am. A 22, 2766–2773 (2005). [CrossRef]
6. J. A. Quiroga, Manuel Servín, J. Luis Marroquín, and Daniel Crespo “Estimation of the orientation term of the general quadrature transform from a single n-dimensional fringe pattern,” J. Opt. Soc. Am. A 22, 439–444 (2005). [CrossRef]
7. D. Crespo, J. A. Quiroga, and J. A. Gomez-Pedrero, “Fast algorithm for estimation of the orientation term of a general quadrature transform with application to demodulation of an n-dimensional fringe pattern,” Appl. Opt. 43,, 6139–6146 (2004). [CrossRef]
8. B. Ströbel, “Processing of interferometric phase maps as complex-valued phasor images,” Appl. Opt. 35, 2192–2198 (1996). [CrossRef] [PubMed]
9. Ng TW., “Photoelastic stress analysis using an object steep-loading method,” Exp. Mech. 37, 137–141 (1997). [CrossRef]
10. J. Villa, J. A. Quiroga, and J. A. Gómez-Pedrero, “Measurement of retardation in digital photoelasticity by load stepping using a sinusoidal least-squares fitting,” Opt. Las. Eng. 41, 127–137 (2004). [CrossRef]
