Single shot interferogram analysis for optical metrology Mahendra Pratap Singh,1,2 Mandeep Singh,1 and Kedar Khare1,* 1
Department of Physics, Indian Institute of Technology Delhi, Hauz Khas, New Delhi 110016, India 2
Instruments Research and Development Establishment, Dehradun 248008, India *Corresponding author: [email protected] Received 12 May 2014; revised 28 August 2014; accepted 28 August 2014; posted 3 September 2014 (Doc. ID 211927); published 3 October 2014
We report a novel constrained optimization method for single shot interferogram analysis. The unknown test wavefront is estimated as a minimum L2-norm squared solution whose phase is constrained to the space spanned by a finite number of Zernike polynomials. Using a single frame from standard phase shifting datasets, we demonstrate that our approach provides a phase map that matches with that generated using phase shifting algorithms to within λ∕100 rms error. Our simulations and experimental results suggest the possibility of a simplified low-cost high quality optical metrology system for performing routine metrology tests involving smooth surface profiles. © 2014 Optical Society of America OCIS codes: (120.3940) Metrology; (120.3180) Interferometry; (120.2650) Fringe analysis; (100.3190) Inverse problems. http://dx.doi.org/10.1364/AO.53.006713
1. Introduction
Interferometry is a widely used tool for many applications such as optical testing and quality control. Typical interferometric methods for optical metrology currently use an in-line Fizeau configuration where the wavefront reflected from an unknown surface interferes with the wavefront reflected from a null surface [1,2] and the measurement of the unknown optic profile can be performed by estimating the phase from the interference data. Several approaches for phase estimation exist in the literature of which the most commonly used approaches are the phase shifting and the Fourier transform methods. In the phase shifting method [3], multiple interferogram frames with known phase shifts are recorded and are then processed digitally for phase estimation. While the resolution of phase map obtained with the phase shifting method is typically same as the digital array detector used for recording the 1559-128X/14/296713-06$15.00/0 © 2014 Optical Society of America
interferograms, the requirement of multiple frames implies that the phase shifting method is not suitable when the metrology system is under the influence of external disturbances and vibrations. A single shot approach for interferometry is most desirable in harsh environments, e.g., real-time metrology of a part under manufacturing/polishing process, since a single short exposure interferogram frame can freeze the effects of vibrations and other ambient disturbances. The phase estimation from a single closed-fringe interferogram is, however, a difficult problem in general. Historically, the most prominent single shot technique is the Fourier transform method (FTM) [4–6] which requires introduction of carrier frequency in the interference pattern. The carrier-frequency structure of the interference pattern, however, imposes a resolution limitation in that the achievable wavefront resolution is typically lower than the resolution capability of the digital detector array used for recording the interferogram. Additionally, FTM cannot be applied to interferograms containing closed fringes, as there is no clear separation of interference terms in the Fourier 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6713
domain. Another important class of single shot methods is based on the regularized phase tracking (RPT) method [7]. The RPT method requires that the lowfrequency dc part of the interference pattern is removed and, further, it aims at estimating the phase map as a solution to a constrained optimization problem. The considerations of data fidelity and a constraint of locally linear phase model are used in the optimization process. Several improvements over the original RPT technique have been considered in the literature; for example, fringe follower phase tracking [8], phase tracking with Bayesian estimation theory with prior Markov random field model [9], and phase tracking without normalization [10]. RPT-based methods are powerful tools for fringe pattern analysis except for the sign ambiguity in estimated phase map that is inherent in the optimization model and the computational complexity due to the requirement of solving local model fitting for every pixel. Fringe demodulation problem via multidimensional quadrature transforms has also received some attention [11,12]. A sequential quadrature and phase tracker has also been developed [13]. The quadrature based methods typically require some form of local frequency and orientation estimate as an additional step. Tracking of arccosine values of a processed interferogram pattern guided by the local frequency estimates has also been demonstrated for single shot interferogram analysis [14]. In this paper, we present a method to analyze a single interferogram that is based on our recent work on single shot high resolution digital holographic imaging [15,16]. In this work we demonstrate a constrained optimization approach for the demodulation of a hologram that can recover high resolution information even when the dc and the cross terms in the hologram overlap with each other. While our experiments in [15,16] used off-axis configuration, we have suitably modified our approach for the metrology problem where interference data often consists of closed fringes. While the approach presented here is similar in spirit to the RPT method, it differs from RPT in details of the modelling. The paper is organized as follows. In Section 2, we present our new approach to phase estimation from a single interferogram. In Section 3, simulation as well as experimental results obtained using this approach are presented. In Section 4 we provide simulation results on the phase estimation error in the presence of artificially introduced noise. Finally, in Section 5 we provide concluding remarks. 2. Optimization Approach to Phase Estimation from Single Shot Interferogram
The recorded intensity pattern for interference of an unknown test wavefront, W W 0 expiΦw , and a known reference beam, R R0 expiΦR , can be represented as H data jWj2 jRj2 WR W R: 6714
(1)
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
All the symbols in the above equation denote functions of spatial variables x; y in the interferogram plane. The knowledge of the phase of the reference beam is assumed here as it represents the prior knowledge about the null surface. Also, the amplitude, jRj, can be recorded easily as a calibration step, if required. We model the problem of determining the phase, Φw , of the unknown test wavefront as a constrained optimization problem as follows: Minimize CW; W ‖H data − jWj2 jRj2 WR W R‖22 with the constraint that
Φw
N X
an Z n :
(2)
n1
Here, the notation ‖…‖22 represents the L2-norm squared error, Zn represents the nth Zernike polynomial ordered as per the convention in [17], an represents the coefficient of expansion, and N represents the highest index of the Zernike polynomial to be used for a given problem at hand. The minimization of the L2-norm squared error term results in a solution that fits in the least-square sense with the measured data and the constraint effectively restricts the unknown phase map to the space spanned by a finite number, N, of Zernike terms. We observe that the RPT method [7] uses a smoothness constraint on the estimated phase map. One major difference between the method described here and the RPT technique is that we do not remove the constant or dc terms from the interferogram as a preprocessing step. The update step, however, involves the complex valued reference beam, R, so that the appropriate complex valued wavefront, Wx; y, is obtained as a solution to the optimization problem. Further, the phase function in this work is not assumed to be locally linear or quadratic, as in the RPT method, but the smoothness is modeled globally in terms of Zernike polynomials. In our recent work on digital holography [15,16] where Fresnel zone hologram configurations were used, a smoothness constraint was utilized as this constraint is suitable for a diffracted field. We remark that the restriction of the phase map to be a linear combination of a finite number of Zernike terms also implicitly enforces a smoothness constraint. Further, Zernike expansion is convenient as it is a standard method used in metrology. We solve the above constrained optimization problem using an iterative gradient descent scheme [18] combined with a projection operator acting on the phase function. We note that the cost function, CW; W , is a function of both W and W . Assuming that the samples of the wavefront function, Wx; y, are written in the form of a column vector, the variation of CW; W may be represented as
δC ∇W CT δW ∇W CT δW ;
the minimum norm solution above is straightforward. In particular, the coefficients may be formally represented as
2 Re∇W CT δW; 2 Re∇W CH δW:
(3)
The superscripts, T and H, above stand for the transpose and Hermitian conjugate operations, respectively, and the last step follows from the fact that CW; W is a real valued function, so that ∇W C ∇W C . From the Cauchy–Schwartz inequality we observe that the direction of steepest descent is achieved if δW is in the same direction as ∇W C [19], as given by ∇W CW;W −H data − jWj2 jRj2 WR W R · W R: (4) Combining the L2-squared error reduction and the Zernike expansion constraint gives an iterative algorithm as follows, U n1 W n − t∇W CWW n W n1 PZ U n1 ; (5) Here, U n1 represents the L2-norm reducing update. The step size, ‘t’, is selected in each iteration by a standard backtracking line search [20]. Typically, one may start with a particular value of ‘t’ (say t 1) and evaluate the cost function, CW; W , associated with the guess solution W U n1 . If required, the value of t may be reduced in some regular manner until, for some t, the cost function has a lower value as compared to CW; W evaluated at W W n . The symbol, PZ , in the above equation represents projection of the phase of U n1 onto the space of the N Zernike polynomials. In particular, application of PZ on a complex valued image involves three steps: 1. Find the argument of the image pixel by pixel to get a phase map, Φ0w argU n1 :
(6)
2. Unwrap the phase map so obtained, Φw unwrapΦ0w :
(7)
3. Find the best L2-norm fit for the unwrapped phase map in terms of the first N Zernike polynomials. This step finds a set of coefficients, an , with n 1; 2; 3; …; N such that the L2-norm error, represented as ‖Φw −
N X n1
a n Zn ‖ 2 ;
is minimized. Since Zernike polynomials form an orthogonal set of functions over a circle, finding
an Zn ; Φw ∕Zn ; Zn ;
(8)
where (,) denotes the scalar product, ZZ f ; g dxdyf x; y gx; y:
(9)
The integration above is performed over the circular aperture used for defining the Zernike polynomials. We note that the operation represented by PZ is applied to the phase part of U n1, as in Eq. (5). This updated phase solution along with amplitude of U n1 constitutes the next guess solution W n1 . We observe that while the L2-norm reducing iteration decreases the cost function, CW; W , the action of the operator, PZ , will tend to increase its numerical value. When these two parts of the iteration act in opposing directions, the overall solution shows very small overall change. The iteration may, therefore, be stopped using some threshold, ε, on the change in the solution in successive iterations, ‖W n1 − W n ‖2 < ε: ‖W n ‖2
(10)
While phase unwrapping may be considered an extra step, several standard methods, such as the Goldstein method [21], are available that may be utilized for this purpose. We observe that the two parts of the iterative procedure result in a solution that is smooth to the extent as represented by Zernike expansion and also matches the interference data. While we have not given any convergence proofs that are out of the scope of this paper, we believe that representing the phase part of the wavefront as a linear combination of Zernike polynomials is a strong constraint that is helpful in restricting the search space for the solution and avoiding problems such as stagnation due to local minima. Based on the form of the constraint that restricts the phase function to be a linear combination of finite number, N, of Zernike polynomials, it is clear that the present method is limited to smooth surfaces and will not include cases involving high frequency features such as phase steps. We believe that in order to address the phase estimation problem for such high frequency cases, simply adding number of Zernike terms may not be the best solution. Instead, a modification of the constraint will be required to include Total Variation or L1-norm based penalties. We will not include these cases in the present paper, but report on such generalization of the constraint in the future. 3. Numerical Simulations and Experimental Results
In this section we describe our simulations and experiments in order to test the iterative phase 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6715
estimation method described in Section 2. For simulations, we assumed a test wavefront in the form of a polynomial given by Wx; y expfikβ1 xy β2 x2 y2 β3 x3 y3 g: (11) The reference beam (from null surface) was assumed to be a quadratic wavefront given by Rx; y expikαx2 y2 :
(12)
The numerical values of parameters in Eqs. (11) and (12) are: β1 6.2λ, β2 8λ, β3 2.2λ, α 4λ, and k denotes the wave number, 2π∕λ. The interferogram and the simulated test surface profile are shown in Figs. 1(a) and 1(b), respectively. Figures 1(c) and 1(d) show the recovered phase map and the interferogram generated using the recovered phase map, respectively. The result was obtained in four iterations starting with a zero image as the first guess solution. The number of Zernike terms used for modeling is N 36. The reason for using 36 Zernike polynomials is that this set provides radial polynomials up to the 8th degree and angular terms up to 7θ. In most routine applications of interferometric surface testing this is the typical order that is used to describe surface profiles [22,23]. The relative squared L2-norm error, as in Eq. (2), is plotted in Fig. 2(a) and is seen to reduce to 0.7%. The error between the test phase [Eq. (11)] and the recovered phase is shown in Fig. 2(b). We observe that the rms error between the test and recovered wavefronts is 0.003 waves after four iterations. The relative change in the solution between successive iterations was observed to
Fig. 1. Simulation for single shot interferogram analysis using the proposed iterative method. (a) Interferogram, (b) test phase map (vertical axis unit is waves), (c) recovered phase map from single interferogram frame, (d) interferogram generated using the recovered phase map. 6716
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
Fig. 2. Recovery error for proposed iteration procedure. (a) Relative L2-norm error between the interferogram data and interferogram model [see Eq. (2)], (b) error map between the known test phase and recovered phase map.
be less than 0.1% at four iterations. Our preliminary code written in MATLAB requires less than 30 s to complete 4–5 iterations on a laptop computer with a 1.5 GHz processor. Approximately 80% of this time is taken by the phase unwrapping algorithm. The computational efficiency of the method can be improved if we can reduce the time taken for the phase unwrapping process. This can be done by using a faster algorithm for phase unwrapping. In order to test our algorithm further we used experimental datasets from a commercial 9-step phase shifting Fizeau interferometer shown schematically in Fig. 3. A number of datasets were recorded that correspond to various types of imaging lenses, and other optical elements involving curved and flat surfaces. Spherical reference surfaces were used to test the curved surfaces in the Fizeau interferometer setup. The performance of the optimization approach is comparable in all the cases. As an illustration, we
Fig. 3. Schematic diagram of Fizeau interferometer setup.
show one of the datasets for a curved surface, in this paper, that would typically be tested using a multistep phase shifting approach. Nine interferograms are recorded with relative phase shifts of π∕2 between consecutive frames that are achieved using a calibrated PZT stepper. The nine frames may be used to recover the phase map from the interferogram using the phase shifting formula [24], −4I 2 12I 4 − 12I 6 4I 8 : arctan I 1 − 8I 3 14I 5 − 8I 7 I 9 (13)
Φtest − Φnull
Here, the intensity patterns corresponding to the nine interferograms are denoted by symbols I n with n 1; 2; …; 9. With the prior knowledge of the null wavefront phase one may recover the unknown test object phase. The null wavefront, if not known, can also be assumed to have tilt and quadratic terms that can be estimated approximately from the fringe pattern center and the number of closed fringes. The final estimate of the test wavefront will then be with respect to this choice of null surface. The first of the nine interferogram frames was also used in our iterative procedure in order to determine the test wavefront phase. Figures 4(a)–4(c) show three of the nine phase shifted interferograms. The recovered test object phase maps using the phase shifting method (with 9 frames) and the proposed iterative method (with single frame) are shown in Figs. 5(a) and 5(b), respectively. Figures 5(c) and 5(d) show 3D surface maps corresponding to phase maps shown in Figs. 5(a) and 5(b), respectively. From numerical comparison we observe that the rms error between the phase maps estimated from the 9-step phase shifting method and the single step iterative method is equal to 0.008 waves showing excellent agreement between the two methods.
Fig. 5. Comparison of single shot iterative method with experimental data from a phase shifting Fizeau interferometer. (a) Phase map obtained using 9-step phase shifting method and (b) phase map obtained using a single frame out of the phase shifting data. (c), (d) 3D surface plots corresponding to phase maps in (a), (b) respectively. The rms error between the two recoveries is 0.008 waves.
approach. One of the four interferograms (corresponding to zero reference phase) was then used for phase estimation using the iterative procedure described in Section 2. Figure 6(a) shows one of the four interferograms for a simulated test object with artificial noise. Figures 6(b) and 6(c) show the phase maps estimated using the phase shifting method and the iterative optimization procedure,
4. Simulation Results with Artificially Introduced Noise
In order to understand the noise sensitivity of our proposed iterative method we used the simulations with test surface as in Section 3, Fig. 1. Uniform random noise patterns in the range (−0.5, 0.5) radians were added to the simulated test wavefront to generate four phase shifted interferograms (corresponding to reference phase shift of 0, π∕2, π, 3π∕2). The four interferograms were first used to determine the test object phase using the standard phase shifting
Fig. 4. Three of the nine phase shifted interferogram frames recorded on the Fizeau interferometer setup shown in Fig. 3.
Fig. 6. Simulations with artificially introduced noise in range [−0.5, 0.5] radians. (a) One of the four simulated noisy interferograms; (b), (c) are estimated phase maps obtained from the noisy interferogram with the phase shifting approach, and proposed iterative approach; (d), (e) are errors between the known ground truth and recovered phase maps using the phase shifting and the iterative method. [Colorbar is in waves with values in range (−0.02, 0.018).] 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6717
respectively. Figures 6(d) and 6(e) show the error between the known ground truth test object and the phase maps obtained by the phase shifting and optimization methods, respectively. We observe that the rms phase errors are similar in two cases and approximately equal to 0.002 waves. The error map for the solution obtained using the optimization method is much smoother, as expected, due to the inherent smoothness property associated with this method. This illustration shows that the single-shot iterative optimization method has similar rms noise performance compared to the multishot phase shifting method. 5. Conclusion
We have demonstrated an optimization approach to phase estimation from a single interferogram frame. Both simulation and experimental results were described and they show excellent phase map recovery for this single-shot approach. The performance of this iterative method is observed to be within λ∕100 rms error with respect to the standard multishot phase shifting approaches, as we have shown with experimental data and simulated noisy data. The method, if implemented with dedicated processing hardware, may provide real time feedback on optical metrology to the polishing machine during the manufacturing process. The proposed method is general enough and can be used in any standard interferometer configuration without much difficulty. We believe that the proposed approach can enable low cost high quality optical metrology systems that are simple to set up. The phase estimation method as presented in this work is limited to smooth surfaces and we will report on possible generalization of the constraint term to include high frequency features, such as steps, in future work. Financial support in order to pursue this allied work is acknowledged from the Department of Biotechnology, Ministry of Science and Technology of India (BT/PR8008/MED/32/306/2013). The authors are also grateful to Dr. A. K. Gupta, Director, IRDE Dehradun, for allowing the use of the experimental facility for data collection. References 1. D. Malacara, M. Servin, and Z. Malacara, Interferogram Analysis for Optical Testing, 2nd ed. (CRC Press, 2005). 2. J. Bruning, D. Herriot, J. Gallagher, D. Rosenfeld, A. White, and D. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt 13, 2693–2703 (1974). 3. K. Creath, “Phase shifting interferometry techniques,” in Progress in Optics, E. Wolf, ed. (Elsevier, 1988), Vol. 26, pp. 357–373.
6718
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
4. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. A 72, 156–160 (1982). 5. T. Kreis, “Digital holographic interference-phase measurement using the Fourier transform method,” J. Opt. Soc. Am. A 3, 847–855 (1986). 6. M. Takeda, “Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: an overview,” Ind. Metrol. 1, 79–99 (1990). 7. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36, 4540–4548 (1997). 8. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Fringe-follower regularized phase tracker for demodulation of closed-fringe interferograms,” J. Opt. Soc. Am. A 18, 689–695 (2001). 9. J. L. Marroquin, L. Rivera, S. Botello, R. Rodriguez-Vera, and M. Servin, “Regularization methods for processing fringe pattern images,” Appl. Opt. 38, 788–794 (1999). 10. R. Legarda-Saenz, W. Osten, and W. Juptner, “Improvement of the regularized phase tracking technique for the processing of nonnormalized fringe patterns,” Appl. Opt. 41, 5519–5526 (2002). 11. 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). 12. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General ndimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20, 925– 934 (2003). 13. M. Servin, J. L. Marroquín, and J. A. Quiroga, “Regularized quadrature and phase tracking from a single closed-fringe interferogram,” J. Opt. Soc. Am. A 21, 411–419 (2004). 14. Q. Kemao and S. H. Soon, “Sequential demodulation of a single fringe pattern guided by local frequencies,” Opt. Lett. 32, 127–129 (2007). 15. K. Khare, P. T. Samsheer Ali, and J. Joseph, “Single shot high resolution digital holography,” Opt. Express 21, 2581–2591 (2013). 16. P. T. Samsheerali, K. Khare, and J. Joseph, “Quantitative phase imaging using single shot digital holography,” Opt. Commun. 319, 85–89 (2014). 17. J. Wyant and K. Creath, “Basic wavefront aberration theory for optical metrology,” in Applied Optics and Optical Engineering, R. R. Shannon and J. Wyant, eds. (Academic, 1992), Vol. XI, pp. 1–53. 18. M. Bertero and P. Boccacchi, Introduction to Inverse Problems in Imaging (IOP, 1998). 19. D. H. Brandwood, “A complex gradient operator and its application in adaptive array theory,” Proc. IRE 130, 11–16 (1983). 20. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University, 2004). 21. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping (Wiley, 1998). 22. C. J. Evans, “PVr—a robust amplitude parameter for optical surface specification,” Opt. Eng. 48, 043605 (2009). 23. C. J. Evans, R. E. Parks, P. J. Sullivan, and J. S. Taylor, “Visualization of surface figure by the use of Zernike polynomials,” Appl. Opt. 34, 7815–7819 (1995). 24. K. Creath, “Choosing a phase-measurement algorithm for measurement of coated LIGO optics,” Proc. SPIE 4101, 46– 55 (2000).
Department of Physics, Indian Institute of Technology Delhi, Hauz Khas, New Delhi 110016, India 2
Instruments Research and Development Establishment, Dehradun 248008, India *Corresponding author: [email protected] Received 12 May 2014; revised 28 August 2014; accepted 28 August 2014; posted 3 September 2014 (Doc. ID 211927); published 3 October 2014
We report a novel constrained optimization method for single shot interferogram analysis. The unknown test wavefront is estimated as a minimum L2-norm squared solution whose phase is constrained to the space spanned by a finite number of Zernike polynomials. Using a single frame from standard phase shifting datasets, we demonstrate that our approach provides a phase map that matches with that generated using phase shifting algorithms to within λ∕100 rms error. Our simulations and experimental results suggest the possibility of a simplified low-cost high quality optical metrology system for performing routine metrology tests involving smooth surface profiles. © 2014 Optical Society of America OCIS codes: (120.3940) Metrology; (120.3180) Interferometry; (120.2650) Fringe analysis; (100.3190) Inverse problems. http://dx.doi.org/10.1364/AO.53.006713
1. Introduction
Interferometry is a widely used tool for many applications such as optical testing and quality control. Typical interferometric methods for optical metrology currently use an in-line Fizeau configuration where the wavefront reflected from an unknown surface interferes with the wavefront reflected from a null surface [1,2] and the measurement of the unknown optic profile can be performed by estimating the phase from the interference data. Several approaches for phase estimation exist in the literature of which the most commonly used approaches are the phase shifting and the Fourier transform methods. In the phase shifting method [3], multiple interferogram frames with known phase shifts are recorded and are then processed digitally for phase estimation. While the resolution of phase map obtained with the phase shifting method is typically same as the digital array detector used for recording the 1559-128X/14/296713-06$15.00/0 © 2014 Optical Society of America
interferograms, the requirement of multiple frames implies that the phase shifting method is not suitable when the metrology system is under the influence of external disturbances and vibrations. A single shot approach for interferometry is most desirable in harsh environments, e.g., real-time metrology of a part under manufacturing/polishing process, since a single short exposure interferogram frame can freeze the effects of vibrations and other ambient disturbances. The phase estimation from a single closed-fringe interferogram is, however, a difficult problem in general. Historically, the most prominent single shot technique is the Fourier transform method (FTM) [4–6] which requires introduction of carrier frequency in the interference pattern. The carrier-frequency structure of the interference pattern, however, imposes a resolution limitation in that the achievable wavefront resolution is typically lower than the resolution capability of the digital detector array used for recording the interferogram. Additionally, FTM cannot be applied to interferograms containing closed fringes, as there is no clear separation of interference terms in the Fourier 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6713
domain. Another important class of single shot methods is based on the regularized phase tracking (RPT) method [7]. The RPT method requires that the lowfrequency dc part of the interference pattern is removed and, further, it aims at estimating the phase map as a solution to a constrained optimization problem. The considerations of data fidelity and a constraint of locally linear phase model are used in the optimization process. Several improvements over the original RPT technique have been considered in the literature; for example, fringe follower phase tracking [8], phase tracking with Bayesian estimation theory with prior Markov random field model [9], and phase tracking without normalization [10]. RPT-based methods are powerful tools for fringe pattern analysis except for the sign ambiguity in estimated phase map that is inherent in the optimization model and the computational complexity due to the requirement of solving local model fitting for every pixel. Fringe demodulation problem via multidimensional quadrature transforms has also received some attention [11,12]. A sequential quadrature and phase tracker has also been developed [13]. The quadrature based methods typically require some form of local frequency and orientation estimate as an additional step. Tracking of arccosine values of a processed interferogram pattern guided by the local frequency estimates has also been demonstrated for single shot interferogram analysis [14]. In this paper, we present a method to analyze a single interferogram that is based on our recent work on single shot high resolution digital holographic imaging [15,16]. In this work we demonstrate a constrained optimization approach for the demodulation of a hologram that can recover high resolution information even when the dc and the cross terms in the hologram overlap with each other. While our experiments in [15,16] used off-axis configuration, we have suitably modified our approach for the metrology problem where interference data often consists of closed fringes. While the approach presented here is similar in spirit to the RPT method, it differs from RPT in details of the modelling. The paper is organized as follows. In Section 2, we present our new approach to phase estimation from a single interferogram. In Section 3, simulation as well as experimental results obtained using this approach are presented. In Section 4 we provide simulation results on the phase estimation error in the presence of artificially introduced noise. Finally, in Section 5 we provide concluding remarks. 2. Optimization Approach to Phase Estimation from Single Shot Interferogram
The recorded intensity pattern for interference of an unknown test wavefront, W W 0 expiΦw , and a known reference beam, R R0 expiΦR , can be represented as H data jWj2 jRj2 WR W R: 6714
(1)
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
All the symbols in the above equation denote functions of spatial variables x; y in the interferogram plane. The knowledge of the phase of the reference beam is assumed here as it represents the prior knowledge about the null surface. Also, the amplitude, jRj, can be recorded easily as a calibration step, if required. We model the problem of determining the phase, Φw , of the unknown test wavefront as a constrained optimization problem as follows: Minimize CW; W ‖H data − jWj2 jRj2 WR W R‖22 with the constraint that
Φw
N X
an Z n :
(2)
n1
Here, the notation ‖…‖22 represents the L2-norm squared error, Zn represents the nth Zernike polynomial ordered as per the convention in [17], an represents the coefficient of expansion, and N represents the highest index of the Zernike polynomial to be used for a given problem at hand. The minimization of the L2-norm squared error term results in a solution that fits in the least-square sense with the measured data and the constraint effectively restricts the unknown phase map to the space spanned by a finite number, N, of Zernike terms. We observe that the RPT method [7] uses a smoothness constraint on the estimated phase map. One major difference between the method described here and the RPT technique is that we do not remove the constant or dc terms from the interferogram as a preprocessing step. The update step, however, involves the complex valued reference beam, R, so that the appropriate complex valued wavefront, Wx; y, is obtained as a solution to the optimization problem. Further, the phase function in this work is not assumed to be locally linear or quadratic, as in the RPT method, but the smoothness is modeled globally in terms of Zernike polynomials. In our recent work on digital holography [15,16] where Fresnel zone hologram configurations were used, a smoothness constraint was utilized as this constraint is suitable for a diffracted field. We remark that the restriction of the phase map to be a linear combination of a finite number of Zernike terms also implicitly enforces a smoothness constraint. Further, Zernike expansion is convenient as it is a standard method used in metrology. We solve the above constrained optimization problem using an iterative gradient descent scheme [18] combined with a projection operator acting on the phase function. We note that the cost function, CW; W , is a function of both W and W . Assuming that the samples of the wavefront function, Wx; y, are written in the form of a column vector, the variation of CW; W may be represented as
δC ∇W CT δW ∇W CT δW ;
the minimum norm solution above is straightforward. In particular, the coefficients may be formally represented as
2 Re∇W CT δW; 2 Re∇W CH δW:
(3)
The superscripts, T and H, above stand for the transpose and Hermitian conjugate operations, respectively, and the last step follows from the fact that CW; W is a real valued function, so that ∇W C ∇W C . From the Cauchy–Schwartz inequality we observe that the direction of steepest descent is achieved if δW is in the same direction as ∇W C [19], as given by ∇W CW;W −H data − jWj2 jRj2 WR W R · W R: (4) Combining the L2-squared error reduction and the Zernike expansion constraint gives an iterative algorithm as follows, U n1 W n − t∇W CWW n W n1 PZ U n1 ; (5) Here, U n1 represents the L2-norm reducing update. The step size, ‘t’, is selected in each iteration by a standard backtracking line search [20]. Typically, one may start with a particular value of ‘t’ (say t 1) and evaluate the cost function, CW; W , associated with the guess solution W U n1 . If required, the value of t may be reduced in some regular manner until, for some t, the cost function has a lower value as compared to CW; W evaluated at W W n . The symbol, PZ , in the above equation represents projection of the phase of U n1 onto the space of the N Zernike polynomials. In particular, application of PZ on a complex valued image involves three steps: 1. Find the argument of the image pixel by pixel to get a phase map, Φ0w argU n1 :
(6)
2. Unwrap the phase map so obtained, Φw unwrapΦ0w :
(7)
3. Find the best L2-norm fit for the unwrapped phase map in terms of the first N Zernike polynomials. This step finds a set of coefficients, an , with n 1; 2; 3; …; N such that the L2-norm error, represented as ‖Φw −
N X n1
a n Zn ‖ 2 ;
is minimized. Since Zernike polynomials form an orthogonal set of functions over a circle, finding
an Zn ; Φw ∕Zn ; Zn ;
(8)
where (,) denotes the scalar product, ZZ f ; g dxdyf x; y gx; y:
(9)
The integration above is performed over the circular aperture used for defining the Zernike polynomials. We note that the operation represented by PZ is applied to the phase part of U n1, as in Eq. (5). This updated phase solution along with amplitude of U n1 constitutes the next guess solution W n1 . We observe that while the L2-norm reducing iteration decreases the cost function, CW; W , the action of the operator, PZ , will tend to increase its numerical value. When these two parts of the iteration act in opposing directions, the overall solution shows very small overall change. The iteration may, therefore, be stopped using some threshold, ε, on the change in the solution in successive iterations, ‖W n1 − W n ‖2 < ε: ‖W n ‖2
(10)
While phase unwrapping may be considered an extra step, several standard methods, such as the Goldstein method [21], are available that may be utilized for this purpose. We observe that the two parts of the iterative procedure result in a solution that is smooth to the extent as represented by Zernike expansion and also matches the interference data. While we have not given any convergence proofs that are out of the scope of this paper, we believe that representing the phase part of the wavefront as a linear combination of Zernike polynomials is a strong constraint that is helpful in restricting the search space for the solution and avoiding problems such as stagnation due to local minima. Based on the form of the constraint that restricts the phase function to be a linear combination of finite number, N, of Zernike polynomials, it is clear that the present method is limited to smooth surfaces and will not include cases involving high frequency features such as phase steps. We believe that in order to address the phase estimation problem for such high frequency cases, simply adding number of Zernike terms may not be the best solution. Instead, a modification of the constraint will be required to include Total Variation or L1-norm based penalties. We will not include these cases in the present paper, but report on such generalization of the constraint in the future. 3. Numerical Simulations and Experimental Results
In this section we describe our simulations and experiments in order to test the iterative phase 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6715
estimation method described in Section 2. For simulations, we assumed a test wavefront in the form of a polynomial given by Wx; y expfikβ1 xy β2 x2 y2 β3 x3 y3 g: (11) The reference beam (from null surface) was assumed to be a quadratic wavefront given by Rx; y expikαx2 y2 :
(12)
The numerical values of parameters in Eqs. (11) and (12) are: β1 6.2λ, β2 8λ, β3 2.2λ, α 4λ, and k denotes the wave number, 2π∕λ. The interferogram and the simulated test surface profile are shown in Figs. 1(a) and 1(b), respectively. Figures 1(c) and 1(d) show the recovered phase map and the interferogram generated using the recovered phase map, respectively. The result was obtained in four iterations starting with a zero image as the first guess solution. The number of Zernike terms used for modeling is N 36. The reason for using 36 Zernike polynomials is that this set provides radial polynomials up to the 8th degree and angular terms up to 7θ. In most routine applications of interferometric surface testing this is the typical order that is used to describe surface profiles [22,23]. The relative squared L2-norm error, as in Eq. (2), is plotted in Fig. 2(a) and is seen to reduce to 0.7%. The error between the test phase [Eq. (11)] and the recovered phase is shown in Fig. 2(b). We observe that the rms error between the test and recovered wavefronts is 0.003 waves after four iterations. The relative change in the solution between successive iterations was observed to
Fig. 1. Simulation for single shot interferogram analysis using the proposed iterative method. (a) Interferogram, (b) test phase map (vertical axis unit is waves), (c) recovered phase map from single interferogram frame, (d) interferogram generated using the recovered phase map. 6716
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
Fig. 2. Recovery error for proposed iteration procedure. (a) Relative L2-norm error between the interferogram data and interferogram model [see Eq. (2)], (b) error map between the known test phase and recovered phase map.
be less than 0.1% at four iterations. Our preliminary code written in MATLAB requires less than 30 s to complete 4–5 iterations on a laptop computer with a 1.5 GHz processor. Approximately 80% of this time is taken by the phase unwrapping algorithm. The computational efficiency of the method can be improved if we can reduce the time taken for the phase unwrapping process. This can be done by using a faster algorithm for phase unwrapping. In order to test our algorithm further we used experimental datasets from a commercial 9-step phase shifting Fizeau interferometer shown schematically in Fig. 3. A number of datasets were recorded that correspond to various types of imaging lenses, and other optical elements involving curved and flat surfaces. Spherical reference surfaces were used to test the curved surfaces in the Fizeau interferometer setup. The performance of the optimization approach is comparable in all the cases. As an illustration, we
Fig. 3. Schematic diagram of Fizeau interferometer setup.
show one of the datasets for a curved surface, in this paper, that would typically be tested using a multistep phase shifting approach. Nine interferograms are recorded with relative phase shifts of π∕2 between consecutive frames that are achieved using a calibrated PZT stepper. The nine frames may be used to recover the phase map from the interferogram using the phase shifting formula [24], −4I 2 12I 4 − 12I 6 4I 8 : arctan I 1 − 8I 3 14I 5 − 8I 7 I 9 (13)
Φtest − Φnull
Here, the intensity patterns corresponding to the nine interferograms are denoted by symbols I n with n 1; 2; …; 9. With the prior knowledge of the null wavefront phase one may recover the unknown test object phase. The null wavefront, if not known, can also be assumed to have tilt and quadratic terms that can be estimated approximately from the fringe pattern center and the number of closed fringes. The final estimate of the test wavefront will then be with respect to this choice of null surface. The first of the nine interferogram frames was also used in our iterative procedure in order to determine the test wavefront phase. Figures 4(a)–4(c) show three of the nine phase shifted interferograms. The recovered test object phase maps using the phase shifting method (with 9 frames) and the proposed iterative method (with single frame) are shown in Figs. 5(a) and 5(b), respectively. Figures 5(c) and 5(d) show 3D surface maps corresponding to phase maps shown in Figs. 5(a) and 5(b), respectively. From numerical comparison we observe that the rms error between the phase maps estimated from the 9-step phase shifting method and the single step iterative method is equal to 0.008 waves showing excellent agreement between the two methods.
Fig. 5. Comparison of single shot iterative method with experimental data from a phase shifting Fizeau interferometer. (a) Phase map obtained using 9-step phase shifting method and (b) phase map obtained using a single frame out of the phase shifting data. (c), (d) 3D surface plots corresponding to phase maps in (a), (b) respectively. The rms error between the two recoveries is 0.008 waves.
approach. One of the four interferograms (corresponding to zero reference phase) was then used for phase estimation using the iterative procedure described in Section 2. Figure 6(a) shows one of the four interferograms for a simulated test object with artificial noise. Figures 6(b) and 6(c) show the phase maps estimated using the phase shifting method and the iterative optimization procedure,
4. Simulation Results with Artificially Introduced Noise
In order to understand the noise sensitivity of our proposed iterative method we used the simulations with test surface as in Section 3, Fig. 1. Uniform random noise patterns in the range (−0.5, 0.5) radians were added to the simulated test wavefront to generate four phase shifted interferograms (corresponding to reference phase shift of 0, π∕2, π, 3π∕2). The four interferograms were first used to determine the test object phase using the standard phase shifting
Fig. 4. Three of the nine phase shifted interferogram frames recorded on the Fizeau interferometer setup shown in Fig. 3.
Fig. 6. Simulations with artificially introduced noise in range [−0.5, 0.5] radians. (a) One of the four simulated noisy interferograms; (b), (c) are estimated phase maps obtained from the noisy interferogram with the phase shifting approach, and proposed iterative approach; (d), (e) are errors between the known ground truth and recovered phase maps using the phase shifting and the iterative method. [Colorbar is in waves with values in range (−0.02, 0.018).] 10 October 2014 / Vol. 53, No. 29 / APPLIED OPTICS
6717
respectively. Figures 6(d) and 6(e) show the error between the known ground truth test object and the phase maps obtained by the phase shifting and optimization methods, respectively. We observe that the rms phase errors are similar in two cases and approximately equal to 0.002 waves. The error map for the solution obtained using the optimization method is much smoother, as expected, due to the inherent smoothness property associated with this method. This illustration shows that the single-shot iterative optimization method has similar rms noise performance compared to the multishot phase shifting method. 5. Conclusion
We have demonstrated an optimization approach to phase estimation from a single interferogram frame. Both simulation and experimental results were described and they show excellent phase map recovery for this single-shot approach. The performance of this iterative method is observed to be within λ∕100 rms error with respect to the standard multishot phase shifting approaches, as we have shown with experimental data and simulated noisy data. The method, if implemented with dedicated processing hardware, may provide real time feedback on optical metrology to the polishing machine during the manufacturing process. The proposed method is general enough and can be used in any standard interferometer configuration without much difficulty. We believe that the proposed approach can enable low cost high quality optical metrology systems that are simple to set up. The phase estimation method as presented in this work is limited to smooth surfaces and we will report on possible generalization of the constraint term to include high frequency features, such as steps, in future work. Financial support in order to pursue this allied work is acknowledged from the Department of Biotechnology, Ministry of Science and Technology of India (BT/PR8008/MED/32/306/2013). The authors are also grateful to Dr. A. K. Gupta, Director, IRDE Dehradun, for allowing the use of the experimental facility for data collection. References 1. D. Malacara, M. Servin, and Z. Malacara, Interferogram Analysis for Optical Testing, 2nd ed. (CRC Press, 2005). 2. J. Bruning, D. Herriot, J. Gallagher, D. Rosenfeld, A. White, and D. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt 13, 2693–2703 (1974). 3. K. Creath, “Phase shifting interferometry techniques,” in Progress in Optics, E. Wolf, ed. (Elsevier, 1988), Vol. 26, pp. 357–373.
6718
APPLIED OPTICS / Vol. 53, No. 29 / 10 October 2014
4. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. A 72, 156–160 (1982). 5. T. Kreis, “Digital holographic interference-phase measurement using the Fourier transform method,” J. Opt. Soc. Am. A 3, 847–855 (1986). 6. M. Takeda, “Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: an overview,” Ind. Metrol. 1, 79–99 (1990). 7. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36, 4540–4548 (1997). 8. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Fringe-follower regularized phase tracker for demodulation of closed-fringe interferograms,” J. Opt. Soc. Am. A 18, 689–695 (2001). 9. J. L. Marroquin, L. Rivera, S. Botello, R. Rodriguez-Vera, and M. Servin, “Regularization methods for processing fringe pattern images,” Appl. Opt. 38, 788–794 (1999). 10. R. Legarda-Saenz, W. Osten, and W. Juptner, “Improvement of the regularized phase tracking technique for the processing of nonnormalized fringe patterns,” Appl. Opt. 41, 5519–5526 (2002). 11. 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). 12. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General ndimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20, 925– 934 (2003). 13. M. Servin, J. L. Marroquín, and J. A. Quiroga, “Regularized quadrature and phase tracking from a single closed-fringe interferogram,” J. Opt. Soc. Am. A 21, 411–419 (2004). 14. Q. Kemao and S. H. Soon, “Sequential demodulation of a single fringe pattern guided by local frequencies,” Opt. Lett. 32, 127–129 (2007). 15. K. Khare, P. T. Samsheer Ali, and J. Joseph, “Single shot high resolution digital holography,” Opt. Express 21, 2581–2591 (2013). 16. P. T. Samsheerali, K. Khare, and J. Joseph, “Quantitative phase imaging using single shot digital holography,” Opt. Commun. 319, 85–89 (2014). 17. J. Wyant and K. Creath, “Basic wavefront aberration theory for optical metrology,” in Applied Optics and Optical Engineering, R. R. Shannon and J. Wyant, eds. (Academic, 1992), Vol. XI, pp. 1–53. 18. M. Bertero and P. Boccacchi, Introduction to Inverse Problems in Imaging (IOP, 1998). 19. D. H. Brandwood, “A complex gradient operator and its application in adaptive array theory,” Proc. IRE 130, 11–16 (1983). 20. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University, 2004). 21. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping (Wiley, 1998). 22. C. J. Evans, “PVr—a robust amplitude parameter for optical surface specification,” Opt. Eng. 48, 043605 (2009). 23. C. J. Evans, R. E. Parks, P. J. Sullivan, and J. S. Taylor, “Visualization of surface figure by the use of Zernike polynomials,” Appl. Opt. 34, 7815–7819 (1995). 24. K. Creath, “Choosing a phase-measurement algorithm for measurement of coated LIGO optics,” Proc. SPIE 4101, 46– 55 (2000).
