- Email: [email protected]
Robust wavefront estimation using multiple directional derivatives in moiré deflectometry
ARTICLE IN PRESS
Optics and Lasers in Engineering 45 (2007) 915–921 www.elsevier.com/locate/optlaseng
Robust wavefront estimation using multiple directional derivatives in moire´ deflectometry$ Ricardo Legarda-Saenz Facultad de Matematicas, Universidad Autonoma de Yucatan, Apartado Postal 172, 97110 Merida, Yucatan, Mexico Received 8 January 2007; accepted 15 April 2007 Available online 6 June 2007
Abstract In this paper, the advantages of two recently proposed techniques to implement an automatic processing of moire´ patterns are discussed; this process initiates with the demodulation of a large number of experimental deflectograms and ends when a robust wavefront reconstruction is obtained. Experimental deflectograms were used to show through one example the application of these techniques. r 2007 Elsevier Ltd. All rights reserved. PACS: 42.30.Ms; 42.15.Dp; 42.30.Va Keywords: Moire´ deflectometry; Gradient field reconstruction; Fringe processing
1. Introduction Moire´ deflectometry is a well-known technique to measure the ray deflections when travelling through a phase object, based on the moire´ and Talbot effects [1,2]. The experimental setup of the moire´ deflectometry consists of a collimated light source and a pair of Ronchi gratings separated along the axis of the beam propagation. A moire´ pattern, or deflectogram, is formed by superimposing the shadow of the front grating projected on the rear grating. For the deflection measurement, the phase object under analysis is mounted in the course of the collimated beam, and the moire´ pattern is modified due to the ray deflections produced by optical properties of the phase object. In this way, moire´ deflectometry provides the lateral gradients of the optical path difference [1]. Traditionally, the changes in the moire´ pattern are defined as [1] dx ¼
DCðx; yÞ , y
(1)
$ This work was supported by CONACYT Grant 49739, and SEP Grant PROMEP-UADY-PTC-52. Tel./fax: +52 999 942 3140. E-mail address: [email protected]
0143-8166/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2007.04.004
where D is the distance between gratings in the experimental setup, y is the angle between gratings lines, and Cðx; yÞ are the ray deflections measured by the moire´ deflectometry. This deflection can be related with the normal direction of the emerging wavefront; however, the expression given in Eq. (1) is a paraxial approximation of the wavefront gradients, and only using a prior knowledge like geometry and position of the phase object, the gradient field of the propagated wavefront can be estimated properly [3]. A non-paraxial approximation of the changes in the moire´ pattern can be obtained as follows [3]: assuming that the deflections analysis can be treated by purely geometrical optics, consider a point P1 in the front grating with coordinates ðx1 ; y1 ; z1 Þ. If this point is projected to the plane z2 ¼ z1 þ D, then the coordinates of the projected point P01 are defined as ðx1 þ nx ; y1 þ ny ; z1 þ DÞ, and the direction of this projection is defined as cos a ¼
nx ; jP1 P01 j
cos b ¼
ny D , 0 ; and cos g ¼ jP1 P01 j jP1 P1 j (2)
where cos2 a þ cos2 b þ cos2 g ¼ 1, and the terms nx and ny represent the changes of the position. These terms are
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
916
expressed in the following way: cos a , nx ¼ jP1 P01 j cos a ¼ D cos g ny ¼ jP1 P01 j cos b ¼ D
cos b . cos g
(3)
(4)
As it was mentioned previously, the projection direction of the first grating can be considered as the normal direction of the wavefront emerging of the phase object. Thus, Eqs. (3) and (4) are written in terms of the gradients of the wavefront as nx ¼ D
cos a qf ¼D , cos g qx
(5)
ny ¼ D
cos b qf ¼D , cos g qy
(6)
where f is the propagated wavefront. On the other hand, the position of the Ronchi gratings used in the experimental setup can be expressed using the indicial equations [1]. For a perfect collimated beam, the gratings in the experimental setup are defined as x cosðjÞ þ y sinðjÞ ¼ r V1 ¼ mp,
(7)
x cosðj þ yÞ þ y sinðj þ yÞ ¼ r V2 ¼ np,
(8)
where r ¼ ðx; yÞ, V1 ¼ ðcos j; sin jÞ and V2 ¼ ðcosðj þ yÞ; sinðj þ yÞÞ are the normal vectors to the gratings, and they describe the orientations of each grating in the setup, m ¼ 0; 1; 2; . . . ; n ¼ 0; 1; 2; . . . ; and p is the grating period. The superposition of both gratings described by Eqs. (7) and (8) will produce a deflectogram with a combination of low- and high-frequency terms; the resultant beat frequency of the superposition is the moire´ pattern, which is defined as r ðV1 V2 Þ ¼ ðm nÞp.
(9)
Any modification in the gradient field of the wavefront will produce changes in the moire´ pattern. These changes will be reflected on position changes of the projected grating, that is [3] ðx þ nx Þ cosðjÞ þ ðy þ ny Þ sinðjÞ ¼ r V1 þ D,f V1 ¼ mp, x cosðj þ yÞ þ y sinðj þ yÞ ¼ r V2 ¼ np,
ð10Þ (11)
where ,f ¼ ðqf=qx; qf=qyÞ is the wavefront gradients, and the resultant moire´ pattern is described as r ðV1 V2 Þ þ D,f V1 ¼ ðm nÞp.
(12)
As it can be observed in Eqs. (9) and (12), the changes of the moire´ pattern represent the directional derivative of the wavefront defined by the orientation of the first grating, that is, the direction of V1 . This implies that it is necessary at least two deflectograms with different orientations of the first grating to obtain the gradient field of the propagated wavefront [1,3,4]. Once the gradient field was obtained, the
wavefront can be reconstructed by numerical integration [4–6]. Several approaches can be founded in literature concerning with the automatic processing of moire´ patterns (for example [7,8]). However, most of approaches require intensive computational efforts, making the complete automatic evaluation of moire´ deflectograms difficult to implement. Recently, several new results on fringe processing and numerical integration of derivatives have been published in literature. In some of them, very useful features that help into the demodulation of deflectograms and the robust wavefront reconstruction obtained from moire´ deflectometry experiment are found. In this paper, the advantages of two of those techniques in the automatic processing of moire´ patterns, that is, the demodulation of a large number of experimental deflectograms and the robust wavefront reconstruction are discussed. Section 2 describes the techniques used to demodulate the deflectograms and integrate the derivative fields, including the considerations to implement an automatic processing of a sequence of deflectograms. Section 3 shows the results obtained to process experimental deflectograms, and the conclusions are given in Section 4. 2. Techniques used to process moire´ patterns 2.1. Quadrature transform for deflectograms demodulation The first step of the procedure is to extract the directional derivative from the deflectogram. The fringe pattern obtained from moire´ deflectometry experiment is defined by [2] IðrÞmoire ¼ aðrÞ þ bðrÞ 2p ðr ðV1 V2 Þ þ D,f V1 Þ , cos p
ð13Þ
where aðrÞ is the background illumination, and bðrÞ is the amplitude modulation. Examples of these fringe patterns are shown in Fig. 1. As it can be observed in Eq. (13), the phase information of this fringe pattern is given by two terms: the carrier frequency described by the moire´ pattern, Eq. (9), and the spatial variations of the frequency given by the directional derivative described in Eq. (12). The demodulation process of a single fringe pattern defined like Eq. (13) is wellknown: the Fourier transform of the fringe pattern is computed, then the phase information is isolated by filtering in the frequency plane, and finally the phase is recovered by computing the inverse Fourier transform [2,9]. To have this method working, the spatial variations of the frequency must be smooth compared to the carrier frequency, and the information in the frequency plane has to be well separated. The conditions imposed at the Fourier-based demodulation method do not allow to implement an automatic demodulation of a sequence of deflectograms because it is
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
917
Fig. 1. Experimental moire´ patterns. For the left column j ¼ 0 , for the middle column j ¼ 60 , and for the right column j ¼ 120 ; for all the columns y ¼ 4 . The top row shows the deflectograms without phase object, and the bottom row shows the distorted deflectograms due to a progressive multifocal ophthalmic lens.
not possible to use the same filtering function for each orientation of the gratings due to the different values of the directional derivative. Instead, one has to select the appropriate filter to isolate the desired frequency information for each orientation, and, consequently, the processing of a sequence of deflectograms using Fourier methods will be manual, and becomes a time-consuming and errorprone process. To illustrate this situation, an example of processing deflectograms is shown in Fig. 2: panels (a)–(c) of Fig. 2 show the Fourier spectra of the deflectograms shown in panels (d)–(f) of Fig. 1, respectively. In panels (d)–(f) of the Fig. 2, the pass-band filter used to isolate the phase information of panels (a)–(c) of the same figure are shown, respectively. The resultant wrapped phase, after the filtering process, are shown in panels (g)–(i) of Fig. 2. As it can be appreciated, the filters were designed explicitly for each frequency spectra. An alternative to the problems mentioned previously is to estimate the quadrature term of the fringe patterns [10,11]. First, it is assumed that the background illumination aðrÞ is filtered resulting in the following fringe pattern gðrÞ ¼ bðrÞ cosðcðrÞÞ,
(14)
where cðrÞ ¼
2p ðr ðV1 V2 Þ þ D,f V1 Þ. p
The quadrature estimation consists of transforming the fringe pattern shown in Eq. (14) into its quadrature term defined by ^ ¼ bðrÞ sinðcðrÞÞ, gðrÞ
and obtaining the following fringe complex pattern: ^ ¼ bðrÞ exp½icðrÞ, g igðrÞ (15) pffiffiffiffiffiffiffi where i ¼ 1. The bi-dimensional quadrature transform is defined as [10,11] q 1 Ffgg i Q2 fgg ¼ nc ðrÞF jqj ¼ bðrÞ sinðcðrÞÞ,
ð16Þ
where Ffg indicates the Fourier transform, q ¼ ðu1 ; u2 Þ is the vectorial position on the frequency domain, and nc ðrÞ is the unit vector to the isophase contours defined as nc ðrÞ ¼
2 X rcðrÞ ¼ cð2pÞ ek , jrcðrÞj k¼1 k
where ek is the orthonormal base on the space domain, cð2pÞ are the director cosines modulus 2p. In principle, k there is no way to obtain nc ðrÞ; however, for the particular case of fringe pattern with carrier frequency, the term nc ðrÞ is equal to the normal direction of the fringe pattern [10,11]. Once the quadrature term was computed, the wrapped phase is obtained with Q2 fgg W fcg ¼ arctan , (17) g where W is the wrapping operator [12]. Then the carrier frequency is subtracted from the phase information, and the resultant displacement due to the object
ARTICLE IN PRESS 918
R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
Fig. 2. Example of demodulating deflectograms using Fourier methods. The top row shows frequency spectra of deflectograms, the middle row shows the pass-band filters, and the bottom row shows the resultant wrapped phase. See text for detailed explanation.
under analysis can be unwrapped using well-known techniques [12]. In addition to the well-known computation advantages of the discrete Fourier transform, the quadrature transform shown in Eq. (16) does not need an accurate estimation of the term nc ðrÞ for a correct phase demodulation because the errors in the orientation estimation cause only second-order errors on the demodulated phase [10]. This property results an important feature that is used to implement the proposed automatic processing of a sequence of deflectograms. Herein, every deflectogram acquired in the experiment is demodulated using Eq. (16) where the carrier-frequency orientation is assigned to the term nc ðrÞ of the quadrature transform. Despite the possible errors generated by a large variations in the fringepattern orientation, which are not so relevant as it was mentioned previously, the assignation to nc ðrÞ allows to automatically process one deflectogram, given that the orientation is the only input parameter necessary to apply the quadrature transform; this parameter is obtained from the experimental setup. Moreover, using a large number of directional derivative fields in the reconstruction, the overall error-contribution of the demodulation could be
reduced even more, as it occurs with the reconstruction technique described in the following section. 2.2. Wavefront reconstruction using multiple directional derivatives The classical procedure for the wavefront reconstruction is to acquire two orthogonal derivatives, usually parallel to each axis of reference, and integrate this gradient field [5,6]. However, the drawback of these methods is its sensitivity to noise: any noisy data or outliers will propagate through the estimation such that the resultant reconstruction depends completely on the reliability of the estimated gradient field. Recently, a technique that uses a large number of directional derivative has been shown to reduce the contribution of this noise in the integration process [4]. This technique is based on acquiring N different directional derivatives. In each case, the directional derivative is defined as Fk ¼
qf qf cos jk þ sin jk ¼ ,f Vk1 , qx qy
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
where k ¼ 2; 3; . . . ; N is the index of the acquired directional derivative field, and the vector Vk1 is given by the orientation of the first grating in the experimental setup.
919
Once the directional derivative fields have been obtained from the deflectograms, a quadratic cost function is used to integrate the derivative information. This function is
Fig. 3. Resultant demodulation of the moire´ patterns shown in Fig. 1 using the quadrature transform defined in Eq. (16).
Fig. 4. Directional derivatives obtained by subtraction of the wrapped phases shown in Fig. 3.
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
920
defined as [4] X X ½ðfxþ1;y fx;y Þ cos jk UðfÞ ¼ k
ðxþ1;yÞ; ðx;yþ1Þ; ðx;yÞ2R
þ ðfx;yþ1 fx;y Þ sin jk Fkx;y 2 ,
ð18Þ
where R ¼ fðx; yÞ : Fkx;y is availableg. The reconstructed ^ is computed as the minimizer of the cost wavefront f function. The cost function defined in Eq. (18) has two principal computational advantages [4]: first, the minimization does not depend on the number of directional derivative field used in the estimation, and second, the solution of this function become a linear system, so that it can be solved by any well-known technique [13,14]. 3. Experiments and discussion of results The techniques described in Section 2 were used to implement an automatic processing of a sequence of deflectograms obtained from a moire´ deflectometry experiment. The phase object under analysis was a progressive multifocal ophthalmic lens from Essilor International [15]. The experimental setup consisted of an expanded collimated He–Ne laser source and a pair of Ronchi gratings of 150 lines/in separated by one Talbot-distance, and both were mounted over rotary stages. Two set of deflectograms were acquired: one set was taken with the phase object rotating the gratings from 0 to 174 in steps of 6 , and the second set was taken only rotating the gratings in the same way for the first set, but without the phase object to correct the systematic experimental setup errors. Examples of the deflectograms acquired are shown in Fig. 1. Once the fringe patterns were acquired, every pair of deflectograms (with and without the object phase information) are processed using the previously described techniques in the following way: first the discrete Fourier transform of both deflectograms are computed [16], then the background illumination is removed using a simple frequency filter, and the quadrature transform is applied using the experimental orientation of the gratings. An example of the resultant demodulation is shown in Fig. 3. Once the phase information is obtained, the carrier frequency are subtracted from phase information and the resultant directional derivative field are unwrapped using a flood technique [17]. Examples of the wrapped directional derivative fields obtained from the subtraction of the carrier frequency are shown in Fig. 4, where panel (a) was obtained from the subtraction of panels (a) and (d) of Fig. 3; panels (b) and (c) of Fig. 4 were obtained by subtraction of the corresponding panels shown in Fig. 3. Once all available pairs of deflectograms are demodulated, the wavefront emerging from the phase object under analysis is reconstructed using the cost function defined in Eq. (18), where the linear system was solved by gradient conjugate algorithm [14].
Fig. 5. Measurement of a progressive ophthalmic lens. Panel (a) shows the recovered shape of the surface. Panel (b) shows the surface height in two horizontal positions indicated in panel (a).
The overall time employed to process sequentially 30 pairs of deflectograms to obtain the wavefront was 174 s on a Intel Pentium-M 1.40 GHz computer with one Gigabyte of memory, and the reconstructed wavefront is shown in Fig. 5. Comparing the reconstructed surface with the technical specifications [15], the estimated heights of the surface are close to the given data. 4. Conclusions Herein, the advantages of two recently proposed techniques to implement an automatic processing of a large number of moire´ patterns were discussed. The process initiated with the demodulation of experimental deflectograms, where the quadrature transform defined in Eq. (16) was used. Taking advantage of the second-error property of this quadrature transform, a large number of experimental deflectograms were sequentially demodulated using
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
only the orientation of the carrier parameter as an input parameter from the experimental setup. As a complement of the demodulation process, a cost function defined in Eq. (18) was used to reconstruct the wavefront using the previously obtained directional derivatives fields without any computational cost. The performance of both techniques on the wavefront reconstruction was proved by a fast and accurate wavefront reconstruction using a large number of experimental deflectograms, as it can be observed in the results described in the Section 3. An extra advantage of the processing based on the above described techniques is its feasibility to be implemented on a dedicated hardware for processing in real-time. References [1] Kafri O, Glatt I. The physics of moire metrology. New York: WileyInterscience; 1990. [2] Patorski K, Kujawinska M. Handbook of the moire fringe technique. Amsterdam: Elsevier; 1993. [3] Legarda-Saenz R, Rodriguez-Vera R, Rivera M. Nonparaxial method for computing the gradient field of a wavefront using moire deflectometry. Opt Commun 1999;160:214–8. [4] Legarda-Saenz R, Rivera M, Rodriguez-Vera R, Trujillo-Schiaffino G. Robust wavefront estimation from multiple directional derivatives. Opt Lett 2000;25:1089–91.
921
[5] Fried DL. Least-squares fitting a wave-front distortion estimate to an array of phase-difference measurements. J Opt Soc Am 1977;67:370–5. [6] Hudgin RH. Optimal wave-front estimation. J Opt Soc Am 1977;67:378–82. [7] Canabal H, Quiroga JA, Bernabeu E. Automatic processing in moire deflectometry by local fringe direction calculation. Appl Opt 1998;37:5894–901. [8] Quiroga JA, Crespo D, Bernabeu E. Fourier transform method for automatic processing of moire deflectometry. Opt Eng 1999;38:974–82. [9] Ga˚svik KJ. Optical metrology. 3rd ed. Chichester: Wiley; 2002. [10] Larkin KG, Bone DJ, Oldfield MA. Natural demodulation of twodimensional fringe patterns: I. General background of the spiral phase quadrature transform. J Opt Soc Am A 2001;18:1862–70. [11] Servin M, Quiroga JA, Marroquin JL. General n-dimensional quadrature transform and its application to interferogram demodulation. J Opt Soc Am A 2003;20:925–34. [12] Ghigila DC, Pritt MD. Two-dimensional phase unwrapping. Theory, algorithms, and software. New York: Wiley; 1998. [13] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. In: C: The art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [14] Golub GH, Van Loan CF. Matrix computations. 3rd ed. Baltimore: The Johns Hopkins University Press; 1996. [15] Essilor International. Progressive multifocal ophthalmic lens, US patent 5,272,495. [16] Frigo M, Johnson SG. The design and implementation of FFTW3. Proc IEEE 2005;93:216–31. [17] Strobel B. Processing of interferometric phase maps as complexvalued phasor images. Appl Opt 1996;35:2192–8.
Optics and Lasers in Engineering 45 (2007) 915–921 www.elsevier.com/locate/optlaseng
Robust wavefront estimation using multiple directional derivatives in moire´ deflectometry$ Ricardo Legarda-Saenz Facultad de Matematicas, Universidad Autonoma de Yucatan, Apartado Postal 172, 97110 Merida, Yucatan, Mexico Received 8 January 2007; accepted 15 April 2007 Available online 6 June 2007
Abstract In this paper, the advantages of two recently proposed techniques to implement an automatic processing of moire´ patterns are discussed; this process initiates with the demodulation of a large number of experimental deflectograms and ends when a robust wavefront reconstruction is obtained. Experimental deflectograms were used to show through one example the application of these techniques. r 2007 Elsevier Ltd. All rights reserved. PACS: 42.30.Ms; 42.15.Dp; 42.30.Va Keywords: Moire´ deflectometry; Gradient field reconstruction; Fringe processing
1. Introduction Moire´ deflectometry is a well-known technique to measure the ray deflections when travelling through a phase object, based on the moire´ and Talbot effects [1,2]. The experimental setup of the moire´ deflectometry consists of a collimated light source and a pair of Ronchi gratings separated along the axis of the beam propagation. A moire´ pattern, or deflectogram, is formed by superimposing the shadow of the front grating projected on the rear grating. For the deflection measurement, the phase object under analysis is mounted in the course of the collimated beam, and the moire´ pattern is modified due to the ray deflections produced by optical properties of the phase object. In this way, moire´ deflectometry provides the lateral gradients of the optical path difference [1]. Traditionally, the changes in the moire´ pattern are defined as [1] dx ¼
DCðx; yÞ , y
(1)
$ This work was supported by CONACYT Grant 49739, and SEP Grant PROMEP-UADY-PTC-52. Tel./fax: +52 999 942 3140. E-mail address: [email protected]
0143-8166/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2007.04.004
where D is the distance between gratings in the experimental setup, y is the angle between gratings lines, and Cðx; yÞ are the ray deflections measured by the moire´ deflectometry. This deflection can be related with the normal direction of the emerging wavefront; however, the expression given in Eq. (1) is a paraxial approximation of the wavefront gradients, and only using a prior knowledge like geometry and position of the phase object, the gradient field of the propagated wavefront can be estimated properly [3]. A non-paraxial approximation of the changes in the moire´ pattern can be obtained as follows [3]: assuming that the deflections analysis can be treated by purely geometrical optics, consider a point P1 in the front grating with coordinates ðx1 ; y1 ; z1 Þ. If this point is projected to the plane z2 ¼ z1 þ D, then the coordinates of the projected point P01 are defined as ðx1 þ nx ; y1 þ ny ; z1 þ DÞ, and the direction of this projection is defined as cos a ¼
nx ; jP1 P01 j
cos b ¼
ny D , 0 ; and cos g ¼ jP1 P01 j jP1 P1 j (2)
where cos2 a þ cos2 b þ cos2 g ¼ 1, and the terms nx and ny represent the changes of the position. These terms are
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
916
expressed in the following way: cos a , nx ¼ jP1 P01 j cos a ¼ D cos g ny ¼ jP1 P01 j cos b ¼ D
cos b . cos g
(3)
(4)
As it was mentioned previously, the projection direction of the first grating can be considered as the normal direction of the wavefront emerging of the phase object. Thus, Eqs. (3) and (4) are written in terms of the gradients of the wavefront as nx ¼ D
cos a qf ¼D , cos g qx
(5)
ny ¼ D
cos b qf ¼D , cos g qy
(6)
where f is the propagated wavefront. On the other hand, the position of the Ronchi gratings used in the experimental setup can be expressed using the indicial equations [1]. For a perfect collimated beam, the gratings in the experimental setup are defined as x cosðjÞ þ y sinðjÞ ¼ r V1 ¼ mp,
(7)
x cosðj þ yÞ þ y sinðj þ yÞ ¼ r V2 ¼ np,
(8)
where r ¼ ðx; yÞ, V1 ¼ ðcos j; sin jÞ and V2 ¼ ðcosðj þ yÞ; sinðj þ yÞÞ are the normal vectors to the gratings, and they describe the orientations of each grating in the setup, m ¼ 0; 1; 2; . . . ; n ¼ 0; 1; 2; . . . ; and p is the grating period. The superposition of both gratings described by Eqs. (7) and (8) will produce a deflectogram with a combination of low- and high-frequency terms; the resultant beat frequency of the superposition is the moire´ pattern, which is defined as r ðV1 V2 Þ ¼ ðm nÞp.
(9)
Any modification in the gradient field of the wavefront will produce changes in the moire´ pattern. These changes will be reflected on position changes of the projected grating, that is [3] ðx þ nx Þ cosðjÞ þ ðy þ ny Þ sinðjÞ ¼ r V1 þ D,f V1 ¼ mp, x cosðj þ yÞ þ y sinðj þ yÞ ¼ r V2 ¼ np,
ð10Þ (11)
where ,f ¼ ðqf=qx; qf=qyÞ is the wavefront gradients, and the resultant moire´ pattern is described as r ðV1 V2 Þ þ D,f V1 ¼ ðm nÞp.
(12)
As it can be observed in Eqs. (9) and (12), the changes of the moire´ pattern represent the directional derivative of the wavefront defined by the orientation of the first grating, that is, the direction of V1 . This implies that it is necessary at least two deflectograms with different orientations of the first grating to obtain the gradient field of the propagated wavefront [1,3,4]. Once the gradient field was obtained, the
wavefront can be reconstructed by numerical integration [4–6]. Several approaches can be founded in literature concerning with the automatic processing of moire´ patterns (for example [7,8]). However, most of approaches require intensive computational efforts, making the complete automatic evaluation of moire´ deflectograms difficult to implement. Recently, several new results on fringe processing and numerical integration of derivatives have been published in literature. In some of them, very useful features that help into the demodulation of deflectograms and the robust wavefront reconstruction obtained from moire´ deflectometry experiment are found. In this paper, the advantages of two of those techniques in the automatic processing of moire´ patterns, that is, the demodulation of a large number of experimental deflectograms and the robust wavefront reconstruction are discussed. Section 2 describes the techniques used to demodulate the deflectograms and integrate the derivative fields, including the considerations to implement an automatic processing of a sequence of deflectograms. Section 3 shows the results obtained to process experimental deflectograms, and the conclusions are given in Section 4. 2. Techniques used to process moire´ patterns 2.1. Quadrature transform for deflectograms demodulation The first step of the procedure is to extract the directional derivative from the deflectogram. The fringe pattern obtained from moire´ deflectometry experiment is defined by [2] IðrÞmoire ¼ aðrÞ þ bðrÞ 2p ðr ðV1 V2 Þ þ D,f V1 Þ , cos p
ð13Þ
where aðrÞ is the background illumination, and bðrÞ is the amplitude modulation. Examples of these fringe patterns are shown in Fig. 1. As it can be observed in Eq. (13), the phase information of this fringe pattern is given by two terms: the carrier frequency described by the moire´ pattern, Eq. (9), and the spatial variations of the frequency given by the directional derivative described in Eq. (12). The demodulation process of a single fringe pattern defined like Eq. (13) is wellknown: the Fourier transform of the fringe pattern is computed, then the phase information is isolated by filtering in the frequency plane, and finally the phase is recovered by computing the inverse Fourier transform [2,9]. To have this method working, the spatial variations of the frequency must be smooth compared to the carrier frequency, and the information in the frequency plane has to be well separated. The conditions imposed at the Fourier-based demodulation method do not allow to implement an automatic demodulation of a sequence of deflectograms because it is
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
917
Fig. 1. Experimental moire´ patterns. For the left column j ¼ 0 , for the middle column j ¼ 60 , and for the right column j ¼ 120 ; for all the columns y ¼ 4 . The top row shows the deflectograms without phase object, and the bottom row shows the distorted deflectograms due to a progressive multifocal ophthalmic lens.
not possible to use the same filtering function for each orientation of the gratings due to the different values of the directional derivative. Instead, one has to select the appropriate filter to isolate the desired frequency information for each orientation, and, consequently, the processing of a sequence of deflectograms using Fourier methods will be manual, and becomes a time-consuming and errorprone process. To illustrate this situation, an example of processing deflectograms is shown in Fig. 2: panels (a)–(c) of Fig. 2 show the Fourier spectra of the deflectograms shown in panels (d)–(f) of Fig. 1, respectively. In panels (d)–(f) of the Fig. 2, the pass-band filter used to isolate the phase information of panels (a)–(c) of the same figure are shown, respectively. The resultant wrapped phase, after the filtering process, are shown in panels (g)–(i) of Fig. 2. As it can be appreciated, the filters were designed explicitly for each frequency spectra. An alternative to the problems mentioned previously is to estimate the quadrature term of the fringe patterns [10,11]. First, it is assumed that the background illumination aðrÞ is filtered resulting in the following fringe pattern gðrÞ ¼ bðrÞ cosðcðrÞÞ,
(14)
where cðrÞ ¼
2p ðr ðV1 V2 Þ þ D,f V1 Þ. p
The quadrature estimation consists of transforming the fringe pattern shown in Eq. (14) into its quadrature term defined by ^ ¼ bðrÞ sinðcðrÞÞ, gðrÞ
and obtaining the following fringe complex pattern: ^ ¼ bðrÞ exp½icðrÞ, g igðrÞ (15) pffiffiffiffiffiffiffi where i ¼ 1. The bi-dimensional quadrature transform is defined as [10,11] q 1 Ffgg i Q2 fgg ¼ nc ðrÞF jqj ¼ bðrÞ sinðcðrÞÞ,
ð16Þ
where Ffg indicates the Fourier transform, q ¼ ðu1 ; u2 Þ is the vectorial position on the frequency domain, and nc ðrÞ is the unit vector to the isophase contours defined as nc ðrÞ ¼
2 X rcðrÞ ¼ cð2pÞ ek , jrcðrÞj k¼1 k
where ek is the orthonormal base on the space domain, cð2pÞ are the director cosines modulus 2p. In principle, k there is no way to obtain nc ðrÞ; however, for the particular case of fringe pattern with carrier frequency, the term nc ðrÞ is equal to the normal direction of the fringe pattern [10,11]. Once the quadrature term was computed, the wrapped phase is obtained with Q2 fgg W fcg ¼ arctan , (17) g where W is the wrapping operator [12]. Then the carrier frequency is subtracted from the phase information, and the resultant displacement due to the object
ARTICLE IN PRESS 918
R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
Fig. 2. Example of demodulating deflectograms using Fourier methods. The top row shows frequency spectra of deflectograms, the middle row shows the pass-band filters, and the bottom row shows the resultant wrapped phase. See text for detailed explanation.
under analysis can be unwrapped using well-known techniques [12]. In addition to the well-known computation advantages of the discrete Fourier transform, the quadrature transform shown in Eq. (16) does not need an accurate estimation of the term nc ðrÞ for a correct phase demodulation because the errors in the orientation estimation cause only second-order errors on the demodulated phase [10]. This property results an important feature that is used to implement the proposed automatic processing of a sequence of deflectograms. Herein, every deflectogram acquired in the experiment is demodulated using Eq. (16) where the carrier-frequency orientation is assigned to the term nc ðrÞ of the quadrature transform. Despite the possible errors generated by a large variations in the fringepattern orientation, which are not so relevant as it was mentioned previously, the assignation to nc ðrÞ allows to automatically process one deflectogram, given that the orientation is the only input parameter necessary to apply the quadrature transform; this parameter is obtained from the experimental setup. Moreover, using a large number of directional derivative fields in the reconstruction, the overall error-contribution of the demodulation could be
reduced even more, as it occurs with the reconstruction technique described in the following section. 2.2. Wavefront reconstruction using multiple directional derivatives The classical procedure for the wavefront reconstruction is to acquire two orthogonal derivatives, usually parallel to each axis of reference, and integrate this gradient field [5,6]. However, the drawback of these methods is its sensitivity to noise: any noisy data or outliers will propagate through the estimation such that the resultant reconstruction depends completely on the reliability of the estimated gradient field. Recently, a technique that uses a large number of directional derivative has been shown to reduce the contribution of this noise in the integration process [4]. This technique is based on acquiring N different directional derivatives. In each case, the directional derivative is defined as Fk ¼
qf qf cos jk þ sin jk ¼ ,f Vk1 , qx qy
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
where k ¼ 2; 3; . . . ; N is the index of the acquired directional derivative field, and the vector Vk1 is given by the orientation of the first grating in the experimental setup.
919
Once the directional derivative fields have been obtained from the deflectograms, a quadratic cost function is used to integrate the derivative information. This function is
Fig. 3. Resultant demodulation of the moire´ patterns shown in Fig. 1 using the quadrature transform defined in Eq. (16).
Fig. 4. Directional derivatives obtained by subtraction of the wrapped phases shown in Fig. 3.
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
920
defined as [4] X X ½ðfxþ1;y fx;y Þ cos jk UðfÞ ¼ k
ðxþ1;yÞ; ðx;yþ1Þ; ðx;yÞ2R
þ ðfx;yþ1 fx;y Þ sin jk Fkx;y 2 ,
ð18Þ
where R ¼ fðx; yÞ : Fkx;y is availableg. The reconstructed ^ is computed as the minimizer of the cost wavefront f function. The cost function defined in Eq. (18) has two principal computational advantages [4]: first, the minimization does not depend on the number of directional derivative field used in the estimation, and second, the solution of this function become a linear system, so that it can be solved by any well-known technique [13,14]. 3. Experiments and discussion of results The techniques described in Section 2 were used to implement an automatic processing of a sequence of deflectograms obtained from a moire´ deflectometry experiment. The phase object under analysis was a progressive multifocal ophthalmic lens from Essilor International [15]. The experimental setup consisted of an expanded collimated He–Ne laser source and a pair of Ronchi gratings of 150 lines/in separated by one Talbot-distance, and both were mounted over rotary stages. Two set of deflectograms were acquired: one set was taken with the phase object rotating the gratings from 0 to 174 in steps of 6 , and the second set was taken only rotating the gratings in the same way for the first set, but without the phase object to correct the systematic experimental setup errors. Examples of the deflectograms acquired are shown in Fig. 1. Once the fringe patterns were acquired, every pair of deflectograms (with and without the object phase information) are processed using the previously described techniques in the following way: first the discrete Fourier transform of both deflectograms are computed [16], then the background illumination is removed using a simple frequency filter, and the quadrature transform is applied using the experimental orientation of the gratings. An example of the resultant demodulation is shown in Fig. 3. Once the phase information is obtained, the carrier frequency are subtracted from phase information and the resultant directional derivative field are unwrapped using a flood technique [17]. Examples of the wrapped directional derivative fields obtained from the subtraction of the carrier frequency are shown in Fig. 4, where panel (a) was obtained from the subtraction of panels (a) and (d) of Fig. 3; panels (b) and (c) of Fig. 4 were obtained by subtraction of the corresponding panels shown in Fig. 3. Once all available pairs of deflectograms are demodulated, the wavefront emerging from the phase object under analysis is reconstructed using the cost function defined in Eq. (18), where the linear system was solved by gradient conjugate algorithm [14].
Fig. 5. Measurement of a progressive ophthalmic lens. Panel (a) shows the recovered shape of the surface. Panel (b) shows the surface height in two horizontal positions indicated in panel (a).
The overall time employed to process sequentially 30 pairs of deflectograms to obtain the wavefront was 174 s on a Intel Pentium-M 1.40 GHz computer with one Gigabyte of memory, and the reconstructed wavefront is shown in Fig. 5. Comparing the reconstructed surface with the technical specifications [15], the estimated heights of the surface are close to the given data. 4. Conclusions Herein, the advantages of two recently proposed techniques to implement an automatic processing of a large number of moire´ patterns were discussed. The process initiated with the demodulation of experimental deflectograms, where the quadrature transform defined in Eq. (16) was used. Taking advantage of the second-error property of this quadrature transform, a large number of experimental deflectograms were sequentially demodulated using
ARTICLE IN PRESS R. Legarda-Saenz / Optics and Lasers in Engineering 45 (2007) 915–921
only the orientation of the carrier parameter as an input parameter from the experimental setup. As a complement of the demodulation process, a cost function defined in Eq. (18) was used to reconstruct the wavefront using the previously obtained directional derivatives fields without any computational cost. The performance of both techniques on the wavefront reconstruction was proved by a fast and accurate wavefront reconstruction using a large number of experimental deflectograms, as it can be observed in the results described in the Section 3. An extra advantage of the processing based on the above described techniques is its feasibility to be implemented on a dedicated hardware for processing in real-time. References [1] Kafri O, Glatt I. The physics of moire metrology. New York: WileyInterscience; 1990. [2] Patorski K, Kujawinska M. Handbook of the moire fringe technique. Amsterdam: Elsevier; 1993. [3] Legarda-Saenz R, Rodriguez-Vera R, Rivera M. Nonparaxial method for computing the gradient field of a wavefront using moire deflectometry. Opt Commun 1999;160:214–8. [4] Legarda-Saenz R, Rivera M, Rodriguez-Vera R, Trujillo-Schiaffino G. Robust wavefront estimation from multiple directional derivatives. Opt Lett 2000;25:1089–91.
921
[5] Fried DL. Least-squares fitting a wave-front distortion estimate to an array of phase-difference measurements. J Opt Soc Am 1977;67:370–5. [6] Hudgin RH. Optimal wave-front estimation. J Opt Soc Am 1977;67:378–82. [7] Canabal H, Quiroga JA, Bernabeu E. Automatic processing in moire deflectometry by local fringe direction calculation. Appl Opt 1998;37:5894–901. [8] Quiroga JA, Crespo D, Bernabeu E. Fourier transform method for automatic processing of moire deflectometry. Opt Eng 1999;38:974–82. [9] Ga˚svik KJ. Optical metrology. 3rd ed. Chichester: Wiley; 2002. [10] Larkin KG, Bone DJ, Oldfield MA. Natural demodulation of twodimensional fringe patterns: I. General background of the spiral phase quadrature transform. J Opt Soc Am A 2001;18:1862–70. [11] Servin M, Quiroga JA, Marroquin JL. General n-dimensional quadrature transform and its application to interferogram demodulation. J Opt Soc Am A 2003;20:925–34. [12] Ghigila DC, Pritt MD. Two-dimensional phase unwrapping. Theory, algorithms, and software. New York: Wiley; 1998. [13] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. In: C: The art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [14] Golub GH, Van Loan CF. Matrix computations. 3rd ed. Baltimore: The Johns Hopkins University Press; 1996. [15] Essilor International. Progressive multifocal ophthalmic lens, US patent 5,272,495. [16] Frigo M, Johnson SG. The design and implementation of FFTW3. Proc IEEE 2005;93:216–31. [17] Strobel B. Processing of interferometric phase maps as complexvalued phasor images. Appl Opt 1996;35:2192–8.