Accuracy of the reconstruction of the crystalline lens gradient index with optimization methods from Ray Tracing and Optical Coherence Tomography data

The accuracy of the reconstruction of the Gradient Refractive Index (GRIN) of the crystalline lens from optimization methods was evaluated. Different input data, including direction cosines of deflected rays, ray impacts obtained from ray tracing and optical path differences from Optical Coherence Tomography (OCT) were studied. Three different GRIN models were analyzed. The experimental errors of the different experimental input data were estimated from comparisons of measurements and simulations using artificial lenses of known geometries. The experimental errors in the surfaces shape measurement and the influence of the number of rays were also considered. OCT-based input data produced the most accurate GRIN reconstructions. We found that optimization methods (combining global and local search algorithms) allow GRIN reconstructions with acceptable accuracies for moderate noise level. ©2011 Optical Society of America OCIS codes: (330.0330) Vision, color and visual optics; (080.2710) Inhomogeneous optics media; (170.4500) Optical coherence tomography. References and links 1. L. F. Garner and G. Smith, “Changes in equivalent and gradient refractive index of the crystalline lens with accommodation,” Optom. Vis. Sci. 74(2), 114–119 (1997). 2. G. Smith, D. A. Atchison, and B. K. Pierscionek, “Modeling the power of the aging human eye,” J. Opt. Soc. Am. A 9(12), 2111–2117 (1992). 3. B. A. Moffat, D. A. Atchison, and J. M. Pope, “Age-related changes in refractive index distribution and power of the human lens as measured by magnetic resonance micro-imaging in vitro,” Vision Res. 42(13), 1683–1693 (2002). 4. A. de Castro, D. Siedlecki, D. Borja, S. Uhlhorn, J. M. Parel, F. Manns, and S. Marcos, “Age-dependent variation of the gradient index profile in human crystalline lenses,” J. Mod. Opt. (2011). 5. A. Huggert, The iso-indicial surfaces of the human crystalline lens (Acta Ophthalmologica supplementum, Stockholm, 1948). 6. R. Weale, “The ageing eye,” in The Lens, H. K. Lewis, ed. (London, 1963). 7. B. K. Pierscionek, “Refractive index contours in the human lens,” Exp. Eye Res. 64(6), 887–893 (1997). 8. M. C. Campbell and A. Hughes, “An analytic, gradient index schematic lens and eye for the rat which predicts aberrations for finite pupils,” Vision Res. 21(7), 1129–1148 (1981). 9. D. Y. Chan, J. P. Ennis, B. K. Pierscionek, and G. Smith, “Determination and modeling of the 3-D gradient refractive indices in crystalline lenses,” Appl. Opt. 27(5), 926–931 (1988). 10. G. Beliakov and D. Y. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. I. Planar systems,” Appl. Opt. 36(22), 5303–5309 (1997). 11. E. Acosta, D. Vazquez, L. Garner, and G. Smith, “Tomographic method for measurement of the gradient refractive index of the crystalline lens. I. The spherical fish lens,” J. Opt. Soc. Am. A 22(3), 424–433 (2005). 12. D. Vazquez, E. Acosta, G. Smith, and L. Garner, “Tomographic method for measurement of the gradient refractive index of the crystalline lens. II. The rotationally symmetrical lens,” J. Opt. Soc. Am. A 23(10), 2551– 2565 (2006). #149271 $15.00 USD Received 16 Jun 2011; revised 6 Aug 2011; accepted 30 Aug 2011; published 19 Sep 2011 (C) 2011 OSA 26 September 2011 / Vol. 19, No. 20 / OPTICS EXPRESS 19265 13. C. E. Jones, D. A. Atchison, R. Meder, and J. M. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. 45(18), 2352–2366 (2005). 14. S. Kasthurirangan, E. L. Markwell, D. A. Atchison, and J. M. Pope, “In vivo study of changes in refractive index distribution in the human crystalline lens with age and accommodation,” Invest. Ophthalmol. Vis. Sci. 49(6), 2531–2540 (2008). 15. O. Pomerantzeff, M. Pankratov, G. J. Wang, and P. Dufault, “Wide-angle optical model of the eye,” Am. J. Optom. Physiol. Opt. 61(3), 166–176 (1984). 16. L. F. Garner, G. Smith, S. Yao, and R. C. Augusteyn, “Gradient refractive index of the crystalline lens of the Black Oreo Dory (Allocyttus Niger): comparison of magnetic resonance imaging (MRI) and laser ray-trace methods,” Vision Res. 41(8), 973–979 (2001). 17. S. Barbero, A. Glasser, C. Clark, and S. Marcos, “Accuracy and possibilities for evaluating the lens gradientindex using a ray tracing tomography global optimization strategy,” Invest. Ophthalmol. Vis. Sci. 45, 1723 (2004). 18. S. Ortiz, S. Barbero, and S. Marcos, “Computer simulations of optical coherence tomography A-scans: what can we learn about refractive index distribution?” Invest. Ophthalmol. Vis. Sci. 45, 2781 (2004). 19. Y. Verma, K. D. Rao, M. K. Suresh, H. S. Patel, and P. K. Gupta, “Measurement of gradient refractive index profile of crystalline lens of fisheye in vivo using optical coherence tomography,” Appl. Phys. B 87(4), 607–610 (2007). 20. A. de Castro, S. Ortiz, E. Gambra, D. Siedlecki, and S. Marcos, “Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging,” Opt. Express 18(21), 21905–21917 (2010). 21. O. Pomerantzeff, H. Fish, J. Govignon, and C. L. Schepens, “Wide angle optical model of the human eye,” Ann. Ophthalmol. 3(8), 815–819 (1971). 22. I. H. Al-Ahdali and M. A. El-Messiery, “Examination of the effect of the fibrous structure of a lens on the optical characteristics of the human eye: a computer-simulated model,” Appl. Opt. 34(25), 5738–5745 (1995). 23. A. Popiolek Masajada, “Numerical study of the influence of the shell structure of the crystalline lens on the refractive properties of the human eye,” Ophthalmic Physiol. Opt. 19(1), 41–49 (1999). 24. D. O. Mutti, K. Zadnik, and A. J. Adams, “The equivalent refractive index of the crystalline lens in childhood,” Vision Res. 35(11), 1565–1573 (1995). 25. Y.-J. Liu, Z.-Q. Wang, L.-P. Song, and G.-G. Mu, “An anatomically accurate eye model with a shell-structure lens,” Optik Internat. J. Light Electron. Opt. 116(6), 241–246 (2005). 26. C. E. Campbell, “Nested shell optical model of the lens of the human eye,” J. Opt. Soc. Am. A 27(11), 2432– 2441 (2010). 27. A. Gullstrand, Handbuch der Physiologitschen Optik, 3 ed. 1909. J. P. Southall trans., (Optical Society of America, 1924). 28. J. W. Blaker, “Toward an adaptive model of the human eye,” J. Opt. Soc. Am. 70(2), 220–223 (1980). 29. G. Smith, B. K. Pierscionek, and D. A. Atchison, “The optical modelling of the human lens,” Ophthalmic Physiol. Opt. 11(4), 359–369 (1991). 30. B. K. Pierscionek and D. Y. Chan, “Refractive index gradient of human lenses,” Optom. Vis. Sci. 66(12), 822– 829 (1989). 31. R. Navarro, F. Palos, and L. M. González, “Adaptive model of the gradient index of the human lens. II. Optics of the accommodating aging lens,” J. Opt. Soc. Am. A 24(9), 2911–2920 (2007). 32. F. Manns, A. Ho, D. Borja, and J. M. Parel, “Comparison of Uniform and Gradient Paraxial Models of the Crystalline Lens,” Invest. Ophthalmol. Vis. Sci. 51, 789 (2010). 33. H. L. Liou and N. A. Brennan, “Anatomically accurate, finite model eye for optical modeling,” J. Opt. Soc. Am. A 14(8), 1684–1695 (1997). 34. M. A. Rama, M. V. Pérez, C. Bao, M. T. Flores-Arias, and C. Gómez-Reino, “Gradient-index crystalline lens model: A new method for determining the paraxial properties by the axial and field rays,” Opt. Commun. 249(46), 595–609 (2005). 35. J. A. Díaz, C. Pizarro, and J. Arasa, “Single dispersive gradient-index profile for the aging human lens,” J. Opt. Soc. Am. A 25(1), 250–261 (2008). 36. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24(8), 2157–2174 (2007). 37. A. V. Goncharov, M. Nowakowski, M. T. Sheehan, and C. Dainty, “Reconstruction of the optical system of the human eye with reverse ray-tracing,” Opt. Express 16(3), 1692–1703 (2008). 38. I. Grulkowski, M. Gora, M. Szkulmowski, I. Gorczynska, D. Szlag, S. Marcos, A. Kowalczyk, and M. Wojtkowski, “Anterior segment imaging with Spectral OCT system using a high-speed CMOS camera,” Opt. Express 17(6), 4842–4858 (2009). 39. S. Ortiz, D. Siedlecki, L. Remon, and S. Marcos, “Optical coherence tomography for quantitative surface topography,” Appl. Opt. 48(35), 6708–6715 (2009). 40. S. Ortiz, D. Siedlecki, I. Grulkowski, L. Remon, D. Pascual, M. Wojtkowski, and S. Marcos, “Optical distortion correction in optical coherence tomography for quantitative ocular anterior segment by three-dimensional imaging,” Opt. Express 18(3), 2782–2796 (2010). 41. O. N. Stavroudis, “Ray tracing,” in The optics of rays, wavefronts and caustics, N. Y. A. Press, ed. (1972). #149271 $15.00 USD Received 16 Jun 2011; revised 6 Aug 2011; accepted 30 Aug 2011; published 19 Sep 2011 (C) 2011 OSA 26 September 2011 / Vol. 19, No. 20 / OPTICS EXPRESS 19266 42. A. Sharma, D. V. Kumar, and A. K. Ghatak, “Tracing rays through graded-index media: a new method,” Appl. Opt. 21(6), 984–987 (1982). 43. B. D. Stone and G. W. Forbes, “Optimal interpolants for Runge-Kutta ray tracing in inhomogeneous media,” J. Opt. Soc. Am. A 7(2), 248–254 (1990). 44. J. H. Holland, Adaptation in natural and artificial systems (The University of Michigan press, 1975). 45. J. A. Nelder and R. Mead, “A Simplex Method for Function Minimization,” Comput. J. 7, 308–313 (1965). 46. J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the Nelder-Mead simplex method in low dimensions,” SIAM J. Optim. 9(1), 112–147 (1998). 47. D. Vasiljevic, Classical and evolutionary algorithms in the optimization of optical systems (Kluwer Academc Publishers Boston/Dordrecht/London, 2002), p. 296. 48. S. R. Uhlhorn, D. Borja, F. Manns, and J.-M. Parel, “Refractive index measurement of the isolated crystalline lens using optical coherence tomography,” Vision Res. 48(27), 2732–2738 (2008). 49. D. A. Atchison and G. Smith, “Chromatic dispersions of the ocular media of human eyes,” J. Opt. Soc. Am. A 22(1), 29–37 (2005). 50. E. Moreno-Barriuso and R. Navarro, “Laser Ray Tracing versus Hartmann-Sha


Introduction
It is well known that the crystalline lens of the eye in many species shows a gradient index (GRIN) distribution, providing the lens with advantageous optical properties, such as increasing its refractive power with relatively low indices, and reducing the ocular spherical aberration and the light losses generated by reflections on the optical surfaces.In the human eye, it is known that the GRIN distribution changes with accommodation [1] and age [2][3][4].Several methods have been used to estimate the GRIN, in most cases based on in-vitro measurements.
Early estimates of the GRIN structure were extracted from measurements on frozen slices of crystalline lenses with an Abbe refractometer [5,6].Also, a reflectometric fiber optic sensor was used to measure the refractive index along the equatorial and sagittal planes [7].These and other destructive methods, performed in human and non-human lenses, triggered discussions regarding the structure of the GRIN distribution and the actual values of the refractive index at the surface and at the nucleus.Destructive methods intrinsically may alter the GRIN structure under test.As an alternative, several non-destructive methods have been proposed: integral inversion methods, Magnetic Resonance Imaging (MRI)-based methods and optimization methods.
Lateral view of rays traced through a lens allows measuring the deflection of each ray refracted by the lens.Using these experimental data, and under certain assumptions, it is possible to relate mathematically the deflection of the rays with the GRIN distribution.By means of an integral inversion, it is then possible to obtain an estimation of the GRIN.Campbell [8] reconstructed the GRIN of a rat lens from lateral ray tracing data, although the method required matching the external media and surface refractive indices, as well as using a GRIN model with iso-indicial concentric surfaces.Chan [9] and Beliakov [10] expanded the method to any possible rotationally symmetric GRIN, but the method required knowledge of the path of the rays inside the lens.With the aim of overcoming these drawbacks, Acosta and Vázquez [11,12] presented a new algorithm based on ray deflection after (but not inside) the lens, which allowed the reconstruction of mono or bi-polynomial GRIN profiles of crystalline lenses.The algorithm involved experimental ray tracing through the lens at several angles (tomographic method) and was applied to porcine lenses.
The magnetic resonance imaging (MRI) technique, assumes a correlation between the transverse spin-spin relaxation time data and the refractive index [13], and requires careful calibrations.Although the technique has produced the first available data of age-related changes of the human crystalline lens profile [13], and has even been applied in-vivo [14], it has also been contested due to the relatively low-resolution and assumptions involved.
Optimization methods are based on finding the optimal parameters of a specific lens model to fit certain experimental data.The first optimization methods, proposed by Pomerantzeff [15], used as experimental data the focal length and the spherical aberration of the eye, and a shell model lens composed of a large number of layers.The high number of unknown variables (thickness, shapes, and refractive indices of the different layers of the crystalline) made the optimization problem ill-defined, as multiple solutions were possible [8].Garner et al. [16] used optimization methods based on lateral ray tracing data to estimate the GRIN structure of a spherical fish lens.Barbero [17] studied the possibilities of using global search optimization for finding the optimal parameters in non-spherical lenses.
The potential of using optical path differences (OPD) to extract information of the GRIN was studied by Ortiz et al. [18], and the first data on a simple spherical fish lens using OPD from OCT were provided by Verma et al. [19].Local search algorithms were sufficient in this case, due to the simplicity and symmetric properties of the GRIN models used.Recently, we have presented a method for the three-dimensional reconstruction of the GRIN in a porcine lens in vitro from OCT images and a global search algorithm [20].We have also applied the method to 2-D OCT images of human lenses of different ages [4].
In addition to the experimental input data and the reconstruction method itself, a critical aspect in the optimization method is the crystalline lens model assumed (function and number of variables).Several works described the crystalline lens using a shell model, where the different layers have different but homogeneous refractive indices [21][22][23][24][25][26].However, such discontinuous structure of the refractive index distribution would produce a discontinuity in the wave aberration [23].The first description of the crystalline lens (accommodated and unaccommodated) with a continuous function was proposed by Gullstrand [27], and later extended by Blaker [28] and by Smith et al. [29], who noted the need for higher order polynomials to describe experimental findings [30].Exponential functions [2,13,31,32], parabolic [33] or other polynomial expressions [34,35] have been proposed to describe the decrease of refractive index from the center to the surface of the lens.Recently, Goncharov and Dainty [36] presented an eye model with a GRIN, which successfully fitted experimental aberration data both on and off-axis.This model was personalized for an individual eye in a later work [37].
In this study, we compared the accuracy and robustness of optimization methods to reconstruct the GRIN using different input experimental data.Three realizations of the Goncharov model, with increasing number of variables were tested.The evaluation was performed based on computer simulations, and the noise expected in each input data estimated from experimental data collected on artificial lenses with a ray tracing setup and OCT imaging system respectively.The effect of measurement errors and experimental limitations of each technique are addressed.

Methods
Experimental measurements on artificial lenses provided real error data of the input information to the reconstruction algorithms.The experimental data were collected using two different methods: a lateral/transverse ray tracing system (developed specifically for this study) and an OCT imaging system, described in an earlier publication [38].
The experimental data acquisition was simulated computationally on human crystalline lenses, assuming the Goncharov lens models.Ray tracing data (direction cosine of the deflected ray, intercept of the outgoing ray with the posterior lens surface and impact on a plane after the lens) and OCT imaging data (OPD) of the lenses were simulated.The merit functions and optimization tools were applied on the simulated data (using the errors obtained experimentally) to reconstruct the GRIN.To quantify the accuracy of the reconstruction, the root mean square (RMS) error of the difference between the reconstructed GRIN and the nominal GRIN was used (GRIN RMS difference).For the computation of GRIN RMS difference the refractive index was evaluated for both the reconstructed and nominal GRIN on a grid of points (10-µm steps) over a 6-mm pupil.The metric is dimensionless and increases with increasing discrepancy of the reconstructed GRIN with respect to the nominal GRIN

Estimation of experimental errors
a. Experimental ray tracing system A Laser Ray Tracing system that allows collection of the rays' deflections (on a lateral viewing camera) as well as rays' impacts behind the lens (on a second transverse viewing camera) was developed.Figure 1 shows a schematic view of the system and examples of collected images.A He-Ne laser (633 nm) was used for illumination, and an x-y galvanometer scanner (Cambridge Technology) deflected the laser to scan the pupil at different locations.The ray beams were delivered sequentially, sampling the lens (in 100-µm steps) two dimensionally across an 8-mm pupil.Measurement time was around 2 minutes per full scan.In the first mode (lateral viewing, Fig. 1a), the lens was placed in a chamber and filled in with distilled water.A CCD camera (Toshiba Teli, 640x480 pixels, 7.4 µm pixel size) focused at the meridional plane of the lens, recorded images of the refracted rays (Fig. 1b).In the second mode (transverse viewing, Fig. 1c) the lens was placed in air and the impacts of the outgoing rays were captured on a plane after the lens, perpendicular to the optical axis of the set-up.The centroids of the spots were obtained from the spot images captured directly onto a bare CCD (Qimaging Retiga 1300, 6.7 µm pixel size, Fig. 1d).
Data were collected on a homogeneous index lens (KPX088, Newport, f = 71 mm, D = 25.4 mm) and a gradient index lens (GPX-30-60, LightPath, f = 75 mm, D = 30 mm) to estimate the error in the measurement of direction cosines and impacts.The GRIN lens had a continuous axial decreasing variation of refractive index from 1.74 (anterior surface) to 1.67 (posterior surface), and a central thickness of 6 mm.The nominal GRIN profile was described by a 7 th order polynomial function.Data from 80 equally spaced rays sampling an 8-mm pupil were analyzed.
Image processing algorithms were developed (MatLab, MathWorks, Inc.) to estimate the direction cosines (lateral viewing mode), and the ray impacts (transverse viewing mode) from the experimental measurements (Figs.1b and 1c).To calculate the deflection of the rays, the image was thresholded and segmented, and the detected rays were fit by a first order polynomial.The position of impacts was estimated from the centroid of the spots.
The differences between the experimental results and simulated direction cosines (in the first mode) or impacts in a plane after the lens (in the second mode) were calculated.The error in the measurement was estimated in terms of standard deviation of these differences in five repeated measurements.

b. Experimental optical coherence tomography system
Optical Coherence Tomography images were obtained using a custom-developed high resolution spectral domain OCT which uses a SLD diode of 840 nm with a FWHM of 50nm.The system was developed in collaboration with Copernicus University and has been described in detail elsewhere [20,[38][39][40].The system allowed three-dimensional imaging of a transversal area of 15x15mm of the samples, and image quality range of 4 mm in depth, with a scanning pattern of 300x300 A-Scans.The axial resolution of the system was 3.42 um in air.
Two homogeneous lenses (EO45-448 and EO45-705, Edmund Optics, focal length = 60 and 72 mm, respectively, lens diameter = 12.5 mm) were imaged with the Optical Coherence Tomography System three dimensionally to estimate the error in the determination of OPD data through the lens.The OPD data were extracted in an 8-mm pupil range.
Image processing algorithms were developed to segment the anterior and posterior surfaces of the lens images.Fan distortion [39] arising from the architecture of the scanning system was corrected.The optical distance between the estimated ray impacts on the first and second surface of the lens represents the OPD of each ray.
The deviation in the measurement of the OPD of each ray was calculated.The experimental error was estimated as the standard deviation of these deviations in five repeated measurements.In addition, the same measurements were used to extract the error in the lens shape measurements with the OCT system.The error was estimated as the mean percent error of the radii of curvature from spherical fits to the segmented anterior surfaces in 5 repeated measurements on each lens.

Ray tracing and GRIN-reconstruction algorithms a. Ray tracing procedure
We implemented computationally a ray tracing procedure in Matlab (MathWorks, Natic MA) for ray propagation through homogeneous and GRIN media.Two major components of the algorithm include the computation of the intersections of the rays with the lens surfaces, and the calculation of the ray path inside the GRIN medium.
The lens surfaces were described by conics.The intersection points of rays and conic surfaces were calculated using the equations derived by Stavroudis [41].Our ray tracing algorithm through GRIN is based on Sharma algorithm [42].This algorithm consists essentially of applying the Runge-Kutta method to solve numerically the ray equation.The algorithm estimates a series of points through which the ray passes, together with the ray direction at those points.These points define a set of intervals, the width of which is set by the step size.The step size can be adjusted for the desired trade-off between accuracy and computation time.Following Stone et al. [43], we introduced an improvement to the original Sharma algorithm for a better estimation of the intersection of the outgoing ray with the lens surface.
The OPD of each ray was calculated by evaluating the optical path at each Sharma interval.The simplest case would involve assuming a linear trajectory inside each interval.However this approach implies reducing the step size, hence increasing the computation time.We observed that using Hermite polynomial interpolation to estimate the trajectories inside each intervals permits to increase the accuracy while reducing the computation time (see section 3.1).

b. Reconstruction algorithm
The GRIN reconstruction algorithm is based on an optimization method where the variables to optimize are the parameters that define the lens GRIN distribution.A merit function (MF) was constructed to quantify the agreement between the experimental data (ray tracing data or OPD from OCT) and simulated ray tracing for a given set of parameters.The MF was defined as the RMS of this difference.The optimization procedure minimizes this merit function.
The optimization procedure comprises two sequential optimization algorithms.First, a global search algorithm was used to find the global minimum avoiding local minimum trapping.Second, starting from the best solution found by the global search algorithm, a local search algorithm was used to descend to the minimum, because of its faster convergence compared to global search.
For the global optimization we used a genetic algorithm implemented in MatLab which we had applied in previous studies [4,20].Genetic algorithms are based on evolutionary strategies [44].Evolutionary rules are applied to generations of solutions, elite percentage was set to 5%, and mutation ratio was 80%.
For the local optimization, we implemented the Nelder-Mead [45] simplex algorithm.We applied the standard values [46] of reflection (ro = 1), expansion (X = 2), contraction (gamma = 1/2) and shrinkage (sigma = 1/2) to the simplex.The algorithm requires a starting simplex.While MatLab proposes to increase the variables a 5% we have used a smaller simplex to avoid missing the global minimum.

c. Gradient index models
The crystalline lens was described using three different GRIN models, proposed by Goncharov and Dainty [36,37], of increasing level of complexity (2, 3 and 4 variables).Specifically, we used the models for a 20 years old crystalline lens, including both the GRIN parameters and lens surfaces shapes as provided by the authors [36].
The first and second GRIN models are defined by two 4 th order polynomials in radial and axial direction describing the GRIN distribution in anterior and posterior regions of the lens respectively.In the first model (unbalanced model, G20U), the refractive index is constant over the lens surfaces and the GRIN distribution is determined only by two variables: surface and nucleus refractive index.In the second model (balanced model, G20B), the last posterior iso-indicial surface may not be coincident with the posterior lens surface and GRIN distribution is described by three variables: surface index, nucleus refractive index and radius of curvature of the posterior iso-indicial surface.The third model (symmetrical model, G20S) assumes a GRIN distribution with one polynomial described by four variables: surface and nucleus refractive indices, center position and a parameter setting the decay of the refractive index along the radial coordinate ─ this parameter could be related to the refractive index 3mm off the optical axis, at the meridional plane.In this model, the marginal iso-indicial surfaces are more curved than the external surfaces.
As we observed that the merit function of the 2-variable model can be optimally minimized simply using the local search algorithm, we did not use the global search algorithm in this model.For the 3-variable and 4-variable models, the genetic algorithm was composed of 5 generations containing 200 individuals each, and 20 generations of 600 individuals respectively.Constraints -needed in presence of large experimental errors-were applied to the range of refractive indices in the nucleus and surface, as well as the nucleus lens position, according to biologically plausible descriptions of the GRIN distribution.These constraints were included in the merit function using penalty terms [47].For all models, the refractive index was constrained to have a value between 1.355 and 1.44, which fits most of refractive index values in human crystalline lenses found in the literature.For the 4-variable model, the position of the lens nucleus was constrained to range from 1.1 to 1.8 mm.This constraint is justified by the observation that values outside of this range can generate two refractive GRIN maxima.Finally, the refractive index at 3 mm off the optical axis (radial coordinate) was constrained between 1.35 and the nucleus refractive index value.This constraint produces a decrease of the GRIN profile along the radial direction, which is in accordance to the morphology of the crystalline lens.

Simulations
In all the simulations we assumed that the experimental data are limited to rays parallel to the optical axis, and within a 6-mm pupil diameter.Unless otherwise noted, 120 rays were traced, i. e. one ray each 50 µm across the 6 mm pupil.
The reconstruction algorithm was tested assuming simulated input data from five different configurations (three from Laser Ray Tracing, and the other two from OCT): (1) Ray deflections of rays outgoing the lens, as used in previous studies to reconstruct GRIN distributions [8,16].(2) Ray deflections and the intercepts of the outgoing rays with the posterior lens surface, as was used in the tomography method proposed by Acosta and Vázquez [11,12].(3) Ray impacts on a CCD placed after the lens, as shown in Fig. 1c.(4) OPD of each ray intercepting the posterior lens surface, as used in some OCT-based reconstruction methods [4,19], and (5) OPD of each ray that intercepts the posterior lens surface as well as the cuvette surface holding the lens (which can be imaged with the OCT in in-vitro measurements [48]), as used in a previous publication [20].
If a broad band source is used, as it is the case in configurations 4 & 5, the experimental data are associated to the group refractive index.However, as the ray tracing algorithm implicitly uses a phase refractive index, we assume that the refractive index reconstructed with the algorithm is equivalent to the group refractive index.For clarity purposes, in the simulations of this study we used the same refractive index for the five configurations.
To explore the limits of the search algorithm in each configuration, we have reconstructed the GRIN distribution assuming no error in the experimental data, and low noise level (range of Gaussian errors ranging from 10 −6 to 10 −3 ) added to the data.The aim of this simulation was to evaluate the reconstruction accuracy of the algorithm itself.
The reconstruction algorithm was then studied when expected experimental noise level was applied to the input data (direction cosines, impacts or OPDs).Noise was simulated introducing Gaussian error with a standard deviation calculated with the experimental ray tracing and OCT imaging described in section 2.1.
As the lens geometry represents additional input information for the GRIN reconstruction algorithms, we also studied the influence of the experimental errors of the lens shape measurements.The simulation of the ray tracing data was performed on lens surfaces where Gaussian noise (standard deviation given by the error measurement estimated as described in section 2.1) was introduced to the surfaces radii of curvature.As before, the experimental error was also introduced to the input data (direction cosines, intercepts of the outgoing rays with posterior surface, impacts or OPDs).The GRIN reconstruction was evaluated using the nominal lens surface radii.
Finally, we studied the influence of the number of rays (ranging from 6 to 1200) in the reconstruction.

Ray tracing algorithm: validation and influence of the step-size
The performance of the custom ray tracing routines was evaluated against a well established ray tracing software (ZEMAX, Bellevue WA) on a simulated lens with homogeneous refractive index and a simulated GRIN lens (polynomial distribution, Goncharov model 20S).
For these examples, a bundle of rays was traced parallel to the optical axis through a 6-mm pupil.In the ray tracing, several data were computed: the outgoing ray intersection points with the posterior lens surface and a plane surface located 10 mm behind the lens and the OPDs of each ray at the posterior surface of the lens.
The differences between both procedures were within computer machine precision for the homogeneous lens.In the case of the GRIN lens the differences (step-size of 10 µm) were below 10 −4 µm for ray intersections and OPDs at the posterior surface of the lens.These differences are probably due to differences in the ray tracing algorithm (although a more detailed comparison is not possible, as the ZEMAX algorithms are proprietary).
We also studied the trade-off between the precision in the computation of the ray trajectories inside the GRIN and the computation time imposed by the Sharma step-size.Figure 2 shows the differences in the estimation of the ray tracing data (for all the experimental data under study) for different step-sizes compared to the smallest step size (10 −4 µm).We found that a 10-µm step size appeared to be a good trade-off, as the estimated errors were less than 10 −9 µm for intersections, 10 −14 for cosine data, and less than 10 −6 µm for OPD, with a computation time less than 0.2 seconds on an Intel Xeon processor at 3GHz.

Limits of the reconstruction algorithm based on laser ray tracing and OCT input data
Figure 3 shows the RMS of the difference between the nominal and reconstructed GRIN (GRIN RMS difference) for the three GRIN models studied and for the five proposed experimental configurations.The simulations were conducted assuming no input error and with added Gaussian errors (10 −6 to 10 −3 ) to the input data (for all conditions).In absence of experimental errors, the GRIN was reconstructed with high accuracy (GRIN RMS difference<10 −8 ) indicating that the global minimum can be theoretically retrieved for all the GRIN models and all conditions under study.As expected, the presence of experimental errors increases the reconstruction error.Also the reconstruction accuracy decreases with the complexity (number of variables) of the lens model.

Estimation of the experimental errors
The experimental set-ups described in section 2.1 (Laser Ray Tracing system and OCT) were used with artificial lenses to obtain estimates of the experimental errors of the input data.
Figure 4 shows the differences, for all rays, between the experimental and the simulated direction cosines of the rays deflected by the lens (a), impacts on a plane after the lens (b) and OPD (c) data for the artificial lenses under study.The standard deviation of the differences in each case was: 1.08 ± 0.49•10 −3 for the homogeneous lens and 0.92 ± 0.14•10 −3 for the GRIN lens (direction cosines); 6.0 ± −2.2 µm for the homogeneous and 6.9 ± −0.3 µm for the GRIN lens (ray impacts); and 2.40 ± 0.17 µm and 2.74 ± 0.18 µm (OPDs) for the two lenses under test.

Influence of experimental errors on the GRIN reconstruction
The reconstruction algorithm was evaluated for different sets of simulated input data representing the 5 proposed experimental configurations (direction cosines of deflected rays, direction cosines and intercept of the outgoing rays with the posterior lens surface, impacts on a plane after the lens, optical path up to the posterior surface of the lens, and optical path up to the posterior surface of the lens and to a plane after the lens) with three different levels of experimental error levels and for the three Goncharov models.
For illustration purposes, we defined three different error levels: errors within the expected order of magnitude of the experimental measurements (error level R), lower (error level L) and higher (error level H) than the expected experimental error.
The simulated errors of the ray direction cosines were 0.5, 1 and 2 • 10 −3 for the errors levels L, R and H respectively; the simulated errors of the intersection points (intercept of the outgoing ray with posterior lens surface and impact on a plane after the lens) were 3, 6 and 12 µm for the errors levels L, R and H respectively; and the simulated errors of the OPD 1.5, 3 and 6 µm for the errors levels L, R and H respectively Figure 5 shows the mean value and the standard deviation of the GRIN RMS difference of the nominal and the reconstructed GRIN for 50 realizations of the reconstruction algorithm.
We found that, for many conditions, the reconstruction error increases with GRIN model complexity.For the three GRIN models, the best reconstructions were achieved using the input OPD data from OCT (with the lowest error occurring for the condition that uses OPD data of both the posterior lens surface and cuvette).Increasing the error of the input data (L, R and H) increases the reconstruction error.Whereas all configurations induced low reconstruction errors (GRIN RMS difference<0.005)for the two-variable Goncharov model, for three and four-variable Goncharov models the errors were higher when using the cosine and impacts as input data.The error of the reconstructed refractive index is of the same order of magnitude in both the axial and meridional planes, although the error is slightly higher along the axial coordinate.Also, the error tends to increase towards the surface for most configurations.

Influence of the surface shape measurement errors in the GRIN reconstruction
We found discrepancies in the measured radius of curvature of 1.02 ± 1.16% and 1.05 ± 0.54% for the two measured lenses.These results are in good agreement with a previous study [39], where the error was found to be around 1%.
We estimated the relative contribution to the GRIN reconstruction induced by the surface shape measurement error (estimated 1%) and the experimental input data errors (level R).Comparisons were made between the results of the simulations with and without errors in the surface shape measurement.Figure 7 shows the GRIN RMS difference for a 4-variable GRIN model with error in both the surface and input data, relative to the reconstructed GRIN with only input data errors.We found that the relative contributions of the surface errors in the reconstructed GRIN were significant (p<0.01) in the configurations using direction cosines and intercept of the outgoing rays with posterior lens surface points (35%), OPD up to the posterior surface of the lens (20%) and OPD up to the posterior surface of the lens and cuvette (40%).

Influence of the ray sampling density on the reconstructed GRIN
We studied the effects of changing the number of input data in the reconstruction algorithm, or equivalently, increasing the corresponding number of rays in the laser ray tracing sampling pattern or increasing the number of A-scans in a cross-section of the OCT images.
While the pupil diameter was kept constant (6 mm), 7 different ray sampling densities were studied ranging from 6 to 1200 total rays, i.e. ray separation distance between 1 mm and 5 µm respectively.All simulations were performed for an error level R, and for the 4-variable Goncharov model, for the 5 proposed experimental configurations.
Figure 8 shows the GRIN RMS difference as a function of number of rays traced.The reconstruction improves as the number of rays increases, although beyond 100 rays the changes are minor.The impact of the number of rays on the accuracy of the reconstruction is lower for the OCT-based OPD input data than for the laser ray tracing-based input data.

Discussion
We have presented a GRIN reconstruction method based on optimization techniques applied to laser ray tracing (ray direction cosines and impacts) and optical coherence tomography input data.We have studied the GRIN reconstruction accuracy and the influence of the experimental errors of the different input data.
In the absence of experimental errors the algorithm has been proved to be sufficiently robust to reconstruct GRIN data for different GRIN models.In the presence of the estimated experimental errors (obtained using custom-developed laser ray tracing and OCT experimental devices), OCT-based data allowed highly accurate reconstructions (GRIN RMS difference<0.005) of the GRIN.Also, the use of OCT-based input data appears less susceptible to errors in the surface radii of curvature, and also on the number of rays.This indicates that OPD is a suitable method for retrieving the GRIN of the crystalline lens.The superiority of this method over the rest proposed in this study (using direction cosines, intercept of the outgoing rays with posterior lens surface or impacts in a plane after the lens), might result from the more direct information on the GRIN that the OPD accumulates at the end of the ray trace.The accuracy of the reconstruction increased when incorporating additional information -at least for the levels of error studied-to the input data: adding the intercept of the outgoing ray with the posterior lens surface in addition to the direction cosines in the ray tracing procedure (configuration 2 versus 1), and adding cuvette OPD in addition to the lens OPD in OCT (configuration 5 versus 4).
To test the flexibility of the reconstruction algorithm, we used three Goncharov models [36] with increasing number of variables.The accuracy of the reconstruction decreased slightly when increasing the complexity of the model.The reconstruction algorithm can be used with other GRIN models.In a previous work [20], using a different 4-variable exponential GRIN model [32] we found similar accuracy (GRIN RMS difference of 4•10 −3 ) using OCT-based input data (posterior lens and cuvette distortions), as in configuration 5 of the current study, and an estimated error of 5 µm in the OPD.
A common problem to all optimization techniques is local minimum trapping, and the possible existence of multiple solutions.Since the merit function used is not described analytically, the absence of local minima or the existence of multiple solutions cannot be proved mathematically.However, the use of a global instead of a local search algorithm prevents local minimum trapping.The large set of multiple experimental input data and relatively low number of unknown variables reduces the possibility for several solutions with the same Merit Function value.
The reconstruction algorithm and the error analysis for the various experimental approaches presented can be compared to those of previous studies.Acosta and Vazquez [12] proposed a reconstruction algorithm based on laser ray tracing at different orientations to reconstruct a non-spherical GRIN.They simulated a Gaussian error of less than 1 µm.With data from orientations up to 80°, the RMS of the GRIN reconstructed was below 10 −4 for a GRIN described by a 9-variable single polynomial and around 10 −3 for a GRIN described by two polynomial expressions (strong GRIN).The RMS reconstruction error, if only orientations up to 10° were available, was >10 −3 in both cases.In comparison, our reconstruction method, when using the same simulated input error, achieved a reconstruction error <5•10 −3 for the 4-variable GRIN model, which is comparable to the results by Vazquez et al (with input data up to 10°).It is worth noting that our reconstruction algorithm could be extended for laser ray tracing data at different orientations.Although we did not explore that possibility in this work, it is likely that the accuracy would increase further by incorporating input data obtained at different orientations.Verma et al. [19] used OCT-based input data and an optimization algorithm (nonlinear least squares fitting) to reconstruct a spherical GRIN lens.Simulating a Gaussian noise of 11 µm they obtained a maximum error of around 0.013.We reproduced their simulations with our local search algorithm finding similar results (maximum error 0.010 and GRIN RMS difference of 5•10 −3 ).
There are other possible sources of experimental error, not addressed here, such as decentration or tilt of the lens, or inaccurate positioning of the ray entrance, which would affect the outcomes.However, several of these errors can be minimized when obtaining threedimensional input data (particularly in the OCT based technique, where the surfaces of the lens are visualized directly).As it was mentioned in section 2.3, the refractive index retrieved from the reconstruction algorithm in configurations 4 & 5 is assumed to be equivalent to the group refractive index.The resulting index can be converted from the group to the phase refractive index at the central wavelength of the OCT light source and subsequently to a different wavelength using available data of the chromatic dispersion of the crystalline lens [48,49].
The techniques described here are designed for estimations of the GRIN of the crystalline lens in vitro.However, the estimation of the GRIN distribution of the lens in vivo would be of major interest to understand the contributions of GRIN to the optical properties of the eye, and its influence in the optical changes with accommodation and aging, particularly considering that the shape of the isolated crystalline lens differs from its unaccommodated state in vivo.To date, all the laser ray tracing GRIN reconstruction methods have proposed imaging the lateral ray deflections, which are unavailable in vivo.In contrast, we have demonstrated that relatively good reconstruction performances (GRIN RMS difference<0.01)can be obtained when using transverse imaging of the impacts after the lens.These impacts can actually be available in vivo (with additional contributions of the cornea) in a double-pass configuration as routinely shown in laser ray tracing (LRT) measurements of the wave aberrations in the eye [50,51], although the algorithm also requires knowledge of the posterior shape of the crystalline lens.Nevertheless, while the OCT-imaging-based technique described here works with in vitro samples, it is likely that this technique, in combination with an in vivo LRT-based technique, operating at wide angles if possible, would provide sufficient information for attempting a reconstruction of the GRIN in vivo through optimization techniques.

Fig. 1 .
Fig. 1.(a) Experimental setup for lateral ray tracing system (b) Corresponding ray images (integrated image of five rays) (c) Single-pass ray tracing configuration to measure the spot diagram in a plane after the lens.(d).Corresponding spot images for five rays.

Fig. 2 .
Fig. 2. (a) Estimated error in simulated data (lateral deviations, impacts or OPD) as a function of the step size in the Sharma algorithm.(b) Computational time to ray trace and to estimate the OPD as a function of the step size.The OPD was calculated either using straight segments or a Hermite polynomial interpolation.Simulations were performed with Goncharov crystalline lens model 20S.

Fig. 3 .
Fig. 3. Difference between nominal and reconstructed GRIN (GRIN RMS difference) as a function of the experimental error of the input data.Data are the mean across 50 repetitions, and the error bars represent standard deviations.

Fig. 4 .
Fig. 4. Difference between experimental and simulated data: Mean and standard deviation for five repeated measurements.(a) Direction cosines of the outgoing rays, (b) Impacts on a plane after the lens and (c) OPD through the lens.

Fig. 5 .
Fig. 5. Difference between reconstructed and nominal GRIN (GRIN RMS difference) for the 5 proposed experimental configurations and three different levels (L R and H) of added experimental error in the input data (L: lower than expected errors; R: realistic errors; H: higher than expected errors).Data are the average GRIN RMS difference and the error bars the standard deviation of 50 realizations of the GRIN reconstruction algorithm.(a), (b) and (c) are the results for three different Goncharov models with increasing complexity (20U, 20B and 20S).

Figure 6
Figure6shows the spatial distribution of the reconstruction error in the different conditions, across the axial and meridional axes for the 4-variable GRIN model and the error level R.The error of the reconstructed refractive index is of the same order of magnitude in both the axial and meridional planes, although the error is slightly higher along the axial coordinate.Also, the error tends to increase towards the surface for most configurations.

Fig. 6 .
Fig. 6.Axial (a) and Meridional (b) deviation from the nominal GRIN profile for the proposed experimental configurations, for realistic input data error levels (R) and the 4-variable Goncharov model.Data represent the mean value and the error bars the standard deviations of 50 realizations of the reconstruction algorithm.

Fig. 7 .
Fig. 7. GRIN RMS difference of the reconstructed and the nominal GRIN data, with realistic error level (R) in the input data without error in the surface shape (solid colors), and with a random deviation of the surface geometry of 1% (light colors).The bars represent average data and the error bars, the standard deviation of 50 realizations of the reconstruction algorithm.

Fig. 8 .
Fig. 8. GRIN RMS difference of the reconstruction for the 5 proposed experimental configurations and with realistic error level R, versus number of rays.Pupil radius was set to 6mm pupil.Data represents mean value, and the error bars, the standard deviation of 50 repetitions of the reconstruction algorithm.