EUROGRAPHICS 2007 / D. Cohen-Or and P. Slavík(Guest Editors)

Volume 26 (2007), Number 3

Interactive Simulation of the Human Eye Depth of Field and
Its Correction by Spectacle Lenses
Masanori Kakimoto1 Tomoaki Tatsukawa1 Yukiteru Mukai1 and Tomoyuki Nishita2
1 SGI Japan, Ltd., Japan

2 The University of Tokyo, Japan

Abstract
This paper describes a fast rendering algorithm for verification of spectacle lens design. Our method simulates
refraction corrections of astigmatism as well as myopia or presbyopia. Refraction and defocus are the main issues
in the simulation. For refraction, our proposed method uses per-vertex basis ray tracing which warps the environment map and produces a real-time refracted image which is subjectively as good as ray tracing. Conventional
defocus simulation was previously done by distribution ray tracing and a real-time solution was impossible. We
introduce the concept of a blur field, which we use to displace every vertex according to its position. The blurring
information is precomputed as a set of field values distributed to voxels which are formed by evenly subdividing theperspective projected space. The field values can be determined by tracing a wavefront from each voxel throughthe lens and the eye, and by evaluating the spread of light at the retina considering the best human accommodation effort. The blur field is stored as texture data and referred to by the vertex shader that displaces each vertex.
With an interactive frame rate, blending the multiple rendering results produces a blurred image comparable to
distribution ray tracing output.
Categories and Subject Descriptors (according to ACM CCS): I.3.3 [Computer Graphics]: Picture/Image Generation

1. Introduction
Reflection and refraction are fundamental themes in computer graphics. Ray tracing is a common solution to thesesubjects. Considering today’s boost in computing power, ray
tracing is still not fast enough to achieve a real-time result for practical applications. Meanwhile, approximations
for reflection and refraction based on environment mappingwith a GPU are used for major applications, especially forgames. For applications of lens system design, such approximation methods are not used since they lack the accuracy
that is required for the design processes. Ray tracing is stilla major tool to simulate refraction for this purpose.
Spectacle lens designers, like most designers of other industrial products, use CAD applications to form the surface shape of their product. Their design task is a repeated processes of shape improvement and verification. The
turnaround speed is getting better because of more sophisticated interactive tools enabling intuitive design of surface
shapes. Unlike other industrial design areas, the turnaroundof lens design is limited because ray tracing is in its verification stage.
c The Eurographics Association and Blackwell Publishing 2007. Published by Blackwell°
Publishing, 9600 Garsington Road, Oxford OX4 2DQ, UK and 350 Main Street, Malden,
MA 02148, USA.

A naked eye view
(30.6fps)

A view through a
A view through aconcave lens (20.1fps) convex lens (20.1fps)

Figure 1: Simulation results of a patient model with bothmyopia and presbyopia.

In spectacle lens verification this is especially noticeable
because distribution ray tracing, which is several times or
an order of magnitude slower than regular ray tracing, ismandatory for the simulation of a defocused image. Our

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

interactive display of human vision corrected by spectaclelenses dramatically changes the design process.
This paper presents a method of human eyesight simulation targeted especially for refraction correction with spectacle lenses. We introduce an interactive algorithm for reflection and refraction, which we call Per-vertex Ray Tracing.
Recently, there are several other research works pursuingreal-time plausible reflection and/or refraction. An advantage of ours is that the image is accurate at every vertex on
the reflective or refractive object. Our algorithm enables theuser to display 3D scenes through the lens in real-time witha sufficient refraction accuracy to simulate the distortion andcolor aberration of the lens.
Defocusing is a more critical issue for spectacle lens simulation. Many interactive applications employ a defocus effect that is often referred to as depth of field (DOF). Much
work has been done for defocus effect with various camera models or image based models. However, these models lack human eye characteristics and are not applicablefor DOF correction by spectacle lenses. There have been asmall number of research activities on such correction. No
previous work has been found which successfully realized
real-time or interactive simulation of the human visual DOF
phenomenon and its correction by lenses.
Two major obstacles are the performance of distribution
ray tracing and the performance of human eye simulation.
We solved the first issue with Per-vertex Ray Tracing andthe second issue by precomputing the result of human eye
and spectacle lens simulation. The result is stored as spatially distributed, volumetric information, which we call a
blur field. In the rendering runtime each vertex of the objects in the scene is displaced according to its location in theblur field.
The blur field concept alone can be combined with distribution ray tracing to simulate eyesight or even other DOF
models, provided that the primitives are polygon based and aroutine is defined to displace each vertex instead of jittering
the viewpoint.
This paper is organized as follows: In Section 2 we review
some previous work on interactive refraction and human visual simulation. Section 3 describes a real-time refraction
technique which we call per-vertex ray tracing. Section 4 introduces our human visual simulation model and the conceptof blur field. Section 5 summarizes results and performanceand Section 6 gives our conclusion with some discussion.

2. Previous Work
Since our proposed method uses real-time refraction, we first
introduce related research on this topic. The second half of
this section presents some previous work on the simulationof the human visual system and refraction correction.

2.1. Real-time Rendering with Refraction
Environment mapping [BN76] is a common technique toapproximate reflection. Since it can easily be hardwareaccelerated [HS93] [VF94], environment mapping became
widely accepted for reflection effects. Real-time refraction mapping techniques were derived [LKM01] [Wym05]
but they are limited to refraction for distant objects,
which is a critical limitation of environment mapping.
Hakura [HSL01] used a layered map to support local objects.
Ohbuchi [Ohb03] realized a real-time method using refraction rays. Both methods map the environment objects to individual image planes and the accuracy is limited. To improve
the accuracy of the local environment mapping, distance impostors [SKALP05] were introduced. Using the depth channel of the environment map, they approximated the ray hit
points by iterative computing. Their method is fast and accurate for large-planar objects but the accuracy decreases forcomplex environment objects.
Unlike other real-time refraction methods, ours computes
per vertex ray-scene intersections to obtain accurate 3D coordinates. Although the performance penalty is higher than
distance impostors, it still runs in 30fps for practical scenes.
More importantly for distortion evaluation, except for a caseof occlusion discrepancy, which virtually does not occur inlens simulation, our method guarantees the refraction accuracy at each mesh point of the refractive surface.

2.2. Simulation of Correction of Human Visual Systems
A number of researchers have proposed DOF models [PC81] [Shi94] [KMH95]. Distribution ray tracing [CPC84] is one of the most popular DOF techniques foroff-line rendering. Accumulation buffer [HA90] is a hardware mechanism which realizes real-time DOF.
Some research work has been done on DOF taking thehuman eye model into account. Santamaria et al. [SAB87]
measured the human eye as a lens system to acquire a point
spread function and tried to reconstruct the image on theretina. Mostafawy [MKL97] introduced a thick lens modelof a human eye and simulated refractive surgery effects.
Loos et al. [LSS98] used wavefront tracing to evaluate thehuman eye accommodation and visualized refractive correction by progressive lenses. Both of the above use distribution ray tracing and real-time processing is impossible.
Barsky [Bar04] advocated the concept of Vision Realistic
Rendering and introduced a DOF technique using aberration
information measured from a patient’s optical system.
Our proposed method produces a similar type of result
to Loos’ or Barsky’s but the performance is enhanced byseveral orders of magnitude, for example 20fps vs 0.002fps.
c The Eurographics Association and Blackwell Publishing 2007.
°

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

Per-vertex ray tracing

Cube mapping

Figure 3: An example of reflection in a car door mirror model. In the central two images a wireframe model of the mirror isoverlaid to show the vertex positions. Each red sphere has been located at the 3D intersection of each reflection vector and thescene, and is reflected in the mirror. The intersections are reflected exactly at the corresponding mirror vertices in the imagesof per-vertex ray tracing while they are not in the cube mapping images.
3. Per-vertex Ray Tracing
The basic idea is to warp the environment map so that eachvertex reflects the accurate scene intersection point of theextension of its reflection vector.
3.1. Algorithm for Each Frame
The processing step in each frame is as follows: The item
numbers correspond to those appearing in Figure 2.
1. Compute reflection (refraction) vector at each vertex of
the reflective (refractive) objects sharing an environmentmap.
2. Compute a virtual viewpoint, form a view volume, andrender an intermediate image to be used as an environment map.
3. For each vertex of the reflective (refractive) objects:a. Trace the extension of the reflection (refraction) vector to find the nearest intersection with the scene.
b. Find the line that connects the intersecting point andthe virtual viewpoint.
c. Find the point where the line intersects with the projection plane for the environment map.
d. Assign the coordinates of the point in the environment
map to the texture coordinates of the relevant vertex.
4. Render the scene and the reflective (refractive) objects
from the original viewpoint.

(3.b) (3.c)(3.d)

reflective object

(2)

(3.a)

virtual viewpoint
vertex(4)
objectsin the scene

original
viewpoint(1)
projection (nearclipping) plane forthe environment map

Figure 2: Processing steps of the per-vertex ray tracing.

c The Eurographics Association and Blackwell Publishing 2007.
°

In step 2 a least-squares approach is used to quickly findthe virtual viewpoint that best represents the reflection (refraction) rays from the vertices.
With this algorithm the accurate intersection point is reflected at each vertex of the reflective object. The only exception is occlusion inconsistency potentially caused by thedifference between the viewing ray from each vertex and thatfrom the virtual viewpoint. Its possibility should, however,be minimized by the least-squares method. In addition, thepossibility is even more negligible for spectacle lens simulation, in which case the distributions of the vertex positions
and the refraction vector directions are moderate and limited.
For refractive lens objects, the process is identical to thatin Figure 2 except that each refraction vector is computedwith double refraction at the back and front faces of the lens.

3.2. Results of reflection and refraction
Figure 3 presents a sample result of reflection using Pervertex Ray Tracing with a cube mapping reference. While
local reflections yield apparent errors in the cube mapping
example, the per-vertex ray tracing results in an exact matching at each vertex.
A rendering example of double-refraction through a lensis shown in Figure 4.
This is also an example of color aberration of the lens.
The right image of Figure 4 is a close-up of the edge ofthe lens, where the color difference is apparent. The color
aberration simulation is accomplished by rendering the lens
three times, each with one of the three separate sets of texture coordinates computed by iterating step 3 described in
Section 3.1 using different refraction indices (and thus different sets of refraction rays) for red, green, or blue. Eachof the above lens rendering tasks is carried out with appropriate color mask settings to make a color-shifted refraction
image over the lens surface. Note that all three tasks sharethe first pass result, i.e., the environment map, and thus theextra computation for color aberration is minimized.

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses
ophthal- view
view
lens
mologic
volume volume geoparameters (2nd pass) (1st pass) metry

virtualviewpointinterior originalscene viewgraph point

textureFigure 4: Color aberration exhibited near the edge of a convex lens. RGB’s are separately rendered for the lens surface
using different refractive indices (R:1.60, G:1.61, B:1.62).
4. Defocus Using Blur Fields
The basic idea for defocus is to displace the objects on a
per vertex basis, according to the vertex location relative tothe lens and the eye. We introduce a concept of blur field,
a precomputed 3D spatial distribution of shape information
which defines how a specific place should be blurred.
We define the blur field in the normalized device coordinate space, by subdividing in x, y, and z directions. Eachsubdivided box, which we call a voxel, contains a blur fieldvalue as a result of the precomputation.
We assume that the scene objects consist of polygon meshmodels. Since the defocus operation is carried out per vertex,
each model needs to be tessellated with sufficient granularity
to produce a desired quality.
This section describes both the precomputation process(Sections 4.1-4.4) and the displacement operation (Section
4.5) in the rendering loop. The outline of the process is described in Figure 5 and Figure 6.
4.1. Blur Field
4.1.1. Use of the Normalized Device Coordinate System
The view volume for the refraction map is used as the space
for the blur field. After the projection transformation andperspective division, the view volume is mapped to the normalized device coordinate system (NDC). The blur field isdefined in the NDC and each sampled value is stored ina voxel. The voxels are formed by evenly subdividing thewhole NDC space.
An advantage of using NDC is its compatibility with theperspective projection. Places closer to the viewpoint in theworld coordinate space are more precisely subdivided thanplaces farther from the viewpoint. This is a similar concept
to that of perspective shadow maps [SD02].
Another major advantage is rendering performance. Since
the vertex shader always transforms all vertices to the clipping space, which is very close to the NDC, the determination of the voxel in which each vertex resides will require a
simple arithmetic rather than a 4 × 4 matrix transformation.
Section 4.5 describes the details.

wavefront
tracing (precomp.)
without lens

wavefront
tracing (precomp.)
with lens

1st
pass

vertex shader input

2nd
pass
n

n

blur field(corrected)

blur field(naked eye)Figure 5: The whole process of the method. The n is thenumber of jittered images to blend.
voxelswave source
wavefront
tracingz
defocusevaluation
at the retina
view volume forthe 1st pass(a) Precomputation

1st pass result
mapped to thelensblur fieldvalue

zvirtualviewpoint(b) 1st pass renderingzview volume for the
2nd pass with anaked eye blur field(c) 2nd pass rendering

Figure 6: The outline of the defocus processing flow. As
a result of the precomputation, an ellipsoid shape information is stored for each voxel. The 1st pass consists of several separate rendering passes applying, in turn, a displacement vector within the ellipsoid to every vertex in the scene.
Their results are blended onto the accumulation or auxiliary buffer. The blurred image is mapped to the lens surface
in the 2nd pass using the per-vertex ray tracing technique.
Note that we assume that the lens position relative to theviewpoint does not change during the rendering loop.
4.1.2. Voxel information
Each voxel contains region information which defines thelimit of the displacement of the vertices residing in the voxel.
The region is determined by an ellipsoid whose size and orientation are computed by wavefront tracing. The wave emanates from a point in the voxel through the spectacle lensand the eye lens, and finally onto the retina. Details are presented in Sections 4.3 and 4.4.
4.2. Finding Light Ray Paths
To trace a wavefront, we need to know the path along which
the ray from the center of the voxel travels through the spectacle lens into the eye, as indicated in Figure 6(a). It is not atrivial problem to find the light paths often described as Fermat’s Principle. We solved this using an image processing
approach.
c The Eurographics Association and Blackwell Publishing 2007.
°

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

First, a 2D (x-y) array of voxels at a certain depth (z) is
rendered as a 2D array of rectangles through the lens with
per-vertex ray tracing, using the same viewpoint position as
in rendering the current frame (2nd pass). The background is
set to black and the color of each rectangle is made unique.
A simple image processing of the rendered result canidentify the pixels corresponding to the centers of the voxels.
Each of these pixels corresponds to a point on the projection
plane, in the world coordinate system. A ray from the viewpoint to such a point should travel through the lens beingrefracted twice, and then reach to the center of the voxel.
This double-refracted ray is a very good approximation of
the path obeying Fermat’s Principle, connecting the viewpoint and the center of the voxel. In addition, in most cases
for lenses, the double-refracted ray strikes through all voxels of different depths but at the same x-y position. Thus, therefracted light ray path can be shared among such voxels asline up along the z direction. To compute a blur field value
for a voxel, a wave along this path is traced back from anappropriate passing point in the voxel. In this paper, we callthis starting point for each voxel a wave source.
4.3. Wavefront Tracing
Wavefront tracing [Kne64] [Sta72] is a general tool tosimulate traveling light rays. It was applied to produce
caustics from reflective surfaces [MH92] [Elb94]. Loos etal. [LSS98] used it to evaluate the human eye accommodation for each pixel. We used wavefront tracing more directly
than Loos et al. in order not only to evaluate accommodation but also to make a defocus spread pattern on the retina
by deriving new wavefront operations, which are describedin Section 4.4.
The wavefront at a point can be represented by the direction of the wave movement, i.e., a normal vector N, two
principal directions e1 and e2 which are perpendicular witheach other, and corresponding principal curvatures κ1 andκ2 . Note that a curvature is negative when the center of curvature is in the opposite direction of the wave movement.
In spreading or transferring by a distance d, a wavefront
is simply transformed as follows.
κ1
κ2
κ10 =
, κ20 =.
(1)
1 − d κ1
1 − d κ2
Here, the vectors N, e1 , and e2 do not change.
The transformation for refraction is derived from Snell’s
law and is fully described in the references [Kne64] [Sta72].
We present only the result here as Mitchell [MH92] andLoos [LSS98] described. The direction vector N is transformed as N0 = µ N + γ N(s) , where N(s) is the unit normal
vector of the refractive surface and µ = nn12 , γ = −µ cos ϕ +cos ϕ 0 . n1 and n2 are the indices of refraction of the media
before and behind the surface, respectively. ϕ and ϕ 0 denotethe angle of incidence and the angle of refraction, respectively. When a unit vector u ∝ N × N(s) is given, which isc The Eurographics Association and Blackwell Publishing 2007.
°

tangent both to the incident wavefront and to the surface,curvatures with respect to u can be obtained as follows. Using another unit vector v = N × u and the angle ω betweenu and the first principal direction e1 , i.e.,u = cos ω e1 − sin ω e2 , v = sin ω e1 + cos ω e2 ,(2)Euler’s Formulas

κu = κ1 cos2 ω + κ2 sin2 ωκv = κ1 sin2 ω + κ2 cos2 ω(3)

κuv = (κ1 − κ2 ) cos ω sin ωdetermine the directional curvatures κu , κv and the torsion
κuv on the incident wave. By refraction, curvatures of the
wavefront are transformed as:(s)κu0 = µκu + γκucos2 ϕ1(s)+ γκu
cos2 ϕ 0
cos2 ϕ 0
cos ϕ
(s) 1
0
= µκuv
κuv+ γκuv,cos ϕ 0cos ϕ 0κv0 = µκv

(s)

(s)

(4)

(s)

where κu , κv , and κuv are, respectively, the directional
curvatures and the torsion of the refractive surface at the incident point. They can be computed using Equation 3 similarly to the wavefront. To represent the wavefront after thetransformation, we need to find the principal curvatures anddirections by inverting Equation 3.

κ1 = κu cos2 ω + 2κuv cos ω sin ω + κv sin2 ωκ2 = κu sin2 ω − 2κuv cos ω sin ω + κv cos2 ω2κuv
tan 2ω =.
κu − κv(5)

Table 1 presents the usage of wavefront tracing operations in
the proposed method, showing which of the operations takes
place in which turn during the computation of the blur field
for each voxel.

Table 1: The process of precomputation for each voxel.
These operations take place along the double-refracted ray
from the wave source of each voxel to the viewpoint.
place or media
wave source
air
front
spectacle
glass
lens
back
air
eye

pupil
retinavoxel

operation
wavefront creationtransferrefraction
transferrefraction
transfer
accommodation
vergence
conoid sectioning
scaling to NDC

equation(1)(3) (4) (5)(1)(3) (4) (5)(1)(13) (15)(3) (7) (5)(8)(9)

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

y4.4. Eye Model
To simulate vision with astigmatic eyes, Loos et al. [LSS98]
used distribution ray tracing. We continued to use wavefront
tracing for the region within the eye by introducing new
wavefront operations which are a transformation by a thinlens and a sectioning. Our model transforms the wavefront
at the eye lens to form a revolved conoid shape and sections
the shape at the retina to acquire correct defocus or a point
spread shape which will be stored as a value of the blur field.
We assume that the eye as a lens system is a thin toric lens,i.e., having orientation dependent refractive powers which
are uniform over the lens surface. From an ophthalmological perspective, our eye model can treat regular or obliqueastigmatism in addition to myopia, hyperopia, and presbyopia.

e1ω1 κ 2
x

κ1

y

y

e2
x ω1

(e)

z

(a) incident wave

x

P1

ω1′

z

κ 2′
z

P2

(b) refractive powerof the eyeκ 1′(c) refracted wave

Figure 7: Wavefront transformation by an astigmatic eyewith principal refractive powers P1 and P2 .
In optical science, the basic vergence formula is

4.4.1. Wavefront Transformation by a Thin Toric Lens

V + P = V 0,

Wavefront transformation in the eye is a complex phenomenon. In addition to approximating the eye as a thin toriclens, we assumed that the incident light ray direction is perpendicular to the lens. This assumption makes sense since
the eyeball rotates to the region of interest so that the ray
converges at the central retinal fovea. This subsection introduces a simple transformation of the incident wavefrontand derives a refracted wavefront which will converge at theretina.

where V is the vergence of the incident light wave and canbe denoted as V = nκ by definition, with the refractive index
n of the media of the incident side and the curvature κ at theincident point, P is the refractive power of the lens, and V 0
is the vergence of the outgoing light. We extended this equation for directional curvature and derived formulas for the
wavefront transformation by an eye as a toric lens, yielding

After refractions by the spectacle lens and a translation to
the surface of the cornea, we obtain a wavefront with principal curvatures κ1 and κ2 and the principal directions e1
and e2 , respectively. As illustrated in Figure 7, we define aneye lens coordinate system as a right-hand coordinate system with its origin at the incident point and the z axis alignedwith the wavefront normal vector. The x axis is chosen to beboth horizontal and tangent to the wavefront and the lens.
For convenience the direction of positive x is determinedsuch that the y axis points upward. Since any tangent vector
on the incident wavefront will be a common tangent vector
with the lens, we can choose the x axis as a common tangent,
which was denoted as u in Section 4.3, in order to represent
the wavefront for the transformation by the thin toric lens.
Let ω1 be the angle between the x axis and the first principal direction e1 , and we can apply Equation 3, substituting
ω = −ω1 , to obtain the directional curvatures κx , κy , andκxy .
For the wavefront transformation at the surface of theeye lens, we employ the vergence formula instead of Snell’sLaw. Here, we introduce directional refractive powers of thetoric lens as a set of input parameters of the astigmatic eye.
Given the principal refractive powers P1 and P2 , and theangle between the x axis and the first principal direction(e)ω1 , we find the directional refractive powers Px (in the x-z
plane), Py (in the y-z plane), and Pxy by applying Equation 3(e)substituting κ1 = P1 , κ2 = P2 , and ω = −ω1 .

κx0 = µκx +

(6)

Py
Pxy
Px
0
, κy0 = µκy + , κxy, (7)= µκxy +n2n2n2where µ = nn12 , n1 is the refractive index of the air, and n2 isthat of vitreous humor of the eye.
We find the final wavefront to converge at the retina by
0 , obusing Equation 5, with κu = κx0 , κv = κy0 , and κuv = κxytaining the principal curvatures κ10 and κ20 , and the angle ofthe first principal direction ω10 = −ω .
4.4.2. Evaluation of Light Spread at the Retina
In general the light wave as derived above, when refracted
by a toric lens, does not converge to a point. After incoming through a circular pupil of a certain size, the transferring light forms a conic-like shape known as Sturm’s conoid,
yielding two focal lines being perpendicular with each otheras illustrated in Figure 8. Two main issues here are (1) the
analysis of Sturm’s Conoid and (2) the human eye accommodation. This subsection describes how to find, given a pupil
and a static light wave forming a conoid, the shape of thelight spread projected onto the retina, whereas the next subsection discusses the accommodation.
For simpler observation of the Sturm’s conoid, we ignore
for a while the rotation of the principal directions about zaxis, i.e., we assume ω10 = 0. To simulate the spread of the
projected light spot, we take a cross-section of the conoid
with a plane perpendicular to z-axis. By setting z as a constant parameter ze for the distance between the cornea andthe retina, we obtain an implicit equation of an ellipsoidc The Eurographics Association and Blackwell Publishing 2007.
°

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

pupil

y
A

ω1′
x
O

rp

κ 2′
κ 1′

circle of posterior
anterior
focal line
least
focal
confusion
line

b
B

C

P

z

1
a
κ 1′

ze

1
κ 2′

zc

Figure 9 illustrates how sh is computed. A small horizontal length hr at the retina should be mapped to a length in
the normalized device coordinate system (NDC) hNDC . The
scaling factor sh can be computed by hNDC /hr . The verticalscaling factor sv is computed similarly.

interval of Sturm

refracted wave

Figure 8: Sturm’s conoid formed by a toric lens.
the wave source of a voxel

side of theview volume

plane perpendicular to the viewing direction

viewing
direction
spectacle lensRpsin θ ′
(D p − Dl )ρ
hlp =
sin θ

θ′

h p = hl + hlp

θ

pupil
(viewpoint)
retina

Dp

Dl
zr

hr

(a length at the retina)

ρ≡

hr
zr

hNDC =

hl = Dl ρ

sh =

hp

Note that every variable in the right-hand side of Equation 9 varies by the individual voxel. The result of the multiplication of the three 2 × 2 matrices is stored in each voxel.
4.4.3. Accommodation of the Eye
It is the common assumption that the eye tries to accommodate so as to align the circle of least confusion on the retina.
Let zr be the distance between the thin lens of our eye modeland the retina, and the above condition is satisfied by letting
a = −b in Equation 8. Thus we obtain

Rp

hNDC
hr

Figure 9: Scaling factor to convert a horizontal length hr at
the retina to h p at the wave source of a voxel and to hNDC in
the normalized device coordinate. The refraction at the spectacle lens is taken into account. The lens is approximated asa thin lens in this case.
22
S(x, y) = ax2 + by2 = 1, wherea = r p (1 − κ10 ze ),Let sh and sv be the scaling factors from the retina to
the NDC for horizontal and vertical directions respectively,
and we obtain the final blur field value Bvx ,vy ,vz at a voxel(vx , vy , vz ) as:
µ
¶µ
¶µ
¶
a 0
s
0
cos ω10 − sin ω10
Bvx ,vy ,vz = h. (9)sin ω10cos ω10
0 b
0 sv

zr =

2
.
κ10 + κ20

We can expand the above formula by using Equations 3, 7,
4, and 5 in turn and using the fact n1 = 1.0, to yieldzr =b = r p (1 − κ20 ze ),(8)and r p is the radius of the pupil. The above a and b can be derived from the proportional relationship among the sides ofseveral of the right-angled triangles along z-axis (e.g. 4AOC
and 4BPC in Figure 8). The parameters κ10 and κ20 should bewhat the system had determined taking into account the accommodation as will be explained in Subsection 4.4.3.
At the retina, where ze = zr , the ellipse shape described bya, b, and ω10 represents the point spread of the light source.
When the light wave from the wave source of a voxel istraced and eventually results in an ellipsoid at the retina, wetreat the shape as an approximated defocusing pattern which
any point light source in the voxel should produce, and register the information for the voxel. Instead of registering a, b,
and ω10 , an equivalent 2 × 2 matrix which should transforma unit circle into the oriented ellipsoid is stored as the voxel
value. This precomputed matrix for each voxel is the entityof the blur field.
To determine the scaling factor of the matrix, we need toconsider a geometrical relationship between the retina andthe normalized device coordinate system (NDC), the scaleof which is equivalent to that of the voxel space.
c The Eurographics Association and Blackwell Publishing 2007.
°

(10)

n2
,
κ̄ + P̄

(11)

where κ̄ = κ1 +2 κ2 and P̄ = x 2 y . Optically κ̄ is the averagecurvature of the incident wavefront to the eye and P̄ is theaverage refractive power of the eye.
P +P

Let P̄min be the standard non-astigmatic refractive powerof the eye in minimum (relaxed) accommodation, and let ∆P
be the accommodation refractive power limited by the maximum accommodation ability ∆max , i.e., ∆P ∈ [0, ∆max ], thenP̄ = P̄min + ∆P.
From Equations 11 and 12, we obtain
µ
µ
¶¶
n2
∆P = max 0, min− κ̄ − P̄min , ∆max
.
zr

(12)

(13)

The meaning of the above accommodation formula is as follows. When the mean curvature of the incident light κ̄ , which
is usually negative, is smaller (−κ̄ is greater), i.e., the lightwave comes from a closer location to the eye, the accommodation effort ∆P becomes greater, and vice versa. A standard
value for ∆max is presented in Table 2 but it declines throughout life.
Let A be the user-defined astigmatism, i.e., A = P1,min −P2,min , where P1,min and P2,min are the principal refractive

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

powers of the eye in minimum accommodation and it is natP+Pural to assume P̄min = 1,min 2 2,min , thenAA
, P2,min = P̄min − ,(14)22and, by definition, we obtain the principal refractive powersof the eye taking account of accommodation,
P1,min = P̄min +P1 = P1,min + ∆P, P2 = P2,min + ∆P.

(15)

The wavefront tracer will receive the above two values with(e)
the user-defined parameter ω1 in order to carry on tracing
as was described in the second half of Subsection 4.4.1.
Table 2 shows the standard parameters of the Gullstrand
and Helmholtz’s schematic eye [GvH09] which we used in
the implementation. Table 3 indicates which of the parameters should be controlled to simulate which of the the humanvision symptoms.
Table 2: Values of the Gullstrand and Helmholtz’s schematiceye. The refractive power of maximum accommodation
should be reduced according to the age.
parameter
minimalaccommodation
refractive maximal
power
accommodation
ability
distance between front
cornea and retina
refractive air
vitreous humor
index

symbol

value

P̄min

58.64 (D)∆max

11.93 (D)zr

22.785 (mm)

n1
n2

1.000
1.336

vector is a 2D vector from the center of a unit circle to a randomly generated point within the circle. In the vertex shader,the seed is transformed into a 2D displacement vector using the blur field value Bvx ,vy ,vz (Equation 9) which has been
looked up with the NDC coordinates of the vertex. The NDC
can be obtained by dividing the x, y, and z by w coordinate.
The shader adds the 2D displacement vector to the NDC x
and y values and gets them back to the homogeneous coordinates just before exiting.
The blur field Bvx ,vy ,vz is supplied to the vertex shader as
either a 3D texture or a 2D texture tiled with sub-textures.
A typical size of a blur field is 32×32 voxels for x and y
directions and 128 for z (depth) direction. Each field value
is a 2 × 2 matrix and can be represented as a 4-component
texel.
A vertex shader code for the displacement using OpenGL
Shading Language [RK06] is shown below.
uniform sampler3D blurField;
uniform vec2 seed;void main(void) {// emulate fixed function lighting......
......
gl_Position = ftransform();// displace the vertex in NDC
const float w = gl_Position.w;
vec3 ndc = vec3(gl_Position.x/w,
gl_Position.y/w, gl_Position.z/w);
vec3 texCoord = vec3(ndc.x/2.0 + 0.5,
ndc.y/2.0 + 0.5, ndc.z/2.0 + 0.5);
vec4 mat = texture3D(blurField, texCoord);
vec3 disp = vec3(mat.x*seed.x + mat.y*seed.y,mat.z*seed.x + mat.w*seed.y, 0.0);
ndc += disp;
gl_Position = vec4(ndc.x*w, ndc.y*w, ndc.z*w, w);}

The unit D is diopter (m−1 ).

Table 3: Simulation targets and their control parameters.
target symptom
refractive
myopiaaxial
refractive
hyperopia
axial
presbyopia
regularastigmatism
oblique
daylight or night visionparameter control
P̄min > 58.64 (D)zr > 22.785 (mm)P̄min < 58.64 (D)zr < 22.785 (mm)∆max < 11.93 (D)(e)
A > 0, ω1 = 0(e)
A > 0, ω1 6= 00.75 < r p < 4 (mm)4.5. GPU Implementation
In the rendering main loop, the scene is iteratively drawn as
many times as the user-defined number of samplings for defocus and the results are blended to produce a blurred image.
The trick to apply the blur field to each vertex is a simple
process.
For each sampling, a displacement seed vector is generated and fed into the vertex shader. A displacement seed

5. Results and Performance
The performance is affected by the overall number of polygons in the scene. Especially, the mesh complexity of therefractive lens has an impact because the rays are traced on a
per-vertex basis. We used lens models of about 1,000-vertex
meshed polygons, which yielded images practically as accurate as ray tracing. All performance numbers were measured
using 3.2GHz Pentium4 CPU with NVIDIA GPU GeForce
7800GTX graphics, with the 1024 × 768 pixel output image
size. Please see the attached movies for the real-time captured interactive operations.
Rendering double-refraction through a lens object
achieves 30-40fps for scenes of approximately 10,000 polygons. With color aberration, the throughput declines to 15fps(Figure 4).
The depth of field processing consists of multiple rendering and the performance is affected by the number of samplings. We used ten-sampling blurred image, which is sufficient for lens design verification in our experience. Withoutlens refraction, the frame rate for a blurred image in Figure 1
c The Eurographics Association and Blackwell Publishing 2007.
°

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses

Figure 10: A virtual eyeglass simulator using polarized
stereo eyeglasses and projectors with polarized filters.
left is around 30fps. Rendering of refraction correction images by a spectacle lens took 7fps to over 20fps dependingon the scene complexity, window size, and the number ofsamplings. Each of the middle and right images in Figure 1
took 20.1fps at 10-sampling with 1024 × 768 window. For
the Figure 1 images, the pupil diameter was set to be 2cm
intentionally to exaggerate the defocusing results. The blur
field values are in proportion to the pupil size r p , which is
obvious from Equations 8 and 9.
Figure 10 shows snapshots of a binocular, stereoscopic
visualization system using this method. The rendering takes
just twice the normal application since the dual images are
rendered independently with different input conditions.

Progressive lensNaked eye reference

Figure 11: Correction by a prototype progressive lens.

Figure 11 is a result of the verification of a progressive
lens under design. While the upper half of the lens is pretty
good at adjusting the myopia, a major portion of the lower
half except the center fails to correct the presbyopia, showing the difficulty of the design of progressive lenses.
An astigmatism example is presented in Figure 12 witha partial blur field voxel information visualized as a set of
ellipsoids. An observation of a whole blur field is presented
in Figure 13.
The precomputation to generate a blur field using the
wavefront tracing takes 300 milliseconds for a 8 × 8 × 128voxel blur field, 1 second for 16 × 16 × 128, and 7 seconds
for a 32 × 32 × 256 blur field.
6. Conclusion
We have proposed a fast depth of field rendering method taking into account the properties of the human eyesight. Our
method has the following characteristics.
• Near real-time frame rate with refraction and defocusing• A precise human-eye model capable of simulating various eyesight syndrome and their correction by spectaclelenses• Precomputed spatially distributed defocus information
stored as a blur field
The blur field is a key and novel idea for the fast rendering.
It is easy to implement on programmable GPUs and the data
set is simple and compact. It could be utilized by other applications which require flexible, space dependent displacec The Eurographics Association and Blackwell Publishing 2007.
°(a) Without DOF(b) DOF with a blur field

Figure 12: An astigmatic naked eye simulation result. The
newspaper object is placed so close (6.1cm at the lower-leftcorners and 10.7cm at the upper-right) that the eye is unable
to accommodate. A partial blur field (voxels of almost thesame depth as the newspaper) is visualized as an ellipsoidfor each voxel. Note that every blur field ellipsoid is also
blurred exactly by its own spread range. Also note that theellipse shapes differ by their distances from the viewpoint,
exhibiting the sections of the Sturm’s conoid.
ment of objects in the scene. Blur fields can be combined
even with distribution ray tracing if the programmers add avertex displacement routine.
A limitation of the blur field technique is that the objects
need be tessellated into meshed polygons with some sufficient granularity comparable to the blur field voxel size. A
solution to this might be to use the blur field in the fragment
shader, which is a future work.
Although our method is currently applied in a spectaclelens design company, applications to the image content generation industry are also possible, for example, as a special
effect for realistic human vision.

M. Kakimoto et al. / Interactive Simulation of the Human Eye Depth of Field and Its Correction by Spectacle Lenses[KMH95] KOLB C., M ITCHELL D., H ANRAHAN P.: A
Realistic Camera Model for Computer Graphics. In Proc.
SIGGRAPH ’95 (1995), pp. 317–324.
[Kne64] K NEISLY J. A.: Local curvature of wavefronts
in an optical system. Journal of the Optical Society ofAmerica 54, 2 (1964), 229–235.
[LKM01] L INDHOLM E., K ILGARD M. J., M ORETON
H.: A user programmable vertex engine. In Proc. SIGGRAPH 2001 (2001), pp. 149–158.

Figure 13: An observation of an 8 × 8 × 128 voxel blur field(top-right) and its applied result (bottom-right). Each of thesmall red circles represents a blur field value of a voxel. The
circles are so closely formed near the viewpoint that they
appear to be a block. The size of a circle is equal to the sizeof the blurring at the place, as shown in Figure 12. Strong
myopia and presbyopia are simulated, and the focus is only
on the object where the blur field values are zeros.
Acknowledgments
We thank Naoya Hirai, Masaki Sawai, and Tomofumi Teratani of SGI Japan for making the interior scenes and thedemo movies, and Takeshi Haga for the advice on the shaderprogram implementation. Thanks also go to Gordon Jolley
for proofreading the drafts.
References[Bar04] BARSKY B. A.: Vision-realistic rendering: simulation of the scanned foveal image from wavefront data of
human subjects. In Proc. APGV ’04 (2004), pp. 73–81.

[LSS98] L OOS J., S LUSALLEK P., S EIDEL H.-P.: Using
wavefront tracing for the visualization and optimizationof progressive lenses. Computer Graphics Forum (Proc.
Eurographics 1998) 17, 3 (1998), 255–263.
[MH92] M ITCHELL D., H ANRAHAN P.: Illuminationfrom curved reflectors. In Proc. SIGGRAPH 1992 (1992),
pp. 283–291.
[MKL97]

M OSTAFAWY S., K ERMANI O., L U H.:
Virtual Eye: Retinal Image
Visualization of the Human Eye. IEEE CG&A 17,
1 (January-February 1997), 8–12.
BATSCHOWSKI[Ohb03] O HBUCHI E.: A real-time refraction renderer forvolume objects using a polygon-rendering scheme. In
Proc. Computer Graphics International 2003 (CGI’03)(2003), pp. 190–195.
[PC81] P OTMESIL M., C HAKRAVARTY I.: A lens andaperture camera model for synthetic image generation. In
Proc. SIGGRAPH ’81 (1981), pp. 297–305.
[RK06] ROST R., K ESSENICH J. M.: OpenGL Shading
Language, 2nd Edition. Addison-Wesley, 2006.

[BN76] B LINN J. F., N EWELL M. E.: Texture and reflection in computer generated images. Communications of
the ACM 19, 10 (1976), 542–547.

[SAB87] S ANTAMARIA J., A RTAL P., B ESCOS J.: Determination of the point-spread function of human eyes
using a hybrid optical-digital method. Optical Society ofAmerica A 4, 6 (1987), 1109–1114.

[CPC84] C OOK R. L., P ORTER T., C ARPENTER L.: Distributed Ray Tracing. In Proc. SIGGRAPH ’84 (1984),
pp. 137–146.

[SD02] S TAMMINGER M., D RETTAKIS G.: Perspective
shadow maps. In Proc. SIGGRAPH ’02 (2002), pp. 557–
562.

[Elb94] E LBER G.: Low cost illumination computation
using an approximation of light wavefronts. In Proc. SIGGRAPH ’94 (1994), pp. 335–342.

[Shi94] S HINYA M.: Post-filtering for depth of field simulation with ray distribution buffer. In Graphics Interface’94 (1994), pp. 59–66.

[GvH09] G ULLSTRAND A., VON H ELMHOLTZ H.:
Handbuch der Physiologischen Optik. 1909.

[SKALP05] S ZIRMAY-K ALOS L., A SZODI B., L AZANYII., P REMECZ M.: Approximate Ray-Tracing on the GPUwith Distance Impostors. Computer Graphics Forum(Proc. Eurographics 2005) 24, 3 (2005), 685–704.

[HA90] H AEBERLI P., A KELEY K.: The accumulationbuffer: hardware support for high-quality rendering. In
Proc. SIGGRAPH ’90 (1990), pp. 309–318.
[HS93] H AEBERLI P., S EGAL M.: Texture mapping as A
fundamental drawing primitive. In Fourth Eurographics
Workshop on Rendering (1993), pp. 259–266.
[HSL01] H AKURA Z. S., S NYDER J. M., L ENGYEL
J. E.: Parameterized environment maps. In SI3D ’01:
Proceedings of the 2001 symposium on Interactive 3Dgraphics (2001), pp. 203–208.

[Sta72] S TAVROUDIS O. N.: The Optics of Rays, Wavefronts, and Caustics. Academic Press, 1972.
[VF94] VOORHIES D., F ORAN J.: Reflection vector shading hardware. In Proc. SIGGRAPH ’94 (1994), pp. 163–
166.
[Wym05] W YMAN C.: An approximate image-space approach for interactive refraction. In Proc. SIGGRAPH
2005 (2005), pp. 1050–1053.
c The Eurographics Association and Blackwell Publishing 2007.
°