Optical coherence tomography (OCT)^{1},^{2} has become a prominent and very useful technique in the studies of anatomy of the eye and its optical properties as well as for clinical diagnosis. Recently, the potential for quantification of OCT anterior segment images has been emphasized, with the report of new image analysis algorithms (including denoising and segmentation),^{3}^{–}^{6} the correction of motion artifacts,^{7}^{–}^{9} and the possibility for extracting biometry and topography from those images. In particular, the availability of algorithms to correct for distortions introduced by the scanning architecture of the systems (fan distortion)^{6},^{10}^{–}^{14} and for distortions arising from refraction by preceding surfaces (optical distortions)^{10},^{11},^{15} opens the possibility for full OCT-based topography of not only the anterior surface of the cornea but also internal surfaces of the ocular components.

On the other hand, optical distortion can be successfully used for gathering additional information on the sample. In fact, the presence of optical distortion in the images of the posterior surface of the crystalline lens [viewed through the anterior surface and the crystalline lens gradient index (GRIN)] has been used to reconstruct the GRIN distribution of the crystalline lens in porcine^{16} and human lenses.^{17} OCT had been applied to reconstruct the GRIN in simple spherical fish lenses^{18} and theoretically envisioned as a tool to provide information of the GRIN structure of the human lens.^{19} However, the first demonstration of the OCT capabilities to provide three-dimensional profiles of the GRIN in mammal lenses and applications in age-related studies in the human lens are very recent. The GRIN distribution in the crystalline lens has been described with various models. It seems well accepted that the index of refraction of the lens varies from a value of 1.41 to 1.43 in the core, to 1.38 in the surface, more monotonically in the young human lens, and with a plateau shape and an abrupt change toward the periphery in the old lens.^{20}^{–}^{22} There are few reports in the literature of the GRIN distribution in the intact lens, from non-destructive methods,^{23},^{24} and only three in lenses in vivo^{21},^{22},^{25} using magnetic resonance imaging. Although magnetic resonance imaging offers the opportunity to obtain non-invasive measurements of the refractive index, it is a complex and costly technique that cannot be readily applied on a large scale and suffers from limited resolution. It must also rely on an indirect calibration which introduces some uncertainty.

The standard processing of OCT anterior segment images involves simple division of the optical path length signals by the tissue group index of refraction. A further sophistication of the method involves ray tracing through the tissue and correction for the distortion produced by refraction, considering the deflection of the rays after the optical surfaces but assuming a constant refractive index in the cornea, in the lens, and in the aqueous humor. The new quantitative tools developed for anterior segment imaging of the eye aim at producing fully quantitative images of the optical structures in vivo.^{26} Fan and optical distortion correction algorithms have made it possible to estimate the true shapes of all the ocular surfaces from the anterior segment in vivo, although the posterior surface of the crystalline lens is affected by the GRIN distribution in the lens.

We have previously reported on the impact of the GRIN on the visualization and quantification of the posterior shape of crystalline lenses in vitro from OCT images.^{19} Comparisons with the actual lens posterior shape were possible, because the crystalline lenses were also imaged with the posterior surface up. Those comparisons suggested that, although the presence of GRIN seemed to have a minor influence on the estimated radius of curvature, the peripheral shape of the lens is misestimated.

In this study, we have further developed the optical distortion correction algorithms by incorporating a model of the GRIN distribution in the crystalline lens. With an appropriate model of the GRIN distribution, the new algorithm would enable the in vivo estimation of the shape of the entire crystalline lens and to assess its changes during aging or accommodation, providing a better understanding of these processes. Interestingly, optical imaging techniques used to date to image the crystalline lens in vivo and to obtain quantitative parameters of its shape (Purkinje and Scheimpflug imaging) did not take into account the presence of GRIN in the lens. To our knowledge, these methods assumed a constant equivalent index of refraction even when provided with distortion correction algorithms.^{27}

#### METHODS

We propose a method for correction of optical distortion through the crystalline lens, applicable to estimating the undistorted posterior shape of the crystalline lens and lens thickness in vivo. The method is applied to a set of nine human lenses in vitro, as this allowed us a direct comparison of the reconstructed posterior lens surface and the actual posterior lens surface (also imaged with the posterior lens up) and the lens thickness. Images were obtained two-dimensionally with a custom-built time-domain OCT system, and the algorithm was based on thru-GRIN ray tracing using Sharma's computation scheme,^{28} which is basically a numerical solution of the ray equation based on the Runge-Kutta method. The Sharma algorithm allows the calculation of the optical path through the GRIN medium, step by step, by retrieving the coordinates of intermediate points through which the ray passes, together with directional cosines of the ray at those points. Fig. 1 depicts a schematic diagram of the whole GRIN optical distortion procedure.

##### OCT Images

Cross-sectional images (B-scans) of human lenses in vitro were captured using a custom-built time-domain OCT system provided with a telecentric beam delivery system which produces a flat scan field. The system is described in detail in previous works.^{29} In brief, 20 A-lines/s are acquired, with an axial resolution of 12 μm and an axial length of 10 mm in air. The light source comes from a superluminescent diode with a 825 nm central wavelength and a 25 nm bandwidth.

Full two-dimensional images of the lenses, immersed in a cuvette with preservation medium, were obtained with a resolution of 500 lines/B-scan and a lateral range of 20 mm. Images of the lenses were taken in two positions: (1) anterior surface facing the OCT beam (anterior-up position) and (2) posterior surface facing the OCT beam (posterior-up position), with a special care taken to ensure that the OCT tomograms were captured on the same meridians for the two lens orientations (although some variability may arise from small differences in alignment).

Denoising and segmentation programs were used to extract the edges of the lens, for the two orientations, as well as the edges of the cuvette. The anterior lens surface shape was calculated by dividing the height of the points by the refractive index of the preservation medium at the measurement wavelength (n = 1.345), because rays are parallel to the optical axis. Lens thickness and the posterior lens shape were also obtained, although these data were only used as reference for comparison with data obtained only from the lens imaged in one orientation. Lens thickness (and average refractive index on axis) was obtained from the distortion produced on the cuvette holding the in vitro crystalline lens and preservation medium following a procedure previously described by Uhlhorn et al.^{29} The actual posterior lens surface was obtained from the lens imaged in the posterior-up position following identical procedures described for the anterior lens surface.

##### Human Lens Samples

The human eyes were obtained from the Florida Lions Eye Bank and used in compliance with the guidelines of the Declaration of Helsinki for research involving the use of human tissue. The mean age of the donor was 45 ± 20 years (ranged from 6 to 72 years), and eyes were received within 48 h postmortem, with measurements performed within an hour of lens extraction from the globe. The lenses used in this study are part of a set of lenses previously used in a previous study^{17} for the purpose of GRIN reconstruction, where the lens handling protocols are described in detail.

##### GRIN Distortion Correction Algorithm

In the “natural” orientation (anterior-up position), the posterior surface of the lens is distorted by the deflection of the rays by the anterior surface and by the presence of a GRIN distribution in the lens. To retrieve the information on its shape, we developed an iterative algorithm based on Sharma's ray tracing algorithm in GRIN medium,^{28} where the parameters responsible for posterior lens shape, and in consequence, for the GRIN distribution, are changed between iterations. The algorithm consists of the following steps:

1. An anatomically plausible model of the GRIN distribution is assumed. The GRIN distribution is described by parameters that are independent of the lens shape (central and peripheral refractive indices, average refractive index, and GRIN profile factor)^{17} and dependent of the lens shape (radius R_{a} and conic constant k_{a} of the lens anterior surface, radius R_{p} and conic constant k_{p} of the posterior lens surface, and thickness t). A more extensive description of the GRIN model used is provided in the next subsection.

2. The actual shape of the anterior surface (radius R_{a} and conic constant k_{a}) can be estimated directly from the measurement, as described above.

3. The optical path inside the crystalline lens is estimated directly from the OCT measurement.

4. An initial posterior lens surface defined by R_{p}, k_{p}, and t is assumed. For the purpose of this calculation, the initial posterior lens surface was computed applying optical distortion correction algorithms assuming a constant refractive index.

5. Rays are traced through the GRIN model using Sharma's algorithm^{28} without calculating the intercept between the ray and the posterior surface. If the surface is posterior to the measured optical path distance (OPD), then the ray tracing is stopped before the ray hits the surface. Otherwise, the model is extended so that the estimated OPD coincides with that measured with OCT at every location. As the Sharma algorithm is based on the Runge-Kutta method, it provides a discrete value which does not provide an exact match of the optical path length within the lens. However, we verified that this effect was negligible for sufficiently small discrete step sizes. A total of 400 rays, within a 4-mm pupil, were used in the ray tracing. Unless otherwise noted, the Sharma step was set to 1 μm.

6. The estimated locations in step (5) are fitted by a conic (new R_{p}, new k_{p} and new t) and the values of R_{p}, k_{p} and t parameters are substituted by the values of the new ones. The algorithm returns to step (4).

7. The iterations are stopped when the surface resulting from two consecutive iterations are comparable [i.e., <0.1 μm difference in terms of root mean square (RMS) metric].

##### GRIN Model

For the purposes of this study, we used a GRIN model introduced by Manns et al.,^{30} where the GRIN is described by means of power coefficient p from the nucleus (having the refractive index of n_{n}) to the surface (having the refractive index of n_{s}):

where ρ_{s} is the distance from the center of the lens to the surface at angle θ and p is the power coefficient of the GRIN distribution.

The shape-independent parameters of this model were obtained from a previous study, where the same OCT data were used to estimate the GRIN distribution from the shape of the distorted posterior surface. In our calculations, we used both exact parameters of GRIN distribution estimated for each eye (in a previous study^{17}) and age-related fits to the experimental data. The surface and refractive indices at the nucleus (n_{n}) and surface (n_{s}) were obtained by linear fits to the data. A statistical analysis revealed that this was approach was not statistically significantly different than the mean (n_{n} = 1.425 and n_{s} = 1.381). The power coefficient (p) was fitted exponentially as a function of age (see Fig. 2 for details).

##### Data Analysis

The accuracy of the lens posterior surface shape and lens thickness obtained from the application of the developed refraction/GRIN distortion correction algorithm on OCT images was compared to those obtained from two other approaches for posterior lens shape and lens thickness estimations. (a) Division of heights of the points of distorted surface by a homogeneous index. This approach does not take into account refraction by the anterior lens surface and the presence of a GRIN distribution in the lens. This approach is followed widely in OCT imaging. (b) Application of optical (refraction) correction algorithms as proposed by Ortiz et al.,^{15} considering the refraction at the anterior lens surface but assuming a homogeneous refractive index. The algorithm works by calculating the refraction at every point of the anterior surface and estimating the locations where the estimated OPD coincides with the measured one. For (a) and (b), the homogenous average refractive indices were obtained from the study of Uhlhorn et al.,^{29} where the group (not phase) refractive index was obtained at 825 nm and an age-dependent expression is provided.

The accuracy of the posterior shape and lens thickness correction methods [the standard methods (a) and (b) above] as well as the new method—using individual or fitted parameters in the GRIN model—was given as the differences with respect to the actual posterior lens shape and lens thickness obtained from the crystalline lens.

A one-way analysis of variance was used to test for significance of differences between the correction results.

#### RESULTS

##### Accuracy in the Reconstruction of the Posterior Lens Shape

Fig. 3 shows the shifts in the reconstructed lens surface radii of curvature (upper panels) and conic constant (lower panels) for individual subjects (a, c) and average across subjects (b, d) for the different reconstruction methods. Data are given relative to the actual parameters on the posterior lens shape (obtained from posterior-up measurements). Although there are great differences across subjects, the largest shift in radius of curvature, on average, occurs for the simple division method, while the refraction distortion correction (with a homogeneous index) and GRIN distortion correction (for both fitted and measured data^{17}) provide similar estimates of radius and conic constant in the posterior surface. The conic constant tends to be similarly retrieved with all methods.

Although an independent analysis of radius of curvature and conic constant did not reveal significant differences across reconstruction methods [F(3,32) = 0.056, p = 0.982 for radius of curvature and F(3,32) = 0.998, p = 0.960 for conic constant], a comparison of the overall surface shape showed differences across methods. Fig. 4 shows an example of difference maps (reconstructed-actual) for the four different reconstruction methods for a 41-year-old eye. When differences are expressed in terms of RMS and peak to valley differences, the GRIN distortion correction method (particularly with the actual GRIN parameters) showed significantly higher accuracy [F(3,32) = 3.260, p = 0.034 and F(3,32) = 3.212, p = 0.036 for RMS and peak to valley analysis, respectively], with the lowest accuracy found for the simple division by refractive index method (Fig. 5).

Figure 4 Image Tools |
Figure 5 Image Tools |

##### Accuracy in the Estimates of Lens Thickness

In particular, the GRIN distortion correction method produced significantly better estimates of the crystalline lens thickness [F(3,32) = 3.983, p = 0.032] compared with the other methods (Fig. 6). The simple division and refraction (with homogeneous refractive index) provided, as expected, identical thickness estimates.

##### Influence of Sharma Step Size

Decreasing the step of the Sharma ray tracing algorithm increased the accuracy of the reconstruction (at the expense of increasing computational time). Fig. 7 shows the average difference between the reconstructed and actual posterior surface (in terms of RMS or peak to valley) as a function of the ray tracing step for both implementations of the GRIN distortion correction algorithm. Although in all cases the accuracy is very high (<30 μm), reducing the iteration step increased accuracy.

##### Convergence

The reconstruction algorithm is characterized by a good convergence. For all the processed lenses, no more than 10 iterations were needed to reach the final parameters in the GRIN distortion correction procedure. In isolated cases (6-, 31-, and 48-year-old lenses), the reconstruction algorithm provided two slightly different local minima, resulting in two slightly different sets of values of radius of curvature and conic constant but the same value of lens thickness. Although for these minima the differences of RMS error were small (up to 1‰), the results described above referred to data from the solution which provides the smaller RMS error between the actual and reconstructed shapes.

#### DISCUSSION

In this study, we proposed an iterative method of optical distortion correction in OCT images of the crystalline lens incorporating the GRIN distribution inside the crystalline lens medium. A comparison with other existing methods of posterior lens shape reconstruction shows slight improvements in the shape reconstruction for averaged input parameters of the GRIN model distribution and a significant improvement when optimal parameters for the GRIN (i.e., those corresponding to the same lens) were used. Interestingly, the method has proved to reconstruct with high accuracy lens thickness, with knowledge of the anterior surface shape and distorted posterior lens shape only (and a model for GRIN distribution). The accuracy of the reconstruction (6 μm in the optimized GRIN reconstruction method) is similar to the reported accuracy of shadow photography (12 μm), broadly used for biometric measurements of the crystalline lenses in vitro,^{31} as tested on calibrated spheres. Our findings show that with an appropriate GRIN model and refraction distortion correction from preceding surfaces, the GRIN distortion correction algorithms, one could reach similar accuracies in vivo. In fact, such accuracy is similar to the repeatability of in vivo lens thickness measurement by means of commercial OCT (8 μm).^{32} It needs to be noted that according to the results presented in the previous section, the thickness estimation has significant contribution to the RMS errors between the actual and reconstructed surfaces.

No significant differences were found across methods in the retrieved radius of curvature and conic constant, although the surface shape is generally best retrieved with the new reconstruction algorithm. This may be in part due to inaccuracies of the surface fitting. Urs et al.^{33} studied different methods for fitting the contours of isolated lens images and found RMS errors in the fits, ranging from 11 to 70 μm, when using 10th order polynomial one lens curve fitting in the posterior lens surface. These fitting errors are of the order of magnitude of the RMS errors of our reconstruction. The use of conic surfaces, with less fitting parameters, is likely prone to higher inaccuracies. In fact, we have shown that comparisons of surfaces using a separate analysis of radius and asphericity may estimate incorrectly the statistical significance in the differences between surfaces, as various combinations of radius and asphericity may describe with similar accuracy the same noisy surface.^{34} For example, the ranges of (correlated) radius and asphericity which described similarly Scheimpflug posterior corneal elevation data were close to 0.2 mm and 0.6, respectively, only slightly lower than the accuracy found for those parameters in this study.

Our reconstruction algorithm is suited for in vivo OCT images of the crystalline lens. The shape of the posterior lens surface measured posterior-up and the distorted shape of the cuvette, which can only be measured in vitro, have been used in this study only for comparison purposes. They are not required as inputs in the proposed correction algorithm. However, a limitation of the study is the general lack of GRIN distribution parametric data. Furthermore, as the OCT technique uses a low coherent light source, which is characterized by broadband wavelengths in the near infrared, group refractive indices in this bandwidth are needed, taking into account the chromatic dispersion of the lens medium. In this study, we used GRIN data that came from our previous work,^{17} which included the set of lenses evaluated in this study. As a result, it is not surprising that the best reconstruction was achieved with the optimal shape parameters of the GRIN in the individual lenses. As more data for GRIN distribution in larger sample become available, it is likely that the GRIN parameters obtained from fitting represent more robustly the population data. It should be noted that the GRIN distribution of isolated lenses will represent more closely the GRIN profile in a maximally accommodated lens, which should be considered when extrapolating GRIN models to perform reconstructions from OCT in vivo measurements.^{22} GRIN distribution estimates in lenses in vitro under simulated accommodation (i.e., with an artificial stretching system) may allow a more direct application of GRIN models to measurements in vivo.^{35}

In this study, we have described a possible implementation of the algorithm, but its core is quite flexible. Basically, instead of Sharma's algorithm, any of the numerical ray tracing and optical path estimation procedures^{36}^{–}^{38} could be adapted to different GRIN models. The only condition for the GRIN model is that it is anatomically plausible, with the isoindicial surfaces functionally related to the external shape of the lens. Potentially, weights could be introduced in the optimization of the GRIN shape parameters, in step (6) of the algorithm, while still guaranteeing the convergence of the algorithm.

We have described an implementation of the algorithm in two dimension, but it could be easily extended three dimensionally. In a previous work, we demonstrated computationally a significantly higher accuracy in the reconstruction of the posterior surface when optical distortion correction algorithms were applied in three dimension.^{15} For example, in a computer eye model with homogeneous index of refraction (and simulated conic surfaces), there was no significant difference in the reconstruction of the posterior lens asphericity by applying a simple division by the refractive index of refraction or two-dimensional refraction distortion correction (with discrepancies of about 50% from the nominal value). However, three-dimensional refraction distortion correction allowed retrieval of the asphericity within 0.3% error.

The algorithm for optical distortion correction through GRIN proposed in the current study will therefore show full potential on three-dimensional images of the crystalline lens in vivo,^{26} particularly as GRIN models of the crystalline lens based on larger populations than current data become available.

Damian Siedlecki

Institute of Physics

Wroclaw University of Technology

Wybrzeze Wyspianskiego 27

50370 Wroclaw, Poland

e-mail: damian.siedlecki@pwr.wroc.pl

##### ACKNOWLEDGMENTS

We thank Bianca Maceo and Raksha Urs for assistance with data processing.

The study was supported in part by Spanish Ministry of Science and Innovation Grant FIS2008-02065 and FIS2011-25637 (Marcos), EUROHORCs-ESF EURYI-05-102-ES (Marcos), ERC-2011-AdG-294099 (Marcos), CSIC JAE-Program (Siedlecki, de Castro, Gambra), National Eye Institute Grants 2R01EY14225, 5F31EY15395 (NRSA Individual Predoctoral Fellowship, Borja), P30EY14801 (Center Grant); the Florida Lions Eye Bank; an unrestricted grant from Research to Prevent Blindness, the Vision Cooperative Research Centre, Sydney, New South Wales, Australia, supported by Polish Ministry of Science and Higher Education, the Australian Federal Government through the Cooperative Research Centres Programme and the Henri and Flore Lesieur Foundation (Parel).