Fully automated laser ray tracing system to measure changes in the crystalline lens GRIN profile

Measuring the lens gradient refractive index (GRIN) accurately and reliably has proven an extremely challenging technical problem. A fully automated laser ray tracing (LRT) system was built to address this issue. The LRT system captures images of multiple laser projections before and after traversing through an ex vivo lens. These LRT images, combined with accurate measurements of the lens geometry, are used to calculate the lens GRIN profile. Mathematically, this is an ill-conditioned problem; hence, it is essential to apply biologically relevant constraints to produce a feasible solution. The lens GRIN measurements were compared with previously published data. Our GRIN retrieval algorithm produces fast and accurate measurements of the lens GRIN profile. Experiments to study the optics of physiologically perturbed lenses are the future direction of this research. © 2017 Optical Society of America OCIS codes: (110.2760) Gradient-index lenses; (140.0140) Lasers and laser optics; (170.0110) Imaging systems. References and links 1. R. P. Hemenger, L. F. Garner, and C. S. Ooi, “Change with age of the refractive index gradient of the human ocular lens,” Invest. Ophthalmol. Vis. Sci. 36(3), 703–707 (1995). 2. E. Vaghefi, A. Kim, and P. J. Donaldson, “Active Maintenance of the Gradient of Refractive Index Is Required to Sustain the Optical Properties of the Lens,” Invest. Ophthalmol. Vis. Sci. 56(12), 7195–7208 (2015). 3. P. J. Donaldson, A. C. Grey, B. Maceo Heilman, J. C. Lim, and E. Vaghefi, “The physiological optics of the lens,” Prog. Retin. Eye Res. 56, e1–e24 (2017). 4. A. de Castro, S. Ortiz, E. Gambra, D. Siedlecki, and S. Marcos, “Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging,” Opt. Express 18(21), 21905–21917 (2010). 5. 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). 6. Y. Verma, K. D. Rao, M. K. Suresh, H. S. Patel, and P. K. Gupta, “Measurement of gradient refractive index profile of crystalline lens of fisheye in vivo using optical coherence tomography,” Appl. Phys. B 87(4), 607–610 (2007). 7. C. E. Jones, D. A. Atchison, R. Meder, and J. M. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. 45(18), 2352–2366 (2005). 8. 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). 9. O. Pomerantzeff, M. Pankratov, G. J. Wang, and P. Dufault, “Wide-angle optical model of the eye,” Am. J. Optom. Physiol. Opt. 61(3), 166–176 (1984). 10. O. Pomerantzeff, H. Fish, J. Govignon, and C. L. Schepens, “Wide-angle Optical Model of the Eye,” Opt. Acta (Lond.) 19, 387–388 (1972). 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. E. Acosta, D. Vazquez, L. Garner, and G. Smith, “Tomographic method for measurement of the gradient refractive index of the crystalline lens. I. The spherical fish lens,” J. Opt. Soc. Am. A 22(3), 424–433 (2005). Vol. 8, No. 11 | 1 Nov 2017 | BIOMEDICAL OPTICS EXPRESS 4947


Introduction
The crystalline lens is a critical component of the eye, and is essential for focusing light onto the retina. The lens compensates for the corneal positive spherical aberration through an inherent negative spherical aberration generated by its gradient of refractive index (GRIN). The existence of the lens GRIN also increases its possible equivalent optical power, enabling flatter lens geometries such as the human lens [1]. A precise knowledge of this GRIN will provide a better understanding of the optical properties of the lens and its contribution to the overall optical quality of the eye. Furthermore, an accurate measurement of the lens GRIN profile, and how it is altered when lens homeostasis is disrupted, will provide insight into the relationship between the optical properties and physiology of the crystalline lens [2], [3]. Such information is essential for understanding age-related changes to the optics of the lens, which lead to refractive errors, presbyopia and cataract.
There are several non-destructive methods that can be used to measure the lens GRIN. These methods include: optical interferometry [4-6], magnetic resonance imaging (MRI) [2], [7,8], using measured ocular and visual parameters [9], [10] and laser ray tracing (LRT) [11][12][13][14][15]. LRT has become the standard for non-destructive GRIN measurements of the in vitro lens. LRT images the amount of deflection of incoming laser beams produced by the lens, when a collection of incident laser beams is shone at it. The earliest LRT algorithms were developed for spherical lenses and required the matching of the surface index to the external refractive index (RI) [16,17]. They also required models which were overly constraining; isoindicial surfaces had to be concentric with and parallel to the surface of the lens [16,17]. More recent methods have overcome these limitations. However, new difficulties were also introduced, demanding substantially more mathematical computation and exhibiting significant difficulty in retrieving a unique solution [18]. Furthermore, a contradiction arose in the choice of mathematical models used. As more complex models had to be used to describe the geometry and GRIN of the lens, the increased complexity amplified the likelihood of falling into erroneous local minimum solutions [12,13].
Genetic algorithms based on optical interferometry have previously been used to overcome this limitation, but are computationally expensive (taking about 45 minutes), especially in high-dimensional problems [4,18]. Moreover, complex models exponentially increase the solution search space size and there is usually no clear convergence criterion for the search algorithm. The 'better' solution is relative to the previous answer. Subsequently, the solution may not always be the desired global optimum [19]. Recent work by Vazquez et al. retrieved the in vitro lens GRIN using a tomographic algorithm [12,13]. The only assumption made is that the lens GRIN profile is rotationally symmetric about its optical axis. The technique allows for a generalized application to non-spherical crystalline lenses. In addition, the method allows for a non-constraining GRIN model (9 coefficients to describe the RI in the lens domain) and does not require RI matching between the surrounding solution and the lens surface.
In this study, we improved upon previous LRT techniques. Specifically, (1) we reduced the computational load of the solution search while avoiding improper local minima answers by using measurements of lens geometry and optics in conjunction with a constraining approach, (2) we have a clear global convergence and termination criteria, and (3) we ensured that our results are experimentally and computationally repeatable. In addition to these improvements, our LRT system is fully automated and has lens homeostatic control systems (temperature and pH). In this communication, we describe our methodology to obtain fast and accurate measurements of the lens GRIN profile using a LRT system.

LRT system
A custom dual chamber LRT system was designed and built to measure the optical changes in ex vivo bovine crystalline lenses over time. Previous studies have evaluated the lens optics or lens physiology independently [3,20]. Yet, it has been established that the optics and physiology of the lens are linked in an actively maintained equilibrium. This link is thought to be through lens's microcirculation system, which generates a circulating flux of ions and water throughout the lens [3,20]. Regulation of the steady state water content ensures that both the lens geometry and RI, which determines its overall refractive properties, are maintained [3]. Thus, measuring the water content, geometry and RI in the aging human lens is of great value and critical towards understanding age-related refractive error shifts, presbyopia and cataract.
The LRT system delivers laser beams sequentially across the crystalline lens and captures images of the individual rays passing through the lens using a front-facing camera. This imaging is performed for multiple delivery angles of the laser, known as projections. Generally, a greater number of projections and images per projection leads to a more accurate GRIN measurement, but this lengthens the duration of the experiment and significantly increases the size of the data set (detailed in Section 2.3.2).
The LRT system is controlled via custom-written MATLAB (R2015b, Mathworks, Natick, MA, USA) code for automated laser delivery and acquisition of ray trace images at every projection. The fully automated LRT system can acquire data from two independent lenses (one lens under normal conditions, the second lens under physiologically perturbed conditions) at pre-programmed time intervals. This allows us to examine changes in lens optics over time and hence study the relationship between the lens physiology and optics in real-time. A SolidWorks (SolidWorks Corp. Waltham, Massachusetts, United States) assembly of the LRT system is depicted in Fig. 1(A). The LRT system is mounted on an optical table (Newport Corporation, VW-3046-OPT-OT).

Laser delivery
The LRT system ( Fig. 1(B)) uses a 514 nm, 5 mW laser (Coherent StingRay Laser Diode Module) with high line uniformity. This wavelength was chosen for its high transmission and low absorption. A precision laser cut 0.2 mm aperture is positioned in front of the laser diode to restrict the width of the laser beam. The laser (Fig. 1, A9) is mounted on a motorized rotational stage (T-RS60A, Zaber Technologies, Vancouver, BC) that pivots about the crystalline lens to allow the rays to pass through the lens at various projections. The laser is positioned close to the specimen tank as shown in Fig. 1(C) to allow for the delivery of large projection angles (~ ± 40°).

Two-dimensional motorized linear stage
The rotational stage (Fig. 1, A10) is mounted on a vertical pillar that is fixed on two adjacently concatenated, motorized, one-dimensional linear stages (T-LSM200 and T-LSM025, Zaber Technologies, Vancouver, BC) (Figs. 1(A) and 1(C)). The 25 mm linear stage (Fig. 1, A11) allows the user to precisely position the laser beam in the meridional plane of the lens. This stage is used for fine adjustment of the beam to ensure proper alignment, detailed in Section 2.3.2. The 200 mm linear stage (Fig. 1, A12) allows the laser to move along the meridional plane of the lens, enabling the delivery of rays at multiple positions along the lens. This stage is used to perform the ray tracing scans.

Cameras
A side-facing digital SLR camera (Canon EOS 1100D, Tokyo, Japan with Canon 50 mm f1.8 lens) (Fig. 1, B3) is used to check the alignment of the laser beam along the meridional plane of the crystalline lens while the second front-facing digital SLR camera (Fig. 1, B4) is used to record images of the individual rays passing through the lens. The front-facing camera is also used to capture images of the lens geometry ( Fig. 1(D)). A diffuse light source (Dell model 1908FP LCD screen) (Fig. 1, B1) positioned behind the lens tank is programmatically controlled to switch on when capturing the lens geometry image, enabling visualization of the lens edges. The diffuse light creates a dark border around the lens producing high contrast images for segmentation. The LCD screen is then turned off for the acquisition of the laser projection images to allow for visualization of the laser beam.
The camera parameters are controlled through a custom-written Windows batch file, which was integrated into MATLAB. The batch file contains code that calls commands in the gPhoto2 software application based on a corresponding camera library (libgphoto2) [21]. The controlled parameters include ISO, shutter speed and aperture of the camera which determined the output image quality (i.e. high sharpness, high contrast, low noise, high brightness). For our experiments, 1600 ISO and 1/10(s) shutter speed were used, with a total scan time of ~35 minutes. For all experiments, a small aperture (f 2.8) was used as this yielded maximum image sharpness.
The cameras were corrected for optical distortion using published methods [22,23] which calibrate a camera using a planar pattern shown at a minimum of two different orientations. For our calibration, we fixed the camera location and moved a planar pattern (checkerboard pattern 8*10 with 18 mm square sizes, generated in MATLAB). A total of 10 photos of the checkerboard were captured at various tilts and distances for each individual camera. The images were then loaded into MATLAB's Single Camera Calibration toolbox and the radial distortion parameters (k1, k2, k3) were calculated. The radial distortion of any point denoted by x distorted and y distorted is related to the undistorted locations x and y through the following equations: (1 ) distorted y y kr kr kr where x and y are the undistorted pixel locations, k1, k2 and k3 are the radial distortion coefficients of a lens and r 2 = x 2 + y 2 .

Specimen tanks and minor components
The specimen tanks ( Fig. 2(A)) consist of a custom-made glass tank and a custom lens holder machined from a thin stainless steel sheet (Marine grade 304, 0.55 mm thickness, Steel & Tube NZ). The thin metal plate ensures minimal obstruction of the lens surface to measure the lens geometry. A 12 mm hole was machined in the center of the plate for placement of the crystalline lens. Threaded screws are positioned at the four corners of the plate to allow for fine height adjustment of the lens within the tank. The specimen tank contains the lens and is filled with the respective solution, either artificial aqueous humor (control) or perturbation media.
The two specimen tanks are placed inside of a larger glass tank which is positioned on a platform made from a 5 mm clear Perspex sheet (Fig. 1, A4). The platform rests atop 4 mounting posts and post brackets (Ø1.5" Mounting Post 1/4"-20 Taps L = 8, Ø1.5" Mounting Post Bracket, 1/4"-20 Taps, ThorLabs) ( Fig. 1, A5-A8). The larger outer tank is filled with water which is heated by a fully submersible heater (AquaOne, 55w heater) (Fig. 2, B1). The temperature of the water bath is monitored and controlled at 34°C. The LRT system is capable of imaging two lenses during an experiment (as seen in Fig. 1(D)), allowing us to measure changes in the optical properties of ex vivo bovine lenses under controlled and physiologically perturbed conditions [2].

Parameter retrieval from LRT images
This section describes the methods used to retrieve the input parameters for the tomographic GRIN retrieval, namely the lens geometrical data (Fig. 3) and laser pathway information (Fig.  4) [13].

Geometry retrieval
The unprocessed lens geometry image is corrected for distortion (see Section 2.1.3) using the MATLAB computer vision system toolbox. The camera calibration is performed using MATLAB's camera calibration toolbox and radial distortion (3-coefficient model) is corrected for in all images prior to post-processing. This ensures that the segmented ray paths and geometry have minimal errors in the ray deflection angles and lens shape.
Experimentally, the user can only approximate the horizontal placement of the lens, but in most cases, a slight tilt is present about its optical axis which is difficult to visually assess ( Fig. 3(A)). To quantitatively determine the amount of the tilt, the geometry of the lens is segmented using the growcut edge detection algorithm [24]. The rotation angle is then calculated by fitting a second-order polynomial to the lens anterior and posterior brushed surface data points, independently. The intersection points of the two polynomials are then found, which provides a quantification of the lens tilt. A normalized unit vector is calculated using the two intersection points (Eq. (3)) and used to calculate the rotation angle (Eq. (4)): where p1 is the point 1 [x, z] and p2 is point 2 [x, z], and where u is the horizontal unit vector [1,0] and v is given in Eq.
(3). The rotation is then performed using the rotation matrix: where, γ is the rotation angle. This results in an undistorted, horizontal image of the lens and ray tracing images. Once corrected for distortion and symmetry (Fig. 3(B)), the lens geometry is obtained by fitting the aspheric lens equation (Eq. (6)) which yields the radii of curvature and conic factor (Figs. 3(C) and 3(D)): where z is the optical axis, x is the equatorial axis, R is the radius of curvature (mm), and k is the conic factor. Since the nature of the aspheric lens equation is symmetric, the image of the lens must be rotated to be symmetric about the lens optical axis to ensure an accurate fit, as described previously.

Ray pathway segmentation
Custom-written MATLAB code is used to segment the ray trace images by detecting the centroids of the individual rays ( Fig. 4(A)). The original RGB images are first converted to grayscale. A median filter is then applied to the grayscale image to reduce image noise. The centroids of the ray are then determined by calculating the average of two values. The first value is the maximum intensity in the row of the image while the second value is the center of the largest change in gradients for the left and right boundaries of the ray. The centroids are subsequently fit using robust least squares (Figs. 4(B) and 4(C)).
To represent the path of laser beam before and after passing through the lens, four nominal points are required (i.e. A, B, C and D) shown in Fig. 4(D). Point A represents the points of intersection between an arbitrary orthogonal (to incident rays) plane and incident rays, B represents the intersection points of the incident rays with the posterior surface of the lens, C represents the intersection points of the refracting rays with the anterior surface of the lens, and D represents the points of intersection between an arbitrary orthogonal (to the optical axis) plane and the refracting rays exiting the lens. The A data points are determined by finding the average slope of all the incoming rays and then solving each incident ray's equation simultaneously and the orthogonal direction equation (z = mx + c, where m = 1/average incident ray slope, c is determined by choosing an arbitrary position of x and z).
Points B and C are obtained by simultaneously solving the ray and lens surface equations. The D points are determined by solving each refracted ray's equations and the orthogonal direction equation (y = c, where c is just at an arbitrary z point greater than the thickness of the lens) [13].

Lens preparation and positioning
Fresh bovine eyes were obtained from a local abattoir (Auckland Meat Processors, Auckland, New Zealand) and immediately transferred to the laboratory. Bovine eyes were dissected to obtain isolated crystalline lenses for organ culture experiments by the removal of the ciliary muscle, zonular fibers and vitreous. The lenses were incubated in isotonic artificial aqueous humor (ISO-AAH; NaCl 125 mM; KC1 4.5 mM; MgCl2 0.5 mM; CaCl2 2 mM; NaHCO3 10 mM; glucose 5 mM; Sucrose 20mM; buffered with 10 mM HEPES to pH 7.1). For the purposes of this study, only the control lenses maintained under physiological conditions will be discussed. A total of 4 bovine lenses were used.
The lens was handled with a custom-made glass scoop and positioned on the lens holder with its anterior surface facing upwards for least obstruction of the lens geometry. The lens was visually adjusted to minimize tilt. Two drops of milk (Meadow Fresh Original) were added to increase the scatter of the solution and therefore increase the visibility of the laser beam.

Laser alignment & running the LRT program
Once the lenses were placed in the LRT system, the "start" and "end" positions for the laser scan were selected manually by the user, prior to start of the experiment. The position of the laser beam along the lens was adjusted using the 200 mm linear stage. The "start" position is the position where the beam is ~1 mm from the equatorial edge. This is defined as the beginning position of the scan. The "end" position is the position where the beam is ~1 mm from the opposing equatorial edge and is defined as the last position of the scan. Thus, the LRT system scans the entirety of the lens with the exception of the outermost 2 mm in order to avoid scattering of the beam at the lens edge. The number of rays in the scan is defined by the user and is independent of the "start" and "stop" position. The MATLAB code requests the total number of laser rays, number of projections, corresponding projection angles, time between each scan sequence and the total scan time.
At the start of the experiment, the side-facing camera was used to capture images of the laser beam at three positions along each lens (the user defined "start", "end", and the halfway point between the "start" and "end" positions) to ensure that the laser was passing through the meridional plane. When the laser is properly positioned, there is minimal deflection of the outgoing ray and no lateral shift of the beam at the three positions. The alignment images were immediately segmented to extract the individual ray paths and calculate the alignment error (difference in angle between the incoming and outgoing rays). The 25 mm linear stage was adjusted by the user, if needed, to minimize the error (i.e. correct for meridional offset) without physically moving the lens. This process was repeated for each delivery angle (e.g. 0°, 10°, etc.) to ensure that the laser was properly aligned on the lens at all of the projections and that the alignment error was within a user-defined tolerance (< 0.5°).
Following the alignment verification procedure, the LRT system proceeded to deliver parallel equally-spaced rays across the lens at the designated projection angles. In previous experiments using porcine lenses [13], it was found that 71 rays retrieves the RI to two decimal places for a ~10 mm lens. To maintain a similar ray interval, we delivered 151 rays across each bovine lens. Vazquez et al. [25] showed that increasing the number of rays had little influence in improving the accuracy of the GRIN retrieval. However, imaging at more projection angles produced a better sampling of the RI distribution.
An image of the lens geometry was then recorded with the front-facing camera to obtain the lens shape measurement. A lens geometry image was captured at the start of each userdefined scan interval to monitor changes in the lens shape over time. Once the LRT experiment was complete, an image of a spherical lens (Edmund Optics, 10.0mm Diameter, N-BK7 Lens) was captured for calibration. Since the position of the crystalline lens relative to the front-facing camera varies slightly between experiments, an image of a spherical lens is captured and used for image calibration. This was achieved by image segmenting the spherical lens and finding the number of pixels in the optical axis (z) and equatorial axis (x).
The number of pixels was then divided by the lens diameter to obtain the z and x scaling factors (in pixels per millimeter) along both axes. All subsequent segmentation was converted to mm by dividing the measured number of pixels by the scaling factor.

Constraining of the inverse problem
The images recorded during the LRT experiment were saved and stored for post-processing. Following the geometry retrieval and ray pathway segmentation detailed in Section 2.2, the relevant parameters were extracted from the experimental data and the lens optical measurements were then associated to the GRIN [13]. More specifically, the experimental exit deflection angles measured by LRT were converted to experimental optical path length (OPL) measurements (Eq. (7)). Using theoretical ray tracing [26] through the lens where the first iteration is assumed to be homogenous, a theoretical OPL was calculated and minimized against the experimental OPL. The minimization process led to a sequence of solutions that were iteratively more accurate than the previous. The image segmentation (Section 2.2) generated a system of equations which had a general form of Ax = b, where A is composed of the line integrals of the rays through the lens from theoretical ray tracing, b is the experimental OPL calculated by the exit deflection angles (Eq. (7)) and x is the vector of GRIN coefficients.
Analysis of the line integral matrix A, indicated our system of equations was illconditioned, meaning the mathematical solution was highly sensitive to the values (or more precisely the errors in those values) of the coefficient matrix, A, or the right-hand side vector, b, in the inverse problem Ax = b. This is common in inverse problems (A −1 Ax = A −1 b) and is usually difficult to deal with mathematically. In our case, the system of equations was an over-determined matrix (more equations than unknown variables), where the matrix condition (a measure of the severity of the ill-conditioned problem), was large ( 8 ( 1 10 ) > × ). This was much greater than unity (1) which is assigned to a well-conditioned inverse problem. When using classical methods that solve inverse problems, such as orthogonal triangular decomposition and unconstrained least squares, significant discrepancies arose between the experimental solution and the expected solution due to experimental error. Hence, we mathematically constrained our solutions space to only obtain solutions that are biologically relevant. In addition, this had the effect of preventing the problem from falling into erroneous local minimum solutions, as they were excluded from the solution space.
The following sections detail the techniques that were applied to retrieve the GRIN distribution considering the lens geometry. Briefly, this involves (1) smoothing the exit deflection angle values, smoothing of errors in the calculated optical path length (Section 2.4.1) and, (2) applying a series of non-conflicting and non-restrictive constraints on the solution space (Section 2.4.2).

Smoothing constraints
As the GRIN retrieval algorithm uses only the intersection points along the lens surface (B and C, Fig. 4(D)) and the exit deflection angles (calculated using points C and D, Fig. 4(D)), these are the parameters that were smoothed [13]. The smoothing was achieved using a ninthorder polynomial with robust bi-square weights. The bi-square weights method minimizes a weighted sum of squares, where the weight given to each data point depends on how far the point is from the fitted line. Points farther from the line were assigned a reduced weight, and points further from the line than expected got zero weight. This method was preferred over the original cubic spline method due to errors in the experimental measurements and errors in the image segmentation. The use of polynomials was also beneficial as concerns such as the over-fitting of data, commonly known as Runge's phenomenon, were avoided [27].
In addition, custom-written MATLAB code was developed to remove "stray" rays in the image post-processing using two methods. A stray ray constitutes any ray that does not show continuity with the adjacent rays and are generally caused by air bubbles, attached vitreous and/or the lens sutures. The code allows the user to visually assess the incoming and outgoing rays and remove any rays that are refracting irregularly. The second method for removing stray rays is an objective approach, where the rays which are overlapping or cross at the anterior refracting surface are excluded. The stray rays were removed by implementing zeros in the matrix A and vector b in the system of equations Ax = b. The remaining rays were then used to calculate the optical path lengths: where S(BC) is the optical path length, α is the exit deflection angle, and x is the equatorial axis. The optical path lengths were then smoothed using a robust fit as the deflection angle parameter, α , and the x-coordinate parameter contained experimental error. This was justified by the continuous RI that has been found in previous studies [28][29][30].

GRIN constraints
Next, physiologically relevant GRIN constraints were mathematically applied to limit the solution search space of our minimization problem. To implement these constraints, the MATLAB built-in function lsqlin was used. Specifically, Ax = b describes the minimization problem min||Ax-b|| and A ineq. x ≤ b ineq describes the inequality constraints (Eq. (8).
Here, A is the matrix containing information about the line integrals, b is the vector of optical path lengths, x is the vector of GRIN coefficients solved for, and A ineq and b ineq contains all the constraints discussed in this section. These constraints were imposed simultaneously with the bipolynomial equality constraints (A eq. x = b eq ) when solving for the bipolynomial model [13].
Physiologically, there are three sensible constraints that do not restrict the solution (Fig.  5): (1) the GRIN profiles along both the equatorial and optical axes observes a second-order profile [28][29][30] and the RI values must decrease from the lens core to the outer cortex. This was the primary constraint used in our algorithm, as it was the least restrictive and found to be the most effective in modelling fiber cell layers. This constraint was implemented as a second-derivative of the RI in the Z direction < 0 (i.e. the profile must be concave down in the Z direction) and the second-derivative of the RI in the X direction is < 0 (Fig. 5(A)).
(2) The RI at the surface of the lens must be within a physiologically appropriate range, and the RI values of the surface points must also not vary significantly for GRIN continuity. This was deemed justified due to the lens capsule having nearly the same RI due to similar water/protein ratio. The constraint was implemented as the surface RI values must not have a difference greater than 0.005 of an index value. Furthermore, to exclude more biologically irrelevant solutions from the solution space, the RI values at the surface of the bovine lens were set to < 1.39, as no values higher than this has been observed in literature [28][29][30] (Fig. 5(B)).
(3) The final constraint limited the maximum and minimum RI values of the lens domain. For a normal bovine lens, the maximum index is approximately 1.45 and the minimum is approximately 1.36. These values indicate the approximate range observed in previous MRI and LRT experiments [2, [28][29][30]. This was constrained so that the maximum and minimum values of all RI in the lens domain were to be within an appropriate biological range (e.g. 1.34 < n internal < 1.47). In the absence of any knowledge of the index range, a larger range could be set. However, it must be noted that with a larger range, the solution space increases and thus the number of possible solutions. We noticed that although relevant, in most cases this constraint was not required, as the previous two constraints produced values within a physiologically acceptable range and usually produce biologically plausible solutions. These three constraints were applied simultaneously by concatenation of each of the separate constraints. The size of the constraints matrix depended on the lens domain size and the number of points chosen in the domain. Limited by computer memory (8 GB RAM), we applied the constraints to 10% of the points across the lens by evenly sampling along the optical axis (Z direction). Using more points would linearly increase the amount of computer memory required due to an increase in the size of the matrix A. No difference in RMS error was observed when using 5% or 10% of the total number of domain points in our simulations. However, we chose to use a larger number of points as computation times did not significantly increase. The size of the total constraint matrix size was approximately 2 × 10 6 rows by, nine model coefficients + number of projections, columns.

Results
Using the LRT system described above, we imaged four ex vivo bovine lenses in AAH and retrieved their individual shape and GRIN profiles. Our measurements were highly reproducible as normal (control) bovine lenses have similar GRIN profiles. The RI values along the equatorial axis of the control bovine lenses were extracted from calculated GRIN maps and compared with previous MRI work [2].

LRT repeatability
The repeatability of the technique for retrieving the optical parameters of the lens was assessed by performing multiple, independent, sequential scans. A single bovine lens was positioned in the LRT system and imaged every hour for t = 0 to 4 hours (5 data sets in total). The lens remained in the same position for all of the scans. The lens shape parameters (anterior and posterior radii of curvature, conic factors, and thickness) and the minimum, maximum and average RI values were analyzed for the five scans. The results are shown in Table 1.
The results indicate high reproducibility in the RI measurements with a maximum SD of 0.01. Relatively larger variances were observed in the conic parameters due to the high sensitivity of the parameter and interdependence with the radius of curvature parameter. Both radii of curvature and lens thickness parameters were shown to also be highly repeatable, and the small variances could be attributed to image segmentation noise and slight changes in the lens shape over time.
where n is the refractive index, r is the radial distance from the center, p is the equatorial radius. The Pierscionek et al. study took into consideration the size and age differences in the lens measured. Figure 6 shows an overlay of the RI of the bovine lenses retrieved using our LRT system and the RI distribution calculated by the above model. Using different mathematics and wavelengths, but the same ray tracing technique, the Pierscionek's model and our LRT system produce similar results. In nearly all regions of the lens, the RI difference is < 0.02 across the equatorial plane (Fig. 6). The green plot indicates the GRIN profile measured using MRI. Good consistency is observed in the cortex regions of the lens but larger discrepancies in the RI are observed at the core region. Accurate retrieval at the lens core using MRI is difficult due to low signals.

Comparison with MRI
We compared our retrieved GRIN profile of a bovine lens with those from the MRI study by Vaghefi et al. (2015) [2]. Comparing these results, LRT obtained values of ~1.370 at the surface of the lens and ~1.440 peak RI (Fig. 6), while MRI produced ~1.380 and ~1.435 respectively. The difference in RI between the lens surface and peak RI is relatively similar (1.440-1.370 = 0.070 for LRT, 1.435-1.380 = 0.055 for MRI). However, a large mismatch in the absolute RI at the core of the lens is observed (Fig. 6).

Regional differences
An analysis of the regional differences in GRIN was conducted to quantify and better visualize differences in the GRIN profile by extracting the mean RI in three different regions: the outer cortex, inner cortex and core. The outer cortex is defined as r/a = −1 to −0.75 and r/a = 0.75 to 1, the inner cortex is defined as r/a = −0.75 to −0.5 and r/a = 0.5 to 0.75, and the core is defined as r/a = −0.5 to 0 and r/a = 0 to 0.5, where r/a is the equatorial axis of the lens normalized from −1 to 1. The results are summarized in Table 2.
The regional difference analysis shows the greatest similarity of RI values in the inner cortex of the lens. The ray tracing based methods are similar at the core; however, the MRI measurements appear to be considerably lower in comparison. Some discrepancy in all three methods is observed in the outer cortex region. Table 2. Mean RI in different regions of the lens using three different measurement methods. The outer cortex is defined as r/a = −1 to −0.75 and r/a = 0.75 to 1, the inner cortex is defined as r/a = −0.75 to −0.5 and r/a = 0.5 to 0.75, the core is defined as r/a = −0.5 to 0 and r/a = 0 to 0.5.

Discussion
A custom-built, fully-automated dual chamber LRT system was developed to allow us to quantify the optical changes occurring in ex vivo bovine lenses over time. A novel mathematical method was implemented to ensure that the inverse ill-conditioned system of equations could be solved, whilst retrieving a biologically valid solution. Our GRIN retrieval algorithm required approximately five iterations to sufficiently retrieve the GRIN distribution, corresponding to a total time of ~10 minutes. The largest number of iterations for convergence (RMS (current iteration -previous iteration) < 1*10 −5 ) observed was eight. Furthermore, this method is more robust to experimental errors due to the implementation of biological constraints. The GRIN retrieval algorithm is inherently more robust to experimental errors because it excludes all biologically non-relevant solutions.
The GRIN data obtained from our LRT system were compared to published LRT and MRI data [2,28]. There is a small difference in the equatorial RI (<0.02) between the published equatorial GRIN profile and our retrieved equatorial profile at the lens surface (r/a = 1), which is likely due to the difference in the models used for fitting. Pierscionek et al. (1989) used a second-order model that does not have a rotationally symmetric profile about the optical axis (Eq. (7). Conversely, our GRIN model is rotationally symmetric about the lens optical axis [13]. Another reason for the observed difference is an inherent setback of LRT, as it is difficult to scan the outermost region of the lens (close to the equator) because the rays exhibit more scattering than refraction. As a result, when our algorithm solves for the lens GRIN, the outer cortex RI values are extrapolated using the central region of the lens.
When comparing our results with MRI, the largest discrepancy in the GRIN profile was at the core region of the lens. We suggest this difference was attributed to: (1) a difference in the measurement modalities, (2) inherent noise in MRI modality and (3) a limitation of the MRI measurement modality, where in general low levels of signals are acquired at the core of the lens. The discrepancy in the outer cortex may be due to the partial volume effect, or partial volume artifact which is caused when the surface voxels contain both lens crystalline and the surrounding solution; thus, the signal is an average RI. Nevertheless, comparison with both previous works showed a difference of < 0.02 of a RI, indicating good comparability.
The primary limitation of our GRIN retrieval technique is that it requires information about the refractive properties of the scanned sample to minimize the effect of noise in the data. This information includes the maximum and minimum possible RIs of the lens and the RI range at the lens surface. Based on the quality of the data set, these limits may not be required. However, this estimation is essential to find a good corresponding solution with noisy data sets. In our testing, implementation of the first and second constraints ensured that the RI distribution was biologically appropriate in most instances. The third constraint often had no influence over the final solution. However, the third constraint was still implemented to ensure that the optimization process would always find a biologically valid solution, even with extremely noisy experimental data sets.
Future work will involve using the developed LRT system and GRIN retrieval algorithm to investigate how the lens optical properties change under different physiological perturbations [31][32][33]. This information is critical towards understanding the age-related physiological and optical changes that occur in the lens, which lead to refractive error, presbyopia and cataract. The LRT system will be used to quantitatively investigate the effects of pharmacological modulations of lens physiology on its optics, which could be important in preventing age-dependent changes in vision quality.

Funding
The study was supported in part by Whitaker International Scholar Award (BMH), New Zealand Royal Society Marsden project #3706540, New Zealand Optometry and Vision Research Foundation #3700559, New Zealand Optometry and Vision Research Foundation #3707089, Maurice and Phyllis Paykel Trust #3707893.