The dipole diffusion approximation, which approximates the volumetric source distribution using a dipole (i.e. two point sources), was originally developed in medical physics community. The idea was proposed by Eason [Jensen et al. 2001] for modeling the back-scattering of light by blood. Farrell et al. [Farrell et al. 1992] used a single dipole to represent the incident source distribution for the noninvasive determination of tissue optical properties in vivo. Jensen et al. [Jensen et al. 2001]
then introduced the dipole diffusion approximation to computer graphics community for modeling the subsurface light transport.
The dipole diffusion approximation consists of positioning two point sources near the surface to approximate an incoming light (see Figure 2.3). One point source, the positive real light source, is locate at the distance zr beneath the surface, and the other one, the negative virtual light source, is located above the surface at a distance zv.
Figure 2.3: An incoming ray is transformed into a dipole source for the diffusion approximation [Jensen et al. 2001] [Poirier 2003].
By using the dipole diffusion approximation to solve diffusion equation, we can get the following expression for the radiant exitanceM (x0) at surface location x
xi 0
where is the albedo, describing the relative
importance of scattering versus absorption); σ = 3σaσt′is the effective transport
coefficient; sr = r2 +z2r is the distance from x to the positive real light source; 0 2
2 v
v r z
s = + is the distance from x to the negative virtual light source; 0
3 )
distance from the dipole source to the surface (shown in Figure 2.3). The mean-free path l
t
lu
σ′
= 1 is the average distance at which the light is scattered:
u . Finally, the
boundary condition for mismatched interfaces is taken into account by the A term which is computed as
dr
approximated from the relative index of refraction η by [Jensen 2001]:
η η
By using Equation 2.12, the subsurface illuminance then could be computed as:
∫
∈where Ê(x ) = Fi dt(η)E(x ) and E(x ) is the irradiance at point xi i i. The diffuse Fresnel
which is the diffuse BSSRDF defined as the ratio of radiant exitance to incident flux [Jensen et al. 2001]. Figure 2.4 shows the graph of Rd (which has exponentially decreasing property that we can employ in Chapter 3).
Finally, since the diffusion approximation already includes a diffuse Fresnel transmittance, the diffuse radiance L is computed as:
( )
η πwhere F is the Fresnel transmittance. t
Alternatively, we could omit the fresnel terms and assume a diffuse radiance:
π
Figure 2.4: The graph of Rd.
The dipole diffusion approximation introduced by Jensen et al. [Jensen et al.
2001] is widely used by the following researchers in the computer graphics
community. However, we find that this model could possibly be improved by using the equation proposed by [Kienle and Patterson 1997] from the medical physics community. They expressed the reflectance as the integral of the radiance over the backward hemisphere and considerably reduced the error in deriving the optical coefficients of translucent materials. By integrating their model into [Jensen et al.
2001], we expect that the measured parameters for translucent materials could be more accurate than [Jensen et al. 2001] and the final appearance could be more plausible as well.
Chapter 3
A Caching Technique for Rendering Translucent Materials
Our development of a caching technique for rendering translucent materials is based on the following observations:
Although the diffusion approximation with a hierarchical integration technique [Jensen et al. 2002] is a very effective way of approximating multiple scattering, it typically requires more than two hundred sample points per pixel to evaluate surbsurface illuminance for each pixel.
The subsurface illuminance tends to change slowly because the subsurface
scattering diffuses the incident light, blurs the effect of small geometric details on the surface, and softens the overall looks. The smooth appearance is distinct especially for high translucent materials and objects with small sizes.
The resulting subsurface illuminance value is view-independent because the subsurface scattering effect is only dependent on object geometry, object material properties, and how the object is illuminated.
It appears that for the sake of efficiency, we should not recalculate subsurface
previously computed values at nearby surfaces. The size of the small set of computed values should be independent of image size, thus high resolution images could be produced efficiently. Also, since subsurface illuminance does not depend on
viewpoints, these computed values could be reused for many images, which is useful for animation production.
In this chapter, we propose a caching technique for rendering translucent materials.
The framework of irradiance caching [Ward et al. 1988] is adopted. To do this, we derive an exact solution to calculate the gradient of subsurface illuminance and propose a split-disk model to determine the spacing of samples. Finally, we integrate the caching scheme into Jensen’s hierarchical evaluation method [Jensen et al. 2002].
3.1 Dipole Diffusion Approximation as a Convolution
Process
Recall that the subsurface illuminance function (Equation 2.14) is
A i BSSRDF expressed as:
( ) ( )
⎥′ is the effective transport
coefficient; sr = r2 +zr2 is the distance from p to the positive real light source; o
distance from the dipole source to the surface (shown in Figure 2.3). The mean-free
path
t
u σ′
l = 1
is the average distance at which the light is scattered. Finally, the boundary condition for mismatched interfaces is taken into account by the K term that is computed as
Fdr
1− Fdr
K +
=1
, where the diffuse Fresnel term Fdr is approximated from the relative index of refraction η by:
η
In the assumption of semi-infinite plane-parallel medium, becomes a function of distance between p and p only. By replacing parameter p and pi o i o with the
and expressing the vector parameter in terms of scalar values, we can rewrite Rd as follows:
Note that Rs is a three-dimensional radial function with its value decaying exponentially with the distance.
The original equation (Equation 3.1) integrates pi over the surface A. It seems that the integration is a two-dimensional process. But actually the integration is performed in three-dimensional space because pi is a three-dimensional point. So we change the integral domain and rewrite the Equation 3.1 in a form of three
dimensions:
. (3.4)
X Y Z
Substituting the Equation 3.3 into Equation 3.4 yields:
.
X Y Z
Because the R is a symmetric function, we can change the sign of the parameter: s
dxdydz
Obviously, the resultant equation is in the form of the following three-dimensional convolution of two functions:
Rs
E
S = ˆ ⊗ .
3.2 Calculating the Subsurface Illuminance Gradient
Given the subsurface illuminance function (Equation 3.1), it is not clear how to calculate the gradient of the subsurface illuminance. With the reformulated
convolution form in Equation 3.5, the gradient could be derived straightforwardly.
Recall that the derivative of a convolution function is:
dx
The gradient of the subsurface illuminance is then derived as follows:
) we could calculate the gradient of subsurface illuminace. However, since Ê does not have an analytic form, it is impossible to calculate analytically. Therefore, in our implementation, we choose Equation 3.7 to calculate the gradient of subsurface illuminance. Note that
∇Eˆ ∇Rs
⎥⎦
And the problem is reduced to evaluating the integral of convolution:
dxdydz
Although the integrals still do not have an analytic solution, by exploiting the
properties of (see Figure 3.1), we could use some integration techniques such as Monte Carlo and quadrature methods to get a good approximation.
Rs
Skimmilk B channel Skimmilk G channel Skimmilk R channel
x Rs
∂ Figure 3.1: The graph of∂ .
Notice that if we choose
to calculate the gradient of subsurface illuminance, it would be infeasible to get a good approximation even using the Monte Carlo method because the could not be easily sampled.
∇Eˆ
3.3 Applying the Gradient to Interpolation
Once we calculate the gradient of subsurface illuminance, we could use the gradient to interpolate the subsurface illuminance more accurately. We use the same weighted average as proposed in Ward’s caching technique [Ward et al. 1992] to interpolate the subsurface illuminance value:
[ ]
p is the position of the point to be computed, Pk is the position of cache k,
wk(p) is the weight of cache k with respect to p, C is the set of valid caches, {cache k: wk(p) > 1/a}, Sk is the computed subsurface illuminance of cache k,
Sk
∇ is the computed gradient of subsurface illuminace of cache k, and a is a user-specified error bound.
The next problem is how to determine the spacing of samples, i.e., how to
determine the weight of each sample, or how to estimate the error of each sample. The
simplest and theoretically most accurate solution is directly using the inner product of the offset (p− pk) and ∇Skderived in last section as our error estimateε
(assuming that the error due to interpolation is proportional to the estimated directional change of subsurface illuminance), i.e.,
S z p
p S y p S x p S
S x y z =∆ ⋅∇
∂
∆ ∂
∂ +
∆ ∂
∂ +
∆ ∂
=
∆
ε ∝ .
Unfortunately, this will lead to bias the calculation. Since gradient is a very local property, areas that just happen to have small subsurface illuminance gradient would be sampled at low density, even though there could still be sudden changes in the subsurface illuminance value due to nearby surfaces. A possible solution is to use some approximation models to capture the largest expected gradient in determining the sample density so that we could not miss anything relevant.
To estimate the largest expected gradient, we introduce a split-disk model analogous to the split-sphere model proposed by Ward et al. [Ward et al. 1988]. The split-disk model, based on the assumption that the geometry is locally flat, relates the subsurface illuminance gradient to the variance V of the irradiance values within nearby surfaces. It assumes that a surface element is located at the center of a disk which approximates nearby surfaces (see Figure 3.2). The radius of the disk, R, is heuristically determined according to the material scattering property. Half of the disk is totally bright with constant irradiance K and the other half is totally dark with constant irradiance of zero. Because the variance of the irradiance values within the disk is V, we can conclude that K=2V. The split disk has the largest expected gradient possible for surfaces with variance V.
An approximate bound to the change of subsurface illuminance in the split disk,ε , is given by the first order Taylor expansion of the function S of one variable:
u u S u
u ∂
− ∂
≤ ( )
)
( 0
ε ,
where uo is the center of the disk and u is some other point on the disk. Note u and uo
areboth one-dimensional value because we only care about the distance between two points.
Sector of circle A
Figure 3.2: The split-disk model. A surface element is located at the center of a half-dark disk.
To derive
we firstly consider a surface element moving from uo to u, the change of S could be computed as twice of an integral over sector of circle A plus an integral over right triangle B (see Figure 3.2):
)
Unfortunately, we cannot find an analytic solution of the integral describing
subsurface scattering over right triangle ∆SB. Inspired by [Mertens et al. 2003] where a semi-analytical integration method is derived to solve the integral describing
subsurface scattering over an arbitrary triangle, we can approximate ∆S by an integral
∆u R
θR
Right triangle B
E = 2V E = 0
uo u
over four sectors of circle, i.e.,
(r) (Equation 3.2) yields Then replacing Rd
( ) ( )
By rearranging terms, we can obtain
(
σ) (
σ)
θ2
Finally, we arrive at
⎟⎟⎠
and the error estimate can be computed as
(
0)
We can extend our approximation to more complicated geometries by replacing u with vector-derived values:
⎟⎟⎠
As in [Ward et al. 1988], the inverse of the error estimate
⎟⎟⎠
⎜⎜ ⎞
⎝
⎛ − + −
− ′ − r − r − v − Rv
v z v R r z r i
i i
i
R e e z
R e e z
V R p
p p σ σ σ σ
π ε ( ) = 2α
= p
w 1 1
)
( (3.10)
is then used as our weight.
Substituting Equation 3.10 into Equation 3.9, the subsurface illuminance of some point of interest then can be computed by interpolating nearby caches. Additionally, the interpretation of the subsurface illuminance gradient allows us to perform the calculation separately for each sample wavelength or to combine into a spectral average. The latter was chosen for our implementation to save storage costs and evaluation time, although this may cause some blurring artifacts when the difference of scattering properties within each sample wavelength is too large.
3.4 A Three-Pass Technique for Rendering Translucent
Materials
To integrate our model into Jensen’s hierarchical evaluation method [Jensen et al.
2002], we use a three-pass approach, in which the first pass consists of computing the irradiance at selected points on the surface, the second pass generates all the necessary cache samples whose values including subsurface illuminance, gradient of subsurface illuminance, and variance of irradiance over nearby surface are computed by using the precomputed irradiance values, and the last pass re-uses the caches to produce the final image via interpolation.
Note that the second pass in our three-pass approach only generates caches. It doesn’t use caches to interpolate any value. The interpolation using caches is done in the third pass. The reason why we do not combine generating caches and re-using caches into a single pass is discussed in section 3.5.
Pass 1: Sampling Irradiance
To solve Equation 3.1, firstly, we need to sample the irradiance function E(x). There are a number of methods to generate sampling positions on the surface. In [Jensen et al. 2002], Turk’s point repulsion algorithm [Turk 1992] is used to obtain a uniform sampling of points on a polygon mesh. However, a uniform sampling seems irrelevant in their hierarchical approach as each sample point is weighted by the area associated with it. Instead of using the Turk’s point repulsion algorithm, which is complicated to implement, we use a very simple method to obtain the sampling positions on the polygon mesh. We directly use the centroid of each face as our sample point and assign the area of the face as the area associated with this sample point. If a model is too coarse and results in low-frequency noise in final image, we subdivide the model until the noise disappears.
For each sample point, we store the position, the area associated with the point, and a computed irradiance estimate. Since we focus on caching technique in this thesis, we do not use any rendering technique that accounts for global illumination (such as photon mapping and distributed ray racing) to compute the irradiance. We simply sum up irradiance contributions from each light source for evaluating direct illumination on each sample point.
As stated in [Jensen et al. 2002], the irradiance samples described in above section should be stored in a hierarchical structure so that,by clustering distant samples, we can exploit the exponential shaped fall-off of Rd, thus facilitating the evaluation of diffusion approximation. There are a number of hierarchical structures that we can use to store irradiance samples. Here we choose an octree as proposed by Jensen et al. [Jensen et al. 2001] in our implementation. Each node in the octree stores some information representing all irradiance samples inside the voxel associated with the node: the total flux Φ, the total area A, and the average position P. Note that the leaf node could have up to 8 irradiance samples for efficiency issue. The Φ, A, and P are computed as follows:
∑ ∑
=
j j j
i A
P P A
∑
Φ=
Φi j , Ai =
∑
Aj , ,where node j is node i’s child and if node i is a leaf node,
where irradiance sample k is within node i.
Furthermore, because our split-disk model needs to compute the variance of irradiance samples over nearby surface and the variance is derived as
,
, (3.11)
where I is the irradiance distribution over nearby surface.
E[I]can be computed as
∑ ∑
where is the irradiance of nearby sample i scaled by Fresnel term, and Ai is associated area of nearby sample i.
For the sake of E[I2], which is computed as
we store additional information M in each node: i
∑
where irradiance sample k is within node i.
Pass 2: Generating All the Necessary Caches
Before we state how to find where all the necessary caches should be generated, we present how to compute the values stored in each cache. The values we stored in each cache are subsurface illuminance S, gradient of subsurface illuminance , and variance V of irradiance over nearby surface. All these values are computed by using
∇S
the precomputed irradiance values (distributed in Pass 1) stored in a hierarchical structure.
To compute the subsurface illuminance and the gradient of subsurface illuminance, we need to solve the integral in Equation 3.1 and Equation 3.8,
respectively. If we directly sum the contribution from all the irradiance samples, i.e., solve the integral in Equation 3.1 and Equation 3.8 using uniform sampling, the computation would be too costly. Thus Jensen et al. proposed to use a hierarchical integration technique [Jensen et al. 2002], i.e., find a set of points N (see Figure 3.3) that could effectively represent all the irradiance samples to compute the integral.
Pre-computed irradiance samples (distributed in Pass 1)
Representative points of a cluster of pre-computed irradiance samples An interested point
The set N
Figure 3.3: The set N.
Finding the set N
To find the set N with respect to a point p, we follow the algorithm proposed by Jensen et al. [Jensen et al. 2002], and begin by traversing the octree from the root. For
each node, we first check whether the node is a leaf. If it is, we include all the
irradiance samples within it to set N; otherwise, we then check whether p is inside the voxel associated with the node. If p is inside the voxel, we continue to traverse down to the children; if not, we then compare the solid angle with respect to the voxel with a user-specified value m, which controls the error. If the solid angle is larger than m, we traverse down to the children; otherwise, we add the representative point of the node to the set N. The pseudo code is listed below:
Traverse_octree(location p, node i) {
IF node i is a leaf node
Add all irradiance samples within voxel i to set N ELSE IF point p is inside voxel i
FOR each child j of node i Traverse_octree(p, j)
ELSE IF the solid angle of voxel i > m FOR each child j of node i
Traverse_octree(p, j) ELSE
Add representative point of node i to set N }
After we find the set N with respect to a point p, the subsurface illuminance at p is computed by summing up the contribution from point i in N:
∫
As for evaluating the gradient of subsurface illuminance at p, ideally, we should use some integration technique exploiting the properties of ∇ ( Figure 3.1) to solve Rs Equation 3.8. Also, the solution to the integral should be evaluated very fast, or the cost for computing the gradient will cancel out the gain from the interpolation and extrapolation. Fortunately, we don’t need to devise some sophisticated integration technique to solve the integral. We could directly use the set N as our sample point to hierarchically solve the integral, just as we use the set N to solve Equation 3.1, though the plots of Rs (Figure 2.4) and (Figure 3.1) are not quite the same. The
gradient of subsurface illuminance could be computed as follows:
Rs
Note that intuitively, it seems that we should use vector Pi - p as our parameter to calculate the gradient, but we use vector p - Pi instead. This is due to the sign change during the derivation of the convolution process (Equation 3.5).
Another value we have to compute is the variance V of irradiance over nearby surface (see Equation 3.11, 3.12, and 3.13):
,
where I is the irradiance samples within nearby surface. Note we again use the set N to calculate the variance.
Determining where the caches should be generated
To determine where all the necessary caches should be generated, we firstly use ray casting to find a set of visible points X. For each point xi in X, we check if there is any previously computed cache at nearby surface that could be used for interpolation, i.e., any cache k with Wk(x ) > 1/a. If any, we leave xi i to next pass; otherwise, we generate a new cache at point x , evaluate S(x ) and i i ∇ associated with the cache. Note that checking each point x in X in different order could change the i
resulting distribution of the caches. Here we choose bottom-up scan-line order in our implementation.
As stated is [Ward et al. 1988], each previously computed cache is only valid for interpolation in some finite space. The “valid domain” of each cache is computed as follows:
which implies
⎟⎟⎠
By using this “valid domain”, we could store these computed caches in an efficient data structure so that nearby caches can be located fast. As in [Ward et al.
1988], we choose an octree to store our computed caches. We begin by using the object’s bound box as the root of the octree. When a new cache is generated at some location, the octree is subdivided on demand to contain this value. Each cache value,
1988], we choose an octree to store our computed caches. We begin by using the object’s bound box as the root of the octree. When a new cache is generated at some location, the octree is subdivided on demand to contain this value. Each cache value,