1 April 2011 Wavefront reconstruction using multiple directional derivatives and Fourier transform
Author Affiliations +
We present a Fourier-based regularized method for reconstructing the wavefront from multiple directional derivatives. This method is robust to noise, and is specially suited for deflectometry measurement.
Legarda-Saenz and Espinosa-Romero: Wavefront reconstruction using multiple directional derivatives and Fourier transform

The deflectometry technique is used for a surface measurement where the local slopes are optically measured and the surface is reconstructed using an integration procedure. This approach has several advantages over a direct height measurement.1, 2 The resultant measurement of this technique is the directional derivative of a wavefront W, which is defined by


[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} {\phi}^{k}({\bf r}) = \triangledown W \cdot {\bf v}_{k} {+} \eta^{k} = \cos \varphi_{k} \frac{\partial W( {\bf r})}{\partial x} + \sin \varphi_{k} \frac{\partial W({\bf r})}{\partial y} + {\eta}^{k},\!\!\!\!\!\!\nonumber\\ \end{eqnarray} \end{document} φk(r)=W·vk+ηk=cosϕkW(r)x+sinϕkW(r)y+ηk,
where r = (x, y) is a position in an M × N rectangular grid of pixels, vk = (cos φ, sin φ) is the measurement direction of a sensor, k = 1, …, K with K ⩾ 2, and ηk represents the measurement and processing errors.

A common way to integrate the derivatives is to use a global integration method,3, 4 where most of the integration techniques only use two orthogonal directions in the measurement, that is ϕ1 and ϕ2 with v1 = (1, 0) and v2 = (0, 1). Recently, in Ref. 5 a procedure was proposed that uses multiple directional derivatives, which reduces the noise contribution in the reconstruction of the wavefront. This consists of acquiring k-number of directional derivatives with different measurement direction, and the estimated wavefront We is computed as the minimizer of the following cost functional:


[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} U(W) &=& \sum _{k}\sum _{(x,y)\in L} \lbrace \left[W(x+1,y) -W(x,y)\right]\cos \varphi _{k}\nonumber \\[-5pt] && + [W(x,y+1) -W(x,y)]\sin \varphi _{k} - \phi ^{k}(x,y) \rbrace ^{2},\nonumber\\ \end{eqnarray} \end{document} U(W)=k(x,y)L{W(x+1,y)W(x,y)cosϕk+[W(x,y+1)W(x,y)]sinϕkφk(x,y)}2,
where L defines the region with valid data.

One drawback of the above cost functional is the computational time employed to found the minimizer of the functional. An alternative to reduce the processing time is to perform the integration through faster methods, that is to perform the integration in the frequency domain. To our knowledge, in the literature there are only two very similar works related to the integration of directional derivatives in frequency domain;6, 7 both methods search an approximation between the measurement data and the directional derivative model in the least-squares sense. One problem with this least-squares approach is its sensitivity to a high level of noise, even when these methods use a large number of derivatives. This situation is illustrated in Fig. 1: (a) shows a synthetic wavefront used as reference, (b) shows the wavefront reconstructed with the method reported in Ref. 5, and (c) shows the wavefront reconstructed with the method reported in Ref. 7. In all cases, a Gaussian noise with a signal-to-noise ratio equal to −20 db was added to 10 directional derivatives employed in the reconstructions. As one can see, the reconstructed wavefront shows artifacts due to noise despite the large number of derivatives utilized on the methods.

Fig. 1

(a) Synthetic wavefront. Wavefronts reconstructed with (b) Ref. 5, (c) Ref. 7, and (d) the proposed method with λ = 20.


To overcome this drawback and compute a fast reconstruction, we add to the cost functional given in Eq. 2, a term which assumes that the wavefront is continuous.8 Thus, the wavefront We is estimated as the minimizer of our proposed regularized cost functional


[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{eqnarray} J(W) &=& \sum _{k}\int\!\!\int |\triangledown W\cdot {\bf v}^{k} - \phi ^{k} |^{2} {\rm d}x\,{\rm d}y \nonumber \\ && + \lambda \int\!\!\int (|W_{xx}|^{2} + |W_{xy}|^{2} + |W_{yy}|^{2} )\, {\rm d}x\,{\rm d}y, \end{eqnarray} \end{document} J(W)=k|W·vkφk|2dxdy+λ(|Wxx|2+|Wxy|2+|Wyy|2)dxdy,
where the subscripts indicate partial derivatives, and λ is the regularization parameter. The second term of this cost functional penalizes the strong oscillations in the curvature of the reconstructed wavefront through the term λ.9

The minimization of the above cost function is performed by means of Fourier transform theory:10 Taking the Fourier transform and using its properties, we obtain the minimizer of Eq. 3 expressed as


[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} {\mathfrak F}\lbrace W_{e}\rbrace = {\displaystyle {i\;\sum \limits _{k} [ {\mathfrak F}\lbrace \phi ^{k}\rbrace ({\bf q}\cdot {\bf v}_{k})] \over \sum \limits _{k}[({\bf q}\cdot {\bf v}_{k})^{2} ] + \lambda \;||{\bf q}||^{2}}}, \end{equation} \end{document} F{We}=ik[F{φk}(q·vk)]k[(q·vk)2]+λ||q||2,
where [TeX:] $i = \sqrt{-1}$ i=1 , q = (u, v) is the position vector on frequency domain, and [TeX:] $ {\mathfrak F}\lbrace \cdot \rbrace$ F{·} denotes the Fourier transform. This expression becomes the same reported in Ref. 10 for the case of K = 2 with v1 = (1, 0) and v2 = (0, 1), and for the case of λ = 0, Eq. 4 becomes the same expression reported in Ref. 7. For practical purposes, one has to consider that the discrete Fourier transform always needs full fields of valid data. The missing data, for instance due to occlusions in the measurement, has to be extrapolated by use of a Gershberg-type algorithm or an alternative technique.11, 12

The performance of the proposed technique is illustrated by two experiments. All numerical experiments were made in a 2.53 GHz Pentium Dual Core PC with 8 GB of main memory using Ubuntu 9.04 as the operative system. The algorithm used to compute the discrete Fourier transform was the FFTW library.13 The first experiment was a numerical estimation using the synthetic wavefront shown in Fig. 1. As mentioned above, a Gaussian noise with a signal-to-noise ratio equal to −20 db was added to 10 directional derivatives, which were generated by use of Eq. 1; the direction vk was equally spaced from 0° to 180°, where every derivative is a 512 × 512 matrix. The resultant wavefront is shown in Fig. 1 using λ = 20. As one can observe, the artifacts in the reconstructed wavefront are notably reduced. The mean square error (MSE) obtained in each reconstruction is shown on Table 1.

Table 1

Mean square error for the reconstructed wavefronts.

MethodMSETime employed (ms)
Ref. 50.2431579011
Ref. 70.320633472
Eq. 4, λ = 200.147356472

The second experiment was the reconstruction of a progressive lens surface. An example of the fringe pattern used in the experiment is shown in Fig. 2. Five directional derivatives were acquired using Eq. 1 with vk equally spaced from 0° to 180° with steps of 30°, where every image has a resolution of 2048 × 1552 pixels. The fringe patterns were processed using the procedure reported in Ref. 14. The computational time employed for the reconstruction was nearly 879 ms using the proposed technique with λ = 5. In the case of the procedure of Ref. 5, the time employed was 272,275 ms. The resultant reconstructions are shown in Figs. 2 and 2. As one can observe, both reconstructions are similar where the uncertainty of both measurements are about 0.04 mm. However, there is a huge difference in the computational process employed for the reconstruction.

Fig. 2

Experimental results: (a) Fringe pattern of a progressive lens in the deflectometry setup. Wavefronts reconstructed with (b) Ref. 5 and (c) the proposed method with λ = 5. The heights are represented in the same gray-scale.


The performance of the cost functional proposed here, Eq. 3 and their minimizer given in Eq. 4, was proved with two fast and accurate reconstructions using a large number of directional derivatives, as it was shown in the results described above. One extra advantage of the proposed method is its feasibility to be implemented on a dedicated hardware for processing in real-time, which is the aim of future research.


The authors were supported in part by grants of the Consejo Nacional de Ciencia y Tecnologia (CONACyT).


1.  C. Wagner and G. Häusler, “Information theoretical optimization for optical range sensors,” App. Opt. 42, 5418–5426 (2003). 10.1364/AO.42.005418 Google Scholar

2.  M. C. Knauer, J. Kaminski, and G. Häusler, “Phase measuring deflectometry: a new approach to measure specular free-form surfaces,” Proc. SPIE 5457, 366–376 (2004). 10.1117/12.545704 Google Scholar

3.  K. Schlüns and R. Klette, “Local and global integration of discrete vector fields,” in Advances in Computer Vision, pp. 149–158, Springer, Wien (1997). Google Scholar

4.  S. Ettl, J. Kaminski, M. C. Knauer, and G. Häusler, “Shape reconstruction from gradient data,” App. Opt. 47, 2091–2097 (2008) and references herein. 10.1364/AO.47.002091 Google Scholar

5.  R. Legarda-Saenz, M. Rivera, R. Rodriguez-Vera, and G. Trujillo-Schiaffino, “Robust wave-front estimation from multiple directional derivatives,” Opt. Lett. 25, 1089–1091 (2000). 10.1364/OL.25.001089 Google Scholar

6.  R. Bronstein, M. Werman, and S. Peleg, “Surface reconstruction from derivatives,” in Proceedings of 11th IAPR International Conference on Pattern Recognition, pp. 391–394, IEEE, Netherlands (1992). Google Scholar

7.  S. Velghe, J. Primot, N. Guérineau, M. Cohen, and B. Wattellier, “Wave-front reconstruction from multidirectional phase derivatives generated by multilateral shearing interferometers,” Opt. Lett. 30, 245–247 (2005). 10.1364/OL.30.000245 Google Scholar

8.  T. Wei and R. Klette, “Height from gradient using surface curvature and area constraints,” in Third Indian Conference on Computer Vision, Graphics and Image Processing (ICVGIP), pp. 204–210, ICVGIP (2002). Google Scholar

9.  M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging, IOP Publishing, Bristol (1998). Google Scholar

10.  T. Wei and R. Klette, “Regularization method for depth from noisy gradient vector fields,” Tech. Report CITR-TR 115 (2002),  http://hdl.handle.net/2292/2850Google Scholar

11.  F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” App. Opt. 30, 1325–1327 (1991). 10.1364/AO.30.001325 Google Scholar

12.  A. Dubra, C. Paterson, and C. Dainty, “Wave-front reconstruction from shear phase maps by use of the discrete Fourier transform,” App. Opt. 43, 1108–1113 (2004). 10.1364/AO.43.001108 Google Scholar

13.  M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93, 216–231 (2005). 10.1109/JPROC.2004.840301 Google Scholar

14.  R. Legarda-Saenz, “Robust wavefront estimation using multiple directional derivatives in moire deflectometry,” Opt. Las. Eng. 45, 915–921 (2007). 10.1016/j.optlaseng.2007.04.004 Google Scholar

Ricardo Legarda-Sáenz, Arturo Espinosa-Romero, "Wavefront reconstruction using multiple directional derivatives and Fourier transform," Optical Engineering 50(4), 040501 (1 April 2011). https://doi.org/10.1117/1.3560540

Back to Top