Geometry-invariant GRIN lens: finite ray tracing - OSA Publishing

2 downloads 47 Views 1MB Size Report
Mehdi Bahrami∗1 and Alexander V. Goncharov2 ... C. E. Jones, D. A. Atchison, R. Meder, and J. M. Pope, “Refractive index distribution ... A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc.
Geometry-invariant GRIN lens: finite ray tracing Mehdi Bahrami∗1 and Alexander V. Goncharov2 1

Faculty of Science, Engineering and Computing, Kingston University London, Kingston-upon-Thames, UK 2 Applied Optics Group, National University of Ireland, Galway, Ireland ∗ [email protected]

Abstract: The refractive index distribution of the geometry-invariant gradient refractive index lens (GIGL) model is derived as a function of Cartesian coordinates. The adjustable external geometry of the GIGL model aims to mimic the shape of the human and animal crystalline lens. The refractive index distribution is based on an adjustable power-law profile, which provides additional flexibility of the model. An analytical method for layer-by-layer finite ray tracing through the GIGL model is developed and used to calculate aberrations of the GIGL model. The result of the finite ray tracing aberrations of the GIGL model are compared to those obtained with paraxial ray tracing. The derived analytical expression for the refractive index distribution can be employed in the reconstruction processes of the eye using the conventional ray tracing methods. The layer-by-layer finite ray tracing approach would be an asset in ray tracing through a modified GIGL model, where the refractive index distribution cannot be described analytically. Using the layer-by-layer finite ray-tracing method, the potential of the GIGL model in representing continuous as well as shell-like layered structures is illustrated and the results for both cases are presented and analysed. © 2014 Optical Society of America OCIS codes: (110.2760) Gradient-index lenses; (330.7326) Visual optics, modeling; (080.5692) Ray trajectories in inhomogeneous media.

References and links 1. R. K. Luneburg, Mathematical Theory of Optics (Brown University, 1944). 2. S. Nakao, S. Fujimoto, R. Nagata, and K. Iwata, “Model of refractive-index distribution in the rabbit crystalline lens,” J. Opt. Soc. Am. A 58, 1125–1130 (1968). 3. B. K. Pierscionek and D. Y C. Chan, “Refractive index gradient of human lenses,” Optom. Vis. Sci. 66, 822–829 (1989). 4. M. Hoshino, K. Uesugi, N. Yagi, S. Mohri, J. Regini, and B. Pierscionek, “Optical properties of in situ eye lenses measured with X-ray talbot interferometry: a novel measure of growth processes,” PLoS ONE 6, e25140 (2011). 5. B. K. Pierscionek, “Refractive index contours in the human lens,” Exp. Eye Res. 64, 887–893 (1997). 6. C. E. Jones, D. A. Atchison, R. Meder, and J. M. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. 45, 2352–2366 (2005). 7. P. P. Fagerholm, B. T. Philipson, and B. Lindstr¨om, “Normal human lens, the distribution of protein,” Exp. Eye Res. 33, 615–620 (1981). 8. A. Tardieu and M. Delaye, “Eye Lens proteins and transparency: from light transmission theory to solution X-ray structural analysis,” Annu. Rev. Biophys. Biomol. Struct. 17, 47–70 (1998). 9. B. K. Pierscionek and J. W. Regini, “The gradient index lens of the eye: an opto-biological synchrony,” Prog. Retin. Eye Res. 31, 332–349 (2012).

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27797

10. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24, 2157–2174 (2007). 11. J. A. D´ıaz, “ABCD matrix of the human lens gradient-index profile: applicability of the calculation methods,” Appl. Opt. 47, 195–205 (2008). 12. G. Beliakov and D. Y. C. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. I. Planar systems,” Appl. Opt. 36, 5303–5309 (1997). 13. G. Beliakov and D. Y. C. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. II. threedimensional systems with symmetry,” Appl. Opt. 37, 5106–5111 (1998). 14. W. S. Jagger, and P. J. Standsl, “A wide-angle gradient index optical model of the crystalline lens and eye of the octopus,” Vision Res. 39, 2841–2852 (1999). 15. 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, 21905–21917 (2010). 16. H. T. Kasprzak, “New approximation for the whole profile of the human crystalline lens,” Ophthalmic Physiol. Opt. 20, 31–43 (2000). 17. E. A. Hermans, P. J. Pouwels, M. Dubbelman, J. P. Kuijer, R. G. van der Heijde, and R. M. Heethaar, “Constant volume of the human lens and decrease in surface area of the capsular bag during accommodation: an MRI and Scheimpflug study,” Invest. Ophthalmol. Vis. Sci. 50, 281–289 (2009). 18. G. Smith, D. Atchison, D. Iskander, J. Jones, and J. Pope, “Mathematical models for describing the shape of the in vitro unstretched human crystalline lens,” Vision Res. 49, 2442–2452 (2009). 19. R. Urs, A. Ho, F. Manns, and JM. Parel, “Age-dependent Fourier model of the shape of the isolated ex vivo human crystalline lens,” Vision Res. 50, 1041–1047 (2010). 20. M. Bahrami and A. V. Goncharov, “Geometry-invariant gradient refractive index lens: analytical ray tracing” J. Biomed. Opt. 17, 055001 (2012). 21. M. Bahrami and A. V. Goncharov, “Geometry-invariant GRIN lens: iso-dispersive contours,” Biomed. Opt. Express 3, 1684–1700 (2012). 22. M. Bahrami and A. V. Goncharov, “The geometry-Invariant lens computational code,” This is a computable document format (CDF) for the equations presented in [20], http://optics.nuigalway.ie/people/mehdiB/CDF.html, (Oct. 2011). 23. G. Smith and B. K. Pierscionek, “The optical structure of the lens and its contribution to the refractive status of the eye,” Ophthalmic Physiol. Opt. 18, 21–29 (1998). 24. G. Smith, D. A. Atchison, and B. K. Pierscionek, “Modeling the power of the aging human eye,” J. Opt. Soc. Am. A 9, 2111–2117 (1992). 25. R. Navarro, F. Palos, and L. Gonz´alez, “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, 2175–2185 (2007). 26. 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, 2531–2540 (2008). 27. P. L. Ruben, “Aberrations arising from decentrations and tilts,” J. Opt. Soc. Am. 54, 45–46 (1964). 28. J. A. D´ıaz, J. Fern´andez-Dorado, and F. Sorroche, “Role of the human lens gradient-index profile in the compensation of third-order ocular aberrations,” J. Biomed. Opt. 17, 075003 (2012). 29. C. J. Sheil, M. Bahrami, and A. V. Goncharov, “An analytical method for predicting the geometrical and optical properties of the human lens under accommodation,” Biomed. Opt. Express 5, 1649–1663 (2014). 30. J. C. Wyant and K. Creath, “Basic wavefront aberration theory for optical metrology,” in Applied Optics and Optical Engineering, Volume 11, R. R. Shannon and J. C. Wyant, eds. (Academic Press, 1992), pp. 1–53. 31. J. Liang, B. Grimm, S. Goelz, and J. F. Bille, “Objective measurement of wave aberrations of the human eye with the use of a Hartmann-Shack wave-front sensor,” J. Opt. Soc. Am. A 11, 1949–1957 (1994). 32. D. Vazquez, E. Acosta, G. Smith, and L. Garner, “Tomographic method for measurement of the gradient refractive index of the crystalline lens. II. The rotationally symmetrical lens,” J. Opt. Soc. Am. A 23, 2551–2565 (2006).

1.

Introduction

Gradient refractive index (GRIN) lenses show a spatial, gradual variation in their index distribution. One celebrated example of a GRIN lens is the Luneburg lens, for which the external spherical surface and rotationally symmetric radial GRIN profile focus the monochromatic rays from a distant object at the lens back surface perfectly [1]. Although such a lens is difficult to manufacture, the concept shows that GRIN lenses have a great potential to attract optical designers’ attention for finding new solutions. The interest in studying GRIN media extends to biology science, since the crystalline lens of the human and animal eyes is also a GRIN lens [2–4]. The medium of a GRIN structure can be seen as constituted of great number of infinitely thin layers

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27798

with the same refractive index, so called iso-indicial contours [5, 6]. Iso-indicial contours in manufactured and biological GRIN lenses have finite thicknesses. For instant, the iso-indicial layers of the GRIN lens of the eye consist of protein and water molecules networks arranged in fiber layers [7–9]. Historically, in the absence of fast numerical calculations, paraxial optics has been a reliable approach for understanding the function and approximating the imaging quality of an optical system. Despite the improvements in computational hardware and software, which nowadays supports fast numerical ray tracing, GRIN lens models of the eye still employ paraxial-optics approximations and third-order aberration calculations for understanding the lens optical behavior [10,11]. However, the complexity of some models structure may require finite numerical ray tracing to examine the accuracy of the paraxial calculations and to recognize their limitations, specifically where optical effects at wide fields might not agree with the paraxial-optics assumptions. Despite the variety of efforts in presenting approximate analytical equations for the paraxial characteristics of the crystalline lens GRIN models, the approach in doing the finite ray tracing through these lenses has been always the same. This conventional approach employs different algorithms, which all require a well-behaved function of coordinates for the refractive index distribution of the lens [10, 12–15]. Limiting the refractive index distribution of the lens to profiles with analytical expressions might exclude a large family of GRIN media, which probably support more realistic representation for the gradient index in crystalline lens. In addition, employing a continuous, analytical expression for the refractive index of the lens does not accurately model ray tracing through the lenses with discrete profiles found in nature. Another common feature between the GRIN lens models supporting finite ray tracing is their use of conical surfaces. Pure conical surfaces are not the most realistic description for the external shape of the lens. There are many studies, which present more flexible geometries to reproduce the external shape of the crystalline lens [16–19]. Although geometries with a better fit to the empirical measurements are appreciated, yet the ray tracing through them can be time-consuming unless an analytical solution supports the calculation of their ray-surface intersection. In this paper we study finite ray tracing through a rotationally symmetric GRIN lens introduced as the geometry invariant GRIN lens (GIGL) model [20]. The optically well-defined external surfaces of this lens are described as conical surfaces with higher-order asphericity terms. Although this model and its external geometry aim at accurate modeling of the GRIN lens of the human eye, its flexible geometry may be used in animal eye studies, as well as in bio-inspired designs. This monochromatic model can be extended to a chromatic model, which is the topic of [21]. The GIGL model supports analytical paraxial ray tracing and provides closed-form expressions for Seidel aberrations of the lens, which is readily supported with computational codes [22]. In the current paper, we compare the aberration coefficients obtained from the paraxial optics equations of the GIGL model to the aberration coefficients calculated with finite ray tracing. This can help us to find the reliable range of the paraxial optics approximation. For our finite ray tracing calculations, we have developed a ray tracing method, which is based on layer-by-layer ray tracing through the well-defined geometry of the iso-indicial contours of the lens. This method can be also applied to more complicated GRIN structures, which do not support an analytical expression for their refractive index distribution, but still employ the surface representation of the GIGL model. In this study, the main benefit of a layer-by-layer ray-tracing, rather than using the traditional method, can be seen in the ability to compare the results for both cases (continuous and layered structure).

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27799

2.

Lens geometry

The external geometry of the GIGL consists of anterior and posterior surfaces meeting at the equatorial interface. The surfaces can be described as conicoids with an additional higher-order aspheric term. This extra asphericity guaranties that the anterior and posterior iso-indicial contours meet at the equatorial interface without any discontinuity in their first derivatives. This can be seen as the smooth connection between the outer anterior and posterior curves at the point (Zc , ρm ) in Fig. 1. The geometry of the iso-indicial contours in the GIGL is obtained by scaling the external shape of the lens. The internal iso-indicial contours radius of curvature r equals Rζ , where R is the external surface radius of curvature and ζ is the normalized distance from the center of the lens (point O). The same definition is used to connect t, the intercept of the internal iso-indicial contour on the optical z-axis measured from the origin O, to the external surface intercept T , such that t = T ζ . The connecting point of the anterior and posterior curves of the iso-indicial contours is located at zc . Similarly, it can be shown that zc = Zc ζ . For the mathematical description of the contours we have   2ra (ta + z) − (1 + ka)(ta + z)2 + ba (ta + z)3 −ta ≤ z < zc 2 (1) ρ =  2r p (t p − z) − (1 + k p)(t p − z)2 + b p(t p − z)3 zc ≤ z ≤ t p , where ρ is the distance from the optical z-axis, k is the conic constant of the surface, b is the higher-order asphericity constant, and subscripts a and p indicate the anterior and posterior parameters of the lens, respectively. It can be shown that the internal higher-order asphericity constant b = B/ζ , where B is the higher-order asphericity constant of the external surface. The mathematical description of the GIGL model and the derivation of the parameters are explained in detail in [20].

3.

Refractive index profile

The refractive index distribution of the geometry-invariant GRIN lens model is based on the power-law profile, which was originally proposed by Pierscionek [3] and later was formulated and supported by several studies [6, 23–26]. In the GIGL model [20], this GRIN profile is presented as: (2) n(ζ ) = nc + (ns − nc )(ζ 2 ) p ; where nc is the refractive index of the center of the lens, ns is the refractive index of the external surface of the lens, and the exponent p is a constant describing the steepness of the GRIN profile. To find the refractive index of the lens as a function of the lens external geometry at any point of the GRIN structure, we replace the internal coefficients in Eqs. (1) with their equivalent descriptions as   2Ra ζ (Ta ζ + z) − (1 + ka)(Ta ζ + z)2 + (Ba /ζ )(Ta ζ + z)3 −Ta ζ ≤ z < Zc ζ 2 ρ = (3)  2R p ζ (Tp ζ − z) − (1 + k p)(Tp ζ − z)2 + (B p /ζ )(Tp ζ − z)3 Zc ζ ≤ z ≤ Tp ζ .

By changing the normalized distance ζ from 1 to 0 in Eqs. (3), we can cover all iso-indicial contours from the surface to the center of the lens. Equationsp (3) should be solved for ζ . The solution is a function of the lens Cartesian coordinates ρ = x2 + y2 and z. Substituting the solution ζ (ρ , z) in Eq. (2) we can find the refractive index of the lens as a function of its coordinates.

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27800

ρ (Zc, ρm)

Ra

ρm

-R p

-rp

ra

O -T a

-t a

Zc zc

tp

Tp

z

Fig. 1. The geometry-invariant GRIN lens model geometry and its iso-indicial contours described in Eqs. (1). The Blue lines indicate the anterior part of the lens and the red lines specify the posterior part. The dashed lines depict the equatorial interface of the lens. In a 3-D presentation the equatorial interface is in the shape of a cone.

Using the first line of Eqs. (3) and bringing all terms to the left hand side of the equation give

ρ 2 − 2Raζ (Ta ζ + z) + (1 + ka)(Ta ζ + z)2 − (Ba /ζ )(Ta ζ + z)3 = 0.

(4)

By rearranging Eq. (4) we have

αa0 z3 + (αa1 z2 + ρ 2 )ζ + αa2 z ζ 2 + αa3 ζ 3 = 0

(5)

where the coefficients are defined by the parameters of its external shape:

αa0 = −Ba , αa1 = 1 + ka − 3BaTa , αa2 = −2Ra + 2(1 + ka)Ta − 3BaTa2 , and αa3 = −2RaTa + (1 + ka)Ta2 − Ba Ta3 . Similarly, the second line of Eq. (3) gives

α p0 z3 + (α p1 z2 + ρ 2)ζ + α p2 z ζ 2 + α p3 ζ 3 = 0,

(6)

where

α p0 = B p , α p1 = 1 + k p − 3B pTp , #220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27801

α p2 = 2R p − 2(1 + k p)Tp + 3B pTp2 , and α p3 = −2R pTp + (1 + k p)Tp2 − B p Tp3 .

Equations (5) and (6) are cubic equations of ζ . Solving a cubic equation can be done analytically in several ways. Here, to avoid complex numbers in the final ζ function coefficients, we shall use so called the trigonometric method. 3.1. Trigonometric method To use the trigonometric method in solving standard cubic equation described as

α0 + α1 ζ + α2 ζ 2 + α3 ζ 3 = 0,

(7)

first, the equation should be reduced to a monic trinomial. This can be done by substituting

ζ =χ−

α2 , 3α3

(8)

which after rearrangement reduces the standard cubic Eq. (7) to

β0 + β1 χ + χ 3 = 0,

(9)

where

β0 =

2α23 − 9α3 α2 α1 + 27α32 α0 , and 27α33

3α3 α1 − α22 β1 = . 2 p 3α3 Now we substitute χ by 2 cos(θ ) −β1 /3 in Eq. (9) and after rearranging we get s −3 3 β 0 = 0. 4 cos3 (θ ) − 3 cos(θ ) − 2β1 β1

(10)

Using the trigonometric identity: 4 cos3 (θ ) − 3 cos(θ ) − cos(3θ ) = 0,

(11)

and then comparing Eqs. (10) and (11) we find that 3β0 cos(3θ ) = 2β1 and

β1 χ j = 2 − cos 3 r

1 arccos 3

3β0 2β1

s

s

−3 β1

−3 , β1 !

2π −j 3

(12)

!

for j = 0, 1, 2.

(13)

Equation (8) can be used to convert this final answer χ to ζ . As long as the coefficients in Eq. (7) are real numbers the coefficients β0 and β1 in Eq. (13) will be real numbers. This is a great advantage in implementing the solutions in low-level programing languages, e.g. C. It is worth mentioning that χ can still be a complex number if the argument of the arccos function exceeds the interval (−1, +1) and/or the square root of the coefficients returns imaginary numbers.

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27802

3.2. Three dimensional GRIN profile Equation (13) does not show any complex coefficients for a real root, and numerical substitution shows that our normalized, real root for ζ (ρ , z) in Eqs (4) and (6) occurs for j = 0, therefore     s s 2 3 β − 9α α z ρ 2 3 α 4 z2 βa1 − 12αa3 ρ 2 1 z a0 a2 a3 a3  − αa2 z cos arccos ζ= 2 2 2 2 2 3 6αa3 (3αa3 ρ − z βa1 ) z βa1 − 3αa3ρ 3αa3 9αa3

(14) for the anterior hemisphere, and     s s 2 3 2 3 α 4 z2 β p1 − 12α p3 ρ 2 1 β α α ρ z − 9 z p3 p0 p2 p3  − α p2 z cos arccos ζ= 2 2 2 2 2 3 6α p3 (3α p3 ρ − z β p1 ) z β p1 − 3α p3ρ 3α p3 9α p3

(15)

for the posterior hemisphere of the GIGL model, where 2 − 3α α , βa1 = αa2 a1 a3 3 − 9α α α + 27α α 2 , βa0 = 2αa2 a1 a2 a3 a0 a3 2 − 3α α , and β p1 = α p2 p1 p3 3 − 9α α α + 27α α 2 . β p0 = 2α p2 p1 p2 p3 p0 p3

For instance, for Ra = R p = Ta = Tp = R and ka = k p = 0 we get Ba = B p = 0, αa0 = α p0 = 0, αa1 = α p1 = 1, αa2 = α p2 = 0, and αa3 = α p3 = −R2 . Hence considering ρ 2 = x2 + y2 , Eqs. (14) and (15) will be reduced to p x2 + y2 + z2 0 ≤ x2 + y2 + z2 ≤ R2 . (16) ζ (x, y, z) = R Equation (16) provides the normalized distance ζ from the origin for any points of a sphere with a radius R. Here the iso-indicial contours are concentric spheres with radius Rζ . Substituting Eq. (14) and (15) in Eq. (2) gives us the refractive index distribution of the GIGL model as a function of the lens coordinates with corresponding parameters as p n(x, y, z) = nc + (ns − nc ) ζ 2 (x, y, z) . (17)

Employing the analytical expression of the refractive index profile is the conventional way of ray tracing through a GRIN lens. The derived refractive index expression used in optimization algorithm is invaluable for reconstructing the lens based on in vivo wavefront or time of flight measurements. Since the lens of the eye is in fact consisting of small fibers with finite thicknesses, that would be beneficial to get a better understanding of discrete structures effect on the lens optical properties. Due to the optically well-behaved structure of the GIGL model, this continuous lens model can be reduced to a shell model. In a shell model the conventional method of GRIN lens ray tracing is not applicable. In the following a method for layer-by-layer ray tracing through the lens is developed to provide a comparison between the optical properties of the lens with different numbers of the shells.

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27803

(a)

(b)

ρ

n 1.415 1.404 1.392 1.381 1.369

z

-1.0

-0.5

0

ζ 0.5

1.0

Fig. 2. A typical geometry-invariant GRIN lens with its iso-indicial contours (a) and the corresponding normalized refractive index profile along radial direction from the center of the lens (b).

4.

Layer-by-layer ray tracing method

Figure 2 depicts a typical GRIN lens with its iso-indicial contours and the corresponding refractive index profile in an arbitrary radial direction, where ζ represents the normalized coordinate in 3-D. The first step in ray tracing is to find the ray-surface intersection and then to transfer the ray to this intersection point. 4.1. Finding ray-contour intersection Figure 3 depicts the schematic refraction of a ray passing through iso-indicial layers of the GRIN lens. The coordinates of the ray at the point P0 is known from the previous steps of ray tracing and the ray with direction cosines of (L, M, N) is emerging toward the next iso-indicial layer interface, where the refractive index changes from n0 to n1 and the ray is refracted to its new direction cosines of (L′ , M ′ , N ′ ). The main task is finding the distance between the point P0 = (x0 , y0 , z0 ) and the next intersection point P1 = (x1 , y1 , z1 ). This distance is indicated as a line segment ∆ in Fig. 3. It is clear from Fig. 3 that   x1 = x0 + L∆, y = y0 + M∆, (18)  1 z1 = z0 + N∆. If P1 is located in the anterior part of the GIGL, considering ρ 2 = x2 + y2 , the top equation in Eqs. (1) can be written as x21 + y21 = 2ra (ta + z1 ) − (1 + ka)(ta + z1 )2 + ba(ta + z1 )3 .

(19)

The coordinates x1 , y1 , and z1 in Eq. (19) can be substituted by Eqs. (18), and after rearrangement, we get δ0 + δ1 ∆ + δ2∆2 + δ3 ∆3 = 0, (20) #220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27804

ų0

ůŏŰŏ

ű





ůƢ ų1 ŏŰƢŏű

Ƣ

Ƒ0

Ƒ1

z

Fig. 3. Schematic refraction of a ray passing through iso-indicial layers of the GRIN lens.

where

δ0 = x20 + y20 − (ta + z0 ) {2ra + (ta + z0 ) [−1 − ka + ba(ta + z0 )]} , δ1 = 2(Lx0 + My0 ) + N {−2ra − (ta + z0 ) [−2 − 2ka + 3ba(ta + z0 )]} , δ2 = L2 + M 2 + N 2 [1 + ka − 3ba(ta + z0 )] , and δ3 = −ba N 3 . The solution of Eq. (20) is discussed in Section 3.1. A similar equation can be derived for the layers of the posterior hemisphere of the GIGL model. The optical path of the ray traveling the distance ∆ in Fig. 3 is n0 ∆. Accumulating the contribution of the layers to the optical path is used to calculate the wavefront shape of the beam and the wavefront aberration coefficients. To find the direction cosines of the ray after the refraction, we need the normal vector of the contour at the intersection point. It is worth mentioning that due to finite number of shells, there might be a possibility that a ray passing through n0 shell might cross it again in the posterior segment of the lens without entering the next inner shell of refractive index n1 . To recognise this situation one needs to check that the solution for ray intercept with the inner shell yields the coordinates of P1 within the anterior lens segment. If the ray solution indicates that the point P1 is located in the posterior lens segment, then the ray should be traced to the posterior intercept point of the current shell of refractive index n0 . The mathematical condition for point P1 (x, y, z) to be within the anterior segment of the lens is simply z < zc , where zc = ζ Zc , and ζ is the normalised distance for the shell of refractive index n as defined in Eq. (2). 4.2. Iso-indicial contours normal vector The normal vector of a surface described as F(x, y, z) = 0 is given by

∂F ∂F ∂F i+ j+ k ∂x ∂y ∂z N = s      . ∂F 2 ∂F 2 ∂F 2 + + ∂x ∂y ∂z

(21)

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27805

(a)

(b)

GRIN Lens

GRIN Lens Cornea

Cornea

Retina Iris

Retina Iris

Fig. 4. Numerical ray tracing for an aligned lens (a), and a tilted (8 degrees) and decentered (2 mm) lens (b). The tilt and decenteration are extremely exaggerated to make the effect more visible in the demonstration.

The anterior part of Eqs. (1) can be rearranged as x2 + y2 − 2ra (ta + z) + (1 + ka)(ta + z)2 − ba (ta + z)3 = 0.

(22)

Therefore the normal vector of the anterior iso-indicial contours of the lens is 2x i + 2y j + [−2ra + 2(1 + ka)(ta + z) − 3ba(ta + z)2 ] k N= p . 4(x2 + y2) + {2ra + (ta + z)[−2 − 2ka + 3ba(ta + z)]}2

(23)

The same approach can be used for the posterior iso-indicial contours of the lens. Knowing the intersection point coordinates, the incident ray direction cosines r, and the normal vector of the surface N, one can use the vector form of the Snell’s law to find the refracted ray direction cosines r′ : n0 (r × N) = n1 (r′ × N). (24) It is trivial that the calculations above can be also carried out for the outer surface of the lens. The concept is also valid for misaligned lenses. For instance, Fig. 4 shows the resultant ray tracing for a centered lens (a), and a tilted (8 degrees) and decentered (2 mm) lens (b).The rays refract through the cornea and the lens which focus the light at the retina. The tilt and decenteration in Fig. 4 are larger than typical values found in the eye and are chosen here just to exaggerate the effect. The ray paths shown in Fig. 4 has been generated by ray tracing in Zemax software program with a user defined GRIN lens surface. 4.3. Examining the ray tracing code The nature of the GIGL model in terms of its external shape and its refractive index distribution is very different from the built-in surfaces and GRIN media in optical design software. Yet for a symmetric geometry the complex GIGL model can be reduced to built-in GRIN profiles with standard conic surfaces. To check the ray tracing and aberration calculation of the new ray tracing code, a spherical lens with radius of R and Gradient 4 refractive index profile in Zemax optical design software is used. We consider a parabolic refractive index profile for this GRIN lens, defined as:  ns − nc  2 x + y2 + (z′ − R)2 , (25) n(x, y, z′ ) = nc + R2 #220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27806

where z′ is defined as z + R to adjust the difference between the coordinate origin convention in Zemax and the GIGL model definitions. Such a spherical lens is identical to the profile derived in Eq. (16) if p = 1 in Eq. (17). Comparing the wavefront coefficients and ray tracing results from the built-in GRIN lens and the GIGL ray tracing code helped us to verify that our code performed correctly.

5.

Seidel and Zernike aberration coefficients

The GIGL model supports analytical paraxial ray tracing leading to closed-form expressions for Seidel aberration coefficients [21]. As we explained before in the introductory article of the GIGL model [20], despite the advantages of the Seidel coefficients in understanding the nature of aberration compensation in the crystalline lens, the third-order aberration calculations do not support tilted or decentered elements, which is the case in the human eye. There have been efforts to modify Seidel theory and use it for tilted or decentered elements [27]. Although this theory is employed in one of the GRIN lens models of the eye [28], tilt and decenter in the eye are not the main limitation for the classical Seidel theory. The external shape of the lenses modeling the crystalline lens can be divided into two groups: models based on pure conic surfaces and models with more advanced geometry using higher order asphericity. This higher asphericity is not taken into account in the Seidel aberration coefficients. For models, with their internal GRIN structure mimicking the lens external shape, the higher-order asphericity terms have their own contribution to the refractive index profile. This means that even the paraxial rays experience the higher-order asphericity in the internal structure of the lens. This can be seen clearly in Fig. 1, where the internal iso-indicial contours are the scaled replica of the external shape of the lens. It is possible to approximately account for B coefficient in paraxial ray tracing by using a modified version of conicoid of revolution [29]. One well-accepted alternative to Seidel polynomials is using Zernike expansion. Zernike coefficients calculation is based on real ray tracing and can account for tilt and decenter in the lens. In the following, the connection between Seidel and Zernike aberration coefficients is discussed and the discrepancy between the two paraxial and real ray-tracing aberration coefficient calculations for the GIGL model is analyzed with numerical examples. A wavefront can be expanded based on Seidel polynomials as W (ρ , θ ) = W040 ρ 4 + W131 H¯ ρ 3 cos θ + W222 H¯ 2 ρ 2 cos2 θ + . . .,

(26)

where ρ and θ are polar coordinates, H¯ is the fractional image height, which varies between 0 to 1, and W040 , W131 , and W222 are wavefront coefficients. Here, we limit our calculation to the widest field of view, where normalized field angle H¯ = 1. The wavefront coefficients can be found using Seidel aberration coefficients as: W040 =

SI SII SIII ; W131 = ; W222 = , 8 2 2

(27)

where SI , SII and SIII are Seidel sums for the third-order spherical aberration, coma and astigmatism. The wavefront can be also expanded based on Zernike polynomials. With no defocus at the image surface, by fitting the wavefront to the first eleven Zernike terms and limiting the filed of view to y direction we have [30]: √ √ √ (28) W040 = 6 5 Z11 ; W131 = 3 8 Z7 ; W222 = −2 6 Z6 , where Z11 , Z7 , and Z6 are standard Zernike coefficients. The aberration calculations for the GIGL model in [20,22] are derived for the Seidel aberration coefficients. These calculations are #220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27807

limited to the paraxial region, where rays “see” the iso-indicial contours of the lens as standard conic surfaces. Finite ray tracing in optical design software Zemax provides Zernike aberration coefficients and allows us to calculate exact wavefront characteristics for non-paraxial regions of the lens. A comparison between the lens aberrations calculated based on the Seidel and Zernike coefficients can provide valuable insight into the accuracy of paraxial calculations for wider apertures and field of views. For this purpose we define the quantity √ SI − 6 5 Z11 × 100%, ∆W040 = 8 SI 8

(29)

as the difference (in percentage) between the Spherical aberration calculated from paraxial Seidel coefficients and Zernike coefficients. Similarly, the quantities ∆W131 and ∆W222 , respectively stand for the deviation in coma and in astigmatism calculations of the lens. This comparison is presented in Table 1 for different apertures and field of views. The optical characteristics of the lens used for numerical calculations in this paper can be summarized as: Ta = 2.1 mm, Tp = 1.4 mm, Ra = 11.0 mm, R p = 7.5 mm, ka = −2.0, k p = −3.0, p = 3.0, ns = 1.370, and nc = 1.415. The lens is surrounded by a medium with a refractive index of 1.336 and the incident rays are not refracted, for instance by the cornea, prior to the lens. Spherical aberration error in Table 1 for the aperture of 0.5 mm indicates that the contribution Table 1. The difference (in percentage) between the aberrations calculated from paraxial Seidel coefficients and Zernike coefficients for spherical aberration, ∆W131 , coma, ∆W131 , and astigmatism, ∆W222 , see Eq. (29).

∆W040 (%)

∆W131 (%)

∆W222 (%)

4.4

10.7

8.3

7.5

5.1

7.3

10.7

21.8

0.5

29.0

30.5

35.3

9.6

10.1

11.9

93.1

94.3

103.4

2.0

87.8

88.7

91.4

45.6

45.9

46.7

255.4

257.6

264.6

4.0

1.0

4.0

8.0

4.0

8.0

1.0 4.0 8.0 1.0 Full field of view (degree)

Pupil diameter (mm)

2.3

of the higher-order asphericity to the error originated from the internal iso-indicial contours is negligible for the paraxial region, but is rapidly increasing for wider pupils. Note that the Zernike aberration coefficients estimated by exact ray tracing remain small at 2 mm pupil diameter, and should not affect the image quality. The Seidel spherical aberration is constant across the field, but as can be seen, Table 1 shows a field dependency in the spherical aberration error. The apparent aperture for the oblique rays is decentered, such that these rays will be more affected by higher-order asphericity in comparison to the paraxial rays. This also means that tilt and decenter of the lens increases the error in the third-order calculations. The data in Table 1 confirms that the image quality calculated from Seidel aberration theory shows a noticeable difference from the exact ray tracing calculations. Therefore the reliability of the classical or the modified version of the Seidel aberration theory in non-conic aspheric GRIN profiles is limited to paraxial apertures and small fields.

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27808

Number of shells 1

10

100

1,000

10,000

100,000

Effective focal length (mm)

51.260 51.460 51.660 51.860 52.060 52.260 52.460 52.660 52.860 53.060 53.260

Fig. 5. The impact of the number of shells on effective focal length.

6.

Discrete refractive index profiles

Although the GRIN lenses with analytical expression for their refractive index profile represent a continuous change in the refractive index across the lens, in numerical ray-racing we need to decide on the number of refractions through the GRIN structure. As we discussed, the refractions take place at the iso-indicial contours, where the refractive index changes. Based on the layer-by-layer ray tracing approach a finite number of refractions reduces the continuous GIGL model to a shell-model lens. The number that we assume for these shells has a direct effect on the optical characteristics of the lens. In Fig. 5 we illustrated the change in the effective focal length (EFL) of the lens with the increasing number of shells. The lens optical characteristics are the same as the one described in Sec. 5. The diagram starts at a two-shell lens, where only one contour separates the refractive indices of surface and the center of the lens. At very large number of shells the diagram gets closer and closer to the EFL value calculated for the continuous profile [22]. Increasing the number of shells from 2 to a continuous case the lens shows a 2 mm change in its EFL; the more number of shells is used, the more optically powerful lens becomes. For a lens with 300 shells the deviation from the continuous lens EFL becomes smaller than 0.1 mm. The change in the GRIN structure of the crystalline lens in human and animal eyes is not continuous. The smallest iso-indicial unit in this medium is constituted from proteins-water network, which their regular orientation provides a transparent GRIN lens. The structure of these networks and their dimension decide on the thickness and the density of iso-indicial shells. A better understanding in optical performance of this structure can help us to chose the number of shells in our numerical calculations to achieve a more accurate aberration prediction in the eye.

7.

Discussion and Conclusion

Instead of starting from an analytical expression for the gradient refractive index of the lens, we started from a well-defined geometry for the iso-indicial contours of the lens and derived

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27809

the refractive index distribution from this geometry. The complexity of this derived analytical expression for the index profile of the GIGL model shows that such a geometry for the interior structure of the lens could not be achieved by a clever guess for the refractive index profile. The conventional method for ray tracing through a GRIN medium requires the refractive index expression and its derivatives within the lens. We developed our ray tracing method based on a layer-by-layer approach, which requires a well-defined geometry for the iso-indicial contours. This guaranties that any future GRIN lens model with iso-indicial contours description similar to the external geometry of the GIGL model can be numerically ray-traced, regardless of availability of an analytical expression for its refractive index profile. In addition, the effects originated from the discrete refractive index profiles can be also studied in such an approach. The finite ray tracing through the GIGL model is used for the aberration coefficient calculation and these coefficients are numerically compared to the ones, which are derived based on paraxial optics approximations. Although for small fields and apertures the approximate coefficients show less than 10 percent error, this deviation for wider angels and apertures is dramatically increased. The source of this difference is higher asphericity in the iso-indicial contours geometry, which is not taken into account in paraxial optics. This indicates that more realistic models of the crystalline lens one should check the accuracy of the paraxial calculations by finite numerical ray tracing. The main advantage of the geometry invariant GRIN lens model used here is the ability to calculate aberration coefficients analytically (with both paraxial and exact ray-tracing). However this requires that the power coefficient p remains constant throughout the lens, yet real crystalline lenses might show non-constant power coefficient depending on the direction of propagation within the lens. It certainly needs to be accepted that all models of biological systems, including the one in this paper, aim to represent the physiological (natural) state, but cannot entirely replicate it. It is worth mentioning that the derived analytical expression for the GIGL model refractive index profile and also the discussed ray tracing code present a great tool for the efforts in reconstruction of the crystalline lens. The GIGL model is not limited to a plane equatorial interface, or does not show sharp tips at its iso-indicial contours. These all can give a more realistic structure of the GRIN lens in the various reconstructions based on ocular aberration measurements [31], magnetic resonance imaging [6, 26], optical tomography imaging [32], optical coherence tomography imaging [15], or refractive index measurements by X-ray Talbot interferometry [4]. Acknowledgments The authors gratefully acknowledge Professor Chris Dainty for his support, we also wish to thank referees for their constructive comments.

#220583 - $15.00 USD Received 7 Aug 2014; revised 22 Oct 2014; accepted 22 Oct 2014; published 3 Nov 2014 (C) 2014 OSA 17 November 2014 | Vol. 22, No. 23 | DOI:10.1364/OE.22.027797 | OPTICS EXPRESS 27810