Abstract
We introduce a new type of lens with two gradients of refractive index (GRIN) and of curvature (GRCU) of iso-indicial surfaces, i.e., GRINCU. The inner structure of the lens resembles that of an onion. Each layer is a meniscus lens with infinitesimal thickness, which coincides with an iso-indicial surface characterized by a conicoid shape and a constant refractive index. The internal distribution automatically adapts to the external geometry. Here, we consider the simplest case of a constant gradient of the curvature radius –G, which indicates a linear decrease as we move along the optical axis. The formulation of this type of lens is presented, including its generalization to nonrotationally symmetric conicoid surfaces. The formulation is then applied to model the crystalline lens; the code corresponding to the numerical computation of the 3D refractive index distribution as well as its gradient is provided as a supplementary file. Finally, we confirmed a refractive power increase of nearly 14% when G changes from 0 to 3.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Gradient index (GRIN) lenses are made of inhomogeneous materials, and inside the lens, light propagates through inhomogeneous media [1,2]. Inhomogeneous media are common both in nature and in optical technology. Two paradigmatic examples are air, which always presents some level of inhomogeneity associated with atmospheric turbulence, and the human crystalline lens, which is a classic example of a GRIN lens. They eyes of most animal species contain some type of GRIN optics. One of the first GRIN lens designs was Maxwell’s fisheye dating from 1854 [1]. Advances in theory, design of index distributions, materials, fabrication and applications are important, and GRIN lenses are currently widely used.
Among the wide variety of functions or profiles of refractive index proposed in the literature, we are interested in an onion-type lens. In the continuous case, this lens is made of an infinite number of layers of differential thickness; each layer acts as a thin (infinitesimal) lens. For the differential onion layer to act as a standard optical surface, we want the refractive index n to be constant within that layer. Such layers are iso-indicial surfaces (IISs) with index n. This surface separates the anterior and posterior media with a difference of the refractive index. If the layer, or IIS, is spherical with radius R, then we can compute its optical power as $P = {{d n} / R} = Cd n$, where $C = {1 / R}$ is the curvature and $d n$ is the differential refractive index. This expression from first-order (paraxial) optics suggests that the gradient index (proportional to $d n$) and the curvature C of the IIS are equally important and have a mutual multiplicative effect. The idea of having optical media with two gradients of index and curvature (GRIN/GRCU or GRINCU) was explored in a previous work [3]; the results obtained for the first-order optics approximation confirmed a strong multiplicative effect on the power of the GRINCU lens. Of course, most types of GRIN distributions have an implicit gradient of curvature (the GRCU can be analyzed by computing the curvatures of the IISs), but to the best of our knowledge, this is the first time that the optical effect of the GRCU has been explicitly modeled and studied [3].
In this context, there are two goals of the present work. In the first part, we develop the formulation of the index distribution of a generic GRIN/GRCU lens composed of differential onion layers. To guarantee analytical, treatable expressions, we restrict ourselves to the particular case when the IISs are revolution conicoids, and the gradient of their apical curvature radius is a constant parameter –G [3]. To this end, we generalize and upgrade a formulation developed before to model the crystalline lens [4].
Our second goal is to apply that generic formulation to model the crystalline lens. The model is adaptive [4,5], geometry-invariant [6] or external-surface-following [7] in the sense that the internal distribution automatically adapts to the external geometry.
This type of model is especially useful for modeling the shape of and optical changes in the lens with both age and accommodation. The inhomogeneity of the crystalline lens has long been known [8]; various sources in the extended literature have reported on experimental data [9–13], models [14–16] and finite ray tracing implementations [17–19] of the GRIN human lens.
2. General formulation
Here, we consider a singlet lens with two surfaces, the anterior and posterior surfaces. The signs of the curvature radii of both surfaces are usually used to classify the different types. When both signs are equal (+, +) or (-, -), the lens is a (positive or negative) meniscus, and when the two signs are opposite, the lens is either biconvex (+, -) or biconcave (-, +) lens.
2.1 Meniscus lens
For our formulation, the simplest case is the meniscus. We use cylindrical coordinates $({\omega ,\theta ,z} )$, but in the case of rotational symmetry around the optical axis, the formulation becomes 2-dimensional (2D) $({\omega ,z} )$, where ${\omega ^2} = {x^2} + {y^2}$ is the squared distance to the Z axis and z is the coordinate along the optical axis. In the meniscus, we consider the inner distribution of the refractive index as a function of the normalized distance along the optical axis ${z_0}$ to the anterior surface.
When the anterior ${R_a} = R({{z_0} = 0} )$ and posterior ${R_p} = R({{z_0} = t} )$ radii are different, there must be a gradient of curvature radii of the IISs to allow a smooth transition between these two external surfaces. In a previous study [3], a constant gradient -G was introduced so that the curvature radius of the IIS varied in a linear way:
Then, it is straightforward to show that $G = {{({{R_a} - {R_p}} )} / t}$ to match the values of both anterior and posterior radii of the meniscus. This means that the GRIN meniscus has two gradients of refractive index GRIN and curvature radius GRCU. Below, we use the acronym GRINCU to refer to this new type of two gradients lens.2.2 Iso-indicial surfaces
For the particular case of a rotationally-symmetric concentric distribution, we showed in Eq. (3) and Eq. (3) of a previous publication [4] that the conical iso-indicial surfaces (IISs) are given by the expression ${r^2} = {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 f}} \right.}\!\lower0.7ex\hbox{$f$}}\left( {\frac{{{{({z - \Delta } )}^2}}}{{{a^2}}} + \varepsilon \frac{{{\omega^2}}}{{{b^2}}} - O} \right)$ in cylindrical coordinates z, ω$({{\omega^2} = {x^2} + {y^2}} )$; a, b are the conic semiaxes; $\varepsilon ={+} 1$ for ellipses and $\varepsilon ={-} 1$ for hyperbolas; $\Delta = \varepsilon a$ means that the correct (left or right) branch of the conicoid intersects the optical axis at point $({z = {z_0},\omega = 0} )$. The variable r is a normalized (between 0 and 1) non-Euclidean distance to the center; the normalization is done by applying an offset O and a normalization factor f. For each value of r we have a different IIS. This generic expression applies to both anterior (Eq. (3)) and posterior (Eq. (3)) parts of the lens simply by assigning the corresponding values to the parameters. By simple operations, we can separate the standard expression of the conicoid on the left side:
To adapt the IIS expression to the present parametrization, we particularize Eq. (3) to the point of intersection with the optical axis $({z = {z_0},\omega = 0} )$, ${\raise0.7ex\hbox{${{{({{z_0} - \Delta } )}^2}}$} \!\mathord{\left/ {\vphantom {{{{({{z_0} - \Delta } )}^2}} {{a^2}}}} \right.}\!\lower0.7ex\hbox{${{a^2}}$}} = f{r^2} + O$, and solve for ${z_0}$:
Now, we can replace $f{r^2} + O$ with its expression in Eq. (3). In this way, we are replacing ${r^2}$ with the parameter ${z_0}$ in the expression of a generic IIS that intersects the optical axis at $({z = {z_0},\omega = 0} )$:Now, let us generalize this formulation to nonconcentric IISs. The simplest way to do this may involve introducing the curvature radius gradient parameter G. Below, we keep $Q({{z_0}} )= Q$ for the sake of simplicity. This means that we restrict ourselves to a meniscus with an external geometry determined by four independent variables (degrees of freedom) t, Ra, Rp and $Q = {Q_a} = {Q_p}$. In the case that ${Q_a} \ne {Q_p}$, one could introduce a constant gradient $- {{({{Q_a} - {Q_p}} )} / t}$ so that $Q({{z_0}} )= {Q_a} - ({{Q_a} - {Q_p}} ){{{z_0}} / t}$. Nevertheless, such generalization to three gradients GRIN/GRCU/GRCC (index/curvature/conic constant) is beyond the scope of the present study and will be the subject of future work.
With the gradient curvature parameter G, the geometry of the internal IIS is determined either by R or Q:
2.3 General conicoid
Nonrotationally symmetric optical elements such as spherocylindrical lenses are commonly used in visual optics, among other applications. Here, we consider a more general 3-axis conicoid, which is a natural generalization of the rotationally symmetric 2-axis case. The equation for the anterior surface ($f{r^2} + O = 1$) is
Therefore, with these definitions, the parameters of the rotationally symmetric case are equal to the mean of the maximum and minimum values. Finally, the surface can be rotated at any angle $\theta $ around the optical axis. If we apply the corresponding 2D rotation matrix ${\mathbf x^{\prime} = }{{\mathbf U}_{\mathbf \theta }}{\mathbf x}$ to the vector of coordinates ${\mathbf x} = {[{x,y} ]^T}$ we arrive at the expression equivalent to Eq. (3) for the generic sphero-cylindrical IIS.
2.4 Biconvex or biconcave lens
The formulation for these lenses is somewhat more complex than that for the meniscus, as the anterior and posterior radii have opposite signs. Not only the external surface but also each ISS has a positive and a negative side or “hemisphere”. Here, it is useful to consider the anterior and posterior sides independently and find the natural central interface of the ISS. This interface is an ellipse or a circumference for the rotationally symmetric case [4]. The axial thickness of the lens is the sum of the anterior and posterior thicknesses $t = {t_a} + {t_p}$, and the anterior and posterior refractive index distributions are functions of the normalized axial distances to their corresponding anterior or posterior external surface:
Now, we can obtain the equation for the central interface surface by replacing ${z_0}_a$ and ${z_0}_p$ with their expressions (Eq. (10)). However, it is straightforward to show that it is a fourth-order equation instead of a second-order conicoid as in the former case of concentric IISs [4].
We implemented two alternative analytical [21] and numerical solutions for the quartic equation. The detailed formulation of the analytical solution can be found in [22], and we include our custom implementation as we show in Code 1 (Ref. [23]). Nevertheless, we implemented a more efficient solution for practical numerical implementation based on transforming Eq. (15) in the associated inequalities for the arguments ${s_a} = \frac{{{t_a} - {z_0}_a({{z_i},{\omega_i}} )}}{{{t_a}}}$ and ${s_p} = \frac{{{z_0}_p({{z_i},{\omega_i}} )- {t_a}}}{{{t_p}}}$. For any point inside the lens, we have three possibilities: (1) ${s_a} = {s_p}$ means that we are at the interface; (2) ${s_a} > {s_p}$ means that the point belongs to the anterior part; and (3) ${s_a} < {s_p}$ means that the point belongs to the posterior part. In addition, a point is inside the lens if both arguments are less than or equal to 1: ${s_a} \le 1$ and ${s_p} \le 1$.
3. GRINCU model of the crystalline lens
In this section, we propose a new general (adaptive [4] or geometrically invariant [6] or external surface-following [7]) model of the crystalline lens based on the above formulation. Basically, we reformulate a previous model [4,24], particularizing the general expression in Eq. (11) for that particular case. That index distribution model was based on a power law, with a noninteger degree (exponent) 2p; the index distribution for the anterior side ${n_a}$ was monotonically increasing between a minimum value at the anterior surface, and conversely, ${n_p}$ was monotonically decreasing between a minimum value at the external surfaces ${n_s} = {n_a}(0 )= {n_p}(0 )$ and a maximum value at the central interface surface ${n_c} = {n_a}(1 )= {n_p}(1 )$.
3.1 Numerical implementation
We implemented two versions of the model in MATLAB (The MathWorks, Inc.) scripts, which permitted us to validate the formulation and was a necessary step toward future implementation of finite ray tracing through the lens model. For the general 3-dimensional (3D) model, we computed both $n({x,y,z} )$ and its gradient $\nabla n({x,y,z} )= ({{{\partial n} / {\partial x,}}{{\partial n} / {\partial y,}}{{\partial n} / {\partial z}}} )$, that is, the three partial derivatives in a 3D sampling grid. In the present formulation, we mainly used cylindrical coordinates so that we applied the chain rule to compute the partial derivatives with respect to the Cartesian coordinates from the partial derivatives with respect to the cylindrical coordinates. Since the calculations of these derivatives are standard and the final expressions are long, they are not included here in the text but are included in the code (in MATLAB notation).
We used a sampling resolution of 0.05 mm. This means that for typical lens dimensions (9 mm diameter and 4 mm thickness), we have 2,592,000 sampling points (with finer resolutions, the computation slows down dramatically in normal PC computers with the risk of memory overflow). In this 3D implementation, we avoided solving the quartic equation in 3D because it strongly increases the computing time. Instead, we included the 2D (rotationally symmetric case) analytical solution in the supplementary script because it is not difficult to generalize it to 3D and then obtain the complete central interface surface. Nevertheless, we applied the numerical solution (much faster) based on the inequalities to the arguments (as explained above) to classify every point among four different options:
- - if either ${s_a} > 1$ or ${s_p} > 1$, then the point is outside the lens; otherwise,
- - if ${s_a} = {s_p}$ the point belongs to the central interface; otherwise,
- - if ${s_a} > {s_p}$ the point belongs to the anterior part; otherwise
- - if ${s_a} < {s_p}$ the point belongs to the posterior part.
3.2 Results
The resulting models for a crystalline lens with different curvature gradients G, ranging from 0 to 3, are compared in Fig. 2 for the rotationally symmetric case. The upper row shows color maps of the refractive index distributions; the black line represents the central interface computed by the analytical method. The lower row shows IISs. The parameters used to compute these models are adapted from Ref. 19 for the unaccommodated lens of a 30-year-old: Ra = 10.96 mm, Qa = -3; Rp = -5.51 mm; t = 3.638 mm, ta = 0.6t; ns = 1.3726, nc = 1.4301, p = 3.741.
Figure 3 shows a transverse slice for the plane z = ta for a nonrotationally asymmetric case (with G = 2) with Da = 1.096 and θ = 30°. Note the elliptical shape with the meridian of maximum radius oriented along the 30° meridian.
The paraxial power of the lens models is represented as a function of the curvature gradient parameter G in Fig. 4. This approximated power was computed by the method described in Ref. 3 (it is included in the last part of the supplementary MATLAB script) using na = 1.3374 for the surrounding aqueous and vitreous humors.
Finally, the three partial derivatives of the gradient $\nabla n = ({{{\partial n} / {\partial x,}}{{\partial n} / {\partial y,}}{{\partial n} / {\partial z}}} )$ are shown in Fig. 5 for G = 2. The computation of the gradient is included here, as it is essential in ray tracing computations. Note the discontinuity in the three partial derivatives at the central interface, where the gradient reverses it sign from positive to negative values. This discontinuity could be avoided in different ways, as for instance using higher order models for the ISSs [6]. However experimental evidence suggests that such discontinuity could be present in real lenses [10].
4. Discussion and conclusions
A generalization of a previous adaptive or external-surface-following type of GRIN lens [4] has been proposed. The first and most important new feature was to include a curvature radius gradient parameter –G of the IIS, which allows us direct and explicit control in the formulation. In this way, we formulated a new type of lens (GRINCU) characterized by two gradients of refractive index (GRIN) and of the curvature of the IIS (GRCU). We believe that the multiplicative effect of the two gradients has a high impact on the optical performance of this type of lens. In fact, we found a mutual multiplicative effect of both gradients on the lens power [3]. As a result, in the case of biconvex (or biconcave) lenses, the second new feature is that the central interface surface between the two anterior and posterior sides of the lens is no longer a plane or a conicoid but a fourth-order surface to guarantee the continuity of the refractive index distribution at the interface. We implemented two alternative analytical and numerical solutions to solve the fourth-order equation to find the surface. As in the former model [4], this does not guarantee continuity in the gradient $\nabla n$. In fact, we can clearly see these discontinuities (sign reversals) in the partial derivatives toward the peripheral parts of the interface (Fig. 5). Therefore, one must be cautious when implementing ray tracing through these GRINCU lenses. One of the simplest solutions is to consider the anterior and posterior parts separately (as a sort of cemented doublet) for ray tracing [4]. Bahrami and Goncharov [6] added a third-order term to the second-order conicoid IIS, which was adjusted to avoid discontinuity in the derivatives.
The third relevant new feature is formulation and implementation of a general nonrotationally symmetric case to include the case of toric-like (spherocylindrical) lenses with two different curvatures along the principal meridians. This is an important feature when modeling biological lenses, such as the crystalline lens, and potentially important as well, for example, when designing intraocular lenses.
The last part of this study consisted of applying this formulation to implement a GRINCU nonrotationally symmetric model of the crystalline lens. We confirmed that the refractive power (diopters, D) increased as we increased the curvature radius parameter G. The power increased by approximately 13.8% when changing from G = 0 (24.25 D) to G = 3 (27.6 D).
To the best of our knowledge, this is the first case of GRINCU lens formulation with an explicit control parameter of the inner curvature gradient. This is the main novelty of the present model. Sheil and Goncharov [6] introduced a new parameter m in their AVOCADO lens model. This parameter was used to control a nonlinear scaling of both curvature radii R and conic constant (shape factor) Q of the IIS and hence the index gradient along the radial variable ω. The main difference is that the gradient $\nabla R ={-} G$ is constant in our model and Q = constant (versus two different nonlinear gradients for both R and Q in Ref. 6). Our goal was to obtain the simplest mathematical formulation, that is, second-order conicoid surfaces: linear scaling (constant gradient) of R and constant Q. There are two main reasons for choosing the simplest approach compatible with experimental data and evidence.
On the one hand, simple (linear or second-order) models strongly facilitate both geometrical and physical (optical) interpretation. On the other hand, the experimental measurement of the refractive index distribution of the lens is extremely difficult. A variety of non-destructive methods have been developed to measure the distribution of refractive index within the crystalline lens. They are based on optical fibers [9], magnetic resonance imaging (MRI) [10], optical coherence tomography (OCT) [12], X-ray interferometry [13] or laser ray tracing [25]. In most cases numerical models are needed for understanding the raw data obtained with the different techniques. For this reason we need to use the models with the least sophistication (maximum simplicity) to analyze those data. The same principle applies to the design of GRIN lenses. In this sense, we believe that further sophistication should be added only when required. For example, experimental evidence suggests that the human lens may exhibit toricity to different extents on its surfaces [12], which motivated us to generalize the results to nonrotationally symmetric conicoids (Subsection 2.3). Further work will include adding a third gradient of conic constant Q, implementing finite ray tracing, and applying the GRINCU lens to develop accommodating age-dependent optical models of the human eye.
Funding
Ministerio de Ciencia, Innovación y Universidades (PID2019-107058RB-I00).
Disclosures
The authors declare no conflicts of interest.
Data availability
No data were generated or analyzed in the presented research.
References
1. E. W. Marchand, Gradient Index Optics (Academic, 1978).
2. C. Gomez-Reino, M.V. Perez, and C. Bao, Gradient-Index Optics: Fundamentals and Applications (Springer-Verlag, 2002). [CrossRef]
3. R. Navarro and N. López-Gil, “Impact of internal curvature gradient on the power and accommodation of the crystalline lens,” Optica 4(3), 334–340 (2017). [CrossRef]
4. R. Navarro, F. Palos, and L. González, “Adaptive model of the gradient index of the human lens. I. formulation and model of aging ex vivo lenses,” J. Opt. Soc. Am. A 24(8), 2175–2185 (2007). [CrossRef]
5. J. W Blaker, “Toward an adaptive model of the human eye,” J. Opt. Soc. Am. 70(2), 220–223 (1980). [CrossRef]
6. M. Bahrami and A. V. Goncharov, “Geometry-invariant gradient refractive index lens: analytical ray tracing,” J. Biomed. Opt. 17(5), 055001 (2012). [CrossRef]
7. Q. Li and F. Fang, “Physiology-like crystalline lens modelling for children,” Opt. Express 28(18), 27155–27180 (2020). [CrossRef]
8. A. Gullstrand, “Appendix II,” in Handbuch der Physiologischen Optik, Helmholtz, 3rd ed. English translation edited by J. P. Southall (Optical Society of America, 1962) Vol. 1, 351–352.
9. B. K. Pierscionek, “Refractive index contours in the human lens,” Exp. Eye Res. 64(6), 887–893 (1997). [CrossRef]
10. C. Jones, D. Atchison, R. Meder, and J. 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). [CrossRef]
11. 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). [CrossRef]
12. 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). [CrossRef]
13. B. Pierscionek, M. Bahrami, M. Hoshino, K. Uesugi, J. Regini, and N. Yagi, “The eye lens: age-related trends and individual variations in refractive index and shape parameters,” Oncotarget 6(31), 30532–30544 (2015). [CrossRef]
14. G. Smith, D. A. Atchison, and B. K. Pierscionek, “Modeling the power of the aging human eye,” J. Opt. Soc. Am. A 9(12), 2111–2117 (1992). [CrossRef]
15. J. A. Díaz, C. Pizarro, and J. Arasa, “Single dispersive gradient-index profile for the aging human lens,” J. Opt. Soc. Am. A 25(1), 250–261 (2008). [CrossRef]
16. C. Sheil and A. Goncharov, “Accommodating volume-constant age-dependent optical (AVOCADO) model of the crystalline GRIN lens,” Biomed. Opt. Express 7(5), 1985–1999 (2016). [CrossRef]
17. H.-L. Liou and N. A. Brennan, “Anatomically accurate, finite model eye for optical modeling,” J. Opt. Soc. Am. A 14(8), 1684–1695 (1997). [CrossRef]
18. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24(8), 2157–2174 (2007). [CrossRef]
19. R. Navarro, “Adaptive model of the aging emmetropic eye and its changes with accommodation,” J. Vision 14(13), 21 (2014). [CrossRef]
20. W. Smith, Modern Optical Engineering (McGraw Hill, 2007).
21. M. Abramowitz and I. A. Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (Dover, 1972), chap. 3.
22. E. W. Weisstein, “Quartic Equation,” in MathWorld–A Wolfram Web Resource. https://mathworld.wolfram.com/QuarticEquation.html.
23. R. Navarro and S. Baquedano, “Numerical implementation of GRINCU model of the crystalline lens,” figshare (2021), https://doi.org/10.6084/m9.figshare.14865105
24. R. Navarro, F. Palos, and L. González, “Adaptive model of the gradient index of the human lens, II: optics of the accommodating aging lens,” J. Opt. Soc. Am. A 24(9), 2911 (2007). [CrossRef]
25. C. Qiu, B. Maceo Heilman, J. Kaipio, P. Donaldson, and E Vaghefi, “Fully automated laser ray tracing system to measure changes in the crystalline lens GRIN profile,” Biomed. Opt. Express 8(11), 4947–4964 (2017). [CrossRef]