A model of the entrance pupil of the human eye

The aperture stop of the iris is subject to refraction by the cornea, and thus an outside observer sees a virtual image: the “entrance pupil” of the eye. When viewed off-axis, the entrance pupil has an elliptical form. The precise appearance of the entrance pupil is a consequence of the anatomical and optical properties of the eye, and the relative positions of the eye and the observer. This paper presents a ray traced model eye that provides the parameters of the entrance pupil ellipse for an observer at an arbitrary location. The model is able to reproduce empirical measurements of the shape of the entrance pupil with good accuracy. I demonstrate that accurate specification of the entrance pupil of a stationary eye requires modeling of corneal refraction, the misalignment of the visual and optical axes, and the non-circularity of the aperture stop. The model, including a three-dimensional ray tracing function through quadric surfaces, is implemented in open-source MATLAB code.

www.nature.com/scientificreports www.nature.com/scientificreports/ that the radius of curvature at the vertex of the cornea varies as a function of spherical ametropia, and provided the parameters for a rotationally-symmetric, prolate ellipsoid, with the radius of curvature (but not asphericity) varying linearly with the refractive error of the modeled eye. I calculated the proportional change in the Navarro values that would correspond to the effect of ametropia described by Atchison. This yields the following expression for the radii of the corneal front surface ellipsoid: . . ⋅ − . cornea front radii SR SR _ _ ( ) [14 26 10 43 10 27] (1 0 0028 ) (1) where SR is spherical refractive error in diopters, and the radii are given (here and subsequently) in the order of axial, horizontal, vertical. The apex of the corneal ellipsoid is on the optical axis. Atchison 10 adopted parameters for the back surface of the cornea that do not vary by ametropia. Navarro and colleagues 9 do not provide posterior cornea parameters. Therefore, I scale the parameters provided by Atchison 10 to be proportional to the axial corneal radius specified by Navarro, yielding: . . cornea back radii _ _ [13 7716 9 3027 9 3027] (2) The center of the back cornea ellipsoid is positioned so that there is 0.55 mm of corneal thickness between the front and back surface of the cornea at the apex, following Atchison 10 .
Iris and aperture stop. The iris is modeled as a plane perpendicular to the optical axis (i.e., zero "iris angle"), positioned 3.9 mm posterior to the anterior surface of the cornea, at the location of the anterior surface of the lens for the cycloplegic eye 11 .
The inner rim of the iris is the source of the image of the border of the entrance pupil. There is evidence that the entrance pupil is decentered nasally with respect to the optical axis, and that the center shifts temporally with dilation (reviewed in 12 ). The current model does not attempt to incorporate these measurements into the properties of the aperture stop. Instead, the aperture stop is modeled as centered and fixed on the optical axis.
The stop is, however, modeled as slightly non-circular. Wyatt 13 measured the average (across subject) ellipse parameters for the entrance pupil under dim and bright light conditions (with the visual axis of the eye aligned with camera axis). The pupil was found to be elliptical in both circumstances, with a horizontal long axis in bright light, and a vertical long axis in dim light. The elliptical eccentricity was less for the constricted as compared to the www.nature.com/scientificreports www.nature.com/scientificreports/ dilated pupil (non-linear eccentricity ε = 0.12 horizontal; 0.21 vertical; where the ratio of the minor to the major axis of an ellipse is equal to ε − 1 2 ). In the current model, I limit the eccentricity of the dilated entrance pupil to a slightly lower value (0.18) to provide a better fit to the measurements of Mathur and colleagues 6 (discussed below). Accounting for the refraction of the cornea (using ray tracing simulations), I calculated the ellipse parameters (area, eccentricity, and tilt) for the aperture stop that would produce the entrance pupil appearance observed by Wyatt 13 at the two light levels. A hyperbolic tangent function of stop radius is then used to transition between the elliptical eccentricity of the constricted stop-through circular at an intermediate size-to the eccentricity of the dilated stop. Consistent with Wyatt 13 , the tilt of the elliptical stop is horizontal (θ = 0) for all stop sizes below the circular intermediate point, and then vertical with a slight tilt of the top towards the nasal field for all larger stops.
The model therefore provides a stop that smoothly transitions in elliptical eccentricity and horizontal-to-vertical orientation (passing through circular) as the radius of the aperture varies from small to large following these functions: where r is the radius of the aperture stop in mm, ε is the non-linear eccentricity of the aperture ellipse, and θ is the tilt of the aperture ellipse (horizontal = 0; vertical = pi/2). The sign of ε is used to determine θ, although the absolute value of ε is used to generate the stop ellipse.
Although not required to determine the entrance pupil, the boundary of the visible iris is modeled for comparison with empirical images of the eye. The horizontal visible iris diameter (HVID) has been reported to be 11.8 mm 14 . Bernd Bruckner of the company Appenzeller Kontaktlinsen AG measured HVID using photokeratoscopy in 461 people 15 and kindly supplied me with the individual subject measurements. These data are well fit with a Gaussian distribution and yield a mean visible iris radius of 5.91 mm and a standard deviation of 0.28 mm. After accounting for corneal magnification, the radius of the iris in the model is set to 5.57 mm.
Lens. The crystalline lens is modeled as a set of quadric surfaces. The anterior and posterior surfaces of the lens are modeled as one-half of a two-sheeted hyperboloid. The radii of these surfaces, their dependence upon age and the accommodative state of the eye, and thickness of the anterior and posterior portions of the lens, are taken from Navarro and colleagues 16 . The interior of the lens is modeled with a gradient of refractive indices. This is realized as a set of ellipsoidal iso-indicial surfaces; 40 surfaces were used in the results presented here. The particular parameters are taken from Atchison 10 , which were in turn derived from prior studies 17,18 . The lens is modeled as center-aligned with the optical axis. The axial position of the lens center is taken from Atchison 10 .
The accommodative state of the model eye is set in terms of the distance of the point of best focus, established by adjusting the parameter "d" in the equations of Navarro and colleagues 16 .
Vitreous chamber and retina surface. Atchison 10 provides the radii of curvature and asphericities of a biconic model of the vitreous chamber, which is decentred and tilted with respect to the visual axis. With increasing negative refractive error (myopia), the vitreous chamber changes size, elongating in the axial direction and to a lesser extent in the vertical and horizontal directions. I model the vitreous chamber as an ellipsoid that is centered on and aligned with the optical axis, deriving the radii, and their dependence upon spherical refractive error, from equations 22-24 provided by Atchison 10 .
Visual and line-of-sight axes. Our goal is to describe the appearance of the entrance pupil. Often, the pupil is observed while the subject under study has their eye directed at a fixation point. To model this, I define the location of the fovea on the retina and the visual and line-of-sight axes.
The visual axis is a nodal ray 19 of the eye that intersects the fovea. The angle α is the angle between the visual axis of the eye and the optical axis 12 . The visual axis is generally displaced nasally and superiorly within the visual field relative to the optical axis. The horizontal displacement of the visual axis is 5-6°2 0 , vertical displacement is on the order of 2-3°2 1 . There is individual variation in these values.
Tabernero and colleagues 20 proposed that individual differences in the axial length of the eye are related to individual differences in α. Their model considers that, as the axial length of the eye increases, the axial distance of the fovea from the nodal points of the eye increases more than does the horizontal or vertical distance. This causes the angle between the visual and optical axes to grow smaller as the eye grows longer. In the current model, I adopt this theoretically motivated relationship with the form: www.nature.com/scientificreports www.nature.com/scientificreports/ When it is necessary to model the eye as focused upon a particular point in space, I obtain the line of sight for the model, which is the ray that connects the fovea and the fixation point via the center of the entrance pupil.
Perspective projection and ray tracing. Given the model eye, I then specify the appearance of the pupil in an image in terms of the parameters of an ellipse. This is accomplished by perspective projection of the eye to the sensor of a simulated camera, and then fitting an ellipse to points on the border of the pupil image. This projection requires us to specify the intrinsic properties of a camera and its position in space.
Perspective projection. A pinhole camera model is defined by an intrinsic matrix. This may be empirically measured for a particular camera using a resectioning approach (also termed "camera calibration") 22,23 . If provided, the projection will model radial lens distortion 24 , although the present simulations assume an ideal camera. The position of the camera relative to the eye is specified by a translation vector in world units, and the torsional rotation of the camera about its optical axis.
With these elements in place, the model defines points on the boundary of the aperture stop that are separated by equally spaced angles with respect to the stop center; these points are slightly unequally spaced in terms of their arc-length distance on the stop border, but this is of little practical consequence. The size of the aperture stop is specified as a radius in millimeters. As the stop is elliptical, radius is defined in terms of a circle of equivalent area. Five boundary points uniquely specify an ellipse and six or more are required to measure the deviation of the shape of the pupil boundary from an elliptical form. Sixteen points were used in this simulation, matching the number of points used by Mathur and colleagues to measure the entrance pupil ellipse 6 . The boundary points are projected to the image plane and then fit with an ellipse 25 , yielding the ellipse parameters.
Ray tracing. As viewed by the camera, the stop boundary points are subject to corneal refraction, and thus are virtual images. The cornea causes the imaged pupil to appear larger than its actual size, and displaces and distorts the shape of the pupil when the eye is viewed off-axis 3,7 . A similar and additional effect is introduced if the eye is observed through corrective lenses worn by the subject (e.g., contacts or spectacles).
To account for these effects, the projection incorporates an inverse ray tracing solution. I model the cornea of the eye and any corrective lenses as quadric surfaces. The set of surfaces (with their radii, rotations, translations, and refractive indices) constitute the optical system. The routines provide refractive indices for the aqueous humor and cornea, as well as values for materials that are present in the optical path as artificial corrective lenses (e.g., hydrogel, polycarbonate, CR-39). The specific refractive index given the wavelength used to image the system is derived using the Cauchy equation, with parameters taken principally from Navarro and colleagues 16 . A wavelength of 550 nm was used in the current simulations.
For the specified optical system, I implement 3D tracing for skew rays through a system of quadric surfaces 26 . Given a boundary point on the stop, the angles with which the ray departs from the optical axis of the eye, and the spatial translation of the camera, we can compute how closely the ray passes the pinhole aperture of the camera (Fig. 2). We search for the ray that exactly intersects the camera aperture and thus will uniquely be present in the resulting image (assuming a pinhole). The point of intersection of this ray at the last surface of the optical system is treated as the virtual image location of this point. In practice, the median ray trace solution finds a ray that passes within 0.0001 mm of the modeled camera pinhole. While greater precision could be obtained with longer computation time, this provides little practical benefit for model accuracy.

entrance pupil of a stationary eye
The ellipse parameters of the entrance pupil produced by the current model may be compared to empirical measurements. Mathur and colleagues 6 obtained images of the pupil of the right eye from 30 subjects. Each subject fixated a target 3 meters distant. A camera was positioned 100 mm away from the corneal apex and was moved to viewing angle positions in the temporal and nasal visual field within the horizontal plane. A viewing angle of zero was assigned to the position at which the line-of-sight of the subject was aligned with the optical axis of the camera. At each angular viewing position, the camera was translated slightly to be centered on the entrance pupil. An image of the pupil was obtained and then the border of the pupil marked by hand with 16 points, and these points fit with an ellipse. The pupil of the subject was pharmacologically dilated (and the eye rendered cycloplegic) with 1% cyclopentolate. The subjects had a range of spherical refractive errors, but the central tendency of the population was to slight myopia (−0.823 D). Figure 3A illustrates the geometry of the measurement.
The measurements of Mathur and colleagues 6 can be simulated in the current model. I created a right eye with a spherical refractive error of −0.823 D, and assumed an entrance pupil diameter of 6 mm to account for the effect of cyclopentolate dilation 27 (accounting for corneal magnification, this corresponds to a 5.3 mm diameter aperture stop). I set the accommodation of the eye such that the point of best focus was at 3 meters and rotated the eye so that the line-of-sight was aligned with the optical axis of the camera at a visual field position of zero. I then positioned the simulated camera as described in Mathur and colleagues 6 . Figure 3B presents renderings of the model eye as viewed from the nasal and temporal visual field. I evaluated the model for viewing angles between −75 and 65. As is described below, the ray tracing performance of the model declined rapidly beyond these limits.
Pupil diameter ratio. At each viewing angle Mathur and colleagues measured the ratio of the minor to the major axis of the entrance pupil, which they termed the pupil diameter ratio (PDR; Fig. 3C). These measurements were combined across subjects and accurately fit by a parameterized cosine function. Figure 4A shows the resulting fit to the empirical data. The PDR decreases as a function of viewing angle, as the entrance pupil is tilted away from the camera and the horizontal axis is foreshortened. The observed entrance pupil follows a cosine function that is "both decentered by a few degrees and flatter by about 12% than the cosine of the viewing angle 6  We may examine how components of the current model contribute to the replication of the empirical result. The simulation was repeated with successive addition of model components. The simplest model assumes a circular stop, no corneal refraction, and alignment of the optical axis and line-of-sight of the eye. Under these circumstances the entrance pupil behaves exactly as a function of the cosine of the viewing angle (Fig. 5A). Next, we rotate the eye so that the line-of-sight is aligned with the optical axis of the camera when the viewing angle is zero (Fig. 5B). Modeling of α (the misalignment of the visual and optical axes) decenters the PDR function. Addition of refraction by the cornea (Fig. 5C) flattens the PDR function, as the entrance pupil now demonstrates magnification along the horizontal axis with greater viewing angles. Finally, the model is augmented with an elliptical, tilted stop (Fig. 5D). Because the dilated stop is itself taller than it is wide, the PDR function is reduced overall, bringing the model into final alignment with the empirical result.
Variation in the pupil diameter ratio function. Mathur and colleagues 6 found that the empirical PDR function could be fit with great accuracy by the expression: where ϕ is the viewing angle in degrees, D is the peak PDR, β is the viewing angle at which the peak PDR is observed (e.g., the "turning point" of the function), and 180 E is the half-period of the fit in degrees. When fit to the output of the current model, we recover the parameters D = 0.99, β = −5.36, and E = 1.14, with an R 2 = 0.99. These parameters are quite close to those found 6 in the fit to empirical data (0.99, −5.30, and 1.12). We may consider how biometric variation in the model changes the observed entrance pupil in terms of these parameters. The value D reflects the eccentricity and tilt of the stop ellipse (Eq. 3) and the vertical component of α. Wyatt 13 found individual variation in pupil shape, and by extension likely variation in the shape of the aperture stop in the iris. While variation in ellipticity of the stop will produce variations in D, overall this model component had a rather small effect upon the simulated PDR (Fig. 5).
The value β is determined by the horizontal component of the misalignment of the visual and optical axes of the eye (angle α). This angle can vary substantially by subject and may have some relation to spherical refractive error. Figure 6A illustrates the PDR function associated with two different β values. Mathur and colleagues 6 obtained β for each of their 30 subjects and related this to subject variation in spherical refractive error. They found a negative correlation between these values, although this did not strongly account for the variance to the optical axis of the eye. The optical system is composed of the aqueous humor, the back and front surfaces of the cornea, and the air. We consider the set of rays that might originate from the edge of the aperture stop. Each of these rays depart from the stop border at some angle with respect to the optical axis of the eye. We can trace these rays through the optical system. Each ray will pass at a varying distance from the pinhole aperture of the camera (blue+). For a pinhole camera, only the ray that strikes the camera stop exactly (i.e., at a distance of zero) will contribute to the image. We therefore search across values of initial angle to find the ray that minimizes the intersection point distance. In this example system, a ray that departs the iris stop border at an angle of −30° with respect to the optical axis of the eye strikes the pin-hole of the camera after undergoing refraction by the cornea. Inset bottom left. The path of the ray after it exits the last surface of the optical system (in this case, the front surface of the cornea) creates a virtual image of the border point of the aperture stop that is displaced (asterisk) from the veridical location. (2019) 9:9360 | https://doi.org/10.1038/s41598-019-45827-3 www.nature.com/scientificreports www.nature.com/scientificreports/ (R 2 = 0.28). Figure 6B replots their data (gray points) and linear fit (black line). In the current model, the value assigned to α based upon spherical refractive error is made to follow the formula proposed by Tabernero and colleagues 20 , producing the dependence of turning point upon spherical refractive error shown (red line). The current model also poorly accounts for individual variation in the value of β (R 2 = 0.16).
The value E varies the "width" of the PDR function (Fig. 7A). In their simulation, Fedtke and collagues 7 found that increasing the modeled aperture stop slightly widened the PDR function (i.e., the parameter E increased). I find the same effect in the current model (Fig. 7B, top). Changes in the depth of the iris relative to the corneal surface had a minimal effect, suggesting that changes in accommodation (apart from a change in stop size) will have little effect upon the shape of the entrance pupil. I examined as well the effect of corneal radius as measured in the horizontal plane (Fig. 7B, bottom). A modest change in the width of the PDR function is found for variation in the assumed corneal power. Given that the corneal surface is astigmatic, this indicates that there will be a small but measurable difference in the appearance of the entrance pupil along the vertical and horizontal dimensions.
Tilt of the entrance pupil ellipse. Mathur   The iris border and pupil border points are virtual images, having been subjected to refraction from the cornea. Iris border points are missing in some views as these points encountered total internal reflection during ray tracing. (C) The pupil diameter ratio (PDR) is the ratio of the minor to the major axis of an ellipse fit to the entrance pupil. θ is the tilt of the major axis of the pupil ellipse. This value of θ is the empirically observed tilt of the entrance pupil ellipse, which is not to be confused with the specified value θ that is the modeled tilt of the elliptical aperture stop in Eq. 3. The PDR as a function of viewing angle that was observed by Mathur and colleagues in 30 subjects. Note that the peak PDR is slightly less than one, indicating that the pupil always had an elliptical shape with a vertical major axis in these measurements. Also, the peak value is shifted away from a viewing angle of zero. www.nature.com/scientificreports www.nature.com/scientificreports/ I calculated this component for the current model and found that the result compares quite well with the empirical measurements (Fig. 8). This tilt component of the entrance pupil ellipse can be appreciated in the rendered eye depictions (Fig. 3B), with the top of the pupil tilted towards the nasal field when the eye is viewed from the temporal field, and vice-a-versa. This tilt is almost entirely the result of modeling the vertical displacement of the visual axis from the optical axis; there is a small additional effect of the tilt of the dilated aperture stop.
Error in ellipse fitting to entrance pupil boundary. When observed in near profile, the shape of the human entrance pupil departs from an elliptical form, as is evident in photographs 4,6 . At these extreme angles, points on the border of the aperture stop of the iris are no longer present in the entrance pupil image seen by an observer, either because of vignetting by the corneal limbus, or because rays from those points encounter total internal reflection while passing though the cornea. The entrance pupil is defined instead by the nearest adjacent point from within the aperture stop that does successfully reach the observer, resulting in a distorted appearance of the pupil border.
I evaluated in the current model the extent to which points from the border of the aperture stop are lost in the entrance pupil image as a function of viewing angle (Fig. 9A). Beyond −75 and 65, the model has increasing difficulty finding valid ray paths for points on the border of the aperture stop, consistent with these points encountering total internal reflection at the corneal surface. www.nature.com/scientificreports www.nature.com/scientificreports/ The current model fits the available pupil border points with an ellipse. I examined the performance of this ellipse fit to the pupil as a function of viewing angle by measuring the root mean squared distance of the entrance pupil border points from the fitted ellipse, expressed in millimeters (Fig. 9B). For most viewing angles, the absolute ellipse fit error is modest relative to the 6 mm diameter of the entrance pupil. However, beyond −75 and 65, the accuracy of the ellipse fit rapidly degrades. This indicates that even the pupil border points that are present in the image at these large viewing angles depart in arrangement from an elliptical form.

Discussion
I find that a model eye is able to reproduce empirical measurements of the shape of the entrance pupil with good accuracy. More generally, the current study provides a forward model of the appearance of the pupil as observed by a camera at an arbitrary position.
Similar to the current study, Fedtke and colleagues 7 had the goal of using a ray traced model eye to examine the properties of the entrance pupil as compared to empirical measures. The current study builds upon this prior work in a few respects. First, the anatomical accuracy of the anterior chamber of the model eye is advanced, incorporating an elliptical aperture stop and a tri-axial ellipsoid cornea. Second, the misalignment of the visual and optical axes is modeled. These refinements both improve the match to empirical observations 6 and support the generalization of the model to arbitrary pupil sizes and to viewing angles apart from the horizontal plane.
The model presented here provides the parameters of an ellipse fit to the border of the entrance pupil. I find that an ellipse provides an accurate description of the modeled entrance pupil for viewing angles between −75 and 65 and, within this range, the model agrees well with the empirical measurements of Mathur and colleagues 6 . The current model, however, is not an accurate description of the pupil beyond these angles. When viewed from the far periphery, the entrance pupil departs from an elliptical shape. While the current model reports that some points on the border of the iris stop are not visible to the observer, it does not attempt to construct the distorted pupil boundary. To capture this appearance, the current model would ideally be improved to account for the vignette effect of the corneal limbus, the deviation of the peripheral cornea from an ellipsoidal form 28 , and the departure of the observer from a pinhole camera model.  www.nature.com/scientificreports www.nature.com/scientificreports/ performing ray tracing from points throughout the plane of the aperture stop, instead of only from the stop border.
While the goal in the current study was to account for the average appearance of the entrance pupil of the eye in a population, the model readily accommodates custom biometric measurements for any particular eye under study. Both the profile of the front surface of the cornea and the angle α have an effect upon the resulting entrance pupil appearance. Both have also been proposed to vary with spherical refractive error 10,20 , although this factor explains a relatively small amount of the individual difference. In the current work, I implemented the Tabernero and colleagues model of variation in α related to spherical refractive error 20 . This model component did not perform well in describing the empirical variation in turning point observed by Mathur and colleagues 6 . If α and spherical refractive error are in fact related, a larger set of empirical measurements, including cases of more extreme ametropia, will be needed to define the functional form.
The model is easily extended to account for rotation of the eye. There is an extensive literature in which anatomically-inspired models of the eye are used to inform gaze and pupil tracking systems (e.g. 3,29 ). These applications have a particular focus upon determining the location of the first Purkinje image and the center of the entrance pupil, as these may be used to deduce gaze angle. While not studied here, these values are readily derived from the current model.
A feature of the current model is that it is implemented in open-source MATLAB code, as opposed to within proprietary optical software. The ray tracing operations of the model are conducted using compiled functions. As a result, the routines execute quickly (e.g., <10 msecs to compute the parameters of the entrance pupil ellipse on a typical laptop computer) and can thus form the basis of an "inverse model" that searches across values of eye rotation and stop radius to best fit an observed eye tracking image.
The biometric parameters that define the model were set almost entirely by reference to prior empirical measurements of the central tendency of studied populations. In an effort to best fit the pupil diameter ratio function of Mathur and colleagues 6 , however, two model parameters were tuned "by hand": the maximum, dilated entrance pupil ellipticity was held to a value slightly smaller than previously reported 13 ; and the value for vertical α was set to 2.5° in the emmetropic eye. Because these parameter values are biologically plausible and within the range of prior empirical measurements, I am optimistic that the model will generalize well to novel observations of entrance pupil appearance.