BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR93-1395
ENTRY:: 1994-03-17
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Towards Accurate and Efficient Volume Rendering
AUTHOR:: Novins, Kevin L.
DATE:: October 1993
PAGES:: 181
COPYRIGHT:: Kevin Lawrence Novins 1994 - All Rights Reserved
ABSTRACT::
This thesis is concerned with improvements to algorithms for volume 
rendering; a technique that provides scientists with the means for visual 
exploration of three-dimensional data. Despite its numerous successes, and 
its increasing use within the scientific community, state-of-the-art volume 
rendering algorithms have many shortcomings. Difficulties include: ensuring 
the accuracy of the rendered images, producing images with modest 
computational resources, and rendering the diverse types of data that are 
currently being produced.

The work in this thesis was motivated by the demands of an ongoing
visualization project in four-dimensional cardiac visualization. We present 
solutions to some key problems in ensuring accuracy and in producing 
algorithms that can scale to handle large datasets. Although the theoretical 
work in this thesis applies to arbitrary data topologies, our implementations 
have assumed that the data is defined by sample points on a regular 
rectilinear grid.

In the area of accuracy, we focus on the error that is introduced during 
volume projection. This phase of the volume rendering process involves the 
evaluation of the emission-absorption volume rendering line integral. This 
thesis presents four techniques for controlled precision volume integration. 
These schema depart from existing approaches in that they provide error 
bounds along with the solutions they generate. In each case,  the error 
analysis leads to an algorithm for evaluating the integral to any specified 
tolerance.

Our investigations into efficiency issues have resulted in two advances. 
First, an adaptive error bracketing scheme is presented that builds on 
the controlled precision volume integration methods. Using adaptive error 
bracketing, the solution for a viewing ray is continually refined until a 
user-specified error tolerance is met. The algorithm allows processing of the 
data without imposing a strict front-to-back or back-to-front evaluation 
order. Second, a suite of tools are presented that can be used to efficiently 
compute perspective projections of volume data. These include a paging 
strategy that is useful when a dataset is too large to fit into RAM memory 
and a ray splitting technique for adaptive supersampling. The latter 
technique ensures that all data features contribute to the final image while 
avoiding overcomputation in regions close to the eyepoint.
END:: CORNELLCS//TR93-1395
BODY::
Towards Accurate and Efficient
Volume Rendering
Kevin Lawrence Novins
Ph.D Thesis
TR 93?i395
October 1993
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
TOWAR?DS ACCURATE AND EF'F'ICIENT
VOLUME RENDERING
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Kevin Lawrence Novins
January 1994
o Kevin Lawrence Novins 1994
ALL RIGHTS RESERVED
TOWARDS ACCURATE AND EFFICIENT
VOLUME RENDERING
kevin Lawrence Novins, Ph.D.
Cornell University 1994
This thesis is concerned with improvements to algorithms for volume rendering; a
technique that provides scientists with the means for visual exploration of three-
dimensional data. Despite its numerous successes, and its increasing use within
the scientific community, state-of-the-art volume rendering algorithms have many
shortcomings. Difficulties include: ensuring the accuracy of the rendered images,
producing images with modest computational resources, and rendering the diverse
types of data that are currently being produced.
The work in this thesis was motivated by the demands of an ongoing visu-
alization project in four-dimensional cardiac visualization. We present solutions
to some key problems in ensuring accuracy and in producing algorithms that can
scale to handle large datasets. Although the theoretical work in this thesis applies
to arbitrary data topologies, our implementations have assumed that the data is
defined by sample points on a regular rectilinear grid.
In the area of accuracy, we focus on the error that is introduced during volume
projection. This phase of the volume rendering process involves the evaluation of
the emission-absorption volume rendering line integral. This thesis presents four
techniques for contrvtled prec?s?on volurne integration. These schema depart from
existing approaches in that they provide error bounds along with the solutions
they generate. In each case, the error analysis leads to an algorithm for evaluating
the integral to any specified tolerance.
Our investigations into efficiency issues have resulted in two advances. First, an
adaptive error bracketing scheme is presented that builds on the controlled precision
volume integration methods. Using adaptive error bracketing, the solution for
a viewing ray is continually refined until a user-specified error tolerance is met.
The algorithm allows processing of the data without imposing a strict front-to-
back or back-to-front evaluation order. Second, a suite of tools are presented
that can be used to efficiently compute perspective projections of volume data.
These include a paging strategy that is useful when a dataset is too large to fit
into R?AM memory and a ray splitting technique for adaptive supersampling. The
latter technique ensures that all data features contribute to the final image while
avoiding overcomputation in regions close to the eyepoint.
BIOGRAPHICAL SKETCH
Kevin L. Novins was born on December 2, 1963. He attended public schools in
Westchester County, New York and graduated from Ardsley High School in 1981.
He studied computer science at Harvard University, and received the Bachelor
of Arts degree in 1985. After graduation from college he performed research in
medical imaging at the North Carolina Memorial Hospital. He entered Cornell
University in the fall of 1987 and received the Master of Science degree in 1990
and the Doctor of Philosophy degree in 1994.
iii
ACKNOWLED GEMENTS
First, I would like to thank the Chair of my thesis committee, Donald P Greenberg,
for making this work possible. He helped me to set my goals and worked to
establish and maintain the research collaborations that were necessary for carrying
them through. I appreciated his guidance, enthusiasm, dedication to quality, and
unwavering support. I also appreciated his teaching me how to present my research
properly. I am glad for all the times that we worked closely together. Don has
created a wonderful environment at the Program of Computer Graphics and I am
grateful that he gave me the opportunity to work and study there.
Much of this research could not have been performed without the help of Jim
Arvo, who was an advisor and a friend throughout the process. It was Jim who
first stimulated my interest in performing an error analysis of volume rendering.
He assisted in high-level planning, in working out implementation issues, and in
settling some of the final production details. Jim also helped me to sort out
my thoughts when I got frustrated, and provided encouragement during research
slowdowns. I am grateful to Jim for his contributions to this thesis and for all that
he has taught me.
iv
Another friend and advisor was Fran?ois Sillion. Francois was instrumental in
teaching me how to perform graduate-level research, and helped to get me started
on this dissertation. I would like to thank him for helping me to become more of
a scientist and less of a hacker.
I would also like to thank David Salesin, who was a colleague in developing the
adaptive error bracketing algorithm. His enthusiasm was key to the success of that
project, and he helped a great deal to move me along towards graduation.
The collaborative research project in cardiac visualization that I discuss in
Chapter 6 was the impetus for this research. It was a pleasure to work with
Dr. Richard Devereux of the Cornell Medical Center, Dr. Riccardo Pini of the
University of Firenze, and Steven Leavitt of Hewlett-Packard. They convinced me
of the importance of the project and provided me with the education that I needed
in the medical and engineering aspects of ultrasonic imaging.
I would like to thank Daniel Huttenlocher and Anthony Reeves for serving on
my thesis committee.
Several of the results presented here have been published in the form of research
papers. I am grateful to my co-authors for their contributions, particularly since
sections of text and figures have been incorporated into this thesis. Specifically, I
would like to thank: flancois Sillion and Donald Greenberg, my co-authors for "An
Efficient Method for Volume Rendering Using Perspective Projection" [NSG9O];
Jim Arvo, my co-author for "Controlled-Precision Volume Rendering" [NA92];
and Jim Arvo and David Salesin, my co-authors for "Adaptive Error Bracketing
for Controlled Precision Volume Rendering" [NA5921.
v
Program of Computer Graphics staff members provided a great deal of support
for me and for this research. Hurf Sheldon was particularly helpful in ensuring
that I always worked in the best possible environment. Jim Ferwerda generously
provided invaluable assistance with many aspects, including the preparation of
slides and videotapes. I would also like to thank Ellen French, Fran Brown, Janet
Brown, and Emil Ghinger for their help.
The CT pelvis dataset used in several of the figures was provided by the De-
partment of Radiation Oncology at the North Carolina Memorial Hospital. The
cardiac dataset of Figures 2.7 and 6.8 were provided by the Cornell Medical Cen-
ter. The electron density dataset of Figure 5.26 was provided by Duncan McRee
of the Scripps Clinic and obtained as part of the Chapel Hill volume rendering test
dataset. The photograph in Figure 3.5 of Andrew James Monks was provided by
Michael Monks and Marietta Christie.
This work was partially supported by the NSF grants "Interactive Computer
Graphics Input and Display Techniques" (CCR-8617880) and "Visualization for
Supercomputing: A Graphics Workstation Approach" (AC-8715478). It was also
supported by the NSF/DARPA Science and Technology Center for Computer
Graphics and Scientific Visualization (ASC-8920219). The work was performed
on equipment donated by the Hewlett-Packard Corporation.
I would like to thank my friends and lab-mates for sharing wonderful times,
and for providing essential support during some not-so-wonderful times.
Finally, I would like to thank my parents, Rochelle and Jay Novins, and my
brother and sister, Doug and Kara, for their constant love and support.
vi
TABLE OF CONTENTS
1 Introduction
2 Volume Illumination Models
2.1			Emission-Absorption Model
2.2			Other			Volume Illumination Models
2.2.1			Maximum-Intensity Projection
2.2.2			Attenuation Projection
2.2.3			Low-Albedo Models.
2.2.4			Global Illumination
2.3			Discussion
3 Volume Rendering Algorithms
3.1			Mapping.
3.1.1			Classification
3.1.2			Absorption			.
3.1.3			Emission
3.2			Reconstruction
3.3			Projection Algorithms
3.3.1 View Transformations
3.3.2 Partitioning the Integral
3.3.3 Ray Casting Algoritlims
3.3.4 Data-Space Algorithms
3.3.5 Frequency Domain Algorithms
3.4			Summary			. . .
4 Controlled-Precision Volume Rendering
4.1 Controlled-Precision Volume Integration
4.1.1 The Integrand and its Derivatives
4.1.2 Piecewise Constant Bounds
4.1.3 Trapezoid and Simpson's Rule
vii
6
7
10
11
12
13
16
19
21
23
24
29
30
35
41
41
42
45
47
55
58
60
62
64
67
70
4.2
4.3
4.4
4.1.4 Power Series
4.1.5 Discussion and Results
Controlled-Precision Ray Casting
Integration over pixel area .
Summary and Future Directions
5 Efficiency Issues
5.1 Efficient Perspective Projection
5.1.1			An Optimal Paging Strategy
5.1.2			Sampling Issues			. . .
5.1.3			Results
5.1.4			Discussion . . .			. . .
5.2			Efficient Controlled-Precision Ray Casting
5.2.1			Overview . .			. .
5.2.2			Selection
5.2.3			Refinement . .			. .
5.2.4			Updating
5.2.5			Initializing the BSP tree
5.2.6			Results
5.2.7			Discussion . .			. .
5.3 Summary and Future Directions
6 Imaging for Cardiology
6.1 The goals of the project
6.2 Data Collection .
6.3 Current Limitations
6.4 Visualization of the data
6.5 Discussion
7 Conclusion and Future Directions
A Nomenclature
B Power Series Convergence Proof
Bibliography
viii
74
82
86
90
92
94
95
98
106
115
118
122
123
124
128
130
130
133
137
140
145
145
147
153
154
159
166
169
172
175
LIST OF TABLES
3.1 Classification via thresholding.
4.1 Integration rules and their restrictions.
A.1 Notation for interval arithmetic
ix
25
85
171
LIST OF FIGURES
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
4.1
4.2
4.3
4.4
The volume illumination problem.
Image of the pelvis dataset computed using the emission-absorption
model
Image of pelvis dataset computed using maximum-intensity projec-
tion
Image of the pelvis dataset computed using attenuation projection.
The low-albedo approximation.
Computing the low-albedo model using the emission-absorption
model
The effects of shadowing. . . .			. .			. . .
Geometry of Figure 2.7.
A volume rendering pipeline.
Geometry of the phase angle.
Angles for Phong illumination.
The effect of gradient shading.
A sample reconstruction problem.
Ideal reconstruction
Nearest-neighbor interpolation.
Linear interpolation. . .
Gaussian interpolation. . .
Two-pass 2D rotation. .
The generic footprint table. . .
Splatting
Decomposition of a cell into subcells.
The Fourier slice theorem.
The volume rendering integrand.
The piecewise constant approximation.
The trapezoid rule
Pseudo-code for the trapezoid rule algorithm.
x
7
10
12
13
15
16
17
18
22
31
34
36
37
38
39
40
41
49
51
53
54
56
66
71
73
74
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
5.12
5.13
5.14
5.15
5.16
5.17
5.18
5.19
5.20
5.21
5.22
5.23
5.24
5.25
5.26
6.1
6.2
6.3
Simpson's rule.			. . .			. .			. .
The power series approximation.
Pseudo-code for computing 1n (a, b).			. .			.
Pseudo-code for computing a remainder term.
Pseudo-code for the complete power series. .
Performance on 10000 test integrals. .
The Snyder and Barr ray/grid intersection algorithm.
Pseudo-code for controlled-precision ray casting.
Controlled-precision ray casting. . . .
Comparison of orthographic and perspective views			.
Computing perspective by predistorting the volume.
Quadratic versus cubic growth rate.
Ray-oriented versus data-oriented processing.
Pseudo-code for initializing the slab queue.			. . .			.
Pseudo-code for processing the slab queue.
Slab processing order. .			. . .			.
Undersampling in perspective			projection.
Uniform supersampling. . . 			. .			. .			.			.
Ray splitting
Adaptive supersampling			. .			.			. .			.			.
Orthographic versus perspective views of CT data.
Uniform vs. adaptive supersampling			. .			. .
Frames from a data fly-through.
Processing order for depth of field.
Lens effects. .			. .			.
The main loop of the adaptive error bracketing algorithm.
Selecting a segment using the span tree.			. . .			.
The selection process for the adaptive error bracketing algorithm.
Updating a segment using the span tree.
The update process for the adaptive error bracketing algorithm.
The front-to-back algorithm. .
The geometry of the dataset in Figure 5.24.			. . .
Adaptive error bracketing: The grid dataset.
Adaptive error bracketing: The CT pelvis dataset. .
Adaptive error bracketing: The SOD enzyme dataset.
Gathering 2D ultrasound data. .
The ultrasound system. . . .			. .
Coordinate axes for spherical and cylindrical grids.
xi
75
77
79
81
82
83
87
89
91
96
97
98
99
100
102
105
108
109
111
112
116
117
119
120
121
125
128
129
131
132
134
136
142
143
144
149
151
152
6.4
6.5
6.6
6.7
6.8
Gathering 3D ultrasound data.
Resampling a cylindrical grid onto a rectilinear grid.
Splat shapes for cylindrical data.
Ray-casting onto a cylindrical-coordinate grid
Comparison of approaches to rendering cylindrical-coordinate data.
155
158
160
161
163
xii
Chapter 1
Introduction
Over the past ten years, great advances have been made in scanning techniques that
produce volume data, such as computed tomography, magnetic resonance imaging
and ultrasonic imaging. These developments have brought large quantities of data
into the hands of physicians and scientists, much of which can only be interpreted
visually. This thesis is concerned with providing tools for the visual exploration of
their data. To be useful, the tools must be intuitive, provide accurate renditions,
and present the results in a timely fashion with modest computational resources.
The work presented in this thesis relates to a subfield of scientific visualization
called volume rendering. Volume rendering algorithms can be distinguished from
other approaches in that images are created directly from a three-dimensional array
of raw data, without fitting a surface as an intermediate step. For the display of
data that is not surface-oriented, such as the output of a tomographic scanner or
the result of a fluid simulation of gaseous phenomena, volume visualization has
2
been the method of choice. It has been successfully applied to data from medicine,
chemistry, astronomy, meteorology, and many other disciplines.
Despite these numerous successes, and its increasing use within the scientific
community, state-of-the art volume rendering algorithms have many shortcomings
as a visualization tool. Developments are needed in the following areas:
Accuracy. Current algorithms employ many approximations for which the impact
on the final image is not known. Error is introduced during data acquisition, recon-
struction from point samples, projection onto screen-space, and display. Volume
rendering algorithms cannot currently be relied upon in fields where high accuracy
is essential.
Efficiency. Volume rendering requires a great deal of calculation, disk storage and
random-access memory. The algorithms become impractical when used with high
resolution dynamic datasets. Efficiency improvements are essential for ensuring
that visualization can be performed outside of special-purpose computing centers.
Ideally, algorithms should be fast enough and computing resources powerful enough
to provide scientists with low-cost tools for interactive manipulation of the images.
Arbitrary data formats. Volume rendering algorithms are currently inappropri-
ate for a large body of scientific data. The technique has only gained acceptance for
visualizing scalar-valued three-dimensional fields. Additionally, most algorithms
impose the limitation that data samples lie on a regular rectilinear grid.
The above difficulties become clear when applying volume rendering techniques
to "real-world" scientific data. The work in this thesis was motivated by problems
that arose as part of an ongoing project in cardiac imaging. In collaboration with
3
researchers at the the Cornell Medical Center, ESAOTE Biomedica, and Hewlett-
Packard, we are developing techniques for dynamic, three-dimensional imaging of
cardiac structures via ultrasound. This method has great potential since it is both
non-invasive and low-cost, making it feasible for diagnosis in a doctor's office. The
project places several new demands on existing volume rendering methods. Since
the application is medical, we need to rely on the images for accurate diagnosis.
Since time effectively adds a fourth dimension to the already large dataset, the
storage demands are excessive, particularly for real-time playback. Finally, the
data sample points do not lie on a regular rectilinear grid.
For this thesis, we concentrated on the issues of accuracy and efficiency as they
relate to volume rendering. While our solutions are not a panacea, they do provide
important advances in each of these areas. Although the theoretical work in this
thesis applies to arbitrary data topologies, the implementations have assumed that
the data is defined by sample points on a rectilinear grid.
In the area of accuracy, we focus on the error that is introduced during vol-
ume projection. This phase of the volume rendering process is concerned with the
evaluation of the emission-absorption volume rendering integral. We present tech-
niques for evaluating this integral to any specified precision. Tools for computing
the integral along the intersection of a viewing ray with a voxel are extended to
consider the intersection of a ray with multiple voxels. Having solved the latter
problem, we can compute point samples of the illumination function across the
image plane to any specified precision. Our approaches represent a departure from
existing volume rendering algorithms in that they provide quantitative assessment
4
and allow the user control over the quality of the solutions.
Our investigations into efficiency have resulted in two advances. The first is an
adaptive error bracketing scheme that builds on the controlled-precision volume
projection methods. Using adaptive error bracketing, the solution for a single
viewing ray is continually refined until a user-specified error tolerance is met. This
avoids overcomputation. For the first time, this algorithm allows processing of
the data without imposing a strict front-to-back or back-to-front evaluation order.
Our studies have shown that adaptive error bracketing saves computational effort
when working with smooth datasets or when computing low-accuracy images.
We also present efficient volume rendering algorithms for computing perspec-
tive projections. Perspective projections are desirable because they provide a depth
cue that is easy to comprehend. Previous research in volume rendering has con-
centrated on the less intuitive parallel projection model since it has been both
easier to compute and provides a regular sampling of the scene. Our results ad-
dress both these problems. We present an optimal paging strategy that is useful
when computing perspective projections of datasets that are too large to fit into
RAM memory. An adaptive supersampling solution to the problem of perspective
ray divergence is also presented. This ensures that all features in the back of the
dataset are accounted for, while avoiding overcomputation in regions towards the
front of the dataset.
The organization of the thesis is as follows. Chapter 2 presents a review of
theoretical models for volume illumination. This includes a discussion of the
emission-absorption model, upon which the work in this thesis is based. Chap-
5
ter 3 contains an overview of the volume rendering process and describes existing
algorithms. Chapter 4 presents our results in approximating the volume projection
problem to any specified precision. Chapter 5 is devoted to efficient computation
of volume-rendered images. It includes techniques for perspective projection and
the adaptive error bracketing algorithm for controlled-precision volume projection.
In Chapter 6 we present the results-to-date of our cardiac imaging project.
The thesis concludes in Chapter 7 with a summary and a discussion of future
research directions.
Chapter 2
Volume Illumination Models
The process of rendering is essentially a simulation of light interacting with an
environment. In rendering for visualization, the goal is to communicate scientific
information. The computational nature of the simulation affords the programmer
the luxury of choosing the rules that will govern the simulation, thus allowing a
great deal of control in devising ways to meet this goal. The existence of a set of
rules is important, since the abstract nature of the visualization means that errors
can only be measured in terms of deviation from the assumed model. This is in con-
trast to methods for realistic image synthesis where the rules are the laws of physics
governing the transfer of light energy, and errors can be measured in terms of devi-
ation from experimental data. In this chapter, we present the em?ss?on-abso?ption
illumination model that has been widely adopted [Luc92,NH92,VPR?+92i and is
the basis for the work in this thesis. We also discuss other illumination models
that have been used in the context of volume rendering.
6
7
The basic illumination problem is depicted in Figure 2.1. Consider a ray of
light passing through a medium between points a and b. We would like to know
the intensity of light at wavelength A that is incident on point a in the direction
of ray % passing through the material between a and b. We denote this quantity
i\(a, b). Color images are generated by computing i\(a, b) at several wavelengths.
In the following discussion, we assume that I?(a, b) is computed independently for
each wavelength and will drop the A subscripts for notational convenience.
Figure 2.1: The volume illumination problem.
2.1 Emission-Absorption Model
A simple model of volume illumination was originally developed as a model of
radiative transfer [Cha6O,Sob63]. It is based on a physical model of an emitting
and absorbing medium, and does not consider scattering effects. The work in this
thesis is based on this model.
The model requires that two nonnegative functions are defined across the vol-
ume. The absorption function, c(x, y, z), determines the fraction of light that is
absorbed by the medium along a differential ray segment passing through (x, y, z).
8
It is normally assumed to be wavelength-independent. The emittance function,
denoted e(x, y, z) determines the wavelength-dependent emission of the medium
along a differential ray segment passing through (w, ij, z).
In the following discussion, the variables t, u, and v will be distance parameters
along the ray ab. The functions c and ? will correspondingly be reparameterized
by distance along the ray. By fixing b, the endpoint of the ray, we can write the
differential equation for I(t, b), the intensity of light at point t along the ray ab
due to material between t and b:
--H )I(t, b) --H?(J) I([, b) t ?(t). (2.1)
This equation states that the change in intensity, as we move towards the eyepoint,
across a differential ray length dt is equal to the light emitted minus the fraction
of incoming light that is absorbed. This differential equation can be solved using
an integrating factor [RB74j. We sketch the solution process below.
We begin by applying the integrating factor,
e?f??(v)dv,
to both sides of the differential equation (2.1),
ftl(t, b) ?(t) --H ?(t) K(t) I(t, b) --H --H?(t) ?(t).
Since,
--H--H --H?(t) e?fbta?(v)dv
we may rewrite Equation (2.2),
d
dtI(t? b) ?(t) ? K'(t) I(t, b)			--H?(t) ?(t)
(2.2)
9
d
dt [I(t, b) ?(?)1			--He(t) ?(t).
Integrating both sides of the equation from b to t yields,
I(t, b) k(t)			--H			?(u) ?(u)du + 1b,
where 1b is a constant of integration. We continue by isolating I(t, b) on the left-
hand side,
Since
we may continue,
I(t, b)			--H			?(u) ?(u)du ?--H? (t) + 1b k--H' (t).
k 1(t) --H e?ftbQ(v)dv,
I(t, b)			--H			E(u) e?f??a(v) dv e?ftbQ(v) dvdu + 1b e?ftbQ(v) dv
--H --H Ibt c(u) e?j??Q(v) dv du + 1b e?ftbQ(v) dv
--H fb ?(u) e?f7Q(v) dv du + 1b e??bQ(v) dv
(2.3)
(2.4)
If we solve Equation (2.4) for t a, we obtain the expression for I(a, b):
I(a, b)			$?e--Hi;?a(v) dv e(u) du + 1b e b?(v)dv			(2.5)
The constant 1b is the intensity of backlighting in the ray direction impinging at
point b. Equation 2.5 summarizes the volume illumination model used in this
thesis. We will normally assume that there is no backlighting. In this case, we
may drop the second term and write:
I(a, b) --H $6e?f?Q(v) dv ?(u) du.			(2.6)
10
The right side of Equation (2.6) shall henceforth be referred to as the volnmc
?endering integral. An image computed using the emission-absorption model is
given in Figure 2.2.
Figure 2.2: Image of the pelvis dataset computed using the emission-absorption
model.
2.2 Other Volume Illumination Models
The emission-absorption model is one of many volume illumination models that
appear in the literature. Its main attractions are its physical basis, simplicity and
its generality. In this section, we review other models in current use and compare
their performance with that of the emission-absorption model.
11
2.2.1 Maximum-Intensity Projection
The simplest projection model for volume rendering is that of mar?murn-intensity
projection [Sab88,VPR+921. Using maximum-intensity projection, the intensity
of the ray between points a and b is proportional to the highest density, p(t),
encountered between points a and b. Thus
I(a,b) k mbax p(t),
where k is an arbitrary scaling constant. In this form, maximum-intensity projec-
tion produces only gray-scale images. However, color is often added to ameliorate
the fact that the technique does not provide any depth information via occlusion.
For example, each pixel can be colored based on the distance of the peak value
from the eyepoint [Sab88].
Maximum-intensity projection is useful in application domains where it is im-
portant to locate extremal values in the data. Such a situation occurs in medical
angiography, where the goal is to detect arterial blockage. The blood stream is
injected with a contrast agent, thus making it the densest object in the body, and
a tomographic scan is performed. By employing maximum-intensity projection in
rendering the resulting data, a physician can ensure that the view of the blood
stream is not occluded by less dense structures.
This model's major drawback is that it presents a narrow view of the volume;
discarding all information except that which involves the most dense material in
the dataset. Thus, the technique remains of interest in only a few highly specialized
fields, and even in these, it is used as only one part of a suite of visualization modes.
12
It cannot be considered a truly general visualization technique. An image of the
pelvis dataset computed using maximum-intensity projection is given in Figure 2.3.
Figure 2.3: Image of pelvis dataset computed using maximum-intensity projection.
2.2.2 Attenuation Projection
Using attenuation projection [Lev92], only the absorption of the media is modeled.
We can obtain an expression for I(a, b) by starting with the emission-absorption
model of Equation (2.5) and setting the emission coefficient, 6(u), to zero yielding,
I(a,b)			Ibe?f?Q(V)dV			(2.7)
Note that in order to use this model, the backlight intensity, 1b must be nonzero.
Equation (2.7) shows that using the attenuation projection model, only the sum of
the density along the ray is considered; the order in which the material is encoun-
tered is irrelevant. As with maximum-intensity projection, there is no occlusion.
13
The effect of attenuation projection is similar to that of an X-ray. In fact, simu-
lation of X-rays is one important application of this model [SNC9OJ. An example
of attenuation projection is given in Figure 2.4. Because of its lack of shading
and depth cues, this model is not completely satisfactory for communicating 3D
information. Nonetheless, it holds attraction due to its easy evaluation relative to
those models that include emissive effects.
Figure 2.4: Image of the pelvis dataset computed using attenuation projection.
2.2.3 Low-Albedo Models
Another popular class of volume illumination models is based on the concept of the
volume media as a varying density field of low-albedo particles that reflect light
[Bli82,NV84,Sab88]. Albedo, or reflectivity, is the fraction of light impinging on a
particle that is reflected. The fundamental assumption of this class of models is
that since particles are of low albedo, all light paths from light source to eyepoint
14
which include more than one reflection off a particle can be ignored. If a particle
reflects at most a fraction, ?, of the incoming light, a light path which bounces off
n particles reflects at most ?? of the incoming light. Therefore, mathematically,
the low-albedo approximation states that
0, n>2.
When using these models, absorption is typically related to particle density, de-
noted p(u) for any given ray. Emission is determined by computing a one-bounce
light reflection from the light source, off the particle, to the eyepoint. The light
reflection model varies from simple phase functions that depend only on the angle
between the light source and the eyepoint [Bli82,KV84] to more complex models
that also consider the orientation of the particles [Lev88,DCH88,Sab88]. The re-
flectance function is denoted ?(u). Choice of ? will be covered in more detail in
Chapter 3. Sometimes, shadowing of the light source by material in the volume is
accounted for by a function, V(u), which computes the optical distance from point
u along the ray to the light source [KV84].
The low-albedo expression for I(a, b) is
I(a, b)			$be?LP(v)dv V(u) ?(u) p(u) du.
This expression can be adapted for use with more than one primary light source
by summing expressions for I(a, b) over all light sources.
The generality of the emission-absorption model allows it to capture the effects
of the above low-albedo model, even though the emission-absorption model does
15
not explicitly consider light reflection. If we define,
p(u)			(2.8)
V(u)?A(u)p(u),			(2.9)
we see that the low-albedo model is equivalent to the emission-absorption model.
Equation 2.9 effectively moves the calculation of the light reflection into a pre-
process that sets the emission coefficients at each point in the volume. Figure
2.5 illustrates how light is reflected off particles using the low-albedo model. Fig-
ure 2.6 shows how this is simulated using the emission-absorption model. In this
fashion, the effects of any illumination model that allows only a fixed number of
reflections in between the light source and the eyepoint may be captured using the
emission-absorption model.
Figure 2.5: The low-albedo methods model a single reflection of light off of the
particles in the medium.
16
(a)
Figure 2.6: Computing the low-albedo model using the emission-absorption model.
During a preprocess (a) the light reflection model is computed, yielding emission
values throughout the volume. Afterwards, (b) the emission-absorption model is
run using these emission values.
2.2.4 Global Illumination
Although the emission-absorption and low-albedo models are physically-based,
they express only a small range of natural lighting phenomena such as occlu-
sion, simple shading, and optical depth. The next step in increased realism is to
account for multiple scattering in the illumination model. Such models are called
global illumination models because the illumination at any given volume element
depends on the illumination of every other volume element.
Global illumination has been given little attention by researchers in scientific
visualization; however, the problem has been studied from the point of view of
realistic image synthesis. Najiya and Von Herzen found that multiple scattering
was essential for realistic rendering of clouds, since water droplets are a high-albedo
medium [NV84]. Researchers in radiosity have used volume techniques to produce
17
photorealistic scenes that include smoke and fog [RTS7,Rus8Si. Unfortunately, the
extensive computational and memory requirements of these algorithrns make them
practical for only very small volumes (currently on the order of 32 x 32 x 32 voxels)
and clearly not sufficient for rendering modern scientific datasets.
(a)			(b)
Figure 2.7: Comparison of image (a) with shadowing and (b) without shadowing.
An additional reason for the lack of interest in global illumination by the visu-
alization community is the lack of clear evidence that more realistic lighting effects
lead to improved image comprehension. In order to get a sense of the value of
increased realism in volume-rendered images, we performed a comparison between
images created with and without shadowing of the light source by the volume. We
rendered several animation sequences of a 64 x 64 x 64 dataset of cardiac ultrasound
data, and tested several light source positions. Some rotation sequences were com-
puted with the light source fixed relative to the object and others were computed
with the light source fixed relative to the eyepoint. Frames computed with and
18
without shadowing from the latter type of sequence is shown in Figure 2.7. There
are two lights, one 45 degrees to the left and one 45 degrees to the right of the
viewer. They are also inclined 45 degrees above the viewer. The relative positions
of viewer, volume environment, and lights are shown in Figure 2.8.
(a)			(b)
Figure 2.8: Geometry of Figure 2.7: (a) side view; (b) top view.
In all our experiments, we found that shadows did not improve comprehension
of the scene, and frequently degraded comprehension. Two specific problems were:
1. Shadows obscured regions of the dataset. The effects of texture and shading
in attenuated areas were minimized, thereby making perception of surface
orientation difficult.
2. Unfamiliarity with the object meant that shadows appeared in unexpected
places. Viewers often incorrectly interpreted shadowed areas as `?dirty vox-
els".
Other researchers have incorporated shadows in volume rendering systems [Lev89,
SN89,ASN92]. Their results demonstrate that shadows are effective when the
19
shadows are cast by approximately convex objects onto orthogonal planes. In
such cases, shadows give cues about the positioning of the volume objects in space
relative to the planes as well as providing silhouettes of the object. None of these
researchers demonstrate the use of shadows for highly concave scenes, such as the
ventricle of the heart shown in Figure 2.7. More study is necessary, but our results
show that improved realism does not guarantee improved comprehension.
Recent work by Nrueger [Nru9O] suggests that scattering terms provided by
global illumination models may be useful for datafeature enhancement. Using this
scheme, scattering is used not to produce photorealistic images but instead to give
the user textural cues about statistical properties of the volume data.
2.3 Discussion
In summary, the emission-absorption model was chosen for the work in this thesis
because
o+ it is physically-based and can model familiar lighting effects such as occlusion
shading, and optical depth;
o+ it is general enough to subsume two of the more popular volume illumination
models: attenuation projection and the low-albedo approximations;
o+ it gives the user great flexibility by allowing emission and absorption to be
decoupled; and
o+ its simple form is conducive to error analysis.
20
The only major drawback of the model is that it cannot capture all the nuances
of global illumination. However, we hope that understanding the accuracy and
efficiency issues for the simpler model will act as a first step towards evaluating
more complex models. In the meantime, more work is necessary to establish that
the added complexity and expense of the global illumination is worthwhile in terms
of visualization applications.
Having presented the model of illumination upon which this work is based, we
continue in Chapter 3 by exploring existing algorithms for evaluating the model
and rendering images.
Chapter 3
Volume Rendering Algorithms
In the previous chapter, we described ideal models of light transfer in an envi-
ronment. In this chapter we shall see how these models are incorporated into a
complete volume rendering algorithm. This chapter is organized around the three
basic steps of volume rendering: rnapp?ng, reconstruction, and projection. These
steps can be envisioned as part of a rendering pipeline, as illustrated in Figure 3.1.
The input to this pipeline is an array of raw data values at a discrete number of
sample points in space. During mapping, raw data samples are interpreted to pro-
duce samples of the parameters that are required by the illumination model. For
the emission-absorption model, these parameters will be ? and e. The reconstruc-
tion phase generates 3D scalar fields from the point samples of the illumination
model parameters. Finally, during the projection phase, the illumination model is
evaluated and a 2D image of the 3D data is produced. This image is then sampled
for display on limited-resolution devices.
21
22
Raw Data Samples
Mapping
cc Samples
nL
Reconstmction
T
cc
E Samples
JL
Reconstmction
t
E
Projection
T
Rendered Image
Figure 3.1: Pipeline showing the three basic stages of volume rendering.
23
Although the phases of mapping, reconstruction and projection are common
to all algorithms, it is important to keep in in mind that this pipeline is only a
conceptualization of the process. No volume rendering algorithm operates in such
discrete stages, and many perform these steps in a different order.
3.1 Mapping
The input to a volume rendering algorithm is usually produced by scientific scan-
ning equipment or by numerical simulations. The types of data that are input
to volume rendering algorithms are extraordinarily diverse. For example, three-
dimensional medical data available with current scanning technology include:
o+ Radiographic density (computed tomography)
o+ Relaxation constants (magnetic resonance imaging)
o+ Interfaces between media with different acoustical impedances (ultrasound)
o+ Fluid flow (Doppler imaging)
o+ Metabolic activity (positron emission tomography)
These data values generally do not correspond to the input parameters required
by volume illumination models. Even when there is some similarity, the direct
mapping does not necessarily communicate the desired information in the dataset.
For example, radiographic density measured by a CT scanner could be used as the
density constant p for particle-based low-albedo models (Section 2.2.3). However,
if the physician only wishes to examine bone anatomy, the dense soft tissue would
24
obstruct the view. Instead, a mapping which only assigns density to bone matter
would be appropriate. Thus, mapping of raw data to rendering parameters is
another type of abstraction that is used in volume rendering.
The emission-absorption model requires that raw data be interpreted and map-
ped to emission (e) and absorption (c) values. The generality of the emission-
absorption model allows a great deal of flexibility in this mapping. The most
common approach has been to first classify each voxel according to what materials
it contains, and then to assign emission and absorption values based on the desired
visual properties of each material.
3.1.1 Classification
The problem of classification is essentially one of object recognition. Raw data
values are examined and each voxel is marked with the material or materials that
it contains. It is closely related to the field of computer vision, where video pixels
are examined in order to recognize objects in a projected scene. The techniques
used for classification depend heavily on the type of input data. Nevertheless
three general strategies have proved to be most popular for scientific visualization.
The simplest of these methods is called thresholding. This technique divides
the raw data values into discrete ranges that correspond to different materials.
The boundaries of these ranges are called thresholds. In order to classify the
data values, ?, into N materials, we sort the materials in order of increasing raw
data range, .1... mN. We then define N + 1 thresholds to.. 1N and perform the
classification indicated in Table 3.1.
25
Table 3.1: Classification via thresholding.
Raw Data Range Material Type
to ? ? ti
--H? ? ? t2
tN-i ? < tN			mN
For each material m?, thresholding makes a binary decision at each sample
point: "Does this sample point represent material m??" This is called b?nar?
classification. There are two problems with binary classification. The first is that
it creates sharp boundaries between different materials, resulting in a misleading
blocky appearance in the final image if the material boundaries are not aligned
with the sample points. Furthermore, it is susceptible to classification errors due
to the "partial volume effect".
The partial volume effect occurs in measured data, where samples do not rep-
resent values at infinitesimal points in space. Rather, they represent the values
across small finite volumes. A sample volume may therefore contain more than
one type of material. In such a domain, a binary classification scheme is bound
to make inappropriate "all-or-none" decisions. The decision it makes is not even
guaranteed to correspond to one of the materials present in the sample volume. For
example, if materials with non-adjacent raw data ranges occupy the same sample
volume, the resulting average data value might indicate to the threshold scheme
26
that the sample represents a material whose data range falls between that of the
two component materials.
An improvement on thresholding, presented by Levoy [Lev88j, compensates for
the partial volume effect. Let ? be the the raw data value output by the scanner
if a voxel is completely occupied by material rn?. His method assumes that
o+ at most two materials share a sample volume; and
o+ if m? and m? share a sample volume then there does not exist another material
mk with a ? value between ? and ?
The technique classifies a raw data value ? by first determining the materials m?
and m? such that
The volume represented by ? is then classified as being
--H? x 100
percent material Tnj, with the remainder classified as being rnj. These perfect
conditions tend not to occur in real datasets that have more than two types of ma-
terials. However, with data such as medical CT scans, the conditions are approxi-
mately held and impressive visual results have been achieved using this technique
despite the inaccuracies.
A generalization of Levoy's method is based on a technique known as maximum-
likelihood classification. This heuristic technique requires that probability density
functions Pj(?) are given a prior? for each material ??. The value of Pj(?) is the
27
probability the raw data value ? will be output by the scanner for a voxel that
contains some amount of material m?. A major drawback of this technique is the
fact that Fj (?) are rarely known in advance, and must be estimated. Assuming
that the Pj (?) functions are known, the probability that data value ? will represent
matenM ntj is estimated by
The maximum-likelihood algorithm classifies the data value as being the material
with the highest p?(?) value. Researchers at PIXAR [DCH88] adapted this scheme
to compensate for the partial volume effect by classifying each data value as a
mixture of materials. The fraction of each component material is given by p?(?).
Since
--H--H1,
the entire sample volu'?e is accounted for with this scheme. In practice, the ab-
straction of mapping of probabilities to fractional occupancy works quite well.
The maximum-likelihood classification scheme is easy to evaluate, and has the ad-
ditional benefit of generalizing to multidimensional input. However, it is limited
to domains where values of Pi(?) are available in advance, and, even in the best of
conditions, is only an estimate.
Another classification scheme presented by Levoy [Lev88] is useful for visual-
izing data isosurfaces. This problem is more difficult than simply marking sample
points with values exactly on the surface, since the isovalue might fall in between
sample points, leaving holes in the rendered image. Levoy's method avoids this
28
pitfall by giving the isosurface a thickness, thus classifying voxels that are close to
the isovalue as partially containing the isosurface. Distance from the isosurface is
estimated by using the magnitude of the gradient as an approximate rate of change
in the neighborhood around a sample point. Using this estimate, the distance of
a sample at (x, y, z) from an isosurface with value ??? is given by
?(w,y,z)--H?50
V?(w,y,z)I
Techniques for computing gradients will be covered in Section 3.1.3. Levoy suggests
that a voxel's fraction of isosurface fall off linearly between d 0 and d r, where
r is the desired thickness of the isosurface. This technique has been demonstrated
to produce realistic isosurfaces of electron density maps [Lev88i.
The automatic classification schemes listed above have been employed to pro-
duce many successful volume-rendered images. However, their performance at
object-recognition does not approach human abilities. For example, in medical
imaging, a radiologist can easily identify the boundaries of different types of soft
tissue on CT scans, such as the liver and the kidney. This type of distinction is not
possible using the techniques presented in this chapter. Often, such classifications
must be performed manually, with trained users identifying objects on individ-
ual slices of data. More sophisticated automatic techniques such as model-based
recognition are being explored in fields such as robotics and computer vision. They
employ domain-specific knowledge to aid in classification. In a medical imaging
system, this knowledge might include the expected shape of the kidney. It is likely
that model-based recognition will find future application in scientific visualization,
29
however to date they have not been widely adopted. This is partly due to the fact
that most scientists are interested in seeing their data in essentially `?raw" form
with a minimum of automatic interpretation.
3.1.2 Absorption
Once the volume is classified, the mapping of raw data to absorption may be per-
formed based on the desired transparency of each material. Since absorption is
measured per unit length, it is also necessary to consider the expected thickness of
the material. The relationship between transparency, T, absorption, ?, and mate-
rial thickness, d can be derived starting with Equation (2.7), which we reproduce
here:
I(a, b) --H j????Q(V) dv
This equation can be used to compute transparency by assuming a unit backlight
value, 1b,
T(a,b) = e?LbQ(v)dv			(3.1)
We then introduce the absorption coefficient and material thickness,
T = e?fod?dv
and solve the integral:
T =
Solving for ? yields an expression for the attenuation coefficient, given the desired
transparency and material thickness,
?=? log(T),
d
30
If the classification of a sample point yields more than one material, the absorption
coefficient for the sample point can be set to the weighted average of the absorp-
tion coefficients of the component materials. The weights are determined by the
fractional occupancy of each material.
3.1.3 Emission
The emission function is typically determined by two factors: self-emission and
reflected light. It therefore corresponds to the "source term" found in the radiative
transfer literature. Selftemission is an inherent glow of the material. This could be
used to highlight muscle tissue by making it emit red light. Self-emission is useful
for generating colored fog-like effects, but gives no surface information.
Reflected light is determined by a light reflection model. The simplest light
reflection models are based on phase functions. Phase functions depend only on
the phase angle, defined as the angle between the viewer and the light source as
measured from the reflecting particle. For example, if the particle is spherical,
significantly larger than a wavelength of light, and the surface is diffuse, the refiec-
tivity will depend on the fraction of the illuminated hemisphere of the particle that
is visible from the eyepoint. Thus, the phase function will decrease monotonically
with increasing phase angle. Figure 3.2 shows views of an illuminated particle with
two different phase angles. To relate this to a familiar effect, note that the shapes
in Figure 3.2 are identical to those associated with phases of the moon. Phase
functions have been used successfully to render images of clouds and fog [BliS2,
KV84,RTS7]. Their usefulness is limited for visualization, since shading depends
31
only on the relative positions of the particle, light source and eyepoint. It does not
reveal any information about the particular dataset.
(p
(a)
(p
(b)
Figure 3.2: View of a particle with two different phase angles, y: (a) 90 degrees,
(b) 45 degrees. Above is a top view of the scene, showing the phase angle. Below
is the view of the particle from the eyepoint.
An important advance in the power of volume rendering came when gradient
shading was introduced [HB86,Lev8S,DCH88] Surface reflection models, which
had been used in computer graphics for years, shade objects in such a way that
humans intuitively discern the orientation of their surfaces. These lighting models
always depend on the normal to the surface in question. The problem for volume
rendering is that since there are no explicitly-defined surfaces, there is no inherent
concept of surface normal. The solution has been to use the gradient of the volume
32
data as a substitute for the surface normal. The gradient will always point in the
direction of maximum change in the data. Thus, if a boundary exists in the data
the gradient will point perpendicular to it.
Since the gradient of the data is not known independently, it is computed
using a discrete approximation. The two most common techniques are forward-
differencing [DCH88],
and central differencing [Lev88j,
G(x, y, z) --H V(X?Y+?YAZy)?V(X?Y?Z)
Az
\?+Ax?y,z)--Hv(x--HAx,y?z)
Ax
C(x, ?, z) --H ?21			v(x?Y+AY?z)?A(X?Y?AY?z)
v(x?y,z??z)--Hv(x?y,z--H?z?
Az
On a regular rectangular grid, computation is simplified by computing the gradient
only at grid points, and choosing Ax, Ay and Az to correspond to the sample
spacing along each dimension.
Two problems exist in using the gradient as a surface normal in shading. The
first is that light reflection models assume that the surface normals point outward
from the ?`front" of the surface. When using the gradient as a surface normal, this
direction is automatically determined to be towards the side with higher density.
This is not necessarily the side that the user considers the front. The problem can
be avoided by using a light reflection model that reflects equally on both sides of
33
the surface, or by manually inverting the gradient directions in regions where an
error is made.
A more thorny problem is that the gradient operator acts as a high-pass filter,
which serves to emphasize noise in the image. The visual result is that the normals
jiggle and do not give the impression of a smooth surface. Note that the noisy
gradient is not always a disadvantage; it adds texture to surfaces which, along
with a perspective projection can be used as a depth cue. Two approaches to
minimizing noise in the gradient have been proposed. First, the gradient can
be taken on a smoothed version of the dataset [DC1188], while the unsmoothed
version is used for the other steps in rendering. Another approach is to compute
the gradient over a larger domain than used in forward-differencing and central-
differencing [HBP+90].
There is some debate in the literature about whether the gradient should be
taken on the classified data or on the raw data values. Using classified data ensures
that the shading is consistent with the objects that are displayed [DCH88]. Levoy
argues that by using the raw data values, errors in classification can be minimized
since proper shading based on the "real" values is always present [Lev88i. This
issue remains unsettled. For the images presented in this thesis, we use the gradient
of the raw data, however, use of the classified gradient would not change the
algorithms or analysis contained herein.
Once the gradient is computed, any of a number of local shading models may
be used. The most popular is the Phong shading model [BT75], which has the
34
general form,
b ka + kdcos(0) + k?cos?(?),
where 0 is the angle between the vector pointing towards the eyepoint, and the
surface normal, and ? is the angle between the vectors pointing towards the view-
point and towards the direction of mirror reflection. These vectors and angles are
illustrated in Figure 3.3. The Phong model is the sum of three terms, an ambient
term, a diffuse term, and a specular term. There are also four constants: ka, kd,
k8 and n. The ambient term is simply the constant ka, and acts as an approxima-
tion to the effects of multiple scatter. The diffuse term models ideal Lambertian
reflection and is weighted by the constant kd. The effect of a specular highlight
is approximated by the third term, which is weighted by the constant k8. The
sharpness of the highlight is controlled by the constant n.
N
R
Figure 3.3: The angles used in Phong illumination, 0 and ?, are shown in the above
diagram. The surface normal is N and the direction of the maximum specular
highlight is 1?. After Foley et. al. [FvFH9O]
The Phong lighting model is not physically-based, however it does produce
35
realistic shading effects. The diffuse term and particularly the specular term are
highly sensitive to deviations in surface normal and thus communicate the orien-
tation of the surface quite well. This model also lends itself to easy evaluation,
since the cosine terms may be evaluated using dot products.
Research has led to more realistic local illumination models than the Phong
model [CT82,HTSG91l. Although there are no technical barriers to their applica-
tion, the use of these models has not been reported in the volume rendering lit-
erature. Increased computational cost and lack of interest in photorealism among
the visualization community are the likely causes for researchers remaining with
the less sophisticated model.
As a demonstration of the effectiveness of applying local illumination models
based on the gradient, we present in Figure 3.4 an image computed using gradient
shading and the same image computed using a phase function.
3.2 Reconstruction
Although it is possible to produce volume data that is defined by mathematical
functions in three dimensions, the overwhelming majority of volume data is pre-
sented as sample points in 3D space. If these datasets were projected without
any type of reconstruction, the resulting image would appear as a scatterplot. In-
stead, we wish to reconstruct the input to the renderer so that the final image
gives the impression of objects that occupy volumes of space. Construction and
reconstruction of functions defined on a continuous domain from point samples is
a well-studied problem in the field of signal processing. What follows is a brief
36
(a)			(b)
Figure 3.4: Comparison of image (a) with gradient shading and (b) with shading
by a phase function.
overview of the techniques used in volume rendering and their benefits and costs.
More detailed information about the reconstruction problem, can be found in an
introductory signal processing text, such as the one by Oppenheim and Schafer
[0575].
Signal processing theory states that a function may be perfectly reconstructed
from its samples if it is band-limited and if it is sampled at twice the maximum
frequency of the data. There are two practical difficulties in applying this ideal
reconstruction to sampled volume data. The first is that the ideal reconstruction
filter is a sinc function, which has infinite extent. Thus, to reconstruct any point
in 3D space, we must consider the values of all sample points, not just those in a
local neighborhood. This quickly becomes computationally intractable. Another
difficulty is that real-life datasets tend not to represent band-limited scenes. Any
37
scene that contains a sharp boundary between two media will represent a step
function in 3D which has infinite extent in the frequency domain. Use of ideal
reconstruction on a dataset which is not band-limited will cause ringing artifacts.
The effects of the different filters that we will present in this chapter will be
demonstrated using a one-dimensional reconstruction problem. The problem is the
reconstruction of a scanline of a 2D image that has been intentionally undersampled
to reflect the kind of difficulties that are faced when visualizing real-world scientific
data. Figure 3.5 shows the original image, and a plot of the selected scanline
showing the samples that were taken. Note that the original image was digitally
scanned, so it is also made up of a discrete number of samples. The plot of
the scanline in Figure 3.5b was created by joining the sample points with linear
segments. Figure 3.6 shows the result of reconstruction with a sinc filter applied
to these samples. Note that the band-limited output of this filter causes ringing
in areas where the actual function is flat.
25
200
150
loo
2			4			b			?			10 12 14
(a)
Figure 3.5: Sample reconstruction problem: (a) 2-D digitized photograph indicat-
ing scanline to be reconstructed. (b) Raw scanline data indicating subsampling
points.
38
10
(a)			(b)
Figure 3.6: R?econstruction with a sinc filter: (a) Filter shape in the x interval
L2,21, (b) reconstructed scanline.
Nearest-neighbor interpolation is the most basic form of reconstruction used
in volume rendering [SSW86, Sab88j. This method states that the value of the
function at any point in three-dimensional space is equal to the value of the nearest
sample. Nearest-neighbor interpolation has the advantage of being very easy to
evaluate. However, the reconstructed signal is distorted, with sharp discontinuities
between cells. This gives images computed using nearest-neighbor interpolation a
distracting, blocky appearance. Figure 3.7 shows the effect on the sampled one-
dimensional data. Another feature of nearest-neighbor interpolation is that the
only output values it produces are those that exist in the raw data --H no values
are synthesized. This is an attractive feature to those more concerned with the
accuracy of the output values than with the spatial positioning of those values.
The degree of continuity between cells can be increased by employing trilinear
interpolation. This thre&dimensional analog of linear interpolation operates on
hexahedral cells with sample values at the corners. A trilinear function has the
39
a ? ? 10 12 la
(a)			(b)
Figure 3.7: Sample scene reconstructed with nearest-neighbor interpolation: (a)
Filter shape, (b) reconstructed scanline.
form:
f(r, y, z) C0 + C1r + Qy + c3z + c4ry + c5yz + c6rz + c7wyz.
The coefficients ..... c7 are computed by solving a linear system. The requisite
eight equations may be written by substituting the r,y,z and f values at the eight
corners into the above equation. If the grid-spacing is regular, the matrix may be
inverted as a preprocess [UKS8j. Trilinear interpolation represents a good trade-
off between smoothness of the output function and computational cost. It is by
far the most popular form of interpolation used in current volume rendering al-
gorithms [Lev88,UK88,WG91]. However, it has some drawbacks. It is limited by
its requirement that there exist hexahedral cells, therefore it is not appropriate for
irregularly-sampled data. Furthermore, the derivative discontinuity at cell bound-
aries can result in perceptual anomalies such as Mach-banding [Rat72]. The result
of linear interpolation is shown in Figure 3.8.
Piecewise constant and trilinear interpolation are both forms of a general in-
terpolation strategy called polynomial spline interpolation. These result in the
40
(a)			(b)
Figure 3.8: Reconstruction with linear interpolation: (a) Filter shape, (b) recon-
structed scanline.
definition of a separate polynomial function for each cell. The order of continuity
between cells can be increased by increasing the degree of the polynomial function
within each cell. To obtain a C1 output function, tricubic interpolation may be
used. In this case, the function can be defined across a hexahedral cell by con-
sidering all the sample points in a 3x3x3 cell neighborhood, for a total of 64 data
values. The use of this computationally intensive technique has not been reported
in the volume rendering literature.
Another approach is to use a Gaussian filter kernel [Wes89,Wes9O,LH911. In
its ideal form, a Gaussian filter kernel has infinite extent In practice, however, it
is clipped so that it extends over only a few sample points. Gaussian interpolation
has the drawback that its output value at the sample points may differ from the
input values at these points. Furthermore, the extreme smoothing that takes place
tends to give the dataset a blurry appearance. Figure 3.9 demonstrates Gaussian
interpolation on the one dimensional example.
41
o+;W??'
10 12
(a)
Figure 3.9: Sample scene reconstructed with Gaussian interpolation: (a) Filter
shape, (b) reconstructed scanline.
3.3 Projection Algorithms
During the projection phase of volume rendering, the three-dimensional scene is
mapped onto a two-dimensional image plane. In volume rendering, projections are
defined by a view transformation and an illumination model. The view transfor-
mation determines the set of points in space that are mapped to a single point on
the image. The illumination model (several of which were discussed in Chapter 2)
determines the color that results from this many-to-one mapping.
3.3.1 View Transformations
The two most popular view transformations are orthographic projection and per-
spective projection. Both require the definition of an image plane. Using an
orthographic projection, the set of points that share an image-space mapping are
those that lie along a common line perpendicular to the image plane. Using per-
spective projection, a center of projection must be defined. The set of points that
share an image-space mapping are those which lie along a common line that passes
42
through the center of projection. Further information about projective geometry
may be found in introductory computer graphics texts, such as the one by Foley
et. al. [FvFH9O].
Perspective projections are desirable since they provide familiar perceptual
cues. However, orthographic projections are usually simpler to compute because
the sets of points that form an equivalence class after projection lie along parallel
lines. This parallel nature of the projection enables the programmer to uniformly
sample world space by uniformly sampling screen space. Further discussion of
projection will be presented in Chapter 5.
Regardless of the choice of view transformation, a two-dimensional reconstruc-
tion step must be added to the rendering process. The illumination function is
defined at all points on the image plane. However, the image must be stored in
limited memory and must be displayed on devices with only a finite number of
rectangular picture elements called pixels. It is assumed that each of these pixels
displays a single constant value across the entire pixel area, although in reality
there is some variation. Typically, the images will be stored as point samples that
correspond to the function's values at pixel centers. Alternatively, the point sam-
ples may correspond to an average of function values across pixel area. Further
information about generating pixel values may be found in Foley et. al. [FvFH9O].
3.3.2 Partitioning the Integral
Before discussing specific approaches to computing projections, we present a tech-
nique for partitioning the volume rendering integral so that solutions to the integral
43
for individual ray segments may be combined to form a solution for an entire ray.
Partitioning in this way is employed in nearly all existing volume projection algo-
rithms.
To make the partition, we first rewrite the integral using the transparency
function, T(a, b) that was first defined in Equation (3.1),
T(a, b) = c?LbQ(v) dv			(3.2)
Recall that this function has a physical interpretation; it is the transparency or
optical depth of the material between points a and b along the ray. The volume
rendering integral can now be rewritten.
I(a,b) = f?T(a,u)?(u)du.			(3.3)
Observe first that transparency obeys a simple multiplicative rule. In particu-
lar, for any point r in the ray segment [a, b], we have
T(a,b) = T(a,x)T(w,b),
as is easily verified from equation (3.2).
We can therefore discretize the ray into n contiguous segments [d0, d1], [d1, d2],
[dn?i,dn] at distances a = d0 < d1 <; . ? = b from the eye point, and
rewrite equation (3.3) as follows:
I(d0,dn)
T(d0 ,u)e(u)du
ffi f??1T(d0 di?i)T(di?i, u) ?(u) du
44
liT(d0,di?1)f??1 T(dj?1, u) ?(u) du
T(d0, dj?1) I(dj?1, di).
Finally, by introducing the notation 1k I(dk-1, dk) and Tk T(dk?l, dk) we
can rewrite this equation more succinctly:
J(d0,dn)			--H
LT(do,di?1)
i--Hi
--H z `? ii Tj,			(3.4)
z=1 j=1
where, H, the product operator is defined,
fi Q a1 x a2 x ... x a?.
i--H--Hi
Equation (3.4) expresses the overall intensity of the ray in terms of the intensities
and transparencies of the discrete segments.
Equation (3.4) can also be expressed in terms of image compositing operators.
These general-purpose operators are used to overlay semitransparent images. In
order to use them, it is necessary to store tuples that associate intensity values and
transparency values with one another. We shall make use of two of the operators,
over and under defined by Porter and Duff [PD84]. Given pixels A and B, with
color and transparency values [Ac, At] and [Bc, Btj the operators are defined,
A over B [Ac t (Bc X At), At x Bt]
A under B [(Ac x Bt) t Bc, At x BtJ,
and correspond to image overlays and underlays. Our notation is slightly different
from Porter and Duff's since we use transparency in our tuples rather than opacity.
45
If we define Ci to be a tuple containing both the intensity and opacity value for
the segment i,
Cj fli, Ti],
then Equation (3.4) may now be written:
[I(d0,dn), T(d0,dn)] = Ci over (C2 over
This formulation is significant, since fast hardware implementations of image com-
positing operators are available on most graphics workstations.
3.3.3 Ray Casting Algorithms
The first category of projection algorithms that we will consider are the ray casting
algorithms [Lev88, Sab88, UN88]. These algorithms work by point-sampling the
image plane. At each of these point samples, the illumination is determined by
firing a ray through all object space points that are mapped to the image-space
sample point. The ray is broken into small regions for which approximations of
I? and Tj are computed. The results are composited using the techniques of the
previous section.
Ray casting algorithms generate approximations of I? and Tj using quadrature
rules that will be described in Chapter 4. Evaluation of these rules require that
either
and e are reconstructed at sample points along each ray, or
o+ an analytic expression for the reconstructed signal in generated for each re-
gion.
46
In either case, the cost of reconstruction dominates the computation.
Research in ray casting has stressed efficiency improvements. In standard ray
tracing, Arvo and Kirk [Gla89] found three basic categories of improvements. Two
of these categories have direct analogies in volume ray casting, as noted by Levoy
[Lev89]. Algorithms in the first category seek to reduce the number of samples
of the image plane that are taken. For example, adaptive supersampling uses a
heuristic to determine regions of the image plane where the illumination function
is slowly varying [Whi8O,Lev9Oa]. The second category of efficiency improvements
seek to reduce the time spent computing the illumination value along each ray.
Octrees can help by enabling the algorithm to skip over empty regions or make
coarse approximations in regions that experience little variation of emission and
absorption [Lev9Oa, DH92j. Early ray-termination algorithms can be used when
rays are traced from front to back. In this case, once a user-specified amount of
dense material has been encountered, the remainder of the ray is ignored since
any structures encountered would largely be occluded. In standard ray tracing,
most speedup techniques introduce no new inaccuracies, however this is not the
case in volume ray casting. The main drawback to all the above techniques is that
they introduce approximation error into the final image without characterizing its
magnitude.
There are several drawbacks to the ray casting approach. Since rays are nor-
mally traced on&at-a-time, references to memory are not well-localized, causing
paging problems with large datasets. Furthermore, the algorithm operates in screen
space, and typically does not adjust based on the resolution of the dataset. Thus,
47
the scene tends to be oversampled in some instances, and undersampled in others.
We address both these problems in Chapter 5. Additionally, ray casting does not
lend itself to distributed-memory parallel implementations, since each processor
requires access to the entire dataset.
Despite these disadvantages, ray tracing remains one of the most popular vol-
ume rendering algorithms. It is conceptually quite simple, and easy to implement.
Ray casting software can be tuned to change projection type or reconstruction fil-
ter with a minimum of effort, By varying the number of samples taken per pixel or
the quadrature rule that is employed along each ray, the user has the opportunity
to trade rendering time for image quality.
3.3.4 Data-Space Algorithms
Whereas ray-casting algorithms have a main loop that iterates once per image-
space sample, data-space algorithms iterate once per data sample. Ray casting
algorithms ask how each data sample affects the current image plane sample. Data-
space algorithms ask how each image-space sample point is affected by the current
data-space sample.
The Pixar Algorithm
The approach introduced by Duff el. al [DCH88] is optimized to employ image
compositing hardware. The basic algorithm requires that
o+ Data is stored in a regular, rectangular grid;
48
o+ An orthographic projection along one of the principal axes of the dataset is
desired; and
o+ Screen resolution is equal to the 2D resolution of the dataset along the axis
of projection.
However, the algorithm is more general than the above restrictions imply. Many of
the preconditions may be satisfied by employing additional resampling operations.
Irregular data may be resampled onto a regular grid. Views that are not aligned
with data axes may be computed by rotating the data within the grid so that
the view may be achieved along one of the principal axes. These rotations may
be performed efficiently using the three-pass technique described by Hanrahan
[Han90]. During each pass, a 1D shearing operation is performed along one of
the axes, as shown in 2D in Figure 3.10. Perspective views may be computed
by predistorting the dataset, so that orthographic projections of the predistorted
dataset exhibit familiar perspective depth effects,
Once the data has been manipulated so that all the preconditions apply, the
data is processed and stored into a single thre&dimensional array of (r, g, b, ?)
values, where r, g, and b are determined by computing the value of L for a ray
passing through the center of each voxel along the view direction. The ? values are
determined by computing the corresponding Tj value. Then each two-dimensional
slab of the data can be treated as an image with an associated alpha channel and
composited using the over and under operators defined in Section 3.3.2.
50
poor choice for both very small and very large datasets.
Splatting
Splatting is a data-space technique that seeks to exploit the coherence that is
present in certain types of reconstructed data. Assuming the samples are regularly-
spaced, the 3D functions of emission and absorption are simply sums of scaled,
translated copies of the reconstruction kernel. The main concept in splatting,
which was introduced by Westover [Wes89,Wes9O], is that each of these scaled and
translated copies of the reconstruction kernel can be projected individually, while
summation (or compositing) occurs in screen space. The coherence that splatting
exploits arises when computing an orthographic view. In this case, the size and
shape of the projected filter kernel does not vary with 3D position, so it can be
projected once as a preprocess.
We now consider the steps involved in the splatting algorithm. The recon-
struction filter, usually a Gaussian [Wes90], is integrated in depth, thus creating a
generic screen-space footprint that is stored in a lookup-table at subpixel resolu-
tion. This process is illustrated in Figure 3.11. The data samples are considered
from front-to-back or from back-to-front. The footprint is centered at the data
point's screen-space projection, and its values are scaled based on the emission
and absorption values at the data point. Finally it is composited into a final image
using over or under operations.
This algorithm differs from other approaches to volume rendering in that the
processes of reconstruction and projection are reordered. Unfortunately, the emis-
51
3D Spherical
Filter Kernel
Projected 2D
Footprint
(a)
Figure 3.11: A preprocess to the splatting algorithm is the creation of a generic
footprint table. (a) A 3D reconstruction filter over a spherical domain is projected
down to a 2D footprint. (b) The generic footprint table stores the integral of
the reconstruction filter at a discrete number of sample points. After Westover
[Wes9Oi.
sion-absorption and low-albedo models of volume rendering are non-linear in c
and 6, hence reconstruction and projection are not associative in the general case.
They are associative when filter kernels do not overlap in space, as noted by Laur
and Hanrahan [LH91], but such filters cannot ensure C1 or higher continuity.
Splatting can be seen as an approximation to the emission-absorption model.
Filters that overlap in space are approximated by those that are spaced so that they
do not. This is illustrated in Figure 3.12. In Figure 3.12a, a function along a ray
is shown, along with the scaled and translated filter kernels which are its addends.
This function can be considered emission or absorption. The filter kernels are
clipped Gaussians. Note that these curves overlap with one another along the
52
ray axis. Therefore, although the ray encounters most of curve A before curve
B, portions of curve A extend into the domain of curve B and vice-versa. The
splatting algorithm renders each of these curves as if they do not overlap along
the ray axis, as shown in Figure 3.12b. This approximation can result in values of
illumination that are either too bright or too dark, depending on the data. Studies
of the magnitude of this error have not been published.
Images produced by splatting tend not to have the same level of detail that is
present in ray-traced images. It is not clear how much of this loss of detail is due to
the use of a wide Gaussian filter kernel, and how much is due to the approximation
error described above.
The initial implementations of splatting required regular, rectangular datasets
and orthographic projections. At the cost of adding additional approximation
error, Laur and Hanrahan introduced modifications to handle rectangular hierar-
chical data and perspective projections [LH91]. A further algorithmic modifica-
tion they introduced was the approximation of the splats by Gouraud-interpolated
(r, g, b, o?) polygons. This enables the use of hardware polygon engines for rendering
the splats.
Splatting is attractive because it is easy to implement, runs quickly, is easily
parallelized, and can be optimized to use rendering hardware. However, it is
also an approximate technique for which the magnitude of the error is unknown.
Furthermore, this error cannot be minimized by modifying the input parameters.
Splatting is best used to create a "quick-preview" visualization engine, useful in
situations where where interactive feedback is required.
53
6666 B
6			t
(a)
I			"4
I			44
c			-.
(b)
Figure 3.12: Errors due to the reordering of reconstruction and projection in splat-
ting. (a) A function defined by truncated Gaussian interpolation (dashed line) is
the sum of several truncated Gaussians (solid lines A, B and C). (b) The splat-
ting algorithm does not allow filters to overlap in depth. To compensate, the
truncated Gaussians must be translated relative to one another so that they no
longer intersect. Recovery of the function is no longer possible.
54
Cell-by-cell processing
Instead of considering data samples individually, the cell-by-cell processing algo-
rithms consider spatial regions whose value is interpolated from the same set of
data samples. In the case of trilinear interpolation, the regions are rectilinear
cells, with data samples at each of the eight corners. The cells are considered in
front-to-back or back-to-front order, while the image is accumulated in a buffer.
In contrast to splatting, there is no overlap in depth between the cells. Thus, the
problem of splatting error described in the previous section is avoided.
Figure 3.13: Decomposition of a cell into subcells. The projected rectilinear cell
at left is into divided into the seven subcells sharing the same front and back face,
an shown in the exploded view on the right. After Shirley and Tuchman [ST9O]
The cell-by-cell processing techniques operate by dividing the cell into regions
called subeciTh that share the same front and back face as shown in Figure 3.13.
The illumination function is assumed to vary slowly within each subcell. Upson
and Neeler's method [UN88] is a scanline algorithm that divides each scanline into
spans along subcell boundaries. The techniques of Shirley and Tuchman [ST9OJ,
and Wilhelms and Van Gelder [WG91] consider the entire screen-space projection
55
of a subcell at a time. They make the approximation that the illumination function
varies bilinearly within each projected subcell. Illumination values are computed
at each region's extremal points to determine the bilinear interpolation coefficients.
Shirley and Tuchman's technique divides projections of tetrahedral cells into tri-
angles, while the approach employed by Wilhelms and Van Gelder divides projec-
tions of rectilinear cells into quadrilaterals and triangles. Projections using either
method can be computed using graphics hardware that renders Gouraud-shaded
(r, g, b, ?) polygons. Wilhelms and Van Gelder make the additional observation
that when computing an orthographic view, the number, size and shape of pro-
jected subcells does not vary. Thus the decomposition into subcells need only be
computed once. They call this technique coherent projection.
Cell-by-cell techniques derive their efficiency from the fact that the illumina-
tion function varies slowly within each subcell. They operate best when each
cell projects onto many pixels. Therefore, they should be considered when screen
resolution is higher than data resolution. The main drawback is that, with the
exception of some experimental results [WG91], analysis of the error introduced
by a piecewise bilinear approximation to the illumination function has not been
performed. When using graphics hardware to render the subcells, error is also
introduced due to the use of fixed-precision arithmetic.
3.3.5 ?requency Domain Algorithms
All the above algorithms for volume projection operate in the spatial domain.
Recently, frequency-domain techniques have been introduced [Mal9O,Lev92,Mal93,
56
4
TL93]. These algorithms are based upon the Fourier slice theorem. This theorem
relates projections of a datset to its Fourier transform, and is illustrated using a
2D example in Figure 3.14. We first define P0(t) to be a parallel projection of the
2D dataset, as shown in Figure 3.14a. The Fourier slice theorem states that the
Fourier transform of P0(t) is equal to a slice through the Fourier transform of the
data, inclined at an angle 0 and passing through the origin, as shown in Figure
3.14b. Further details of the Fourier slice theorem may be found in textbooks such
as the one by Nak and Slaney [NS88].
x
v
Se(w)
(a)
Figure 3.14: The Fourier slice theorem. (a) A parallel projection at angle 0 in the
spatial domain, P0(t) , is the inverse Fourier transform of So(w), a slice through
the origin of the Frequency domain at angle 0 (b)
Fourier volume rendering employs the three-dimensional Fourier slice theorem.
The steps involved in computing a 2D projection, Po,?(s, 1) of a sampled function
57
f(x, y, z) using the Fourier volume rendering technique are:
1.			Compute the discrete Fourier transform F(u, v,			? [f(w, y,
2. Let So,?(q, ?) be a regular resampling of F(u, v, w) at screen resolution, along
the plane passing through the origin at angles 0, ?.
3. Let Po?(s,t) be ?--H` [So,?(q,r)].
Note that the use of this technique does not obviate the need for a reconstruction
phase. Although the 3D function is never reconstructed in the spatial domain, it
is reconstructed in the frequency domain when computing So??(s, t) in step 2.
The model of projection that is assumed by the Fourier slice theorem does not
correspond to any of the volume illumination models presented in Chapter 2. It is
simply the integral of the density function along the ray, i.e.,
I(a, b)			p(t)dt.
This expression may easily be transformed into the attenuation-only model pre-
sented in Equation 2.7 by taking the negative exponential of the result. As noted
in Chapter 2, this model is extremely limited in application since it does not ad-
mit shading or occlusion. Thus, the projection model is a major drawback of
the technique. Levoy has investigated techniques for improving the appearance of
Fourier-rendered images by introducing depth-cueing and linear approximations
to diffuse shading. Totsuka and Levoy implement these effects more efficiently by
filtering in the frequency domain [TL93j.
If we restrict our focus to computing the attenuation-only model, Fourier vol-
ume rendering is attractive because of the improvement in worst-case asymptotic
58
complexity over both image-space and data-space algorithms. Assume a volume
dataset with n3 samples, and an image with m2 pixels. The worst-case asymp-
totic complexity for data-space algorithms is 0(n3) and O(m2n) for image-space
algorithms. Fourier rendering techniques require that the dataset first be trans-
formed using a 3D DFT; an O(n3 log n) preprocess that is performed only once per
dataset. Projections are then be computed by computing 0(m2) interpolations in
the frequency domain, and computing an inverse 2D DFT in O(m2 log m) time.
If we consider the initial Fourier transiorm as an independent preprocess, then
Fourier volume rendering can be considered an O(m2 log m) process.
There are several drawbacks to using the technique as well. The speedups that
are theoretically possible with Fourier volume rendering are difficult to achieve in
practice, due to the cost of reconstruction in the frequency domain [Lev92,TL93i.
The technique is also quite stringent in its preconditions: The density function
must be sampled on a regular, rectangular grid and an orthographic projection
model must be used. The technique is relatively new, and further development of
the algorithm may overcome these disadvantages.
3.4 Summary
The work presented in this thesis builds upon the techniques presented in this
chapter. The methods we explore assume that the inputs to the renderer are
point-sampled functions of emission and absorption. Thus we do not assume any
particular strategy for mapping. For the sample images, we used either maximum-
likelihood or isosurface classification, along with gradient shading based on the
59
Phong model. For reconstruction, we use trilinear interpolation, although our
methods will be general enough to handle any polynomial spline interpolation.
The projection algorithms that we develop are extensions to ray casting. We felt
that this was the best choice, because of its flexibility, and the ability to control
error.
Chapter 4
Controlled-Precision Volume
Rendering
The process of volume rendering scientific data involves many approximations. In
some cases, the approximations are introduced intentionally, such as when the pro-
grammer seeks to increase processing speed. Other approximations are inherent
in the model of illumination; for example there is no known analytic solution for
the volume rendering integral of Equation (2.6). Approximations in the rendering
process can lead to image artifacts that are not faithful to the data being ren-
dered. At present, there are no tools to evaluate these artifacts quantitatively,
so their presence and severity are evaluated subjectively. Often, the possibility of
inaccuracy in rendering is not considered at all.
In this chapter, we introduce the notion of controlled-precision volume render-
ing. Algorithms with controlled precision differ from previous approaches in that
60
61
they include both quantitative assessment and user control over the quality of the
rendered image. Measures of quality are in terms of maximum deviation from
the mathematically exact solution to the rendering model. Controlled-precision
volume rendering has several advantages over existing algorithms:
o+ Using controlled-precision rendering, images can be computed that are guar-
anteed to have no visible rendering artifacts. Artifacts are image features
that would not be present in the mathematically exact solution.
o+ Controlled-precision methods can efficiently compute images by halting com-
putation when the image meets the desired precision.
o+ Controlled-precision methods can be used to compute reference solutions that
can serve to evaluate other volume rendering algorithms quantitatively.
Guarantees of accuracy are essential in fields where scientists draw critical con-
clusions about their data based on the images that are produced. For example,
in medical imaging, physicians must be able to determine if a small, anomalous
dark spot on an image is indicative of disease or is merely a rendering artifact. In
fact, physicians should not have to worry at all about the presence of rendering
artifacts. Since controlled-precision techniques can provide reliable solutions to
any specified accuracy, they can solve these problems by ensuring that images are
computed to the accuracy of the display device.
In the first section of this chapter we present controlled-precision techniques for
evaluating the model of volume illumination along a ray passing through a single
voxel. In the second section, we extend these results to include intersections of rays
62
with an entire volume. In the final section, we address the problem of accurate
computation of pixel values.
4.1 Controlled-Precision Volume Integration
We now turn our attention to controlled-precision methods for solving the volume
rendering integral, introduced as Equation (2.6) in Chapter 2.
I(a, b)			f6e?iJQ(v) dv 6(t) dt.			(4.1)
For convenience, we define f(t) to be the integrand of Equation (4.1),
f(t) --H e?LtQ(v)dv 6(t).			(4.2)
Note that each ray has a different function of ? and 6 and hence a different f(t).
No general solution to fab f(t) dt using a finite number of arithmetic operations
and trigonometric functions is known. Thus, we must rely on some form of numer-
ical quadrature. Although many standard techniques exist, previous approaches to
solving this integral for volume rendering have not included error bounds. Many
volume rendering algorithms evaluate the integral using a piecewise constant ap-
proximation the domain of integration is broken into small segments where 6
and a are approximated by constant values, k( and ka [Lev88,Sab88,WG91]. Under
this restriction, the integral of each of the segments reduces to a simple arithmetic
expression,
$bf(t) dt
--H $b??LtkQdvk?dt
I(a,b)
63
--H k?			dt
--H kka? F?k?[b?a] --H
--H e?k?[b?ai1
(4.3)
kQ
Note that the piecewise constant functions for which the technique is named are
approximations to 0 and e. The resulting approximation to the integrand is not
piecewise constant.
If o(v) is symbolically integrable, a similar arithmetic expression may be found
by dividing the domain of integration into small regions where ? is assumed to be
related to c by a constant kQ? [MHC90,WG91]:
J(a,b) =
f(t)dt
fbe?LtQ(v)dv ka?o(t) dt
= koE f6e--HJJQ(v)dvo(t) dt
--H			kQ?			--H e?a(v)dv1
= kQ?			--H e?LbQ(v)dv1
The approach employed by Williams and Max [WM92] approximates 0 and e in
small regions by piecewise linear segments. The solution of the integral for each
segment includes an error function that can be approximated to any specified
precision. Another common technique that is applicable when 0 is symbolically
integrable is to solve the inner integral of Equation (4.1) directly, and to approx-
imate the outer integral using a deterministic quadrature formula [NV84,Uk88,
64
Kru9Oi, or via Monte Carlo methods [I?N89j. While all these strategies are known
to increase in accuracy as the number of segments or samples increases, no precise
error bounds have been given in the literature.
We present four techniques for computing approximate solutions to the volume
rendering integral with bounded error. The first method, which we call piecewise
constant bounds propagates bounds on d and e using interval arithmetic. The
remaining three methods are the trapezoid rule, Simpson's rule, and a power series
method. In all cases, by bounding the magnitudes of the derivatives of f, we
are able to obtain remainder terms that in turn bound the approximation error.
The convergence properties of these algorithms vary considerably. The method
of piecewise constant bounds is practical only for low-precision solutions. The
trapezoid rule and Simpson's rule are most efficient for low- and medium-precision
solutions. The power series method efficiently computes very high precision results
making it the method of choice for computation of highly accurate solutions. Such
results may be used as reference standards by which the output of other algorithms
may be compared. Note that much of the work in this section is based on standard
numerical techniques. The details of standard methods are included for reasons of
completeness.
4.1.1 The Integrand and its Derivatives
We begin by exploring the integrand and its derivatives. In this section, we focus
our attention on individual voxels and defer discussion of collections of voxels
until Section 4.2. We assume reconstruction by polynomial spline interpolation.
65
The functions ?(v) and 6(t) are therefore polynomial. Since a is a polynomial, the
exponent of Equation (4.2) evaluates to a polynomial, and we rewrite the integrand
in a more convenient form
f(t) --H eP(t)q(t)			(4.4)
where p and q are polynomials. In order to give a sense of the typical shapes of
these functions, graphs of p(t), q(t) and f(t) for two functions a(v) and 6(t) are
given in Figure 4.1.
A prerequisite for computing solutions to the volume rendering integral with
error bounds will be bounding the derivatives of f(t) for all tin an interval [a, b].
In addition, the power series method will require exact values of the derivatives of
J(t) at t --H 0. For this reason, we begin our discussion by presenting a formula for
the nth derivative of f, denoted f(n), in terms of the derivatives of p and q.
Theorem 4.1 (Derivatives of the Integrand) Let f(t) --H eP(t)q(t), where p
andq are n t?mes d?fferentiable. Then
f(n) (t) --H n! ? q(i)(t)An?j(t)
j--H--HO
where An?j(t) ? defined recu??ve1y by
(4.5)
A0(t)			--H ep(t)			(4.6)
An(t)			(4.7)
1 p(k)(t)An?k(t)
k--H1
Proof: The theorem follows from the fact that
An(t)			ep(t)
6(3
cp(t) e?f?Q(v) dv			q(t)
0.8			0.8 N
0.6
0.4			0.4
0.2			0.2
0
0.2			0.4 0.6 0.8			1
I			I
f(t) eP(t)q(I)
0
0.8
0.6
N
0.2
0
0.2			0.4			0.6			0.8
0.2			0.4			0.6			0.8
Figure 4.1: Graphs of the volume rendering integrand. (top) A ray is passed
through a rectilinear voxel with samples of ? and 6 at the corners. (middle) This
gives rise to functions of attenuation, ep(t), and emission, q(I) along the ray. (bot-
tom) These two functions are multiplied to obtain the volume rendering integrand,
f(t) (bottom). The solid and dashed lines present these functions for two different
input voxels.
67
and from the formula for the ?th derivative of a product,
fg  f(i)g(fl?i) fl
i--HO			i
Equation (4.5) may then be verified by induction. ?
Although equations (4.5)-(4.7) are formulated to give values of the derivatives
of f at a point, they can be modified to give bounds on the derivatives of f over
an interval. The only change necessary is to replace instances of derivatives of p
and q with bounds on those derivatives over the interval. Derivatives of p and
q are polynomials, so bounds can be obtained using the following rule. Given a
polynomial ?(t) of degree D with
t ?t t... t pDtD
an upper bound on ? for t ? [a, b] is given by
1% t P1c +...+ PDcD
(4.8)
(4.9)
where c --H max(?a , Ib). This is a coarse bound but has the advantage of requiring
no root finding.
4.1.2 Piecewise Constant Bounds
In Equation (4.3), we showed that the volume rendering integral reduces to a simple
arithmetic expression in ray segments where & and ? are constants. The method of
piecewise constant bounds relies on the fact that this formula can also be used to
form upper and lower bounds on I(a, b), given piecewise constant bounds on
and e(u).
68
We begin by dividing the domain of integration, [a, bl,into n subregions, [d0, d1j,
[d1, d2], ..., [dn--H1, dn] as in Section 3.3.2. Let ?k be an interval [?k7, ?t1 sud? that
o < ? ?(t) ? ?k+ for all tin [dk?l,dk]. We solve for an interval Tk [Thk,i71
bracketing the exact transparency Tk in the kth subregion using standard interval
arithmetic [Moo66l.
Tk			Ldk?1 Qkdv
--H e?Qk (dk--Hdk?1)
(4.10)
The rules of interval arithmetic are quite simple to define in this case, since the
two operations, integration and negative exponentiation, are monotonic. To get the
upper bound on the interval for Tk, we use the lower bound on ?k and vice-versa.
These rules correspond to the intuitive notion that transparency and absorption
are inversely related.
Similarly, if we bracket the emission function ? within the segment [dk--H1, dk]
by an interval Ek --H [?7?, ?k+?l' then it is also possible to solve for 1k [1k?? 1k+]' an
interval containing the exact intensity in the kth subregion, 1k, as follows:
1k = fd?kl e?Ltk?l Qk dv Ekdt
--H Ek(l--HTk)/oLk.			(4.11)
Once again, the interval arithmetic corresponds to an intuitive notion. Maximum
intensity is obtained using maximum emission and minimum absorption. Min-
imum intensity is obtained using minimum emission and maximum absorption.
Mathematically, the computation of 1k is a little more complicated than for Tk
69
due to the possibility of a zero in the denominator. Equation (4.11) can be broken
into three cases:
Ek(1 --H Tk)/Qk			if ?k # 0 and ?k+ # ?
1k= [6k(1--HI+k)/0?+k %(dk--Hdk1)1 ff??=0ando?+#O
Ek(dk --H dk-I)			if ?k = 0 and ??+ = 0.
We can then bound the overall intensity for the ray with an interval
I(a,b) = k--H(a,b),I+(a,b)1
(4.12)
by rewriting Equation (3.4) with interval arithmetic:
n			i--Hi
I(a,b) = Z `? II 7)
?=1			j--H--H1
Finally, by introducing the notation for the width and midpoint of an interval,
we may state:
II(a,b)			I??(a,b) --H I?(a,b)
I (1+ (a b) t I?(a, b))
I(a, b)			=			2
(4.13)
I(a,b)--HI(a,b) ?2I I(a, b)
This means that the midpoint of the interval may be used an as approximation to
I(a, b), with a margin of error of I?21 I(a, b)
Up until now, we have assumed that across any segment, upper and lower
bounds on the values of o?(v) and 6(t) are available. If they are not, cheap bounds
may be computed by sampling the function at the midpoint of the segment, and
70
using bounds on the first derivative to determine the maximum range of the func-
tion within the segment. Using this method carries an additional benefit. The
sample points may be used to compute a standard piecewise constant approxima-
tion to the intensity, which we call Iffia, b). We have found that Ifta, b) serves as a
better estimate to I(a, b) than I(a, b), even though mathematically the margin of
error is larger.
The method of piecewise constant bounds is valuable because it enables us
to compute a reliable error term for the piecewise constant approximation. In
Figure 4.2, we show the volume rendering integrand, and the approximation to
the integrand which is used for the piecewise constant technique. Recall that the
name of the method refers to how ?(v) and ?([) are estimated. The approximation
to the integrand itself is not piecewise constant.
The method of piecewise constant bounds can be used to compute an approx-
imation and remainder term once a partitioning of the domain of integration has
been established. It is not truly a controlled-precision technique, since it currently
lacks an algorithm for determining in advance a partitioning, given a desired level
of accuracy. For situations where an image must be computed to a prescribed
precision, the following methods are more suitable.
4.1.3 Trapezoid and Simpson's Rule
The trapezoid rule and Simpson's rule are two instances of Newton-Cotes quadra-
ture formulse [SB8O,Act7O]. These rules generate a polynomial spline approximant
to the integrand, f(t). The integral of the spline is easily computed, and is used
71
0.8
0.2
Integrand
Piecewise Constant, n=2
Piecewise Constant, n=8
0			I			I ________ I			I _________
0			0.2			0.4			0.6			0.8
Figure 4.2: The piecewise constant approximation. The true integrand is compared
to the approximation used for the piecewise constant method with two and eight
segments. In this graph, the horizontal axis is distance along the ray, t.
as an approximation to lab f(t) dt. The spline is formed by evaluating f(t) at n t 1
evenly-spaced sample points, a, a t h, a + ..... . , a + nh b. If a linear spline
is used, the resulting Newton-Cotes formula is the trapezoid rule. If a quadratic
spline is used, the result is Simpson's rule. In the following subsections, we present
these formul? along with their corresponding remainder terms.
72
Trapezoid Rule
Using the trapezoid rule, we may write:
I(a,b) --H $bf([)dt
where T, the trapezoid approximation, is defined by
f(a)+f(b) n-1
2 +Sf(a+ih)
i=1
(4.14)
T--H--Hh			(4.15)
and TT, the remainder, is defined by
b--Ha
--H			ftl(c) h2			(4.16)
12
for some c between a and b. Approximations to a sample volume rendering inte-
grand using the trapezoid rule with n 2 and n --H 6 are shown in Figure 4.3.
Because c is unknown, ?T cannot be computed exactly. However, using the results
from Section 4.1.1 we can compute an upper bound on fit (c) . Substituting this
bound into Equation (4.16) yields a bound on the magnitude of TT.
To obtain a solution to the volume rendering integral to within a given tolerance
we set ?? ? in Equation (4.16) and solve for the appropriate step size It.
This approach is perhaps the simplest method for controlled-precision integration,
requiring very little computation in addition to the basic trapezoid rule. Pseudo
code is given in Figure 4.4.
73
0 I
0.2
0			0.2			0.4			0.6			0.8
I(a,b) --H fbf(t)dt			(4.17)
--H S?rs			(4.18)
where S, the Simpson's rule approximation, is defined by
h
3?3
Figure 4.3: The trapezoid rule. The true integrand is compared to the approxi-
mation used for the trapezoid rule with two and four segments. In this graph, the
horizontal axis is distance along the ray, t.
Simpson's Rule
Using Simpson's rule, we may write:
0.8
Integrand
Trapezoid Rule, n=2
Trapezoid Rule, n=4
74
TrapezoidRule(p,q,a,b,?)
[Given polynornials p and q, compute the volume rendering integral
from a to b within tolerance 6 using the trapezoid rule.
P ?Upper bound on f"(c) for c C [a, b]
?			126/ ((b --H a)P)			[h'is largest step size that will ensure 6]
n ?			--H a)/h'i			[Divide range into regular steps of size ?
h ?			--H a)/n			[h is the step size for this partition]
Compute T using Equation (4.15)
return T
Figure 4.4: Pseudo-code for the trapezoid rule algorithm.
2
and r?, the remainder, is defined by
tf(a)+f(b)
b--Ha
--H 180 f(4)(c)h4			(4.19)
for some c between a and b.
Here n must be even. Approximations to the integrand using Simpson's rule
with n 2 and n --H 4 are shown in Figure 4.5. As with the trapezoid rule,
Equation (4.19) may be used to compute a bound on the approximation error for
a given step size h, or to compute the step size necessary to ensure a given accuracy.
The pseudo-code in Figure 4.4 may be used with only slight modification.
4.1.4 Power Series
We have developed another method for controlled precision volume integration
based on a power series expansion. It is obtained from Taylor's theorem coupled
75
0			I			I			I			I
0.2
with Lagrange's form of the remainder:
0			0.2			0.4			0.6			0.8			1
Figure 4.5: Simpson's rule. The true integrand is compared to the approximation
used for Simpson's rule with two segments. In this graph, the horizontal axis is
distance along the ray, t.
0.6
n f(k)(0) k
fl			= kA			k!
(4.20)
Theorem 4.2 (Taylor-Lagrange) Let f be a function that is C?+' on the in-
terval [0,1]. Then f(t) can be expressed as the sum of an n-term power series Pn(t)
and a remainder term 1?n(1),
f(t) = Pn(t) + ?n (1)
0.8
Integrand
Simpson's Rule, n=2
where
(4.21)
76
1?n(t)
+ 1)!
fof some c ? [0, ti.
(4.22)
Proofs of this theorem may be found in most introductory calculus texts, for ex-
ample the one by Thomas and Finney [TF79i. Approximations to the volume
rendering integrand using a two term and 10-term power series are shown in Fig-
ure 4.6.
We apply Theorem 4.2 to the problem of approximating I(a, b) by letting f be
the volume rendering integrand and writing
I(a, b)			L?(a, b) + fb Rn(t)dt.
where the n-term approximation 1n (a, b) is given by
I?(a,b) $6 Pn(t)dt.
Because Rn is bounded above for all n, and
?l\\? I?n(t) 0
(4.23)
for all t, as shown in Appendix B, the integral of Rn in Equation (4.23) must also go
to zero. Thus, as n increases, I?(a, b) provides an arbitrarily close approximation
to I(a, b). We derive the series and remainder terms for the volume rendering
integral in the following subsections.
77
0.8
Integrand
Power Series, n=3
Power Series, n=8
0.6
0.2
0o+			I
0			0.2			0.4			0.6			0.8
Figure 4.6: The power series approximation. The true integrand is compared
with its three-term and eight-term Taylor series approximation. In this graph, the
horizontal axis is distance along the ray, t.
The Series Term
The first step in forming the series term is to substitute the expression for j(n),
derived from Theorem 4.1, into Equation (4.21), giving
k q(i)(o)Ak?j(O) tk			(4.24)
k--H--HO j--H--HO
A0 (0)			(4.25)
An (0)			(4.26)
1 p(k)(0)An?k(0)
?k--H1			(k--H1)!
where
78
Since p and q are polynomials, these equations may be simplified. Let p and q be
defined by
p(t) --H potpit$... +PNtN
M
--H--H qo t qi t.... t qM t
For any polynomial A defined as in Equation (4.8), we have
n!P? ff n<D
0			otherwise.
Therefore, Equations (4.24)-(4.26) yield
n min(k,M)
0
ePo
Q Ak%(0) 1k			(4.27)
A0 (0)			(4.28)
1 rnin(n?N)
An (0)			--H ? kPkAn?k(0).			(4.29)
k--H1
Finally, integrating Pn(t) we obtain
min(k?M)			--H			(4.30)
In(a,b)--H--HZ ? q)A??j(O)
k--H--H0 j--H0			k+1
which approximates I(a, b). Pseudo code for computing I?(a, b) is given in Figure
4.7. Since M and N are fixed by the degrees of p and q, the complexity of this al-
gorithm is 0(n). The accuracy of this approximation is provided by the remainder
ten?.
The Remainder Term
In this section we derive a bound on Rn (t), the remainder term introduced in
Theorem 4.2, when the integrand f is given by Equation (4.4).
79
Approximate?I(p, q, a, b, n)
[Given polynomials p and q, compute the volume rendering
integral from a to b via an n-term power series.1
A0 ? e?0
`4 ?q0A0(b--Ha)
for k ? 1 to n do
Ak ? Gompute?A(p, k, A0,.. . , Ak?1)
min(k,M)
? QAk (bk+l --H ak??1) /k t I
j--H--H0
14 ? In t
end for
returnin
Gompute?A(p, k, A0,... Ak?l)
and terms A0... Ak?1, compute Ak as
[Given pol?nomial p
in Equation (4.29).]
Ak?1 min(k?N)
k
return Ak
IPj Ak?j
Figure 4.7: Pseudo-code for computing `4(a, b).
80
As outlined in Section 4.1.1, a bound on f(?)(t) for t ? [a, bj can be obtained
from Equations (4.5)-(4.7). We begin by substituting this bound for f(n?l) into
Equation (4.22), giving
min(n???1?M) ?q(i)An+1?j
Rn(t)			t?+1			(4.31)
where q(i) is the upper bound on the jth derivative of q from Equation (4.9). An+i?j
is an upper bound on An+i?j computed via the recurrence relation
A0 --H			c?			(4.32)
1 min(n,N) P(k)An?k
--H			(4.33)
k--H1 (k --H 1)!
The remainder term for the fl-term power series is formed by integrating l??(t)
from a to b,
min(n+1,M) ?q(i)An+i?j
fb 1?n(t)dt --H
Finally, using Equation (4.23),
--H a?+2 min(n+1?M) ?q(i)An+i?j
II(Qb)--Hi4(Qb)I ? nt2
Pseudocode to compute I(a, b) --H A?(a, b)I is given in Figure 4.8.
The form of the expression above does not provide a means of computing a pm-
ori the number of terms necessary to compute I(a, b) to attain a given accuracy ?.
Therefore, to use the power series, we accumulate terms while simultaneously eval-
uating the remainder term. The loop exits when the desired tolerance is achieved.
Termination is guaranteed by the convergence proof in Appendix B. The loop
may be coded efficiently by precomputing bounds on the derivatives of p and q.
81
Compute?Remainde?p, q, a, b, n)
[Given polynomials p and q, compute a nemaThder te?m,
bounding I(a, b) --H i4(a, b) .J
Compute p, bounds on the derivatives of p [See Section ?.1.1.?
Compute q, bounds on the derivatives of q
for k ? 1 to (n + 1) do
-Ak ? Compute?A?Bound(p?, k, A0,... , Ak?1)
end for
b?+2 --H a????2 rnin(n??1?M) ?q(i)-An+i?j
_			z
r m			__
n+2
return T
Gompute?A?Bound(p, k, A04k-1)
[Given bounds on polynomial p and te?ms-A0... Ak?1,
compute -Ak as in Equation (4???)?1
.?l minfflk,N) Pj Ak?j
kk o+?
return Ak
Figure 4.8: Pseudo-code for computing a remainder term.
82
PowerSeries(v, q, a b, ?)
[Given polynomials p and q, compute the volume rendernng
integral from a to b with tolerance ? using a power series.]
Compute p, bounds on the derivatives of p [See Section 4.1.1.]
Compute q, bounds on the derivatives of q
A0 ? e?0
A0
14 ?q0A0(b--Ha)
k?1
repeat
Ak ? Compute?A(p, k, ..... . , Ak?l)
Ak ? Compute?A?Bound(p, k, ..... . , Ak?1)
min(k?M)
In ? In + ? qiAk_ - (bk+l --H ak+l) /k + 1
j--H0
r ? b?+2 --H a?+2 min(n?1?M) q(i)An+1?j
n+2			j--H--H0
until r ?
return I
Figure 4.9: Pseudo-code for the complete power series.
Iterative techniques are then used to evaluate An(t) and An, as shown in pseudo
code in Figure 4.9.
4.1.5 Discussion and Results
In order to evaluate the relative performance of the controlled-precision integration
methods, each of the approaches were tested using a random distribution of rays
passing through a single voxel whose emission and absorption values at each corner
were randomly distributed over the interval [0, 1]. The time to compute 10,000
83
such integrals to a specified precision was measured. The results are summarized
in a Figure 4.10, where computation time is plotted against precision. Note that
precision is measured in bits, so the horizontal scale is logarithmic.
800
700
? 600
500
0
? 400
0
? 300
Simpson's Rule
Trapezoid Rule
200____
4			8			24			28
12			16			20
Bits of Precision
32
Figure 4.10: Performance on 10000 test integrals.
Figure 4.10 shows that the three quadrature methods exhibit very different
convergence properties, and that each method has a range where its performance
is superior. For low- and medium-precision results, trapezoid and Simpson's rule
are faster than the power series, although the relative difference in computation
time is small. Trapezoid rule convergence is only quadratic, while Simpson's rule
is fourth order. However, the performance of the trapezoid rule is superior at low
precision because it needs to bound only the second derivative of the integrand,
in contrast to Simpson's, which requires bounds on the fourth derivative. The
84
power series requires several terms before beginning to converge, making it less
efficient for lower precisions. Convergence is exponential, however, which enables
it to compute very high precision solutions that would not be practical with the
other methods.
Since most computer graphics displays have an 8-bit dynamic range, one might
conclude from the graph in Figure 4.10 that Simpson's rule is the method of choice,
This is not always the case, however. First, the precision reported in Figure 4.10 is
absolute, whereas most volume-rendered images are scaled relative to the brightest
pixel value in the image. Relative scaling will shift the curves along the precision
axis, thus changing the location of the crossover points. Secondly, the precision
reported in Figure 4.10 is for a single voxel. As described in the next section
when the precision along an entire ray is considered, the errors from each voxel
accumulate; thus, the tolerance for a single voxel may have to be considerably
tighter than the tolerance for the ray itself.
The methods described in this section compute guaranteed bounds on the er-
ror, but not necessarily the tightest bounds possible. When a given solution is
significantly more accurate than can be proven with these tools, computation time
is wasted in refining the results until they are provably below tolerance. There is a
tradeoff in utilizing tighter bounds: they are more expensive to compute, but they
may allow the quadrature algorithm to compute a result to the same precision
with less work. Investigation of other methods of bounding these derivatives is
warranted.
Further work is needed on the applicability of other quadrature rules, such as
85
Table 4.1: Quadrature rules and their information requirements. The function f(t)
is the volume rendering integrand, e?j7Q(v) d? 6(t).
Piecewise			Trapezoid,			Power
Constant			Simpson's			Series
Conditions			Bounds			Rules			Method
Bounds on a'(v), 6(t)			X
f(t) Computable			x
f(n) (0) Computable			x
Bounds on f(?)(t)			x			x
Gaussian quadrature, to the volume rendering integral. The power series could
also be converted to a continued-fractions format, in hopes of achieving better
convergence properties [PFTVSS].
While the derivations we have presented are specific to polynomial ? and 6, the
underlying quadrature rules are more general. For example, Newton-Cotes meth-
ods require only evaluation of the integrand at points and bounds on derivatives in
the domain of integration. A summary of the requirements for each algorithm is
given in Table 4.1. This work should be extended to exploring more general forms
of c? and 6.
86
4.2 Controlled-Precision Ray Casting
Since we are assuming polynomial spline interpolation, the controlled-precision
quadrature techniques described above can only be used for computing values of
I(a, b) within a single cell. The trapezoid rule, Simpson's rule and the power
series method require bounds on derivatives to compute remainder terms. These
bounds are not available across cell boundaries because the interpolated function
is discontinuous. In this section, we consider how these quadrature rules can be
used in a simple ray-casting volume rendering algorithm to compute I(a, b) across
multiple cells to any specified precision. In Chapter 5 we will develop a more
efficient hierarchical algorithm.
We begin by partitioning the ray along its intersections with the voxel grid.
Where these intersections correspond to discontinuities of ? and ?, this strategy is
akin to the discontinuity meshing algorithms for radiosity [LTG92]. The algorithm
we use for partitioning the ray is the one developed by Snyder and Barr for ray-
casting an axis-aligned rectilinear grid [SB87J. The grid spacing must remain
constant along each axis, but may vary between the axes. Their algorithm is
based on the following observations:
1. All intersection points occur on an x, y, or z plane of the volume grid.
2. A ray intersects each of these planes at regular intervals.
Let t be the distance parameter along the ray. During the initialization phase, the t
values of the first intersection of the ray with the bounding x, y and z planes of the
data are computed and stored as tx, Jy, tz. The values Ax, Ay and Az are set to the
87
distance along the ray between data slices along each of the axes. The algorithm
then iterates by choosing the smallest of fix, ty, [zi as the next intersection with
grid. That value is then incremented by the appropriate A value. The algorithm
therefore identifies all voxels in a ray's path and the distances along the ray to
each of these intersections. The algorithm requires no multiplications or divisions
in the inner loop. An iteration of the this algorithm is shown in Figure 4.11.
Figure 4.11: Tracing a ray through a grid using the Snyder and Barr algorithm
[SB87]. The ray enters the cell at t0. The next intersection points with a grid
plane along each of the axes are tx, ty and tz. The smallest of these values, in this
case tx, is the ray's next grid intersection. The algorithm proceeds to the next cell
and increments tx by Ax.
Recall from Section 4.1.2 that if intervals of intensity, 1k, and transparency,
Tk, are given for each segment k, then bounds on I(a, b) can be computed using
interval arithmetic,
n			i--Hi
I(a,b)			?			Ii			fi 7).			(4.34)
i--Hi			j=1
We can obtain 1k using the quadrature techniques from the previous section by
adding and subtracting the remainder terms from the estimates. Since absorption
88
is polynomial, the transparency integral reduces to a simple arithmetic expres-
sion. Neglecting roundoff error, we can assume that Tk can be computed exactly.
Therefore we write it as an ordinary scalar, Tk. If the user wishes to consider
roundoff error, the more sophisticated approach detailed in Chapter 5 can be used
to properly handle intervals of Tk.
To compute a ray with an error tolerance of ?, we may distribute the error
evenly across all the segments allocating a tolerance
fl
to each of the additive terms in Equation (4.34),
i--Hi
`2117).
j=1
Since each of the Tk's are in the range [0,1], the value of the product will be a
fraction between zero and one. Thus, the error tolerance for the ?th term may be
met by computing 1k with a tolerance of
74%..74-i
which will be larger than or equal to ?. A front-to-back algorithm that implements
this strategy is shown in Figure 4.12. Line (1) is an application of a quadrature
rule such as the ones introduced in Section 4.1. The choice of quadrature rule can
be based on the value of ? to obtain the fastest convergence for the given tolerance.
The algorithm developed in this section was used to compute the 256 x 256
image of the pelvis dataset shown in Figure 4.13a. Each ray was computed to the
nearest grey level. For comparison, Figure 4.13b shows the same scene rendered
89
CaSt?Ray(rffla,6,?)
[Compute I(a, 6) along a ray ? to with?n ??.?
T = 1 [Initialize transparency]
I = 0 [Initial?ze ?ntens?ty]
Partition T into N ray--Hvoxel intersections
for each ray--Hvoxel intersection i from front to back do
Compute Tj			[Equation (i.2)]
Compute I? to error tolerance ?			[Section (4.1)]
I ? I t j? x T
T? Tx Tj
end for
return I
Figure 4.12: Pseudo-code for controlled-precision ray casting,
90
using the piecewise constant approximation. Piecewise constant regions were equal
in length to half the voxel grid spacing. The images appear quite similar. Fig-
ure 4.13c is a subtraction and rescaling of the Figures 4.13a and 4.13b. Neutral
grey indicates an exact match to 8-bit precision. The brightest and darkest areas
indicate a difference of 20 grey levels a significant disparity. When the step size
of the piecewise constant algorithm was set to one-tenth of the grid spacing, the
two images agreed to 8 bits of precision. A radial beating artifact in the output
of the piecewise constant algorithm is enhanced in Figure 4. 13c. This is aliasing
caused by undersampling along the ray axis. The coherence of the artifact arises
from sample points that are placed at uniform distances from the eyepoint.
4.3 Integration over pixel area
The approach described in the previous section correctly computes the value of
I(a, b) along individual rays to any specified precision. Recall from Section 3.3.1
that ray casting point-samples the illumination function across the image plane.
By itself, guaranteeing the accuracy of each ray does not guarantee the accuracy
of the final image. This is due to the fact that we have no bound on the variation
of the illumination between the sample points. Error due to point-sampling of the
image plane can be severe. If an object falls between rays, it will not contribute
to the final image at all.
Our preliminary investigations into bounding screen-space error concentrated
on obtaining bounds on the illumination function across screen-space disks of a
given radius ?. A practical solution to this problem enables a controlled-precision
91
(a)
(c)
Figure 4.13: Controlled-precision ray casting. Image (a) was computed using con-
trolled-precision techniques to within one grey level at each ray. Image (b) was
computed using the piecewise constant approximation with a step size of one-half
the maximum grid spacing. Image (c) is a subtraction of images (a) and (b),
highlighting the errors in image (b)
92
algorithm, since a screen-space region could be approximated with smaller and
smaller disks until a specified error tolerance is reached. Our initial approach
operated by sampling absorption and emission at the center of the disk and con-
sidering the maximum variation of the functions across the requisite volume. The
resulting algorithm is impractical, since the bounds computed in this manner are
far too coarse.
Obtaining tight bounds for the variation of the illumination function is difficult
because there are no guarantees that the function will be slow-varying. For exam-
ple, two rays in any finite neighborhood may intersect a completely different set
of voxels. A tempting approach would be to consider all possible discontinuities of
the illumination function by projecting the voxel grid onto the screen. However,
this technique would be impractical because it would subdivide the screen into
0(n3) individual regions, where n is the linear resolution of the data.
4.4 Summary and Future Directions
The work presented in this chapter represents a step toward more rigorous error
analysis of volume rendering techniques. NVe first presented controlled-precision
techniques for computing the emission-absorption model along rays that pass
through a single interpolation cell. These results were extended to form an algo-
rithm that computes illumination along rays passing through collections of interpo-
lation cells. A theoretically valid algorithm has also been presented for controlled-
precision calculation of pixel values.
There remains much work to be done in this area. A practical algorithm for
93
bounding error in screen-space is an important next step. There are also other
sources of error in volume-rendered images that we have not touched upon. These
include reconstruction error and uncertainty in measured data values. Errors in-
troduced by approximations such as fixed-precision arithmetic and splatting have
yet to be explored.
Throughout this chapter, we have been assuming an objective error metric,
one that measures deviation from an ideal image. While these metrics are easy
to fit into a mathematical framework, their correspondence to how we perceive
error in images is unknown. Study of perceptual effects of rendering error on the
interpretation of volume-rendered images is warranted.
The abstract nature of volume rendering poses some difficulty in performing
any type of error analysis. Cue abstraction is the process of mapping raw data
values to functions of emission and absorption. This process is currently under
user control. Clearly errors can occur during this assignment phase. But since we
have no ideal model of how mapping should take place, we are currently unable to
quantify errors introduced at this stage.
Chapter 5
Efficiency Issues
Volume rendering is often used in tim&critical applications, such as medical imag-
ing, that require fast imaging of data. Even in situations that don't require im-
mediate access to final results, fast-feedback images are often helpful when testing
rendering parameters such as camera position, lighting and color. While improve-
ments in processor speed and memory chip capacity have resulted in a continual
increase in the power of graphics workstations, a simultaneous explosion in dataset
size has more than kept pace. Advances in scanning technology result in an 0(n3)
increase in data size for every 0(n) improvement in linear resolution. Algorithmic
improvement is necessary to continue to provide users with reasonable rates of
feedback.
While we wish to provide fast feedback, we don't want to compromise on our
goals of presenting information clearly and accurately. Feedback that confuses or
misleads the user is certainly not an asset. In this chapter, we consider issues
94
95
regarding efficient computation of volume-rendered images In the first section, we
address techniques for fast computation of perspective projections of volume data.
This is important since many of the faster algorithms presented in Section 3.3
could only compute orthographic projections. In the second section we introduce
an adaptive error bracketing scheme that attempts to compute ray values to any
specified precision using minimum computation.
5.1 Efficient Perspective Projection
The majority of images produced using volume rendering have employed ortho-
graphic projection. While quite effective in many circumstances, orthographic
views have, in a sense, held scientists at an infinite distance from their data. Or-
thographic projections often result in depth ambiguities and reversals of apparent
depth as is shown in the surface rendering of Figure 5.1. Since our visual system
is tuned to comprehend the perspective views projected onto our retinas, scene
geometry and depth are more readily apparent in simulated perspective views. In
addition, data fiythroughs and lens effects become available to aid in the compre-
hension of spatial relationships within the data when perspective is used.
The widespread use of perspective projection has been prevented by a lack of
efficient algorithms. The fastest rendering algorithms such as splatting and coher-
ent projection exploit the coherence of parallel viewing rays. Frequency-domain
rendering techniques are currently constrained to the parallel projection model
used by the Fourier slice theorem. R?ay-casting approaches have been hampered
by the need for random access to large datasets and by the nonuniform sampling
96
Figure 5.1: A comparison of a perspective (top) and an orthographic (bottom)
view of two equal-sized spheres connected with a cylinder. In the orthographic
view, lack of perspective depth cuing causes size and orientation ambiguity.
rate inherent in perspective ray divergence. In this section, we present an opti-
mal paging strategy and an adaptive supersampling technique that address each
of these problems. Used in concert, they result in an algorithm that computes
perspective projections with modest requirements in terms of both computation
time and resident memory.
Several researchers have explored the use of perspective in volume rendering.
However, to date none have described efficient, single pass slice-by-slice processing
techniques. Most ray casting implementations have assumed random access to the
volume dataset [SSW86,UN88,SabS8,Web9O]. If the volume data has n voxels on
a side, this approach requires 0(n3) storage space, which quickly outstrips RAM
capacity for high resolution arrays. Random access to large virtual memory spaces
97
(a)			(b)
Figure 5.2: Computing perspective by predistorting the volume. The perspec-
tive rays, (a), can be computed by firing parallel rays, (b),into an appropriately
distorted version of the dataset.
causes an unacceptable amount of disk paging. Another approach, Upson and
Keeler's cell-by-cell processing method, assumes that the data can be accessed
a slice-at-a-time along any of the three coordinate axes [UN88]. This presents a
difficulty because volume data is usually stored so that slice data may be efficiently
accessed on only one axis.
Perspective views may also be realized by casting parallel rays into an ap-
propriately transformed copy of the volume grid [DCH88,HBP+89,Han9O]. This
approach is illustrated in Figure 5.2. Aside from being expensive, the resampling
step is prone to error. In fact, even under ideal conditions, where the original
volume is sampled above the Nyquist frequency, error can be introduced using this
method. This is because the data is resampled onto a grid with the same resolution
as original grid. The shrinking of parts of the volume that are distant from the
eyepoint, an effect that can be observed in Figure 5.2, will cause an increase in
the frequency of the 3D signal. This in turn raises the signal's Nyquist frequency,
possibly rendering the original sampling rate inadequate [Han9O].
98
0
100000
10000			I			-			---
1000
100
0.1
512			1024			1536			2048
Data Resolution
Figure 5.3: Storage required for a 3D volume compared to storage required for two
2D slices. Even at the highest data resolutions, the 2D slices will fit into RAM
memory on current workstations. A 4-byte floating-point representation for each
sample of ot and 6 is assumed. Note that the vertical axis is logarithmic.
5.1.1 An Optimal Paging Strategy
In this section, we present an algorithm that reduces the resident memory re-
quirements of a ray-casting volume renderer from 0(n3) to 0(n2). The increase
in efficiency is due to processing two 2D slices of data at a time, rather than the
entire dataset. The significance of this improvement may be observed in Figure 5.3
which compares an 0(n3) growth curve to an 0(n2) growth curve. Assuming 32-
bit floating point representations for each sample of emission and absorption, a
2048 x 2048 x 2048 dataset requires 68 g?gabytes of data. Two 2D slices of this
same dataset will require only 67 megabytes of data storage that is easily found
99
on modern workstations. We present the new algorithm in two parts. First, a
basic scheme to allow ray tracing on a slice-by-slice basis is introduced. Next, we
describe a processing order that requires that each slice of the dataset is loaded
into memory only once.
t
A
A
A
A
(a)
Figure 5.4: Ray-oriented versus data-oriented processing. In the diagrams above,
ray segments with the same shade dot are processed during the same iteration.
(a) Using ray-oriented processing, all segments on each ray are processed at the
same time. (b) Using data-oriented processing, all segments in the same region are
processed at the same time.
Slice by Slice Ray Tracing
A typical ray tracer has an outer loop that iterates once for each ray. Each ray is
traced until it terminates, yielding a final pixel value. Processing then concentrates
on the next ray. This scheme, illustrated in Figure 5.4a, is extremely inefficient
with regard to memory access because each ray requires access to sample points
from across the entire dataset. If the data is too large to fit into RAM, the
processing of the segments that comprise a ray is likely to cause virtual memory
loo
page faults. The faults will be repeated for each ray, and each part of the dataset
may be loaded and re-loaded into RAM thousands of times.
To increase efficiency the algorithm is inverted, so that it is driven by the
data currently in memory. This technique, illustrated in Figure 5.4b, performs
as many of the necessary operations on the data currently loaded into memory
before swapping it out of RAM. We decided to base our paging strategy on two-
dimensional slices of data. As we have already established, the growth rate of n2
is slow enough that slices can be accommodated in typical memory configurations.
At the same time, slices are large enough so that the overhead involved in reading
a slice of data from a disk is relatively insignificant. Because we rely on trilinear
interpolation at sample points, at least two slices of data are required in memory
at one time. We call the volume bounded by a pair of adjacent slices a slab.
Higher-order interpolation can be used at the cost of maintaining thicker slabs.
In?t?alize?SlabQueue(r, V, s)
[Gene?ate an ?n?t?al slab queue, g?ven a set of rays, r, a volume V,
and a set of slabs, 5
SlabQueue ? empty slab queue
for each ray rj in r do
? slab that contains the first intersection of r? with V.
Enqueue ray r? onto SlabQeueue[s?]
end for
return SlabQueue
Figure 5.5: Pseudo-code for initializing the slab queue.
In order to implement the data-oriented strategy, we developed a data structure
101
called a slab queue. A slab queue has an entry for each slab of data. Each entry is
a queue of rays. These are rays that require the data in the slab for the processing
of their next segment. The slab queue structure is initialized by tracing each
ray to its first intersection with the volume array. The slab that contains this
intersection point is determined, and each ray is placed on the appropriate queue.
This algorithm is shown in pseudo-code in Figure 5.5. The existence of the slab
queue enables the following data-oriented processing strategy, given in pseudo-code
in Figure 5.6. Whenever a slab is loaded into memory, it activates each ray on its
queue in turn. The ray processes the segments in its path until the next segment
lies in another slab, or until it exits the volume. In the former case, the slab
containing the next segment is determined, and the ray is placed on the proper
slab queue.
This algorithm requires the ability to stop and restart the tracing of a ray. We
therefore introduce another data structure that contains all information pertinent
to a ray's processing. A pointer to this data structure is what actually is stored in
the slab queue. At a minimum, this structure must contain:
o+ For each wavelength of interest, the intensity accumulated along the ray up
to the current segment;
o+ the transparency accumulated along the ray up to the current segment;
o+ the direction of the ray (the origin is simply the eyepoint).
o+ a pointer to the next ray structure on the slab queue, or an indicator that
this ray is the last on the queue;
102
P?ocess?SlabQueue(SlabQueue)
[Given an initialized slab queue, compute a volume rendering]
while SlabQueue is not empty do
Choose a slab s? to load into memory
for each ray r? in SlabQueue[s?] do
repeat
Evaluate illumination for current ray segment [Section 11]
? slab containing r?'s next segment
until s? ?
if s? is within the volume then
Enqueue r? onto SlabQueue[s?]
end for
end while
Figure 5.6: Pseudo-code for processing the slab queue.
Additional fields are useful for storing information that speed ray operations. For
a piecewise constant algorithm that takes evenly-spaced sample points along the
ray, such as the one presented by Levoy [Lev88], retaining the current sample
point reduces ray reactivation overhead. Storing the last sample point that is in
the volume speeds testing of the termination condition.
In order to apply the quadrature rules given in Chapter 4, the exact boundaries
of the ray segments must be computed. For this purpose we use the ray/grid
intersection of Snyder and Barr that was described in Section 4.2 [SB87]. For this
algorithm, it is helpful to store distance, t, from ray origin to the start of the
current segment and the t values of the next data slice along each axis: tx, ty, tz.
The exact amount of storage will depend on the bits of precision used to store
103
numeric values and on the optional fields that are used. The current implementa-
tion for computing single-precision grayscale images uses 40 bytes for evenly-spaced
sample points and 56 bytes for exact ray/grid intersection. Total storage cost for
the ray data structures and one slab of emission and absorption values is 56rn2t8n2
bytes of storage, where n? is the screen resolution and n is the data resolution. This
compares to 8n3 bytes of storage when the entire dataset is stored in memory. If
we assume screen resolution equal to data resolution, the method presented here
is more memory efficient for n > 8.
The slab queues are maintained as linked lists through the pointer field in the
ray data structure. The ray queuing and dequeueing operations require only two
pointer reassignments each, since objects are added and removed from the front of
the queue. Thus, the overhead for maintaining the queues is negligible.
Slab Processing Order
Now that we have established a method for data-driven evaluation of ray segments,
we wish to determine a slab processing order that allows all rays to be computed
while loading each slab into memory only once. We derive the processing order for
front-to-back processing of ray segments. The reverse order is optimal for back-to-
front processing. The order we define takes advantage of the fact that all viewing
rays share a common origin --H the eyepoint. Assume without loss of generality
that the data may be accessed as slices perpendicular to the axis. Let k be the
number of slabs, and let the slabs be indexed by increasing z value, Si, S2,???, 8k
We denote the slab that contains the coordinate of the eyepoint s?. If no slab
104
contains the z coordinate of the eyepoint, 8 is the slab whose z range is closest the
the z coordinate of the eyepoint. Regardless of view direction, all rays emanating
from the eyepoint must fall in to one of three categories:
1. The ray may stay in the Sz slab.
2. The ray may diverge in the +z direction from the eyepoint. In this case it
requires access to the slabs 8?, 8z+1,?? 8k in increasing order.
3. The ray may diverge in the --Hz direction from the eyepoint. In this case it
requires access to the slabs 8?, 8z--H1,?? si in decreasing order.
Rays in any of the categories may therefore be processed during a single pass
through the dataset. We first load slabs s? through 8k in increasing order and then
slabs 8z-1 through 8i in decreasing order. Alternatively the slabs may be loaded
starting with s? and decreasing first. In either case, the order is optimal in terms of
paging since each data point is loaded only once. Several examples using different
eyepoints and processing orders are shown in Figure 5.7.
Combining the data-driven processing technique with the slab processing or-
der, we obtain a complete algorithm for ray-tracing volume datasets using per-
spective projection that requires only o(rn2) storage, makes a single pass through
the dataset, and requires only a negligible amount of computation time over that
required by algorithms requiring 0(n3) storage. Adaptive refinement may be
achieved by displaying the image after each slab is processed. The effectiveness
of the intermediate results depends on the orientation of the slabs relative to the
image plane.
105
Sz
(a)
Sz
(b)
Sz
(c)
Figure 5.7: Slabs are processed starting from slab 8z, the slab containing the
z-coordinate of the eyepoint, and continuing outward. This order is indicated by
the arrows at the bottom of each diagram. In (a) and (b), the z-coordinate of the
eyepoint is the same. This leads to identical processing orders even though the
view directions are different. In (c), the eyepoint lies outside the z extent of the
slabs. The slab closest to the eyepoint becomes s?.
106
5.1.2 Sampling Issues
We now turn our attention to sampling problems inherent in perspective projection.
There are two issues:
o+ Using perspective projection, objects that are more distant from the eyepoint
occupy less screen space than similar-sized objects that are close to the eye-
point. Therefore, there is a loss of object-space resolution toward the back
of the scene being viewed.
o+ Perspective rays diverge, causing the 3D sampling rate to change with depth.
Undersampling of the volume data is likely to occur.
The first issue is inherent in a perspective projection. This effect is familiar in
conventional rendering and photography as well as in our own visual system. We
do not consider this a serious problem. If increased resolution is necessary for
objects near the back of the dataset, the user may move the eyepoint closer to the
object of interest, rotate the eyepoint to look at the back of the dataset, or simply
increase image resolution.
Of greater concern is the potential for undersampling the volume data due to
perspective ray divergence. This can result in severe aliasing, which is visually
unappealing, and more seriously, the erroneous omission of small objects from the
image entirely. To better understand this latter case, consider a 60-degree field
of view looking directly onto a face of the volume bounding box. One ray is shot
through the center of each pixel on a screen equal in resolution to a slice of the
dataset. The first slice of the dataset will be sampled at the same rate as the
107
screen. However, this sampling rate will decrease for each successive slice. By the
time the last slice is reached, the sampling rate is half the screen sampling rate.
This example is illustrated in 2D in Figure 5.8.
One way to control undersampling is to increase the number of rays that are
fired per pixel. Such supersampling was the basis for the controlled-precision image
rendering algorithm presented in Section 4.3. With appropriate error tolerances
and sufficient computational resources, this technique can compute perspective
images with no perceptible rendering artifact. Uniform supersampling of screen
space is not, however, an efficient method for increasing the accuracy of perspective
images. As shown in Figure 5.9, when we fire enough rays to adequately sample
the regions of the dataset that are distant from the eyepoint, the regions close to
the eyepoint are ove?ampled, thus wasting computation time.
An alternative approach to maintaining adequate sampling rates with perspec-
tive projection would be to use multiresolution filterings of the data as in the
3D-mipmap algorithms [Wil83,Lev90a]. In these schema, instead of firing addi-
tional rays to obtain the average value across a larger region, samples are taken
from successively blurred versions of the volume data. This has the positive ef-
fect of reducing the number of samples taken during ray tracing. In particular, it
reduces computation effort in distant areas which do not contribute much to the
final image. However, it adds an extra level of interpolation between the raw data
and the ray sampling process that is likely to introduce error. It also adds some
complexity to the ray tracer which will have to pass through several sets of slabs.
108
Eyepoint
Screen
Figure 5.8: Divergence of perspective rays causes nonuniform sampling of the
volume data. In this case, although the first slice of data is sampled once per
cell, the last slice is sampled only every two voxels. The shaded cells are not
encountered by any rays. Any data contained in these cells will not contribute to
the final image.
109
Eyepoint
Screen
Figure 5,9: A simple solution to the problem illustrated in Figure 5.8 is to use pixel
supersampling. Enough rays are fired through each pixel so that the most distant
voxels are sampled adequately. This results in wasted computation as voxels close
to the eyepoint are oversampled.
110
We have developed a supersampling approach that does not waste processing
time by oversampling regions near the eyepoint. The method originates supersam-
pling rays within the volume, rather than at the eye point. The basic scheme is
illustrated in Figure 5.10. A ray is traced through the center of a screen pixel,
and into the volume. Processing continues until the ray reaches a point where the
projected area of the pixel reaches a user-specified threshold. The ray is then split
into four supersampling rays that originate at the centers of subpixels, extrapo-
lated to the distance of the split point. The splitting continues recursively for the
supersampling rays as the projected area of their subpixels reaches the threshold.
The process stops when a ray exits the volume. The algorithm results in a tree of
rays spawned from each initial viewing ray. The overall effect is shown in 2D in
Figure 5.11.
For each ray, values of intensity, I(a, b), and transparency, T(a, b) are computed
using the quadrature rules from Section 4.1. After the splitting is complete, a pixel
value must be determined based on the supersampling ray values, which are stored
in the tree of rays. At each level of the tree, the four child rays are composited to
yield a single value of intensity, I, and transparency, T. These composite values
will be combined with those obtained along the parent ray using the methods
from Section 3.3.2. Composition of the ray values of the children is analogous to
the process of obtaining a single pixel value from ordinary supersampling. In the
current implementation, we use a box filter, however any filter could be used in its
place.
111
fy.			1.
Screen Pixel
Projected Pixe!
(a)
Supersampling
Projected			Rays
Subpixels
(b)
Figure 5.10: The ray splitting algoritlim. (a) A ray is traced until the projected
area of the screen pixel reaches a user-specified threshold. (b) At this point, the
ray is halted. Further processing is done along four supersampling rays. This
prevents the sampling rate from dropping too low.
112
Screen
Eyepoint
v
Figure 5.11: Using adaptive supersampling, each ray is split when divergence
causes sampling below a threshold value. The resulting tree of rays are then corn-
posited to yield final pixel values. Using the method, the benefits of supersampling
distant regions can be attained without the cost of oversampling nearby regions.
113
The size of the projected pixel that will induce splitting is determined by con-
sidering the smallest-sized object that the end-user wishes to be considered in the
final image. By setting the threshold projected pixel area to be less than the square
of the smallest voxel grid spacing, all voxels within the viewing pyramid can be
guaranteed to contribute to the final image. Reducing the threshold value will re-
sult in increased accuracy at the cost of more supersampling rays. Unfortunately,
it is not currently possible to determine in advance the splitting threshold that is
necessary to ensure that final pixel values are computed to within a certain error
tolerance.
This adaptive supersampling scheme works well in conjunction with the slice-
by-slice processing algorithm presented above. The ray split points do not depend
on the data values, so the start and end points of each ray in the tree may be
computed as a preprocess. Each of these rays is then queued separately during
initialization. The slice processing order remains unchanged.
The computation of pixel values from the ray trees is performed during a post-
process. In order to implement this phase of the algorithm, it is necessary to amend
the ray data structure described in Section 5.1.1 to include pointers to the ray's
children. If the child rays are allocated as a contiguous block of memory, only one
pointer needs to be stored.
Our adaptive supersampling method is attractive because it automatically com-
putes the number of child rays necessary to maintain a given sampling rate in a
manner that is independent of screen resolution. Through ray splitting, computa-
tional effort is not wasted by oversampling the region close to the eyepoint. While
114
supersampling is not always necessary for good image quality, this technique can be
implemented quite easily as an option to any ray tracer, thus providing a guarantee
that all voxels are taken into account when such precision is absolutely necessary,
115
5.1.3 Results
Figures 5.12 through 5.16 show 256 x 256 renderings of a 119 x 112 x 73 CT study
of a human pelvis. The viewpoint is on the left side and slightly toward the back of
the body. Image 5.12 compares a perspective view with an orthographic view. In
both cases, one ray was cast per screen pixel. Note that in the orthographic view,
the size of the femurs (leg bones) do not vary with distance, so the orientation
of the pelvis is difficult to discern. The perspective view was rendered using the
slice-by-slice processing technique presented here. It took 18 seconds to render on
a HP750 workstation. Profiling showed that 2.5 per cent of the total time was
spent performing the ray enqueue and dequeue operations that were described in
Section 5.1.1. This percentage could be substantially reduced if the calls were
implemented as macros rather than as procedure calls.
Figure 5.13 compares uniform supersampling to the ray-splitting adaptive su-
persampling technique. Figure 5.13a was computed using uniform supersampling
to ensure that at least two rays pass through each voxel. Approximately 12 mil-
lion samples of emission and absorption were taken inside the volume array. Fig-
ure 5.13b was computed using adaptive supersampling to ensure the same min-
imum sampling rate while taking only about 6 million samples of emission and
absorption inside the volume array. Figure 5.13c is a subtraction of images a and
Perspective projection also enabled us to experiment with alternate forms of
presenting the volume data. One approach was to attain data fiy-throughs by
moving the eyepoint within the volume. Three frames of a fly-through animation
116
(a)
(c)
Figure 5.12: (a) A 256 x 256 orthographic view of a 192 x 113 x 73 CT study of a
human pelvis. Note that the orientation of the pelvis is difficult to discern. (b) A
256 x 256 perspective view of the same dataset.
117
(a)			(b)
(c)
Figure 5.13: Uniform versus adaptive supersampling. (a) Uniformly supersampled
image that ensures at least 2 rays pass through each voxel. (b) Adaptively super-
sampled image that ensures at least two rays pass through each voxel. Computa-
tion time was approximately halved by using the adaptive ray-splitting technique.
(c) Subtraction of the two images.
118
sequence are presented in Figure 5.14,
We also experimented with depth of field effects [Pc82,cPc84l. If the user
defines an area of interest within the volume data, we found that one may improve
image comprehension by focusing a virtual camera on that region. The result is
better depth perception because objects are blurred as distance from the focal
region increases. The blurred remnants of distant objects help to maintain global
perspective, but do not interfere with the comprehension of objects in the focal
plane. One method for computing depth of field effects is is to distribute the rays
over a finite aperture to blur the image outside the focal region [CPC84]. Because
this technique effectively involves the jittering of the eyepoint, the slice processing
order must be modified to take into account the minimum and maximum z extents
of the eyepoint. This necessitates making two passes through the slices within this
extent, as illustrated in Figure 5.15. Three images with different aperture sizes are
presented in Figure 5.16. The images are focused on the base of the spine using a
simulated normal lens at f/2, J/5.6 and f/11
5.1.4 Discussion
There are several aspects of our research in efficient perspective projection that
warrant further research:
Alternate processing orders. The slice processing order presented here is op-
timized so that each slice of memory is loaded only once. One could imagine
alternative slab processing orders that would provide more useful intermediate
images for adaptive refinement. This technique could be useful in time-critical
119
(a)
(c)
Figure 5.14: Using perspective projection, the eyepoint may be moved inside the
volume. The images (a), (b), and (c) show successive frames of an animated
fly-through of the CT dataset.
120
Processing
order
Jt			__
w			r
extended eyepoint,? I
(lens surface) /, I			\ \
??i I
I			?
I			?
z
I			/			11			I			?
I,			?I			I			I
IIII? I
I I
Figure 5.15: When jittering samples over a lens surface, the processing order must
be modified to account for the shifting ray origins. Slabs between the minimum
and maximum z extents of the lens must be processed twice.
computing applications. If each ray were given an importance weight, processing
could be directed towards the slab with the most weight in its queue. A ray would
be assigned a high weight if, for example, it is known to pass through a region of
interest.
Parallelism. The slice-by-slice processing technique could be modified to take
advantage of shared memory coarse-grained parallelism by distributing the rays
on a given slab queue to several processors. If memory is not shared, all processors
must load the current slab into memory. Thus parallelism should only be used if
121
f/2			f/5.6
f/11
Figure 5.16: Images computed by focusing a virtual camera on the base of the
spine and using different f settings.
122
the slab queue is longer than a certain threshold.
Alternate data topologies. The ray tracing algorithm assumes that the data
is stored as parallel rectilinear slabs. Other regular data structures should be
amenable to the paging technique. Such data structures include octrees and data
sampled in cylindrical and spherical coordinate grids.
Controlled-precision techniques. The ray-splitting algorithm should be ex-
tended so that the ray split points can be automatically determined, given an
error tolerance for the final pixel value.
5.2 Efficient Controlled-Precision Ray Casting
We now consider the efficient computation of rays to any specified precision. Recall
that the technique presented in Section 4.2 had three major shortcomings:
1. The error tolerance for the ray was distributed evenly across all segments,
regardless of the ease of approximating that segment.
2. All voxels along a ray's path had to be visited.
3. The algorithm assumed that transparency could be computed exactly, and
that no error was introduced due to floating-point arithmetic.
Other researchers in volume rendering have investigated methods for trading
speed for accuracy in the computation of ray intensities. We described the most
popular of these approaches, early ray termination, in Section 3.3.3. Early ray
termination avoids computation in areas of the volume that are largely occluded
[Lev9Oa,DH92]. Danskin and Hanrahan introduce a more sophisticated approach
123
whereby coarse approximations to the volume rendering integral are employed in
segments where such approximations will introduce little overall error [DH92]. Such
segments include those that are mostly occluded; those that are empty or nearly
empty; and those that are largely homogeneous. The method employs an octree
organization of the data. At each node, the range of values encountered among all
child nodes is stored, as well as the average value. The range information is used
in finding near-empty and homogeneous regions. The average values are used in
the computation of coarse approximations to the integral across a segment.
This previous work is lacking in that the approximations that are introduced do
not include error bounds. Thus, although speed can be traded for accuracy, there
is no indication of how accurate the solution is at any given stage. Furthermore,
these techniques are limited to a strict front-to-back traversal of the dataset. The
method we develop here, which we call adaptive erwr bracketing addresses both
these concerns.
5.2.1 Overview
As in the method described by Danskin and Hanrahan, [DH921, we assume a
hierarchical organization of the data. In our implementation, we use a BSP tree
[FNN8Oj instead of an octree. The hierarchical structure is view-independent, and
can be built once for each dataset and stored as part of the data.
Our method is based on the method of piecewise constant bounds, thus we use
concepts and notation from Section 4.1.2. Associated with each segment k are two
intervals 1k and Tk, which respectively bracket the intensity and transparency of
124
the segment. We actually store amended intervals that are three-tuples consisting
of rn?n?rnum, estimated, and maximum contributions. We denote the estimated
value of an amended interval A by A. When doing arithmetic, the minimum and
maximum values are treated as ordinary intervals [Moo66], whereas the estimated
values are treated as scalars. For example, given two amended intervals A
[A-,A%A+]andB[B-,B,B+],wehaveA?B [A---HB+, A--HB, A+--HB-].
For each ray, the algorithm maintains a set 8 of contiguous segments f[dj?1, dj]?
along the ray. Initially, the set contains just a single segment whose endpoints give
the intersection of the ray with the entire volume. The algorithm proceeds itera-
tively, at each step selecting the segment from 8 that is estimated to contribute
the most error, then refining that segment either by splitting it in two, or else by
computing its contribution more precisely and, finally, updating the approxima-
tion to I(a, b) based on the new 8 and associated bounds. The process is repeated
until the approximation to I(a, b) is computed to within a given error tolerance,
?. This solution strategy is shown in pseudo code in Figure 5.17, where V is the
given volume dataset. In the following sections, we look at each these steps?select,
refine, and update in more detail.
5.2.2 Selection
A good selection strategy would be to choose the segment of 8 that contributes
the most error to the overall ray intensity 1. Unfortunately, this information is
not normally available. Note that finding the exact error contribution of each
segment is equivalent to the problem of computing the ray intensity with no error.
125
RenderVolume(Nffi ?)
[Render volume V to error tolerance ?
for each ray R do
let d0, d1 be the entry and exit points of R through V
initialize 8 ? fW0,d1]J.
initialize I using d1 --H d0 and global bounds on a and e.
while ?I >?do
1 [Select.] Choose a segment [dk?1, dk] from 8 to refine.
2. [Refine.] improve the bounds 1k and Tk for [dk-1, dk],
either by applying a quadrature rule or by
subdividing and bounding the smaller segments.
3. [Update.] Compute the new value of I
end while
output I for this ray.
end for
Figure 5.17: The main loop of the adaptive error bracketing algorithm.
126
Our selection strategy is based on an error estzrnate that provides a heuristic for
choosing the best segment to refine. It is important to note that the quality of the
error estimate will only effect the selection process. The algorithm is guaranteed
to maintain reliable bounds on I(a, b) at all times.
In principle, the error contributed by a single segment k to the overall ray inten-
sity can be described as the sum of two components. The first is error contributed
directly due to inaccuracy in the estimated intensity contribution of the segment,
1k The second is error contributed indirectly due to inaccuracy in the estimated
transparency of the segment, Tk. This latter error contribution is indirect because
the error is only introduced when the transparency estimate is multiplied by the
intensity of the combined segments behind it. The sum of these two components
may be written as
Ek = Tfront(k)( I%--HJ? + Tk--HTk Iback(k)),
where Tfront(k) is the accumulated transparency of the segments in front of k, and
Iback(k) is the intensity of the combined segments behind k. We do not know
either of these quantities exactly, but from Equation (3.4) we can estimate them
as follows:
Tfront(k) = k--H1			n			1
llTi,			fback(k) = Z			II Tj
i--H--Hi			?--Hk+1			j=k+1
Finally, we can bound the quantities 1% --H 1k and Tk --H Tk to arrive at a practical
estimate Ek for the error contribution of segment k:
Ek = Tfront(k) (11k I + ITkI Iback(k))
(5.1)
127
Computing and updating this error estimate efficiently for each segment k is
complicated by the fact that refining a single segment changes the error estimate
for all the others. We therefore employ a binary tree structure, which we call a
span tre?, to aid in propagating these updates efficiently. Each leaf of the tree
represents a single segment k of the ray. The intervals 1k and Tk are at the leaves.
Each internal node represents the union of all the segments below it, and stores an
intensity and transparency interval for this union. Note that the actual intensity
I(a, b) of the entire ray is contained in the I interval of the span tree's root.
Using this structure, error estimates for any segment k can be computed while
traversing the tree from the root to k. Initially, the quantities Tfront(k) and Iback(k)
are set to complete transparency and zero intensity, respectively. At each node, the
error formula from Equation 5.1 is evaluated for the front and back child segments.
The algorithm will descend either toward the front or toward the back of the ray,
depending on which segment contains the higher error estimate. If descending
towards the front, Iback(k) is updated to include the intensity of the back child. If
descending towards the back, Tfront(k) is updated to include the transparency of
the front child. This algorithm is illustrated in Figure 5.18 and shown in pseudo
code in Figure 5.19. The procedure selects a "promising" leaf segment in O(log n)
time. Note that the search described here does not guarantee that the leaf segment
with the highest estimated error will be chosen. However, a more careful selection
would require visiting all the leaf nodes. In practice, we have found that the more
careful 0(n) search is not generally worth the extra selection cost.
128
Figure 5.18: In order to avoid a sequential seard? for the segment with the highest
estimated error, a binary tree structure called a span tree is employed. Amended
intervals for intensity and transparency are stored at each node. Selection starts
by examining the root node. At each stage, processing is directed towards the
child node with the higher error, until a segment is located. This process is not
guaranteed to find the segment with the highest estimated error, however it runs
in only O(log n) time.
5.2.3 Refinement
The purpose of refinement is to obtain tighter bounds on the intensity and trans-
parency within the selected segment, thereby improving the approximation of the
overall intensity.
R?efinement may be performed in many different ways. We chose to employ a
two-tier strategy. To refine ray segments that span multiple voxels, we split the
segment in two. After splitting, we obtain new bounds on the functions of ab-
sorption and emission across each of the new segments. Since the new segments
are smaller, we expect the bounds to be tighter, thus resulting in a better approx-
imation. This strategy, which corresponds to the method of piecewise constant
bounds, could be applied below the voxel level at the cost of slow convergence. We
chose instead to introduce a second strategy that we apply when a segment spans
129
Select(root)
[Select a segment to refine, given the span tree rootj
return RecursiveSelect(root, 1, 0)
RecursiveSelect(s, Tfront, 1back)
if 5 is a leaf then
return s
else
[Descend tree via front child or back child]
f ? front child of 5
b ? back child of 5
1mid			? 1b + Tbiback
Tmid			? TfrontTf
[Evaluate Equation (5.1) for each child]
Ef ? Tfront( I `?			+ ITj			1mid)
Eb ? T?id( 1b + I Tb I 1back)
if Ef > Eb then
else
return RecursiveSelect(f, Tfront, 1mid)
return RecursiveSelect(b, Tmid, 1back)
end if
end if
Figure 5.19: The selection process for the adaptive error bracketing algorithm.
130
a single voxel. This strategy is to compute the integral along the segment directly
using one of the other controlled-precision quadrature rules from Chapter 4.
If refinement is done by subdivision, a splitting point must be chosen. In our
implementation this point is chosen according to the splitting planes of the BSP
tree. If refinement is done by application of a quadrature rule, a specific rule and
error tolerance must be chosen. We used the power series to compute the integral
to the highest-possible precision. This way, we know that we will not have to
revisit the node for further refinement.
5.2.4 Updating
Once we have refined a segment, we need to update the values of I and T in the
span tree. The update carries a dual purpose. First, it computes a new value
of I(a, b) that can be used to check if the current solution is within tolerance.
Second, it updates values that will be used in computing error estimates for the
next invocation of the selection code. The update is accomplished by a single
bottom-up traversal of the tree beginning with the refined span. Each time we
move to a parent node, we recompute I and T at that node by compositing the
intensity and transparency information from the children. The complexity of this
process is O(log n), where n is the number of spans. The algorithm for update is
illustrated in Figure 5.20 and given in pseudo code in Figure 5.21.
5.2.5 Initializing the BSP tree
Our implementation uses an axis-aligned BSP tree as a simple hierarchical or-
ganization of the voxel data. This data structure can be built once and stored
131
Figure 5.20: Because intermediate results are stored at the span tree nodes, up-
dates may be processed efficiently. During a single leaf-to-root traversal, revised
estimates and error bounds are propagated using pairwise compositing.
along with the data. Each node of the BSP tree stores amended intervals Q and
E that bracket the absorption and emission coefficients for all the voxels repre-
sented by that node. The estimated values, ? and 6% are computed by averaging
the respective coefficients. The bounds and estimates may be propagated during
a bottom-up traversal of the hierarchy [LH91]. If the original dataset has N vox-
els, the resulting BSP tree will have 2N --H 1 nodes. Each node requires a fixed
amount of storage, so the BSP tree requires only a constant factor more storage
than the original dataset. Unfortunately, this constant factor can become quite
large. Our current implementation stores 36 bytes per node, thus a dataset with
128 x 128 x 128 voxels requires 144 megabytes of storage. Although we have not
yet implemented this strategy, the BSP tree planes could also be used to subdi-
vide screen space. This information could be used to determine the parts of the
data that are required for the processing of an individual ray. This would reduce
resident memory requirements.
132
Update(s, Tnew, Inew)
[Update the bonnds associated with leaf node s?i
Ts <7Tnew
1s ? 1new
[Pwpagate the change up to the root.]
while s has a parent do
p ? parent of 8
f ? front child of p
b ? back child of p
[Compose T and I of p's children]
Tp ? TjT6
Ip ? Ij $ TfIb
s?p
end while
Figure 5.21: The update process for the adaptive error bracketing algorithrn.
133
Since we split the ray along BSP tree planes, a segment always corresponds
directly to a ray intersection with a single BSP tree node. Therefore, the intervals
Tk and 1k for a new segment k can be computed from the Q and E values stored
in the BSP tree, according to equations (4.10) and (4.11). Note that splitting
along BSP tree planes does not restrict ray direction. However, it does mean that
segments are not necessarily split directly in half during refinement.
We apply a quadrature rule when the ray segment to be refined passes though a
cell in the lowest level of the BSP tree. In this case, the BSP tree is queried for the
values of ? and ? at the corners of the cell. These values are used in determining
the coefficients of the emission and absorption polynomials, ?(t) and ?(t), along
the ray segment.
5.2.6 Results
The adaptive error bracketing scheme attempts to compute each ray to a speci-
fied accuracy using the smallest number of refinement operations. The cost of the
strategy comes in the selection and refinement operations that enable the determi-
nation of the most promising segment to refine. Selection and update each require
a traversal of the span tree, imposing an O(log n) cost at each stage, where n is the
number of segments under consideration. This computational cost is only worth-
while if it pays for itself in terms of a sufficient reduction in the overall number of
refinement operations.
To evaluate this tradeoff, we implemented another controlled-precision ray eval-
uation algorithm that implements both early termination and machinery to skip
134
Figure 5.22: Using the front-to-back algorithm, the frontmost part of the ray,
shown in black, is evaluated to the highest possible precision. Empty regions
are identified by the BSP tree and skipped, since they don't contribute to the
final result. While processing front-to-back, a coarse bound on the contribution
of the back part of the ray is maintained. Processing halts when the maximum
contribution of the back part of the ray falls below the user-specified error tolerance.
over large empty regions. In this manner, it is similar to the scheme presented by
Levoy [Lev9Oa], however our scheme allows refinement to a precise error bound.
The ray is refined front-to-back; thus, selection and update take constant time.
The frontmost region is subdivided until it is empty, or until it is the size of a
single voxel, in which case a quadrature rule is applied. Since we do not have the
luxury of revisiting a segment using a front-to-back scheme, the error tolerance is
set to machine precision. In order to maintain an error bound on I(a, b) during
each stage of the algorithm, we compute the maximum possible emission along the
back part of the ray that has not been evaluated. This is illustrated in Figure 5.22.
Termination occurs when the product of the accumulated transparency and this
maximum remaining emission falls below the error tolerance.
We tested the adaptive error bracketing scheme (AEB) described in this section
against the front-to-back method (FTB) outlined in the preceding paragraph using
three different datasets and a range of error tolerances. We measured performance
in terms of both total computation time and total number of quadrature rules
employed. We ran all our tests on an HP 750, a 70 MIPS machine.
The results are shown in Figures 5.24, 5.25, and 5.26. For each dataset, we
135
present a four-by-six array of images. The top three rows show the components
of the amended intensity interval for the image. From top to bottom, they are
the minimum, estimated, and maximum intensity respectively. The bottom row
shows the difference of maximum and minimum intensities for each pixel. The
six columns show the images produced for increasingly tight error bounds. These
bounds are listed below each column, both in terms of absolute error in ray intensity
(?) and relative error. Relative error is the ratio ?/I??, where 1w is the intensity
of the brightest pixel in the highest-quality image. We also show graphs of total
CPU time, and total number of quadratures below the images for each of the two
methods AEB and FTB.
The first example, shown in Figure 5.24, is a 64 x 64 x 64 dataset depicting
a uniformly emissive grid behind a cloud of absorptive material that decreases
linearly in density from left to right. The geometry of the dataset is shown in
Figure 5.23. The slow variation of cloud absorption allows for accurate bracketing
by the AEB algorithm without many applications of the quadrature rule. Since
the cloud is not completely opaque, the FTB algorithm cannot employ early ter-
mination, thus making it 125% to 2090% slower, as the relative error ranges from
0.2% to 170%.
The second example, shown in Figure 5.25, is representative of isosurface data.
It is a 192 x 112 x 73 CT dataset of a human pelvis. In this example, transparency
drops sharply upon intersection with the bone surface, reaching negligible values
after penetrating only a few voxels. In this domain, the FTB algorithm quickly
locates the first intersection with the bone surface and applies the quadrature rule
136
AbsorptiveCloud
II
Grid
Figure 5.23: The geometry of the dataset in Figure 5.24 (plan view). The grid,
indicated at right, is partially obscured by an absorptive cloud whose density
increases linearly as indicated with shading, The location ofthe view position is
shown at bottom left.
where it is most effective. The AEB algorithm, on the other hand, tries to find
regions where the bracketed error is acceptable, but ultimately applies quadrature
rules to many of the same voxels on the bone surface. On this high-gradient
dataset, the extra cost ofdeciding where to refine causes AEB to run between 85%
and 98% slower, as the relative error ranges from 0.5% to 30%.
The third example, shown in Figure 5.26, is a 97 x 97 x 116 electron density
map ofan SOD enzyme. Electron density is mapped directly to absorption ?, and
emission ? is obtained via gradient shading. Absorption is low enough that rays
penetrate deep into the volume, but high enough that few rays have any trans-
parency left when they exit the dataset. For low accuracy results (within 134%)
the FTB algorithm is slower than AEB. However, the benefit ofusing AEB grad-
ually decreases as accuracy increases, and for the highest-precision images (within
0.5%), AEB is 104% slower even though it uses only 77% ofthe quadratures.
The examples suggest that adaptive bracketing is very effective in reducing the
137
number of quadratures necessary to compute an accurate image. However, fewer
quadratures does not always translate into an overall speed-up, since the extra
complexity in choosing where to refine imposes additional cost. In the limit, when
the highest accuracy rendering is required, there is no advantage in trying to be
careful about where to refine, since any method would inevitably have to perform
all of the same quadratures along the entire ray. Thus, the adaptive bracketing
method pays off most when a lower-accuracy rendering is desired.
As demonstrated by the first example, the adaptive bracketing scheme has a
great advantage over the front-to-back approach for `?cloudy" or low-gradient data.
5.2.7 Discussion
There are several aspects of this research that warrant further research:
Other refinement strategies. The select--Hrefine--Hupdate scheme, which is central
to adaptive error bracketing, can make use of any tools that refine error estimates
within a segment. Currently, we have implemented only a simple two-tier scheme:
segments are split down to the voxel level, after which highly accurate quadrature is
applied. Possibilities for improvement include limiting the depth of the span tree,
thus applying quadrature over larger regions. The addition of other refinement
tools with intermediate costs and benefits might also be useful. For example, a
cheaper, lower-accuracy quadrature formula like the trapezoid rule might be one
such tool. There is a tradeoff since overall cost would be higher if the same segment
must be refined multiple times. We would also like to explore incorporating some
notion of the relative costs and benefits of the various refinement tools into the
138
selection strategy.
Appropriate error tolerances. The error bounds provided by adaptive bracket-
ing are absolute. This type of error bound is useful if results are desired in precise
units that are well-defined in advance. Unfortunately, this is rarely the case in
volume rendering. Instead, once an image is created it is typically scaled to make
effective use of the dynamic range of the display. This scaling also affects the ab-
solute error. We would therefore like to investigate an iterative strategy whereby
the absolute error tolerance is refined along with the image in order to achieve a
prescribed "relative" error bound.
Bounding other errors. There are many sources of error that are currently not
addressed but that may in fact be well-suited to the interval arithmetic approach
described here. For example, if bounds could be placed on the data acquisition
and reconstruction errors, this information could be maintained and propagated
quite naturally in our scheme. The roundoff error introduced due to floating-point
or fixed-precision arithmetic could also be accounted for in this way.
Irregular data. So long as we can bracket the contribution of any given sample to
a segment along a ray, we can apply our technique to achieve any given accuracy.
Furthermore, the BSP tree is a general-purpose data structure that does not require
rectilinear cells. Thus, irregularly-spaced data could be handled directly.
Exploiting Coherence. The current algorithm starts with a new span tree for
each ray. It should be possible to get a rough idea of where refinement will have to
occur by using information gathered while processing nearby rays. One approach
would be to use the structure of a neighboring span tree as an initial guess. This
139
type of approach was employed successfully for maintaining active edge lists for
polygon scan-conversion by Watkins [Wat7O].
Progressive refinement. Ideally, the effort invested in making a fast but low-
quality image could be put to good use when making a highly-accurate rendition
of the same image. The idea of progressive refinement is to generate approximate
images as intermediate results on the way to the final image [BFG86]. The pro-
gressive technique of Levoy [Lev9Ob] refines solutions by continually adding rays
and interpolating the regions between them. An approach based on adaptive error
bracketing would instead refine all rays simultaneously. However, such a scheme
would require storing a span tree for every ray, which would increase the overall
storage requirement to size O(rn? ? n?p2), where n? is the linear resolution of the
volume data, and p is the linear resolution of the display. This increased storage
requirement may be prohibitive for very high-resolution data or displays.
Interactive exploration. The span trees store intermediate results at each node,
thus enabling fast updates if the emission and absorption functions were modified.
If the span tree were stored for every ray, an interactive exploration tool could be
developed that would allow real-time manipulation of these data parameters.
Adding color. Color is extremely valuable in visualizing precisely those datasets
for which adaptive bracketing is most effective. Color can be added in a straight-
forward way to our algorithm without increasing storage requirements since the
BSP nodes need only bound either the greatest contribution in any one channel,
or else a single luminance contribution.
Loose bounds. As we discussed in Chapter 4, one problem with the current
140
controlled-precision techniques is that the actual error is often substantially smaller
than the computed bounds indicate. For example, in Figure 5.25, all the estimated
images appear identical, even though the pixel values of the leftmost image are not
guaranteed to be within 170% of the correct values. Since adaptive error bracketing
is based on the principle of doing the least work in refining a solution as possible,
its effectiveness is diminished by a failure to recognize the true quality of solution
at any given stage.
5.3 Summary and Future Directions
The algorithms presented in this chapter aim towards efficient solutions to the
problems of intuitive display and accurate imaging for volume rendering. We have
presented an algorithm for perspective projection that does not require random
access to the dataset and that avoids both undersampling and oversampling the
data. We also presented an adaptive error bracketing algorithm that attempts
to compute guaranteed-accuracy solutions with the fewest refinement operations
possible. It is more efficient than a conventional front-to-back algorithm for low-
precision solutions and for low-gradient data.
These techniques have been successful in that their computation and storage re-
quirements are easily met by today's workstation-level machines. However, neither
of the approaches is currently valid for rendering at interactive rates. The ability to
interactively explore the dataset is important for improving comprehension. Even
when computing still images, it is often helpful to generate fast preview images
to test rendering parameters. Currently, users seeking interactive feedback must
141
abandon precise error bounds and accurate perspective projections. An important
future research topic is to lift these restrictions.
142
2500
Figure 5.24: Adaptive error bracketing. Results using the grid dataset
1000
500
0
6:
2000
1500
1			1			1			1			1			1
4			16			64			256			1024 4096
Minimum
Quadratures (Thousands)
__ AEB
FTB
Estimated
Maximum
Maximum
Minus
Minimum
1			1			1			1			1
4			16			64			256			1024
% of Max.			170%			43%			11%			2.7%			.7%
CPU Seconds
250?			AEB
FTB
200
150
`00
50
0
6:
1			1			1			1			1			1
4			16			64			256			1024 4096
4096
.2%
143
Minimum
Maximum1
Maximum
Minus
Minimum
1			1			1			1			1
16			64			256			1024			4096
% of Max.			478%			119%			30%			7.5%			1.9%
CPU Seconds
150 AEB
FTB
125
loo
75
50
25
0
1			1			1			1			1			1
16			64			256			1024 4096 16384
16384
0.5%
Quadratures (Thousands)
400[AEB
FIB
300
200
loo
0
1			1			1			1			1			I
16			64			256 1024			4096 16384
Figure 5.25: Adaptive error bracketing. Results using the CT pelvis dataset
144
Minimum
Estimated ?
Maximum
Maximum
Minus
Minimum
1			1			1			1			1			1
16			64			256			1024			4096			16384
% of Max.			538%			134%			34%			8.5%			2.1%			0.5%
CPU Seconds
500
400
300
200
loo
0
1			1			1			1
16			64			256			1024 4096 16384
Quadratures (Thousands)
__ AEB
2000F			FTB
1500
1000
500
0
1			1			1			1			1			1
16			64			256 1024			4096 16384
Figure 5.26: Adaptive error bracketing. Results using the SOD enzyme dataset
Chapter 6
Imaging for Cardiology
In this chapter, we present the preliminary results of our ongoing project in ren-
dering cardiac ultrasound data. This case-study will illustrate the difficulties en-
countered in providing reliable imaging of the types of datasets that are generated
by modern scientific scanning equipment.
6.1 The goals of the project
Our continuing project is oriented towards the development of a low-cost, dynamic,
four-dimensional imaging system for the human heart in vivo. We are interested in
developing non-invasive diagnostic tools that are practical for use in a physician's
office, rather than at a special-purpose imaging center. The project is performed
in collaboration with researchers at the Cornell Medical Center, Hewlett-Packard,
The University of Firenze, and ESAOTE Biomedica. Our portion of the project
is limited to the display of the 3D data.
145
146
The overall medical goals of the project include the display of:
o+ muscle wall motion, useful for locating "dead" muscle tissue;
o+ valve action, for early detection of diseases such as Marfan's syndrome; and
o+ blood-flow patterns, in order to detect malfunctioning valves.
Without the capability for displaying dynamic, three-dimensional projections of
the heart, many of these goals are not realizable. For example, it is practically im-
possible to understand the twisting motion of an abnormal valve flap by examining
two-dimensional slices.
The project goals determined both the data acquisition method and the vi-
sualization platform that were used. Current clinical choices for gathering 3D
cardiac data include x-ray computed tomography (CT), magnetic resonance imag-
ing (MRI), and ultrasound. X-rays result in harmful ionizing radiation, thus CT is
considered an invasive technique. MRI systems are non-invasive and have recently
been overcoming early problems in providing the spatial and temporal resolution
necessary for diagnostic cardiology. However, these systems are currently too large
and too expensive for use as a casual diagnostic tool. We chose to base our work
on ultrasonic imaging. Ultrasound devices are small and inexpensive. Patients are
scanned using an external hand-held probe. The ultrasonic waves that are used in
normal diagnostic practice are not considered harmful.
For visualization, we required a mass-produced system with high-end graph-
ics capabilities. Although supercomputers are attractive for their computational
power and large RAM caches, we found them to be impractical for our purposes.
147
Personal computers are widely available, but currently lack both high-end graphics
and the resources necessary for handling large datasets. Today's graphics worksta-
tions represent a good tradeoff between the two extremes of supercomputers and
personal computers, and we decided to design our software for these machines.
Workstations are small enough and inexpensive enough to meet our clinical goals.
Since the power of personal computers continues to grow exponentially, we expect
that our software will be feasible for execution on personal computers in just a few
years, thus further reducing the cost of the total system.
6.2 Data Collection
We begin by reviewing the data-collection process. Ultrasonic imaging works in a
manner similar to that of radar. An ultrasonic transducer converts electrical energy
into ultrasonic energy. It produces an ultrasonic beam that passes through the
body. At each interface between structures with differing acoustical impedances,
some of the energy is reflected back towards the transducer in the form of an
echo. Upon receiving the echo, the transducer converts the acoustical signal into
electrical pulses. The time between emission of the beam and reception of the echo
is used to infer the depth of the interface. More detailed information concerning
the physics of ultrasound may be found in review articles, such as the one by
Havlice [11T79].
The process is extended to two dimensions by sequentially reorienting the ul-
trasound beams at regular angular intervals across a thirty- to ninety-degree arc.
The sweep can be performed either by mechanically rocking the transducer or
148
by electronically "beam-steering" the waveform generated by a phased array of
transducers. In either case, the system must pause to gather echoes before emit-
ting a beam at the next angular position, It is desirable that the entire sweep be
performed in a thirtieth of a second, in order to allow real-time video display of
the fan-shaped 2D slice. The speed of ultrasound in the body is approximately
1480 m/s. Assuming instantaneous reorientation of the transducer, up to 308 an-
gular samples are possible across a ninety-degree arc scanned to a depth of 16
cm. Currently available hardware attains 120 angular samples. This corresponds
to sampling the plane along scanlines separated by three-quarters of a degree in-
crements. The process of gathering this fan-shaped slice of data is illustrated in
Figure 6.1.
The simple extension of this scheme to three dimensions is not feasible since
the time required to sweep across even a small-sized volume is too large relative to
the speed of motion of the heart. Current technology does not allow us to gather
data from multiple beams simultaneously due to the interference that would result
from multiple ultrasonic waves. Thus we cannot reconstruct 3D static images by
taking a snapshot at a single point in time. Instead, we rely on the fact that
heart motion is cyclic. We reconstruct static images of the heart by combining
information obtained over several cardiac cycles. First, we take a planar snapshot
from a known position and orientation. Then we wait for the heart to return to the
same position and take another snapshot from a different orientation. This process
is repeated until enough planes of data are available. Then the planes are merged
to form a single 3D dataset. The timing of the snapshots is determined by electrical
149
(a)			(b)
(3
(c)			(d)
Figure 6.1: Gathering 2D ultrasound data. (a) Inferring depth from the time be-
tween emission and echo, an ultrasound transducer is able to gather data along
a line that extends into the body. (b) This line may be moved by mechanically
rocking the transducer (as shown) or by beam-steering a phased-array of transduc-
ers. (c) The transducer is swept over a 90 degree arc, gathering data in acquisition
lines from across a 2D plane. Data samples are indexed in polar coordinates. The
sweep occurs within 1/30th sec. (d) A 2D cardiac ultrasound image of a normal
volunteer.
150
signals from the autonomic nervous system that control the heart. These may be
measured outside the body by an electrocardiogram (EKG) device. We use the
EKG impulse as a reference frame for each heart cycle, and take our snapshots at
regular intervals after the pulse. The gating is performed electronically. Note that
this technique is not useful for diagnosing diseases that result in acyclic irregularity
of cardiac function.
The choice of planes that are used to gather three-dimensional data is lim-
ited by the fact that ultrasound waves are attenuated sharply by bone, shadowing
any structures at greater depths. Thus, the rib cage presents difficulties for the
acquisition of ultrasound data about the human heart. Our colleagues at the Uni-
versity of Firenze, ESAOTE Biomedica and New York Hospital developed an ultra-
sound probe with a rotating transducer for the acquisition of 3D data [PMM+90,
PCM+92]. The new probe can be held fixed while internally, the transducer rotates
about its central axis. By sweeping the fan-shaped scan profile across a full arc,
we are able to gather data within a cone-shaped region. The radius of the cone
is small enough in the near field so that the ultrasound probe may be manually
aligned with a gap between two ribs with no intersection of the scanning region
with the rib cage. A photograph of the ultrasound system during data acquisition
is given in Figure 6.2.
Since we are interested in studying the dynamic aspects of heart function, we
collect data at several points in the heart cycle. Because rotation of the transducer
is not instantaneous, it is more efficient to gather data from all time steps at a
single orientation before rotating the transducer. Data is collected during every
151
Figure 6.2: The ultrasound system during data collection. The hand-held probe
contains a rotating transducer that is controlled by the apparatus at center right.
The probe connects to a standard 2D ultrasound device at lower right.
152
0
(a)			(b)
Figure 6.3: Coordinate axes relative to the transducer for spherical and cylindrical
ultrasound data. (a) Spherical coordinate (r, 0, ?) axes. (b) Cylindrical coordinate
axes.
other heart cycle; the intervening cycles are used for rotating the transducer. In
order to obtain 62 radial slices for every time-step, 124 heart cycles are necessary,
requiring 75 to 123 seconds. After the scan, the imaging proceeds without further
participation by the patient.
With the addition of the transducer rotation, data is gathered at regularly-
sampled points in a spherical coordinate system. The full dataset is indexed by
(r, 0, ?, t), where r is distance along the acquisition line, 0 is the angle of the rocking
of the transducer and ? is the angle of rotation of the transducer. These axes are
shown in Figure 6.3a. The time-step within the heart cycle is given by t.
153
6.3 Current Limitations
Unfortunately, at the present time, it is not possible for us to obtain digital data
in the native spherical coordinate system. The bandwidth and fast RAM require-
ments for moving and storing the mass of data accumulated by the transducer over
short time intervals are too high for practical implementation. While we await de-
velopments in this area, we have been able to gather 2D slices of data after they
have been scan-converted for display by the ultrasound device. Conversion of the
data to video standard requires a resampling onto a 2D rectangular grid. This
step is performed by hardware inside the ultrasound device that employs a bilin-
ear interpolation algorithm [LL8O]. The scan-converted data is recorded on video
tape, which is later digitized for further processing. With scan-conversion of each
2D slice, the data arrives as regular samples in a cylindrical (r, y, ?) coordinate
system, with the axes defined as in Figure 6.3b. The resolution of the dataset is
as follows. Each 2D slice is made up of 512 x 512 8-bit pixels. Each 3D dataset
contains 62 slices. The heart cycle is sampled at 15 time steps. Thus, the total
size of the dataset is
512 x 512 x 62 x 15 243, 793, 920 bytes 232MB.
A great deal of information is lost using the current acquisition system. Con-
sider the many resampling steps that exist in the pipeline:
0 Analog ultrasound signals are converted to digital signals.
o+ The digital signals are resampled onto a rectangular grid.
154
o+ The gridded data is converted to analog form for 2D display.
o+ The 2D display is recorded on video tape.
o+ The tape is played back into a frame grabber, where the video frames are
captured and converted to digital data.
The sum of these steps causes a significant amount of degradation that cannot be
accounted for during subsequent processing.
6.4 Visualization of the data
Once the raw data has been digitized, the individual slices must be sorted according
to time index. For each time index, each slice is placed in its proper spatial
orientation, resulting in a single cylindrical (r, y, ?) dataset per time frame. The
process of collecting the 3D data is illustrated in Figure 6.4.
Before projection, it is necessary to convert the raw data to absorption and
emission functions. As described in Section 3.1.1, the first step is to identify the
objects of interest in the raw data. In this case the objects of interest are the
muscle walls and valve structures that make up the heart. The internal textures
of these tissues result in echoes throughout their mass, not just at the interfaces.
Since we are using parasternal views, the heart is the only internal organ present
in the scans. Classification is largely a matter of filtering out weak echoes, caused
by ultrasound interactions with blood cells, and strong "speckle" artifacts caused
by anisotropic reflection of ultrasound waves by structures such as the cartilagi-
nous fibers in the ventricle. We chose to use isosurface techniques to extract the
155
(a)
(b)
(c)
1?r11
Figure 6.4: Gathering 3D ultrasound data. The transducer gathers data at n?
angular positions. Illustration (a) shows the transducer and the fan-shaped slices
resulting from three transducer orientations. At each orientation, data is gathered
for n time intervals. After scan-conversion, the dataset is recorded onto videotape
as ni orientations, ?1, , ?m, with n rectilinear slices per orientation, as shown in
illustration (b). These slices are digitized and then sorted according to time index.
The result is n cylindrical 3D datasets,t1,... , tn, as shown in illustration (c).
156
midrange of echo intensities. The algorithm that we used was a direct implemen-
tation of Levoy's isosurface extraction algorithm [Lev88j that was described in
Section 3.1.1. Our results were acceptable as a first approach, even though noise
and speckle artifact were not completely eliminated. Once classification was com-
pleted, absorption parameters were assigned so that the muscle wall appeared as
a solid mass. Emission characteristics were determined by computing the gradient
of the data and applying the Phong shading model, as described in Section 3.1.3.
The output of the above process is a dataset with two scalar quantities at
each sample point: an absorption coefficient and an emission coefficient. We have
investigated three approaches to projecting this data that arrives as regular samples
across a cylindrical grid.
The first approach was to resample the data onto a regular rectilinear grid.
Resampling enables the use of the wide range of projection algorithms that are
currently available for this data topology, including the controlled-precision tech-
niques presented in this thesis. We chose to use a ray-casting algorithm and a
perspective camera model. The expensive conversion of sample points from cylin-
drical to rectilinear coordinates is performed only once. Problems arise using this
method because the original dataset samples Cartesian space nonuniformly. If the
rectilinear grid samples the core of the dataset adequately, it oversamples the outer
regions, imposing high storage costs. Conversely, if the outer regions are sampled
adequately the core is undersampled, potentially resulting in the omission of small
details. This problem is presented in Figure 6.5. An additional drawback to this
approach is that the original data is not sampled above the Nyquist frequency,
157
therefore resampling is guaranteed to introduce error. We are currently unable to
quantify the magnitude of this error.
A second approach was based on the splatting technique that we described in
Section 3.3.4. We developed a new splatting algorithm that operates in cylindrical
coordinates. Splatting is difficult to apply to cylindrical data since the shape of the
filter kernels must vary over space to account to the varying sampling rate. Ideally
the spherical filter kernels that are used for rectilinear data need to be stretched
and bent along the curvature of the cylindrical grid. We chose to approximate
this transformation by merely stretching the kernels along the axis perpendicu-
lar to each radial line. This results in kernels with the same ellipsoid shape that
has been studied in the splatting literature [Wes9O,LH91]. An illustration com-
paring the coverage of a 2D polar-coordinate plane using uniform circular filter
kernels and elliptical filter kernels is illustrated in Figure 6.6. Once the proper
orientation and stretching of a kernel is determined, it may be projected using the
hierarchical splatting technique and a hardware z-buffer [LH91]. This technique
is an improvement over the resampling method in that it operates in the native
coordinate system of the data. However, it spends most of its computational en-
ergy computing the splats for the multitude of sample points at the core. This
is unfortunate since the screen-space contribution of these sample points is mini-
mal. The largest drawback of the approach is the cascade of approximations that
are employed. In raw form, splatting introduces error due to the improper treat-
ment of filter kernel overlap. Hierarchical splatting approximates the raw splats by
Gouraud-shaded polygons and approximates image compositing operations with
158
..... .......
.o+			..			...
(a)
(b)
? (b)
(a)
Figure 6.5: Resampling a cylindrical grid onto a rectilinear grid. In this example
the cylindrical grid has 128 radial lines, each containing 256 samples. The spacing
between rectilinear grid lines was chosen to be equal to the spacing between samples
along a single radial line. Two portions of a horizontal planar slice of the dataset
are expanded. (a) Note that near the core, many cylindrical samples fall within
each rectilinear cell. (b) A short distance away, most cells contain no cylindrical
samples.
159
fixed-precision arithmetic. Our modification of the technique adds further inaccu-
racy due to the approximation of the ideal cylindrical splat shape.
Our third approach was to use a ray-casting algorithm that was adapted to
operate in cylindrical coordinates. We took evenly-spaced sample points along
the ray and computed the volume rendering integral using the piecewise-constant
approximation. Thus the algorithm was similar in spirit to the one presented
by Levoy [Lev88]. Each sample point along the ray was converted from (r, y, z)
coordinates to (r, y, ?). The sample points at the corners of the enclosing cell
were located, and trilinear interpolation along each of the cylindrical-coordinate
axes was used. The drawback of this approach is that choosing an acceptable
uniform sampling rate along the ray is difficult. The situation is analogous to that
encountered when resampling the data onto a regular grid, as shown in Figure 6.7.
In this instance, because of concerns about computational expense, we chose to
undersample the core. While it is desirable to spend less effort computing the
contribution of the core, this method ignores large numbers of samples entirely,
thus introducing the possibility of gross inaccuracy.
6.5 Discussion
Our studies with these three algorithms revealed the many difficulties in apply-
ing current-technology volume rendering algorithms solving to real-world scientific
visualization problems. The techniques produced images that, when viewed indi-
vidually, appeared plausible. However, they differed markedly from one another.
Images produced by each of the methods with the same dataset is given in Fig-
160
(a)
(b)
Figure 6.6: Comparison of uniform circular splats to nonuinform elliptical splats
in polar coordinates. (a) Uniform splats do not account for the varying spatial
sampling rate. Filters close to the coordinate system origin overlap while gaps
appear between filters in more distant regions. (b) The circular splats are stretched
along the axis perpendicular to the radial lines. A more uniform coverage of the
plane is attained even though ideally, the splats would be bent to conform to
curvature of the radial axis as well. For the purpose of illustration, we chose to
stretch the filters until they "just touch" one another. In practice, a uniform
fraction of overlap is maintained.
161
Figure 6.7: Ray-casting onto a cylindrical-coordinate grid. When evenly-spaced
sample points are taken, some voxels at the core may not be sampled at all, as
shown in this diagram. The cost of using a regular sampling rate that ensures that
all voxels in the ray's path are considered is prohibitively high.
162
ure 6.8. It was the disparity that is evident in side-by-side comparisons that led
us to begin an error analysis of the volume rendering process.
There are many sources of error in the process used in our cardiac imaging
system. They include:
Acquisition Error. The data that arrives for rendering contains inaccuracies,
including limited spatial and temporal resolution, noise, inaccurate positioning of
the transducer, and lossy scan-conversion.
Reconstruction Error. The volume rendering algorithms operate on 3D scalar
fields. Thus, a reconstruction filter must be applied to the discrete data samples
that are produced during acquisition. Due to practical considerations, ideal fil-
ters cannot be used, thus error is introduced during this phase. The process is
complicated by the fact that data does not arrive as samples on a rectilinear grid.
Classification error. The raw data does not correspond to the parameters of
emission and absorption that are required by the rendering model. The raw data
must be interpreted to produce these values. Incorrect classification will result in
misleading images.
Rendering Error. Once the functions of emission and absorption are determined,
the volume rendering model must be evaluated. Since analytic expressions for the
resulting integrals are not known, approximation error is introduced.
The controlled-precision methods that were introduced in this thesis address
the problems of rendering error. We presented techniques for computing individual
rays to any specified precision, given polynomial or piecewise-polynomial functions
of absorption and emission. We also developed an algorithm that, in principle, will
163
(a)			(b)
(c)
Figure 6.8: Comparison of three methods for rendering cylindrical-coordinate data:
(a) Ray casting with the original grid; (b) Ray casting after resampling onto a
regular grid; (c) Cylindrical splatting. The variation in appearance is due to the
fact that each method employs a different set of approximations.
164
compute accurate reconstructions of screen space.
This work needs to be extended so that the controlled-precision algorithms can
operate directly on cylindrical and spherical coordinate data. This will eliminate
a source of resampling error. In order to apply the techniques, it will be necessary
to derive analytic expressions for the interpolated functions of emission and ab-
sorption as they pass through the cells of a curvilinear-coordinate grid. While the
evaluation of the resulting expressions is likely to be expensive, the images that are
produced can be presented with known accuracy. Using a cylindrical-coordinate
version of the adaptive error bracketing algorithm might help to alleviate some of
the problems of computational expense. We expect that such an algorithm would
be able to use coarse approximations within the core while maintaining reliable
error bounds for the final values. Thus, an efficiency that was not attainable with
any of our existing rendering techniques would become available.
As the above listing implies, there are many more sources of error that need to
be investigated. The logical next step would be to quantify the error due to non-
ideal reconstruction. This analysis should be performed so that it applies to data
that is not regularly sampled in any coordinate system. This will be useful, for
example, since we know that the hand-held transducer does not remain perfectly
fixed in space during a two-minute scanning session. Variation in the position
and orientation of the transducer will result in data that arrives on arbitrarily
positioned planes. Precise position and orientation information will be necessary
before we can begin to consider this effect in our analysis.
Although more research is needed in all aspects of our cardiac imaging project,
165
we are encouraged by our results thus far. Our basic goal of a low-cost, non-
invasive diagnostic imaging tool still appears to be attainable. Although many of
the difficulties we face are unique to our project, we find that similar obstacles
arise when imaging any real-world volume data. Our results thus far will help
scientists by providing guarantees of accuracy of the computed projections. Our
results will also provide a starting point for further error analysis of the volume
rendering process.
Chapter 7
Conclusion and Future Directions
Three major contributions to the field of volume rendering have been presented:
Controlled-precision volume projection. We developed techniques for com-
puting the intensity accumulated along individual viewing rays to any specified
precision. We presented four approaches for solving the emission-absorption inte-
gral across individual ray-voxel intersections; piecewise-constant bounds, trapezoid
rule, Simpson's rule, and a power series method that is practical for computing
high-precision solutions. These results were extended to compute the emission-
absorption integral to any specified precision for rays that pass through multiple
voxels. Our approaches represent a departure from existing volume rendering algo-
rithms in that they allow quantitative assessment and user control over the quality
of solutions to the line integral that are produced.
Adaptive error bracketing. Adaptive error bracketing addresses the problem
of efficiently computing point-samples of screen-space to any specified precision.
166
167
Using this method, computation is directed toward large sources of error, and
easy-to-compute approximations are used whenever possible. The approach allows
continual refinement of the solution until a user-specified error tolerance is met,
thus avoiding overcomputation while maintaining guaranteed error bounds. The
algorithm allows processing of the data without imposing a strict front-to-back or
back-to-front evaluation order. Our studies have shown that adaptive error brack-
eting saves significant computational effort when working with smooth datasets or
when computing low-accuracy images.
Efficient perspective projection. We also presented efficient volume rendering
algorithms for computing perspective projections. An optimal paging strategy was
introduced that makes ray-casting practical in situations where the dataset is too
large to fit into RAM memory. A ray-splitting approach to adaptive supersampling
was also presented. This technique ensures that all data features are accounted for
in the final image, while at the same time avoiding overcomputation in regions of
the dataset that are near the eyepoint.
In the introduction to this thesis, we outlined several shortcomings of the state-
of-the art in volume rendering. While the above contributions represent an im-
provement in terms of both accuracy and efficiency, there remains much work to
be done:
Accuracy. This thesis provides tools for controlling one source of error in the
volume rendering process: approximating the value of the illumination function
at point-samples. An important next step is to extend these results to include
controlling projection error across entire areas of screen-space. In addition, there
168
are many other sources of error in volume rendered images that need to be ad-
dressed. These include non-ideal reconstruction of functions in object-space and
in screen-space, and the propagation of errors that are introduced during data
acquisition.
Efficiency. The techniques presented in this thesis can reduce computational and
storage requirements so that volume-rendered images may be produced in several
minutes on workstation-level machines. An important next step is to provide
methods providing scientists with interactive feedback. Since the adaptive error
bracketing algorithm allows progressive refinement along individual rays, we believe
it will serve as a good starting point for this type of research.
Arbitrary data formats. The main research contributions of this thesis all as-
sume a restricted data format; only scalar data organized on a regular rectilinear
grid was considered. While the theoretical work in controlled-precision projection
remains valid for any scalar field, there remain many unresolved implementation
issues. These include proper reconstruction of non-rectilinear data and the inter-
section of a ray with irregularly-shaped voxels. Generalization of the algorithms
to the visualization of multiple scalar and vector fields await further research in
techniques for the effective display of such information.
Clearly, significant research is necessary to transform volume rendering into a
reliable general-purpose visualization tool. We hope that the contributions of this
thesis represent a step in the direction of solving this larger problem.
Appendix A
Nomenclature
This chapter is a summary of the nomenclature that is defined in the thesis. It is
included only for reference.
Symbol Definition
Absoprtion
Emission
I(a, b) Intensity along a ray between points a and
b
tc(t)			Integrating Factor
p			Density
Page Introduced
7
8
7
8
11
169
170
Symbol Definition
V(v) Optical distance from a point on a ray to
the light source
?(u) Reflectance at a point on a ray
Raw data value
T(a, b) Transparency along a ray between points
a and &
f(t) The volume rendering integrand
Pn(t) Power series approximant
Rn (t) Power series remainder
Page Introduced
14
14
24
29
62
75
75
Interval Arithmetic
As used in this thesis, intervals are contigous regions of the real axis The interval
bounded by scalars a and b is denoted by the bracketed pair [a, bi. Variables
representing intervals are denoted by variables in bold letters, ?. e. Q, E. Additional
notation regarding intervals is given in Table A.1.
171
Table A.1: Notation for interval arithmetic. The definitions below apply to the
interval x			[a, b]
Notation Meaning			Definition
Lower bound on x			a
Upper bound on x			b
x			Midpoint of x			-21(a +
I x			Width of x			b --H a
Appe?dix B
Power Series Co?vergence Proof
In this section we prove that lim??? Rn(t) 0 for all 1. Equation (4.31) expresses
Rn(t) as the sum of a finite number of terms Tn?1 where
Tn?i(t)
Since j is fixed, these terms are of the form
Kit?An?j
where K1 is the constant ?q(i)/j!. It will be sufficient to show that
lim Tn 0.
NVe first show that tT1(t),T2(t),.. .J is bounded for all t, from which conver-
gence to zero will follow easily.
Substituting for A from Equation (4.33) yields
K N p?(k)Ah-k
n?JkZ (kYl)!
172
173
for n --H 1 > N+j. Letting
we have
Now letting
we have
N ______
K2 K1 ??,,aiX (k --H 1)!'
N
? K2 ? An?j?k
--H J k=1
= K2  )?Wkk
--H J k--H1
K2 N
--H n --H ZTn?ktk.
k--H--H1
K3 A% mNax Itk
k_
K3 N
Tn-i ? --H ?Tn-k
k--H--H1
As n ? 00, Tn-i (t) cannot grow without bound, since when
K2			1
fl--Hi
Tn--Hi (t) is strictly less than the average of its preceding N values, and hence strictly
less than the maximum of these values. Thus after a finite number of terms Tn(t)
is a decreasing sequence. Therefore we may bound the sequence
Tn < maxtTi(t), T2(t), . . .? = B ? 00.
Finally, from Equation (B.4) it follows that
T i<K3NB			K4
fi
174
and
lim K4
_ --0.
n?oo fl
This completes the theorem. E
Note that Tn, and hence 1?n is bounded over all j and all t E [a, b]. This
bound allows us to interchange limiting operations and integration,
$6nl? R? J% $6m
which we use implicitly in Chapter 4.
BIBLIOGRAPHY
[Act7Ol Forman 5. Acton. Numerneal Methods That Work. Harper & Row,
New York, 1970.
[ASK92] Ricardo 5. Avila, Lisa M. Sobierajski, and Arie E. Kaufman. Towards
a comprehensive volume visualization system. In Proceedings: Visual-
ization `92, pages 13--H20. IEEE Computer Society, 1992.
[BFG86] Larry Bergman, Henry Fuchs, and Eric Grant. Image rendering by
adaptive refinement. Computer Graphics, 20(4):29--H38, August 1986.
[B1i82] James F. Blinn. Light reflection functions for simulation of clouds and
dusty surfaces. Computer Graphics, 16(3):21--H29, July 1982.
[BT75] Phong Bui-Tong. Illumination for computer-generated pictures. Corn-
munications of the ACM, pages 311--H317, June 1975.
[Cha60] 5. Chandrasekar. Radiative Transfer. Dover Publications, New York,
1960.
[CPC84i Robert L. Cook, Thomas Porter, and Loren Carpenter. Distributed
ray tracing. Computer Graphics, 18(3):137--H145, July 1984.
[CT82] Robert L. Cook and Kenneth E. Torrance. A reflectance model for
computer graphics. ACM Transactions of Graphics, 1(1):7--H24, 1982.
[DCH88] Robert A. Drebin, Loren Carpenter, and Pat Hanrahan. Volume ren-
dering. Computer Graphics, 22(4):65--H74, August 1988.
John Danskin and Pat Hanrahan. Fast algorithms for volume ray
tracing. In 1992 Boston Workshop on Volume Visualization, pages
91-98. ACM SIGGRAPH, 1992.
[DH92]
175
176
[FKN80]
[FvFH90j
H. Fuchs, Z. M. Kedem, and B. F. Nayloy. On visible surface gener-
ation by a priori tree structures. Computer Graphics, 14(3):124--H133,
July 1980.
James D. Foley, Andries van Dam, Steven K. Feiner, and John F.
Hughes. Computer Graphics: Principles and Practice, Second Edi-
tion. Systems Programming Series. Addison-Wesley, Reading, Mas-
sachusetts, 1990.
[Gla89] Andrew 5. Glassner, editor. An Introduction to Ra? Tracing. Aca-
demic Press, San Diego, California, 1989.
[Han90] Pat Hanrahan. Three-pass affine transforms for volume rendering.
Computer Graphics, 24(5):71--H78, November 1990.
[HB86]
[HBP+ 891
[HBP+901
[HT79]
[HT5G911
Karl Heinz Ho?hne and Ralph Bernstein. Shading 3D-images from CT
using gray-level gradients. IEEE Transactions on Medical Imaging,
MI-5(1):45--H47, March 1986.
Karl Heinz Ho?hne, Michael Bomans, Andreas Pommert, Martin
Riemer, Carsten Schiers, Ulf Tiede, and Gunnar Wiebecke. 3D vi-
sualization of tomographic data using the generalized voxel model. In
Proceedings of the 1989 Chapel Hill Workshop on Volume Visualiza-
tion, pages 51--H57, May 1989.
Karl Heinz H5hne, Michael Bomans, Andreas Pommert, Martin
Riemer, Carsten Schiers, Ulf Tiede, and Gunnar Wiebecke. 3D visu-
alization of tomographic data using the generalized voxel model. The
Visual Computer, 6(1): 28--H36, 1990.
James F. Havlice and Jon C. Taenzer. Medical ultrasonic imaging: An
overview of principles and instrumentation. Proceedings of the IEEE,
67(4):620--H641, April 1979.
Xiao Dong He, Kenneth E. Torrance, Francois X. Sillion, and Don-
ald P. Greenberg. A comprehensive physical model for light reflection.
Computer Graphics, 25(4), August 1991.
[KK89] James T. Kajiya and Timothy L. Kay. Rendering fur with three di-
mensional textures. Computer Graphics, 23(3):271--H280, July 1989.
[Kru90] Wolfgang Krueger. Volume rendering and data feature enhancement.
Computer Graphics, 24(5):21--H26, November 1990.
177
[kS88] Avinash C. Nak and Malcom Slaney. Principles of Computerized To-
mographic Imaging. IEEE Press, New York, 1988.
[NV84] James T. Kajiya and Brian P. Von Herzen. Ray tracing volume densi-
ties. Computer Graphics, 18(3):165--H173, July 1984.
[Lev88] Marc Levoy. Display of surfaces from volume data. IEEE Computer
Graphics and Applications, 8(3):29--H37, May 1988.
[Lev89]
Marc Levoy. Display of Surfaces from Volume Data. Ph.D. disser-
tation, The University of North Carolina at Chaptel Hill, 1989. Also
published as Department of Computer Science Technical Report TR89-
022.
[Lev9Oa] Marc Levoy. Efficient ray tracing of volume data. Transactions on
Graphics, 9(3):245--H261, July 1990.
[Lev90b] Marc Levoy. Volume rendering by adaptive refinement. The Visual
Computer, 6(1):2--H7, February 1990.
[Lev92]
[LH91]
[LL8oj
[LTG92]
Marc Levoy. Volume rendering using the Fourier projection-slice theo-
rem. Technical Report CSL-TR-92-521, Stanford University Computer
Systems Laboratory, April 1992.
David Laur and Pat Hanrahan. Hierarchical splatting: A progres-
sive refinement algorithm for volume rendering. Computer Graphics,
25(4):285?288, July 1991.
H. G. Larsen and 5. C. Leavitt. An image display algorithm for use
in real-time sector scanners with digital scan converters. In B. R.
McAvoy, editor, 1980 Ultrasonics Symposium Proceedings, pages 763
765. IEEE Group on Sonics and Ultrasonics, November 1980.
Dani Lischinski, Fillipo Tampieri, and Donald P. Greenberg. Discon-
tinuity meshing for accurate radiosity. IEEE Computer Graphics and
Applications, 12(6):25--H39, November 1992.
[Luc92] Bruce Lucas. A scientific visualization renderer. In Proceedings: Vi-
sualization `92, pages 227--H233. IEEE Computer Society, 1992.
[Mal9O] Tom Malzbender. Fourier volume rendering. In Proceedings of the
1990 Hewlett-Packard Graphics Symposium, 1990.
178
[Mal93l Tom Malzbender. Fourier volume rendering. A CM Trunsactions on
Cruphics, 12(3), July 1993.
[MHC9oj
Nelson Max, Pat Hanrahan, and Roger Crawfis. Area and volume
coherence for efficient visualization of 3d scalar functions. Computer
Crnph?cs, 24(5):27--H33, November 1990.
[Moo66] Ramon E. Moore. Interval Analysis. Prentice-Hall, 1966,
[NA92] Kevin L. Novins and James R. Arvo. Controlled precision volume
integration. In 1992 Boston Workshop on Volume Visualization, pages
83-89. ACM SIGGRAPH, 1992.
[NAS92]
[NH92]
[NSG90]
[0575]
[PC82]
[PCM+92]
[PD84]
Kevin L. Novins, James R. Arvo, and David H. Salesin. Adaptive
error bracketing for controlled-precision volume rendering. Technical
Report TR92-1312, Cornell University Department of Computer Sci-
ence, Ithaca, New York, November 1992.
Paul Ning and Lambertus Hesselink. Vector quantization for volume
rendering. In 1992 Boston Workshop on Volume Visual&ation, pages
69-74. ACM SIGGRAPH, 1992.
Kevin L. Novins, Irancois X. Sillion, and Donald P. Greenberg. An effi-
cient method for volume rendering using perspective projection. Com-
puter Graphics, 24(5):95--H102, November 1990.
Alan V. Oppenheim and Ronald W. Schafer. Digital Signal Processing.
Prentice-Hall, Englewood Cliffs, New Jersey, 1975.
M. Potmesil and I. Chakravarty. Synthetic image generation with a
lens and aperture camera model. A CM Transactions on Graphics,
1(2):85--H108, April 1982.
R. Pini, M. Costi, G.A. Mensah, L. Masotti, K. L. Novins, D. P. Green-
berg, B. Greppi, M. Cerofolini, and R.B. Devereux. Computed tomog-
raphy of the heart by ultrasound. In Computers in Carliology 1992,
pages 17--H20. IEEE Computer Society, 1992.
Thomas Porter and Tom Duff. Compositing digital images. Computer
Graphics, 18(3):253 259, July 1984.
William II. Press, Brian Flannery, Saul A. Teukolsky, and
William T. Vetterling. Numerical Recipies in C: The Art of Scien-
tific Computing. Cambridge University Press, Cambridge, 1988.
[PFTV88]
179
[PMM+901
Riccardo Pini, Elisabetta Monnini, Leonardo Masoti, Kevin L.
Novins, Donald P. Greenberg, Barbara Greppi, Marino Cerofolini, and
Richard B. Devereux. Echocardiographic three-dimensional visual-
ization of the heart. In Karl Heinz H5hne et al., editor, SD Imag-
ing in Medicine: Algorithms, Systems, Applications, pages 197?215.
Springer-Verlag, Berlin, 1990. NATO ASI Series, Volume F 60.
[Rat72] Floyd Ratliff. Contour and contrast. Scientific American, 226(6):91--H
101, June 1972.
[RB74] Earl D. Rainville and Phillip E. Bedient. Elementary Differential
Equations. Macmillan, New York, 1974. Fifth Edition.
[RT87]
[Rus88]
Holly E. Rushmeier and Kenneth E. Torrance. The zonal method for
calculating light intensities in the presence of a participating mediu'?.
Computer Graphics, 21(4):293--H302, July 1987.
Holly Edith Rushmeier. Realistic Image Synthesis for Scenes with Ra-
diatively Participating Media. Ph.D. dissertation, Cornell University,
Ithaca, New York, May 1988.
[Sab88] Paolo Sabella. A rendering algorithm for visualizing 3D scalar fields.
Computer Graphics, 22(4):51--H57, August 1988.
[SB80]			Introduction to Numerical Analysis.
1980. Translated by R. Bartels, W.
[SB87]
[SN89]
J. Stoer and R. Bulirsch.
Springer-Verlag, New York,
Gautschi, and C. Witzgall.
John M. Snyder and Alan H. Barr. Ray tracing complex models con-
taining surface tessellations. Computer Graphics, 21(4):119--H128, July
1987.
Peter Shirley and Henry Neeman. Volume visualization at the cen-
ter for supercomputing research and development. In Proceedings of
the 1989 Chapel Hill Workshop on Volume Visualization, pages 17--H20,
May 1989.
George W. Sherouse, Kevin L. Novins, and Edward L. Chaney. Com-
putation of digitally reconstructed radiographs for use in radiotherapy
treatment design. International Journal of Radiation Oncology, Biol-
ogy and Physics, 18(3):651--H658, March 1990.
[SNC90]
180
?Sob63]			V.
[SSW86l
[ST9O]
[TF79]
[TL93]
[UK88]
wPR?+?92i
[Wat70]
V. Sobolev. A Treatise on Radiative Trensfer. D. Van Nos-
trand Company, Princeton, New Jersey, 1963. Translated by 5.1.
Gaposchkin.
Daniel 5. Schlusselberg, Wade K. Smith, and Donald J. Woodward.
Three-dimensional display of medical image volumes. In Preceedings of
NCCA `86, pages 114--H123. National Computer Graphics Association,
March 1986.
Peter Shirley and Allan Tuchman. A polygonal approximation to
direct scalar volume rendering. Computer Craphics, 24(5):63--H70,
November 1990.
George B. Thomas, Jr. and Ross L. Finney. Calculus and Analytic
Geometry. Addison-Wesley, New York, fifth edition, 1979.
Takashi Totsuka and Marc Levoy. Frequency domain volume render-
ing. Computer Graphics, pages 271--H278, August 1993.
Craig Upson and Michael Neeler. V-buffer: Visible volume rendering.
Computer Graphics, 22(4):59--H64, August 1988.
Dirk Vandermeulen, Peter Plets, Steven Ramakers, Paul Suetens, and
Guy Marchal. Integrated visualization of brain anatomy and cerebral
blood vessels. In 1992 Boston Workshop on Volume Visualization,
pages 39--H46. ACM SIGGRAPH, 1992.
G. 5. Watkins. A real time visible surface algorithm. Ph.D. The-
sis UTEC-CSc-70-101, Computer Science Department, University of
Utah, June 1970. NTIS AD-762 004.
[Web90] Robert E. Webber. Ray tracing voxel data via biquadratic local surface
interpolation. The Visual Computem 6(1):8--H15, 1990.
[Wes89l Lee Westover. Interactive volume rendering. In Proceedings of the
Chapel Hill Workshop on Volume Visualization, May 1989.
[Wes9O] Lee Westover. Footprint evaluation for volume rendering. Computer
Graphics, 24(4):367 376, August 1990.
[WG91] Jane Wilhelms and Allen Van Gelder. A coherent projection approach
for direct volume rendering. Computer Graphics, 25(4):275--H284, July
1991.
181
[Whi8O] Turner Whitted. An improved illumination model for shaded display.
Commun?cat?ons of the ACM, 23(6):343--H349, June 1980.
[WilS3] L. Williams. Pyramidal parametrics. Compute? Graphics, 17(3): 1--H11,
July 1983.
Peter L. Williams and Nelson Max. A volume density optical model.
In 1992 Boston Workshop on Volume Visualization, pages 61 68. ACM
SIGGRAPH, 1992.
[WM92]
