IEEE/EG International Symposium on Volume Graphics (2010)

R. Westermann and G. Kindlmann (Editors)

Advanced Light Material Interaction

for Direct Volume Rendering

Florian Lindemann and Timo Ropinski

Visualization and Computer Graphics Research Group (VisCG),

Department of Computer Science, University of Münster, Germany

Abstract

In this paper we present a heuristic approach for simulating advanced light material interactions in the context
of interactive volume rendering. In contrast to previous work, we are able to incorporate complex material func-
tions, which allow to simulate reflectance and scattering. We exploit a common representation of these material
properties based on spherical harmonic basis functions, to combine the achieved reflectance and scattering ef-
fects with natural lighting conditions, i. e., incorporating colored area light sources. To achieve these goals, we
introduce a modified SH projection technique, which is not just tailored at a single material category, but adapts
to the present material. Thus, reflecting and scattering materials as assigned trough the transfer function can be
captured in a unified approach. We will describe the required extensions to the standard volume rendering integral
and present an approximation which allows to realize the material effects in order to achieve interactive frame
rates. By exploiting a combination of CPU and GPU processing, we are able to modify material properties and
can change the illumination conditions interactively. We will demonstrate the outcome of the proposed approach
based on renderings of real-world data sets and report the achieved computation times.
Categories and Subject Descriptors (according to ACM CCS): I.3.3 [Computer Graphics]: Picture/Image
Generation—Line and curve generation

1. Introduction ever, only by considering complex reflectance functions in

combination with area light sources, realistic illumination
In recent years sophisticated volumetric illumination mod- effects can be simulated. By not considering advanced ma-
els for direct volume rendering (DVR) have been proposed. terial properties as well as area light sources, the quality of
Interactive techniques are available, that support ambient oc- DVR images is rather limited compared to polygonal com-
clusion [RMSD∗ 08, SPH∗ 09, LHY09], shadows [HKSB06] puter graphics.
as well as scattering [KPH∗ 03, Sal07]. While all these tech-
niques brought realism in the are of volume graphics to a Besides the material properties, light material interactions
new level, an essential part needed to allow even more re- depend on the light source geometry, i. e., its extend and di-
alistic image synthesis has only attracted little attention in rection, and in some cases also the viewing direction. There-
the area of DVR so far: light material interaction. The ap- fore, changing the camera or the light source requires a re-
pearance of an illuminated object depends on the combina- computation. In this paper we introduce an approach for in-
tion of material properties, lighting properties as well as the tegrating different material types into the volume render-
light transport. While these three factors are often consid- ing process and illuminate them using area light sources.
ered when illuminating polygonal objects, volume rendering With the presented approach, it becomes possible to simulate
techniques so far mainly focus on modeling the light trans- realistic materials of different kinds interactively. We will
port. Regarding the material and the lighting properties only demonstrate how to incorporate complex reflectance and
little work has been done, i. e., materials are usually reduced phase functions as well as area light sources. By transform-
to assigning an emissive color through the transfer function ing all these properties into a common basis representation,
and a point light source is used to illuminate the scene. How- we are able to perform all computations required for render-

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

ing on the GPU, and are able to achieve interactive frame tation, Ljung et al. exploit a ray-casting approach to deter-
rates, by still allowing to change the light properties and to a mine the ambient illumination [LHY09]. Both of these ap-
certain degree also the material properties. The key idea is to proaches can be applied to DVR as well as isosurface ren-
exploit a spherical harmonic (SH) basis function representa- dering. In contrast, Wyman et al. [WPHS06] and Beason
tion for the material properties as well as the light sources. et al. [BGB∗ 06] focus on a precomputation, which allows
This allows us to integrate colored area light sources, and to render isosurfaces under static illumination conditions.
simulate their interaction with advanced materials interac- More recently, Banks and Beacon have demonstrated how to
tively, i. e., the light direction as well as the camera position decouple illumination sampling from isosurface generation,
can be changed. To our knowledge, this is the first attempt, in order to achieve sophisticated isosurface illumination re-
which incorporates colored area light sources into DVR. Fur- sults [BB09].
thermore, we are not aware of other techniques dealing with
SH basis functions have been used in many areas of polyg-
advanced material properties in the area of interactive vol-
onal rendering. They allow to approximate spherical func-
ume rendering, in the sense that complex reflectance func-
tions by exploiting a set of SH basis functions. Among
tions as well as scattering functions can be used. Addition-
other applications, this concept has been exploited to cap-
ally, we exploit an adapted SH projection approach, which
ture BRDFs [WAT92], as well as to allow the integration of
allows to handle different material categories, while allow-
low frequency shadowing effects [SKS02]. It could also be
ing an easy integration of conventional SH lighting [SKS02].
shown, that SH basis functions can be used for other render-
ing effects, as caustics or scattering [KH84]. Ritschel has
2. Related Work applied SHs in the field of volume rendering to simulate
low-frequency shadowing [Rit07]. Similar as in this paper,
In the past, several interactive volumetric illumination mod- a GPU-based approach is exploited in order to compute the
els have been proposed. Many of these techniques are based required SH coefficients. The ray traversal with increasing
on the optical models initially derived by Max [Max95]. step size as proposed by Ritschel, has also been integrated
One important approach in this area is the scattering tech- into our implementation.
nique presented by Kniss et al. [KPH∗ 03]. By restricting
themselves to a forward scattering phase function, the au- As can be seen, all cited volume rendering papers rather
thors are able to achieve convincing results at interactive target light transport, than light material interaction. How-
frame rates. The idea of constraining the phase function ever, we believe that this area needs to be covered in order
has been picked up by other researchers. Schott et al. ex- to generate more convincing renderings. To our knowledge
tend this idea in order to integrate directional occlusion ef- only the style transfer function approach for illustrative vi-
fects into a slice-based volume renderer [SPH∗ 09]. Patel et sualization by Bruckner and Groeller [BG07] deals with ad-
al. have adapted the concept to not only incorporate occlu- vanced material properties in the area of interactive volume
sion information, but to also simulate light emitting struc- rendering so far.
tures [PBVG10]. However, since the techniques by Schott
et al. as well as Patel et al. are based on a slice-wise front- 3. Light Material Interaction
to-back compositing, only headlights are possible. Ropinski
et al. also use cone-shaped phase functions in order to sim- In the following, we will incorporate advanced material
ulate scattering effects together with an interactive volume properties as well as colored area light sources into inter-
processing technique, were light is swept through an illu- active volume rendering. We will refer to reflectance and
mination volume [RDRS10]. While all these techniques are phase functions as material functions, and will call materials
capable to produce convincing results, they are rather strong having significant differences as belonging to other material
approximations, since scattering is assumed to be directed categories. When considering material and lighting proper-
along one major axis, and not over the whole sphere. This is ~ o ) which
ties together with light transport, the radiance L(~x, ω
especially problematic, when dealing with area light sources leaves from position~x inside a volume in direction ω ~ o can be
having multiple spatial components. An alternative approach defined as:
has been introduced by Rezk-Salama [Sal07]. He proposes
the use of Monte-Carlo integration to allow scattering effects Z
on a limited set of surfaces. ~ o ) = EM (~x, ω
L(~x, ω ~o) + ~ o ) · Li (~x, ω
~ i, ω
M(~x, ω ~ i )d ω
~ i.
Ω
Besides the integration of scattering effects, researchers
have also targeted the simulation of diffuse interreflections. Here, L(~x, ω~ o ) is the sum of the emission EM (~x, ω
~ o ) and
As with the related work regarding scattering, we cover only the local light material interaction, which is expressed by
those publications targeting interactive volume rendering. the integral over the unit sphere Ω. To compute the local
Ropinski et al. use a set of local histograms in order to cap- light material interaction, the material function M(~x, ω ~ i, ω
~o)
ture the vicinity of a voxel [RMSD∗ 08]. While their tech- ~ i ) coming from
is modulated by the incident intensity Li (~x, ω
nique is based on a rather expensive histogram precompu- ~ i . Since we deal with volume data, where the
direction ω

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

notion of a surface does not inherently exist, we do not Intuitively, the Ylm (ϑ , ϕ) can be thought of having two
weight this modulation by the clamped cosine term of the independent components to span the sphere. Complex ex-
surface normal and ωi . While this weighting could be eas- ponential functions with a single argument of period 2π are
ily achieved by using the gradient 5τ( f (~x)), it would not used to define the variations along the latitude, while asso-
compromise with the definition of the material function ciated Lengendre polynomials are used to capture variations
M(~x, ω~ i, ω~ o ), which is in contrast to surface reflectance func- along the longitudes. To allow an efficient computation and
tions not defined over the hemisphere, but over the whole representation of f , the representation can be bandlimited:
sphere. So far, only local effects, such as emission and
local light material interactions, i. e., out-scattering or re-
L l
flectance, have been incorporated. To integrate environmen-
f˜(ϑ , ϕ) = ∑ ∑ clmYlm (ϑ , ϕ)
tal effects, as in-scattering or color bleeding, the unmod- l=0 m=−l
ified incident intensity Li (~x, ω ~ i ) has to be modulated with
~ i ), which is based on the environment of ~x:
Lmi (~x, ω In this bandlimited version, the number of summands is
equal to 1 + 3 + ... + (2L + 1) = (L + 1)2 , which directly
R influences the number of required ci j . For practical rea-
~ o ) = EM (~x, ω
L(~x, ω ~o) + ~o) ·
~ i, ω
Ω M(~x, ω sons we have set L = 3 throughout this paper, resulting in
~ i ) · Li (~x, ω
(Lmi (~x, ω ~ i ))d ω
~ i. (1) (3 + 1)2 = 16 coefficients. To project f and compute the
ci j , we use a least square projection, which minimizes the
In Section 4 we will describe how to compute Lmi in order squared error between our approximation and the original
to integrate diffuse interreflections as well as in-scattering. function f . Thus, the SH coefficients c j can be computed as
In the simplified model described by Equation 1, light is only follows:
scattered when traveling towards~x, such that no light is scat-
tered along the viewing direction. Please note, that we use Z
a slightly modified notion for L as compared to Li and Lmi . cj = Y j (s) f (s) ds. (2)
S
While L is defined based on the direction ω going out from
~x, Li and Lmi are defined based on the direction ω going in Due to the orthonomality of the SH basis functions, it can
to ~x. be shown, that the integration of two SH-projected functions
f and g becomes a sum of the products of their coefficients
With the simplified illumination model specified in Equa-
ci and di :
tion 1, we are able to incorporate complex material functions
specified through M. M, Lmi as well as the incident light Li
are defined over the unit sphere Ω. Obviously, the integration n
over Ω involves a significant computing effort, and cannot I( f˜, g̃) = ∑ ci di . (3)
be made for each frame. However, when assuming that the i
functions are continuous, we can exploit SH basis functions
This property makes SH integration very efficient, since
in order to approximate Equation 1 and achieve interactive
when interpreting c and d as vectors, the integration can be
frame rates. In the following subsection we briefly review
easily computed on the GPU by using a dot product. Thus,
the main properties of SHs, before describing their integra-
instead of reconstructing and integrating the original func-
tion into our approximation presented in Section 4.
tions, just a dot product is be used. Another important SH
property is the rotational invariance. When a rotation is de-
3.1. Spherical Harmonics sired, it can be directly applied to the SH coefficients. This
is very helpful, when using SH representations for lighting,
Providing a detailed derivation of the concepts behind SHs where the light source may rotate around an object.
would be beyond the scope of this paper. Instead, we refer to
the main properties exploited by our approach and point out
reference to more detailed explanations [Gre03]. 4. Realizing Advanced Material Effects

SH basis functions allow to represent a continuous func- While most SH-based material representations for polygonal
tion given on the sphere. An infinite number of SH basis data are tailored to support one specific material category,
functions Ylm (ϑ , ϕ) of varying order l and degree m can be i. e., highly reflective or scattering, in volume rendering the
weighted with SH coefficients ci j , to represent arbitrarily material category may change just by modifying the trans-
complex functions f (ϑ , ϕ) defined over the sphere: fer function. When for instance a lower opacity is assigned
through the transfer function, stronger scattering can be ex-
pected. Therefore, a specialized SH compression method is
∞ l needed in order to support different material categories.
f (ϑ , ϕ) = ∑ ∑ clmYlm (ϑ , ϕ).
l=0 m=−l In our approach, the SH projection adapts to the material

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

function M at each sample position ~x. In particular, we con-

sider the transparency τ( f (~x)) assigned through the transfer
function, as well as the gradient length | 5 τ( f (~x))|, depict-
ing the heterogeneity of the medium. To illustrate our ap-
proach, we would like to first emphasize the two existing
extreme cases. First, when τ( f (~x)) represents high trans-
parency and | 5 τ( f (~x))| is small, we can assume that the
whole unit sphere centered at ~x affects its illumination. In
contrast, when τ( f (~x)) depicts a rather opaque medium and (a) (b)
| 5 τ( f (~x))| is large, only the front facing hemisphere affects
the illumination of ~x, since ~x can be considered as a part of a Figure 1: The Christmas tree data set, with gradient-based
boundary surface [Lev88]. The other possible combinations lighting (a) and with a color bleeding contribution com-
of the extrema of τ( f (~x)) and | 5 τ( f (~x))| lie somewhere pressed using SH basis functions (b).
between on this scale. So loosely speaking, we consider the
whole contribution of the unit sphere, when no surface-like
structure can be assumed depending on τ( f (~x)), and just the 4.1. Color Bleeding Effects
hemisphere, when a surface-like structure is present.
By using the modified Monte-Carlo integration, the compu-
For a low-frequency spherical function f , it can be shown tation of the color bleeding effects can be done easily by per-
that Monte-Carlo integration over a sufficiently high number forming SH projection. This is a straight forward extension
of samples is sufficient to compute Equation 2 and project f to the diffuse shadowed radiance transfer. However, we omit
into SH space. The strength of this technique is the ability a weighting with the clamped cosine term, depending on the
to generate a uniform distributions of samples over the unit angle between 5τ( f (~x)) and the light vector, since a sim-
sphere. However, in the case where we deal with a surface- ilar behavior is already incorporated in our modified sam-
like structure, we need to sample the front facing hemisphere pling approach described above. Thus, from each voxel, rays
only, since the back facing hemisphere is not assumed to are cast in random directions ωi based on our modified sam-
contribute to the illumination of ~x. Therefore, we have cho- pling scheme. The color and intensity changes, which occur
sen to integrate an adaptive sampling, in order to perform the to the ray while traveling through the medium, are modeled
SH projection of the environmental lighting. By considering by evaluating the front-to-back compositing volume render-
the probability density function, the SH coefficient calcula- ~ i ). This is different to
ing integral in order to obtain Lmi (~x, ω
tion given in Equation 2 can be reformulated that it equals the approach by Ritschel, where only absorption and thus no
the expected value E of the random variable r: chromaticity effects are considered [Rit07]. After the ray-
casting for the current sample has been completed, each of
the accumulated color channel values are used as a weight
for the values Ylm (ϑs , ϕs ) of the current sample, which are
Z Z
f (s)Y j (s) added to the total sum of the SH projection. When substi-
cj = Y j (s) f (s)ds = p(r)ds = E(r). tuting the thus computed Lmi in Equation 1, renderings as
S Ω p(r)
shown in Figure 1 (b) can be generated. When comparing
the outcome with conventional gradient-based shading (see
To evaluate this equation, the probability-density function Figure 1 (a)), it can be seen that a better depth separation can
p(r) needs Rto be known. In our case, it has to fulfill the re- be achieved.
quirement Ω p(r)ds = 1. Since in the spherical case this In this approach, we can use two different ways of incor-
1
area of integration equals 4π, setting p(r) to 4π is a good porating material properties. When evaluating the material
choice. In our case, we have to incorporate both a hemi- function M for each sample on the ray, it directly influences
spherical or spherical region, in order to deal with translu- the bleeding color. However, due to the fact that M is con-
cent and surface-like structure properties. Since these prop- sidered in the SH projection, it cannot be exchanged without
erties are continuous, we introduce a below-surface angle α, performing a new SH projection. Alternatively, M could be
which is proportional to the transparency of a voxel and lies neglected during the ray-marching and instead the transfer
between 0 and π2 . During Monte-Carlo sampling, only rays function color could be taken into account. In this case, ma-
are considered if their angle with 5τ( f (~x)) is less or equal terial functions could be exchanged during rendering inter-
to π2 + α. To fulfill the requirements of the probability den- actively, as long as they have a similar hue.
sity, 2π is used as a base value for the hemisphere, and the
additional ring below the surface given by α is included as
4.2. Scattering Effects
2π + 2π · sin α = 2π(sin α + 1). This term is equal to 2π for
α = 0 and to 4π for α = π2 , which corresponds to the whole With the approach described in the previous subsection, we
sphere. are able to incorporate indirect illumination effects by ex-

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

(a) (b) (c) (d)

Figure 2: Application of the interactive scattering technique to different volume data sets: For the hand data set, we have set Rd
to allow a realistic skin appearance (a). For the engine data set Rd has been set to mimic a porcelain-like appearance (b). For
the teapot an Rd resulting in more scattering has been used (c), while the bonsai shows the occurrence of the bleeding colors
due to the use of brown and green (d).

ploiting the computation of Lmi . However, to further increase tion with a volumetric source distribution [JMLH01]. Thus,
the level of realism for certain materials, scattering effects Equation 4 enables us to compute the scattering contribu-
need to be incorporated. Hao and Varshney could show in ~ i . Notice
tion for each voxel ~x and a given light direction ω
the context of subsurface scattering, that although scattering that the color bleeding contribution is also incorporated in
is in principle a global effect, its influence is rather local due Equation 4, since we use Lmi as derived in the previous sub-
to the exponential falloff of light [HV04]. We have adopted section. The used material parameters Rd have been set as
this insight, and exploit an additional scattering pass. Since specified by Jensen et al. [JMLH01]. Although, the param-
we have already computed and SH projected the Lmi for eters have been originally acquired to simulate subsurface
all ~x as described in the previous subsection, we can ex- scattering, we were able to achieve compelling scattering ef-
ploit this information when computing the scattering con- fects by using these parameters.
tribution. While previous interactive approaches have con-
In our implementation, the neighborhood size |N(~x)| is
strained scattering just to a cone angle [KPH∗ 03, RDRS10],
assumed to be constant within the entire volume. Obviously,
we are thus able to perform the scattering over the whole
the resolution and the average size of structures contained in
unit sphere, or the hemisphere depending on the adapted
the data set should be taken into account in order to specify
Monte-Carlo sampling. When assuming, that multiple scat-
|N(~x)|. However, in our case we could always achieve good
tering dominates, the in-scattering reaching~x from its neigh-
results, when setting |N(~x)| to 10% of the maximum vol-
borhood can be computed as follows:
ume dimension. An alternative approach would be to adapt it
based on the material properties of the current region around
1 1 ~x.
0
Lmi ~ i) =
(~x, ω ∑ ( Rd ) · Lmi (~xi , ω
~ i) · (4)
~x ∈N(~x)
π |N(~
x)| The main benefit of this approach is the fact, that we can
i
0
again use SH projection, in order to efficiently store the Lmi
We have specified this equation in its discrete form to de- as the results of this scattering pass. Thus, as with conven-
pict, that this additional scattering pass is computed based tional SH lighting, the light properties can be changed inter-
on the already present SH coefficient volumes, containing actively. Figure 2 shows the application of Equation 1, when
we substitute Lmi with Lmi0 and use different parameters for
the coefficients for Lmi . Since these volumes are by nature
given in a discrete data structure, the ~xi in Equation 4 re- Rd .
fer to actual voxel positions instead of sample positions.
Equation 4 is inspired by the subsurface scattering equa-
4.3. Local Material Effects
tion introduced by Hao and Varshney [HV04]. Since we
are not interested in subsurface scattering, but true volu- Now, that we are able to compute the environmental illu-
metric scattering, we have neglected the cosine weighting mination contribution, i. e., the in-scattering reaching ~x, we
from the original equation. Furthermore, we have dropped can focus on the realization of the local material effects. As
the Fresnel term, and modified the weighting to comply with seen in Equation 1, this is specified by the material func-
the volumetric scattering neighborhood N(~x). π1 Rd describes ~ i, ω
tion M(~x, ω ~ o ). Since M is also a spherical function, we
the diffuse reflectance obtained by a single dipole approxi- can compress it using standard SH projection techniques.
mation for multiple scattering. This approximation achieves Therefore, we cast rays from the the center of a sphere repre-
a comparable accuracy as using the diffusion approxima- senting the material to the outside, and project the hit values

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

signed to additionally incorporate a complex material func-

tion.

5. Implementation
The concepts described in this paper have been integrated
into a GPU-based volume ray-caster implemented in C++
by using OpenGL in combination with GLSL. The coeffi-
(a) (b) 0 are stored in
cients for the SH projections of Lmi and Lmi
Figure 3: In this example we use BRDFs as material func- 3D coefficient volumes, whereby those representing Lmi can
be discarded, after Lmi0 has been computed, since they are
tions M. The BRDFs have been SH projected by using 1250
samples distributed over the hemisphere. not accessed during rendering. To generate the SH coeffi-
cient volumes, we exploit framebuffer objects in conjunction
with multiple render targets. The actual ray marching is per-
formed in a fragment shader, while rendering slice by slice
on the sphere’s surface. Since this can be done once in a into the 3D coefficient texture. While the GLSL specification
preprocessing for each material added to the system, perfor- describes functionality for generating random numbers, this
mance is not an issue when regarding this precomputation. feature is not implemented yet. Therefore, we have decided
However, to simplify the rendering process, we have decided to generate the random samples for the Monte-Carlo sam-
to stick with the same number of SH bands as in the other pling on the CPU in form of normalized cartesian vectors
SH projections. At the same time, by increasing the number and upload them to the shader as a 2D texture. The same is
of samples, we can improve the quality of the approxima- done for the Ylm (ϑ , ϕ) for each of the sample directions. By
tion. Examples where we have used 1250 samples distributes using the random samples, rays are cast within the shader.
over the hemisphere, are shown in Figure 3. In this case, we After a ray has been terminated, its contribution is multiplied
have used data from the MERL BRDF library [MPBM03] to with the values Ylm (ϑs , ϕs ) obtained from the second 2D tex-
describe the material function M for surface-like structures ture and added to the total integration sum. Thus, we obtain
having a high degree of reflectance. The actual representa- all (L + 1)2 SH-coefficients for the current voxel, which are
tion in our implementation is described in Section 5. distributed to the attached render targets.
We have used the same projection technique to represent To speed up the computation, we have integrated early
light probes, in order to integrate natural lighting conditions ray-termination based on a ray opacity threshold of 95%.
depicted as Li . Thus, we have all necessary factors for com- Furthermore, since the ray-casting is performed directly in
puting Equation 1 and can compute the illumination at ~x by volume space [0...1]3 , a sample ray is terminated if one of
evaluating Equation 3. However, we have the problem, that its position-coordinates reaches 0 or 1.
while SH integration can be done easily for two functions
0 Currently multiple render target functionality is limited
(see Equation 3), dealing with three functions, i. e., M, Lmi
to 8 simultaneous output targets. Therefore, we have chosen
and Li , becomes more complex. While existing techniques
L = 3 as the highest band, yielding 16 coefficients per color.
allow to realize the integration of three functions [Sny06],
We generate these 16 coefficients by rendering into 4 slices
we have decided to use a simpler technique due to perfor-
with 4 channels simultaneously. However, to reduce memory
mance issues. Therefore, we approximate the integration by
requirements, we exploit the fact, that human beings are less
using the following equation, which is not physical correct
sensitive to chromaticity variations as compared to varia-
but led to convincing visual results:

R
L(~x, ωo ) = ~ i, ω
Ω M(~x, ω ~ o )Li (~x, ω ~i ·
~ i )d ω
R 0
L
Ω mi (~x, ~
ω )L
i i (~x, ~
ω i )d ~
ω i.

Intuitively, we simplify the triple product to two succes-

sive double products. First, we integrate over M and Li to ob-
tain the local material effects. The result of this integration
is modulated by the environmental illumination described by
(a) (b) (c)
the second integral. Applications of the introduced rendering
technique are shown in Figure 4. Figure 4 (a) shows the ap- Figure 4: Different effects demonstrated to the engine data
plication of conventional gradient-based shading, while in set: the application of M · Li (a), the environmental light con-
Figure 4 (b) the color bleeding contribution through Lmi has tribution Lmi (b), and the the application of M · Lmi · Li (c).
been incorporated and in Figure 4 (c) BRDFs have been as-

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

tions in luminance [Liv02]. Therefore, we have chosen to use Shadows +Scattering +Bleeding
only a 2-band projection for the chromaticity values while 1283 0.9s / 24 fps 21.3s / 17 fps 58.9s / 13 fps
keeping the 4-band projection for the luminance. Thus, the 2562 × 128 2.6s / 20 fps 70.2s / 14 fps 180s / 5 fps
three color-coefficient sets sum up to only 3 · 4 = 12 coef- 2563 3.4s / 6 fps 149s / 3 fps 388s / 1 fps
ficients altogether, which can be handled using three FBO- Table 1: Timings for SH projection and rendering with data
attached textures during the projection, as well as three 3D sets of different voxel sizes. From left to right: shadows only,
volume textures during rendering. Even when using chro- scattered shadows, scattered shadows with color bleeding.
maticity and luminance, this adds up to 4 + 3 = 7 render
targets, which makes this operation possible in a single pass.
While the SH coefficient volumes are generated on the scattered shadows and color bleeding. During all tests as
GPU, the projection of the material functions and the light well as to generate the renderings shown in the paper, we
probes is less time critical and can be performed on the CPU, have used downsampled coefficient volumes of about 33%
where we generate a partitioned 2D texture representation of the original size, and voxels having an opacity of zero
to be passed to the shader during rendering. When using a have been discarded. Assuming that the radiance transfer is
4-band compression, the resulting texture is partitioned as of low-frequency, the downsampling is a simple and effec-
shown in Figure 5. For each ωo , a combination of texels rep- tive optimization. Furthermore, we have used increasing ray
resents the 16 SH coefficients needed for the ωi . The tex- step size [Rit07]. The windows resolution used to compute
ture consists of three columns, having 4 quadrants of size the given frame rates was 400 × 400.
N × N texels, when using N random samples. The first 4 co-
As it can be seen in Table 1, for data sets of 128 × 128 ×
efficients (bands 0 and 1) are saved in the lower left, the
128 voxels, we achieve interactive frame rates for all com-
next 4 in the lower right (the first 4 coefficients of band 2),
pared SH techniques. This allows to inspect the data set for
the next 4 in the upper left (one coefficient of band 2 and
a given transfer function interactively. Furthermore, for the
3 coefficients of band 3) and the final 4 in the upper right
shadowing technique the SH projection is still fast enough to
(the remaining 4 coefficients of band 3). The x-axis of each
allow to change the transfer function. Obviously, the scatter-
quadrant represents ϑ of ωo , ranging from 0 to π, and the y-
ing pass results in an additional performance penalty, such
axis represents ϕ of ωo , ranging from 0 to 2π. Thus, we have
that the transfer function cannot be changed without waiting
all needed coefficients in a GPU-friendly format and the in-
for the SH projection. However, since the reprojection has
tegration as described in Subsection 4.3, can be performed
to be performed only when the transfer function is changed,
during rendering in a fragment shader by using matrix and
the actual rendering can still explored at 14 −17 frames, i. e.,
vector operations.
camera and light parameters can be changed freely. As de-
scribed in Subsection 4.1, the material function M can only
6. Performance Results be exchanged interactively when it is not considered dur-
ing the SH projection. In fact, in many application areas of
Table 1 shows a comparison of the performance of the SH
medicine or computational fluid dynamics, predefined trans-
projection process for different parameters. The tests were
fer functions are used, which could also involve M. Table 1
performed with three different data sets, having a resolution
also shows, that when additionally incorporating the color
of 128 × 128 × 128, 256 × 256 × 128 and 256 × 256 × 256
bleeding during the shadowing pass, both SH projection as
voxel. All tests have been conducted on a Intel Core2 Quad
well as rendering take significantly longer. This delay is due
CPU Q9450, running at 2.66 GHz, with 4 GB of main RAM
to the fact, that three additional SH coefficient volumes have
and an nVidia GeForce 9800 GTX with 512 MB RAM. In
to be generated and accessed during rendering.
the columns the SH projection times as well as the frame
rates are shown for shadows only, scattered shadows and
7. Conclusions and Future Work
In this paper we have demonstrated how to integrate realistic
light material interaction effects into interactive volume ren-
dering. With the proposed technique we are able to exploit
material functions of reasonable frequency. By using a mod-
ified SH projection, which has been adapted to the needs of
volume rendering, we are able to represent these effects us-
ing a common representation. As we have shown, this also
Figure 5: The BRDF coefficient texture for the aluminum- allows the seamless integration of conventional SH lighting
bronze material. One quadrant holds 4 coefficients for each effects. To our knowledge, this is the first application of non
ωo . The x-axis spans the values ϑ = 0...π and the y-axis the cone-shaped phase functions in the area of interactive vol-
values ϕ = 0...2π. ume rendering. To achieve this, we have explained how to
incorporate the desired effects within the volume rendering

c The Eurographics Association 2010.

 F. Lindemann, T. Ropinski / Advanced Light Material Interaction for DVR

integral, and have described our hybrid CPU/GPU imple- [KPH∗ 03] K NISS J., P REMOZE S., H ANSEN C., S HIRLEY P.,
mentation. We believe that the proposed concepts bring the M C P HERSON A.: A model for volume lighting and modeling.
realism of volume rendered images to a new level. IEEE Transactions on Visualization and Computer Graphics 9, 2
(2003), 150–162. 1, 2, 5
While we were able to produce high-quality imagery at [Lev88] L EVOY M.: Display of surfaces from volume data. IEEE
interactive frame rates, there are still several open issues, Computer Graphics and Applications 8, 3 (1988), 29–37. 4
which could potentially be addressed in the future. The main [LHY09] L JUNG P., H ERNELL F., Y NNERMAN A.: Local ambi-
drawback of SH-based techniques is the restriction to ap- ent occlusion in direct volume rendering. IEEE Transactions on
proximate low-frequency functions only. While this is suffi- Visualization and Computer Graphics 15, 2 (2009). 1, 2
cient for most phase functions, for some materials having [Liv02] L IVINGSTON M.: Vision and Art: The Biology of Seeing.
sharp specular highlights this might not be sufficient and Harry N. Abrams, New York, 2002. 7
Gibbs ringing may occur. Therefore, alternative concepts [Max95] M AX N.: Optical models for direct volume rendering.
should be investigated for these scenarios. Such as the SH IEEE Transactions on Visualization and Computer Graphics 1, 2
(1995), 99–108. 2
approximation method by Zafar et al. [ZARM06], or wavelet
compression techniques, which have been proven useful in [MPBM03] M ATUSIK W., P FISTER H., B RAND M., M C M IL -
LAN L.: A data-driven reflectance model. ACM Transactions on
the area of polygonal rendering [SM06], and might be also
Graphics 22, 3 (July 2003), 759–769. 6
used in DVR. Another issue with the proposed technique is
[PBVG10] PATEL D., B RUCKNER S., V IOLA I., G ROELLER E.:
the availability of appropriate material functions. While sev-
Seismic volume visualization for horizon extraction. In IEEE
eral simplified models exist, having a true volumetric cap- Pacific Visualization (PacificVis 2010) (2010). 2
turing would be of great interest. Finally, since the presented
[RDRS10] ROPINSKI T., D ÖRING C., R EZK -S ALAMA C.: Inter-
technique can be considered as approximative, it should be active volumetric lighting simulating scattering and shadowing.
evaluated with respect to perceptual requirements. In IEEE Pacific Visualization (PacificVis 2010) (2010). 2, 5
[Rit07] R ITSCHEL T.: Fast GPU-based Visibility Computation
for Natural Illumination of Volume Data Sets. In Short Paper
Acknowledgments Proceedings of Eurographics 2007 (2007), pp. 17–20. 2, 4, 7
This work was partly supported by grants from Deutsche [RMSD∗ 08] ROPINSKI T., M EYER -S PRADOW J., D IEPEN -
Forschungsgemeinschaft, SFB 656 MoBil (project Z1). We BROCK S., M ENSMANN J., H INRICHS K. H.: Interactive vol-
would like to thank the reviewers for their constructive com- ume rendering with dynamic ambient occlusion and color bleed-
ing. Computer Graphics Forum (Eurographics 2008) 27, 2
ments. The presented concepts were implemented using the (2008), 567–576. 1, 2
Voreen volume rendering engine (http://www.voreen.
[Sal07] S ALAMA C. R.: GPU-based monte-carlo volume raycast-
org).
ing. In Pacific Graphics (PG) (2007). 1, 2
[SKS02] S LOAN P., K AUTZ J., S NYDER J.: Precomputed radi-
References ance transfer for real-time rendering in dynamic, low-frequency
lighting environments. In ACM SIGGRAPH (2002), pp. 527–536.
[BB09] BANKS D. C., B EASON K.: Decoupling illumination 2
from isosurface generation using 4d light transport. IEEE Trans-
actions on Visualization and Computer Graphics 15, 6 (2009), [SM06] S UN W., M UKHERJEE A.: Generalized wavelet prod-
1595–1602. 2 uct integral for rendering dynamic glossy objects. In ACM SIG-
GRAPH (2006), pp. 955–966. 8
[BG07] B RUCKNER S., G RÖLLER M. E.: Style transfer func-
tions for illustrative volume rendering. Computer Graphics Fo- [Sny06] S NYDER J.: Code generation and factoring for fast eval-
rum 26, 3 (2007), 715–724. 2 uation of low-order spherical harmonic products and squares.
Tech. rep., Microsoft Corporation, 2006. 6
[BGB∗ 06] B EASON K. M., G RANT J., BANKS D. C., F UTCH
B., H USSAINI M. Y.: Pre-computed illumination for isosurfaces. [SPH∗ 09] S CHOTT M., P EGORARO V., H ANSEN C.,
In Visualization and Data Analysis ’06 (SPIE Vol. 6060) (2006), B OULANGER K., B OUATOUCH K.: A directional occlu-
pp. 1–11. 2 sion shading model for interactive direct volume rendering.
In Computer Graphics Forum (Eurographics/IEEE VGTC
[Gre03] G REEN R.: Spherical harmonic lighting: The gritty de- Symposium on Visualization 2009 issue) (2009), pp. 855–862. 1,
tails. In Game Deveopers Conference (2003). 3 2
[HKSB06] H ADWIGER M., K RATZ A., S IGG C., B ÜHLER K.: [WAT92] W ESTIN S. H., A RVO J. R., T ORRANCE K. E.: Pre-
GPU-accelerated deep shadow maps for direct volume render- dicting reflectance functions from complex surfaces. ACM SIG-
ing. In ACM SIGGRAPH/EG Conference on Graphics Hardware GRAPH 26, 2 (1992), 255–264. 2
(GH) (2006), pp. 49–52. 1
[WPHS06] W YMAN C., PARKER S., H ANSEN C., S HIRLEY
[HV04] H AO X., VARSHNEY A.: Real-time rendering of translu- P.: Interactive display of isosurfaces with global illumination.
cent meshes. ACM Trans. Graph. 23, 2 (2004), 120–142. 5 IEEE Transactions on Visualization and Computer Graphics 12,
[JMLH01] J ENSEN H., M ARSCHNER S., L EVOY M., H ANRA - 2 (2006), 186–196. 2
HAN P.: A practical model for subsurface light transport. In [ZARM06] Z AFAR N. B., A KESSON J., ROBLE D., M USETH
ACM SIGGRAPH (2001), pp. 511–518. 5 K.: Scattered spherical harmonic approximation for accelerated
[KH84] K AJIYA J. T., H ERZEN B. P. V.: Ray tracing volume volume rendering. In SIGGRAPH ’06: ACM SIGGRAPH 2006
densities. In ACM SIGGRAPH (1984), pp. 165–174. 2 Sketches (2006), p. 148. 8

c The Eurographics Association 2010.