4278

OPTICS LETTERS / Vol. 37, No. 20 / October 15, 2012

Fringe demodulation using the two-dimensional phase differencing operator Gannavarpu Rajshekhar and Pramod Rastogi* Applied Computing and Mechanics Laboratory, Ecole Polytechnique Fédérale de Lausanne, Lausanne 1015, Switzerland *Corresponding author: [email protected] Received June 22, 2012; revised September 4, 2012; accepted September 4, 2012; posted September 5, 2012 (Doc. ID 171145); published October 11, 2012 The Letter proposes a method for phase estimation from a fringe pattern. The proposed method relies on a parametric approach where the phase is locally approximated as a two-dimensional (2D) polynomial, with the ensuing polynomial coefficients as the respective parameters. These coefficients are then estimated using the phase differencing operator. Because of the 2D formulation, the proposed method simultaneously analyzes signal samples along the horizontal and vertical dimensions, which enables robust estimation in the presence of noise. In addition, the method directly provides the desired phase without the requirement of complex unwrapping algorithms. Simulation and experimental results are presented to validate the method’s potential. © 2012 Optical Society of America OCIS code: 120.5050, 090.2880, 120.2650.

Reliable estimation of phase from a fringe pattern is an important problem in fringe analysis [1]. The phaseshifting approach [2] is a popular technique for phase estimation; however, its application is limited by the requirement of multiple frames. For phase estimation using a single frame, several methods have been proposed, including Fourier transform [3], regularized phase tracking [4], wavelet transform [5], windowed Fourier transform (WFT) [6], one-dimensional (1D) or univariate polynomial phase approximation method [7], etc. In the Letter, we propose a fringe analysis method based on two-dimensional (2D) phase differencing operator for phase estimation. The 2D or bivariate polynomial modeling provides the generalization of the corresponding univariate approach with row(column)-by-row (column) operations on a fringe pattern, and offers better robustness against noise. The method is presented in the context of digital holographic interferometry (DHI) [8], though its applicability to other optical techniques is also discussed. The reconstructed interference field in DHI can be expressed as Γ
x; y  a
x; y expjϕ
x; y  η
x; y;

(1)

where a
x; y is the amplitude, assumed to be slowly varying, ϕ
x; y is the phase, and η
x; y represents noise. The real part of Γ
x; y constitutes the fringe pattern. Here, x and y represent pixels along the horizontal(column) and vertical(row) dimensions. For an N x × N y fringe pattern, we have x ∈ 0; N x − 1 and y ∈ 0; N y − 1. In the proposed method, the phase is modeled as a local 2D polynomial. Accordingly, Γ
x; y is divided into multiple rectangular segments or blocks bmn , each of size Lx × Ly , and with m ∈ 1; M along x and n ∈ 1; N along y, as shown in Fig. 1. For a given block bmn , we have

where I  f0 ≤ k; l and 0 ≤ k  l ≤ dg, with d being the degree of the polynomial. Equivalently, the phase is composed of d  1 layers, i.e., ϕmn
x;y  fc
0;0g fc
0;1yc
1;0xg  |


{z


} |













{z













} layer 1

layer 2

fc
0;dyd c
1;d−1xyd−1 c
d−1;1xd−1 yc
d;0xd g; |




















































{z




















































} layer d1

(4) th

where the i layer corresponds to a bivariate polynomial of degree i − 1, and has i coefficients. Consequently, phase retrieval is reduced to a parameter estimation problem, and the required parameters are the coefficients for each layer. To estimate the bivariate coefficients, the phase differencing operator [9,10] is used in the analysis. For the block bmn , the 2D phase differencing operator is given as PDx
P ;y
d−P−1 Γmn
x; y ( ) d−P−1
q  d−P−1 P Y Y P
†
pq 
; fΓmn
x  pτx ; y  qτy g p  q0

p0

Γmn
x; y  amn
x; y expjϕmn
x; y  ηmn
x; y; (2) and ϕmn
x; y 

X

c
k; lxk yl ;

(3)


k;l∈I

0146-9592/12/204278-03$15.00/0

Fig. 1. 2D block representation. © 2012 Optical Society of America

(5)

October 15, 2012 / Vol. 37, No. 20 / OPTICS LETTERS

with
†
pq Γmn

 

Γmn ; p  q even ; Γmn ; p  q odd

(6)

where “” denotes the complex conjugate, P is a parameter that varies in 0; d − 1, τx  Lx ∕
P  1, τy  Ly ∕
d − P; and
nr  n! ∕
r!
n − r!. An important property of the phase differencing operator is that it transforms a polynomial phase signal of degree d into a 2D complex exponential signal. In other words, ignoring the scaling term for the simplicity of analysis, we have

4279

Subsequently, the phase ϕmn is constructed from these coefficients using Eq. (4). Following the above procedure, the phase is estimated for all blocks bmn shown in Fig. 1. The next step is to align the estimated phases between adjacent blocks to obtain a smooth 2D phase distribution. This is achieved by a phase offsetting operation. Consider the phase ϕmn
x; y within an arbitrary block bmn of size Lx × Ly , so that x ∈ 0; Lx − 1 and y ∈ 0; Ly − 1. To ensure phase continuity along x near the end of bmn and the beginning of bm1n , the difference between the phase values of ϕmn at x  Lx and ϕm1n at x  0 should be minimal. Hence, we have L −1

PDx
P ;y
d−P−1 Γmn
x; y ≈ exp j
ωx x  ωy y;

(7)

ϕoffset

y 1 X 
ϕ
L ; y − ϕm1n
0; y; Ly y0 mn x

(13)

where the frequencies are related to the coefficients as ωy ; c
P; d − P 
−1d−1 P!
d − P!τxP τyd−P−1

ϕˆ m1n
x; y  ϕm1n
x; y  ϕoffset : (8)

c
P  1; d − P − 1 

ωx :
−1d−1
P  1!
d − P − 1!τxP τyd−P−1

(9)

The above equations signify that the respective coefficients can be obtained by estimating the spatial frequencies corresponding to the phase differencing operator, i.e., PDx
P ;y
d−P−1 Γmn
x; y. For frequency estimation, the 2D discrete Fourier transform is used:  X X ˆ x ; ωˆ y  arg max PDx
P ;y
d−P−1 Γmn
x; y ω ω1 ;ω2

y

(14)

The above offsetting operation ensures that the phase varies smoothly along x near the edges of the blocks. Subsequently, in a similar manner, the offsetting operation is performed along y to obtain the overall continuous phase distribution. It needs to be noted that the unwrapped phase estimate is directly obtained. To test the performance of the proposed method, complex reconstructed interference field Γ
x; y (512× 512 pixels) was simulated in MATLAB. Additive white

x

2  × exp−j
ω1 x  ω2 y ;

(10)

where ω1 and ω2 denote the frequencies at which the Fourier spectrum is computed. By varying P ∈ 0; d − 1 in Eqs. (8) and (9), the coefficients c
0; d; c
1; d − 1;    ; c
d − 1; 1; c
d; 0 for the layer d  1 in Eq. (4) are computed. Subsequently, the contribution of these coefficients is removed from the phase ϕmn
x; y using a peeling operation as d−1 Γmn
x; y

  X d P d−P :  Γmn
x; y exp −j c
P; d − Px y P0

(11) d−1
x; y to a 2D polynomial The above operation reduces Γmn phase signal of degree d − 1, i.e., with d layers of coefficients. Hence, by successive application of the phase differencing operator followed by a peeling operation, the coefficients corresponding to each layer are computed. Finally, the zero-order coefficient for the block bmn is computed as

 Ly −1 L x −1 X 1 X Γ0mn
x; y : c
0; 0  angle Lx Ly y0 x0 

(12)

Fig. 2. (Color online) (a) Simulated fringe pattern. (b) Estimated phase in radians. Estimation error in radians for (c) proposed method, (d) 1D polynomial phase approximation method, and (e) WFT method.

4280

OPTICS LETTERS / Vol. 37, No. 20 / October 15, 2012 Table 1. Performance at Various Noise Levels

Noise σ η 0.1 0.2 0.3 0.4

RMSE for 2D Method

RMSE for 1D Method

0.0310 0.0328 0.0537 0.0776

0.0701 0.1171 0.3220 0.4663

Gaussian noise was added with zero mean and standard deviation σ η  0.25. The corresponding fringe pattern is shown in Fig. 2(a). To apply the proposed method, the signal was divided into four blocks along each dimension, i.e., M  4 and N  4. Within each block, the phase was approximated as a quartic-order 2D polynomial, i.e., d  4. The estimated phase and the corresponding estimation error in radians are shown in Figs. 2(b) and 2(c). The root mean square error (RMSE) for the proposed method was 0.0477 rad. The required computational time was about 56 s. For comparison, phase was also estimated using the 1D polynomial phase approximation method [7] and WFT method [6]. For the univariate approach, each column was divided into four segments and quartic phase approximation was used. The estimation error in radians is shown in Fig. 2(d), with the corresponding RMSE being 0.1804 rad, and the required computational time being about 6 s. The RMSEs in radians, obtained by the proposed 2D method and the 1D method at varying levels of noise, are shown in Table 1. These results clearly indicate that the proposed 2D approach offers better estimation accuracy in the presence of noise. For the WFT method, the implementation parameters were selected as wxh  wyh  1, wxl  wyl  −1, wxi  wyi  0.05, and σ x  σ y  10. The corresponding estimation error in radians is shown in Fig. 2(e). The RMSE for the WFT method was 0.2801 rad, and the required computational time was about 210 s. These results demonstrate the superior performance of the proposed method. Further, the practical applicability of the proposed method was validated through a DHI experiment. The experimental fringe pattern is shown in Fig. 3(a). The estimated phase in radians using the proposed method is shown in Fig. 3(b). The proposed method can be extended to other optical interferometric techniques too. Applying a normalization operation [11] and real to

Fig. 3. (Color online) (a) Experimental fringe pattern. (b) Estimated phase in radians.

analytic signal conversion [12] in presence of a carrier, we obtain a complex signal similar to the one in Eq. (1), and the proposed method can be subsequently applied. To conclude, an elegant fringe analysis method based on 2D phase differencing operator was proposed for phase estimation. The method offers robust estimation without the need of unwrapping operations. The presented results illustrate the method’s potential. This work was funded by Swiss National Science Foundation under Grant 200020-131900. References 1. G. Rajshekhar and P. Rastogi, Opt. Laser Eng. 50, iii (2012). 2. K. Creath, in Progress in Optics, E. Wolf, ed. (NorthHolland, 1988), p. 349. 3. M. Takeda, H. Ina, and S. Kobayashi, J. Opt. Soc. Am. 72, 156 (1982). 4. M. Servin, J. L. Marroquin, and F. J. Cuevas, Appl. Opt. 36, 4540 (1997). 5. L. R. Watkins, S. M. Tan, and T. H. Barnes, Opt. Lett. 24, 905 (1999). 6. Q. Kemao, Opt. Laser Eng. 45, 304 (2007). 7. S. S. Gorthi and P. Rastogi, Rev. Sci. Instrum. 80, 073109 (2009). 8. U. Schnars and W. P. O. Juptner, Meas. Sci. Tech. 13, R85 (2002). 9. B. Friedlander and J. M. Francos, IEEE Trans. Image Process. 5, 1084 (1996). 10. J. M. Francos and B. Friedlander, Multidimens. Syst. Signal Process. 9, 173 (1998). 11. J. A. Quiroga, J. Antonio Gomez-Pedrero, and A. Garcia-Botella, Opt. Commun. 197, 43 (2001). 12. S. Lawrence Marple, Jr., IEEE Trans. Signal Process. 47, 2600 (1999).