Total variation regularization cost function for demodulating phase discontinuities Ricardo Legarda-Saenz,* Carlos Brito-Loeza, and Arturo Espinosa-Romero Facultad de Matematicas, Universidad Autonoma de Yucatan, Mexico *Corresponding author: [email protected] Received 25 November 2013; revised 5 March 2014; accepted 6 March 2014; posted 10 March 2014 (Doc. ID 201717); published 4 April 2014
We introduce a method based on the minimization of a total variation regularization cost function for computing discontinuous phase maps from fringe patterns. The performance of the method is demonstrated by numerical experiments with both synthetic and real data. © 2014 Optical Society of America OCIS codes: (120.2650) Fringe analysis; (120.3180) Interferometry; (120.3940) Metrology; (120.6650) Surface measurements, figure; (100.3190) Inverse problems. http://dx.doi.org/10.1364/AO.53.002297
1. Introduction
The main goal of fringe analysis techniques is to recover accurately the local modulated phase from one or several fringe patterns [1,2]; such phase is related to some physical quantities, such as shape, deformation, refractive index, and temperature. The basic model for a fringe pattern is given by I x ax bx cosφx ϕx ; where x x1 ; x2 , ax is the background illumination, bx is the amplitude modulation, and ϕx is the phase map to be recovered; the spatial carrier frequency of the fringe pattern is defined by the term φx . One challenging task is the recovery of a discontinuous phase map from a single fringe pattern. In recent years, many methods that successfully estimate the phase from a single pattern [3–6] have been reported. These methods, however, consider phase maps and illumination terms as continuous. Recently, a method was proposed for computing discontinuous phase maps, from a fringe pattern with carrier frequency, based on the minimization of a regularized cost function, which uses a second-order edge-preserving potential [7]. Although this method is able to detect and reconstruct phase discontinuities, its cost functional is 1559-128X/14/112297-05$15.00/0 © 2014 Optical Society of America
not convex; hence convergence to an optimal solution is conditioned to the provided initial phase usually computed by standard methods. In this work we present a method for recovering a discontinuous phase map from a single fringe pattern based on the total variational approach. First we describe the minimization problem based on the total variational approach, where we made some remarks about the convexity of the proposed method. Then we describe the numerical solution of the proposed cost function, which results in a simple algorithm. The resulting algorithm performance is tested by two numerical phase demodulations. 2. Demodulation of Fringe Pattern Using Total Variational Approach A. Variational Formulation
Following the image denoising method proposed by Rudin et al. [8] capable of preserving discontinuities, we propose to estimate the discontinuous phase map, background, and amplitude terms as the solution of the minimization problem defined by min Eϕx ; ax ; bx ≡
ϕx ;ax ;bx
Z Z λ I x − gx 2 dx j∇ϕx jdx 2 Ω Ω Z Z j∇ax jdx j∇bx jdx ; (1) Ω
Ω
10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2297
where gx is the noise input fringe pattern, Ω ⊂ R2 denotes the continuous signal domain, and λ is a positive parameter or Lagrange multiplier to control the contribution of the similarity term. Due to each total variation regularizer in Eq. (1) being convex [9], it is easy to verify that minimization with respect to ax yields the global minimum. Obtaining global minima for bx and ϕx is, on the other hand, not that simple. A numerical analysis carried out on Eq. (1) to study its convexity showed that provided the terms I x and gx are in the range 0; 2, the second-order derivatives with respect to bx and ϕx deliver positive semi-definite functionals. In Appendix A we give insight about the positive definiteness of each term in Eq. (1). In order to minimize Eq. (1), the first-order optimality condition or Euler–Lagrange equation has to be derived for each function ϕx , ax , and bx . In the formal derivation we assume that all these vector fields are smooth enough such that gradients are well defined and the variations have compact support over Ω so that we can use the divergence theorem to get rid of the boundary terms. By following standard methods from calculus of variations, the gradient descent equations can be obtained and are given by ∂ϕx ∇ϕx ∂I ∇· λI x − gx x ; ∂t j∇ϕx j ∂ϕx ∂ax ∇ax ∂I ∇· λI x − gx x ; ∂t j∇ax j ∂ax ∂bx ∇bx ∂I ∇· λI x − gx x ; ∂t j∇bx j ∂bx
(2)
with boundary conditions ∂ϕx 0; ∂ν
∂ax 0; ∂ν
and
∂bx 0; ∂ν
where ν denotes the outer normal to the boundary. B.
Numerical Solution
To solve Eq. (2) we apply a sequential method, i.e., we solve the first equation for ϕx, and the new value is immediately used in the second equation when solving for ax ; the same procedure is applied for bx, and this process is repeated iteratively until reaching convergence to the solution. We now proceed to show how to discretize the differential operator in the first equation of Eq. (2). The same procedure can be applied to the other two equations. For simplicity we drop the subscript x from the function ϕx and let ϕi;j ϕx1 ; x2 to denote the value of a grid function ϕ at point i; j x1 ; x2 defined on the cell-centered grid Ωh equal to all points x1 ; x2 ∈ Ω such that 2i − 1hx1 2j − 1hx2 x1 ; x2 ; 2 2 with 1 ≤ i ≤ m, 1 ≤ j ≤ n. This way Ωh consists of m × n cells of size hx1 × hx2 with hx1 1∕m, hx2 1∕n the grid spacing. 2298
APPLIED OPTICS / Vol. 53, No. 11 / 10 April 2014
To approximate the derivatives, we use forward and backward finite differences defined as follows: δ x1 ϕi;j
ϕi1;j − ϕi;j h x1
and δ x2 ϕi;j
ϕi;j1 − ϕi;j : hx2
The divergence term is approximated using ∇ · V i;j δ−x1 V 1i;j δ−x2 V 2i;j ; where V i;j V 1i;j ; V 2i;j
∇ϕi;j ; j∇ϕji;j
∇ϕi;j δ x1 ϕi;j ; δx2 ϕi;j ; q 2 2 j∇ϕji;j δ x1 ϕi;j δx2 ϕi;j :
The rest of the terms in the equation are implemented by straightforward evaluation at point i; j. Finally, the Neumann’s boundary condition on ∂Ω is treated as ϕi;0 ϕi;1 ;
ϕi;n1 ϕi;n ;
ϕ0;j ϕ1;j ;
ϕm1;j ϕm;j ;
ϕ0;0 ϕ1;1 ;
ϕ0;n1 ϕ1;n ;
ϕm1;0 ϕm;1 ;
ϕm1;n1 ϕm;n :
Examples of numerical implementations similar to Eq. (2) can be found in Refs. [8] and [10]. 3. Numerical Experiments
To illustrate the performance of the proposed method, we carried out two numerical experiments using an Intel Core i5 @ 2.40 GHz laptop with Debian (wheezy) 64-bit and 4 GB of memory. The first experiment was the demodulation of a synthetic fringe pattern of size 250 × 250 pixels generated in similar way to that described in Ref. [7]: 8 a 0.6; bx 0.5; if 90 ≤ x1 ; x2 ≤ 160 > > > x > < ϕ 1.2 0.005x x Ix ; > ax 1.0; bx 0.9; otherwise > > > : ϕx 0.2 0.005x φx
2π x; P
(3)
where x1 ; x2 and P 12.5 are given in pixels. Figure 1(a) shows the synthetic fringe pattern, and Fig. 1(b) shows the synthetic phase term ϕx . The initial values for the demodulation process were ax bx 1.0, ϕx 0.2. The resultant phasemap estimation using the proposed method is shown in Fig. 2(a). The mean squared error (MSE) of the
Fig. 3. (a) Experimental fringe pattern. (b) Demodulated phase term. Fig. 1. (a) Synthetic fringe pattern. (b) Synthetic phase term ϕx .
demodulation was 0.00285, and the time employed to get the solution was 144 s using a simple gradient descent scheme with λ 100 [11–13]. For purpose of comparison, we implemented the method described in Ref. [7]. The resultant phase-map estimation is shown in Fig. 2(b). The MSE of the demodulation was 0.06722, and the time employed by this method was 800 s using the same numerical scheme. As can be observed, the estimation based on Eq. (2) successfully demodulates the fringe pattern with phase discontinuities. Because of the simple definition of the term ϕx , Eq. (3), the proposed cost function will converge to an accurate solution given an arbitrary initial point. On the contrary, the time used by the method reported in Ref. [7] was larger than the time used by our method, Eq. (2), and fails to demodulate the fringe pattern using initial points far away of the optimal point. The second experiment was the processing of a fringe pattern obtained from a holographic interferometry experiment [14], which consisted of the height measurement of a microthin film. The fringe pattern obtained from this experiment, with 640 × 480 pixels, is shown in Fig. 3(a). The fringe pattern was demodulated using Fourier-based quadrature transform [15]; the resultant demodulated phase term is shown in Fig. 3(b). During the demodulation process, the background, modulation, and spatial carrier terms are estimated; then the carrier term is subtracted from the wrapped phase term and the phase term is unwrapped if necessary [16]. These estimations are used as initial values in our method.
Fig. 2. Estimated phase terms using (a) Eq. (2) and (b) method from Ref. [7].
Fig. 4. Estimated phase terms using (a) Eq. (2) and (b) Schwider– Hariharan 41 algorithm [17].
The resultant phase-map estimation is shown in Figs. 4(a) and 5. The time employed for the proposed method was 210 s using a simple gradient descent scheme with λ 100 [11–13]. For purposes of comparison, we show in Figs. 4(b) and 6 the resultant phase-map estimation using a 41 phase shifting method [17].
Fig. 5. Estimated phase terms using Eq. (2).
Fig. 6. Estimated phase terms using Schwider–Hariharan 41 algorithm [17]. 10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2299
The second-order derivatives of F with respect to ax , bx , and ϕx read as follows: d2 F 2 da2x d2 F 2 db2x d2 F 2 dϕ2x Fig. 7. Columns of the estimated phase terms using Eq. (2) and Schwider–Hariharan 4 1 algorithm [17].
As can be observed in Figs. 4 and 7, the estimated phase map obtained by the solution of Eq. (2) is very similar to the one obtained using a well-proven demodulation method. An important point to highlight from the result shown in Fig. 7 is the preservation of the dynamic range of the estimated phase map. In Ref. [18] the authors showed that the change in image intensity (here the dynamic range of the phase) due to total variation regularization is directly proportional to the regularization parameter and inversely proportional to the scale of the image feature; that means total variation regularization causes that smaller-scaled changes in phase be partially or entirely removed while larger-scaled ones be relatively unaffected. The resultant phase map of our method can be seen as the filtered phase term of the phase shifting method, as one can observe in Figs. 5 and 6. Another important aspect of this experiment is the need for using a prior estimate for demodulating real fringe patterns because the optimization process could be trapped on a local minimum due to large phase variations. 4. Conclusion
As can be seen from the resultant phase estimations in the above two experiments, the proposed cost function allows the accurate demodulation of a single fringe pattern with discontinuities. The numerical solution of Eq. (2) results in a very simple algorithm that can be implemented with several computational efficient techniques [10]. An extra advantage of the proposed method is its feasibility to be implemented on dedicated hardware to obtain real-time processing. This will be one aim of our future research. Appendix A
Let us define F as Z λ I − gx 2 dx: F 2 Ω x Z ax bx cosφx ϕx − gx 2 dx: Ω
2300
APPLIED OPTICS / Vol. 53, No. 11 / 10 April 2014
Z Ω
Z Ω
dx > 0;
(A1)
cos2 φx ϕx dx ≥ 0;
(A2)
Z Ω
gx − ax cosφx ϕx
− bx cos2φx ϕx dx.
(A3)
From Eq. (A1), it is clearR that F is convex with respect to ax , and due to Ω j∇ax jdx being convex as well, as proven in [9], minimization with respect to ax yields the global minimum. In a similar fashion, Eq. (A2) reveals that F is at most a semi-positive definite functional. Lately, from Eq. (A3), not much can be said about F with respect to ϕx . However, through a large number of numerical tests we noticed that provided I x and gx are in the range 0; 2, Eq. (A3) is not negative and therefore F is a semi-positive definite functional. Although this observation is by no means conclusive, it is a good practice to set I x and gx to this interval. References 1. D. W. Robinson and G. T. Reid, Interferogram Analysis, Digital Fringe Pattern Measurement Techniques (Taylor & Francis, 1993). 2. G. Rajshekhar and P. Rastogi, “Fringe analysis: premise and perspectives,” Opt. Lasers Eng. 50, iii–x (2012). 3. J. L. Marroquin, M. Rivera, S. Botello, R. Rodriguez-Vera, and M. Servin, “Regularization methods for processing fringepattern images,” Appl. Opt. 38, 788–794 (1999). 4. J. Villa, J. A. Quiroga, and M. Servin, “Improved regularized phase-tracking technique for the processing of squared-grating deflectograms,” Appl. Opt. 39, 502–508 (2000). 5. R. Legarda-Saenz, W. Osten, and W. P. Juptner, “Improvement of the regularized phase tracking technique for the processing of nonnormalized fringe patterns,” Appl. Opt. 41, 5519–5526 (2002). 6. M. Rivera, “Robust phase demodulation of interferograms with open or closed fringes,” J. Opt. Soc. Am. A 22, 1170– 1175 (2005). 7. C. Galvan and M. Rivera, “Second-order robust regularization cost function for detecting and reconstructing phase discontinuities,” Appl. Opt. 45, 353–359 (2006). 8. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). 9. R. Acar and C. R. Vogel, “Analysis of bounded variation penalty methods for ill-posed problems,” Inverse Probl. 10, 1217–1229 (1994). 10. P. Getreuer, “Rudin-Osher-Fatemi total variation denoising using split Bregman,” IPOL 2012, 1–20 (2012). 11. J. Nocedal and S. Wright, Numerical Optimization, 2nd ed. (Springer, 2006). 12. M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging (Taylor & Francis, 1998). 13. H. Liao, F. Li, and M. K. Ng, “Selection of regularization parameter in total variation image restoration,” J. Opt. Soc. Am. A 26, 2311–2320 (2009).
14. T. Kreis, Holographic Interferometry: Principles and Methods (Wiley-VCH, 1996). 15. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20, 925–934 (2003).
16. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley-Interscience, 1998). 17. K. J. Gasvik, Optical Metrology, 3rd ed. (Wiley, 2002). 18. D. Strong and T. Chan, “Edge-preserving and scale-dependent properties of total variation regularization,” Inverse Probl. 19, S165–S187 (2003).
10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2301
We introduce a method based on the minimization of a total variation regularization cost function for computing discontinuous phase maps from fringe patterns. The performance of the method is demonstrated by numerical experiments with both synthetic and real data. © 2014 Optical Society of America OCIS codes: (120.2650) Fringe analysis; (120.3180) Interferometry; (120.3940) Metrology; (120.6650) Surface measurements, figure; (100.3190) Inverse problems. http://dx.doi.org/10.1364/AO.53.002297
1. Introduction
The main goal of fringe analysis techniques is to recover accurately the local modulated phase from one or several fringe patterns [1,2]; such phase is related to some physical quantities, such as shape, deformation, refractive index, and temperature. The basic model for a fringe pattern is given by I x ax bx cosφx ϕx ; where x x1 ; x2 , ax is the background illumination, bx is the amplitude modulation, and ϕx is the phase map to be recovered; the spatial carrier frequency of the fringe pattern is defined by the term φx . One challenging task is the recovery of a discontinuous phase map from a single fringe pattern. In recent years, many methods that successfully estimate the phase from a single pattern [3–6] have been reported. These methods, however, consider phase maps and illumination terms as continuous. Recently, a method was proposed for computing discontinuous phase maps, from a fringe pattern with carrier frequency, based on the minimization of a regularized cost function, which uses a second-order edge-preserving potential [7]. Although this method is able to detect and reconstruct phase discontinuities, its cost functional is 1559-128X/14/112297-05$15.00/0 © 2014 Optical Society of America
not convex; hence convergence to an optimal solution is conditioned to the provided initial phase usually computed by standard methods. In this work we present a method for recovering a discontinuous phase map from a single fringe pattern based on the total variational approach. First we describe the minimization problem based on the total variational approach, where we made some remarks about the convexity of the proposed method. Then we describe the numerical solution of the proposed cost function, which results in a simple algorithm. The resulting algorithm performance is tested by two numerical phase demodulations. 2. Demodulation of Fringe Pattern Using Total Variational Approach A. Variational Formulation
Following the image denoising method proposed by Rudin et al. [8] capable of preserving discontinuities, we propose to estimate the discontinuous phase map, background, and amplitude terms as the solution of the minimization problem defined by min Eϕx ; ax ; bx ≡
ϕx ;ax ;bx
Z Z λ I x − gx 2 dx j∇ϕx jdx 2 Ω Ω Z Z j∇ax jdx j∇bx jdx ; (1) Ω
Ω
10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2297
where gx is the noise input fringe pattern, Ω ⊂ R2 denotes the continuous signal domain, and λ is a positive parameter or Lagrange multiplier to control the contribution of the similarity term. Due to each total variation regularizer in Eq. (1) being convex [9], it is easy to verify that minimization with respect to ax yields the global minimum. Obtaining global minima for bx and ϕx is, on the other hand, not that simple. A numerical analysis carried out on Eq. (1) to study its convexity showed that provided the terms I x and gx are in the range 0; 2, the second-order derivatives with respect to bx and ϕx deliver positive semi-definite functionals. In Appendix A we give insight about the positive definiteness of each term in Eq. (1). In order to minimize Eq. (1), the first-order optimality condition or Euler–Lagrange equation has to be derived for each function ϕx , ax , and bx . In the formal derivation we assume that all these vector fields are smooth enough such that gradients are well defined and the variations have compact support over Ω so that we can use the divergence theorem to get rid of the boundary terms. By following standard methods from calculus of variations, the gradient descent equations can be obtained and are given by ∂ϕx ∇ϕx ∂I ∇· λI x − gx x ; ∂t j∇ϕx j ∂ϕx ∂ax ∇ax ∂I ∇· λI x − gx x ; ∂t j∇ax j ∂ax ∂bx ∇bx ∂I ∇· λI x − gx x ; ∂t j∇bx j ∂bx
(2)
with boundary conditions ∂ϕx 0; ∂ν
∂ax 0; ∂ν
and
∂bx 0; ∂ν
where ν denotes the outer normal to the boundary. B.
Numerical Solution
To solve Eq. (2) we apply a sequential method, i.e., we solve the first equation for ϕx, and the new value is immediately used in the second equation when solving for ax ; the same procedure is applied for bx, and this process is repeated iteratively until reaching convergence to the solution. We now proceed to show how to discretize the differential operator in the first equation of Eq. (2). The same procedure can be applied to the other two equations. For simplicity we drop the subscript x from the function ϕx and let ϕi;j ϕx1 ; x2 to denote the value of a grid function ϕ at point i; j x1 ; x2 defined on the cell-centered grid Ωh equal to all points x1 ; x2 ∈ Ω such that 2i − 1hx1 2j − 1hx2 x1 ; x2 ; 2 2 with 1 ≤ i ≤ m, 1 ≤ j ≤ n. This way Ωh consists of m × n cells of size hx1 × hx2 with hx1 1∕m, hx2 1∕n the grid spacing. 2298
APPLIED OPTICS / Vol. 53, No. 11 / 10 April 2014
To approximate the derivatives, we use forward and backward finite differences defined as follows: δ x1 ϕi;j
ϕi1;j − ϕi;j h x1
and δ x2 ϕi;j
ϕi;j1 − ϕi;j : hx2
The divergence term is approximated using ∇ · V i;j δ−x1 V 1i;j δ−x2 V 2i;j ; where V i;j V 1i;j ; V 2i;j
∇ϕi;j ; j∇ϕji;j
∇ϕi;j δ x1 ϕi;j ; δx2 ϕi;j ; q 2 2 j∇ϕji;j δ x1 ϕi;j δx2 ϕi;j :
The rest of the terms in the equation are implemented by straightforward evaluation at point i; j. Finally, the Neumann’s boundary condition on ∂Ω is treated as ϕi;0 ϕi;1 ;
ϕi;n1 ϕi;n ;
ϕ0;j ϕ1;j ;
ϕm1;j ϕm;j ;
ϕ0;0 ϕ1;1 ;
ϕ0;n1 ϕ1;n ;
ϕm1;0 ϕm;1 ;
ϕm1;n1 ϕm;n :
Examples of numerical implementations similar to Eq. (2) can be found in Refs. [8] and [10]. 3. Numerical Experiments
To illustrate the performance of the proposed method, we carried out two numerical experiments using an Intel Core i5 @ 2.40 GHz laptop with Debian (wheezy) 64-bit and 4 GB of memory. The first experiment was the demodulation of a synthetic fringe pattern of size 250 × 250 pixels generated in similar way to that described in Ref. [7]: 8 a 0.6; bx 0.5; if 90 ≤ x1 ; x2 ≤ 160 > > > x > < ϕ 1.2 0.005x x Ix ; > ax 1.0; bx 0.9; otherwise > > > : ϕx 0.2 0.005x φx
2π x; P
(3)
where x1 ; x2 and P 12.5 are given in pixels. Figure 1(a) shows the synthetic fringe pattern, and Fig. 1(b) shows the synthetic phase term ϕx . The initial values for the demodulation process were ax bx 1.0, ϕx 0.2. The resultant phasemap estimation using the proposed method is shown in Fig. 2(a). The mean squared error (MSE) of the
Fig. 3. (a) Experimental fringe pattern. (b) Demodulated phase term. Fig. 1. (a) Synthetic fringe pattern. (b) Synthetic phase term ϕx .
demodulation was 0.00285, and the time employed to get the solution was 144 s using a simple gradient descent scheme with λ 100 [11–13]. For purpose of comparison, we implemented the method described in Ref. [7]. The resultant phase-map estimation is shown in Fig. 2(b). The MSE of the demodulation was 0.06722, and the time employed by this method was 800 s using the same numerical scheme. As can be observed, the estimation based on Eq. (2) successfully demodulates the fringe pattern with phase discontinuities. Because of the simple definition of the term ϕx , Eq. (3), the proposed cost function will converge to an accurate solution given an arbitrary initial point. On the contrary, the time used by the method reported in Ref. [7] was larger than the time used by our method, Eq. (2), and fails to demodulate the fringe pattern using initial points far away of the optimal point. The second experiment was the processing of a fringe pattern obtained from a holographic interferometry experiment [14], which consisted of the height measurement of a microthin film. The fringe pattern obtained from this experiment, with 640 × 480 pixels, is shown in Fig. 3(a). The fringe pattern was demodulated using Fourier-based quadrature transform [15]; the resultant demodulated phase term is shown in Fig. 3(b). During the demodulation process, the background, modulation, and spatial carrier terms are estimated; then the carrier term is subtracted from the wrapped phase term and the phase term is unwrapped if necessary [16]. These estimations are used as initial values in our method.
Fig. 2. Estimated phase terms using (a) Eq. (2) and (b) method from Ref. [7].
Fig. 4. Estimated phase terms using (a) Eq. (2) and (b) Schwider– Hariharan 41 algorithm [17].
The resultant phase-map estimation is shown in Figs. 4(a) and 5. The time employed for the proposed method was 210 s using a simple gradient descent scheme with λ 100 [11–13]. For purposes of comparison, we show in Figs. 4(b) and 6 the resultant phase-map estimation using a 41 phase shifting method [17].
Fig. 5. Estimated phase terms using Eq. (2).
Fig. 6. Estimated phase terms using Schwider–Hariharan 41 algorithm [17]. 10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2299
The second-order derivatives of F with respect to ax , bx , and ϕx read as follows: d2 F 2 da2x d2 F 2 db2x d2 F 2 dϕ2x Fig. 7. Columns of the estimated phase terms using Eq. (2) and Schwider–Hariharan 4 1 algorithm [17].
As can be observed in Figs. 4 and 7, the estimated phase map obtained by the solution of Eq. (2) is very similar to the one obtained using a well-proven demodulation method. An important point to highlight from the result shown in Fig. 7 is the preservation of the dynamic range of the estimated phase map. In Ref. [18] the authors showed that the change in image intensity (here the dynamic range of the phase) due to total variation regularization is directly proportional to the regularization parameter and inversely proportional to the scale of the image feature; that means total variation regularization causes that smaller-scaled changes in phase be partially or entirely removed while larger-scaled ones be relatively unaffected. The resultant phase map of our method can be seen as the filtered phase term of the phase shifting method, as one can observe in Figs. 5 and 6. Another important aspect of this experiment is the need for using a prior estimate for demodulating real fringe patterns because the optimization process could be trapped on a local minimum due to large phase variations. 4. Conclusion
As can be seen from the resultant phase estimations in the above two experiments, the proposed cost function allows the accurate demodulation of a single fringe pattern with discontinuities. The numerical solution of Eq. (2) results in a very simple algorithm that can be implemented with several computational efficient techniques [10]. An extra advantage of the proposed method is its feasibility to be implemented on dedicated hardware to obtain real-time processing. This will be one aim of our future research. Appendix A
Let us define F as Z λ I − gx 2 dx: F 2 Ω x Z ax bx cosφx ϕx − gx 2 dx: Ω
2300
APPLIED OPTICS / Vol. 53, No. 11 / 10 April 2014
Z Ω
Z Ω
dx > 0;
(A1)
cos2 φx ϕx dx ≥ 0;
(A2)
Z Ω
gx − ax cosφx ϕx
− bx cos2φx ϕx dx.
(A3)
From Eq. (A1), it is clearR that F is convex with respect to ax , and due to Ω j∇ax jdx being convex as well, as proven in [9], minimization with respect to ax yields the global minimum. In a similar fashion, Eq. (A2) reveals that F is at most a semi-positive definite functional. Lately, from Eq. (A3), not much can be said about F with respect to ϕx . However, through a large number of numerical tests we noticed that provided I x and gx are in the range 0; 2, Eq. (A3) is not negative and therefore F is a semi-positive definite functional. Although this observation is by no means conclusive, it is a good practice to set I x and gx to this interval. References 1. D. W. Robinson and G. T. Reid, Interferogram Analysis, Digital Fringe Pattern Measurement Techniques (Taylor & Francis, 1993). 2. G. Rajshekhar and P. Rastogi, “Fringe analysis: premise and perspectives,” Opt. Lasers Eng. 50, iii–x (2012). 3. J. L. Marroquin, M. Rivera, S. Botello, R. Rodriguez-Vera, and M. Servin, “Regularization methods for processing fringepattern images,” Appl. Opt. 38, 788–794 (1999). 4. J. Villa, J. A. Quiroga, and M. Servin, “Improved regularized phase-tracking technique for the processing of squared-grating deflectograms,” Appl. Opt. 39, 502–508 (2000). 5. R. Legarda-Saenz, W. Osten, and W. P. Juptner, “Improvement of the regularized phase tracking technique for the processing of nonnormalized fringe patterns,” Appl. Opt. 41, 5519–5526 (2002). 6. M. Rivera, “Robust phase demodulation of interferograms with open or closed fringes,” J. Opt. Soc. Am. A 22, 1170– 1175 (2005). 7. C. Galvan and M. Rivera, “Second-order robust regularization cost function for detecting and reconstructing phase discontinuities,” Appl. Opt. 45, 353–359 (2006). 8. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). 9. R. Acar and C. R. Vogel, “Analysis of bounded variation penalty methods for ill-posed problems,” Inverse Probl. 10, 1217–1229 (1994). 10. P. Getreuer, “Rudin-Osher-Fatemi total variation denoising using split Bregman,” IPOL 2012, 1–20 (2012). 11. J. Nocedal and S. Wright, Numerical Optimization, 2nd ed. (Springer, 2006). 12. M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging (Taylor & Francis, 1998). 13. H. Liao, F. Li, and M. K. Ng, “Selection of regularization parameter in total variation image restoration,” J. Opt. Soc. Am. A 26, 2311–2320 (2009).
14. T. Kreis, Holographic Interferometry: Principles and Methods (Wiley-VCH, 1996). 15. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20, 925–934 (2003).
16. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley-Interscience, 1998). 17. K. J. Gasvik, Optical Metrology, 3rd ed. (Wiley, 2002). 18. D. Strong and T. Chan, “Edge-preserving and scale-dependent properties of total variation regularization,” Inverse Probl. 19, S165–S187 (2003).
10 April 2014 / Vol. 53, No. 11 / APPLIED OPTICS
2301
