Shack-Hartmann spot dislocation map determination using an optical flow method

We present a robust, dense, and accurate Shack-Hartmann spot dislocation map determination method based on a regularized optical flow algorithm that does not require obtaining the spot centroids. The method is capable to measure in presence of strong noise, background illumination and spot modulating signals, which are typical limiting factors of traditional centroid detection algorithms. Moreover, the proposed approach is able to face cases where some of the reference beam spots have not a corresponding one in the distorted Hartmann diagram, and it can expand the dynamic range of the Shack-Hartmann sensor unwrapping the obtained dense dislocation maps. We have tested the algorithm with both simulations and experimental data obtaining satisfactory results. A complete MATLAB package that can reproduce all the results can be downloaded from [http://goo.gl/XbZVOr]. ©2014 Optical Society of America OCIS codes: (100.0100) Image processing; (120.5050) Phase measurement. References and links 1. R. G. Lane and M. Tallon, “Wave-front reconstruction using a Shack-Hartmann sensor,” Appl. Opt. 31(32), 6902–6908 (1992). 2. D. N. Neal, J. Copland, and D. Neal, “Shack-Hartmann wavefront sensor precision and accuracy,” Proc. SPIE 4779, 148–160 (2002). 3. M. Nicolle, T. Fusco, G. Rousset, and V. Michau, “Improvement of Shack-Hartmann wave-front sensor measurement for extreme adaptive optics,” Opt. Lett. 29(23), 2743–2745 (2004). 4. K. L. Baker and M. M. Moallem, “Iteratively weighted centroiding for Shack-Hartmann wavefront sensors,” Opt. Express 15(8), 5147–5159 (2007). 5. L. A. Poyneer, D. W. Palmer, K. N. LaFortune, and B. Bauman, “Experimental results for correlation-based wavefront sensing,” Proc. SPIE 5894, 58940N (2005). 6. C. Leroux and C. Dainty, “Estimation of centroid positions with a matched-filter algorithm: relevance for aberrometry of the eye,” Opt. Express 18(2), 1197–1206 (2010). 7. J. Vargas, R. Restrepo, J. C. Estrada, C. O. S. Sorzano, Y. Z. Du, and J. M. Carazo, “Shack-Hartmann centroid detection using the spiral phase transform,” Appl. Opt. 51(30), 7362–7367 (2012). 8. J. Pfund, N. Lindlein, and J. Schwider, “Dynamic range expansion of a Shack-Hartmann sensor by use of a modified unwrapping algorithm,” Opt. Lett. 23(13), 995–997 (1998). 9. B. K. P. Horn and B. G. Schunck, “Determining Optical Flow,” Artif. Intell. 17(1-3), 185–203 (1981). 10. J. Vargas, J. A. Quiroga, C. O. S. Sorzano, J. C. Estrada, and M. Servín, “Multiplicative phase-shifting interferometry using optical flow,” Appl. Opt. 51(24), 5903–5908 (2012). 11. J. Vargas, J. A. Quiroga, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Two-step interferometry by a regularized optical flow algorithm,” Opt. Lett. 36(17), 3485–3487 (2011). 12. J. Y. Bouguet, “Pyramidal Implementation of the LK Feature Tracker,” Tech. Rep., Intel Corporation, Microprocessor Research Labs (1999). 13. B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” Proc. of Imaging Understanding Workshop 121–130 (1981). 14. J. Vargas, L. González-Fernandez, J. Antonio Quiroga, and T. Belenguer, “Shack–Hartmann centroid detection method based on high dynamic range imaging and normalization techniques,” Appl. Opt. 49(13), 2409–2416 (2010). 15. J. C. Estrada, J. Vargas, J. M. Flores-Moreno, and J. A. Quiroga, “Windowed phase unwrapping using a firstorder dynamic system following iso-phase contours,” Appl. Opt. 51(31), 7549–7553 (2012). 16. M. A. Navarro, J. C. Estrada, M. Servin, J. A. Quiroga, and J. Vargas, “Fast two-dimensional simultaneous phase unwrapping and low-pass filtering,” Opt. Express 20(3), 2556–2561 (2012). #197550 $15.00 USD Received 12 Sep 2013; revised 18 Nov 2013; accepted 10 Dec 2013; published 14 Jan 2014 (C) 2014 OSA 27 January 2014 | Vol. 22, No. 2 | DOI:10.1364/OE.22.001319 | OPTICS EXPRESS 1319 17. J. Marzat, Y. Dumortier, and A. Ducrot, “Real-time dense and accurate parallel optical flow using CUDA,” WSCG2009 105–111. (2009). 18. R. T. Frankot and R. Chellappa, “A method for enforcing integrability in shape from shading,” IEEE PAMI 10(4), 439–451 (1988).


Introduction
A Shack-Hartmann (SH) wavefront sensor consists on a two-dimensional microlens array focusing on a CCD camera.The measuring principle of a SH sensor consists on determining the local slopes of an incoming wavefront with respect to a reference one, both sampled by the microlens array.Each microlens focuses the incident rays into its focal plane, where a CCD sensor is placed capturing the spot map called Hartmann spot diagram.Typically, the local slopes of the wavefront error are obtained computing the displacements or dislocations   , xy  between corresponding spot centroids of the Hartmann diagrams generated by the distorted and reference beams.The reference Hartmann spot diagram is acquired using a high quality reference beam.The two-dimensional field of partial derivatives of the wavefront error is given from the corresponding centroids dislocations by [1],   where W  is the wavefront error, which corresponds to the difference between the distorted and the reference wavefronts, f is the distance between the microlens and CCD plane, i, j denotes the i th and j th microlens and   , xy  are the displacements or dislocations between corresponding distorted and reference spot centroids.Note that, from a practical point of view, the wavefront error determination consists in three steps.In a first process, the distorted and reference Hartmann diagrams are process in order to obtain the different spot centroids.Secondly, corresponding spots between the distorted and reference Hartmann diagrams are identified, and then, it is determined the field of dislocations   , xy  . Finally, we determine the wavefront error integrating the discrete field of wavefront partial derivatives shown in Eq. (1), by a zonal or modal reconstruction [2].
The location of the centroid is typically determined calculating the spot center of mass.This method is fast, especially for small numbers of pixels per SH spot, and it is extensively used.However, the results obtained from this technique are not accurate in Hartmann patterns affected by noise, background illumination and spot modulating signals.In the past, several approaches have been proposed to accurately obtain the Hartmann centroids in these adverse situations, such as, the weighted centre of gravity [3], the iteratively weighted centroiding [4] and matched filter approaches [5].These methods can obtain accurate results in presence of noise, but they are not appropriate in cases where the Hartmann pattern is affected by background illumination and/or a varying spot modulation or contrast signal [6,7].Recently, in [7] it is presented a Shack-Hartmann centroid detection approach using the spiral phase transform.This method is robust against strong noise, background illumination and spot modulating signals.Observe that any of the algorithms presented above [3][4][5][6][7] will not obtain accurate wavefront error reconstruction results if some of the detected spot centroids are not accurately located as these methods do not introduce any prior knowledge about the smoothness of the wavefront.
The determination of the spot correspondences is typically performed defining for each reference spot a Hartmann diagram region in which the distorted spot must be placed.These regions are defined by the microlens subaperture area.Observe that this procedure defines a limited wavefront error slope dynamic range that the Shack-Hartmann device can measure.However, in [8] it is presented a method to extend this limited dynamic range unwrapping the spot distribution.This approach remains a limitation on the maximum wavefront error curvature.Additionally, the unwrapping process can be challenging in some cases because the low spatial resolution of the dislocation map, especially when some spot centroids have not being detected because the noise and/or a strong background signal.
In this work, we present a method that can obtain dense dislocations maps   , xy  in a single process, without the need of computing the spot centroids.The proposed algorithm is capable of determining dense dislocations maps in adverse situations, such as in presence of strong noise, background illumination and spot modulating signals, as well as, in cases where there is a large number of reference spots that have not corresponding ones in the distorted Hartmann diagram.The proposed method is based on using an optical flow algorithm between the reference and distorted Hartmann diagrams, which imposes prior information about the smoothness of the wavefront error to be detected.This approach calculates the motion between the two Hartmann patterns at every pixel position, which corresponds to the dislocation maps   , in a robust way.Observe that these dense dislocation maps can be easily unwrapped, following a similar procedure than in [8], if needed.
In Section 2, we present the proposed method.Section 3 includes some simulations and in Section 4, we show the experimental results.Finally, in Section 5 the conclusions are drawn.

Proposed method
Optical flow methodologies consist on algorithms that try to calculate the local motion fields between two image frames taken at times t and t + Δt at every pixel position [9].Note that optical flow approaches have been successfully applied previously in optical metrology and specifically in phase demodulation processes to determine a regularized phase direction map [10,11].In this work, we have used an advanced optical flow approach based on a pyramidal implementation of the well-known Lucas-Kanade (LK) optical flow algorithm with iterative refinement [12].
The LK method [13] assumes that the flow or motion between the two images is essentially constant in a local neighborhood of the pixel under consideration, and it solves the basic optical flow equations for all the pixels in that neighborhood, by the least squares criterion.Note that an arbitrary pixel at time t given by   ,, x y t and with intensity   ,, I x y t will move between two consecutive frames to the position   Assuming the movement is small, we can expand the intensity map at tt  as We suppose that the brightness of an object does not change along time, therefore we define the brightness constancy equation as, and Eq. ( 2) can be rewritten using Eq.(3) as 0 where corresponds to the velocity components.As explained before, the LK method assumes that the displacements between the images are small and approximately constant within a neighborhood of the point    (6) This system has usually more equations than unknowns and then is overdetermined.We can obtain x through the Moore-Penrose pseudoinverse of A given by and can be rewritten as .In a practical point of view, in order to avoid instabilities in the inversion of T AA , it is recommended to sum a small number to the diagonal elements and around 0.1.The standard LK optical flow presented above is only valid if the pixel displacements between the images are small in order to assure the first order Taylor expansion shown in Eq. (2).In cases where this first order approximation is not accurate, it is necessary to iterate multiple times on this scheme.Therefore, after the k th iteration, the brightness constancy equation holds Note that we have assumed without lost of generality that A key parameter of the LK optical flow algorithm is the window or neighborhood size.This window size introduces a trade-off between the accuracy and robustness of the approach.The accuracy of the optical flow method relates to the local sub-pixel accuracy that can be achieved, while robustness refers to the sensitivity of the approach respect to disturbances, such as, changes of lighting between the images, noise and outliers, among others.Note that large window sizes provide more robustness, but less accuracy.Additionally, in order to handle large motions between the images, it is necessary to use a large integration window [12].
In order to solve this problematic, in [12] it is presented a pyramidal implementation of the iterative LK algorithm presented above.This method is based on performing a recursive LK optical flow approach over different resolution representations of the input images, called pyramids representations, and obtained downsampling the input images.We can sort these image pyramid representations from coarse to fine resolutions (note that the images with highest resolution are the input images) and then perform the LK approach, beginning with the images with lowest resolution through the ones with highest resolution.This process will provide first a coarse estimation of the motion (velocity) vectors   00 , uv.These vectors can be used as an initial estimation for the next optical flow estimation using the next following higher resolution pyramidal representation images.Note that this process can be used recursively and the brightness constancy equation can be rewritten as the number of pyramidal representation images used.Observe that this advanced optical flow approach based on a pyramidal implementation of the LK optical flow algorithm with iterative refinement provides as clear advantage that the integrating window size can be remained small, assuring high accuracy, while being able to compute large overall pixels displacements with robustness, thanks to the pyramidal decompositions and iterative refinements.
In this work, we have used the presented pyramidal implementation of the LK optical flow algorithm with iterative refinement to compute a dense, accurate and robust measure of the displacement or dislocation maps   , xy  between the distorted and reference Hartmann patterns.In this case,   ,, corresponds to the reference and distorted Hartmann diagrams.In order to assure the brightness constancy assumption shown in Eq. ( 1), we preprocess the Hartmann diagrams by a fast pattern normalization using the spiral phase transform [7] [14], and then, we compute the optical flow approach over this normalized patterns.This normalization performs background suppression and modulation equalization processes.Additionally, there are practical situations where, on one hand, some of the distorted spots have not being detected because the presence of noise and/or a strong background signal.Moreover, the obtained dislocation maps can be highly affected by noise.On the other, in some cases, the wavefront error slope dynamic range exceeds the one that the Shack-Hartmann device in principle can measure.In this situation, the obtained dense can be filtered and smoothed using a Zernike fitting approach.Secondly, these dense displacements maps can be efficiently unwrapped using the methods presented in [15,16] that correspond to fast unwrapping two-dimensional algorithms based on linear recursive filters.These recursive filters can be seen as composed by two terms: a predictor, that is an estimation based on previously unwrapped values, and a corrector, that takes into account wrapped input data to correct the current estimation.These method are robust to noise because their smoothing capabilities.

Simulations
In order to show the effectiveness of the proposed method, we have performed some simulations.In Fig. 1  In Fig. 2, we show the theoretical dislocation maps x  (a) and y  (d), and the displacements obtained from the optical flow approach in (b) and (e), respectively.In Figs.2(c) and 2(d) we show the dislocation maps after performing a Zernike fitting to the displacements maps computed from the optical flow approach.The Root-Mean-Square error (rms) between the recovered and the ground truth displacement maps is of 0.39 and 0.028 px for the  The optical flow approach is performed with three pyramidal representation images and two iterations in each case.Additionally, the Zernike fitting is performed using twenty five Zernike polynomials.Our implementation of the LK optical flow algorithm with pyramidal implementation and iterative refinement requires 50s using a 2.4 GHz laptop and processing with MATLAB.Note that an implementation of this algorithm for real time applications has been proposed in [17], which can be used in applications that require fast processing, if desired.As can be seen from the presented results, the proposed approach can obtain very dense and accurate dislocation maps results.We have compared the accuracy of the proposed method with a robust Shack-Hartmann centroid detection algorithm capable to measure in presence of strong noise, background illumination and spot modulating signals [7].We have also processed the Hartmann patterns shown in Fig. 1 with this algorithm, which has detected the 60% of the centroids.In Fig. 3, we show the captured centroids marked with white squares.The rms between the recovered and the ground truth displacement maps is of 0.7 px for both x  and y  dislocation maps, and the processing time is of 60s.Observe that the proposed optical flow approach is significantly more accurate that the centroid detection algorithm [7], at the same time that it uses all the information contained in the Hartmann patterns and not only the 60% of the centroids.We have also studied the accuracy of the proposed approach for different noise levels.In all cases, we have used the same background, modulation and theoretical dislocation maps  Note that these dislocation maps are dense but appears wrapped because the limited dynamic range of the SH sensor.In order to unwrapping these maps, we have used a similar method that the one presented in [8] but we have used the unwrapping approach presented in [15,16], that is very fast and robust to noise.In Figs.7(b) and 7(e), we show the obtained results and in Figs.7(c) and 7(f), we present the theoretical dislocation maps.The rms between the recovered and the ground truth displacement maps is of 0.33 and 0.33 px for x  and y  respectively.Observe that if we perform the Zernike fitting, the obtained rms are 0.042 and 0.038 px, respectively.

Experimental results
We have also checked our proposed method with experimental data.We have used an Imagine Optics HASO3 commercial Shack-Hartmann wavefront sensor confronted with a Spatial Light Modulator (SLM) from Hamamatsu model X8267.The SH device is equipped with the HASO R-Flex unit that consists on a collimated light source.This collimated beam implies in the SLM and it is modified by a controlled and known local phase shift.Finally, the modified beam enters in the SH sensor and the introduced aberration is measured.We have applied a quadratic phase-shift in the SLM that is shown in Fig. 8  The effective focal length of the HASO3 microlenses is of 5mm and the sampling rate of the CCD camera is approximately 12.5 μm/px.In Fig. 9, we show the obtained reference (a) and distorted (b) Hartmann spot diagrams and in Fig. 10, we present the wavefront error computed by the HASO3 sensor.The obtained rms and peak to valley (pv) are of 0.11 and 0.51 μm, respectively.Finally, in Fig. 11, we show the obtained wavefront error by the proposed method.In this case, the obtained rms and pv are 0.11 and 0.49 μm respectively.We have used the fast method of Frankot and Chellappa [18] for constructing an integrable surface by the dense gradient information obtained from the proposed LK approach.Observe from Figs. 10 and 11 that the results are similar while the proposed approach presents an increase of spatial resolution with respect to the result computed by the HASO3.

Conclusions
In conclusion, we have presented a straightforward, robust, dense, and accurate Shack-Hartmann spot dislocation map determination method based on a regularized optical flow algorithm that does not require obtaining the spot centroids.This approach is capable to measure in presence of strong noise, background illumination and spot modulating signals; that are the typical limiting factors of the traditional centroid detection algorithms.Additionally, the dense dislocation maps obtained by the proposed approach can be easily unwrapped, if needed.We have tested the algorithm with simulations and experimental data obtaining very satisfactory results.A complete MATLAB package that can reproduce all the results can be downloaded from [http://goo.gl/XbZVOr]

,
consideration.Thus, the optical #197550 -$15.00USD Received 12 Sep 2013; revised 18 Nov 2013; accepted 10 Dec 2013; published 14 Jan 2014 (C) 2014 OSA 27 January 2014 | Vol.22, No. 2 | DOI:10.1364/OE.22.001319| OPTICS EXPRESS 1321 flow equation can be assumed to hold for all pixels   T AA is the structure tensor of the image at point   00 , xy .Observe that Eq. (6) is valid only if T AA is invertible, which is verified only if x I and y I are different than zero in the neighborhood of the point   00 , xy
motion vectors using the iterative LK algorithm introduces above.Note that we can assume without of generality that 1 t  and the final motion vectors are given by , we show simulated reference and distorted Hartmann patterns affected by noise, background and modulation signals.The noise is Gaussian and additive with a Noise-to-Signal-Ratio (NSR) of 11%.Note that as we are using white Gaussian noise, the NSR parameter is related with typical readout camera noise.The Hartmann diagram has size of 640 × 640px and contains around 3600 spots.The background and modulation signals correspond both in arbitrary units (a.u).

Fig. 2 .
Fig. 2. Theoretical dislocation maps x  (a) and y  (d) and displacements obtained by the optical flow approach (b) and (e) without Zernike fitting.In (c) and (d) we show the dislocation maps after performing a Zernike fitting.

Fig. 3 .
Fig. 3. Simulated reference (a) and distorted (b) Hartmann pattern shown in Fig. 1 with the detected centroids marked with white squares.
#197550 -$15.00USD Received 12 Sep 2013; revised 18 Nov 2013; accepted 10 Dec 2013; published 14 Jan 2014 (C) 2014 OSA Observe from Figs. 4 and 6 that our proposed method is more robust to noise, and at the same time, uses all the centroids in the Hartmann patterns.Additionally, we have analyzed the performance of the proposed method for different wavefront error dynamic ranges.In Figs.7(a) and 7(d), we show the obtained dislocation map x  and y  by the proposed optical flow algorithm when the simulated wavefront error dynamic range exceeded the, in principle, limited dynamic range of the SH sensor.The background and modulation signals, as as, the noise level used to simulate the respective reference and distorted spot distributions are the same as the used in Fig. 1.

Fig. 7 .
Fig. 7. Obtained dislocation maps x  (a) and y  (d) by the proposed optical flow algorithm when the simulated wavefront error dynamic range exceeded the limited dynamic range of the SH sensor.In (b) and (e) it is shown the obtained results after unwrapping the dislocation maps, and in (c) and (f) the theoretical dislocation maps.

Fig. 9 .
Fig. 9. Obtained reference (a) and distorted (b) Hartmann spot diagrams obtained by the commercial Shack-Hartmann device.