用於球諧函數照明架構之有效取樣與重建入射光強度之方法

國

立

交

通

大

學

資訊工程系

碩

士

論

文

用於球諧函數照明架構之

有效取樣與重建入射光強度之方法

Effective Sampling and Reconstruction of Incident Radiance for

Spherical Harmonic Lighting

研究生：孫毓翔

指導教授：莊榮宏 博士

ày¦×Ðóï€Úx

b[ãø¥˜á ú—]°

@~ß Ò 

¼0>0: v }ÿ

»ñø;
.£G.

`Š

¦×Ðóï€ ×Ç`• *
ÍÈb[£Ý•îΛ3å՛VŽ;cXÝ G]Ùï€ìÝYŇÅ[Œ ݯhÚx|§ãÏÙXCWÝï €[Œ
&ÆèŒ×Íãø¥˜ÝÚx3GÝ]°
ÏÙXWÝá ú—Ðó
ÎDÄ×ÍG§¼í8÷µãø›H9°øÍ#½à×Í«û ÒbnÝÐó•/æ
|ÿÕΛ«N×ÍcFÝá ú—ãyøÍÝ÷µ ]P¬½á ú—‚Ž;
î/æ`XbÝøÍ¡GíÊ¼
á ú—Ý IŽ;P°W¥˜&ÆèŒ×͛VqAΛ¿¢|Cá ú —¼÷µøÍÝ]°
9Í]°|¯¸àï¼×ý0 Â
|-3øÍó ¥˜`² •×¨²
EyNÍøÍ
&ÆK›×Íb[Pš
¸ÿá ú—Ý IŽ;Èb[1D
¬vè>¥˜[£ i

Spherical Harmonic Lighting

Student: Yu-Hsian Sun

Advisor: Dr. Jung-Hong Chuang

Department of Computer Science and Information Engineering

National Chiao Tung University

ABSTRACT

Pre-computed radiance transfer accounting for self-shadowing and interreflection allows ob-jects to be shaded by distant, dynamic, low-frequency lights [SKS02]. In the previous work, incident radiance function is sampled at points uniformly selected by an off-line process, and a radial function is then used to interpolate all samples. If lights are near to objects, or occluding one another, a large number of samples are required. Since interpolation is performed through all samples, locality can be likely to lose. In this thesis, we present a framework of sampling and reconstruction in which objects can be shaded by nearby luminaries. Our work dynami-cally selects sample points according to object geometries and nearby luminaries. The quality of sampling is controlled by a user-specified error bound. Furthermore, the valid domain of samples can be used as error bounds to enhance the interpolation.

ACKNOWLEDGMENT

First of all, I would like to thank my advisor, Professor Jung-Hong Chuang and Mr. Wang-Yeh Lee for their guidance, support, and inspiration throughout my master’s degree. I am grateful to Mr. Shih-Ling Keng and Mr. Ling-Lin Shih for their comments and encouragements. Thanks to my colleagues in CGGM lab: Chih-Chun Chen, Tan-Chi Ho, Min-Sheng Chien, Ren-Hao Jen, Chi-Han Peng, Jong-Hon Lu, Yong-Cheng Chen, Yi-Gin Lin, Roger Hong, Cheng-Li Hou, and Chao-Wei Juan. It is always a great pleasure to study and exchange ideas with everyone in CGGM lab during the past two years. It is good to have my friends Chih-Hao Liang, Chin-Min Lin, Bo-Han Li, and Jin-Jen Huang to stand by my side. Last but not least, I would like to thank my family for their love and support.

List of Figures vi

List of Algorithms viii

List of Tables ix 1 Introduction 1 1.1 Literature Review . . . 2 1.2 Thesis Overview . . . 3 1.2.1 Problem Statement . . . 3 1.2.2 Contributions . . . 4 1.2.3 Thesis Organization . . . 4 2 Literature Review 6 2.1 Fundamentals of Physically-Based Rendering . . . 6

2.1.1 Spherical Coordinates . . . 6

2.1.2 Projected Area and Solid Angle . . . 7

2.1.3 Radiometry . . . 8

2.1.4 Light Reflection . . . 10

2.1.5 Formulations of Local Reflection Integral . . . 11

2.2 Reflection Integral and Spherical Harmonic Lighting . . . 12

2.2.1 Reflection Integral as a Convolution Process . . . 12

2.2.2 Orthonormal Basis Functions and Frequency Domain . . . 13

2.2.3 Fast Evaluation of Reflection Integral . . . 14

2.2.4 Precomputed Radiance Transfer . . . 17

2.2.5 Limitations of Spherical Harmonic Lighting . . . 18

2.3 Caching Techniques for Rendering . . . 20

3 Sampling and Reconstruction of Incident Radiance 24 3.1 Framework Overview . . . 24

3.2 Preprocessing Stage . . . 25

3.2.1 Precomputing Radiance Transfer Coefficients . . . 26

3.2.2 Bounding Volume Hierarchy Creation . . . 28

3.3 Runtime Stage . . . 30

3.3.1 The BVH and Its Use for Sampling . . . 30

3.3.2 The Oracle Function . . . 32

3.3.3 Hierarchy Traversal . . . 35

3.3.4 Incident Radiance Sampling . . . 36

3.3.5 Efficient Reconstruction using Octree . . . 37

4 Results 40 5 Conclusion and Future Work 49 5.1 Summary . . . 49

5.2 Future Work . . . 50

5.2.1 GPU Acceleration . . . 50

5.2.2 Higher-order Interpolator . . . 50

5.2.3 Oracle Improvement . . . 50

5.2.4 Real-time Rendering Applications . . . 51

Bibliography 52

2.1 Notations of spherical coordinates. . . 7

2.2 Three types of BRDF. Left to right: diffuse, specular, glossy [DBB03]. . . 10

2.3 Visualize spherical harmonics of various bands [RH01]. . . 15

2.4 A light probe and its approximation with 9 spherical harmonic coefficients [RH01]. 16 2.5 A directional function and its approximation with spherical harmonics. Left to right: the original function, and the approximation with l = (0, 2, ..., 10). The number of coefficients are 1, 4,...,100, respectively [Gre03]. . . 16

2.6 Approximate various BRDF and lighting condition with 25 coefficients [KSS02]. 17 2.7 Precomputed Radiance Transfer [SKS02]. . . 17

2.8 Reconstruct L(x, ωi) through ICP generated points. . . 19

2.9 An example of irradiance caching [PH04]. . . 21

3.1 Framework Overview. . . 26

3.2 BVH construction for the armadillo model. . . 30

3.3 Approximating an 1D function with discreet samples. . . 32

3.4 If the blockers and the receiver are near, higher variations in irradiance oc-curs [Arv94]. . . 34

3.5 Computing locational codes for nodes and vertices in a 2D quadtree. . . 38

4.1 Lighting configuration in Fig. 4.2. . . 41

4.2 Comparisons of the armadillo model for our methods, the previous work with ICP, and the ground truth. . . 42

4.3 Sampling distribution of our methods and the ICP algorithm. . . 43

4.4 Different lighting configurations for the church model. . . 43

4.5 Comparisons of the church model for our methods with  = 0.007,  = 0.004, the ICP method with 256 samples, 512 samples, and the ground truth. . . 46 4.6 Our methods detect the lighting configurations, and the sample numbers under

the same error bound are adjusted accordingly. Both lighting configurations are rendered with  = 0.001. . . 47 4.7 Results of different error bounds for the Buddha model. . . 48

3.1 Compute tipfor a vertex p . . . 27

3.2 Generate uniform direction samples over S2 . . . 28

3.3 Recursively Build BVH for O . . . 31

3.4 BVH Traversal . . . 35

List of Tables

3.1 Coordinate system of each cubemap face [Ope99]. . . 37

4.1 Error and timing comparisons of the armadillo model in Fig. 4.2. . . 41 4.2 Error and timing comparison of the church model in Fig. 4.5. . . 44 4.3 Sample numbers and costs under different lighting configurations with our method. 45 4.4 Errors and timing statistics for Fig. 4.7. . . 45

1

Introduction

With modern graphics hardware, complex lighting and shading effects seen only in traditional global illumination methods are now possible to be rendered in real time. One possible approach is the precomputed radiance transfer [SKS02], which precomputes complex transfer functions including occlusions and interreflections and encodes this function in textures as spherical har-monic coefficients. These methods allow for fast display of static objects with complex lighting conditions, assuming the distant light that is represented by a single environment map. However, several restrictions limit their applicable areas. For example, the object must not be deformed, the lighting environment can be rotated but not be translated, and the spatial relationship be-tween two objects can not be altered. These methods, however, shed some lights on the way to physically-based, high-quality real-time rendering on graphics hardware.

We present a framework of sampling and reconstruction that allows objects to be shaded by nearby luminaries under the framework of precomputed radiance transfer. In the previous work, incident radiance function are reconstructed at sample points uniformly selected by an off-line process and are interpolated using a radial function over all samples. This method gen-erally requires a large number of samples to fully reconstruct the incident radiance function. Furthermore, how many samples are needed is unknown. We propose a hierarchical sampling scheme that selects sample points according to object geometries and nearby luminaries, result-ing in dynamically scalable image that has better quality compared to previous methods. We also bound the errors with the valid domain of samples to better the interpolation.

1.1 Literature Review 2

1.1

Literature Review

Real-time rendering has lots of important applications and is currently an active research area. Traditional graphics hardware has many restrictions. For example, light sources are limited to point light sources. The lighting model is limited to the Phong model. As graphics hard-ware evolves, we are possible to render a more physically-based, realistic images in real time. To shade a point of interest, one must evaluate the local reflectance integral, which is an in-tegration of the incident radiance function L(x, ωi) and a bidirectional reflectance distribution function(BRDF) that models surface reflectance properties over the entire incident hemisphere. The BRDF is a known function, but incident radiance is not. We can separate the problem of the evaluation of the local reflection integral into two: one is how to compute L(x, ωi) and the other is how to efficiently compute the integral itself.

The incident radiance function L(x, ωi) can be decomposed into direct and indirect nations. Direct lighting can be obtained through the monte carlo integration. Indirect illumi-nation can be obtained by global illumiillumi-nation algorithms such as path tracing [Kaj86], photon mapping [Jen01], or radiosity. Indirect illumination is usually ignored in rendering algorithms based on traditional graphics hardware due to its high computational costs, but can be computed by performing radiosity simulation or approximated by dynamic calculations of ambient occlu-sionsor some multi-pass rendering algorithms. Simulating the global transportation of light requires global knowledge of the scene, which is expensive to derive. Precomputation makes it feasible to add global illumination effects into real-time rendering applications. Daubert et al. precomputes the interreflection between surface mesostructures and parameterize the result with viewing and lighting direction [DKS+03]. By this way the dynamic interreflection on sur-face can be efficiently evaluated with a few texture lookups. Hao et al. precomputes the dipole approximation integral for subsurface scattering at many lighting directions, and store the result with a reference point scheme to compress the data set [HV04]. Sloan et al. proposed a method that pre-calculates self-transfers such as occlusion and interreflection and stores the transfer function as spherical harmonic coefficients for each vertex [SKS02]. Even if we consider only direct illumination, we still have to evaluate the local reflection integral. A physically accurate approach to compute the integral is ray tracing with monte carlo integration, as proposed by Kajiya [Kaj86] and Cook [CPC84]. However since we must gather radiance values from all directions to shade a point of interest, the coherence between different rays is so low that the

number of rays required may be large. This leads to a slow process and is not inherently suitable for graphics hardware. This is why most graphics hardware supports only a finite number of point light sources. We can intuitively think that evaluating the reflection integral with ray trac-ing is to compute the convolution of the incident radiance function and the surface BRDF in the directional domain, which is inevitably a slow process. Ramamoorthi and Hanrahan proposed a framework that evaluates the reflection integral in the frequency domain[RH01, RH04] and justified the work done by Sloan and Kautz et al. [SKS02, KSS02]

1.2

Thesis Overview

1.2.1

Problem Statement

In the framework of Sloan et al. [SKS02], distant light sources are assumed and in consequence the radiance field does not vary along surface positions. Incident radiance function L(x, ωi) on each point x is represented with only a single directional function Lp(ω) in spherical harmonics, where p is usually the center of the shaded object. The assumption does not hold true if the light source is not infinite far away, and such violations may result in incorrect shading. It is indicated that the problem posed by non-distant lights could be solved by a multiple sampling of the radiance function across the surface with ICP algorithms, but no justifications is given. Since what we want is a correct shading result, one must take into account both the variation of the incident radiance function L(x, ωi) and the transfer function T (x, ω) to decide the number and distribution of the samples.

We developed a sampling scheme that deals with non-distant light sources. In general,

L(x, ωi) varies both in positions and directions and can be very complex due to occlusions and complex light sources. As in [SKS02] and [AKDS04], for each sample point xiin space we can compute its incident radiance function Lxi(ω) and represent it in a vector of spherical harmonic

coefficients. In this framework, each L(x, ωi) is point-sampled in discreet points xiand recon-structed from these sample points. The error introduced by the approximation is determined by two factors: the sampling distribution and the reconstruction filter. The work by Annen et al. [AKDS04] is a higher-order reconstruction filter. Higher-order reconstruction filters can improve shading results, but cannot capture the local properties of the approximated function. In the case of the incident radiance function, changes due to occlusions and complex lighting geometries cannot be captured by higher-order reconstruction filters. To capture the local

prop-1.2 Thesis Overview 4

erties of the function, more samples must be placed at proper locations. For example, more samples should be placed on the boundary of shadows, as in the work by Ward et al. [WRC88]. We want to address the problem of sample distribution and reconstruction. Our framework dy-namically decides the sample number and distribution, accounting for the variation of L(x, ω). Since it is possible that the variation of L(x, ω) is too high to be reconstructed in desired com-putational cost, the framework are able to adjust the sample number, given user-specified error bounds or computational costs.

To dynamically distribute samples and to decide the number of samples needed, a bounding volume hierarchy of the object is created as an off-line process. During run time, we evaluate an estimated error when the incident radiance function for each bounding volume hierarchy node is being sampled. If the estimated error is larger than a user-specified error bound, more samples are added and the positions of samples are determined by the hierarchy. The incident radiance function is reconstructed from samples with the same interpolator used in Sloan et al’s work [SKS02, AKDS04]. However, to achieve more efficient reconstruction and to capture local variations of lighting, a valid domain is associated with each sample, and only valid samples are used for interpolation. Moreover, the samples are inserted to an octree for efficient queries, when they are generated.

1.2.2

Contributions

Our contributions can be concluded as follows:

• A hierarchical object-space sampling scheme for spherical harmonic lighting that

sup-ports non-distant lights.

• A spatial partitioning data structure that accelerates the reconstruction of incident

radi-ance functions.

• A mechanism that allows users to control the trade-off between the shading quality and

the rendering cost.

1.2.3

Thesis Organization

In Chapter 2, we introduce notations and the backgrounds required for our work. In Chapter 3, we describe the framework and theoretical details of our work, while Chapter 4 demonstrates

the results, which are compared to previous methods, and in Chapter 5, we discuss possible improvements and the future work.

C H A P T E R

2

Literature Review

This chapter lays the groundwork for physically-based shading on modern graphics hardware. We start from the foundations of physically-based rendering within a traditional realistic image synthesis framework. We then show how spherical harmonics and precomputed radiance trans-fer allow for fast physically-based shading on modern graphics hardware. Finally, we discuss existing techniques for caching radiance fields and their theoretical backgrounds.

2.1

Fundamentals of Physically-Based Rendering

First we introduce the physical quantities for lighting and rendering. Then we discuss the math-ematical formulation that models the interaction between the light and the surface. Finally we discuss how modern graphics hardware is capable of converting simulation results into dis-playable low-dynamic-range pixels.

2.1.1

Spherical Coordinates

We are interested in quantities about light. Lighting differs from each location and each di-rection. To define a location in space we use a Cartesian system. In this thesis we assume a right-handed coordinate system. A point in space, usually denoted as p in this thesis, can be ex-plicitly written as an element in R3_{space, (x, y, z). To denote a vector, we use notations with an} arrow above an English alphabet, like ~s. For pure directional quantities, we simply use the term

direction. In this thesis we denote a direction as ω, and a differential solid angle around a certain

direction as dω. Solid angle will be introduced later. The set of all possible direction, which

* + , / -.

Figure 2.1: Notations of spherical coordinates.

can be imagined as directions from the origin to surface points on a unit sphere, is denoted as

S2. When calculating a point’s outgoing radiance due to reflection, we only get interested the upper hemisphere above the plane formed by the normal at the point of interest. We denote the set of directions of the upper hemisphere as Ω. Though a direction is an unique element in S2 domain, in practice we often parameterize the directional domain into horizontal angle, φ, and vertical angle, θ, as shown in Fig. 2.1. The vertical angle θ is the angle between the direction vector ω and z-axis. The horizontal angle φ is the angle between the projection of ω onto the xy plane and the x-axis. Most directional functions, such as spherical harmonics, can be explicitly defined as a function over θ and φ. To convert between the angle notation (θ, φ) and normalized Euclidean vector notation ω = (x, y, z) we use the following relations:

cosθ = z; θ = cos−1(z); sinθ cos φ = x; φ = tan−1(y

x); sinθ sin φ = y.

2.1.2

Projected Area and Solid Angle

Projected area is defined as the orthogonal projection of a surface of any shape onto a plane

with a unit vector v as normal. The differential form is dAproj = cos(θ)dA, where θ is the angle between the local surface normal and the unit vector v. We can integrate dAprojover the visible surface area to get the total projected area:

Aproj =

Z

A

2.1 Fundamentals of Physically-Based Rendering 8

Solid angle is an extension of plane angle from two dimensions to three dimensions. Recall the

definition of a plane angle is “One radian is the plane angle between two radii of a circle that cuts off on the circumference an arc equal in length to the radius.” And the definition of a solid angle is extended to: “The solid angle subtended by an object from a point P is the area of the projection of the object onto the unit sphere centered at P [SP94]”. Solid angle is expressed in steradians (sr). The solid angle subtended by the whole sphere is 4π sr, that is, the entire area of a unit sphere. The solid angle of an object is thus the area of the projection of the object onto a unit sphere. Note that two objects different in shape can still subtend the same solid angle. We can think of the differential solid angle as representing both a direction and an infinitesimal area on the unit sphere. We use dω to denote the differential solid angle around a certain direction

ω. Sometimes given a point P in space, we need to convert between differential solid angle

and differential surface area. Let r be the distance between P and the differential area dA, and

θ is the angle between normal of dA and the vector from dA to P , the differential solid angle

subtended by dA is:

dω = cos θdA

r2 . (2.1)

If dAsis a small differential surface element on a sphere of radius r, then

dAs= r2sinθ dθdφ,

and

dω = dAs

r2 = sinθdθdφ.

2.1.3

Radiometry

Physically-based rendering is the simulation of light. The basic terminology to describe light is radiometry, a measurement of optical radiation, which is electromagnetic radiation within the frequency range between 3∗1011_{and 3∗10}16_{Hz. This range corresponds to wavelength between} 0.01 and 1000 micrometers (µm) and includes the regions commonly called the ultraviolet, the visible, and the infrared.

The most basic quantity in radiometry is the photon. The energy eλ of a photon with wave-length λ is

eλ =

hc λ ,

where h ≈ 6.63 · 10−34J · sec is Planck’s constant, and c is the speed of light. eλ is measured in joules (J ).

Spectral radiant energy Qλ in nλ photons with wavelength λ is

Qλ = nλeλ = nλ

hc λ.

Radiant energy Q is the energy computed by integrating the spectral radiant energy over all

possible wavelengths:

Q = Z inf

0

Qλdλ.

Radiant flux or radiant power Φ describes the flow of radiant energy per unit time:

Φ = dQ dt ,

where t is measured in second. Radiant flux is often just called flux.

Radiant flux density M or B or E is defined as the differential flux per differential area at

certain surface location x:

M (x) = B(x) = E(x) = dΦ dA,

where M is referred to as radiant existence, which is the flux leaving a surface, B is referred to radiosity, which is exactly the same as radiant existence; and E is referred to as irradiance, which is the flux arriving at a surface.

Radiant intensity I is defined as the differential flux per differential solid angle dω:

I = dΦ dω.

The essential measurement used in the context of rendering is radiance, L, defined as the dif-ferential flux per difdif-ferential projected area per difdif-ferential solid angle:

L = d

2_Φ

cosθdAdw = dE cosθ dω.

Radiance can be described by a five-dimensional function (three for the position and two for the direction) per wavelength, usually written as Lλ(x, ω). It is the most important quantity in radiometry, since it could most closely represent the color of an object. Most light receivers, such as cameras and the human eye, are sensitive to radiance, while the response curve of these sensors may be different. In practice, we often simplify per-wavelength calculation into tristimulus values, LR(x, ω), LG(x, ω), and LB(x, ω). In this thesis, we simply use L(x, ωi) and did not distinguish between them, though in the implementation we have to process these three values, respectively.

2.1 Fundamentals of Physically-Based Rendering 10

2.1.4

Light Reflection

In this section we introduce the theoretical framework used to model the interaction between light and object surfaces. Radiant energy will be scattered or absorbed by surfaces. We consider only reflection here.

Reflection is a behavior that light enters and leaves the material at the same point on the sur-face. The amount of energy been reflected can be modeled by a six dimensional function called

bidirectional reflectance distribution function, or BRDF. BRDF is the ratio between irradiance

and outgoing radiance, usually denoted as ρ [DBB03]:

ρ(x, ωi, ωo) =

dL(x, ωo)

dE(x, ωi)

,

where x is the 2D position on the surface, ωo is the differential solid angle around outgoing direction of radiance, ωi is the differential solid angle around incident direction of irradiance. BRDF becomes a 4D function for a given surface location, and a 2D directional function if we fix both the surface location and outgoing direction. We denote a BRDF with fixed outgoing direction as ρωo(x, ωi) and a BRDF with both outgoing direction and surface location fixed as

ρpωo(ωi). BRDF can be roughly categorized into the following three types, as shown in Fig. 2.2:

diffuse Reflect light uniformly over the entire reflecting hemisphere. The BRDF is simply a

constant, also known as the Lambert’s model.

specular A perfectly mirror surface that reflects incident radiance into mirror directions. The

BRDF is a Dirac delta function.

glossy Most surfaces are neither ideally diffuse nor ideally specular and are called glossy

sur-faces. Their BRDFs are much more complex.

Given a point x on the surface, the incident radiance function L(x, ωi) at the point x, and the BRDF ρ(x, ωi, ωo), we are interested in computing the outgoing radiance L(x, ωo) on a given direction ωo. We start from the definition of BRDF, integrate over the hemisphere of directions on the right side, finally arriving at the local reflection integral [RH04]:

ρλ(x, ωi, ωo) = dL(x, ωo) dE(x, ωi) dL(x, ωo) = ρλ(x, ωi, ωo)dE(x, ωi) L(x, ωo) = Z Ωx ρλ(x, ωi, ωo)dE(x, ωi) L(x, ωo) = Z Ωx ρλ(x, ωi, ωo)L(x, ωi)cos(Nx, ωi)dω, (2.2)

where Nx is the surface normal at the point x, and Ωx is the hemisphere around the surface normal Nx. The reflected radiance field is given by the integral above. It is an important formula for rendering and will be evaluated at each shading point. As we can see, it is an integration of the incident radiance function and the surface BRDF weighted by the projected area over all possible directions in Ω.

2.1.5

Formulations of Local Reflection Integral

The local reflection integral can be rewritten in various forms. In previous section we introduced a form that integrate over the hemisphere. This is called the hemispherical formulation. Another possible form is the area formulation [AKDS04]. In this formulation we consider the surfaces of objects in the scene that contribute to the incoming radiance at a point x. We integrate over all surface points on the scene surface A, taking the occlusion and geometry relationship into account. The occlusion can be modeled as a binary visibility function between two points x and

y [DBB03]: V (x, y) =     

1 if x and y are mutually visible,

0 if x and y are not mutually visible,

∀x, y ∈ A.

The geometry relationship of two surface points x and y, called the geometry term [DBB03], is:

G(x, y) = cos(Nx, Ψ)cos(Ny, Ψ) r2

xy

,

where Nx and Ny are surface normals at point x and y, respectively, and Ψ is the normalized unit vector from x to y. With the geometry term and the visibility function, we can write the

2.2 Reflection Integral and Spherical Harmonic Lighting 12

local reflection integral as:

L(x, ωo) =

Z

A

L(x, ωi)ρ(x, ωi, ωo)V (x, y)G(x, y) dAy.

The area formulation is useful in computing the lighting from nearby geometries, while the hemispherical form is easier to analyze in the frequency domain. We will have further discussion in the next section.

2.2

Reflection Integral and Spherical Harmonic Lighting

Computing the reflection integral in the frequency domain is an efficient approach. To see why and how this works, we will first show that the reflection integral can be viewed as a convolution in the directional domain. Then we will demonstrate that through orthonormal basis functions, we can project functions defined in the directional domain into the frequency domain, where convolution simply becomes a multiplication. The frequency domain analysis of reflection integral is thoroughly elaborated in the work of Ramamoorthi et al. [RH04].

2.2.1

Reflection Integral as a Convolution Process

The definition of convolution of two functions f (x) and g(x) is given as:

f (x) ⊗ g(x) = Z

f (x)g(t − x) dx.

Convolution is usually performed in the entire domains of f (x) and g(x). In the domain of real numbers, we integrate over [−∞, ∞]. In the directional domain, the entire sphere of directions,

S2_{, is considered instead. Moreover, we can generalize the notion of convolution to some other} transformation Rt[RH04], where Rtis a function of t and write

(f (x) ⊗ g(x))(t) = Z

x

f (x)g(Rt(x)) dx.

When Rt is a rotation by angle t, the above formula defines the convolution in the angular domain. In the case of reflection, there is no value on the hemisphere below the surface; hence we integrate only over Ω, the upper hemisphere.

As in [SKS02], we define a transfer function Tp,ωo(ωi) for a given point p on surface and

a give viewing direction ωoas the multiplication of the binary visibility function Vp(ωi) repre-senting the occlusion at point p, the BRDF ρp,ωo(ωi), and the projected solid angle cosθ. For

diffuse surfaces, the BRDF does not depend on the viewing direction, so we can discard the ωo term in the transfer function. Note that Vp(ωi) is defined as the self occlusion of the object:

Vp(ωi) =     

1 if p is not occluded in direction ωi,

0 if p is occluded in direction ωi. Now the transfer function for diffuse surfaces is:

Tp(ωi) = Vp(ωi)ρp(ωi)cosθ. (2.3)

To see why Equation (2.2) can be seen as a convolution, we can rewrite Equation (2.2) as the multiplication of the incident radiance function and the transfer function:

Lp(ωo) =

Z

Ω

Lp(ωi)Tp(R0(ωi))dωi (2.4)

= Lp(ωi) ⊗ Tp(ωi), (2.5)

where p is the point of interest and R0is an identical rotation.

We can now consider Lp(ωo) as a convolution of the incident radiance Lp and the trans-fer function Tp. Intuitively, a convolution in spatial domain will become a multiplication in frequency domain. To transform Lp and Tp into the frequency domain, we need a set of ortho-normal basis functions that span the entire Ω domain.

2.2.2

Orthonormal Basis Functions and Frequency Domain

A pair of functions fi(x) and fj(x) are said to be orthonormal to each other if they are normal-ized and orthogonal to each other. A function fi(x) is normalized if

Z b

a

fi(x)fi(x) dx = 1.

Two functions fi(x) and fj(x) are said to be orthogonal to each other if

Z b a

fi(x)fj(x) dx = 0.

We may say that if a set of function {fi(x)} are orthonormal to one another, then the convolution of any two functions fi and fj in the set is exactly a Kronecker delta function δij:

Z fi(x)fj(x) dx = δij =      1 if i = j, 0 if i 6= j.

2.2 Reflection Integral and Spherical Harmonic Lighting 14

Let {Bi(x)} be a complete set of orthonormal functions. Assume that Bi(x) and f (x) have the same range. The function f (x) can be expanded as

Fi = Z f (x)Bi(x) dx, (2.6) where f (x) =X i FiBi(x). (2.7)

We say that Fi are the coefficients obtained by projecting f (x) onto the basis functions Bi(x). After projecting two functions, f (x) and g(x), onto the same set of basis functions, the convolu-tion of f (x) and g(x) will become the sum of products of their respective coefficients [NRH04], that is, f (x) ⊗ g(x) = Z f (x)g(x)dx = Z X i FiBi(x) ! X j GjBj(x) ! dx =X i X j FiGj Z Bi(x)Bj(x) dx =X i X j FiGjδij =X i FiGi. (2.8)

If we project Lp and Tp in Equation (2.2) onto some basis functions over the hemispherical domain, the integral will become simply a dot product between their coefficients. Such a set of functions as Bi(x) is called spherical harmonic functions.

2.2.3

Fast Evaluation of Reflection Integral

Spherical harmonics is used extensively in computational chemistry and introduced to the con-text of computer graphics by Cabral et al. [CMS87]. They form a complete orthonormal basis set for S2_{domain. Spherical harmonics, denoted as Y}l

mhere, are actually the product of terms of Fourier basis functions cosθeimφ, associated Legendre polynomials P_l|m|, and a normalization factor Km

l , as follows:

Y_lm(θ, φ) = K_lmP_l|m|(cosθ)eimφ.

The indexing variable l ≥ 0 and −l ≤ m ≤ l can be thought as the frequency of the basis function, where m controls the frequency in horizontal directions. Since −l ≤ m ≤ l, we can

Figure 2.3: Visualize spherical harmonics of various bands [RH01].

combine l and m into a single indexing variable i with i = l(l + 1) + m, where i ≥ 0. Though

Yiforms a complete system only when i → ∞, most directional functions can be approximated by a small number of basis functions.

To compute the exact value of a spherical harmonic Yiat certain direction, we must evaluate the associated Legendre polynomial, multiplied by the normalization factor and terms of Fourier basis function. Normalization term is defined as follows [Gre03]:

K_ml = s (2l + 1) 4π (l − |m|)! (l + |m|)!. Pm

l (x) is computed with the following recursion rules [Gre03]:

P_mm = (−1)m(2m − 1)!!(1 − x2)m/2 (2.9)

(l − m)P_ml = x(2l − 1)P_l−1m − (l + m − l)P_l−2m (2.10)

P_m+1m = x(2m + 1)P_mm (2.11)

Both Rule (2.10) and (2.11) can be used to compute a higher band (with larger l) solution from a lower band one. Rule (2.11) introduces more floating point roundoff error and Rule (2.10) is used whenever possible. Rule (2.9) is the termination rule of the recurrence relation, and

k!! means the product of all odd integers less than or equal to k. Fig. 2.3 visualizes spherical

harmonics with l = (0, 1, 2). Since m ≥ −l and m ≤ l, there are nine spherical harmonic func-tions. Positive values are in green and negative values are in blue. To project a known function defined over all directions, f (ω), onto spherical harmonic coefficients, one usually evaluates the integral in Equation (2.6) with numerical methods like the quadrature rule. Fig. 2.4 demon-strates projecting an environmental map that represents incident radiance from the surrounding environment onto spherical harmonic coefficients, and reconstruct the map using Equation (2.7). Note that how sharp lighting details are blurred, since only 9 coefficients are used to approx-imate the environmental map. As the number of coefficients increases, we can preserve more on sharper features. This is shown in Fig. 2.5. At first a small number of coefficients seems to

2.2 Reflection Integral and Spherical Harmonic Lighting 16

Figure 2.4: A light probe and its approximation with 9 spherical harmonic coefficients [RH01].

Figure 2.5: A directional function and its approximation with spherical harmonics. Left to right: the original function, and the approximation with l = (0, 2, ..., 10). The number of coefficients are 1, 4,...,100, respectively [Gre03].

be a very rough approximation, but taking surface material into consideration, it is often suf-ficient. In [RH01], it is shown that for Lambertian surfaces, nine coefficients are sufsuf-ficient. It is proved that for lambertian surfaces, the error introduced by this approximation is negligible in [RH04]. Kautz et al. proposed a method that stores the spherical harmonic coefficients of BRDF ρωo(x, ωi) on a texture[KSS02], parameterized by viewing angle ωo, to represent

arbi-trary BRDF. The result is shown in Fig. 2.6. In their work, 25 spherical harmonic coefficients is sufficient to represent surfaces from perfect lambertian to anisotropic glossy metallic one for each viewing direction. However, their method cannot handle specular surfaces.

Now we know that the reflection integral is a convolution of the incident radiance function and the transfer function. We can evaluate the convolution by projecting both functions onto spherical harmonics to obtain two vectors of spherical harmonic coefficients, and the reflected radiance becomes the dot product of two coefficient vectors as Equation (2.8) shows. If we assume all emitters are infinitely far away, then the incident radiance function does not depend on the position. This is an assumption commonly used in environmental-map based rendering techniques. Now the the incident radiance function can be represented by a vector of spherical

Figure 2.6: Approximate various BRDF and lighting condition with 25 coefficients [KSS02].

Figure 2.7: Precomputed Radiance Transfer [SKS02].

harmonic coefficients, and can be quickly convoluted. Since the transfer function Tωo(x, ωi) is

different at each shading point, we must precomputed them. In the next section, we introduce previous works on precomputing the transfer function.

2.2.4

Precomputed Radiance Transfer

Various representations may encapsulate precomputed or acquired global illumination solu-tions. Light maps store radiosity simulation results in surface diffuse textures and add global il-lumination to interactive 3D environments. Surface light fields [CBCG02] record 4D exiting ra-diance sampled over an object’s surface, but both the emitters and the scene must remain static. Horizon map [Max88] stores precomputed visibility information for surface micro-geometry, which can efficiently render self-shadowing effects. Polynomial texture maps [MGW01] fit acquired or precomputed surface BRDF into polynomial functions of incident light vector, al-lowing for real-time interreflection effects. However, both horizon map and polynomial texture map are limited to point light sources, since for area light sources costly multi-pass rendering is required. As shown in Fig. 2.7, Sloan et al. [SKS02] precompute the transfer function for each vertex, and compute the dot product on graphics hardware, resulting in real-time rendering with effects like occlusions and interreflections on modern graphics hardware. They call T the

2.2 Reflection Integral and Spherical Harmonic Lighting 18

radiance transfer function or simply the transfer function, and we follow their denotation. The precomputed radiance transfer framework proposed by Sloan et al. [SKS02], also known as spherical harmonic lighting [Gre03] that computes Equation (2.2) on frequency domain by projecting both the incident radiance function and the transfer function onto spherical harmonic basis functions. Sloan’s work extends the work of Ramamoorthi et al. [RH01], which consid-ers only diffuse surface and precomputes self-shadowing and self-interreflection. Their goal is to shade a object with an environmental map as distant light sources in real time. They also separate surface BRDFs out from the transfer function, thus allowing dynamic changes of sur-face materials. The framework is further extended by various researchers to address arbitrary BRDFs [KSS02] and compress high-dimensional transfer functions by clustered principal com-ponent analysis [SHHS03]. To prevent the transfer function from being recomputed when the surface BRDF is changed, they further decouple the BRDF from the transfer function. By this way the dot product becomes a matrix multiplication. For details, see [SKS02] and [SHHS03]. Ng et al. [NRH04] use haar basis functions defined on directional domain instead of spherical harmonics to handle all-frequency lighting conditions, including specular reflections.

2.2.5

Limitations of Spherical Harmonic Lighting

Though spherical harmonic light provides an efficient way for representing and computing radi-ance functions, it has its weaknesses and limitations. Spherical harmonics is very inefficient for high-frequency lighting conditions such as point light sources and perfect specular reflections. Hence, the scene must remain static, and only infinitely far environmental light sources can be handled.

Spherical harmonic lighting cannot handle point light sources and specular reflections be-cause of its frequency domain nature. Remember that though {Yi} forms a complete ortho-normal basis system, {Yi} is a infinite set; and in practice we can only compute a finite set of coefficients to represent function f (ω) defined over the directional domain, such as the incident radiance function at some point p. If f (ω) contains high frequency components, a finite set of coefficients can not be able to accurately represent f (ω). Point light sources and specular reflections are Dirac delta functions in directional domain, thus requiring an infinitely large number of coefficients. Therefore spherical harmonic lighting is suitable for slow varying, low-frequency lighting conditions. Luckily most natural lighting environments satisfy this condi-tion. To overcome this problem, we can choose another set of basis functions that can represent

(a) Ground truth (b) Sample location (c) Result of reconstruction

Figure 2.8: Reconstruct L(x, ωi) through ICP generated points.

high-frequency signals more efficiently, as in [NRH04] and [LSSS04].

The second limitation is the nature of precomputation. Since the transfer function are pre-computed over an object, once the object is deformed, the prepre-computed coefficients are in-valid. One possible solution is to encode all possible deformations in advance and use principal component data reduction to relate pose deformation to the change of transfer function coef-ficients [JF03]. However this solution is only applicable to objects that are not deformed too much–rubber models, for example.

The third limitation, also the problem we’d like to address in this thesis, comes from the fact that we assume that emitters are infinitely far away. If there are emitters near the object being shaded, the incident radiance function L(x, ωi) may vary rapidly along x. The previous work [SKS02] simply stated that we may solve this problem through multiple samples over objects and interpolate through linear combinations of coefficients with the following Equa-tion [AKDS04]: Lp = X i wiLi, (2.12) where wi = 1 kp−xikb P j 1 kp−xjkb . (2.13)

In Equation (2.13), x is the location of a sample point, Li is a coefficient of a sample, and

b controls the locality of reconstruction result. Both [SKS02] and [AKDS04] distribute the

sample with iterative closest point algorithm(ICP) [LBG80], which is a vector quantization algorithm also known as LBG because of its inventor. This algorithm distributes samples evenly

2.3 Caching Techniques for Rendering 20

across all vertices, though the samples are unnecessarily fall on mesh surface. While this may be a good choice if the emitter is not close to the object, it does not give any relationship between the quality, the number of samples, and error bounds. As the emitters draw near to the object or occluding one another, very large number of samples are required to capture the local details of the incident radiance. Since L(x, ωi) is reconstructed by interpolating all sample points, reconstruction is slow when there are many samples. The effect of b in Equation (2.13) relies on the distance between samples, which is hard to control and predict. As shown in Fig. 2.8, the armadillo model is illuminated by two distant bluish light sources and one yellow nearby light source. Though 64 samples are used, the lighting near the head of armadillo is missing. Ground truth is obtained by sampling L(x, ωi) at every vertex.

Another problem introduced by non-distant light sources is that transfer function stores only self-occlusion or self-interreflection information. If the emitter intrudes the object’s convex hull, or there are multiple objects close to one another, the precomputed transfer function becomes invalid. In this thesis we address the problem of varying L(x, ωi). The problem of invalid

T (x, ωi), however, may need further investigation.

2.3

Caching Techniques for Rendering

Computing the exact value L(x, ωi), whether or not indirect illumination is included, is a rather expensive operation. There are several ways to accelerate this operation, including precompu-tation and sparse sampling. We discuss sparse sampling here.

The idea behind sparse sampling can be considered as a function approximation from dis-crete sample points. L(x, ω) is an unknown function, but can be evaluated at certain points

(xi, ωi) where xi ∈ R3 and ωi ∈ S2. We can construct another function L0(x, ω) from each point (xi, ωi), trying to minimize the difference between L and L0. The difference between two function f (x) and g(x) is usually measured by Lp norm, as defined in [Wat80]:

normp(f, g) =k f − g kp=

Z

| f (x) − g(x) |p −p

.

The error introduce by the approximation is related with three factors: the number of samples, the filter used to reconstruct the function, and the distribution of samples. Increasing the sample density or simply the sample number leads to better approximation at the expense of higher sampling costs. If function sampling is expensive, then we should avoid evaluating too many

(a) Scene rendered with irradiance caching (b) Distribution of caches

Figure 2.9: An example of irradiance caching [PH04].

samples; and proper sample distribution becomes important. After evaluating samples, to con-struct the approximating function, we must apply some kind of interpolator or filter to existing samples. Possible choices include linear combination and piecewise polynomial approximating functions. In general, high order interpolator gives smooth results, but needs more samples as input and extra computational costs. Several methods exploiting the sparse sampling over spatial domain. We briefly introduce them since they are somehow related to our work. Ward et al. [WRC88] proposed a method that greatly improves the rendering speed for ray tracing diffuse surfaces. The idea is to sparsely sample the incident indirect irradiance function over surfaces. In diffuse environments, irradiance Eindirect(x) due to indirect illumination cannot change rapidly, hence justifing the motivation of sparse sampling. To decide the sample distri-bution, Ward et al. proposed a split sphere model that decides the valid domain of each sample point, which is called irradiance cache in their work. The split sphere model has the largest gra-dient possible for an environment without concentrated sources. During the ray tracing process, at each shading point p (the first hit point of primary ray), they check if there is any usable cache nearby. Nearby usable caches are then used to interpolate the irradiance at p. If there is no usable cache nearby, a new cache is generated at p and stored into an octree constructed from scene geometry for efficient queries. Since each cache generation involves computing irradi-ance E(p) at p which dominates the computational costs of ray tracing diffuse environments, we can control the quality and speed by scaling the valid domain of each cache. The irradiance

2.3 Caching Techniques for Rendering 22

caching is very suitable for rendering scenes containing numerous diffuse objects. Note that cache distribution is decided by the ray tracing order of pixels and the valid domain of each cache, although the caches themselves are stored in R3 space. As shown in Fig. 2.9, caches are placed densely where indirect illumination changes rapidly, e.g., corners of walls. The blocky artifacts on the pillars and walls due to insufficient sampling can be compensated by computing gradients of irradiance during cache sampling and by using the gradients when interpolating irradiance values [WH92].

In [WRC88], they reconstruct the irradiance E(x) at each point of interest x as the weighted sum of irradiance values Ep of caches [PH04]:

E(x) = P pwpEp P wp , where wp =  dmax(Np·Nx)−kx−pk dmax(Np·Nx) 2

, and dmax is a user-specified upper bound of cache’s valid domain.

While the irradiance caches are stored in R3 space separated from the scene geometry, they represent irradiance of each single point on a surface. These information may become invalid as we move away from the sample position, or the surface normal is rotated. Irradiance volume [GSHG98] generalizes the idea of irradiance caching by storing irradiance values for all possible directions of normals at each sample point (densely sampled for around 200 directions) and arranging sample points in a bilevel uniform grid. Also, irradiance volume stores both direct and indirect irradiance, while irradiance cache stores only indirect irradiance. Greger et al. [GSHG98] define a semi-dynamic environment in which most surfaces are stationary, and a few dynamic objects are small relative to the scale of the environment. For example, positioning a piece of furniture in an architectural application would fall into this category. With irradiance volume, not only irradiance on existing surfaces can be interpolated, but also a few dynamic objects moving in this semi-dynamic environment can be shaded by irradiance volume computed from the stationary geometries. However, a few moving objects cannot affect the shading results of the stationary geometries, since the irradiance volume is precomputed and static.

The sample density in an irradiance volume is decided by the number of primitives falling inside a grid cell. If the number of primitives inside a grid cell exceeds some user-specified limit, the cell is further divided into another uniform grid, resulting in a bilevel uniform grid. The sample positions are corner points of grid cells. To shade a point p with normal Np, one

must query the irradiance volume to find its corresponding irradiance. To query the irradiance at p with normal Np, the cell containing p is accessed to find the sample values corresponding to the direction Np. Since the cell has eight corner points, a trilinear interpolation is used to reconstruct the final irradiance from the eight values retrieved from the cell. Unlike Ward et al’s work, irradiance volume gives no error bounds for the interpolation.

Both irradiance caching and irradiance volume require a long computational time, and the sampling is done off-line. On the other hand, our work targets interactive 3D environments under the framework of spherical harmonic lighting. Moreover, since irradiance cache and irradiance volume store irradiance rather than directional radiance, they are limited to diffuse environments. Our work adopts the idea of the valid domain from irradiance caching, but use a different sampling strategy. Irradiance caching places the sample points in the progress of ray tracing, whereas our method distributes our samples as the object hierarchy is traversed. Our method has some similarities with irradiance volume as we both place samples in R3rather than on object surfaces. Unlike irradiance volume, we dynamically determine the sample distribution and the number of samples; and sample the incident radiance function at runtime. Also, we store directional radiance values (in spherical harmonic coefficients) rather than irradiance values, and in consequence, the method can be extended to glossy materials as in [KSS02].

C H A P T E R

3

Sampling and Reconstruction

of Incident Radiance

3.1

Framework Overview

For distant lights, the incident radiance function L(x, ωi) could not change a lot around x and hence vertices can use the same incident radiance coefficients. For non-distant lights, this is no longer true. Our goal is to incorporate non-distant light sources to the spherical harmonic lighting. To prevent from sampling the incident radiance at each vertex, we use the sparse sampling and from which to reconstruct the spherical harmonic coefficients of L(p, ω) for a vertex p. Spare sampling is error-prone, but the error will be decreased when the number of samples increases. The problems are that we have to decide where should the samples be placed and when to stop sampling. To dynamically distribute samples and decide the required number of samples, a bounding volume hierarchy (BVH) of the object is constructed during the preprocessing stage. At runtime we traverse the hierarchy, and at each node of the hierarchy, we sample the incident radiance and compute an estimated error. If the error is lower than a user-specified error bound, we stop traversing the hierarchy and terminate the sampling stage. After stopping traversing the hierarchy, the spherical harmonic coefficient vectors representing the incident radiance of every vertices are then reconstructed from those sparse samples.

Previous methods interpolate all samples for the value at a vertex [SKS02, AKDS04]. If the number of samples is larger than a hundred or more, interpolating through all samples becomes very slow. In our framework, to preserve the locality of lighting and to accelerate

the reconstruction speed, each sample is associated with a valid domain. During reconstructing the coefficient vector for the incident radiance at a vertex, only samples whose valid domain contains the vertex will be considered.

Our framework considers the whole scene geometry as a static object that can be lit by dynamic emitters. Both the object and the emitters are represented by triangle meshes. The scene geometry is denoted as O. Geometry of emitters is denoted as Gl. A emitter can be placed at any position as long as it does not intrude the convex hull of O. The size of an emitter should be large enough, since spherical harmonic lighting prefers area light sources and cannot handle point light sources.

As shown in Fig. 3.1, our framework consists of a preprocessing stage involving O and a run-time stage that involves both O and Gl. The preprocessing stage computes the transfer function for each vertex in O and build an axis-aligned bounding volume hierarchy for O. At runtime the incident radiance function is dynamically sampled and reconstructed. Thus the runtime stage can be further divided into two stages. The first stage is a breadth-first traversal of the BVH aiming to determine the number and the distribution of samples by means of the oracle described in section 3.3.2; and to sample the incident radiance function at the same time. Sampling

L(xi, ω) at a sample position xi is performed with hardware cubemaps and is described in Section 3.3.4. The second stage reconstructs the vector of spherical harmonic coefficients ~lp representing the incident radiance function at every vertex p in O by querying the octree data structure and then evaluates the reflection integral for each vertex using Equation (2.8), which is discussed in Section 3.3.5. Shading is calculated with spherical harmonic lighting. The final displaying is done with graphics hardware.

3.2

Preprocessing Stage

The preprocessing stage takes the scene’s geometry O as the input. No prior knowledge of emitters Gl is required. The result is a vector of spherical harmonic coefficients for transfer function at each vertex p in O, and a BVH of O. The BVH is organized as an binary tree flattened to an array. Each BVH node stores information for sampling and error estimation such as an averaged position of vertices, radius of bounding box, etc.

First we explain in Section 3.2.1 how to compute the transfer function and its spherical harmonic coefficient vector for each vertex, then in Section 3.2.2 we show how to build the

3.2 Preprocessing Stage 26

Figure 3.1: Framework Overview.

BVH.

3.2.1

Precomputing Radiance Transfer Coefficients

We precalculate the transfer function Tp(ω) for each vertex p of O and project Tp(ω) onto spherical harmonic basis functions as in [SKS02]. The result is a series of spherical harmonic coefficients, denoted as tipfor each spherical harmonic index i and each vertex p. We denote the vector formed by these coefficients as a function ~t(p) defined over all vertices or the function value ~tp at the vertex p. The transfer function is defined as the product of local surface BRDF ρ, the projected solid angle cosθ, and the visibility function V (p, ω) formed by O’s self occlusion.

Algorithm 3.1: Compute tipfor a vertex p

input : Vertex p and its normal Npand object geometry O

output: tipfor p Φ ←GenerateUniformDirectionSample(pnumberOfRay); 1 for i ← 0 to 15 do 2 tip ← 0; 3 foreach ω = (θ, φ) ∈ Φ do 4

convert ω into Cartesian coordinate ω0;

5

ray.origin ← p;

6

ray.direction ← ω0_;

7

vis←IntersectionTest(O,ray);

8 tip= tip+ Yi(ω)(Np· ω0) ·vis; 9 tip = 4tipπρp numberOfRay; 10

Now Equation (2.2) at each vertex p of O,

Lp(ωo) = Z Ω Lp(ω)V (p, ω)ρ(ωo, ω)cosθdω, can be rewritten as Lp(ωo) = Z Ω Lp(ω)Tp(ω)dω,

where Tp(ω) is equal to V (p, ω)ρ(ωo, ω)cosθ. Here we assume homogeneous materials across the object, and in our current implementation, the BRDFs are simply lambertian materials, which is a constant value, denoted as ρo, over the entire object. Note that the simple lambertian model can be further extended to arbitrary materials as in [KSS02]. The BRDFs should be wave-length dependent, but in our implementation we approximate them by tristimulus channels; so there will be Tpr, Tpg, Tpb, totally 3 functions to be projected on to spherical harmonic basis. Since the lamebertian model is simply a constant for each channel, we post-multiply the transfer functions with the material’s diffuse reflectivity. We compute V (p, ω) for each p in O by ray-casting 6400 rays in uniformly sampled directions on S2; see details in algorithm 3.2.

3.2 Preprocessing Stage 28

Algorithm 3.2: Generate uniform direction samples over S2

input : An integer√N , square root of the number of samples desired.

output: Φ, A set of directions (θi, φi).

/* ξ1 and ξ2 are uniform random variables. */

for i ← 0 to N do 1 φ ← 2πξ2; 2 θ ← arccos(1 − 2ξ1); 3 Φ ← (θ, φ) ∪ Φ; 4

functions Yi(ω) to yield coefficient tipin band i using the following equation

tip= Z S2 V (p, ω)ρocosθYi(ω) dω = 4π N X V (p, ω)ρocosθYi(ω) = 4πρo N X V (p, ω)cosθYi(ω),

where θ is the angle between the surface normal at p and the ray direction ω. The number of bands effects the accuracy of the reconstructed transfer function ˜Tp from tip. We use 16 bands (l = 4) for the approximation, where each band is represented by an IEEE 32-bit float-ing point number. Note that all tip are computed under the same coordinate system; thus no rotation is required. Some previous works compute the transfer functions under the coordinate system formed by the local tangent plane on the surface because the BRDFs are defined over the tangent plane [KSS02]. In such approaches, rotations of spherical harmonic coefficients are required. Our implementation computes the coefficients under the same world coordinate system, and sacrifices some flexibility for simplicity. Algorithm 3.1 shows this process. Inter-section tests for computing V (p, ω) are computed with an axis-aligned binary space partition tree for acceleration [PH04].

3.2.2

Bounding Volume Hierarchy Creation

The bounding volume hierarchy is a hierarchical structure of bounding volumes. Vertices in close spatial proximity are put in clusters and the clusters are enclosed in bounding volume.

There are three possible strategies to build a BVH [Eri05]: top-down, bottom-up, and in-sertion. Usually a bottom-up or insertion approach may result in a tree that is more balanced

than a top-down approach, but is much slower and complex to construct. We adopt a top-down approach because it is fast, simple, and less error-prone. The degree of children is also an issue. Higher degrees decreases the maximum depth of the hierarchy, but in the mean time increase the cost of hierarchy traversal. As stated in [Eri05], we find that the binary tree structure is a reasonable choice.

For example, since in practice our vertices number will not exceed 100000, the depth of a full binary tree will never exceed 18. Moreover, a binary tree can be flattened into an array without pointers for indexing; that is, for a node ni with array index i, the children of ni are

n2i+1and n2i+2. This is much more cache-friendly, since during a breadth-first search, children with the same parent node will be placed in nearby memory addresses.

The hierarchy is built in a top-down approach. We start from the bounding box of all vertices in O and build an axis-aligned bounding box for the entire object. We then find the longest axis of the bounding box and split the node into two children. The split point is the vertices’ coordinate median on the longest axis. This involves sorting the vertices along the longest axis. Splitting with vertices’ coordinate median gives us a balanced and near full binary tree, but not necessarily a BVH with optimal bounding volume [Eri05]. However, it still tends to produce a hierarchy with vertex clusters that are near to one another because a typical mesh for spherical harmonic lighting is finely and evenly tessellated. Since the vertices are distributed evenly over the object, split according to their coordinate median tends to minimize the bounding volume of nodes. Algorithm 3.3 shows this process. Fig. 3.2 shows the construction process of the BVH for the armadillo model. Note that some bounding boxes are quite large, since our strategy does not guarantee optimum sizes of bounding volumes. In our experiments, however, those large bounding boxes do not lead to visual artifacts.

Two attributes associated with each node are also precomputed during the creation of the hierarchy. One is the radius of the bounding volume of each node, which is a variable used in evaluating the oracle function to be discussed in section 3.3.2. The other is the sample position for the node, which is the position average of vertices inside the node. Since our bounding primitive is axis-aligned bounding box, the radius of the bounding volume for each node is set to be half of the diagonal length of the AABB.

3.3 Runtime Stage 30

(a) The model. (b) Vertices and AABB. (c) Split once.

(d) Split twice. (e) Split 12 times.

Figure 3.2: BVH construction for the armadillo model.

3.3

Runtime Stage

In this section we first explain how the bounding volume hierarchy and the error estimator are used to distribute samples, and then we discuss how to sample the incident radiance in Section 3.3.4. Finally how the octree is used for efficient reconstruction is described in Sec-tion 3.3.5.

3.3.1

The BVH and Its Use for Sampling

For a set of vertices P , the result of shading is the exiting radiance from a vertex p ∈ P to the camera. The color C(p) of the vertex p is the exiting radiance, computed a dot product of the two coefficient vectors for the incident radiance function and the transfer function, respectively, as shown in Equation (2.8). During run time, the incident radiance coefficient vector ~lp is unknown, while the coefficient vector for the transfer function ~tp is precomputed. If we replace

Algorithm 3.3: Recursively Build BVH for O

input : A BVH node and vertices V inside the node output: A binary axis-aligned bounding box tree BVH

node.AABB ←ComputeAABB(V);

1

node.avgPos ←AveragePosition(V);

2

node.d ←DiagonalLength(node.AABB);

3

laxis←GetLongestAxis(node.AABB);

4

(V1, V2) ←SplitAlongMedian(V ,laxis);

5

BuildBVH(V1,node.leftChild);

6

BuildBVH(V2,node.rightChild);

7

the true coefficient vector ~lp of the incident radiance function at every vertex p with a vector of coefficients ~lx sampled at some point x, the error introduced by such approximation ˜C(p) is

kC(p) − ˜C(p)k = k(~t(p) · ~l(p)) − (~t(p) · ~lx)k,

where the function ~t(p) is the vector function representing the coefficient vectors for the transfer function at p ∈ P , and in particular with L1 norm, we have

kC(p) − ˜C(p)k1 = X p∈P   (~tp· ~lp) − (~tp· ~lx)   , (3.1)

As shown in Fig. 3.3 for 1D cases, if we partition the domain P , which is R3 in our case and add more samples to the approximation, we get more accurate results. Previous works uniformly partition the domain using the ICP algorithm [SKS02] or uniform grid [GSHG98], here we can do it better. Since the error kC(p) − C0(p)k1 is defined on vertices rather than the entire R3, we can partition the set of all vertices P rather than uniformly partition R3. A useful partition on P is the bounding volume hierarchy.

We cannot apply Equation (3.1) to estimate the error of each hierarchy node, since it still involves computing ~l(p) at every vertex. Since it is impossible to know the exact value of

~l(p) unless all vertices p are sampled, our method tries to estimate them. The oracle can be

seen as a partition strategy. We don’t care about the exact values of error. Only the relative errors between different hierarchy nodes need to be computed. With relative errors, our method decides whether or not a node should be partitioned.

3.3 Runtime Stage 32

(a) An 1D function (b) Replace with an constant

(c) Partition the domain and adding samples

(d) As samples increase, the error diminishes

Figure 3.3: Approximating an 1D function with discreet samples.

3.3.2

The Oracle Function

Since the transfer function coefficients ~t(p) is the combination of the visibility function and the surface BRDF, the outgoing radiance will never be larger than the incident radiance, as a result of the energy conservation law. The errors in Equation (3.1) are thus never greater than

P P   |~lp| − |~lx|  

, which is further bounded by a looser bound as follows: max

p∈P (|~lp|) ||P ||, (3.2) where ||P || is the number of vertices of the set P . The intuition behind Equation (3.2) is that we will partition a vertex set that has a higher value of |~l(p)|, since discreet sampling introduces larger errors at this region. Note that the norm of vector of spherical harmonic coefficients representing the incident radiance is less than or equal to√i times the L2 norm of the incident radiance function. The magnitude of the incident radiance function at the sample point x can be obtained by computing the L2 norm of the incident radiance function:

kLx(ω)k2 =

s Z

S2

and the norm of the vector |~lx| is   ~lx   = v u u t X i Z S2 Lx(ω)Yi(ω) dω 2 .

By applying the Cauchy inequality to the norm of the vector |~lx|, we can see that it is related to the norm of the incident radiance function and have

  ~lx    2 =X i Z S2 Lx(ω)Yi(ω) dω 2 ≤X i Z S2 Lx(ω)2dω  Z S2 Yi(ω)2dω  = i Z S2 Lx(ω)2dω  = i kLx(ω)k 2 2, which implies |~lx| ≤ √ i kLx(ω)k2.

We can thus use the magnitude of the vector ~lx to adequately estimate the magnitude of the incident radiance function at point x. In addition to partitioning at brighter regions, vertex set P having a large number of vertices should also be partitioned first, since the error could be more appealing.

However, we must perform sampling to know |~l(p)| since ~l(p) is known only at run time. As a consequence, our oracle function is an a posteriori error estimator; that is, the error estimation is done after |~l(p)| is sampled. We made a few assumptions to further approximate the error.

One assumption is that ~l(p) does not vary too much around p so that we can take a single sample x inside set P and use it to approximate maxp∈P(|~l(p)|). Another assumption is that the points are evenly distributed on the object surface; so we can use the radius of the bounding volume of P to approximate ||P ||. Equation (3.2) can now be approximated by:

max

p∈P (|~lp|) ||P || ≈ r |~lx|, where r is the radius of the bounding volume of P .

Furthermore, we take the possibility of variations into account. The error is higher if there are more variations in ~l(p). We approximate the variation with the ratio of the radius of the bounding volume to the average distance to emitters. The intuition behind this choice is shown

3.3 Runtime Stage 34

Figure 3.4: If the blockers and the receiver are near, higher variations in irradiance oc-curs [Arv94].

in Fig. 3.4. If the emitters are far away from our receiver (object), the variation due to nearby blockers could be smaller, as the B curve in Fig. 3.4 shows. The blockers can be considered as a type of emitters with less emission. Also, the errors due to the variation are smaller if the receiver is very small compared to the distance to the emitters. One extreme case is a scene lit by infinitely far emitters; thereby no irradiance variations can occur. In the oracle, the variations are approximated by the ratio of the radius of the bounding box of the vertex set P to the harmonic mean of distances from the sample points to all emitters. The harmonic mean of a set of distances di, where i = 1...N , is [PH04],

Hd= N P N 1 di .

By using the harmonic mean we are able to detect the presence of nearby emitters, as a larger distance has a smaller effect on the harmonic mean. Now we can write our oracle as

r |~lx|

H2 d

, (3.3)

where ~lx is the vector of coefficients representing the incident radiance sampled at point x, r is the radius of the bounding box of the vertex set P computed in the preprocessing stage, and Hd is the harmonic mean of distances from x to emitters. The computation of ~lx and Hd will be described in section 3.3.

The oracle function in Equation (3.3) computes the estimated error between ˜C(p) and C(p).