OSA's Digital Library

## Optics Express

• Editor: C. Martijn de Sterke
• Vol. 18, Iss. 16 — Aug. 2, 2010
• pp: 16567–16572

## N-dimensional regularized fringe direction-estimator

Jesús Villa, Juan Antonio Quiroga, Manuel Servin, Julio Cesar Estrada, and Ismael de la Rosa  »View Author Affiliations

Optics Express, Vol. 18, Issue 16, pp. 16567-16572 (2010)
http://dx.doi.org/10.1364/OE.18.016567

View Full Text Article

Acrobat PDF (723 KB)

 Vol Issue Page

Browse by Journal and Year

Lookup Conference Papers

## Article Tools

Share
Citations
 Select an action... ------------------------ Export Citation in:   ► BibTeX   ► EndNote (RIS)   ► HTML (.html)   ► Plain Text ------------------------ Save to:   ► My Article Collections

### 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 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

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

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

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]

] 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

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]

] which is is defined as

$Qn {f (x)} = nϕ· ∇f ∥ ∇ϕ∥,$
(2)

where

$nϕ= ∇ϕ ∥ ∇ϕ∥= Σ k=1n ( ∂ϕ ∂ xk) ek ∥ ∇ϕ∥= Σ k=1ncos ( αk) ek.$
(3)

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:

$nf= ∇f ∥∇f∥= −sinϕ Σ k=1n ( ∂ϕ ∂xk) ek ∣sinϕ∣∥∇ϕ∥= −sgn (sinϕ) nϕ.$
(4)

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

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]

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

$θ= tan −1 ( ∂f ∂x2 ∂f ∂x1), θ∈[− π2, π2),$
(5)

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

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]

]. Now, the vector field n ϕ is then defined as

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

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

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]

]

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

where

$βk= tan −1 ( ∂f ∂x k+1 ∂f ∂xk)$
(8)

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

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]

]. There are two main inconveniences using the methodology of reference [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]

] 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

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]

] was proposed in reference [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]

], 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

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]

]. For ℝ n consider we compute from n f a set of n − 1 unit vectors d k , where k = 1,2…,n−1, and

$d1= ( d11, d12, ..., d 1n)T,$

$d2= ( d21, d22 ,..., d 2n)T,$

$⋮$

$d n−1= ( d ( n−1)1, d ( n−1)2 ,..., d ( n−1)n)T,$
(9)

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 2a 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 xL 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

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

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

$Ux (p)= Σ x˜∈Γ { Σ k=1 n−1 [p (x)· dk ( x˜)]2+μ ∥p (x) −p ( x˜) ∥2s ( x˜)}.$
(11)
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 represents the coordinates in the region Γ rounding x, p() the already estimated vectors in Γ, s() a Boolean function used to indicate if the site ∊ Γ 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

$Ap=b,$
(12)

where

$A= ( Σ x˜∈Γ { Σ k=1 n−1 d k1 ( x˜)2+μs ( x˜)} Σ x˜∈Γ { Σ k=1 n−1 d k2 ( x˜) d k1( x˜)} ⋯ Σ x˜∈Γ { Σ k=1 n−1 d kn ( x˜) d k1( x˜)} Σ x˜∈Γ { Σ k=1 n−1 d k1 ( x˜) d k2( x˜)} Σ x˜∈Γ { Σ k=1 n−1 d k2 ( x˜)2+μs ( x˜)} ⋯ Σ x˜∈Γ { Σ k=1 n−1 d kn ( x˜) d k2( x˜)} ⋮ ⋮ ⋱ ⋮ Σ x˜∈Γ { Σ k=1 n−1 d k1 ( x˜) d kn( x˜)} Σ x˜∈Γ { Σ k=1 n−1 d k2 ( x˜) d kn( x˜)} ⋯ Σ x˜∈Γ { Σ k=1 n−1 d kn ( x˜)2+μs ( x˜)}),$
(13)
$b= ( μ Σ x˜∈Γ { p1 ( x˜)s ( x˜)} μ Σ x˜∈Γ { p2 ( x˜)s ( x˜)} ⋮ μ Σ x˜∈Γ { pn ( x˜)s ( x˜)}).$
(14)

To estimate the full vector field p(x) in L we start by setting the field s(x) = 0 (xL), 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

B. Ströbel, “Processing of interferometric phase maps as complex-valued phasor images,” Appl. Opt. 35, 2192–2198 (1996). [CrossRef] [PubMed]

] 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

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]

], 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

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]

] that require to solve a non-linear system by means of iterative methods.

The following algorithm describes our method for fringe direction estimation.

1. Compute n f and d k , and set s = 0 for every site in the fringe image. Define the size of Γ.
2. Choose a site in the field n f for the first estimation, use p = n f for such a site and set s = 1.
3. Compute p from Equation 12 in an adjacent neighbor of a previously computed site (for example, following the scanning strategy reported in [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]

]), and set s = 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 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 $κ (z)= 2πNz$ , where z = 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 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

$θ= 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)

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

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]

], 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

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]

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

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]

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

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]

].

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

 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]

OCIS Codes
(100.2650) Image processing : Fringe analysis
(100.5070) Image processing : Phase retrieval
(120.5050) Instrumentation, measurement, and metrology : Phase measurement

ToC Category:
Instrumentation, Measurement, and Metrology

History
Original Manuscript: May 28, 2010
Revised Manuscript: July 7, 2010
Manuscript Accepted: July 7, 2010
Published: July 22, 2010

Citation
Jesús Villa, Juan Antonio Quiroga, Manuel Servin, Julio Cesar Estrada, and Ismael de la Rosa, "N-dimensional regularized fringe direction-estimator," Opt. Express 18, 16567-16572 (2010)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-16-16567

Sort:  Year  |  Journal  |  Reset

### References

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. J . Villa, I . De la Rosa, G . Miramontes, and J. 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, M. Servín, J. Luis Marroquín, and D . 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. T. W. Ng, "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]

### Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

### Figures

 Fig. 1. Fig. 2. Fig. 3.

OSA is a member of CrossRef.