Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging

We present an optimization method to retrieve the gradient index (GRIN) distribution of the in-vitro crystalline lens from optical path difference data extracted from OCT images. Three-dimensional OCT images of the crystalline lens are obtained in two orientations (with the anterior surface up and posterior surface up), allowing to obtain the lens geometry. The GRIN reconstruction method is based on a genetic algorithm that searches for the parameters of a 4-variable GRIN model that best fits the distorted posterior surface of the lens. Computer simulations showed that, for noise of 5 μm in the surface elevations, the GRIN is recovered with an accuracy of 0.003 and 0.010 in the refractive indices of the nucleus and surface of the lens, respectively. The method was applied to retrieve threedimensionally the GRIN of a porcine crystalline lens in vitro. We found a refractive index ranging from 1.362 in the surface to 1.443 in the nucleus of the lens, an axial exponential decay of the GRIN profile of 2.62 and a meridional exponential decay ranging from 3.56 to 5.18. The effect of GRIN on the aberrations of the lens also studied. The estimated spherical aberration of the measured porcine lens was 2.87 μm assuming a homogenous equivalent refractive index, and the presence of GRIN shifted the spherical aberration toward negative values (−0.97 μm), for a 6-mm pupil. ©2010 Optical Society of America OCIS codes: (330.0330) Vision, color and visual optics; (080.2710) Inhomogeneous optics media; (170.4500) Optical coherence tomography; (110.6880) Three-dimensional image acquisition. References and links 1. B. K. Pierscionek, “Presbyopia effect of refractive index,” Clin. Exp. Optom. 73(1), 23–30 (1990). 2. W. S. Jagger, “The optics of the spherical fish lens,” Vision Res. 32(7), 1271–1284 (1992). 3. J. W. Blaker, “Toward an adaptive model of the human eye,” J. Opt. Soc. Am. 70(2), 220–223 (1980). 4. G. Smith, and B. K. Pierscionek, “The optical structure of the lens and its contribution to the refractive status of the eye,” Ophthalmic Physiol. Opt. 18(1), 21–29 (1998). 5. 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). 6. 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). 7. 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). 8. O. Pomerantzeff, H. Fish, J. Govignon, and C. L. Schepens, “Wide-angle optical model of the eye,” Am. J. Optom. Physiol. Opt. 19, 387–388 (1972). 9. 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). 10. B. K. Pierscionek, and D. Y. Chan, “Refractive index gradient of human lenses,” Optom. Vis. Sci. 66(12), 822– 829 (1989). #132469 $15.00 USD Received 29 Jul 2010; revised 7 Sep 2010; accepted 13 Sep 2010; published 30 Sep 2010 (C) 2010 OSA 11 October 2010 / Vol. 18, No. 21 / OPTICS EXPRESS 21905 11. 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). 12. 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). 13. M. C. Campbell, “Measurement of refractive index in an intact crystalline lens,” Vision Res. 24(5), 409–415 (1984). 14. 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). 15. G. Beliakov, and D. Y. Chan, “Analysis of Inhomogeneous Optical Systems by the Use of Ray Tracing. II. Three-dimensional Systems with Symmetry,” Appl. Opt. 37(22), 5106–5111 (1998). 16. E. Acosta, D. Vázquez, 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). 17. D. Vázquez, 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). 18. C. E. Jones, and J. M. Pope, “Measuring optical properties of an eye lens using magnetic resonance imaging,” Magn. Reson. Imaging 22(2), 211–220 (2004). 19. 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). 20. Y. Verma, K. Rao, M. Suresh, H. Patel, and P. Gupta, “Measurement of gradient refractive index profile of crystalline lens of fish eye in vivo using optical coherence tomography,” Appl. Phys. B 87(4), 607–610 (2007). 21. 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). 22. 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). 23. D. Siedlecki, A. de Castro, S. Ortiz, D. Borja, F. Manns, and S. Marcos, “Estimation of contribution of gradient index structure to the amount of posterior surface optical distortion for in-vitro human crystalline lenses imaged by Optical Coherence Tomography,” presented at the Fifth European Meeting on Visual and Physiological Optics (2010). 24. 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). 25. F. Manns, A. Ho, D. Borja, and J. Parel, “Comparison of Uniform and Gradient Paraxial Models of the Crystalline Lens,” Invest. Ophthalmol. Vis. Sci. 51, E-Abstract 789 (2010). 26. A. M. Rosen, D. B. Denham, V. Fernandez, D. Borja, A. Ho, F. Manns, J. M. Parel, and R. C. Augusteyn, “In vitro dimensions and curvatures of human lenses,” Vision Res. 46(6-7), 1002–1009 (2006). 27. J. H. Holland, Adaptation in natural and artificial systems (The University of Michigan Press, 1975). 28. O. N. Stavroudis, The optics of rays, wavefronts and caustics (New York Academic Press, 1972), pp. 81–95. 29. 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). 30. 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). 31. J. A. Nelder, and R. Mead, “A simplex method for function minimization,” Comput. J. 7, 308–313 (1965). 32. M. Dubbelman, and G. L. Van der Heijde, “The shape of the aging human lens: curvature, equivalent refractive index and the lens paradox,” Vision Res. 41(14), 1867–1877 (2001). 33. M. Dubbelman, G. L. van der Heijde, and H. A. Weeber, “The thickness of the aging human lens obtained from corrected Scheimpflug images,” Optom. Vis. Sci. 78(6), 411–416 (2001). 34. C. E. Jones, and J. M. Pope, “Measuring optical properties of an eye lens using magnetic resonance imaging,” Magn. Reson. Imaging 22(2), 211–220 (2004). 35. 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). 36. S. Ortiz, D. Siedlecki, L. Remon, and S. Marcos, “Optical coherence tomography for quantitative surface topography,” Appl. Opt. 48(35), 6708–6715 (2009). 37. R. J. Zawadzki, C. Leisser, R. Leitgeb, M. Pircher, and A. F. Fercher, “Three-dimensional ophthalmic optical coherence tomography with a refraction correction algorithm,” Proc. SPIE 5140, 8 (2003). 38. A. Roorda, and A. Glasser, “Wave aberrations of the isolated crystalline lens,” J. Vis. 4(4), 250–261 (2004). 39. 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). 40. S. G. el-Hage, and F. Berny, “Contribution of the crystalline lens to the spherical aberration of the eye,” J. Opt. Soc. Am. 63(2), 205–211 (1973). 41. J. G. Sivak, and R. O. Kreuzer, “Spherical aberration of the crystalline lens,” Vision Res. 23(1), 59–70 (1983). 42. P. Artal, E. Berrio, A. Guirao, and P. Piers, “Contribution of the cornea and internal surfaces to the change of ocular aberrations with age,” J. Opt. Soc. Am. A 19(1), 137–143 (2002). #132469 $15.00 USD Received 29 Jul 2010; revised 7 Sep 2010; accepted 13 Sep 2010; published 30 Sep 2010 (C) 2010 OSA 11 October 2010 / Vol. 18, No. 21 / OPTICS EXPRESS 21906 43. S. Barbero, S. Marcos, and J. Merayo-Lloves, “Corneal and total optical aberrations in a unilateral aphakic patient,” J. Cataract Refract. Surg. 28(9), 1594–1600 (2002).


Introduction
The understanding of the optical quality of the eye, and its changes with accommodation and age requires knowledge of the geometry and optical properties of the crystalline lens.Due to its complicated structure and its importance in the accommodation process, the lens is an interesting element of the optical system of the eye, however it is less accessible than the cornea through optical techniques.In addition, the posterior lens surface is distorted by the refraction of the preceding optical surfaces and by the non-homogeneous refractive index of the crystalline lens, known to have a gradient distribution (GRIN).The presence of GRIN in the lens plays a critical role in the crystalline lens power [1], and very likely also in its high order aberrations.In some animals, such as in fish, it is well known that the GRIN compensates for a large amount of the spherical aberration induced by the spherical geometry of the lens [2].In humans and other species, although there have been many attempts to incorporate GRIN in schematic model eyes [3][4][5][6], these in general fail to explain the change both in power and spherical aberration with age and accommodation, presumably because of the lack of validations of the GRIN models used, and the use of average geometrical data and ocular aberrations from different literature sources (and not from the same eye).
There have been various reports in the literature on the use of non-destructive techniques to measure GRIN in vitro, or to estimate it through limited optical data in vivo.A Purkinje imaging technique for lens phakometry has been used to estimate the crystalline lens GRIN [7], however the model used was very simple and the method required assumption of the value of the refractive index in the nucleus.Optimization approaches using shell models have also been proposed [8,9].For example, Pomerantzeff attempted to estimate the GRIN structure in vivo by fitting a shell model to the spherical aberration of the eyes of 50 subjects.However, the model included multiple unknown variables such as lens thickness, curvatures and refractive indices of the different layers making the problem largely multidimensional.In addition to the problem of multiple solutions matching the spherical aberration of the eye [10] a local search algorithm was used, which very likely missed a global solution.The use of local search optimization algorithms allowed retrieving the GRIN of in vitro spherical fish lenses [11].In addition, we had proposed the use of global search algorithms to retrieve the GRIN parameters of a commercial artificial GRIN lens [12].
Ray tracing has also been applied to estimate the crystalline lens GRIN structure in isolated lenses in combination with measurements of the lens shape.Campbell [13] applied a ray tracing method to reconstruct the GRIN of a rat lens, but the method needed matching between external media and surface refractive index, as well as use of an elliptical isoindicial surfaces GRIN model.This method was improved by Chan [14] and Beliakov et al. [15] to avoid the use of elliptical isoindicial surfaces.Acosta et al. [16,17] presented an algorithm, also based on ray tracing, using bipolynomial GRIN distribution for the lens.The method required measuring the deflections of a bundle of rays at different angles (tomographic reconstruction).These methods relate mathematically the exit point and deflection of the rays with the GRIN distribution.While computer simulations predict that the parameters of the model can be obtained with high accuracy (even in the presence of simulated experimental noise), the technique is ultimately limited by the precision of the detection of the point where each ray intersects both surfaces of the lens and therefore is only applicable if lateral ray tracing data are available.In addition to optical methods, magnetic resonance imaging (MRI) has been used to estimate the human crystalline lens GRIN distribution in vitro [18], and in fact, results have been even presented in vivo [19].Limitations of the technique include its relatively low resolution, and the assumption of the direct relationship between water content and refractive index of the lens.
To our knowledge, only one study has attempted to use optical coherence tomography (OCT) to reconstruct experimentally the GRIN in a fish lens [20].In an earlier study, we had analyzed theoretically the potential of OCT to retrieve the GRIN in human lenses [21].The OCT measures optical path differences, and its signal carries intrinsically information on the refractive index.The application of an optimization method in the fish lens to reconstruct the GRIN from OCT images is relatively simple, as the spherical shape of the lens allows assuming full knowledge of the surface geometry and a very simple GRIN model which avoids local minimum problems in the optimization process.
The distortion produced by the isolated crystalline lens on the visualization of a plane surface (the base of the cuvette holding the lens) imaged with OCT technique has been previously used to estimate the physical thickness of the lens as well as the average refractive index of human lenses of different ages [22].In a recent study with isolated human lenses we have shown that the distortion of the posterior surface of the crystalline lens is in part produced by the presence of GRIN (which primarily affects the estimates of posterior lens asphericity) [23].The comparison with the actual shape of the posterior lens surface was possible by imaging the lens in two orientations (posterior lens surface up and down).We have proposed a ray-tracing method for the correction of optical distortion in OCT imaging [24], and in particular, in the presence of GRIN [23].The fact that the GRIN produces a significant distortion in the posterior lens surface and in the cuvette holding the lens can be actually used to attempt the reconstruction of the GRIN, provided that the lens can be measured in two orientations, as it occurs in vitro.
In this study we present a reconstruction method of the GRIN based on optimization algorithms applied to crystalline lenses with non-spherical geometry and a GRIN structure similar to that expected to be found in the human crystalline lens.The input data are the optical path differences measured in OCT through the lens, and the shape of the anterior and posterior crystalline lens surfaces (with the posterior directly measured in vitro by flipping the lens), also obtained from OCT.The anatomical GRIN model used is a 4-variable model proposed by Manns et al. [25].The method is demonstrated computationally and tested experimentally on a porcine lens in vitro.Three-dimensional OCT measurements allowed a full 3-D reconstruction of the GRIN structure, showing asymmetries in the GRIN distribution and avoiding potential artifacts arising from measuring GRIN in a single meridian.Knowledge of both the surface geometry and GRIN allowed estimating the relative contribution of each to key optical aberrations (i.e.astigmatism and spherical aberration).

GRIN model
A 4-variable model was used to describe the crystalline lens GRIN structure.The model has been described in detail elsewhere [25].In brief, the anterior and posterior surfaces of the crystalline lens are described by conics, the GRIN follows the surface shapes, with an exponential variation in the refractive index from the nucleus to the surface in both the axial and meridional directions.The center of the GRIN lies is in the intersection of the optical axis and the equatorial plane.We implemented this model, setting the equatorial plane at a distance from the anterior surface vertex equal to 0.41 times the lens thickness [26].The GRIN can be expressed in polar coordinates centered in the lens as where ρ S is the distance from the center of the lens to the surface at angle θ; p(θ) is a monotonic function which expresses the change of the exponential decay value from the axial direction in θ = 0, p 1, to the meridional direction, θ = π / 2, p 2. This function was chosen to depend on the fifth power of θ for this work; n N is the refractive index in the nucleus of the lens, and ∆n the difference between surface and nucleus index

GRIN reconstruction algorithm
We developed an optimization algorithm that searches for the GRIN parameters by minimizing a Merit Function defined as: ), OPD (Surf3) , where MF stands for the Merit Function as a function of the 4-variables of the GRIN (GRIN_coeffs), the RMS refers to the root mean square, OPD calc is the estimated optical path difference for the posterior surface (Surf2) and the cuvette surface (Surf3), and OPD exp the optical path difference measured by OCT.The Merit Function is therefore defined by the sum of RMS of the differences between the estimated and measured optical path differences for the posterior surface of the crystalline lens and the cuvette.The search algorithm is based on a genetic optimization [27].
A ray tracing algorithm was programmed in MatLab (MathWorks, Natick MA) using Stavroudis formulas [28] to implement the vectorial Snell's law, Sharma algorithm [29] for ray tracing in a GRIN structure, and Stone et al.'s method [30] to calculate the optimal interpolant of the last step of Sharma algorithm to locate the intersection of the ray with the posterior surface of the lens.
The algorithm starts with a population of solutions randomly selected within the expected range of solutions.The starting population ranged from 1.36 to 1.39 for the surface refractive index (n S ), 1.38 to 1.42 for the nucleus refactive index (n N ) and 2 to 9 for the optical and meridional exponential decay (p1 and p2).The algorithm was set without constrains, i.e., there were not limits in the search of variable values.A schematic diagram of the rules of the algorithm is presented in Fig. 1.In brief, the Merit Function is evaluated in each solution, and then a second population is created by selecting some of the solutions stochastically (taking into account the merit function value) and recombining (85%) two solutions (i.e.combining their variables) and mutating them (75% of recombined children) by adding gaussian noise.Some of the best solutions of a previous iteration (5%) are directly copied to the next iteration and a few (10%) are chosen randomly.With this mixture of determinism and randomness we assured that the algorithm explored the entire solution space, avoiding local solutions by not converging in the first iterations to a local minimum.In the simulations presented here the population for each iteration was formed by 400 solutions and the process was repeated 10 times.A local search algorithm [31] was then applied starting with the genetic optimization solution.The computation time was less than 45 minutes for each complete optimization using a 10 µm step in the Sharma algorithm.We tested that using different Sharma's algorithm steps below that value does not produce significant changes in the ray tracing algorithm or in the optimization results.Fig. 1.Schematic diagram of a genetic algorithm iteration.In this study, 85% of the population is created by crossover and mutation.The 5% best solutions are directly copied to the next solution (elite) and the rest of the population is randomly selected from previous generations.

Simulations
The performance of the GRIN reconstruction method was studied simulating three human crystalline lenses (20, 40, and 60 years old), with realistic geometry and a GRIN distribution consistent with previous literature.The anterior and posterior lens surface shapes and lens thickness were obtained from Scheimpflug studies [32,33] and the parameters for the GRIN model were obtained from MRI studies [34].The nominal refractive indices were assumed to be 1.378 (surface index) and 1.410 (nucleus index).The axial exponential decay (p1) were chosen to be 3 (young lens), 4 (middle age lens) and 6 (old lens); and the meridional exponential decay (p2), 2(young lens), 4 (middle age lens) and 8 (old lens).These values represent the commonly assumed GRIN profiles of a young human lens (distributed GRIN from the nucleus to the surface in the young lens), and an old human lens (a large plateau and rapid index change toward the surface) [34].Figure 2 presents assumed GRIN distributions for young (A), middle age (B) and old (C) human lenses used in simulations.
The optical path difference of each ray was obtained by simulating the distortion of the posterior surface of the lens and the cuvette produced by refraction and was calculated using the ray tracing algorithm (see section 2.2) over a 6 mm pupil.Gaussian noise was added to simulate inaccuracies in the detection of the surfaces from OCT in absence of other sources of experimental error.The standard deviations of the noise used in the simulations ranged from 0.1 (well below the expected experimental values) to 20 µm (much higher than the OCT resolution), with intermediate values of, 0.1, 1, 5, and 10 µm.Simulations of the GRIN reconstruction using the global search algorithm described in Section 2.2 were performed, with 100 repetitions in each condition.The goodness of the reconstruction was assessed by comparison of the nominal and reconstructed GRIN parameters, and in terms of root mean square (RMS) of the differences between nominal and recovered GRIN in a grid of points over the lens.

Experimental measurements
Experimental measurements were performed on an isolated porcine lens.Enucleated porcine eyes were obtained in a local slaughterhouse, and used within 4 hours post-mortem.The lens was extracted from the cadaver eye and placed on a ring holder in a PMMA chamber filled with preservation medium (BSS plus, Alcon, Fort Worth TX, n = 1.345).The lens was aligned with the OCT system and 3-D OCT images were collected of the anterior and posterior lens surfaces, and the base of the cuvette.The lens was first imaged with the anterior surface up, and then flipped over and imaged with the posterior surface up (similarly to the procedure described [24], but with 3-D rather than 2-D data acquisition).Also, the cuvette was imaged with and without the crystalline lens filled in with preservation medium in both cases.The whole procedure of the measurement took less than one hour, during which the crystalline lens was continuously immersed in the preservation medium.

OCT system
The OCT images were obtained using a custom-developed high resolution spectral domain wich uses a SLD diode of 840 nm as central wavelength and a FWHM of 50 nm, the SNR of the instrument was estimated on 97 dB, the adquisition speed was 25,000 A-Scans/sec resulting in an integration time of 40 µm.More details can be found in previous works [24,35].The system was adapted to the measurement of horizontally resting crystalline lens by inserting a mirror at 45 deg in front of the last lens of the instrument and on top of the cuvette holding the crystalline lens.An unfolding mode [35] was used to expand the axial range of the system up to 10 mm in air.In addition, to ensure optimal signal from both lens surfaces, two 3-D images were obtained per condition at two different focal planes, and then merged using the cuvette as a reference, to produce a full 3-D image of the crystalline lens.The merging algorithm was specifically developed in Matlab for this study.A single 3-D image consisted of 1668 A-scans and 70 B-scans in a lateral range of 12x12 mm with a resolution of 170x7 µm.The axial resolution of the images was 3.42 µm, and the acquisition time of a full 3-D image was 4.5 seconds.Figure 3 shows a full 3-D view of the crystalline lens acquired as described in the text.

OCT image processing
OCT images were corrected for fan distortion, which arises from the scanning architecture, using algorithms developed in a previous study [36].The fan distortion correction and posterior image analysis requires segmentation of the anterior and posterior crystalline lens surface, as well as the base of the cuvette.The detection of anterior and posterior surfaces of the crystalline lens was performed using two different strategies, according to the differences in the reflectivity of the surfaces.For the anterior surface, the segmentation algorithm searches the maximum of intensity in each A-scan, and for the posterior surface the segmentation algorithm is based on edge detection of the thresholded image.The surface of the distorted cuvette was segmented using the first method.
The pair of 3-D full crystalline lens images obtained in two orientations (anterior surface and posterior surface up) was registered, to ensure that comparisons were made for the lens in an identical position.The axis of astigmatism of the entire lens (which should be the same in the two orientations of the lens) was used as a reference for alignment.The axis of astigmatism of each lens surface was computed from the Zernike fits (up to the 6th order, 6mm diameter) to each surface elevation map.The astigmatism axis of the lens was calculated from the astigmatic distortion of the image of the base of the cuvette.Using astigmatism as a signature for rotation leaves an uncertainty of 180 degrees, which is easily solved by carefully monitoring the lens handling while the lens was flipped, and by inspection of irregularities in the image that serve as landmarks.Figure 4 shows the data obtained from three dimensional OCT measurements of the crystalline lens placed with the anterior surface up (A) and the posterior up (B).

Crystalline lens shape and thickness
The crystalline lens shape was obtained from the direct OCT images of the crystalline lens (anterior lens up and posterior lens up), i.e. not subject to optical distortion, and corrected from fan distortion.The 3-D elevation maps were fitted by Zernike polynomials up to 6th order.However, as the GRIN model used is defined assuming conic surfaces, only symmetric Zernike coefficients and astigmatism were used to implement the lens surfaces in the algorithm.
The crystalline lens thickness was calculated from the distortion induced on the cuvette, following a method described in an earlier study [22].

Three-dimensional experimental GRIN reconstruction algorithm
As the crystalline lens surfaces lack from rotational symmetry, it was necessary to generalize the GRIN model of Eq. ( 1) to three dimensions.The generalization was made by selecting 18 meridians (from 0° to 170° in steps of 10°), and applying the optimization algorithm described in 2.2 to a combined merit function for all meridians.The GRIN profile of all meridians had the same values for three of the variables (surface index, nucleus index and the exponential decay in the optical axis, p 1 ).The exponential decay in the meridional axis, p 2 , was left free, to account for the differences of the GRIN across meridians.

Analysis of the optical effects on GRIN
In order to study the influence of the GRIN distribution in the optics of the crystalline lens, aberrations of the entire lens were calculated with the experimental GRIN and with a homogeneous equivalent index (defined as the homogeneous index that makes a lens with the same surface geometry match the focal length of the lens with the GRIN).The geometry of the lens was described by the Zernike fits to the raw data, as described before.A standard computational ray tracing analysis was followed to estimate the lens aberrations.The plane of best focus was calculated searching the plane where the deviations of ray impacts (RMS) from the optical axis reached a minimum.The aberrations were calculated by fitting the optical path differences of rays (at best focus) with respect to central ray to a 6th order Zernike polynomial.Calculations were done for 6-mm pupils and 200 rays in each meridian.

Simulations: Predicted performance of the GRIN reconstruction algorithm
Figure 5 shows the goodness of the GRIN reconstruction in simulated crystalline lenses as a function of the simulated noise in the lens surface detection.Without errors, the exact nominal parameters of the GRIN were obtained.The increase in surface detection noise increased the error in the reconstructed GRIN, as well as the standard deviation of the retrieved parameters (across repetitions).The mean values and standard deviation for all the ages for a 5 µm error in the data are 0.010 ± 0.004 for the surface index, 0.003 ± 0.001 for the nucleus index and 1.4 ± 0.5 and 2.2 ± 0.5 for the exponential decay in optical and meridional axis values.As shown in Fig. 5, no large differences were found in the quality of the reconstruction of the simulated lenses across ages, although for small errors, the reconstruction is worst for the surface index and best for the nucleus index in older lenses.This is possibly due to the larger role of the nucleus in the description of the lens GRIN in old lenses in comparison to young lenses.The error in the global description of the GRIN in terms RMS of is 0.004 ± 0.001, and practically constant across ages.

Experimental results: Reconstruction of the porcine crystalline in 3-D
The GRIN structure in the porcine lens was obtained by applying the global search algorithm, as described in Section 2.3.The optimization search algorithm was run 5 times for the same set of data.
The retrieved refractive index in the nucleus, n N was 1.443, and the index of refraction in the surface, n S was 1.362.The axial exponential decay, p 1 , was 2.62, and the meridional exponential decay, p 2 , varied from 3.56 to 5.18 depending on the meridian angle.
The standard deviation of these parameters across repetitions was 0.0007 and 0.003 for the nucleus and surface index respectively, 0.09 for p1, and 0.30 for p2 (averaged across meridians).
The change in the meridional exponential decay p 2 is shown in Fig. 6, in comparison with the radius of curvature of the anterior and posterior lens surface, as a function of the meridional angle.Changes in these parameters across angle are indicative of astigmatism.The anterior surface of the lens appears more astigmatic than the posterior surface in this eye, with the steepest meridian at about 50 deg.Interestingly, p 2 peaks at this angle, making the crystalline lens less convergent in this meridian, and therefore producing a compensation of the surface astigmatism.In fact, in this particular lens, GRIN overcompensates the astigmatism induced by the surfaces.The actual values of astigmatism are reported in section 3.3.Gradient index structure for different meridians is shown in Fig. 7.

Discussion
We have proposed a new method to reconstruct the GRIN of the crystalline lens in vitro, based on OCT imaging and a global search algorithm.The method is based on the retrieval of the GRIN profile from the optical distortions produced on the posterior surface of the lens by the anterior surface and the GRIN.Previous works proposed methods for compensation of the optical distortion in OCT [24,37].Here, we take advantage of this distortion to reconstruct the GRIN.
Computer simulations show that, provided that the lens surfaces are known within certain accuracy, it is possible to retrieve the GRIN from the optical path difference of each ray, which is the information provided in OCT.For the detection errors and resolution of state-ofthe-art OCT, we predict an accuracy in the GRIN reconstruction of 0.004 RMS.The algorithm is easily implemented in 3-D, and we have achieved, for the first time to the best of our knowledge, a 3-D reconstruction of the gradient index distribution in a porcine lens.
A critical aspect of the GRIN reconstruction is the model used to describe the index structure within the lens.We used a 4-variable GRIN model, which is sufficiently simple to allow an appropriate convergence of the optimization algorithm, while at the same time it is sufficiently flexible to accommodate differences in the variation of the refractive index from the nucleus to the surface (expected to vary strongly with age in humans) and across meridians (important in non-rotationally symmetric lenses such as the porcine lens studied).The GRIN model used is defined for a lens with biconic surfaces shapes, and therefore cannot capture all the potential complexity of surface irregularities and refractive index inhomogeneities suggested by direct measurements of the high order aberrations in porcine crystalline lenses [38].However, our estimates of spherical aberration of the porcine lens (from the geometrical and GRIN data) show negative values in agreement with aberrometry measurements.A 3-variable model such as Goncharov's balanced model [39] worked with similar performance in simulations but the lack of rotational symmetry in the crystalline lens required an additional parameter.The optimization method was not successful with polynomial GRIN models such as the Gullstrand's model, where coupling between coefficients and the high number of variables did not allow an easy relation between changes in the polynomial variables and changes in GRIN profile.Despite the difference in the GRIN model and procedure, we retrieved very similar values for the surface and nucleus refractive indices as those previously reported in porcine crystalline lenses by Vazquez et al using a tomographic ray tracing algorithm in a 2-D GRIN reconstruction [17]: 1.366 and 1.444 using a mono-polynomial model and 1.361 and 1.449 for by-polynomial model), and 1.362 and 1.443 using our OCT method and genetic algorithm in a 3-D GRIN reconstruction of a 4variable GRIN model.
The presence of the GRIN structure in the lens has a strong impact on its optical quality.Numerous studies in the literature both in vivo and in vitro have shown that the spherical aberration in the young human and primate lens tends to be negative, compensating at least in part the positive spherical aberration of the cornea [38,[40][41][42][43].We found that the spherical aberration of the isolated lens would be positive if it had a homogeneous refractive index, while the measured GRIN structure makes the spherical aberration to be negative.The estimated spherical aberration in the presence of GRIN (−0.97 µm) is in reasonable good agreement with reported measurements on an isolated porcine lenses using ray tracing (−0.7 µm) [38].By measuring the lens shape and GRIN we have been able then to identify the relative contribution of the surface shape and GRIN in the spherical aberration.As previous studies in other animals such as fish [2] we have shown a compensatory effect of GRIN also in a more complex lens structure such as the porcine lens.
The measurement of the GRIN structure in 3-D has allowed us to evaluate the influence of GRIN in other, non-rotationally symmetric terms, such as astigmatism.In the lens studied, the presence of GRIN increased the impact of the magnitude of lens astigmatism, as it would be predicted with a homogeneous index, and also plays a role in the axis of the overall lenticular astigmatism axis.
As the method proposed requires knowledge of the posterior lens shape (as well as its distortion as it is imaged through the lens), it is suited for in vitro studies, but cannot be directly applied in vivo.However, while different optimization strategies will be needed to extract the gradient index from OCT images of the crystalline lens in vivo, knowledge of the GRIN structure in vitro have allowed to identify an anatomically plausible GRIN model and can help to limit the search of the solution space in vivo.

Fig. 2 .
Fig. 2. Crystalline lens GRIN tested in the simulations for a young human lens (A), a middle age lens (B) and an old human lens (C).The index of refraction ranged from 1.410 in the nucleus to 1.378 in the surface.The black line represents the undistorted cuvette, the red line represents the distorted posterior surface, and the blue line the distorted cuvette, for a 6-mm pupil.These curves are the input data to the optimization algorithm.

Fig. 3 .
Fig. 3. (Media 1) Three dimensional image of the in-vitro posterior-up crystalline lens measured with OCT.

Fig. 4 .
Fig.4.Three dimensional OCT data from images of the crystalline lens placed with the anterior surface up (left) and the posterior up (right).The blue and red points correspond to the segmented anterior and posterior surfaces of the lens, respectively.The green points correspond to the segmented cuvette surface imaged through the crystalline lens and the black points to the segmented cuvette surface without the crystalline lens.All data are fan distortion corrected, and the surfaces are also distorted due to the presence of preservation media.

Fig. 5 .
Fig. 5. Deviation of the parameters resulting from optimization with different amounts of error in the detection of distorted surfaces.

Fig. 6 .
Fig. 6.Radius of curvature of the anterior (red) and posterior (green) lens surface, in mm. and meridional exponential decay parameter p2 (black) as a function of meridional angle.

Fig. 7 .
Fig. 7. (Media 2) Change of the gradient index of the lens with meridian angle.

Figure 8
Figure 8 shows the simulated ray tracing on the measured porcine lens, for one meridian (0 deg), with the estimated GRIN (A) and with a homogeneous equivalent refractive index (B).The presence of GRIN increases the effective astigmatism of the crystalline lens, from 1.03 D (homogeneous index) to 4.76 D (GRIN).It also shifts the axis of astigmatism from −34 deg (homogeneous index) to 61 deg (GRIN) The presence of GRIN has also a large influence on spherical aberration, shifting the 4th order spherical aberration term from positive values (2.87 µm, homogeneous index) to lower, and negative values (−0.97µm,GRIN).

Fig. 8 .
Fig.8.Ray tracing in one meridian for the crystalline lens with the reconstructed gradient index (A) and with the equivalent refractive index (B).Main differences between them are an increase in astigmatism of the lens (not visible here) and a shift of spherical aberration toward negative values.