1. Introduction

The crystalline lens is a biconic lens with a gradient refractive index (GRIN) that is known to peak near the lens center and decrease toward the periphery. The lens contributes about 30% of the optical power of the eye [1] and, to some extent, compensates for some of the aberrations introduced by the cornea [2]. The gradient index of refraction of the lens is known to change with age [3–7] and accommodation [5,8,9], and is responsible for part of the change in the lens spherical aberration [8,10].

Information on the GRIN of the lens has been derived mostly from measurements on excised lenses using methods that either solve iteratively the ray equation [11–15] or optimize the variables of a model to fit the experimental data [16–18]. To apply these methods, which rely on precise measurements of the external geometry, the lens must be extracted from the ocular globe. This approach naturally limits the applicability in vivo, since the posterior lens surface can only be imaged through the unknown GRIN. The distortion introduced by the GRIN prevents accurate quantification of the posterior lens surface shape, particularly its asphericity [19,20].

Imaging methods, such as Magnetic Resonance Imaging (MRI) [4] or X-ray Talbot interferometry [21], on the other hand, have shown that it is possible to study the crystalline lens refractive index without information about the external lens geometry. X-ray Talbot interferometry can only be applied in vitro [21,22], but MRI has been used both in vitro [4,23] and in vivo [5,24,25]. However, MRI relies on sensitive calibrations to relate the image intensity and the local refractive index using the relationship between the proton transverse relaxation time and the concentration of proteins [26]. In addition, the complexity and cost of MRI imaging systems limit its applicability for routine use in large scale studies.

In the past, our group has presented experimental validation of an optical method applicable in vitro that searches for the best values for the parameters of a 4-variable model that fit the optical path difference between anterior and posterior lens surface collected from OCT imaging [27]. A pivotal step of the method is to image the crystalline lens twice: first with the anterior lens surface facing the OCT beam (anterior-up image), and then with the posterior lens surface facing the OCT beam (posterior-up image). With this method, both lens surfaces are imaged without distortion from the GRIN media, and the external geometry of the lens can be accurately calculated. The distorted posterior lens surface from the anterior-up image is used to reconstruct the GRIN by finding the variables of the GRIN that best fit the distortion by means of a global search algorithm [17,27]. We have applied this method to study changes in the refractive index of human crystalline lenses with age [6,7] and with accommodation in animal models [8], and have quantified its influence on spherical aberration and astigmatism [7,10,28] in both relaxed and accommodated states [8,9].

In the present study, we demonstrate a method to accurately reconstruct the gradient refractive index distribution and posterior lens surface using off-axis OCT images [29] of the lens acquired in the anterior-up orientation, without the need for posterior-up images. This opens up the possibility of reconstructing the GRIN in vivo using OCT. If only on-axis OCT images are acquired, then the shape of the posterior lens surface is coupled to the GRIN distribution and reconstruction is therefore not possible. The acquisition of images at multiple rotations of the OCT delivery system (i.e. off-axis) enables us to uncouple the posterior surface shape and the GRIN.

2. Methods

2.1 Rotational optical coherence tomography system

The experimental system has been described in detail elsewhere [29]. Essentially, we modified the delivery probe of a commercial Spectral Domain Optical Coherence Tomography system (Envisu R4400, Leica-Microsystems GmbH, Wetzlar, Germany) and mounted it on a rotation stage to image the isolated crystalline lens on and off axis (Fig. 1). The rotation stage was motorized to achieve smooth movement of the optical and mechanical systems. The central wavelength of the system was 880 nm and the nominal axial resolution was 8 µm. The system was programed to acquire 2D cross-sectional images (B-scans) at each rotation of the delivery probe (1000 A-lines x 2048 pixels/A-line, 10.0 x 7.4 µm pixel size in air). The scan width was 10 mm at normal incidence. The position of the crystalline lens was adjusted so that the B-scan imaged passed through the lens center by displacing the lens perpendicularly to the image plane until both the anterior and posterior surface reflex were visible. Also, shifting the vertical position of the cuvette allowed to avoid displacements in the image of the lens at different rotations and ensured that the lens was within the image range at all angles. Images were then collected by automatically rotating the OCT beam between −45 and +45 degrees in steps of 5 degrees. There was no need to correct for the fan distortion arising from the scanner architecture [30] because this was negligible in the central area where the rays optical path are measured [29].

 

Fig. 1 Picture of the OCT probe mounted in a rotation stage (A) and (B), schema of the experimental setup (C) and (D), and OCT images for one of the lenses (E) and (F) on axis and with a rotation of 45 degrees, respectively.

Download Full Size | PDF

2.2 Crystalline lenses

We used 9 in vitro crystalline lenses from 7 cynomolgus monkeys (Macaca fascicularis, post-mortem time (i.e. time from eye enucleation to measurement): 20 ± 15 hours, age: 7.4 ± 3.2 years old). All experiments adhered to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Visual Research. The eyes were obtained from the Division of Veterinary Resources at the University of Miami as a part of a tissue-sharing protocol and were used in accordance with Institutional Animal Care and Use Guidelines. Enucleation of the eyes was performed immediately after euthanasia and the samples were wrapped in gauze and stored in a closed container. No animals were euthanized for the sole purpose of this study.

After their arrival at the laboratory, the eyes were either immediately prepared by an ophthalmic surgeon for the measurements, or refrigerated at 4°C. Preparation involved extracting the crystalline lens from the eye and placing it on a rubber ring in a chamber filled with preservation medium (Dulbecco’s Modified Eagle’s Medium, DMEM) [31]. The lenses were imaged in the anterior-up position on- and off-axis with the OCT rotating the imaging beam between −45 and +45 in steps of 5 degrees, and then flipped over and imaged with the posterior lens surface facing the OCT beam on-axis. This allowed us to compare the undistorted posterior lens surface shape to the reconstructed posterior surface in order to assess the error of the method.

2.3 Image analysis and registration

The anterior and posterior surfaces of the crystalline lens were segmented in all the images (on-axis and off-axis) using semiautomatic algorithms written in MatLab (Mathworks, Natick, MA). In the off-axis images, refraction on the interface air-preservation medium was taken into account (Fig. 2). The angle of the beam inside the preservation medium was calculated with Snell’s law, $\sin (\text{α'})={\text{n}}^{\text{-1}}\sin (\text{α}),$ where αʹ is the angle of the refracted rays, n is the refractive index of the preservation medium at the OCT wavelength (n = 1.345) and α is the angle of rotation of the OCT probe. The distance between rays, i.e. the horizontal pixel size in the OCT images, was enlarged due to refraction by$\text{cos(α')/cos(α),}$since two collimated rays separated a distance d will intersect the air-preservation media interface at points separate h = d / cos(α) and will continue traveling inside the preservation media with an angle α’; hence separated inside the medium a distance $\text{d'}=\text{h}\cos (\text{α'})=\text{d}\cos (\text{α'})/\cos (\text{α}).$These equations were used to calculate the angle of incidence and horizontal position of the rays entering the lens through the central 4 mm.

 

Fig. 2 (A) Schematic of the OCT rays refracted at the air-DMEM interface. The rays impact the lens with an angle different than that of the OCT probe, in the example of the figure, $\alpha \text{'}$ = 21.82°. (B) OCT image of the lens with the probe tilted by 30 degrees. (C) Processing steps to calculate the anterior surface apex, which was assumed to be the lens axis. The detected anterior surface was corrected for DMEM distortion, rotated by an angle, $-\alpha \text{'}$, and fitted to a conic to detect the apex of the anterior surface (shown in blue). (D) OCT image of the lens with the detected anterior surface apex superimposed.

Download Full Size | PDF

When the OCT probe was aligned with the lens axis, the apex of the lens anterior surface, i.e. the lens center, was located by fitting the data to a conic surface. This point was calculated from the images obtained when the OCT probe was tilted with the following three steps: 1) correction of the anterior surface of the lens for distortions due to the preservation medium, by dividing the sag of the surface by the DMEM refractive index, 2) rotation of the surface data points by an angle equal to that of the refracted OCT beam, and 3) data fit by a conic surface (Fig. 2). To ensure the accuracy of this procedure, we verified that the radius of curvature and conic of the anterior surface measured from the images with different rotations of the OCT probe were similar. Once the apex was located, the input data to the ray tracing algorithm, i.e. the points of incidence of each ray and their corresponding tilts before the lens surface, were calculated using the equations above. The experimental optical path difference between the anterior and the posterior surface was then calculated from the images.

2.4 Ray tracing and reconstruction algorithm

A custom ray tracing algorithm [17] was used to calculate the optical path of rays passing through the lens at varying angles. The crystalline lens GRIN model [17] is centered on the lens optical axis at a distance from the anterior surface of 0.41 times the thickness of the lens [32]. A power-law equation is used to describe the change in refractive index between the center and the surface in polar coordinates:

Both surfaces of the crystalline lens were fitted to conics using the equation

The reconstruction problem was posed as a minimization problem. We programmed a merit function to calculate the sum of the root mean squares difference between the measured optical path and the optical path calculated with the ray tracing algorithm for each ray and each orientation. The 4 parameters of the GRIN distribution (surface and nucleus refractive index and exponent in the axial and meridional directions) and one (radius) or two (radius and conic) parameters of the posterior surface geometry, were set as variables. The optimization algorithm explored different combinations of the variables to minimize the merit function. To ensure that the problem did not converge in a local minimum, genetic algorithms were used [17]. Due to the stochastic nature of the genetic algorithms, the search was repeated three times, and the convergence of the algorithm was assessed by confirming that the solution of three consecutive repetitions was within three digits for the refractive index and within 2 digits for the radius and conic constant.

2.5 Reconstruction

Three optimization problems were solved separately with increasing number of variables: (a) reconstruction of GRIN using the external lens geometry (4-variable problem), (b) reconstruction of both GRIN and the radius of curvature of the posterior surface (5-variable problem), and (c) reconstruction of the GRIN and the radius and conic constant of the posterior surface (6-variable problem). Problem (a) can be solved experimentally on crystalline lenses when both anterior-up and posterior-up images are acquired to enable direct measurement of both the anterior and posterior lens surfaces. It is equivalent to the in vitro reconstructions presented before [6–8,17,28] except that off axis data were also used for the reconstruction in the present study, but not lens power which was used as an input constrain before. Here only the optical path obtained from the on- and off-axis OCT images is used to reconstruct the GRIN distribution. Problems (b) and (c) do not require acquiring images in more than one position (only anterior-up images), and are used to show the potential of the technique for future in vivo applications.

Simulations were performed to estimate the accuracy of the reconstruction. Optical path difference for rays at different orientations were calculated with a crystalline lens surface shape and a given GRIN distribution. The data was then used as input to reconstruct a) the 4 GRIN parameters b) the 4 GRIN parameters and one (radius) posterior surface parameter, c) the 4 GRIN parameters and two (radius and conic) parameters of the posterior surface. The reconstruction problem was solved after adding Gaussian noise of different variances to the optical path difference to more closely simulate real-life conditions. The simulations were repeated with different Gaussian noise distributions with the same magnitude and the average and standard deviation of the difference with the nominal parameters was calculated.

The reconstructed posterior lens surface shape was compared with the experimental measurements of the same surface derived from posterior-up images. To evaluate the effect of adding one (posterior surface radius) or two (posterior surface radius and conic constant) additional variables to the optimization in the GRIN reconstruction, the reconstructed GRIN distributions were compared.

3. Results

3.1 Computer simulations

The simulations show that off-axis data are essential for reconstructing the posterior lens surface and the GRIN distribution using OCT imaging, when the crystalline lens is imaged from only one orientation. Figure 3 shows the error in the reconstruction with 10 µm of Gaussian noise added to the calculated optical paths for different values of the maximum OCT probe angle. The same number of rays, 240, was used in all reconstructions. If the external geometry was known, the error in the reconstruction of the four GRIN parameters using probe angles of up to 45 degrees was, on average, 0.007, 0.004, 1.1 and 0.9 for the 4 parameters defining the GRIN distribution (surface and nucleus refractive index, and axial and meridional exponent, respectively). When the radius of curvature of the posterior lens surface was reconstructed, the reconstruction error was on average, 0.008, 0.004, 1.5 and 1.0 for the 4 GRIN parameters and 0.058 mm for the posterior surface radius. When both the radius and the conic constant were set as variables, the reconstruction error was, on average, 0.009, 0.004, 1.7 and 1.1 for the 4 GRIN parameters, and 0.089 mm and 0.3 for the radius and conic respectively. The described methodology provided a higher accuracy than when only on-axis information was used. The reconstruction errors in the latter case were: 0.012, 0.005, 1.6, and 1.5 for the nominal GRIN parameters for the 4-variable problem; 0.020, 0.006, 3.4, and 2.6 for the GRIN parameters, and 0.148 mm for the radius of curvature of the posterior surface for the 5-variable problem; and 0.014, 0.005, 2.3, and 2.7 for the GRIN parameters, and 0.142 and 0.27 for the radius of curvature and conic constant of the posterior surface for the 6-variable problem.

 

Fig. 3 Results of the reconstruction error from computer simulations (assuming 10 µm Gaussian noise added to the measured optical path). (A) Reconstruction error for the 4-variable problem: 4 GRIN parameters, using anterior and posterior-up images; (B), (D) Reconstruction error for the 5-variable problem: four GRIN parameters (B), and one posterior surface parameter (lens radius of curvature) (D), using only information from anterior-up images; (C), (E) Reconstruction error for the 6-variable problem: four GRIN parameters (C), and two surface parameters (lens radius of curvature and conic) (E), using only anterior-up images. Error bars represent the standard deviation of 50 repetitions of the reconstruction algorithm with different errors. As the number of orientations used for the reconstruction increased, the accuracy of the reconstruction increased for both the GRIN parameters and the posterior surface shape.

Download Full Size | PDF

3.2 Experimental results

The average difference between the radii of curvature of the anterior surface measured from images at different rotations of the delivery probe for all the lenses was 4% (maximum difference 7.5%). Figures 4 and 5 show the results of the three optimization problems for all experimental lenses. The surface refractive index was found to be very similar across lenses (average ± STD: 1.365 ± 0.004), whereas the nucleus refractive index and the axial and meridional exponents showed higher inter-subject variability (1.451 ± 0.008, 6.2 ± 2.5 and 9.2 ± 1.8 respectively). Including posterior surface parameters in the reconstruction did not significantly change the reconstructed GRIN parameters. When compared to the GRIN distribution reconstructed using the posterior-up image (lens shape), the differences between refractive index values within the central 4 mm in terms of root mean square were always below 0.015, with an average difference of 0.006 and 0.008 for the 5 and 6-variable problem respectively.

 

Fig. 4 Parameters of the reconstructed GRIN for all 9 crystalline lenses: (A) surface refractive index, (B) nucleus refractive index, (C) axial exponent and (D) meridional exponent. Black bars stand for the 4-variable problem (GRIN variables were optimized), gray bars for the 5-variable problem (GRIN and the posterior surface radius were optimized) and white bars for the 6-variable problem (GRIN and two parameters of the posterior surface were optimized).

Download Full Size | PDF

 

Fig. 5 Lens GRIN map reconstructions for the 9 crystalline lenses. Upper row: 4-variable problem (GRIN reconstruction); Middle row: 5-variable problem (GRIN and posterior surface radius of curvature reconstruction); Lower row: 6-variable problem (GRIN and posterior surface radius of curvature and conic constant reconstruction). Differences in the posterior surface shape between the three reconstruction problems are not noticeable. Scale bar represents 1 mm.

Download Full Size | PDF

The reconstructed geometry of the posterior surface of the lens was compared to the corresponding undistorted lens surface shape obtained experimentally (from direct images of the crystalline lens with posterior surface up). Figure 6 shows reconstructed and measured radius of curvature (A), and corresponding Bland-Altman plot (B), when the GRIN parameters and the radius of the posterior surface were set as variables (5-variable problem). The difference between the measured and the reconstructed radius of curvature was below 10% in all cases and was on average 0.13 mm. The root mean square (RMS) of the difference between measured and reconstructed posterior surface shape was 32 ± 27 µm.

 

Fig. 6 Reconstructed posterior lens surface shape in the 5-variable problem. (A) Reconstructed radius of curvature compared with the radius of curvature measured experimentally. (B) Bland-Altman plot.

Download Full Size | PDF

Figure 7 shows the reconstructed and measured posterior surface radius of curvature (A) and conic constant (C), and corresponding Bland-Altman plots (B, D), when the GRIN parameters, the radius of the posterior surface and the conic constant were set as variables (6-variable problem). The difference between the measured and the reconstructed radius of curvature was below 15% for all the lenses and was 0.20 mm, on average, and the average conic constant difference was 0.47. The RMS between the measured and the reconstructed surface was, on average, 38 ± 22 µm.

 

Fig. 7 Reconstructed posterior lens surface shape (radius of curvature and conic constant) in the 6-variable problem. Comparison between the measured and the reconstructed posterior surface radius of curvature (A) and corresponding Bland-Altman plot (B). Comparison between the measured and reconstructed conic constant (C) and corresponding Bland-Altman plot (D).

Download Full Size | PDF

4. Discussion

In the present study, we have demonstrated a method to accurately reconstruct the gradient refractive index and the posterior surface shape exclusively using on and off-axis OCT images of the crystalline lens with its anterior surface facing the OCT beam. The GRIN and the posterior surface were successfully reconstructed in 9 in vitro cynomolgus monkey crystalline lenses with OCT images obtained by tilting the OCT probe from −45 to 45 degrees.

The precision in the reconstruction of the 4 GRIN parameters (surface and nucleus refractive index and axial and meridional exponent) when introducing posterior surface shape as additional variables of the optimization problem, (0.008, 0.004, 1.5, 1.0, and 0.009, 0.004, 1.7, 1.1, for the 5 and 6-variable problem respectively) is comparable to the precision when only the GRIN is reconstructed (0.007, 0.004, 1.1 and 0.9). The fact that the error does not increase for the nucleus refractive index together with the lower error in the reconstructions (0.004 in the nucleus vs 0.008 for the surface, on average) is possibly due to the higher influence of this parameter in the optical path of the rays. The results of the simulations show that only using the information from the off-axis images, the error in the reconstruction of the posterior radius of curvature could be lowered to the third decimal place.

The simulations presented in this manuscript assume that the error in the experimental data can be modeled with a 10 µm Gaussian noise. However, a decrease in the signal to noise ratio in some parts of the image of the surface of the crystalline lens acquired at high rotation angles of the OCT probe, may lead to a slightly higher error in the optical path data obtained from off-axis images. In addition, while care was taken to image the central meridian of the lens, slight errors in the positioning and tilts of the lens may also affect the quality of the GRIN reconstruction.

The convergence of the optimization algorithm was faster when more off-axis information was used. To decrease the average deviation of the variables in the initial population by 66%, more than 40 generations were needed when only on-axis ray tracing was used, while less than 15 generations were sufficient when off-axis information was used as additional input to the genetic algorithm. This result suggests that the off-axis information smooths the solution space, which in turn decreases the number of required generations of the algorithm to converge. In this study, the computational time was about 1 hour for each reconstruction.

The surface, nucleus refractive index and exponents found in the crystalline lenses studied using the four-variable problem (average ± standard deviation 1.365 ± 0.004, 1.451 ± 0.008, 6.2 ± 2.5 and 9.2 ± 1.8 respectively) are in agreement with the values we found in a previous study [8] using a different set of cynomolgus monkey lenses (1.375 ± 0.003 and 1.429 ± 0.003 for surface and nucleus refractive indices and exponents ranging from 2.1 to 9.1) and a technique that used images of the crystalline lens in the two orientations. The higher value of the nucleus index in the current data set is consistent with the fact that we find a slightly lower value of the surface index. The contribution of the lens surfaces to total lens power decreases when the lens surface refractive index decreases. The larger difference between nucleus and surface index compensates for this loss of surface power. Some of the differences between the two data sets could reflect age-related changes in the lens shape and refractive index as well as differences in methodology. On average, the donors in the present study were older than in the previous study (7.4 years versus 5.7 years). In addition, the previous study relied only on on-axis images and used the lens power as one of the input parameters of the reconstruction algorithm. The present study uses both on- and off-axis images but does not take into account lens power. The power of the lenses was calculated with the reconstructed GRIN distributions (average power 52.7 D) and the values are within the range of previous reports on the power of the isolated cynomolgus monkey crystalline lens (from 35 to 60 D) [33]. The change in the power of the crystalline lenses when shifting axially the peak of refractive index between 0.36 and 0.46 times the central thickness was always below 4%.

The use of OCT to reconstruct the crystalline lens GRIN was first reported by Verma et al. [16] on isolated Zebrafish lenses, which have spherical symmetry. In their experiment, the complete lens external geometry was measured using a single OCT image, and the only variables in their optimization problem were those of the GRIN model. The monkey or human crystalline lens is not spherical and even in the case of a homogeneous refractive index lens, the optical path traveled by the rays from the anterior to the posterior surface depends both on the refractive index and the anterior surface shape. Our results show that this dependency can be removed if data from multiple orientation images are used.

A rotational OCT has also been used to characterize GRIN manufactured lenses capturing only one A-scan at each rotation position of the OCT probe [34]. This alternative scanning architecture allows for minimizing the optical distortions as well as increasing the signal to noise ratio, which can be important in some imaging modalities [35].

An added difficulty to studying the crystalline lens in vivo is that the lens is viewed through a cornea with high optical power. Optical distortion correction algorithms are applied on Scheimpflug [36–38] or on Optical Coherence Tomography [39–41] images of the anterior segment acquired in vivo to compensate for refraction by the cornea, and subsequently by the anterior surface of the lens. However, these algorithms typically assume a homogeneous index of refraction for the crystalline lens to calculate the shape of the posterior surface of the lens, since the GRIN is unknown. In a previous computational study using in vitro lens images [19], we showed that approximating the GRIN with a homogeneous refractive index produced an error in the estimation of the radius of curvature within the measurement reproducibility. However, to accurately reconstruct the conic constant of the posterior surface, a correction algorithm that takes into account the actual gradient index was needed [20]. In addition, even if the posterior lens radius of curvature can be estimated in principle by using a homogeneous refractive index for distortion correction, there is some uncertainty because the exact value of the homogeneous refractive index is unknown.

The methodology presented here can be adapted to reconstruct the GRIN in vivo, as off-axis OCT imaging of the anterior segment of the eye is feasible on human subjects. In intact eyes, it would be necessary to take into account refraction by the anterior and posterior corneal surfaces. Previous studies have shown that such corrections are feasible [39–41]. Accurate alignment procedures and high SNR would be needed as input. Also, the method presented here assumes that the thickness of the lens is measured. Although this parameter could be introduced as an additional variable in the optimization, we found that the excessive number of variables limited the adequate convergence of the algorithm. Some methods have been envisioned to measure the average axial group refractive index in vivo [42] and hence could be used to determine the thickness of the lens using OCT, although their feasibility in a wide population needs to be proved. In addition, knowledge of the eye’s aberrations, both on-axis and off-axis can provide additional inputs to help refine the GRIN and posterior surface lens reconstruction. These OCT/Laser Ray Tracing data, currently available in in vitro configurations [29], could also be acquired in humans in vivo, up to 40 deg eccentricities [43–45].