Estimation of the full shape of the crystalline lens in-vivo from OCT images using eigenlenses

: Quantifying the full 3-D shape of the human crystalline lens is important for improving intraocular lens power or sizing calculations in treatments of cataract and presbyopia. In a previous work we described a novel method for the representation of the full shape of the ex vivo crystalline lens called eigenlenses , which proved more compact and accurate than compared state-of-the art methods of crystalline lens shape quantification. Here we demonstrate the use of eigenlenses to estimate the full shape of the crystalline lens in vivo from optical coherence tomography images, where only the information visible through the pupil is available. We compare the performance of eigenlenses with previous methods of full crystalline lens shape estimation, and demonstrate an improvement in repeatability, robustness and use of computational resources. We found that eigenlenses can be used to describe efficiently the crystalline lens full shape changes with accommodation and refractive error


Introduction
The human crystalline lens is a sophisticated optical element in the eye with rather unique properties.In the young eye, the crystalline lens is able to change its shape dynamically to focus near and far objects onto the retina (accommodation) [1].With aging, the crystalline lens loses its accommodation capacity (presbyopia) and transparency (cataract).Standard treatment for the correction of these conditions involve the replacement of the crystalline lens material by an artificial intraocular lens (IOLs) that is placed inside the capsular bag [2].Accurate quantification of the shape of the crystalline lens is important to help selection of the IOL.For instance, selecting the power of standard and of emerging accommodating IOLs requires an accurate estimation of the IOL position after surgery which is related to lens shape [3][4][5].Also, selecting the size of accommodating IOLs can be critical for its correct mechanism of action [5].
Magnetic Resonance Imaging (MRI) [6][7][8][9] and Ultrasound Biomicroscopy (UBM) [10] techniques have been demonstrated to image the full shape of the lens including the peripheral region behind the iris.However, typically the resolution of these techniques is low (approximately 100 µm in MRI and 40 µm in UBM versus <10 µm in Optical Coherence Tomography, OCT) and the MRI acquisition time is long (several minutes in MRI versus miliseconds in UBM and OCT).Furthermore, UBM requires contact with the eye, and MRI is costly and not available in a typical ophthalmology clinic setting.Optical imaging techniques, including Purkinje imaging [11,12], Scheimpflug imaging [13][14][15][16] and OCT [17][18][19][20] have been demonstrated for quantification of crystalline lens geometrical parameters.Optical techniques, however, have two main drawbacks: (i) they suffer from refraction distortions from preceding surfaces, which need to be corrected, and (ii) the iris blocks the incident light and thus only the shape of the lens within the pupil is available.
Refraction correction has been successfully achieved using ray tracing methods [21,22].On the other hand, extrapolation of the surfaces blocked by the iris from data visible through the pupil has often been achieved using the so-called intersection approach, which extrapolates the lens periphery using circle fitting of the anterior and posterior surfaces of the lens, with the equatorial diameter obtained from their intersection.This approach has been reported in some lens studies [23][24][25][26] and is used in some commercial OCT systems such as the CASIA2 (Tomey, Nürenberg, Germany) and the Catalys laser (Johnson & Johnson Vision, Santa Ana, CA, USA).In a previous work on ex vivo lenses we showed that the intersection approach overestimates the equatorial diameter and the volume of the lens, and underestimates (anterior shift) the equatorial plane position [27].
In that earlier study [27] we proposed a 2-region parametric model (2RM) to estimate realistically the full lens shape in vivo from OCT images.In brief, the method was trained and validated using OCT images from ex vivo isolated lenses, in which the full shape of the lens is visible because the iris is removed before the measurements [28].The method was demonstrated in vivo on accommodating young subjects [29], and in cross-sectional studies where the crystalline lens was quantified as a function of age [30] or refractive errors [31].The method has also been used to improve the estimation of the intraocular position in a cataract surgery [3] and to quantify the crystalline lens in a guinea pig myopia model in vivo [32].While the 2RM method allowed us to fully quantify the crystalline lens in human and animal models, the description of the lens required merging a set of parametric surfaces that described anterior, posterior and peripheral parts of the lens, making its practical implementation complicated and requiring additional user interaction to monitor the constructed models.
Recently [33] we have presented a new method for the representation of the full shape of the crystalline lens ex vivo that we called eigenlenses.Briefly, eigenlenses represent the most common "deformation patterns" (with respect to an average lens shape) among the human crystalline lens shapes that can be found in nature.Thus, for example, the first eigenlens (the most "common" deformation) is related with changes in the size of the lens, while the second eigenlens is related with changes in the aspect ratio (i.e., lenses more rounded or "stretched").Eigenlenses representation is compact, capturing the 96% of the variance in shape of the training set of lenses with only 6 coefficients.
In this manuscript we demonstrate the use of eigenlenses to estimate the full shape of the lens in vivo, which is made possible thanks to the compaction ability of eigenlenses.Unlike the 2RM, eigenlenses method provides, by construction, a meaningful shape, is more repeatable and robust, and requires lower computational resources.We also demonstrate eigenlenses ability to describe the full shape crystalline lens changes with accommodation and its dependence on refractive error.

Subjects
Seventeen eyes from 17 young human subjects (6/11 males/females; mean age: 28.6 ± 3 y/o, age range: 22 to 30 y/o; spherical error range: -6.75 to + 0.75 D) were studied.Subjects signed a consent form approved by the Institutional Review Boards after they had been informed on the nature and possible consequences of the study, in accordance with the tenets of the Declaration of Helsinki.

Optical coherence tomography imaging system and experimental protocol
Images of the anterior segment of the eye were acquired with a custom swept-source optical coherence tomography (SS-OCT) system employing a MEMS-based vertical-cavity surfaceemitting laser (VCSEL) swept-source (SL132120, Thorlabs, USA), centered at 1300 nm and a fiber-optics based Mach-Zender interferometer configuration.The system had been described in an earlier publication [34].The acquisition speed was 200000 A-Scans/s.The axial depth range was 15.95 mm in air, sampled by 1920 pixels, resulting in a pixel size of 8.3 µm.The measured axial resolution was 16 µm in air.The transverse field-of-view was 15 mm × 15 mm, sampled by 300 and 150 pixels in the horizontal and vertical (fast and slow) directions, resulting in a pixel size 50 µm × 100 µm, respectively.
Images of the cornea and crystalline lens were acquired for relaxed accommodation (0 D) and 4.5 D of accommodative demand.To induce accommodation and compensate defocus, an external accommodative channel, based on a Badal system mounted on a manual stage and a Digital-Light-Processing (DLP) pico projector (PicoPix, 854 × 480 pixels, Philips NV, Amsterdam, Netherlands; 55 lum), was incorporated into the OCT system through a dichroic mirror.The fixation stimulus consisted of a white Snellen E-letter presented on a black background.The subject's spherical error was corrected using a trial lens placed on a pupil conjugate plane and subject's stabilization was guaranteed using a bite bar.Prior to imaging, the subject's eye pupillary axis was aligned to the optical axis of the instrument following the procedure explained in prior work [19].Repeated 3D-OCT measurements of the right eye of each subject (around five per subject) were acquired for each accommodative demand.

Eigenlenses and eigencenters
The eigenlenses approach describes the shape of any crystalline lens as an average lens plus a linear combination of a K number of eigenlenses ("deformation patterns") [33]: where l n represents the lens n to be represented, l the average lens obtained from a set of ex vivo training lenses, e k the eigenlens k and a k the coefficients that multiply each eigenlens (establishing the contribution of each deformation pattern) to obtain the lens l n .Elements of vectors l n and l represent the distances r from the origin of coordinates to the surface at every pair of sampled elevation (θ) and azimuth (φ) angles in spherical coordinates, i.e., l n = (r θ 1, ϕ 1 , . . ., r θ P, ϕ Q , ), with P the number of elevation angle sampling points, Q the number of azimuth angle sampling points, and M=P × Q the number of points that define the lens [33].Figure 1(a) represents the 6 first eigenlenses, and Fig. 1(b) shows an example of how the representation of a given lens improves as more eigenlenses are included.Figure S1 shows an example of the construction of an ex-vivo lens of age 65 y/o with eigenlenses.
Furthermore, eigenlenses were additionally constructed to represent the optical zone (central area) of the crystalline lens that is visible through the pupil in-vivo.In the current study we introduce the term eigencenters as a denomination of this representation.The concept of the eigencenters is similar to the eigenlenses, but in this case, they are obtained using only the central area of a given diameter in the training ex-vivo lenses, and thus, they represent the variations of the shape of the lens within the pupil.Therefore, the central area of any crystalline lens  can be represented as the average crystalline lens-central area plus a linear combination of K eigencenters: where z n represents the central area of lens n to be represented, z the average shape within the central area from a set of training lenses, g k the eigencenter k and c k the coefficients that multiply each eigencenter to obtain the central area of lens n. Figure 1(c) shows the first 6 eigencenters representing the most common "deformation patterns" in the central area.
In our previous work [33], we determined that K = 6 was the optimal number of eigenlenses/eigencenters to be used, as we found that a higher K did not result in a substantial increase of the representation accuracy (with K = 6 we were able to explain the 96% of the variance in the training set).Therefore, we will also use 6 eigenlenses/eigencenters in the present work.

Training a model for the estimation of the full shape of the crystalline lens from its central area
Eigenlenses can be used to estimate the full shape of the crystalline lens from the information of the optical/central area that is available in-vivo.To this end, a mathematical expression for the estimation of full-shape coefficients a k from the eigencenter coefficients c k was obtained, following the next steps: (i) Eigenlenses representation was calculated in a training set of N=133 isolated lenses, obtaining the corresponding set of a k coefficients {a k , k = 1, . . ., 6} for each of the 133 lenses [33].Furthermore e k and l were stored as they are needed to obtain the full shape of new (previously "unseen") test lenses.
(ii) In vivo conditions were simulated in the training set assuming that only the central part of the lens was visible.Specifically, experiments were done assuming pupil diameters from approximately 3.4 to 6.5 mm (sampling the elevation angle in polar coordinates from π/4 to π/8, see Fig. 3 in [33] for details).From this central part, the eigencenter representation was calculated for each lens, obtaining the corresponding set of c k coefficients {c k , k = 1, . . ., 6}.Furthermore g k and z were stored as they were needed to obtain the eigencenter representation of new test lenses.
(iii) An expression for the estimation of coefficients a k from c k (a k f (c 1 , . . ., c 6 )) was obtained using multivariate linear regression.Specifically, we solved the problem: where A is an N×6 matrix where each column corresponds to a coefficient a k and each row to a training lens, C is an N×(6+1) matrix where each column k corresponds to a coefficient c k (first column is a vector of ones) and each row to a training lens, β is a 7×6 coefficients matrix (to be obtained) and θ is the error matrix [35].The matrix β contains the coefficients needed to estimate each a k from a set of c k eigencenter coefficients.
Figure 2 illustrates the process.From the whole set of 133 training lenses, we obtained and stored e k and l, and, after simulating in vivo conditions, g k and z.Then, for each ex vivo training lens, we obtained a k (representation of the full shape) and c k (representation of the central area) coefficients.Finally, β was estimated solving Eq. (3).With this information, we can estimate the full shape of a new, previously unseen, test lens in vivo as explained in the next section.

Full shape estimation in-vivo using eigenlenses
In this section, we describe the estimation of the crystalline lens full shape from a volume of OCT images of a crystalline lens in vivo using the models previously trained (section 2.4).Specifically, the method consists of three steps: (i) the 3-D shape of the lens is obtained within the pupil; (ii) the eigencenter representation of this 3-D shape is obtained; and (iii) the full shape coefficients and, from them, the full shape of the lens is estimated.

Obtaining the 3-D shape of the lens within the pupil
Obtaining the accurate 3-D shape of a lens within the pupil from the OCT images involved automatic surfaces detection and segmentation, and distortion correction [21,36].The surfaces of interest (anterior and posterior cornea and anterior and posterior crystalline lens) were automatically segmented in each B-scan using thresholding, Canny edge detectors, mathematical morphological operations, and a priori knowledge of the measurements [19,37].Then, each 3-D surface composed of all B-scans segmentations was fit with Zernike polynomials up to the 4-th order, and the resulting smooth surface was used to guide a new refined segmentation in each B-scan.In this case, the edges were searched in the local neighborhood of the "smooth segmentation", obtained (for each B-scan and surface) by evaluating the Zernike representation of the corresponding 3-D surface at the corresponding coordinates of that B-scan [28].
The custom-designed OCT system produced a partially-telecentric scan.The residual "fan distortion" arising from the scanning architecture and the optics of the system was corrected following the process described by Ortiz et al. [36].Optical distortion, due to the refraction of the rays in the different surfaces of interest, was corrected by using 3-D ray tracing [21].The corneal group refractive index was taken as 1.389 [38], the aqueous humor group refractive index as 1.345, and the crystalline lens refractive index was obtained from the age-dependent average refractive index expression derived by Uhlhorn et al. [39].Figure 3 illustrates the process to obtain the 3-D model within the pupil.

Obtaining the eigencenter representation of the 3-D shape
The eigencenters representation of the 3-D shape of the lens within the pupil was obtained by calculating the corresponding coefficients c k .Specifically, the 3-D crystalline lens data were first centered laterally at the anterior lens maximum elevation point and axially at the mid-point between anterior and posterior surfaces, converted into polar coordinates, and interpolated in the specific M = P×Q=10000 defined sampling points, obtaining z.The residual data was obtained by subtracting the mean lens, and projected into the eigencenter basis: T , where c is a vector with the 6 coefficients, c = (c 1 , . . ., c 6 ) T , E is an M×6 matrix where each column is the ordered (starting form the highest eigenvalue) first 6 eigencenters E = (g 1 , . . ., g 6 ), and z is an M×1 vector of the test data; T indicates the transpose operator.

Estimation of the full shape coefficients and the full shape of the lens
From the set of eigencenter coefficients c = (c 1 , . . ., c 6 ) T , full shape coefficients â = (â 1 , . . ., â6 ) are estimated using the β matrix obtained in the training phase (as illustrated in Fig. 2 where â = (â 1 , . . ., â6 ) are the estimated coefficients, c ′ = (1, c 1 , . . ., c 6 ), and β is a 7×6 coefficients matrix.Once we have â, we can obtain the full shape l by projecting into the eigenlenses basis and adding the mean lens: Figure 4 illustrates the process to estimate the full shape from the 3-D model within the pupil.

Full shape quantification
Once the full shape 3-D models were constructed, various biometric parameters were quantified.Specifically, we obtained: (1) lens thickness (LT); (2) lens equatorial diameter (DIA); (3) lens equatorial plane position (EPP); (4) lens surface area (LSA); and (5) lens volume (VOL).The EPP was defined as the distance between the anterior lens apex and the position of the equatorial plane [3].The VOL was estimated by numerically solving the double integration of the anterior and posterior lens surfaces [29].The LSA was estimated from the Delaunay triangulation of the 3-D anterior and posterior lens contours [29].

Data analysis
We investigated if the âk coefficients and the geometrical parameters (DIA, VOL, LSA and EPP) varied significantly with the pupil diameter (PD) from which the full shape was estimated, applying one-way ANOVA with Bonferroni's correction for multiple comparisons.
We analyzed the relationship between eigenlenses coefficients and full shape geometrical parameters by means of linear regression, obtaining the Pearson correlation coefficient (ρ) and the p-value for testing the hypothesis of no correlation (p).
We also analyzed the agreement between the proposed method and our previous method for the full shape estimation of the crystalline lens [27] using Bland-Altman plots [40].To evaluate the repeatability of each method, we calculated the standard deviation across repeated measurements for each subject, and the mean across subjects was obtained.
We compared the mean values of geometrical parameters and eigenlenses coefficients between 0 D and 4.5 D accommodative states using a two-tails paired t-test.For all analyses, statistical significance was defined as a p-value lower than 0.05.Calculations were performed with Matlab software (MathWorks, Natick, MA, USA).

Pupil diameter analysis
The âk coefficients and the geometrical parameters (DIA, VOL, LSA and EPP) estimated with the proposed method were analyzed as a function of the pupil diameter (PD) from which the estimation of the full lens shape is obtained (n=17 subjects at 0 D of accommodation, mean of the repeated measurements).We investigated coefficients and parameters for PDs between 3.4 to 6.5 mm.In cases in which the actual PD of the measurement was higher than the tested PD, we removed the information outside the tested diameter.
We did not find a statistically significant difference of the â1 , â2 , â3 , â4 and â5 coefficients and the reconstructed lens geometrical parameters across pupil diameters (p>0.05,one-way ANOVA with Bonferroni's correction).Furthermore, the values did not differ more than 1% for PD larger or equal to 4.4 mm with respect to the largest PD tested (6.5 mm) for any parameter except for the â1 (that differed <2%) and â3 (that differed <4%).Only the â6 coefficient showed a statistically significant variation with PD.These results are shown in Figure S2.

Relationship between eigenlenses cofficients and full shape geometrical parameters
Figure 6 shows the correlation between eigenlenses coefficients and full-shape geometrical parameters.The first column shows the correlation of DIA (first row), VOL (second row), LSA (third row) and EPP (fourth row) with the first coefficient â1 ; the second column with the second coefficient â2 ; and the third column with the fifth coefficient â5 .The black solid line represents the best linear fitting, which is plotted when the correlation is statistically significant.Correlation coefficient (ρ) and the p-value for testing the hypothesis of no correlation (p) are also shown.
For these calculations we considered a PD of 5.2 mm, which was the closest value to the actual mean PD of the in-vivo measurements.These results are for measurements obtained at 0 D of accommodative demand (n=83, including repeated measurements on the 17 subjects).
In general, VOL and LSA strongly depended on â1 , while DIA depended on a combination of â1 , â2 and â5 .As expected, we did not find correlation of â3 or â4 eigenlens coefficients (asymmetric terms) with any of the geometrical parameters, and of â6 coefficient with any parameter, except for the EPP (ρ =-0.51, p=8•10 −7 , not shown in the figure).Table 1 shows a multiple linear regression analysis using least squares between the crystalline lens geometrical parameters and the first two coefficients â1 and â2 .First row shows the R 2 , second row shows the equation that relates each parameter and the two coefficients and third row shows the Mean Absolute error (MAE) between the estimation using the two coefficients and the actual parameter.LT, DIA, VOL and LSA can be accurately estimated from the first 2 coefficients.On the contrary, EPP estimation improves if the fifth coefficient is included instead of the first coefficient (R 2 = 0.90, EPP = 1.841 + 0.0167â 2 + 0.0194â 5 , MAE = 0.025 mm).

Comparison with our previous 2-region parametric model (2RM) to estimate the full lens shape in vivo
Figure 7 shows the Bland-Altman plots for measurements at 0 D accommodation (n = 17, using the mean of the repeated measurements for each subject).The Y axis represents the difference between the methods calculated as 2RM-proposed eigenlenses.The limits of agreement (LoA) are calculated as 1.96 • SD, where SD is the standard deviation of the difference.Table 2 shows the mean value ± standard deviation across subjects of each method (2RM and proposed eigenlenses) and geometrical parameter.In general, the 2RM method underestimated the size of the lens relative to the proposed eigenlenses method (∼2% lower DIA, 3% lower VOL, 1% lower LSA) and overestimated EPP (7%), especially for smaller lenses.
Table 3 shows the repeatability of each method (2RM and proposed eigenlenses) for each geometrical parameter (mean across subjects of the standard deviation across repeated measurements).The standard deviation across measurements is statistically significantly lower (better

388
Table 3 shows the repeatability of each method (2RM and proposed eigenlenses) for each 389 geometrical parameter (mean across subjects of the standard deviation across repeated 390 measurements).The standard deviation across measurements is statistically significantly lower 391 (better repeatability) with the eigenlenses than with the 2RM method for all geometrical 392 parameters except for the VOL (paired t-test, asterisks in Table 3).repeatability) with the eigenlenses than with the 2RM method for all geometrical parameters except for the VOL (paired t-test, asterisks in Table 3).The mean computational time across lenses was statistically significantly lower (paired t-test) with eigenlenses (4.3 ±0.4 ms) than with the 2RM method (8.1±2.3 s), measured in a 3.5 GHz CPU, 8 GB RAM, using Matlab software.Table 4 shows the mean changes across subjects of DIA, VOL and eigenlenses coefficients (â 1 , â2 , â5 and â6 ), between the two accommodative demand states, and the p-value of a paired t-test to compare if the mean changes with accommodation were statistically significant.As expected, asymmetric terms â3 , â4 (not shown) and VOL did not statistically change with accommodation.

411
Table 4 shows the mean changes across subjects of DIA, VOL and eigenlenses coefficients (a1, 412 a2, a5 and a6), between the two accommodative demand states, and the p-value of a paired t-413 test to compare if the mean changes with accommodation were statistically significant.As 414 expected, asymmetric terms a3, a4 (not shown) and VOL did not statistically change with   Figure 9 shows the geometrical changes with accommodation of the full shape obtained with the eigenlenses method for two different subjects.Blue line represents the non-accommodated state (0 D of accommodation demand); red dashed line shows the accommodated state (4.5 D of accommodation demand).For visualization purposes, the figure shows 2-D profiles calculated as the mean of the meridians that compose the 3-D model (azimuth angles defined in the sampling in spherical coordinates).

S#2 S#11
Fig. 9. Changes of the full shape with accommodation obtained using the eigenlenses method.Blue: non-accommodated state (0 D of accommodation demand); red: accommodated state (4.5 D of accommodation demand).2-D profiles calculated as the mean of the meridians that compose the 3-D model.

Changes of the eigenlenses coefficients and the geometrical parameters with refractive error
Table 5 shows the mean DIA, VOL and eigenlenses coefficients (â 1 , â2 , â5 and â6 ) pooling emmetropic subjects (n=4, mean refractive error= 0.19 ± 0.37 D, spherical error ranging from 0 to 0.75 D) and myopic subjects (refractive error ≤ −1 D, n=13, mean refractive error = -3.36 ± 1.98 D, ranging from -6.75 to -1.25 D) separately.Data are shown for relaxed accommodation (0 D) and for 4.5 D of accommodative demand.We found statistically significant differences between emmetropes and myopes (t-test; p<0.05 significance) in the lens diameter at 4.5 D of accommodative demand and in the coefficients â2 and â5 for both accommodative states (marked by *).

Discussion
In earlier work [33] we presented the novel concepts of eigenlenses and eigencenters for the representation of the full shape and the central area of the crystalline lens, respectively.They were constructed from 3-D OCT images of 133 ex vivo lenses (newborn to 71 y/o), thus capturing compactly and accurately the geometry of lenses of different ages.In this paper, we have shown the application of the eigenlenses and eigencenters method to the full shape estimation of the crystalline lens from in-vivo OCT images.Eigencenter coefficients are obtained from the central part of the lens that is visible through the pupil, and eigenlenses coefficients (and thus, the full shape of the lens) are estimated from them using an estimator previously obtained using the set of ex vivo lenses.The estimation of eigenlenses from eigencenters is made possible thanks to the compaction ability of the representations, that capture the 96% of the variance in the training set with only 6 coefficients and more than 80% with 2 coefficients.This allows to obtain an accurate estimator that only depended on 6 variables, avoiding overfitting problems [33].
The eigenlenses method has several advantages when compared with our previous 2-region parametric model (2RM) [27] for the full shape estimation of the crystalline lens.Unlike in 2RM, where the full shape is obtained by merging a set of parametric surfaces that represent different portions of the lens, eigenlenses always lead to a smooth and realistic shape by construction, because lens deformation patterns (eigenlenses) are added to the average lens (Eq.( 1)), which is inherently smooth.Furthermore, as shown in Table 3, eigenlenses is a more repeatable method than 2RM (e.g., reducing the standard deviation across measurements by a factor of around 1.8 in DIA and 1.4 in VOL), and much more robust: from the total 83 available 3-D models of the central part, 16 full shapes obtained using 2RM were not accurately constructed and thus discarded (i.e., ∼19% led to non-smooth or anatomically impossible shapes, remaining 67 useful models), while no models were discarded due to a wrong construction with eigenlenses.Thus, 2RM required user interaction to discard those models.Finally, to obtain the representation of a new test lens, the eigenlenses method simply finds the projection of the lens on the eigenlenses space and no parameter optimization is required, unlike in the 2RM approach.This reduced by a factor of 2,000 the computational time of a full shape estimation (from 8.1 s to 4.3 ms).High repeatability, robustness, and low computational cost are very desirable properties for the implementation of the method in commercial OCT systems.
Geometrical changes with accommodation can be well described with changes in the value of the eigenlenses coefficients (Table 4).Specifically, accommodation shape changes are driven by coefficient â2 , which increased approximately 1.92 D −1 with accommodation (i.e., lens became more rounded), while â1 slightly decreased 0.6 D −1 , â5 increased 0.4 D −1 , and â6 decreased 0.09 D −1 .There are several studies in the literature where finite element models (FEM) are used to predict the deformation of the crystalline lens in the accommodation process upon application of equatorial forces by the ciliary muscle, assuming some age-dependent lens material properties and specific lens geometries as a function of age (11, 29 and 45 y/o in [41,42]).The combination of the eigenlenses representation (for the description of patient-dependent full shape changes with accommodation) and FEM can be very useful to study the crystalline lens mechanical parameters of a given patient [43], improving the understanding of the interplay between the mechanical, geometrical, and optical factors in the accommodative response.Also, our results on the changes of eigenlenses coefficients with accommodation can be useful in the generation of geometrical models with accommodation, which is challenging with classical descriptions of the crystalline lens using for example radii of curvature and asphericities or other more complex full shape models [44].
Mean DIA obtained with eigenlenses for un-accommodated lenses (8.87 mm) is close to the one reported in [29] based on the 2RM method (8.94 mm) and by Atchison et al. [6] using MRI (9.03 mm); DIA decreased at 4.5 D of accommodative demand to 8.65 mm (eigenlenses method).Atchison et al. reported average 8.71 mm lens diameters (MRI) in young subjects for 4.8-6.9D accommodative demand.For VOL, our mean reported value for unaccommodated lenses was 140 mm 3 , which is slightly lower than our report using 2RM [29] (146 mm 3 ) and MRI-based data (154 mm 3 in [9]).Changes of DIA were statistically significant with accommodation (0.24 mm upon 4.5 D of accommodative demand).Assuming linear changes of DIA with accommodation, our mean rate of change (around -0.05 mm/D) is similar to those reported by Atchison et al. [6] (slopes around 0.048-0.067mm/D) Hermans et al. [7] (0.07 mm/D), and lower than those reported with 2RM [29] (0.13 mm/D) and in Shepard et al. [9] (0.09 mm/D).Nevertheless, note that changes are reported as a function of accommodative demand, and some differences may arise from individual differences in accommodative response.The presence of accommodative lags will induce an underestimation of the rates compared to estimates as a function of the actual accommodative response [29].Changes in VOL with accommodation were not statistically significant, as expected and found in previous literature using OCT [29,45] and MRI [7].
Second and fifth eigenlenses coefficients were significantly different between myopes and emmetropes.Specifically, â2 decreased in myopes (for both, 0 D and 4.5 D states) consistent with more stretched lenses in myopic subjects, and â5 increased (only in the 0 D state).Also, we found that DIA in myopes was statistically significantly higher than in emmetropes (only significant for the 4.5 D state), in good correspondence with our previous report [31].Furthermore, changes of DIA with accommodation (between relaxed and 4.5 D of accommodative demand) were statistically significantly higher in emmetropes (0.34 ± 0.05 mm) than in myopes (0.20 ± 0.09 mm), which could indicate a lower accommodative response in myopes.The ability of the eigenlenses to capture the full shape of the crystalline lens with a small number of coefficients makes them suitable metrics in applications monitoring myopia and accommodation [46], and in studies of the role of the crystalline lens during emmetropization process, myopia development, and myopia treatment [31,47].
In previous work [3], we demonstrated that including some parameters of the full shape of the crystalline lens (e.g., the EPP, LSA and VOL) improved the estimation of the post-operative IOL position in cataract surgery, which is the largest contributor to post-surgical refractive errors [48]).Thanks to their properties and construction, eigenlenses coefficients describe the full shape of the crystalline lens in a more accurate and compact way than other "classical" geometrical parameters (for example, the lens shape is better described with the first two eigenlenses coefficients than with its VOL and EPP).This opens the possibility of improving the IOL position estimation formulas using directly eigenlenses coefficients instead of geometrical parameters as previously shown [3].The same reasoning applies also to accommodative IOLs (A-IOLs) for presbyopia, where prior knowledge of the full shape of the crystalline lens could be critical for the proper sizing and thus the correct operation of the implanted A-IOL [2,49,50].Finally, full shape lens models could be useful for improving IOL design and for customizing accommodating IOLs, by studying how the full shape geometry modulates the optical/mechanical performance of a given IOL or how the capsular bag shape relates with the IOL tilt and decentration.
The gradient of the refractive index (GRIN) of the crystalline lens may affect the posterior lens shape retrieved from OCT images.Nevertheless, we demonstrated a minor influence in the obtained full shape estimation in earlier work [27,29].Eigenlenses have been applied in vivo in this paper, but they are constructed from a set of ex vivo lenses [33], which geometry can be considered fully-accommodated.Cataract and presbyopia applications of eigenlenses correspond to nearly presbyopic or presbyopic eyes, where one expects ex vivo and in vivo lenses to have approximately the same shape, as changes with accommodation will be small.Furthermore, the set of ages of the ex vivo lenses used for the eigenlenses construction ranged from new born to 71 y/o, thus including very different geometries that are learned by the models.Eigenlenses were constructed using lenses from two different populations: 28 lenses from an eyebank in Spain, and 105 lenses from an eyebank in India.In a recent study [51] we did not find population differences in geometry between isolated European and Indian lenses.
In the future, we plan to include in the training set isolated ex vivo young lenses measured in a lens stretcher system that simulates non-accommodated states [52,53] in order to improve the representation of young lenses with accommodation.In addition, coupling of the eigenlenses representation of the lens with Finite Element Models (and applied equatorial forces) on this lens stretcher database, and further applications in vivo can give insights into the mechanism of accommodation and its loss in presbyopia.The full shape coefficients a k are estimated from coefficients c k using multivariate linear regression.In the future, the estimation method could be improved using "classical" machine learning methods for supervised learning, as neural networks or support vector machines (SVMs).Finally, a related interesting future application of the eigenlens framework is the study of the change of eigenlens coefficient/lens parameters as a function of age, and as a function of refractive error, in relation to the accommodative response.

Conclusion
We have shown application of the proposed eigenlenses representation for the estimation of the crystalline lens full shape in vivo from OCT images, where only the central part of the lens visible through the pupil is available.We have demonstrated that the eigenlenses method has higher repeatability, robustness and computational speed compared to the previously proposed 2RM method for the estimation of the full shape of the lens in vivo.Also, we have shown the ability of eigenlenses to capture geometrical changes with accommodation, only requiring a small number (2-4) of coefficients, and analyzed the relationship between eigenlenses coefficients, accommodation and refractive error.Accurate patient-specific modelling of the crystalline lens in vivo with the eigenlenses representation has the potential to improve cataract surgery (more accurate selection of IOL power) and the design of solutions to restore accommodation (accommodative IOLs).

Fig. 1 .
Fig. 1.(a) Representation of the 6 first eigenlenses obtained, for each eigenlens k, as l + a k e k .Three different a k values are illustrated for every eigenlens k (in blue, red and yellow).The average lens l is represented in all subfigures in red (a k = 0) and the different eigenlenses can be interpreted as deformations from this average lens shape.(b) Example of how the representation of a given lens improves as more eigenlenses are included.The black lens shows the lens l n to be represented.Blue, red, yellow and green lenses show the obtained representation as the average lens l, the a 1 = −25, a 2 = −15 and a 5 = −5 coefficients are included, respectively.Asymmetric terms (a 3 and a 4 ) are not included in the example for clarity.(c) Representation of the 6 first eigencenters obtained, for each eigencenter k, as z + c k g k .Three different c k values are illustrated for every eigencenter k (in blue, red and yellow).The average central area z is represented in all subfigures in red (c k = 0) and the different eigencenters can be interpreted as deformations from this average central area.

Fig. 2 .
Fig. 2. Training a model for the estimation of the full shape from the central area.From the training set of ex vivo lenses, e k , l and g k , z are firstly obtained and stored.Then, for each ex vivo training lens, a k (representation of the full shape) and c k (central part) are calculated.From them, a matrix β is obtained using multivariate linear regression, that allows the estimation of the a k from c k in a new unseen lens (and thus, the estimation of its full shape from its central part measured in-vivo).

Fig. 3 .
Fig. 3. 3-D model construction within the pupil, including the automatic surfaces detection in each B-scan and the distortion correction.

Fig. 4 .
Fig. 4. Full shape estimation from the 3-D model within the pupil, including: (1) Change of coordinates: 3-D model within the pupil change of coordinates, from cartesian to spherical; (2) Projection of the residual data in the eigencenters basis to obtain c k coefficients; (3) Prediction of âk coefficients from the c k coefficients; (4) Reconstruction of the full shape from the âk coefficients.

Figure 5
Figure 5 shows an example of a central B-scan acquisition and the 3-D full shape model obtained after applying the proposed method.

4 Fig. 6 .
Fig. 6.Correlation between eigenlenses coefficients and full-shape geometrical parameters.First column: correlation between â1 and DIA, VOL, LSA and EPP.Second column: correlation between â2 and DIA, VOL, LSA and EPP.Third column: correlation between â5 and DIA, VOL, LSA and EPP.

Fig. 7 .
Fig. 7. Bland-Altman plots for the different geometrical parameters; Y axis: difference of the

Fig. 7 .
Fig. 7. Bland-Altman plots for the different geometrical parameters; Y axis: difference of the methods (2RM-proposed); LoA: Limit of Agreement (1.96 • SD, where SD is the standard deviation of the difference).

Figure 8 shows
Figure 8 shows DIA and "symmetric" eigenlenses coefficients (â 1 , â2 , â5 and â6 ) for the 17 subjects in two accommodative demands: 0 D (blue circles) and 4.5 D (red crosses).Table4shows the mean changes across subjects of DIA, VOL and eigenlenses coefficients (â 1 , â2 , â5 and â6 ), between the two accommodative demand states, and the p-value of a paired t-test to compare if the mean changes with accommodation were statistically significant.As expected, asymmetric terms â3 , â4 (not shown) and VOL did not statistically change with accommodation. 404

Figure 9 Fig. 8 .
Figure9shows the geometrical changes with accommodation of the full shape obtained with

Table 4 .
Mean changes of DIA, VOL and eigenlenses coefficients ( â1 , â2 , â5 and â6 ) between two accommodative states (0 and 4.5 D) values correspond to a paired t-test; * indicates statistically significant changes at a significance level of p<0.05.

Table 3 . Repeatability (standard deviation) in the estimation of each geometrical 395 parameter for the 2RM and the eigenlenses methods.
CPU, 8 GB RAM, using Matlab software.