.po 1.125i
.ll 6.5i
.ls 1
.nr tm 6v
.nr bm 7v
.nr pp 9
.nr sp 10
.EQ
delim $$
.EN
.ps 12
.b
.(b C
Volume Rendering using the
Fourier Projection-Slice Theorem
.)b
.r
.ps 10
.(b C
.i "Marc Levoy"
.)b
.ps 9
.(b
.(c
Computer Science Department
Center for Integrated Systems
Stanford University
Stanford, CA 94305-4070
Email: levoy@cs.stanford.edu
.)c
.)b
.ps 10
.2c 0.25i
.b
.(b C
Abstract
.)b
.r
.ps 9
The Fourier projection-slice theorem states that the inverse transform of a
slice extracted from the frequency domain representation of a volume yields a
projection of the volume in a direction perpendicular to the slice. This
theorem allows the generation of attenuation-only renderings of volume data in
$O(N sup 2~ log~ N)$ time for a volume of size $N sup 3$. In this paper, we
show how more realistic renderings can be generated using a class of shading
models whose terms are Fourier projections. Models are derived for rendering
depth cueing by linear attenuation of variable energy emitters and for
rendering directional shading by Lambertian reflection with hemispherical
illumination. While the resulting images do not exhibit the occlusion that is
characteristic of conventional volume rendering, they provide sufficient depth
and shape cues to give a strong illusion that occlusion exists.
.ps 10
.b "Keywords:"
.ps 9
Volume rendering, Fourier projections, Shading models,
Scientific visualization, Medical imaging
.sh 1 "Introduction"
.pp
Volume rendering is a technique for visualizing a sampled scalar or vector
volume of three spatial dimensions by modeling the propagation of light in a
participating medium. In its most general form, parameters of the data are
mapped to physical parameters of the medium such as energy density, absorption
and scattering coefficients, and scattering phase function [Krueger90]. The
integro-differential equation governing light transfer is then solved to obtain
an energy equilibrium for the volume [Siegal81].
.pp
To reduce solution costs, practical volume rendering algorithms assume a low
albedo, allowing them to ignore the effects of interreflection [Levoy88,
Sabella88, Drebin88]. The resulting equation can be evaluated independently
for each viewing ray using an expression of the form
.EQ (1.1)
I~~ =~~
size +3 int from {"-"inf} to inf~
f~(t)~
exp size +4 [ "-"
int from {"-" inf} to t~ g (u)~ du~ size +4 ]~
dt
.EN
where $t$ represents distance along the ray. This expression is typically
evaluated by digital compositing [Porter84], which in the context of volume
rendering is equivalent to numerical integration using the rectangle rule.
Given a cubical volume measuring $N$ voxels on a side, and assuming pixel
spacing equal to voxel spacing, the cost of generating an image using these
algorithms is $O(N sup 3 )$.
.pp
If we simplify the physical model still further by ignoring emission and first
scatter, we produce an image that looks like an X-ray. The expression for each
viewing ray now takes the form
.EQ (1.2)
I~~ =~~
exp size +4 [ "-"
int from {"-" inf} to inf~ g (t)~ dt~ size +4 ] .
.EN
.pp
For an imaging geometry that employs uniformly spaced parallel rays,
expressions of this form can be evaluated efficiently using the Fourier
projection-slice theorem. This theorem allows us to compute integrals over
volumes by extracting slices from a frequency domain representation of the
volume. The application of this theorem to image synthesis has been
independently proposed in [Dunne90] and [Malzbender]. An application to MR
angiography is described in [Napel91].
.pp
The computational advantage of the projection-slice theorem is, at least in
theory, considerable. Given a cubical volume and pixel spacing equal to voxel
spacing as above, the cost per image is $O(N sup 2~ log~ N)$ after an $O(N sup
3~ log~ N)$ preprocessing step to compute the 3D frequency domain
representation. For a volume and image of width $N~ =~ 256$, volume rendering
using digital compositing requires $16M$ operations, while rendering using the
projection-slice theorem requires only $512K$ operations, a speedup of more
than an order of magnitude.
.pp
There are two disadvantages to this approach. From a practical standpoint,
signal processing considerations have thus far prevented researchers from
obtaining speedups matching the theoretical complexity. We expect these
problems to be solvable, and, aside from a brief explanation in section 2, we
do not address them further in this paper. A more fundamental disadvantage of
the Fourier approach is that the resulting images do not exhibit occlusion and
are consequently not as useful as volume renderings for many applications.
Fortunately, occlusion is only one of many cues employed by the human visual
system to determine the shape and spatial relationships of objects. Other cues
include perspective, shading, texture, shadows, atmospheric attenuation,
stereopsis, ocular accommodation, head motion parallax, and the kinetic depth
effect.
.pp
In this paper, we consider techniques for augmenting equation 1.2 to include as
many depth and shape cues as possible while retaining the computational
efficiency of evaluation using the projection-slice theorem. The key
restriction placed on us is that, to avoid recomputing the 3D frequency domain
representation for each view, the integrand must be independent of viewing
direction. If we wish to support rotation of the volume relative to a light
source, the integrand must also be independent of lighting direction. A class
of shading models that satisfies these constraints has the form for each
viewing ray
.EQ (1.3)
I~~ =~~
sum from {i=0} to n~ w sub i~ left (~
size +3 int from {"-"inf} to inf~ f sub i (t)~ dt~
right )
.EN
for scalars $w sub i$ and scalar volumes $f sub i$. As we shall see, a
surprisingly large class of shading models are expressible as linear
combinations of this form where the $f sub i$'s are independent of both viewing
and lighting directions. In particular, we consider the following three
models:
.ba +0.5i
.ip \(bu 2
X-rays with depth cueing
.ip \(bu 2
X-rays with directional shading
.ip \(bu 2
X-rays with depth cueing & directional shading
.ba -0.5i
.pp
The remainder of the paper is organized as follows. Section 2 describes the
Fourier projection-slice theorem and gives an algorithm for using it to
generate arbitrarily oriented X-rays of a sampled scalar volume. In section 3,
we derive formulations for the three shading models listed above plus a few
extensions. Section 4 describes our implementation and results, and section 5
gives conclusions and suggests directions for future research.
.sh 1 "The projection-slice theorem"
.pp
The Fourier projection-slice theorem states that the inverse 2D Fourier
transform of a slice passing through the origin of the 3D Fourier transform of
a function gives a projected image of the function, where the projection is in
a direction perpendicular to the extracted slice. Each point on the projection
is the integral of the function over the infinite ray passing through that
point and normal to the slice [Dudgeon84]. For later use, we define a Fourier
projection operator $PI sub V$ with input function $f$ on $bold R sup 3$ and
output image $I$ on $bold R sup 2$ such that
.EQ (2.1)
I~~ =~~
down 10 {size +2 PI sub V} (f)~~ =~~
down 20 {size +4 F sub 2 sup {"-"1}}~ down 20 size +4 "{"~
down 20 {size +4 F sub 3}~ down 10 size +3 "{"~ f~ down 10 size +3 "}"~
delta sub V~ down 20 size +4 "}" .
.EN
$delta sub V$ restricts the spectrum to a plane passing through the origin and
perpendicular to viewing direction $bold V$, and $down 20 {size +4 F sub N}$
and $down 20 {size +4 F sub N sup {"-"1}}$ are the forward and inverse
N-dimensional Fourier transforms, respectively.
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 1:"$~~$
Algorithm for generating an attenuation-only rendering of a volume dataset
using the Fourier projection-slice theorem.
.xl 3.125i
.)z
.pp
The Fourier projection-slice theorem is the basis for one variant of computed
tomography (CT). Specifically, if the 2D Fourier transforms of a set of X-rays
of an object taken at different angles are inserted into a 3D spectrum such
that each slice passes through the origin and is perpendicular to the
(parallel) X-ray beam used to acquire that slice, the inverse 3D Fourier
transform of the composite spectrum yields a 3D radiometric facsimile of the
object, i.e. a CT dataset [Macovski83].
.pp
By reversing the CT acquisition process, we can generate X-rays from CT data.
A typical implementation is shown in figure 1. In a preprocessing step, we
compute the 3D discrete Fourier transform (DFT) of a scalar volume sampled on
the three-dimensional integer lattice. For real volumes (i.e. having no
imaginary component), a discrete Hartley transform (DHT) [Bracewell86] may be
substituted. For each view, we extract a slice passing through the origin and
compute its inverse 3D DFT. Since an arbitrarily oriented slicing plane does
not typically pass through points on the integer lattice, slice values must be
computed by interpolation from nearby samples.
.pp
For a volume measuring $N$ voxels on a side where $N$ is an integer power of
two, the computation of a 3D spectrum using a Fast Fourier Transform (FFT) or
Fast Hartley Transform (FHT) algorithm requires $O(N sup 3~ log~ N)$
operations. The cost of extracting a 2D slice measuring $N$ samples on a side
from the 3D spectrum is $O(W sup 3~ N sup 2 )$ where $W$ is the nonzero width
of the interpolating filter. To complete the algorithm, we perform a 2D
inverse FFT or FHT, which requires $O(N sup 2~ log~ N)$ operations. Properly
speaking, sampling theory demands that we employ $2x$ oversampling during slice
extraction. We have found that the penalty of not doing so is minor for the
datasets we typically encounter. Assuming $N~ =~ 256$ and $W~ =~ 4$, we have
slice extraction and inverse transformation costs of $16M$ and $2.25M$
operations, respectively. Therefore, although the Fourier transformation is
asymptotically dominant, the slice extraction step usually takes more time.
.pp
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 2:"$~~$
Imperfect reconstruction in the frequency domain produces incompletely
suppressed replicas in the spatial domain. These appear in the image as
noticeable ghosting.
.xl 3.125i
.)z
The data structures required by the algorithm (in addition to the input volume)
include a 3D spectrum measuring $N$ coefficients on a side and a 2D slice that
also measures $N$ coefficients on a side. While the input volume may be
adequately represented by 8 bits per voxel, accurate frequency domain
representations typically require 32-bit floating point coefficients, thereby
increasing memory requirements by a factor of 4. Quantization or coding of
Fourier coefficients such as specified in the JPEG still picture compression
standard is possible, but reconstruction time must be considered.
.pp
The use of finite support interpolation filters in the slice extraction step
can give rise to objectionable artifacts in Fourier projections. Although
these problems are not the focus of this paper, they impact performance
significantly, so a brief discussion of them is warranted.
.pp
Convolution of a sampled function with a filter of finite support causes
incomplete suppression of periodic replicas in the other domain. Usually,
convolution is performed in the spatial domain and the incompletely suppressed
replicas occur in the frequency domain. In the Fourier projection algorithm,
we convolve in the frequency domain, producing incompletely suppressed replicas
in the spatial domain that appear as noticeable ghosting. The problem is
summarized in one lower dimension in figure 2. Two solutions are known for
this problem:
.np
Employ a better (and typically wider) interpolation filter. If we increase the
nonzero filter width by a factor of $k sub 1$, $k sub 1~ > 1~$, we increase
per-view processing time by a factor of ${k sub 1} sup 3$.
.np
Widen the separation between spatial domain replicas by zero padding the input
volume. If we increase the volume width by a factor of $k sub 2$, which
corresponds to surrounding the data on each of its six sides by a margin of
zeros of width $N~ (k sub 2~ "-"~ 1)/2$, we increase per-view processing time
by a factor of ${k sub 2} sup 2~ (log~ {k sub 2}~ +~ log~ N )~ /~ log~ N$. We
also increase memory requirements by a factor of ${k sub 2} sup 3$ and
preprocessing time by a factor of ${k sub 2} sup 3~ (log~ {k sub 2}~ +~ log~ N
)~ /~ log~ N$, but these effects are less important.
.lp
The second solution has theoretically superior per-view performance. The
relative impact of these two solutions on image quality has not been studied,
however. Malzbender reports satisfactory results with a $5 sup 3$ tap filter
and a 10% margin of zeros on each side of the data [Malzbender]. The images in
this paper were made using a 16% margin ($k sub 2~ =~ 4/3$) and a $4 sup 3$ tap
filter.
.sh 1 "Shading models"
.pp
As discussed in the introduction, the class of integrals that can be directly
evaluated using the Fourier projection-slice theorem does not allow the
modeling of emission and scattering in a participating medium. In this
section, we consider several shading models that are linear combinations of
Fourier projections in the spirit of equation 1.3 and which restore some of the
lost visual cues.
.sh 2 "Depth cueing"
.pp
The fraction of light lost to absorption or scattering in a participating
medium is exponentially related to the distance it travels. This attenuation
provides strong perceptual cues concerning the relative depth of objects in the
scene. In this section, we develop an approximate shading model in which
exponential attenuation with locally varying absorption and scattering
coefficients is replaced by linear attenuation with a constant coefficient. It
will be shown that images based on this model can be rendered efficiently using
the projection-slice theorem.
.pp
Let us begin with exponential attenuation. We define a volume of density $rho
(x,y,z)$ and color $C(x,y,z)$ to emit light with an energy of $C rho$ per unit
length and absorb light with an opacity of $tau~ rho$ per unit length where
$tau$ is a constant. We also define a viewing ray $bold P (t)$ parameterized
by length $t$ and two points $bold P (a)$ and $bold P (b)$ lying on the ray.
Ignoring scattering, the total energy $I sub P$ arriving at $bold P (b)$ due to
emission and absorption along that portion of the ray lying between $bold P
(a)$ and $bold P (b)$ is given by
.ba +0.125i
.EQ L (3.1.1)
I sub P~~ =~~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
exp size +4 [ "-"~ int from a to t~ tau~ rho~ ( bold P (u))~ du~ size +4 ]~
dt .
.EN
.ba -0.125i
The notation is adapted from [Max90]. Replacing exponential attenuation with
locally varying coefficients $tau~ rho$ by linear attenuation with constant
coefficient $tau$ gives
.ba +0.125i
.EQ L (3.1.2)
I sub P~ =~~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
size +1 [ tau /2~ +~ tau /2~
( bold D~ ( bold P (t))~ cdot~ bold V ) size +1 ]~
dt
.EN
.ba -0.125i
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 3:"$~~$
Linear depth cueing for a volume centered at the origin is given by $bold D~ =~
m~ bold P~ +~ b$ for constants $m$ and $b$ where $bold P (t)$ is position along
a viewing ray and $bold V$ is the normalized viewing direction.
.xl 3.125i
.)z
where $bold D~ =~ (D sub x ,D sub y ,D sub z )$ is a linear function of
position along the ray, i.e. $bold D~ =~ m~ bold P~ +~ b$ for constants $m$ and
$b$, and $bold V~ =~ (V sub x ,V sub y ,V sub z )$ is the normalized viewing
direction as shown in figure 3. Note that $bold V$ points away from the eye,
not towards it as is common in many textbooks. We have chosen this convention
in order to underscore the similarity between equation 3.1.2 and equation 3.2.3
in section 3.2.
.pp
In the introduction, we stated that in order to avoid recomputing the 3D
frequency domain representation for each view, all integrals over the volume
must be independent of viewing direction. Looking back at the last equation,
we observe that viewing direction $bold V$ is independent of ray parameter $t$.
We may therefore factor it out of the integral, producing
.ba +0.125i
.EQ L (3.1.3)
I sub P~ mark =~~
V sub x~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
size +1 [ tau /2~ +~ tau /2~
(D sub x ( bold P (t)) ) size +1 ]~
dt
.EN C
.EQ L
lineup +~~
V sub y~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
size +1 [ tau /2~ +~ tau /2~
(D sub y ( bold P (t)) ) size +1 ]~
dt
.EN C
.EQ L
lineup +~~
V sub z~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
size +1 [ tau /2~ +~ tau /2~
(D sub z ( bold P (t)) ) size +1 ]~
dt
.EN C
.EQ L
lineup +~~
c sub 1~
size +3 int from a to b~
rho~ ( bold P (t))~ C( bold P (t))~
dt .
.EN
.ba -0.125i
where $c sub 1~ =~ tau /2~ "-"~ tau /2~ ( V sub x + V sub y + V sub z )$.
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 4:"$~~$
Algorithm for generating a depth cued rendering of a volume dataset. During
precomputation, volume $rho~ C$ is multiplied by depth operator $bold D$. For
each view, the extracted slices are multiplied by viewing direction $bold V$.
.xl 3.125i
.)z
.pp
These four integrals are now in a form that can evaluated using the Fourier
projection-slice theorem. Applying the projection operator defined in
equation 2.1, we obtain the following expression for output image $I$
.ba +0.125i
.EQ L (3.1.4)
I~~ mark =~~
V sub x~
PI sub V size +1 [ ( tau /2~ +~ tau /2~ D sub x (x,y,z) )~
rho~ (x,y,z)~ C(x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
V sub y~
PI sub V size +1 [ ( tau /2~ +~ tau /2~ D sub y (x,y,z) )~
rho~ (x,y,z)~ C(x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
V sub z~
PI sub V size +1 [ ( tau /2~ +~ tau /2~ D sub z (x,y,z) )~
rho~ (x,y,z)~ C(x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
c sub 1~
PI sub V size +1 [ rho~ (x,y,z)~ C(x,y,z) size +1 ] .
.EN
.ba -0.125i
.pp
The algorithm is shown in figure 4. Four copies of $rho~ (x,y,z)~ C(x,y,z)$
are created. Each is multiplied by the appropriate depth cueing function and
Fourier transformed, yielding four 3D spectra. The first three spectra
represent volumes that are depth cued along one of object space coordinate
axes. The fourth spectrum represents the original volume. For each view,
slices are extracted from each spectra, inverse Fourier transformed, weighted
by the appropriate component of the current viewing direction as shown in the
figure and summed, producing an image that is depth cued for that viewing
direction.
.pp
To summarize, we have replaced per-voxel shading followed by projection with
projection followed by per-pixel shading. This transformation is allowed
because integration is a linear operator and the current viewing direction is
constant for all voxels. Both formulations produce identical results, but the
latter allows a more efficient implementation using Fourier projections.
.sh 2 "Directional shading"
.pp
In the previous section, we showed that a per-voxel dot product followed by
projection could be replaced with projection followed by a per-pixel dot
product. There are other places where dot products appear in per-voxel shading
models. Most notable is probably the $| bold N~ cdot~ bold L |$ product giving
the irradiance on a surface of orientation $bold N$ produced by a parallel
light source originating from direction $bold L$. Unfortunately, the presence
of the nonlinear absolute value operator prevents evaluation of this dot
product using the Fourier projection-slice theorem. In this section, we
develop an alternative, physically valid, shading model in which illumination
from a parallel light source is replaced by illumination from a hemispherical
source. It will be shown that images based on this model can be rendered
efficiently using the projection-slice theorem.
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 5:"$~~$
Hemispherical illumination of a surface. Irradiance is equal to integral of
the projection of shaded portion of the hemisphere down onto the plane
containing the surface.
.xl 3.125i
.)z
.pp
Let us begin with Lambertian reflection and parallel illumination. We define a
volume of albedo $OMEGA (x,y,z)$ and an orientation field $bold N~ (x,y,z)$.
The latter is typically estimated from the scalar input data using a central
difference operator [Levoy88]. Ignoring attenuation, the energy $I sub P$
arriving at point $bold P (b)$ along viewing ray $bold P (t)$ is given by
.EQ (3.2.1)
I sub P~ =~~
size +3 int from a to b~
I sub s~ OMEGA~ ( bold P (t))~
| bold N~ ( bold P (t))~ cdot~ bold L |~
dt
.EN
where $I sub s$ is the intensity of a parallel light source originating from
direction $bold L$. This expression assumes illumination in both the $bold L$
and $"-" bold L$ directions. One-sided illumination is obtained by replacing
$| bold N~ cdot~ bold L |$ with $max ( bold N~ cdot~ bold L ,~ 0)$.
.pp
Since the absolute value (or max) operator is nonlinear, we cannot factor
lighting direction $bold L$ out of the integral as we did with viewing
direction $bold V$ in section 3.1. Omitting these operators entirely would
cause voxels containing backward-facing normals (relative to the light source)
to produce negative light that would combine under integration with positive
light produced by voxels containing forward-facing normals to yield incorrect
images.
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 6:"$~~$
Algorithm for generating a directionally shaded rendering of a volume dataset.
During precomputation, volume $I sub s~ OMEGA$ is multiplied by orientation
field $bold N$. For each view, the extracted slices are multiplied by lighting
direction $bold L$.
.xl 3.125i
.)z
.pp
To eliminate this problem, we consider an alternative illumination method - an
emissive hemisphere. The irradiance $E sub i$ on a surface having normal
vector $bold N$ illuminated by a hemisphere whose pole points in direction
$bold L$ as shown in figure 5 is equal by Nusselt's analogue (as described in
[Cohen85]) to the projection of the visible portion of the hemisphere down onto
the plane containing the surface, or
.ba +0.5i
.EQ L (3.2.2)
E sub i~~ =~~
I sub s~
size +3 int from {"-" pi over 2~ +~ theta} to {pi over 2}~~~
size +3 int from {"-" pi over 2} to {pi over 2}~
cos phi~ cos sup 2 psi~
d psi~
d phi
.EN C
.EQ L
~~~~~ =~~
pi /2~ +~ pi /2~ cos psi
.EN C
.EQ L
~~~~~ =~~
pi /2~ +~ pi /2~ ( bold N~ cdot~ bold L ) .
.EN
.ba -0.5i
Assuming Lambertian reflection with albedo $OMEGA~ (x,y,z)$, we have
.ba +0.25i
.EQ L (3.2.3)
I sub P~ =~~
size +3 int from a to b~
I sub s~ OMEGA~ ( bold P (t))~
size +1 [ pi /2~ +~ pi /2~
( bold N~ ( bold P (t) )~ cdot~ bold L ) size +1 ]~
dt .
.EN
.ba -0.25i
Although seldom seen in the graphics literature (a notable exception is
[Nishita86]), this shading model is entirely plausible - it models the
appearance of diffuse objects outdoors on a cloudy day (ignoring
interreflection among objects on the ground). The directional shading provided
by this model is less sensitive to surface orientation than shading arising
from parallel illumination, just as our perception of shape is poorer on a
cloudy day than on a sunny day. Nevertheless, the model is a significant
improvement over ignoring scattering effects entirely.
.pp
Continuing our derivation, we observe that equation 3.2.3 is nearly identical
to equation 3.1.2. If we wish to support rotation of the volume relative to
the light source yet avoid recomputing the 3D frequency domain representation
for each view, the integrand must be independent of lighting direction.
Looking at the equation, we see that hemisphere pole direction $bold L$ is
independent of ray parameter $t$ and may be factored out of the integral. As
in the case of depth cueing, this manipulation produces four integrals that can
each be evaluated using the Fourier projection-slice theorem. Repeating the
steps in section 3.1, we obtain the following expression for output image $I$
.ba +0.125i
.EQ L (3.2.4)
I~~ mark =~~
L sub x~
PI sub V size +1 [ ( pi /2~ +~ pi /2~ N sub x (x,y,z) )~
I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
L sub y~
PI sub V size +1 [ ( pi /2~ +~ pi /2~ N sub y (x,y,z) )~
I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
L sub z~
PI sub V size +1 [ ( pi /2~ +~ pi /2~ N sub z (x,y,z) )~
I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
c sub 2~
PI sub V size +1 [ I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN
.ba -0.125i
where $N sub x (x,y,z)~$ is the $X$-component of the orientation field at
location $(x,y,z)$ similarly for $N sub y$ and $N sub z$, and
$c sub 2~ =~ pi /2~ "-"~ pi /2~ ( L sub x + L sub y + L sub z )$.
.pp
The algorithm is shown in figure 6. Four copies of $I sub s~ OMEGA~ (x,y,z)$
are created. Each is multiplied by the appropriate component of $bold N~
(x,y,z)$ and Fourier transformed, yielding four 3D spectra. Because of the
method used to compute $bold N$, the first three spectra essentially represent
volumes that have been gradient-enhanced in one of the object space coordinate
directions. The fourth spectra represents the original volume. For each
change in view or lighting direction, slices are extracted from each spectra,
inverse Fourier transformed, weighted by the appropriate component of the
current lighting direction as shown in the figure and summed, producing an
image that is directionally shaded for that lighting direction.
.pp
As in the case of depth cueing, we have replaced per-voxel shading followed by
projection with projection followed by per-pixel shading because one of the
terms in a dot product - the current lighting direction - is constant for all
voxels.
.sh 2 "Depth cueing & directional shading"
.pp
If depth cueing and directional shading as defined in sections 3.1 and 3.2 can
each be expressed as linear combinations of Fourier projections, it should be
possible to combine them, producing a composite shading model that exhibits
both effects and which is expressible as a linear combination of Fourier
projections.
.pp
Multiplying the integrand in equation 3.2.3 by the linear attenuation factor in
equation 3.1.2, we obtain for a viewing ray
.ba +0.0625i
.EQ L (3.3.1)
I sub P~ =~~
size +3 int from a to b~
I sub s~ OMEGA~ ( bold P (t))~
size +1 [ pi /2~ +~ pi /2~
( bold N~ ( bold P (t)~ cdot~ bold L ) size +1 ]
.EN C
.EQ L
~~~~~ roman "x"~~
size +1 [ tau /2~ +~ tau /2~
( bold D~ ( bold P (t))~ cdot~ bold V ) size +1 ]~
dt .
.EN
.ba -0.0625i
After we factor and apply the Fourier projection operator,
we obtain for the output image
.EQ L (3.3.2)
I~~ mark =~~
tau~ pi /2~ bold {V sub L}~ cdot~
PI sub V size +1 [ ( 1/2~ +~ 1/2~ bold {D sub N}~ (x,y,z) )~
I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
tau~ pi /2~ c sub 3~
PI sub V size +1 [ I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN
for $c sub 3~ =~ 1/2~ "-"~ 1/2~ up 15 {sum from i}~ bold {V sub L} sub i$ where
.ba +0.0625i
.EQ L
bold {V sub L}~~ =~~ ( mark
V sub x L sub x ,~ V sub x L sub y ,~ V sub x L sub z ,~
V sub y L sub x ,~ V sub y L sub y ,~ V sub y L sub z ,~
.EN C
.EQ L
lineup
V sub z L sub x ,~ V sub z L sub y ,~ V sub z L sub z ,
V sub x ,~ V sub y ,~ V sub z ,~
L sub x ,~ L sub y ,~ L sub z )
.EN
.ba -0.0625i
and
.ba +0.0625i
.EQ L
bold {D sub N}~~ =~~ ( mark
D sub x N sub x ,~ D sub x N sub y ,~ D sub x N sub z ,~
D sub y N sub x ,~ D sub y N sub y ,~ D sub y N sub z ,~
.EN C
.EQ L
lineup
D sub z N sub x ,~ D sub z N sub y ,~ D sub z N sub z ,~
D sub x ,~ D sub y ,~ D sub z ,~
N sub x ,~ N sub y ,~ N sub z ) .
.EN
.ba -0.0625i
This requires the extraction and inverse Fourier transformation of 16 slices
from 16 precomputed spectra. While this algorithm is still asymptotically
less expensive than conventional volume rendering, its high constant makes it
slower except for very large volumes. We note that our composite shading
model ignores emission $rho~ C$ in equation 3.1.2. In the context of a
directional shading model, this term can be thought of as representing ambient
illumination. If included, it adds another 3 terms to equation 3.3.2 and hence
another 3 slices that must be extracted.
.sh 2 "Specular reflections and other extensions"
.pp
Shiny surfaces exhibit a greater change in reflected light intensity for small
variations in surface orientation than dull surfaces. A common formulation of
the specular term in Phong's illumination model [Phong75] is $| bold N~ cdot~
bold H | sup n$ for some exponent $n$ where $bold H$ is the vector halfway
between the directions of the light source, $bold L$, and the viewer, $bold V$.
Using hemispherical illumination to avoid the absolute value operator, we have
for a viewing ray
.ba +0.0625i
.EQ L (3.4.1)
I sub P~ =~~
size +3 int from a to b~
I sub s~ OMEGA~ ( bold P (t))~
size +1 [ pi /2~ +~ pi /2~
( bold N~ ( bold P (t)~ cdot~ bold H ) sup n size +1 ]~
dt .
.EN
.ba -0.0625i
After we factor and apply the Fourier projection operator,
we obtain for the output image
.ba +0.0625i
.EQ L (3.4.2)
I~~ mark =~~
bold {H sub n}~ cdot~
PI sub V size +1 [ ( pi /2~ +~ pi /2~ bold {N sub n} (x,y,z) )~
I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN C
.EQ L
lineup +~~
c sub 4~
PI sub V size +1 [ I sub s~ OMEGA~ (x,y,z) size +1 ]
.EN
.ba -0.0625i
for suitable constant $c sub 4$ where $bold {H sub n}$ and $bold {N sub n}$
contain the list of factors of the polynomial $( bold N~ cdot~ bold H ) sup n$.
The number of terms and hence the number of spectra is given by the multinomial
theorem [Knuth73]. For $n~ =~ 10$, a typical value, we require 286 precomputed
spectra, making this extension slower than conventional volume rendering even
for very large volumes.
.pp
For scenes of high depth complexity, linear attenuation may provide an
insufficiently steep intensity falloff to disambiguate depth relationships.
Equation 3.4.1 can be adapted to provide $n$-degree polynomial depth cueing by
replacing $bold N$ and $bold H$ with $bold D$ and $bold V$, respectively, and
adjusting the constants. For example, inverse square law falloff would be
given by $( bold D~ cdot~ bold V ) sup 2$ and would require 10 slices from 10
precomputed spectra.
.sh 1 "Implementation and results"
.pp
Figures 7 through 13 each show the application of a different shading model to
a $96 sup 3$ voxel cube extracted from a CT scan of a human skull mounted in a
lucite head cast. The right eye orbit and a portion of the zygoma (temple) are
visible at center and left, respectively, and the lucite nose is visible at
lower right. To facilitate comparisons, figure 7 is a conventional volume
rendering generated using the techniques described in [Levoy88].
.pp
Figures 8 through 13 were generated using the Fourier projection-slice theorem
as follows. The input data was first edge enhanced and nonlinear windowed
(except for figure 8 as noted below), then surrounded with a 16% margin of
zeros to avoid the ghosting problem discussed in section 2. The resulting $128
sup 3$ voxel cube was weighted according to equations 3.1.4, 3.2.4, or 3.3.2,
and Fourier transformed using a 3D FFT algorithm. For each view, slices of
width equal to the spectra were extracted using a tricubic interpolating
spline, inverse transformed using a 2D FFT algorithm, weighted appropriately
for the current viewing or lighting directions, and summed to produce a $128
sup 2$ image.
.pp
Memory requirements varied from 48 megabytes for equation 3.1.4 to 256
megabytes for equation 3.3.2. Per-view rendering time varied from 6 seconds
for equation 3.1.4 (about 1.5 seconds per extracted slice) to over a minute for
equation 3.3.2 (due mostly to paging). Timings are for a 36 MHz Silicon
Graphics 4D/320 VGX and a sequential C implementation. No attempt was made to
optimize the code.
.pp
By replacing the FFT with an FHT, one can immediately halve the per-view
rendering time and the memory requirements. By optimizing the filter and
trimming the zero padding, another factor of two speedup and some reduction in
memory usage can be expected. For greater speedups, the algorithm can be
parallelized at either a fine grain (within the FFT) or a coarse grain (between
slices). It is also amenable to hardware acceleration using video coprocessors
or video digital signal processors (DSP). These devices frequently operate on
fixed-point representations of frequency spectra, suggesting additional
opportunities for reducing memory costs.
.pp
Figure 8 shows an attenuation-only projection (i.e. an X-ray) generated using a
variant of equation 1.2. The lack of occlusion and shape-from-shading cues is
evident. Figure 9 shows an attenuation-only projection of the input volume
after edge sharpening and nonlinear windowing of the density range. Figures
10 through 13 also include sharpening and windowing. Figure 10 shows a
depth-cued X-ray generated using the shading model described in section 3.1.
The weighting assigned to each point in the volume falls off linearly with
increasing distance from the viewer. Figure 11 shows a directionally shaded
X-ray generated using the shading model described in section 3.2. Each point
in the volume is weighted as if it sat on a diffusely reflecting surface
illuminated by a hemisphere placed on the right side of the picture and shining
leftwards. Figures 12 and 13 show two views that each include both depth
cueing and directional shading using the model described in section 3.3.
Although no occlusion occurs in either of these figures, the shading effects
enhance our ability to perceive shape and discern spatial relationships. When
displayed in motion, the visual impact is sufficiently strong that viewers
often overlook the fact that these are X-rays, not opaque renderings.
.sh 1 "Conclusions and future work"
.pp
We have described a family of algorithms for generating realistically shaded
renderings of volume data using the Fourier projection-slice theorem. While
the images produced using these algorithms do not exhibit occlusion, they
provide sufficient depth and shape information to make them moderately
understandable in static views and very understandable in motion sequences. In
either case, they produce images richer in visual cues than the
attenuation-only projections with which the projection-slice theorem is usually
associated.
.pp
One aspect of the Fourier approach that we have not discussed is how it
impacts interactive data segmentation. The edge sharpening and nonlinear
windowing operators used in figures 9 through 13 are applied before the forward
3D Fourier transform. As a consequence, changing them requires recomputing the
spectra. Depending on the application, it may be undesirable to hardwire data
segmentation in this way.
.pp
Several strategies are available for avoiding this difficulty. If the data
segmentation operator is expressible as a linear combination of basis
functions, it can be handled using the same scheme we have employed for shading
models. For example, if a CT dataset is segmented into fat, muscle, and bone
volumes during preprocessing, yielding three spectra, the weights and colors
assigned to each material can be altered at image generation time by applying
different weights to the slices extracted from each spectrum.
.pp
If the segmentation operator is expressible as a convolution, it may be
possible to implement it equally or more efficiently as multiplication in the
frequency domain. For example, the edge sharpening operator used in figures 9
through 13 could have been implemented as multiplication of each extracted
slice by a high-frequency emphasis function. Using this approach, the amount
and type of enhancement can be specified at image generation time.
.pp
Conversely, if the segmentation operator is expressible as multiplication by a
function of position in the spatial domain, it may be possible to implement it
as convolution in the frequency domain. For an intriguing example of this
strategy, consider the linear depth cueing described in section 3.1. Although
the Fourier transform of this function is infinite in extent, equally effective
depth cueing can be obtained using the first half period of a suitably
oriented cosine function. The Fourier transform of a cosine is two delta
functions. For each view, we determine the placement of these deltas,
preconvolve them with our interpolation filter (which amounts to summing two
shifted copies of the filter function), and use the composite filter during
slice extraction. This technique requires only one spectrum as input, making
it potentially more efficient than the technique described in this paper. We
are currently investigating this approach.
.pp
Looking beyond shading models, sums of Fourier projections can also be used to
compute integrals over irregular spatial regions. Specifically, the integral
of any function over the ray segment connecting two points in the interior of a
volume can be computed as a sum of integrals over spatially disjoint intervals
of the ray that taken together span the segment. Assuming well behaved
functions, each integral can be computed using the projection-slice theorem.
We are currently investigating an algorithm based on recursive subdivision of
volumes into octants that can compute images in $O(N sup 2~ log sup 2~ N)$
time. Possible applications include rendering of spatially clipped volumes for
medical visualization, fast calculation of inter-element attenuation for zonal
radiosity algorithms, and approximate visibility determination for
geometrically defined scenes.
.ps 10
.uh "Acknowledgements"
.ps 9
.lp
This research was supported by the National Science Foundation (NSF), the
National Aeronautics and Space Administration (NASA), and the sponsoring
companies of the Stanford Center for Integrated Systems (CIS). The SGI 4D/320
VGX workstation was donated by Silicon Graphics, Inc. The CT scan was provided
by North Carolina Memorial Hospital. The notion that shading could be factored
with respect to digital compositing was suggested to me by Brice Tebbs.
Discussions with Adam Levinthal were useful in the early stages of this
project, Takashi Totsuka suggested the frequency domain implementation of
depth cueing, and Alice Yu provided a helpful reading of the manuscript.
>>>>> Ackowledgement of Brice Tebbs did not appear in published paper or TR.
.ps 10
.uh "References"
.ps 9
.in 0.35i
.ne 0.5i
.ti 0
[Bracewell86] Bracewell, R.,
.i "The Fourier Transform and Its Applications" ,
McGraw-Hill, 1986.
.ne 0.5i
.ti 0
[Cohen85] Cohen, M.F. and Greenberg, D.P.,
``The Hemi-Cube: A Radiosity Solution for Complex Environments,''
.i "Computer Graphics (Proc. SIGGRAPH '85)" ,
Vol. 19, No. 3, July, 1985, pp. 31-40.
.ne 0.5i
.ti 0
[Drebin88] Drebin, R.A., Carpenter, L., and Hanrahan, P.,
``Volume Rendering,''
.i "Computer Graphics (Proc. SIGGRAPH '88)" ,
Vol. 22, No. 4, August, 1988, pp. 65-74.
.ne 0.5i
.ti 0
[Dudgeon84] Dudgeon, D.E. and Mersereau, R.M.,
.i "Multidimensional Digital Signal Processing" ,
Prentice-Hall, 1984.
.ne 0.5i
.ti 0
[Dunne90] Dunne, S., Napel, S. and Rutt, B.,
``Fast Reprojection of Volume Data,''
.i "Proc. First Conference on Visualization in Biomedical Computing,"
IEEE Computer Society Press, May, 1990, pp. 11-18.
.ne 0.5i
.ti 0
[Hottel67] Hottel, H.C. and Sarofim, A.F.,
.i "Radiative Transfer" ,
McGraw-Hill, 1967.
.ne 0.5i
.ti 0
[Knuth73] Knuth, D.,
.i "The Art of Computer Programming" ,
Addison-Wesley, 1973.
.ne 0.5i
.ti 0
[Krueger90] Krueger, W.,
``Volume Rendering and Data Feature Enhancement,''
.i "Computer Graphics (Proc. San Diego Workshop on Volume Visualization)" ,
Vol 24, No. 5, November, 1990, pp. 21-26.
.ne 0.5i
.ti 0
[Levoy88] Levoy, M.,
``Display of Surfaces from Volume Data,''
.i "IEEE Computer Graphics and Applications" ,
Vol. 8, No. 3, May, 1988, pp. 29-37.
.ne 0.5i
.ti 0
[Macovski83] Macovski, A.,
.i "Medical Imaging Systems" ,
Prentice-Hall, 1983.
.ne 0.5i
.ti 0
[Malzbender] Malzbender, T.,
``Fourier Volume Rendering.''
Submitted for publication.
.ne 0.5i
.ti 0
[Max90] Max, N., Hanrahan, P. and Crawfis, R.,
``Area and Volume Coherence for Efficient Visualization
of 3D Scalar Functions,''
.i "Computer Graphics (Proc. San Diego Workshop on Volume Visualization)" ,
Vol. 24, No. 5, November, 1990, pp. 27-33.
.ne 0.5i
.ti 0
[Napel91] Napel, S., Dunne, S. and Rutt, B.K.,
``Fast Fourier Projection for MR Angiography,''
.i "Magnetic Resonance in Medicine" ,
Vol. 19, 1991, pp. 393-405.
.ne 0.5i
.ti 0
[Nishita86] Nishita, T. and Nakamae, E.,
``Continuous Tone Representation of Three-Dimensional Objects
Illuminated by Sky Light,''
.i "Computer Graphics (Proc. SIGGRAPH '86)" ,
Vol. 20, No. 4, August, 1986, pp. 125-132.
.ne 0.5i
.ti 0
[Phong75] Bui-Tuong, Phong,
``Illumination for Computer-Generated Pictures,''
.i "Communications of the ACM."
Vol. 18, No. 6, June, 1975, pp. 311-317.
.ne 0.5i
.ti 0
[Porter84] Porter, T. and Duff, T.,
``Compositing Digital Images,''
.i "Computer Graphics (Proc. SIGGRAPH '84)" ,
Vol. 18, No. 3, July, 1984, pp. 253-259.
.ne 0.5i
.ti 0
[Rushmeier87] Rushmeier, H.E. and Torrance, K.E.,
``The Zonal Method for Calculating Light Intensities
in the Presence of a Participating Medium,''
.i "Computer Graphics (Proc. SIGGRAPH '87)" ,
Vol. 21, No. 4, July, 1987, pp. 293-302.
.ne 0.5i
.ti 0
[Sabella88] Sabella, P.,
``A Rendering Algorithm for Visualizing 3D Scalar Fields,''
.i "Computer Graphics (Proc. SIGGRAPH '88)" ,
Vol. 22, No. 4, August 1988, pp. 51-58.
.ne 0.5i
.ti 0
[Siegel81] Siegel, R. and Howell, J.R.,
.i "Thermal Radiation Heat Transfer" ,
Hemisphere Publishing, 1981.
.(z L
.sp 2.75i
.xl 2.875i
.ip " " 0.125i
.b "Figure 7:"$~~$
Conventional volume rendering of a CT scan of a human skull mounted in
a lucite head cast.
.xl 3.125i
.)z
.bp
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 8:"$~~$
Attenuation-only projection (i.e. X-ray). All figures on this page were
generated using the Fourier projection-slice theorem.
.xl 3.125i
.)b
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 10:"$~~$
Depth-cued X-ray generated using the shading model described in section 3.1.
.xl 3.125i
.)b
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 12:"$~~$
Depth-cued and directionally shaded X-ray generated using the shading model
described in section 3.3.
.xl 3.125i
.)b
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 9:"$~~$
Attenuation-only projection after edge sharpening and nonlinear windowing.
Figures 10-13 also include these enhancements.
.xl 3.125i
.)b
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 11:"$~~$
Directionally shaded X-ray generated using the shading model described in
section 3.2.
.xl 3.125i
.)b
.(b L
.sp 2.4i
.xl 2.875i
.ip " " 0.125i
.b "Figure 13:"$~~$
Rotated view of depth-cued and directionally shaded X-ray.
.xl 3.125i
.)b