N-dimensional regularized fringe
dir ection-estimator

Jeśus Villa,1,∗ Juan Antonio Quiroga,2 Manuel Serv́ın,3 Julio Cesar
Estrada,3 and Ismael de la Rosa1

1Laboratorio de Procesamiento Digital de Señales, Facultad de Ingenierı́a Eléctrica,
Universidad Aut́onoma de Zacatecas, Avenida Ramón López Velarde 801, Zacatecas 98000,

México.

2Departamento déOptica, Facultad de Ciencias Fı́sicas, Universidad Complutense de
Madrid, Ciudad Universitaria S/N, Madrid 28040, España.

3Centro de Investigaciones eńOptica A.C., Loma del Bosque 115, Col. Lomas del Campestre,
León, Guanajuato 37150, Ḿexico.

∗jvillah@uaz.edu.mx

Abstract: 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 singlen-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

OCIS codes: (100.2650) Fringe analysis; (100.5070) Phase retrieval; (120.5050) Phase
Measurement.

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. A18,1862–1870 (2001).
2. M. Serv́ın, J. A. Quiroga, and J. L. Marroquı́n, “Generaln-dimensional quadrature transform and its application

to interferogram demodulation,” J. Opt. Soc. Am. A20,925–934 (2003).
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).
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. A19,1524–1531 (2002).
5. Jeśus 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. A22, 2766–2773
(2005).

6. J. A. Quiroga, Manuel Servı́n, J. Luis Marroqúın, and Daniel Crespo “Estimation of the orientation term of
the general quadrature transform from a singlen-dimensional fringe pattern,” J. Opt. Soc. Am. A22, 439-444
(2005).

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 ann-dimensional fringe pattern,” Appl. Opt.
43,, 6139–6146 (2004).

8. B. Str̈obel, “Processing of interferometric phase maps as complex-valued phasor images,” Appl. Opt.35,2192–
2198 (1996).

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16567



9. Ng TW., “Photoelastic stress analysis using an object steep-loading method,” Exp. Mech.37,137–141(1997).
10. J. Villa, J. A. Quiroga and J. A. Ǵomez-Pedrero, “Measurement of retardation in digital photoelasticity by load

stepping using a sinusoidal least-squares fitting,” Opt. Las. Eng.41,127–137 (2004).

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

f (x) = a(x)+b(x)cosφ(x), x = (x1,x2, ...,xn) ∈ L, (1)

wherex = (x1,x2, ...,xn) is the coordinate vector in the region of valid dataL, φ(x) the phase
of interest,a(x) the background illumination, andb(x) the contrast. The last two terms are usu-
ally, for convenience, suppressed by means of a normalization procedure such that from now in
the text we may considerf (x) ≈ cosφ(x). The signal in quadraturefq(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 ofn-dimensions, the signal in quadra-
ture can be determined by means of the generaln-dimensional quadrature transform [2] which
is is defined as

Qn{ f (x)} = nφ ·
∇ f

‖∇φ‖
, (2)

where

nφ =
∇φ
‖∇φ‖

=

n
∑

k=1
(∂φ/∂xk)ek

‖∇φ‖
=

n

∑
k=1

cos(αk)ek. (3)

This vector contains the direction cosines that point out to the fringe direction, whereek are
the standard vectors inRn. The key point using this transformation is the determination ofnφ ,
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:

n f =
∇ f

‖∇ f‖
=

−sinφ
n
∑

k=1
(∂φ/∂xk)ek

|sinφ |‖∇φ‖
= −sgn(sinφ)nφ . (4)

This relation indicates that the unit vector fieldn f is parallel tonφ 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
betweenn f andnφ can also be established in the two-dimensional case defining the anglesθ
andα which represent the orientation and direction angles respectively, where

θ = tan−1
(

∂ f/∂x2

∂ f/∂x1

)

, θ ∈
[

−
π
2

,
π
2

)

, (5)

andW(2θ) = W(2α). The symbolW 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

nφ = cos(α)e1 +sin(α)e2. (6)

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16568



For more than two dimensions, however, the relation between the angles ofnφ andn f is not so
obvious. An alternative solution is by determining the fringe pattern quadrature signQS{ f} =
−sgn(sinφ) which can be obtained with [6]

QS{ f} = sgn[cos(βk)]sgn[cos(ϕk)], (7)

where

βk = tan−1
(

∂ f/∂xk+1

∂ f/∂xk

)

(8)

represents orientation angle subtended by the fringe pattern gradient projection on the(k,k+1)
plane with thekth 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. Then-dimensional regularized fringe direction-estimator

As mentioned before,nφ and n f are parallel butn f changes the sign with respect tonφ at
every fringe contour. The key idea of the method presented here is based on our previously
reported work [5]. ForRn consider we compute fromn f a set ofn−1 unit vectorsdk, where
k = 1,2. . . ,n−1, and

d1 =(d11,d12, . . . ,d1n)
T ,

d2 =(d21,d22, . . . ,d2n)
T ,

...

dn−1 =(d(n−1)1,d(n−1)2, . . . ,d(n−1)n)
T ,

(9)

such thatn f anddk form an orthonormal basis forRn. The set of vectorsdk can be obtained
from n f in the following way: When calculating the null space ofn f by means of its QR de-
composition, Q will be formed by a set of orthonormal column vectors, that isQ= (a1 a2 . . .an)
wherea1 andn f are parallel, so that setdk can be selected asd1 = a2, d2 = a3, . . . ,d(n−1) = an.
By observing Figure (1), which is the case forR

3, we note thatnφ ⊥ dk. Oncedk is com-
puted for allx ∈ L the idea is to compute a smooth vector fieldp(x) = (p1, p2, . . . , pn)

T that
points out to the same direction ofnφ (x). The first restriction to construct our estimator is that
dk(x) ⊥ p(x), or which is the same

dk(x) ·p(x) = 0, x∀L. (10)

On the other hand, to avoid abrupt sign changes we most restrictp(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 sitex to be estimated. The vector field
p(x) can be locally estimated by means of a scanning strategy and the following cost function

Ux(p) = ∑̃
x∈Γ

{
n−1

∑
k=1

[p(x) ·dk(x̃)]2 + µ‖p(x)−p(x̃)‖2s(x̃)}. (11)

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16569



nf

fn = -nf

d1

d2

fn = nf

d1

d2

nf d1

d2

fn d1

d2

nf d1

nf d2

nf

fn

(A) (B)

Fig. 1. Relation between vectorsnφ andn 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 thatdk andn f form an orthonormal set.

In this cost functioñx represents the coordinates in the regionΓ roundingx, p(x̃) the already
estimated vectors inΓ, s(x̃) a Boolean function used to indicate if the sitex̃ ∈ Γ has already
been estimated, andµ a regularization parameter that controls the smoothness of the estimated
vector field. To computep in a given sitex we set∇Ux(p) = 0, which represents a simple linear
system ofn equations that can be written in matrix form as

Ap = b, (12)

where

A =















∑
x̃∈Γ

{∑n−1
k=1 dk1(x̃)2+µs(x̃)} ∑

x̃∈Γ
{∑n−1

k=1 dk2(x̃)dk1(x̃)} ... ∑
x̃∈Γ

{∑n−1
k=1 dkn(x̃)dk1(x̃)}

∑
x̃∈Γ

{∑n−1
k=1 dk1(x̃)dk2(x̃)} ∑

x̃∈Γ
{∑n−1

k=1 dk2(x̃)2+µs(x̃)} ... ∑
x̃∈Γ

{∑n−1
k=1 dkn(x̃)dk2(x̃)}

...
...

...
...

∑
x̃∈Γ

{∑n−1
k=1 dk1(x̃)dkn(x̃)} ∑

x̃∈Γ
{∑n−1

k=1 dk2(x̃)dkn(x̃)} ... ∑
x̃∈Γ

{∑n−1
k=1 dkn(x̃)2+µs(x̃)}















, (13)

and

b =















µ ∑
x̃∈Γ

{p1(x̃)s(x̃)}

µ ∑
x̃∈Γ

{p2(x̃)s(x̃)}

...
µ ∑

x̃∈Γ
{pn(x̃)s(x̃)}















. (14)

To estimate the full vector fieldp(x) in L we start by setting the fields(x) = 0 (x∀L),
then the linear system (12) solved for every site inL. 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 estimatep otherwise. In practice we only selectp = n f . Once a sitex is estimated
the corresponding indicators(x) is set equal to 1. As the estimation ofp(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 systemAp = b have to be
solved by means of any direct method, beingA of sizen×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.

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16570



The following algorithm describes our method for fringe direction estimation.

(1) Computen f anddk, and sets= 0 for every site in the fringe image. Define the size ofΓ.

(2) Choose a site in the fieldn f for the first estimation, usep = n f for such a site and set
s= 1.

(3) Computep from Equation 12 in an adjacent neighbor of a previously computed site (for
example, following the scanning strategy reported in [5]), and sets= 1.

(4) 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 thex,y andz directions respectively), andµ = 1. The first experiment was
a simple 100×100 noisy simulated phase-shifted fringe pattern withN = 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 toI(x,y,z) =cos[φ(x,y)+κ(z)]. The function
κ(z) defines the phase shift such thatκ(z) = 2π

N z, wherez= 0,1, . . . ,N−1. The phaseφ(x,y)
was calculated with

φ(x,y) =
√

802− (x−50)2− (y−50)2. (15)

Figure 2 (A) shows some samples of the sequence where thez-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

θ = tan−1
(

p2

p1

)

. (16)

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

ε =
∑‖nφ −p‖2

∑‖nφ‖2 , (17)

wherenφ 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 simu-
lated 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.

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16571



Fig. 2. (A) Sequence of phase-shifted interferograms wherez-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-mapobtained using then-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 singlen-
dimensional fringe pattern. Our proposed theoretical approach was validated presenting some
simulated experiments of three-dimensional fringe patterns. As mentioned, the fringe direc-
tion estimation of an-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 compu-
tational 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 Tecnoloǵı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.

#129221 - $15.00 USD Received 28 May 2010; revised 7 Jul 2010; accepted 7 Jul 2010; published 22 Jul 2010
(C) 2010 OSA 2 August 2010 / Vol. 18,  No. 16 / OPTICS EXPRESS  16572