EUROGRAPHICS 2012 / P. Cignoni, T. Ertl
(Guest Editors)

Volume 31 (2012), Number 2

Real-time Rendering of Dynamic Scenes under All-frequency
Lighting using Integral Spherical Gaussian
Kei Iwasaki1 Wataru Furuya1 Yoshinori Dobashi2 Tomoyuki Nishita3
1 Wakayama University 2 Hokkaido University 3 The University of Tokyo
1 iwasaki@sys.wakayama-u.ac.jp 2 doba@ime.ist.hokudai.ac.jp 3 nis@is.s.u-tokyo.ac.jp

Figure 1: Real-time rendering of fully dynamic scenes with all-frequency lighting and highly specular BRDFs.
Abstract
We propose an efficient rendering method for dynamic scenes under all-frequency environmental lighting. To
render the surfaces of objects illuminated by distant environmental lighting, the triple product of the lighting, thevisibility function and the BRDF is integrated at each shading point on the surfaces. Our method represents theenvironmental lighting and the BRDF with a linear combination of spherical Gaussians, replacing the integral of
the triple product with the sum of the integrals of spherical Gaussians over the visible region of the hemisphere.
We propose a new form of spherical Gaussian, the integral spherical Gaussian, that enables the fast and accurate
integration of spherical Gaussians with various sharpness over the visible region on the hemisphere. The integral
spherical Gaussian simplifies the integration to a sum of four pre-integrated values, which are easily evaluatedon-the-fly. With a combination of a set of spheres to approximate object geometries and the integral spherical
Gaussian, our method can render object surfaces very efficiently. Our GPU implementation demonstrates realtime rendering of dynamic scenes with dynamic viewpoints, lighting, and BRDFs.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism —Color, shading, shadowing, and texture

1. Introduction
Real-time rendering of fully dynamic scenes illuminated byenvironmental lighting, where the viewpoint, lighting, and
BRDFs can be changed and objects can move and deform
at run-time, is beneficial in many applications such as lighting/material design, animation, and games. For photorealistic rendering under environmental lighting, the triple product of the environmental lighting, the BRDF, and the visibility function is integrated. Integrating the triple productis computationally expensive, especially for dynamic scenes
c 2011 The Author(s)Computer Graphics Forum c 2011 The Eurographics Association and Blackwell Publishing Ltd. Published by Blackwell Publishing, 9600 Garsington Road, Oxford OX4 2DQ,
UK and 350 Main Street, Malden, MA 02148, USA.

with highly-specular BRDFs, since the visibility function
is high-frequency and BRDFs have sharp peaks; moreover,both of these vary at run-time due to the dynamic viewpoint,changing BRDFs and objects. To achieve real-time performance, previous methods tackled this problem by restricting
the methods to static scenes [SKS02,NRH03,TS06,GKD07,
WRG∗ 09], low-frequency lighting [RWS∗ 06, NKF09], or
moderately specular reflectance [ADM∗ 08, RGK∗ 08]. Realtime rendering of fully dynamic scenes under all-frequency
lighting still remains a challenging task.

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

To address this, we propose an efficient rendering methodfor dynamic scenes illuminated by all-frequency environmental lighting. Our method uses spherical Gaussians(SGs) that are identical to isotropic Gaussians [GKD07].
SGs are also a type of spherical radial basis functions(SRBFs) [TS06]. For diffuse BRDFs, our method represents
the environmental lighting with a small number of SGs. For
specular BRDFs, our method uses prefiltered environment
maps with SGs of different sharpness, which makes it possible to deal with dynamic, highly-specular BRDFs. Both calculations require the integral of the product of the SGs andthe visibility, which makes the rendering difficult and limits
the previous method [WRG∗ 09] to static scenes. We propose a novel calculation method, integral spherical Gaussian (ISG), for integrating SGs with arbitrary sharpness overthe visible region on the hemisphere at each shading point.
SGs with various sharpness are pre-integrated over different
integration intervals of the elevation and azimuthal angles.
The integrated values are fitted to a simple function of theelevation angle, azimuthal angle, and the sharpness of theSG, enabling direct evaluation on the GPU. Therefore, theintegral spherical Gaussians do not need huge precomputed
data nor precomputational time, which are required for previous methods [NRH03, TS06]. The integral of the SG overthe visible region on the hemisphere is efficiently calculatedusing the pre-integrated SG and a set of spheres approximating the object geometries. Furthermore, our method utilizes
a hierarchical representation of spheres for scalability andlevel of detail. All the computations are performed in parallel on the GPU.
Our method contributes the following:• Real-time rendering of fully dynamic scenes under allfrequency lighting with highly specular BRDFs.
• An efficient calculation method to integrate various SGs
over the visible region on the hemisphere using the integral spherical Gaussian, which is easy to implement andtherefore amenable to GPUs.
• No huge amount of precomputed data for the visibility
function nor the BRDFs is needed, which allows the user
to change the object geometries and BRDFs at run-time.
2. Previous Work
Precomputation based methods: Sloan et al. proposed a
Precomputed Radiance Transfer (PRT) method to render
static scenes illuminated by low-frequency lighting environments [SKS02]. Several methods have been proposed torender all-frequency shadows using Haar wavelets [NRH03,
NRH04, WTL06] or spherical radial basis functions [TS06].
To render highly specular materials, Green et al. proposed areal-time rendering method by approximating BRDFs withisotropic Gaussian functions [GKD07]. Sun et al. developed an interactive relighting method for dynamic BRDFs
with indirect illumination [SZC∗ 07]. Wang et al. proposed a real-time rendering method for dynamic, spatially-varying BRDFs [WRG∗ 09]. Lam et al. proposed a multiscale spherical radial basis function that represents allfrequency lighting, and achieved real-time rendering with its
GPU-implementation [LHLW10]. Although these methods
can render scenes illuminated by environmental lighting inreal-time, they are limited to static scenes.
Kautz et al. proposed an interactive rendering method for
self-shadows of dynamic objects using hemispherical rasterization [KLA04]. Sloan et al. proposed a local deformable
PRT method using zonal harmonics [SLS05]. Ren et al. proposed a real-time rendering method for dynamic scenes under low-frequency lighting [RWS∗ 06]. Sloan et al. extended
this method to calculate the indirect illumination for dynamic scenes [SGNS07]. These two methods approximate
the object shape with a set of spheres and calculate the occlusion due to the set in the Spherical Harmonics (SH) domain using a SH exponentiation technique. Nowrouzezahrai
et al. presented a real-time shadowing method with arbitrary BRDFs for animated models [NKF09]. These methods,
however, cannot be applied to an all-frequency lighting environment.
Zhou et al. proposed shadow fields that record shadowing
effects due to a rigid object around the object and achieved
real-time rendering of dynamic scenes for low-frequency
lighting and interactive rendering for all-frequency lighting [ZHL∗ 05]. Tamura et al. extended precomputed shadow
fields to record shadowing effects adaptively and achieved
real-time rendering on the GPU [TJCN06]. Sun et al. proposed a generalized wavelet product integral method to render dynamic scenes of rigid objects under all-frequency
lighting [SM06]. Xu et al. proposed a real-time rendering
method of dynamic scenes under all-frequency lighting using spherical piecewise constant functions [XJF∗ 08]. Thesemethods, however, cannot deal with deformable objects.
Shadow map based methods: Krivanek et al. proposed
filtered importance sampling and rendered dynamic scenes
using the approximate visibility calculated from shadow
maps [KC08]. This method, however, utilizes a small number of shadow maps, which may cause incorrect visibilities. Annen et al. proposed a real-time rendering methodfor dynamic scenes under all-frequency lighting [ADM∗ 08].
This method approximates the environmental lighting witha set of area light sources, and render shadows from theselight sources using a convolution shadow map method. Thismethod, however, does not integrate the product of the
BRDF and the incident radiance. Ritschel et al. proposed
imperfect shadow maps to calculate the indirect illumination for dynamic scenes [RGK∗ 08]. They also demonstrated
that these imperfect shadow maps can be applied to environment map lighting. Hollander et al. proposed the concept of ’ManyLoDs’ to accelerate the many-view rasterization process [HREB11]. These methods determine the light
positions to render shadow maps by using importance sampling of the environment maps, and ignore the importancec 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussianof the material BRDFs. Therefore, these methods are unsuitable for rendering highly specular materials.

3. Overview
Our method focuses on fast rendering of dynamic scenes
with direct illumination. The outgoing radiance L(x, ω o ) at
point x in the outgoing direction ω o is calculated by the following equation.
Z

L(x, ω o ) =ωi )V (x, ω i )ρ(x,ωωi ,ωωo )max(0,n(x)·ω
ωi )dω
ωi ,
L(ω

S2

(1)
where S2 denotes the unit sphere in R3 , ω i is the incidentωi ) is the distant lighting represented by the endirection, L(ωvironment maps, V (x, ω i ) is the visibility function at x, ρ isthe BRDF, and n(x) is the normal at x. To simplify the notation, we omit x in the following. Our method represents
ωi , ω o ) = kd + ks ρs (ωωi , ω o ) where kd is
the BRDF ρ with ρ(ωωi , ω o ) is a specular term. By suba diffuse term and ks ρs (ωωo ) can be calculated from thestituting this in Eq. (1), L(ωsum of the diffuse component Ld and the specular compoωo ) as follows:
nent Ls (ω

Ld =

Z

S2

ωi )V (ωωi ) max(0, ω i · n)dω
ωi ,
L(ω

ωo ) =
Ls (ω

Z

S2

ωi )V (ωωi )ρs (ωωi , ωo ) max(0, ωi ·n)dωωi .(3)L(ω(2)ωi ) and the coTo calculate Ld , our method approximates L(ωsine term max(0, n · ω i ) with the sum of spherical Gaussians.
ωi ; ξ , η, µ) is represented by
A spherical Gaussian (SG) G(ωthe following equation:(4)

where the unit vector ξ is the lobe axis, η is the lobe sharpness, and µ is the lobe amplitude that consists of RGB components. Our method represents the environmental lighting Lωi ) ≈ ∑K
ωi ; ξ k , ηk , µk ) where K
as the sum of SGs; L(ωk=1 G(ωis the number of SG lobes. The cosine term is also approxωi ; n, ηc , µc ). By substituting
imated by a single SG as G(ωthese terms into Eq. (2), Ld is calculated by the following
equation.
K Z

Ld ≈ ∑

2k=1 S

ωi ; ξ k , ηk , µk )G(ωωi ; n, ηc , µc )V (ωωi )dω
ωi . (5)G(ωThe product of two SGs can be represented with a single SG:
ωi ; ξ, η, µ) = G(ωωi ; ξk , ηk , µk )G(ωωi ; n, ηc , µc ),
G(ω

(6)

where the lobe sharpness η is ||ηk ξ k + ηc n||, the lobe axis
ξ is (ηk ξ k + ηc n)/η, and the amplitude µ is µk µc eη−ηk −ηc .
The diffuse component Ld is calculated by integrating theωi ) and the spherical
product of the visibility function V (ωωi ; ξ , η, µ). For fully dynamic scenes where theGaussian G(ωc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

To calculate the visibility function and integrate the product at run-time, our method approximates the object geometry by a set of spheres. This makes it possible to calculate
the solid angles of occluded regions for the SGs efficiently.
To calculate the integral of the product of the visibility function and the SG, our method extends the summed area table(SAT) [Cro84] (or integral image [VJ01]) over the 2D image
domain to the SG over the 2D domain on a unit sphere. We
call this integral spherical Gaussian (ISG) (Section 4). The
integral of the product of the visibility function and arbitrary
SGs is efficiently calculated using the ISG.
To calculate the specular component Ls , our method represents the BRDF ρ and the cosine term as the sum of SGs.
ωi , ω o ) can be calculated fromAs [WRG∗ 09], the BRDF ρ(ωωi , ω o ) ≈ ∑Jj=1 G(ωωi ; ξ j , η j , µ j ),the sum of the SGs as ρ(ωwhere J is the number of SGs for ρ, and the lobe axis ξ j depends on the outgoing direction ω o . By substituting this intoωo ) can be rewritten as:
Eq. (3), Ls (ωJ Z

ωo )≈ ∑Ls (ωωi )G(ωωi ; ξ j ,η j ,µ j )G(ωωi ; n,ηc ,µc )V (ωωi )dω
ωi .(7)
L(ω

2

j=1 S

ωo ) = kd Ld + ks Ls (ω
ωo ),
L(ω

ωi ·ξξ−1)
ωi ; ξ , η, µ) = µeη(ωG(ω,lighting and viewpoint change, and objects move and deωi ) and G(ωωi ; ξ , η, µ) change and cannot be preform, V (ωcomputed.

Since the product of two SGs is also represented with a sinωi ) = G(ωωi ; ξ ′j , η′j , µ′j ), Ls (ωωo ) can be rewritten
gle SG G j (ωas:J Zωo ) ≈ ∑Ls (ωωi )G j (ωωi )V (ωωi )dω
ωi .
L(ω

2
j=1 S

(8)

Our method approximates the integral of the triple productas:J Zωo ) ≈ ∑
Ls (ω

2
j=1 S

ωi )V (ωωi )dω
ωi
2 G j (ωωi )G j (ωωi )dω
ωi S R
L(ω
.(9)
ω
ω
G
(ω
)dω
2
j
i
i
S
R

ωi )V (ωωi )dω
ωi can be calculated in the
The numerator S2 G j (ωRωi )dω
ωi can besame way as Ld and the denominator S2 G j (ω
R

′

−η j
). The integral of thecalculated analytically as 2π
η′ (1 − ejproduct of L and G is calculated as a 3D function ΓL (ξξ, η):ZS2ωi )G(ωωi ; ξ , η, µ)dωωi = µΓL (ξξ, η).
L(ω

(10)

Our method precomputes ΓL (ξξ, η) and stores the data asprefiltered environment maps for various lobe sharpness
η as [WRG∗ 09]. Therefore, our method can change the
BRDFs (i.e. the SGs) at run-time. Since our method integrates the product of the BRDF approximated by the SGs
and the lighting, our method can represent highly specularmaterials.
In summary, efficient computation of the integral of the
product of the SGs and the visibility function is the key element for both diffuse and specular components to achieve
real-time rendering of fully dynamic scenes. We propose a

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

blocker spheres

patch(a) patchoccludedpatch(b) occluded patch

!"#

!$#

!%#

!&#

!'#

!(#

Figure 3: 3D plots of (a)Ŝ(θ, φ, 100), (b) f (θ, φ, 100) and(c) difference between (a) and (b) , and (d)Ŝ(θ, φ, 1000),(e) f (θ, φ, 1000), and (f) difference between (d) and (e). Relative squared error between (a) and (d) is 0.061%, and
0.0013% between (b) and (e). Maximum error between (a)and (b) is 0.0966 and that between (d) and (e) is 0.0251.

blocker sphere(c) projected sphere(d) occluded angle

Figure 2: Local coordinate system for integral spherical
Gaussian and patch.

for various sharpness η (from 1 to 10000) using the leastsquares method. The polynomials of g and h are describedin the Appendix. Consequently, ISG S(θ, φ, η) is calculated−ηby 2π) f (θ, φ, η).
η (1 − efast calculation method of the integral of the product in thenext Section.

Our method calculates the integral of the product of the
SGs and the visibility function by subtracting the integral of
the SGs over the occluded region Ω̄ on the hemisphere fromthe ISG S(π, π, η) as follows:

4. Integral Spherical GaussianZ

S2

We now focus on the calculation of the integral of the product of the visibility function and the SGs. Our method sets
the local coordinate system so that the axis ξ is set to the zaxis as shown in Fig. 2, and the direction ω is parameterizedby two angles θ ∈ [0, π] and φ ∈ [0, π] as:ω (θ, φ) = (sin θ cos φ, cos θ, sin θ sin φ),(11)

where θ is the angle between ω and the y axis, and φ is theangle between ω and the xy plane. We define the integral
spherical Gaussian S(θ, φ, η) as
S(θ, φ, η) =Z θZ φ00

ω(θ′ , φ′ ); z, η, 1) sin θ′ dθ′ dφ′ . (12)G(ωWe also define the normalized integral spherical GaussianŜ(θ, φ, η) as:Ŝ(θ, φ, η) = S(θ, φ, η)/S(π, π, η) =ωi ; z, η, 1)V (ωωi )dω
ωi =S(π, π, η) −G(ωS(θ, φ, η).(13)(2π/η)(1 − e−η )Fig. 3 shows normalized ISGs for different values of lobe
sharpness η. Our method approximates the normalized ISG
with the following function f (θ, φ, η) consisting of two sigmoid functions:11f (θ, φ, η) =, (14)
1 + e−g(η)(θ−π/2) 1 + e−h(η)(φ−π/2)
where g and h are polynomial functions of lobe sharpness
η. Our method fits g and h to fourth degree polynomialsZωi ; z, η, 1)dω
ωi .
G(ω

Ω̄

To simplify the explanation, we first describe the integralof the SG over the occluded region Ω̄ bounded by four directions ω (θ0 , φ0 ), ω (θ0 , φ1 ), ω (θ1 , φ0 ), and ω (θ1 , φ1 ) as
shown in Fig. 2(a). The area bounded by these four cornerdirections is referred to as the patch P(θ0 , θ1 , φ0 , φ1 ).
Then the integral of the SG over the patch P(θ0 , θ1 , φ0 , φ1 )is efficiently calculated by:Z θ1 Z φ1θ0

φ0ω(θ′ , φ′ ); z, η, 1) sin θ′ dθ′ dφ′
G(ω= S(θ1 , φ1 , η) − S(θ0 , φ1 , η) − S(θ1 , φ0 , η) + S(θ0 , φ0 , η).
4.1. Calculation of Occluded Patches
With our method, the integral of the SG over the occluded
regions due to multiple overlapping blocker spheres can becalculated using ISG. To compute the integral efficiently, ourmethod discretizes the surface on the hemisphere into small
patches adaptively based on the lobe sharpness η. Firstly, we
1
calculate the limiting angle γ that satisfies 1+e−g(η)(γ−π/2)> ε,where the threshold ε is specified by the user. In our implementation, we set ε to 0.01 and obtain good results. Then
our method uniformly subdivides the surface on the hemisphere around the z axis within the limiting angle γ as shownin Fig. 2(b).
For each blocker sphere, the patches overlapped by thec 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

projected blocker sphere are calculated and considered as
occluded. The occluded patches are determined by using thefour angles of the patches. Let cl and rl be the center and theradius of l-th blocker sphere in the local coordinate system(see Fig. 2(c)). We calculate the angles θl and φl for cl , andthe angle δl that corresponds to the radius of the projected
blocker sphere as δl = sin−1 (rl /||cl ||). Then the occluded
region in the θ direction is calculated as θl −δl ≤ θ ≤ θl +δl .
For the m-th discretized angle θm as shown in Fig. 2(d), theoccluded region in the φ direction is√φl − δm ≤ φ ≤ φl + δm ,Figure 4: Comparison of shadows using ISG and SSDF. To
compute ISG and SSDF, a sphere set that approximates thecharacter geometry is used. To render (c) using ray tracing,the character geometry itself is used to compute the visibility.
The environment map is approximated by 9 SGs.

δ2 −(θ −θ )2mll. A patchwhere δm is calculated as δm =sin θmwhose four angles overlap these regions is determined as occluded.

Finally, the ISGs of the occluded patches are calculated−ηand subtracted from 2π), resulting in the integral ofη (1 − ethe product of the visibility function and the SG. For narrow SG lobes, since the limiting angle is very small, thediscretized patches are very small and this provides good
accuracy for discretization. For broader SG lobes, the SGchanges smoothly and large patches are sufficient to approximate occluded regions.

4.2. Discussion
To integrate the product of the SG and the visibility, the previous method [WRG∗ 09] proposed a spherical signed distance function (SSDF) V d (ξξ). The SSDF V d (ξξ) returns thesigned angular distance between the lobe axis ξ and thenearest boundary direction. The visible region is approximated as the region where the angle φ in Fig. 2 is larger thanπ/2 − V d (ξξ). Then the integral of the product of the visibility function and the SG is calculated by integrating the
SG over the visible region φ ≥ π/2 − V d (ξξ). That is, theintegral calculation using the SSDF can be viewed as a special case of our ISG as S(π, π/2 − V d (ξξ), η). Although the
SSDF representation can render realistic shadows in the casewhere the occluded region widely covers the y axis (i.e. θ in
Fig. 2 is nearly π), the SSDF can overestimate the occluded
region by assuming that θ is π. Fig. 4 shows a comparison
between our method using an ISG, the SSDF, and a ray tracing method. As shown in Fig. 4, shadow regions computed
using the SSDF are large compared to the reference solution,and is overall dark due to overestimation of the invisible region. The ISG can calculate shadows more accurately, andcan produce plausible shadows matching the reference solution. Note that computing SSDF is time-consuming even foran object approximated by spheres, especially for the SGs
whose lobe axis direction is invisible (occluded). In this case,
we need to shoot many rays to search for the nearest visible
direction. This process is too expensive to compute at runtime for dynamic scenes.
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

5. Implementation Details
Our method calculates the outgoing radiance kd Ld +ωo ) by approximating the object geometries in theks Ls (ωscene with a set of spheres. The computational cost of the
diffuse component Ld for each point is O(KNs ) and that of
the specular component Ls is O(JNs ), where K is the numberof SGs to represent the lighting L, J is the number of SGs torepresent the BRDF ρ, and Ns is the number of spheres. As
described in [WRG∗ 09], K and J are small enough to produce plausible lighting (e.g. K ≤ 10 and J ≤ 7 for typicalscenes). Therefore the performance of our method mainly
depends on the number of spheres.
To render dynamic scenes efficiently, our method represents the spheres using a bounding volume hierarchy (BVH).
A set of spheres that approximates each object geometry is
prepared using [WZS∗ 06] or [BO04]. These spheres are considered as leaf nodes and parent nodes are constructed thatbound spheres of child nodes. In our implementation, a binary tree structure is used to represent the BVH. The BVH
is constructed using [LGS∗ 09] whose complexity is O(Ns )and therefore this method scales well for many spheres.
To compute the integral of the product of the SG with lobe
axis ξ and the visibility, our method traverses the binary treeof spheres. As described in the previous section, our method
only considers the spheres that intersect a cone defined by
the limiting angle γ centered on the lobe axis ξ as shown in
Fig. 5. From the root node, our method calculates the cut of
the binary tree as follows. For each node, our method first
determines whether the bounding sphere corresponds to thenode intersecting the cone or not. If the sphere does not intersect the cone, the node and its descendant nodes are culled.
If the sphere intersects the cone, the projected sphere in thelocal coordinate system is calculated. If the projected sphereis larger than the threshold specified by the user, our method
traverses the child nodes. Otherwise, the patches that overlap the projected sphere are calculated and considered as occluded.
Although sphere approximation can accelerate the rendering, it is difficult to approximate thin shells (e.g. clothes)using spheres. To handle thin shells, our method projects triangle meshes instead of spheres, calculating the occluded
patches by rasterizing the projected triangle meshes similar
to [KLA04].

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

cone(a) intersection test
between sphere and cone

culledspheres(b) BVH of spheres

Figure 5: Intersection tests between spheres and cones using BVH. Dashed (solid) spheres in (a) represent inner (leaf)
nodes of BVH in (b). Shaded spheres in (a) intersect thecone. Since red dashed sphere does not intersect the cone,the descendant nodes are culled.

Our method utilizes two methods to calculate the outgoingradiance, the per-vertex shadowing method and the per-pixel
shadowing method. For the per-vertex shadowing method,
the integral of the product of the SG and the visibility is
calculated at each vertex. Then the integral is interpolated
and shading is performed at each pixel using a fragment program. For the per-pixel shadowing method, our method calculates the integral at each pixel using CUDA. To do this,our method first renders the scene into the framebuffer object to store the position and normal at each pixel. Then thebuffer is transferred to CUDA and the integration of the SG
is performed at each pixel.
6. Results
We have implemented our rendering algorithm on a standard
PC with an Intel Corei7 Extreme 965 CPU and a GeForceGTX 580 GPU. The statistics of examples in this paper are
listed in Table. 1. The SGs that approximate the environmenωi ) are calculated using [TS06]. Our method
tal lighting L(ωrepresents isotropic BRDFs such as the Blinn-Phong BRDF
model and the Cook-Torrance BRDF model with 1 SG, andthe anisotropic BRDF (Ashikhmin-Shirely BRDF) with 7
SGs. The precomputation time and the data size for ΓL are
10 min and 6MB, respectively.
Fig. 6 shows images of a Buddha model; (a), (b) withshadows and (c) without shadows. As shown in Fig. 6(a), ourmethod can render self-shadows onto the pedestal. Fig. 6(b)shows the shadows due to the Buddha model with specular reflections on the floor. Fig. 7 shows still images of achess scene that consists of 32 chess pieces and a board. The
BRDFs of the chess pieces can be changed from diffuse tohighly specular interactively at run-time. Our method discretizes the hemisphere into 100 to 400 patches in these examples. Fig. 8 shows still images of six deformable characters. Since our method does not require precomputed visibility data, our method can render deformable objects inreal-time. Fig. 9 shows still images of rigid body simulations. The interactions between 8 bunny models and 8 chesspieces are simulated using NVIDIA PhysX. Fig. 10 shows
still images of simulations of interactions between a clothand 4 bunny models. As shown in Fig. 10, our method candeal with thin shells (e.g. the cloth). The shadows due to thecloth are rendered by projecting triangle meshes of the clothand calculating the overlapped occluded patches.
Fig. 11 shows images rendered using (a) our method, (b)a ray tracing method, and (c), (d) shadow map methods.
The BRDF of the bunny model is the Blinn-Phong BRDF
whose glossiness parameter is 1000. As shown in Fig. 11,
the image rendered using our method matches that renderedusing the ray tracing method. Fig. 11(c) is rendered using
30 shadow maps which is equal to the number of shadow
maps used in [ADM∗ 08]. The lighting directions of the
shadow maps are calculated using importance sampling of
the environment maps. As shown in Figs. 11, the specular reflections of a highly specular BRDF cannot be rendered using only 30 shadow maps. Fig. 11(d) is renderedusing 1024 shadow maps which is equal to the number of
shadow maps in [RGK∗ 08], but the specular reflection isstill different from Fig. 11(b). The computational times forFig. 11 are (a) 14.5ms, (b) 1044sec, (c) 23.0ms, and (d)262ms, respectively. Similar to Fig. 11(c), the convolution
shadow map method [ADM∗ 08] does not integrate the product of the BRDF and the incident radiance. Therefore, thismethod fails to capture highly specular reflections. As described in [RGK∗ 08], the imperfect shadow map method canrender the scene (Fig. 11 in [RGK∗ 08]) in 25ms, while therendering time using 1024 classic shadow maps is 405ms.
Therefore, the estimated time to render Fig. 11(d) is about
16.2ms using imperfect shadow maps, comparable to our
method. However, our method can render specular reflections better as shown in Fig. 11 (especially as highlighted in
the red frame).
Fig. 12 shows a comparison between our method and a
ray tracing method. Figs. 12(a) and (d) are rendered by using our method, Figs. 12(b) and (e) are rendered by using a
ray tracing method, and Figs. 12(c) and (f) are difference images scaled by 4. Figs. 12 (b) and (d) are rendered by tracing
6 × 64 × 64 rays per pixel to sample the incident radiance.
The RMS errors in (c) and (f) are 3.2% and 4.6%, respectively. As shown in Fig. 12, our result matches the reference
solution very well and any differences are subtle.
Although our method can render convincing results, italso has several limitations. Firstly, our method handles only
direct illumination. This may be addressed by combining[LWDB10] with our method, but this remains as a future
work. For almost ideally-specular surfaces, the circle boundaries due to the sphere set approximation may be visible.
This may be relaxed by increasing the number of spheres
approximating the geometries at the cost of decreasing therendering performance. However, as shown in Fig. 12, ourc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

Table 1: Rendering performance. #vertices and #spheres
represent the number of vertices and spheres, respectively.
fps represents the rendering frame rate for fully dynamic
scenes (dynamic lighting, objects, BRDFs, viewpoint). The
image resolutions are 640 × 480 except for Fig. 11(400 ×
400). Per-vertex shadowing is used except for Fig. 12 inwhich per-pixel shadowing is performed.
FigureFig. 6
Fig. 7
Fig. 8Fig. 9Fig. 10
Fig. 11
Fig. 12# vertices
57K
101K
27K
66K
25K
18K
18K

# spheres
46
306
192
192
216
53
53

fps
40fps
7-9fps33-47fps
13-60fps
2.5-3fps
69fps
10fps

Figure 7: Chess scene.

Figure 8: Character scene.

[BO04] B RADSHAW G., O’S ULLIVAN C.: Adaptive medial-axis
approximation for sphere-tree construction. ACM Trans. Graph.
23, 1 (2004), 1–26. 5[Cro84] C ROW F. C.: Summed-area tables for texture mapping.
SIGGRAPH Comput. Graph. 18 (1984), 207–212. 3

(a)
(b)
(c)
Figure 6: Buddha scene (a) viewed from front, (b) viewed
from back with shadows, and (c) without shadows.
method can render plausible hard shadows comparable to thereference solution.
7. Conclusions and Future Work
We have proposed a real-time rendering method for fully dynamic scenes with highly specular BRDFs illuminated by
all-frequency lighting. The key contribution of our method isan efficient method for calculating the integral of the product of various SGs and the visibility function. The integralof the product of arbitrary SGs and visibility is efficiently
calculated using the sum of the integral spherical Gaussians
which can be easily evaluated on-the-fly for various SG lobe
sharpness. The integral of the triple product of the lighting,the visibility function and the BRDFs can be efficiently calculated using the integral spherical Gaussian.
In future work, we would like to extend our method to
calculate indirect illumination.
Acknowledgments
This work was supported in part by Grant-in-Aid for Scientific Research(C) No. 20500082, and by the Japan Prize
Foundation.
References[ADM∗ 08] A NNEN T., D ONG Z., M ERTENS T., B EKAERT P.,
S EIDEL H.-P., K AUTZ J.: Real-time all-frequency shadows indynamic scenes. ACM Trans. Graph. 27, 3 (2008), 34:1–34:8. 1,2, 6
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

[GKD07] G REEN P., K AUTZ J., D URAND F.: Efficient reflectance and visibility approximations for environment map rendering. Computer Graphics Forum 26, 3 (2007), 495–502. 1,
2[HREB11] H OLLANDER M., R ITSCHEL T., E ISEMANN E.,
B OUBEKEUR T.: Manylods: Parallel many-view level-of-detailselection for real-time global illumination. Computer Graphics
Forum 30 (2011), 1233–1240. 2[KC08] K RIVANEK J., C OBERT M.: Real-time shading with filtered importance sampling. Computer Graphics Forum 27, 4(2008), 1147–1154. 2[KLA04] K AUTZ J., L EHTINEN J., A ILA T.: Hemispherical rasterization for self-shadowing of dynamic objects. In Eurographics Symposium on Rendering (2004), pp. 179–184. 2, 5
[LGS∗ 09] L AUTERBACH C., G ARLAND M., S ENGUPTA S.,
L UEBKE D., M ANOCHA D.: Fast bvh construction on gpus.
Computer Graphics Forum 28, 2 (2009), 375–384. 5[LHLW10] L AM P.-M., H O T.-Y., L EUNG C.-S., W ONG T.T.: All-frequency lighting with multiscale spherical radial basis functions. IEEE Transactions on Visualization and Computer
Graphics 16, 1 (2010), 43–56. 2[LWDB10] L AURIJSSEN J., WANG R., D UTR Ã L’ P., B ROWN
B.: Fast estimation and rendering of indirect highlights. Computer Graphics Forum 29, 4 (2010), 1305–1313. 6[NKF09] N OWROUZEZAHRAI D., K ALOGERAKIS E., F IUME
E.: Shadowing dynamic scenes with arbitrary brdfs. Computer
Graphics Forum 28, 1 (2009), 249–258. 1, 2
[NRH03] N G R., R AMAMOORTHI R., H ANRAHAN P.: Allfrequency shadows using non-linear wavelet lighting approximation. ACM Trans. Graph. 22, 3 (2003), 376–381. 1, 2
[NRH04] N G R., R AMAMOORTHI R., H ANRAHAN P.: Triple
product wavelet integrals for all-frequency relighting. ACM
Trans. Graph. 23, 3 (2004), 477–487. 2[RGK∗ 08] R ITSCHEL T., G ROSCH T., K IM M. H., S EIDEL H.P., DACHSBACHER C., K AUTZ J.: Imperfect shadow maps

K. Iwasaki, & W. Furuya & Y. Dobashi & T. Nishita / Real-time Rendering of Dynamic Scenes using Integral Spherical Gaussian

Figure 11: Comparison of (a) our method (14.3ms), (b) ray
tracing method (1044sec), (c) 30 shadow maps (23.0ms),
and 1024 shadow maps (262ms). Shadow map resolutionsof (c) and (d) are 256 × 256.

!"#

!$#

!%#

!&#

!'#

!(#

Figure 9: Falling rigid objects scene.

Figure 10: Cloth scene.
for efficient computation of indirect illumination. ACM Trans.
Graph. 27, 5 (2008), 129:1–129:8. 1, 2, 6
[RWS∗ 06] R EN Z., WANG R., S NYDER J., Z HOU K., L IU X.,
S UN B., S LOAN P.-P., BAO H., P ENG Q., G UO B.: Real-time
soft shadows in dynamic scenes using spherical harmonic exponentiation. ACM Trans. Graph. 25, 3 (2006), 977–986. 1, 2
[SGNS07] S LOANP.
P.,
G OVINDARAJUN.
K.,N OWROUZEZAHRAI D., S NYDER J.:
Image-based proxyaccumulation for real-time soft global illumination. In PacificGraphics 2007 (2007), pp. 97–105. 2[SKS02] S LOAN P.-P., K AUTZ J., S NYDER J.: Precomputed
radiance transfer for real-time rendering in dynamic, lowfrequency lighting environments. ACM Trans. Graph. 21, 3(2002), 527–536. 1, 2

Figure 12: Comparison of our method and ray tracing
method. BRDFs of the bunny and floor for (a) and (b) areboth Phong BRDF model (the glossiness parameter is 1000
and ks = 0.5). For (d) and (e), ks of the bunny is set to 0.01.
(c) and (f) are difference images scaled by 4.
Conference on Computer Vision and Pattern Recognition, 2001
(2001), pp. 511–518. 3[WRG∗ 09] WANG J., R EN P., G ONG M., S NYDER J., G UO
B.: All-frequency rendering of dynamic, spatially-varying reflectance. ACM Trans. Graph. 28, 5 (2009), 133:1–133:10. 1, 2,
3, 5[WTL06] WANG R., T RAN J., L UEBKE D.: All-frequency relighting of glossy objects. ACM Trans. Graph. 25, 2 (2006),
293–318. 2[WZS∗ 06] WANG R., Z HOU K., S NYDER J., L IU X., BAO H.,
P ENG Q., G UO B.: Variational sphere set approximation forsolid objects. Vis. Comput. 22, 9 (2006), 612–621. 5[SLS05] S LOAN P.-P., L UNA B., S NYDER J.: Local, deformableprecomputed radiance transfer. ACM Trans. Graph. 24, 3 (2005),
1216–1224. 2[XJF∗ 08] X U K., J IA Y.-T., F U H., H U S., TAI C.-L.: Spherical piecewise constant basis functions for all-frequency precomputed radiance transfer. IEEE Transactions on Visualization andComputer Graphics 14, 2 (2008), 454–467. 2[SM06] S UN W., M UKHERJEE A.: Generalized wavelet product integral for rendering dynamic glossy objects. ACM Trans.
Graph. 25, 3 (2006), 955–966. 2[ZHL∗ 05] Z HOU K., H U Y., L IN S., G UO B., S HUM H.-Y.: Precomputed shadow fields for dynamic scenes. ACM Trans. Graph.
24, 3 (2005), 1196–1201. 2[SZC∗ 07] S UN X., Z HOU K., C HEN Y., L IN S., S HI J., G UO B.:
Interactive relighting with dynamic brdfs. ACM Trans. Graph.
26, 3 (2007), 27:1–27:10. 2[TJCN06] TAMURA N., J OHAN H., C HEN B.-Y., N ISHITA T.: A
practical and fast rendering algorithm for dynamic scenes using
adaptive shadow fields. The Visual Computer (2006), 702–712.
2[TS06] T SAI Y.-T., S HIH Z.-C.: All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation. ACM Trans. Graph. 25, 3 (2006),
967–976. 1, 2, 6

Appendix
The polynomial parameters for g and h are listed in the following:g(η)=g4 (η/100)4 + g3 (η/100)3 + g2 (η/100)2 + g1 η/100 + g0 ,h(η)=h4 (η/100)4 + h3 (η/100)3 + h2 (η/100)2 + h1 η/100 + h0 ,(g4 , g3 , g2 , g1 , g0 ) = (−2.6856e−6 , 7e−4 , −0.0571, 3.9529, 17.6028),(h4 , h3 , h2 , h1 , h0 ) = (−2.6875e−6 , 7e−4 , −0.0592, 3.9900, 17.5003).

[VJ01] V IOLA P., J ONES M.: Rapid object detection using aboosted cascade of simple features. In Proceedings of IEEE
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.