Eigenlenses : a new model for full crystalline lens shape representation and its applications

: The crystalline lens is an important optical element in the eye, responsible for focusing, and which experiences signiﬁcant changes throughout life. The shape of the lens is usually studied only in the optical area (central 4 to 6 mm). However, for a great number of applications, a description of the full shape of the crystalline lens is required. We propose a new method for the representation of the full shape of the crystalline lens, constructed from 3-dimensional optical coherence tomography images of 133 isolated crystalline lenses (0-71 y/o), which we have called eigenlenses . The method is shown to be compact and accurate to describe not only the full shape of the crystalline lens, but also the optical zone in comparison with other methods. We also demonstrate its application to the extrapolation of the full shape of the crystalline lens from in-vivo optical images of the anterior segment of the eye, where only the central part of the lens visible through the pupil is available, and in the generation (synthesis) of realistic full lenses of a given age. The method has critical applications, among others, in improving and evaluating myopia and presbyopia treatments.


Introduction
The crystalline lens of the eye is a fascinating optical structure that, along with the cornea, projects the images of the outside world onto the retina. The crystalline lens is transparent, although its transparency is lost with age, in what is known as a cataract, requiring its replacement by an artificial lens. The crystalline lens is coupled with the cornea to minimize the aberrations in central vision, although this balance is lost with age. Additionally, the crystalline lens shows a gradient index structure, this profile being also age-dependent, and playing an optical role. Very remarkably, the crystalline lens changes its shape to focus near and far objects in response to an accommodative demand, by contraction of the ciliary muscle and action of the zonulae which connect the equatorial region of the crystalline lens with the ciliary body [1].
The cornea is the external part of the eye and therefore its shape can be easily quantified. Nevertheless, the crystalline lens lies behind the cornea and the iris, and its quantification becomes more complicated. Optical techniques, including Purkinje imaging [2,3], Scheimpflug imaging [4][5][6][7][8] and Optical Coherence Tomography (OCT) [9][10][11][12][13][14][15][16], have been used to assess the crystalline lens, although special care needs to be taken to correct the refraction from preceding surfaces to allow an accurate reconstruction of the lens surfaces geometry. In any case, in vivo optical techniques are limited by the iris, which invariably blocks the incident light preventing the imaging of the peripheral zone of the crystalline lens. Non-optical techniques such as Magnetic Resonance Imaging (MRI) [17][18][19][20] or Ultrasound Biomicroscopy (UBM) [21] allow visualization of the full shape of the lens including the peripheral region, but at a much lower resolution and longer acquisition times [17][18][19][20][21], or limited by the need of contact, the need of a very experienced examiner, and other types of distortion (UBM) [21].
There are many applications in which an accurate and compact (with a small number of parameters) description of the crystalline lens full shape is required. We showed recently that the ability to compactly describe the full shape of the crystalline lens is important to improve the accuracy in the selection of an intraocular lens (IOL) that replaces the crystalline lens after a cataract surgery [22]. Also, knowledge of the full shape of the crystalline lens turns critical for the success of methods to restore accommodation by modifying the elasticity of the lens material, refilling the lens capsule with a deformable polymer [23][24][25], or for proper sizing of accommodative introacular lenses [26]. In addition, an accurate description of the full geometry at different ages will improve the biomechanical models of the crystalline lens, which has direct benefit to the computational modelling of the accommodative process and to the design of customized accommodating IOLs [27]. Knowledge of the changes undergone by the full crystalline lens in early ages will give insights on the emmetropization process and the role of the crystalline lens in the development of refractive errors [15,28]. Finally, with emerging pharmacological treatments of myopia and presbyopia, in many ways targeting directly the crystalline lens, it has never been as relevant as now to accurately quantify its full dimension to monitor its changes and to assess the treatment.
Typically, studies describing the crystalline lens geometry use radius of curvature and asphericity (R, Q of the anterior and posterior surfaces) [29] or a Zernike polynomial fit of the surface elevation [12,14], but these models are not aimed at representing the periphery of the lens. The full shape geometry of the lens has been previously modeled using hyperbolic cosine functions [30], generalized or interdependent conic functions [31,32], and polynomial functions [32,33]. However, some of these models are purely theoretical, and not validated against real crystalline lenses [30,31]. While there have been attempts to study the linear dependency of the parameters of the lens model with age, no parameter appeared to correlate significantly with age (one curve approach [33]), or the correlation was weak, and only for a low number of parameters (two curves approach, [33]). This lack of correlation and compactness (full shape models have a high number of parameters, ranging from 8 [30] to 20 [33]) make difficult their use for the interpretation of the changes of the geometry with age, the generation of realistic lenses, or the estimation of geometrical and optical properties from the coefficients of the models. Finally, all the previous full shape models assume rotational symmetry around the optical axis.
In this manuscript we propose a new method for the representation of the full shape of the crystalline lens that we call eigenlenses. Their construction is inspired by the concepts of eigenfaces [34] and active shape models [35]. The proposed method allows representation of any lens as the addition of an average lens (the mean shape over a training set of lenses) and a linear combination of typical deformations from this mean shape (eigenlenses). Eigenlenses lead to a 3-Dimensional (3-D) accurate and compact representation of the full geometry of the crystalline lens and have an easy interpretation (each eigenlens is related with an intuitive change in the lens geometry, as for example their general size or their shape factor). We compare the accuracy and compactness of eigenlenses with other representation methods. We also present applications in the full shape generation of realistic crystalline lenses of any age, and in the estimation of the crystalline lens full shape from images collected in-vivo with optical techniques, in which only the data within the pupilary region are available.

Tissue preparation
We studied 133 isolated crystalline lenses from 112 human donors. Two eye banks were used for this study. From the first eye bank (Banc de Sang i Teixits, Barcelona, BST, Spain) we obtained 28 eye globes from 24 donors (age range 19-71, mean±STD 48±13 y/o) [36] and from the second eye bank (Ramayamma International Eye Bank at L V Prasad Eye Institute, LVPEI, Hyderabad, India) we obtained 105 eye globes from 88 donors (age range 0-56, 26±14 y/o) [37]. After enucleation, a surgeon carefully isolated the crystalline lens from the eye globe and immediately placed it on a custom made lens holder of nylon sutures within a cuvette filled with preservation media (DMEM/F-12 HEPES no phenol red, GIBCO; and BSS, Alcon Laboratories was used in the first and second eye banks respectively). The holder avoids the contact between the lens surfaces and the bottom of the cuvette. The median post mortem time across lenses was 34 hours.
Note that lenses with capsular separation or with cataract were discarded (not measured). The 133 crystalline lenses of the study were selected from a larger group of 157 measured lenses of which lenses with visible deformations or low image quality, which prevented obtaining the full lens shape, were excluded. Also, to prevent the possibility of water uptake, we followed the suggested protocols by Augusteyn et al. [38]. More details on tissue preparation and experimental protocols followed can be found in our previous work [37].
This research followed the tenets of the Declaration of Helsinki and was approved by the Institutional Review Boards of CSIC, BST, and LVPEI.

Optical coherence tomography imaging systems
The lenses from the first eye bank were measured with a custom developed spectral domain OCT (SD-OCT) system [39,40] composed by a superluminescent diode with a central wavelength of 840 nm and a full width half maximum (FWHM) bandwidth of 50 nm. The axial range was 7 mm in air, resulting in an axial pixel size of 3.4 µm with a measured optical axial resolution of 9.2 µm in air. The acquisition speed was 25000 A-scans/s and each 3-D volume was composed of 1668 A-scans/B-scan and 60 B-scans on a 12x12 mm lateral area.
The lenses from the second eye bank were measured with a different SD-OCT system described before [41]. The core of the instrument is a commercial imaging system (ENVISU R4400, Bioptigen Inc.) equipped with a superluminescent diode with a central wavelength of 880 nm and a FWHM of 40 nm. The axial range was 15.18 mm in air, resulting in a pixel size of 7.4 µm and the calculated optical axial resolution was 8.5 µm in air. The acquisition speed was 32000 A-scans/s and each 3-D volume was composed of 600 A-scans/B-scan, and 100 B-scans on a 15x15 mm lateral area.
The crystalline lenses were aligned with the OCT systems to collect images of the full shape of the lens. The lenses were first imaged with the anterior surface facing the OCT beam ("anterior-up" measurement), and then flipped over and imaged again with the posterior surface facing the OCT beam ("posterior-up" measurement). Figure 1 shows an example of the central B-scan from "anterior-up" (Fig. 1(a)) and "posterior-up" (Fig. 1(b)) measurements.

Obtaining the 3-D shape of the lenses
The main steps for obtaining the 3-D shape of the lenses from OCT raw B-scans were: segmentation of the images, distortion correction, tilt removal, and registration.
The surface of the full lens was automatically segmented in each B-scan using thresholding, Canny detectors, morphological operations, and a-priori knowledge of the measurements [10,14]. Then, 3-D data composed of all B-scans were fit with Zernike polynomials up to the 4-th order, and the resulting smooth surface was used to refine the segmentation in each B-scan. This process was repeated applying three different orientations to the images, which allowed improvement of the detection of the central part and the periphery of the lens [37].
Although the OCT systems produced a partially-telecentric scan, the segmented surface elevation was corrected from the residual "fan distortion" arising from the scanning architecture and the optics of the system [39,42]. Vertical distances in the OCT images represent the optical path of the beam. Thus, the 3-D shape of the lenses was obtained by merging the 3-D anterior lens surface segmented from the "anterior-up" measurement and the 3-D posterior lens surface segmented from the "posterior-up" measurement, because these surfaces are not affected by refraction in previous optical surfaces or in the lens gradient index [43]. The distortion due to the presence of the preservation media was corrected by dividing the sag of the surfaces by the group refractive index of the media, n=1.345 for the DMEM (at 840 nm) and BSS (at 880 nm) [44].
The tilt of the distortion-corrected anterior and posterior 3-D surfaces was removed, and both surfaces were registered laterally (in the X-Y plane) using the equator of the lens [9]. Then, to account for possible rotations when flipping the lens over in the measurements, posterior surface was rotated until the equator overlapping between anterior and posterior surfaces was maximized. Axial registration (along the Z axis) was performed by matching the central thickness of the lens (LT), calculated using the deformation in the image of the cuvette [45]. Figure 2 illustrates the process for obtaining the 3-D shape of the lenses.

Eigenlenses construction
From the 3-D shape of the set of N = 133 lenses, eigenlenses were constructed following the next steps: 1. The origin of coordinates was set in the lens center, defined laterally as the (X, Y) center of the lens equator and axially as Z = LT/2.
2. The segmented points of the 3-D shape of the lenses were transformed to spherical coordinates, i.e., (X, Y, Z) −→ (r, θ, ϕ), where r is the radial distance from the origin of coordinates to the surface, and θ and ϕ are the elevation and azimuth angles respectively, as defined in Fig. 3(a).
3. In order to have all the lenses sampled at the same (θ, ϕ) angles, we proceeded as follows: • Azimuth ϕ was sampled at P = 100 uniformly spaced angles in the interval [−π, π], and elevation θ was sampled at Q = 100 uniformly spaced angles in the interval Note that, although this sampling does not result in evenly placed points on the surface of a sphere, we observed that eigenlenses resulting from different sufficiently dense samplings were very similar, and we chose this approach for simplicity.
• For each lens n, the distances from the origin of coordinates to the surface at every pair θ i , ϕ j , r θ i ,ϕ j , were obtained by cubic interpolation from the segmented points of the corresponding n-th 3-D shape. This led to a column vector with P × Q elements describing the shape of the lens n: where T denotes transpose. 4. The average lens (l) was obtained as the mean of all lenses, and the residual data for each lens (∆ n ) was calculated by subtracting the average lens from the lens datā Figure 3(a) illustrates the processes defined above.

5.
We call eigenlenses to the eigenvectors of the covariance matrix of the residual data (i.e., their principal components), that were obtained solving the classical diagonalization problem [46] is the sample covariance matrix of the residual data, e k is the k eigenlens and λ k its corresponding eigenvalue. With this formalism, any crystalline lens can be represented with a certain level of accuracy as the average lens plus a linear combination of eigenlenses: where K is the number of eigenlenses used in the representation and a k is the multiplicative coefficient of eigenlens k. This implies that a given lens can be defined -with a certain level of accuracy-by a set of coefficients {a k , k = 1, . . . , K}.

Eigenlenses for the representation of the optical zone
Eigenlenses were additionally constructed to represent the optical central area of the crystalline lens. We followed the process explained in the previous section but, in this case, defining the lenses only in their central part. To do that, we changed the sampling of the elevation angle θ so that, instead of covering all the 3-D space and thus the full shape of the lens ([−π/2, π/2]), it just covered an angle corresponding to the central area of the lens (e.g., [−π/2, −π/4] ∪ [π/4, π/2] for covering approximately the central area of 4 mm of diameter). Figure 3(b) illustrates the θ sampling process.

Data analysis
To evaluate the accuracy of the model and its generalization ability (the capability of representing an unseen lens) 10-fold cross validation, i.e., training the model with 90% of the lenses and testing it on the remaining 10%, shifting the test set in each fold, was used. The 10-fold cross validation was repeated 100 times, and the root mean squared error (RMSE) was estimated averaging the error (the difference between the actual lens and its representation with K eigenlenses) in the test lenses. Note that the RMSE is calculated in spherical coordinates (approximately orthogonal to the surfaces) to obtain the eigenlenses performance for describing the full shape.
To assess if the RMSE produced by state-of-the-art methods were significantly different from the proposed, a one-way ANOVA with Bonferroni's correction for multiple comparisons test was applied. To compare if the proposed method was significantly different to other representations when using the same number of parameters, paired t-tests were used.
Multiple linear regression analysis was performed to study the prediction of the full shape from the central part of the lens, and the adjusted coefficient of determination R 2 and the p-value for testing the hypothesis of no correlation were obtained.
For all analyses, statistical significance was defined as a p-value lower than 0.05. Calculations were obtained using Matlab software (MathWorks, Natick, MA, USA).

Eigenlenses interpretation and properties
Eigenlenses with high associated eigenvalues represent the main way in which the points of the surface tend to move together (i.e., how the shape varies) across lenses with respect to the average lens, and thus can be thought as deformation patterns with respect to the average lens. The eigenvalue λ k in Eq. (4) represents the amount of variance explained by the corresponding k eigenlens (e k ). Thus, eigenlenses with larger λ k capture the most relevant modes of variation across lenses. In addition, the residual variations represented by the eigenlenses are mutually orthogonal, leading to an efficient representation. Figure 4 shows the shape changes associated to the 6 eigenlenses with the largest λ k . To illustrate these changes, each eigenlens was multiplied by 3 different coefficients a k = −3σ k , 0, 3σ k , where σ k is the standard deviation of the coefficients over all the lenses for the eigenlens k, and added to the average lensl: Visually, the first eigenlens (the most common deformation) is related with changes in the size of the lens; the second eigenlens is related with changes in the aspect ratio (i.e. lenses more rounded or "stretched"); the third and fourth with asymmetric changes in X or Y directions; and the fifth and sixth with fine changes in the shape of the surfaces (related with the asphericity of conicoids [47] or with rotationally symmetric Zernike polynomials). As a consequence of eigenlenses construction and properties, any lens can be represented with high accuracy with a small number of eigenlenses (the ones with largest eigenvalues) plus the average lens. Furthermore, eigenlenses coefficients are not correlated by construction and thus any lens created from a set of coefficients {a k } within the range of values of the training set is realistic.
A graphical user interface (GUI) was developed (available upon request) to illustrate the changes that the first 5 eigenlenses represent and their influence of on the final generated lens (see Visualization 1).

Choosing the number of eigenlenses
The number of eigenlenses K used in the representation should be chosen as a trade-off between obtaining a compact representation (with a small number of eigenlenses) and capturing the relevant geometrical information of the lenses to achieve an accurate representation.
We proposed two metrics for choosing K: percentage of variance (PVar) explained by the set of the first K eigenlenses; and root mean squared error (RMSE), obtained averaging across test lenses (from 100 repetitions of 10-fold cross validation) the error between the actual lens and its representation with K eigenlenses. Standard deviation (STD) of the RMSE (across lenses and folds) and of the PVar (across folds) was also calculated. Note that if PVar=100 or RMSE=0 we can represent all the test lenses without error. Figure 5(a) shows PVar±STD and RMS±STD as a function of K. We chose K = 6 as the optimal value for our purposes, as this was a point where increasing K did not decrease substantially the RMSE (or increased the PVar). We will use these value of K = 6 throughout the rest of the paper.

Model performance for describing the full shape: comparison with other representations
RMSE of the full shape description averaged across the test lenses was compared with the RMSE of three other representations in the state of the art, namely: • Lenses reconstructed from the best sphere fitting of anterior and posterior surfaces of the lens, lens thickness, and posterior surface position defined as the (X,Y) coordinate of its apex (5 parameters).
• Lenses reconstructed from the best conicoid fittings, that comprised the same parameters than sphere fittings plus asphericity values of anterior and posterior lenses (7 parameters) [29,47].
• Zernike approximation of the surfaces, using 6, 15 and 28 terms for representing anterior and posterior surfaces (12, 30 and 56 terms in total respectively) [12,14]. Figure 5(b) shows the RMSE for the representation using the first 6 eigenlenses and for the other three representations. The error with the proposed method was significantly lower (p<0.05, one-way ANOVA with Bonferroni's correction) than with sphere, conicoid, 12 Zernike coefficients (Z12) and 30 Zernike coefficients (Z30) representations (error 2.23, 2.20, 2.39 and 1.51 times lower respectively). This was achieved using similar number of parameters (6) than sphere (5) and conicoid (7), and lower number of parameters than the number of Zernike coefficients used by Z12 and Z30. Only the Z56 representation produced similar RMSE as with the eigenlenses method, but with much higher number of parameters (56 vs 6).

Model performance for describing the optical zone: comparison with other representations
Accuracy in the representation of the optical zone with eigenlenses, constructed as described in Section 2.5, was compared with the accuracy obtained with sphere, conicoid, and Zernike representations described in the previous section. Figure 6 shows the RMSE obtained with the different methods calculated in the central 6 mm. The error using 6 eigenlenses was significantly lower (p<0.05, one-way ANOVA with Bonferroni's correction) than with sphere and conicoid representations. Note that using 12 eigenlenses the mean RMSE (RMSE=0.0085 ± 0.0034 mm) was similar to the one obtained with Z30, and using 17 eigenlenses (RMSE=0.0069 ± 0.0026 mm, shown in the figure with the label Eiglen 17) was similar to the one of Z56, but with lower STD in both cases. Note also that the RMSE was significantly lower (p<0.05, paired samples t-test) with eigenlenses than with Zernike when using the same number of parameters in both representations, independently of the number of parameters chosen (e.g., the RMSE was significantly lower with 30 eigenlenses than with 30 Zernike terms).

Application example 1: lens generation (synthesis)
We used the eigenlenses representation to generate realistic lenses of a given age. Firstly, we estimated from the data the underlying conditional probability distribution of the coefficients given an age: P(a 1 , . . . , a K |age = A), with K = 6 the number of eigenlenses used. We assumed that P(a 1 , . . . , a K |age = A) was distributed as a multivariate normal distribution, and we estimated its mean vector and covariance matrix. To avoid estimating these parameters using only the lenses with a specific age, we used all the data (i.e., all ages) weighting every sample i as a function of the square distance between its age (age i ) and the age of interest (A) using a Gaussian kernel and normalizingŵ whereŵ i is the weight assigned to each specific sample, fulfilling that N i=1ŵ i = 1; and the parameter W, which controls the width of the Gaussian kernel, was set to 5. From the weighted data, we estimated the weighted covariance matrix and mean vector [48]. Note that the weighting process gives more importance in the estimation to samples of ages closer to the age of interest, and the normalization of the weights is required for the estimation [48].
Once the probability distribution function had been estimated, lenses of a given age were generated sampling from P(a 1 , . . . , a K |age = A), obtaining for example typical lenses (mean value), "atypical" lenses that were at a given distance to the mean, or random lenses. Figure 7 illustrates the process. For visualization purposes, in this example we used only 2 eigenlenses. Figure 7(a) shows P(a 1 , a 2 |age = A) for ages A = 5, A = 30 and A = 60 y/o. From each distribution, we obtained a random pair (a 1 , a 2 ), and we generated 3 lenses of 0, 30 and 60 y/o as l =l + k=2 k=1 a k e k (Fig. 7(b)). A GUI showing the lens generation process is available upon request (see Visualization 2).

Application example 2: lens extrapolation from pupil information
When imaging the crystalline lens in-vivo with optical technologies, the iris blocks the incident light and only the information within the pupil (the central part of the lens) is available (see Visualization 3). Eigenlenses were used to estimate the full shape of the crystalline lens from the information available in-vivo, following the next steps: 1. In-vivo conditions were simulated in the training set of isolated lenses, i.e., we assumed that only the central part of the lens was visible. Experiments were done assuming a pupil diameter of approximately 4 and 5 millimeters, sampling the elevation angle in polar coordinates at the interval [−π/2, −π/4] ∪ [π/4, π/2] and at [−π/2, −π/6] ∪ [π/6, π/2] respectively.

2.
Eigenlenses for the representation of the optical zone were obtained as described in Section 2.5. We call the corresponding coefficients {c k , k = 1, . . . , 6}.
3. A parametric expression was found to estimate each coefficient a k , k = 1, . . . , 6, that describe the lens full shape, from a set of coefficients of the central area (i.e., a 1 = f (c 1 , . . . , c 6 ), a 2 = f (c 1 , . . . , c 6 ),. . . , a k = f (c 1 , . . . , c 6 )) with multiple linear regression using least squares. Figure 8 illustrates the process. For each original full lens shape of the training set ( Fig. 8(a), blue lens on the top) its representation was obtained (coefficients {a k , k = 1, . . . , 6}, Fig. 8(a), black lens on the bottom-left). In the same lens, in-vivo conditions were simulated and the c k coefficients representing the central area were calculated ( Fig. 8(a), blue area on the bottom-right). A model for the prediction of the a k from the c k was obtained from all the training lenses. Then, this model was applied to the test lenses ( Fig. 8(b)), for which only the information of the central part (and thus the corresponding c k ) was assumed to be known. Note that the model could also be applied in the same manner to in-vivo measurements.
The goodness of fit of the coefficients prediction was evaluated with the adjusted coefficient of determination R 2 . Table 1 shows the R 2 and p-values for the prediction of the a k coefficients for 4 mm and 5 mm pupils. Note that the model explained a high proportion of the variance when the first five a k coefficients were predicted, but the performance decreased substantially in the prediction of a 6 .
The full shape reconstruction was evaluated in the test lenses, obtaining an average RMSE between the original full lens and the reconstructed from the central area of RMSE=0.074 ± 0.025 mm for a 4 mm pupil and RMSE=0.066 ± 0.020 mm for a 5 mm pupil. A GUI showing the lens extrapolation process is available upon request (see Visualization 4).

Discussion
We proposed a method for obtaining a compact and accurate representation of the full shape of the isolated human crystalline lens constructed from 3-D OCT images in 133 lenses from 112 human donors in a range of ages spanning from new born to 71 y/o. Then, we used the method in two relevant applications: (1) generation of realistic lenses of a given age, and (2) estimation of the full shape of a crystalline lens from the central part of the lens that is visible in in-vivo measurements when using optical techniques. The proposed method represents a change of paradigm in lens representation, usually focused on closed-form mathematical expressions where several parameters are optimized so that the model fits some data. Differently, eigenlenses are constructed on the basis of a training set of lenses, capturing the main shape variations across them. Thus, eigenlenses are very efficient in the representation of test lenses ("new" lenses not previously seen by the model) if their shape is captured in the training set. In other words, an important factor to obtain a general and efficient representation is the availability of a training set that captures the crystalline lens shape variability that is present in nature. This variability occurs primarily because of changes with age, but also because of inter-subject variability for a given age. Our data set of 133 lenses of ages ranging from 0 to 71 y/o captures well this variability. Our model shows a strong generalization ability, with test and train error very close to each other (RMSE=0.045 ± 0.011 mm and RMSE=0.042 ± 0.010 mm respectively).
Thanks to the orthogonality of the eigenlenses, the model is very compact and is able to represent accurately the full shape of crsytalline lenses with a small number of parameters (with 2 eigenlenses the model explains more that the 80 % of the variance of the data). Using 6 eigenlenses, the RMSE was significantly lower than using classical representations such as sphere fitting, conicoids, or Zernike polynomials with 30 coefficients. As these methods are not aimed at the description of the periphery of the lens, eigenlenses were also assessed in the optical zone, obtaining an error significantly lower than all other representations when comparing using the same number of parameters.
Urs et al. [33] represented the full shape of the crystalline lens with 2 different models: one curve approach (OCA), that describes the full shape using only one equation, and two curves approach (TCA). They obtained a RMSE in the order of the one obtained with our method but using a larger number of parameters (10 with the OCA and 20 with the TCA). More importantly, Urs et al. did not find correlation of the model parameters with age for the OCA, or the correlation was weak and only in 3 of the 20 parameters in the TCA (R 2 = 0.51, 0.31 and 0.17), preventing the model to be used properly in applications that entail relating its parameters with lenses of different ages. Other interesting properties of eigenlenses are that (i) they are highly intuitive as shown in Fig. 4, allowing to interpret lens geometrical properties as a whole, and that (ii), to obtain the representation of a new test lens, eigenlenses method simply finds the projection of the lens on the eigenlenses space, and, unlike state of the art parametric approaches, no parameter optimization is required.
Sources of error in the 3-D models construction and thus in the lens geometry may arise from different factors as the preservation media used (DMEM or BSS), tissue water uptake [38], capsule detachments or lens deformations, optical or fan distortions, or the segmentation and registration processes. We treated carefully all those factors, as explained in previous sections. Furthermore, eigenlenses are useful independently of small errors in the measured lens geometry. The representation may slightly lose in efficiency, but keep all its characteristic properties intact.
The age distribution of the lenses is not perfectly uniform (see the age distribution for the Indian and Spanish lenses in [37] and [36] respectively). We did not find statistically significant differences in the RMSE (using the 6 eigenlenses representation) between three different age groups (younger ages group, age<20 y/o, n=37; mid ages group, 20 y/o≤age<40 y/o, n=60; and older ages group, age≥40 y/o, n =36; p<0.05, ANOVA with Bonferroni's correction). Therefore we can conclude that the model is able to represent crystalline lenses of different ages with similar accuracy. We did not find a statistically significant difference with gender or ethnicity of lens thickness or lens volume in our data set (20 % of females, 80% of males; 79 % of Indian, 21 % of Spanish lenses). In any case, the scope of the current study is the presentation of a general model to describe the lens, which we have shown it is capable to capture and describe individual differences across lenses, such as those related to age. Future interesting work may include the investigation of the dependence of eigenlenses coefficients on other variables such as gender or ethnicity. Although we have applied the eigenlenses method to OCT images, it is also applicable to crystalline lens images obtained with other imaging modalities, such as MRI or UBM. Eigenlenses could be implemented in current commercial 2-D devices (assuming rotationally symmetric models), and may be relevant in clinic, for example, for IOL selection or complete personalized anterior segment geometry.
Eigenlenses can be further used in other areas. In a previous study [9] the authors presented a method for estimating the full shape of the crystalline lens from OCT in-vivo measurements that, although it was accurate, it required for human supervision to remove badly constructed models. In the proposed method the lens is estimated as a linear combination of deformations plus the average lens, always leading to a smooth and realistic model, and thus being more reliable and robust, specially in non favorable scenarios (e.g., a very tilted crystalline lens measurement). Note that the RMSE obtained between the actual full shape of a lens and its shape estimated from its central optical zone (RMSE=0.074 ± 0.025 mm and RMSE=0.066 ± 0.020 mm for 4 and 5 mm pupils respectively) is in the order of the RMSE between the actual full shape and its representation using 6 eigenlenses (Fig. 5, RMSE=0.045 ± 0.011 mm), indicating that the estimation of the full shape from the optical zone is accurate. The practical procedure to apply the lens extrapolation method to in-vivo OCT images will incorporate the surfaces segmentation and fan and optical distortion corrections (Fig. 2) to obtain the 3-D shape of the lens visible through the pupil. Note that, in-vivo, the optical distortion produced by the refraction of the rays in the preceding optical surfaces has to be corrected to obtain accurate anterior and posterior lens shapes [42]. Once the actual geometry of the central part of the lens is obtained, the corresponding c k coefficients are calculated and from them the a k are estimated (and thus the full shape) using the prediction model obtained in the training phase (Fig. 8). Estimations of the posterior lens shape geometry may be affected by the presence of a gradient refractive index, although we demonstrated a minor influence in these applications in earlier work [9,10]. Further discussions on the the application of ex-vivo based lens models to in-vivo images can be found in [9,10].
In [22] we presented a method for estimating the IOL position from some quantified parameters of the full crystalline lens, specifically the volume and the equatorial plane position. Using the eigenlenses coefficients may improve the estimation of the post-operative IOL position because a more accurate description of the pre-operative crystalline lens shape is allowed with the same number of parameters. For example, the lens shape is better described with the eigenlenses 1 and 2 than with only its volume and equatorial position.
Another suggestive area of application is the combination of the proposed crystalline lens representation with finite element models (FEM). FEM are used to predict the deformation of the crystalline lens with accommodation upon application of equatorial forces by the ciliary muscle, with age-dependent lens material properties. Generally these models rely on specific lens geometries at three ages (11, 29 and 45 y/o) [49][50][51]. Undoubtedly, the incorporation of realistic full shape geometries at different ages will improve the model and thus the understanding of the geometrical, mechanical, and optical factors in the accommodative response. Therefore, it is useful to be able to generate realistic full lenses of a given age, which is difficult with classical descriptions of the crystalline lens (radii of curvature, asphericities, Zernike coefficients of the surface elevation) or with state of the art full shape models [33]. Also, lens generation can improve the IOL design, by studying how the full shape geometry modulates the optical/mechanical performance of a given IOL. For example, a key step in the design process of an IOL platform is the evaluation of the mechanical stability to forces consistent with tensions produced on the lens haptics by the capsular bag (of different diameters). It is conceivable to simulate this test (and the effect on IOL tilt and decentration) using realistic full lens shape models, such as those provided by the eigenlenses representation.
Finally, the ability of the eigenlenses to capture the main geometrical information of the crystalline lens with a small number of uncorrelated parameters makes them optimal for improving the in-vivo estimation of both the optical and biometric changes with a single clinical instrument. This has applications in the monitoring of refractive error developments, assessment of accommodation, pre-operative testing of cataract patients or evaluation of presbyopic corrections, among others [52]. Also, eigenlenses will allow to better understand the geometric changes of the lens as an unique object than using different parameters for quantifying the anterior and posterior surfaces and the thickness of the lens. This will be very useful to study the role of the crystalline lens during the emmetropization process, myopia development, and myopia treatment.

Conclusion
We have proposed the eigenlenses, a novel crystalline lens full shape representation method that is accurate, compact and intuitive. These properties allowed us to apply it successfully for full shape estimation from in-vivo simulated optical images, where only the central part of the lens is visible, and for realistic lenses generation of a given age. Further work will focus on extending the full shape model to different accommodative states using stretched ex-vivo lenses for the training process, and its application in refraction estimation from in-vivo measurements. Also, further investigations may include coupling eigenlenses with finite element models to describe crystalline lens mechanical properties, and with ray tracing analysis to describe its optical properties.