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
24
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
convert ω into Cartesian coordinate ω0;
5
Now Equation (2.2) at each vertex p of O, Lp(ωo) =
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.
We project the transfer function Tp(ω) at vertex p onto a set of bands i of spherical harmonic
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).
functions Yi(ω) to yield coefficient tipin band i using the following equation tip=
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);
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 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
, which is further bounded by a looser bound as follows:
maxp∈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
Lx(ω)2dω ,
and the norm of the vector |~lx| is
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
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:
maxp∈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|
Hd2 , (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).
Once the estimated error exceeds a user-specified bound, we must traverse to the next level and
Algorithm 3.4: BVH Traversal
input: Desired error bound ε and sample limit N Insert root node intoqueue;
1
evaluate the samples at two children nodes.
3.3.3 Hierarchy Traversal
In order to determine the number and the distribution of samples, the bounding volume hier-archy is traversed; and samples are then positioned at proper nodes. Starting from the root, we compute the incident radiance at the sample position precomputed in the preprocessing stage.The oracle function is then evaluated to check whether the sample is adequate for recon-struction. If the oracle exceeds the user-specified error bound, two children of this node are recursively traversed in order to distribute more samples; otherwise, the recursive traversal ter-minates. Shown in algorithm 3.4, this is performed with a queue. The traversal stops when the oracles evaluated at all nodes being traversed are below the error bound or the number of samples exceeds some user-specified limit. During the traversal, the samples are inserted into an octree for efficient reconstruction, which will be discussed in section 3.3.5. Each sample is associated with a valid domain, aiming to prevent the global reconstruction results from being influenced by local, high-intensity samples. The valid domain is defined as a sphere with radius r centered at the sample’s position.
3.3 Runtime Stage 36
3.3.4 Incident Radiance Sampling
To dynamically sample the incident radiance function at a sample point x, we exploit the use of graphics hardware. We also compute the averaged distance from x to emitters for error estimation at the same time. Emitters Gl are sent down to graphics pipeline six times to render a cubemap. The cubemap now represents the emission of emitters to x. We then project the incident radiance represented by the cubemap onto the spherical harmonic bases. Because cubemapping is formed by perspective projections onto six planes, during the integration we must take into account the non-uniform solid angle associated with each pixel. Let ω(c) be the solid angle of a pixel c(s, t) of a k × k cubemap face with viewing direction R. The relation of the cubemap faces and viewing directions R and their corresponding up vectors T is shown in Table 3.1. The vector ~c from the sample point x to the pixel c(s, t) is
~c = ( s
2k − 1 + 1
k)S + ( t
2k − 1 + 1
k)T + R,
where S is the cross product of R and T [Ope99]. The coefficients li representing the incident radiance at x is computed by the following equation
li =
where ω(c) is the solid angle associated with pixel c(s, t). We use Equation (2.1) to compute the solid angel ω(c) of a pixel c from the area of the pixel:
ω(c) = X
The sample position is the average position of all vertices inside the node and is precom-puted. The number of coefficients for each vector of coefficients ~lx for sample x is also 16.
After the sampling and coefficient computation are done, the vector of coefficients is then used to estimate the error of the node.To estimate the error, we also need the average distance from
R up vector (T )
+X −Y
−X −Y
+Y +Y
−Y −Y
+Z −Y
−Z −Y
Table 3.1: Coordinate system of each cubemap face [Ope99].
the sample to all the surrounding emitters. As we render the emitters into the cubemap, we also render their distances to the sample into another color buffer. In our implementation this is done with OpenGL’s multiple draw buffers and shaders. Then we read the color buffer back to main memory and compute the harmonic mean of distances.
3.3.5 Efficient Reconstruction using Octree
To approximate true incident radiance coefficient vector ~lp0 from a given set of samples {~lx}, we
To approximate true incident radiance coefficient vector ~lp0 from a given set of samples {~lx}, we