Chapter 2 Related Work
2.2 BRDF Importance Sampling
BRDF importance sampling is a technique to reduce the image variance in physically-based rendering. Its concept is to find the distribution based on the representation of BRDF. Shirley demonstrated how to sample the traditional Phong BRDF model efficiently [19], Lafortune also presented importance sampling schemes for the modified Phong model [11]. Ward [22] showed the stochastic sampling method for the BRDF models composed of (elliptical) Gaussian kernels. For more examples, please see Pharr and Humphreys [16].
Lafortune [10] used multiple cosine-lobes for representing the BRDF, he used non-linear fitting algorithm to fit sums of cosine-lobes to an analytical model or to actual measurements. Though this representation is simple and can be applied for Monte Carlo importance sampling efficiently, it is hard to approximate the complex BRDF by using his fitting process.
Lalonde [12] used wavelets to represent the BRDF and proposed an importance sampling scheme by binary searching the tree constructed against the coefficients of each wavelet basis. Matusik [14] also used a wavelets representation of BRDF, and he presented a numerical sampling method based on wavelets analysis.
Lawrence et al. [13] demonstrated an importance sampling method based on a factored representation. They reparameterized the BRDF by using half-angle [18], then used non-negative matrix factorization (NMF) twice to decompose the BRDF data for efficient importance sampling.
Our approach represents the BRDF using scattered SRBFs. Tsai and Shih [20]
have shown this representation performs better than wavelets in data analysis.
Furthermore, we can directly fit scattered SRBFs to the BRDF data without reparameterization. We also propose an efficient importance sampling method based on this representation. In the following chapters, we will describe how we represent the measured BRDF data with scattered SRBFs and how we estimate the probability distribution for importance sampling.
Chapter 3
Spherical Radial Basis Functions
Spherical radial basis function (SRBF) is a special radial basis function (RBF) defined on the unit sphere. SRBFs avoid distortions and the artificial boundaries when they are used to analyze the data on sphere. In this chapter, we first introduce the background of SRBFs, and express the useful properties of SRBFs that had been developed by previous researchers. Then we discuss the scattered SRBFs briefly, and point out the advantages of scattered SRBFs.
3.1 Background of SRBFs
SRBF is recognized as an axis-symmetric reproducing kernel function defined on Sm, the unit sphere embedded in \m+1. It can be taken as a special radial basis function (RBF) defined on the unit sphere, i.e. the kernel functions only depend on the spherical distance between unit vectors (Figure 3.1).
Let η and ξ be two points on Sm, and let θ η ξ( , ) be the geodesic distance between η and ξ on , i.e. the arc length of the great circle joining the η and ξ.
Because kernel functions of SRBF are depending on θ, and SRBF can be expressed in the expansions of Legendre polynomials:
Sm
0 each Legendre polynomial satisfy the following conditions:
Gl
Figure 3.1: 3D plot of a Gaussian SRBF
When all Gl are positive, a spherical function F( )η can be represented in SRBF
Since SRBF has expansions of the Legendre polynomials, there is a very useful property based on the orthogonal property of Legendre polynomials in [-1, 1] called spherical singular integral order-l spherical harmonics on , and
Sm dm l,
Sm dω denotes the differential surface element on Sm.
The spherical singular integral could be transformed into a simpler form in some conditions. Tsai and Shih [20] have proved that the convolution of two Abel-Possion SRBFs kernels, or two Gaussian SRBFs kernels has a mathematically simple form with small m. For example, the definition of Gaussian SRBF kernel is
( )
( ; ) , 0
GGau η ξ λ⋅ =e e−λ λ η ξ⋅ λ> , (5) where λ, the parameter called bandwidth, describes the coverage of the SRBF, and the convolution of two Gaussian SRBFs can be written as
( ) 21 to [20], and there are two references for further about SRBF [15] [6].
3.2 Scattered SRBFs
Distribution of the SRBFs’ centers on the sphere affects the compression efficiency significantly. If we use uniform SRBFs to represent the data with sparse distribution, it would waste lots of basis kernels on the region without any data. On the other hand, using scattered SRBFs, i.e. adapting the center, bandwidth and coefficient of each basis, we can locate the SRBF kernels based on the data distribution on the sphere. Therefore, scattered SRBFs can capture the feature of the original data with much fewer bases than those used in uniform SRBFs. Nevertheless, it may be hard to compute the parameters of each scattered SRBF (i.e. center, bandwidth, and coefficient) from original data. Fortunately, Tsai and Shih [20] present a novel approach to fit data with scattered SRBFs recently. Their approach makes scattered SRBFs more practical for representing data on sphere.
Chapter 4
Sampling Scheme
Our system consists of two major processes, one is the off-line process and the other is the run-time rendering process. In the off-line process, we transform the BRDF measurements into the scattered SRBF representation. In the run-time rendering process, we first determine how many samples should be taken from each SRBF. Then we generate samples based on the probability density function (PDF) calculated from each SRBF. Finally, we combine the sampling results with multiple importance sampling technique presented by Veach and Guibas [21].
The following sections are organized as follows. In Section 4.1, we will introduce the non-uniform and non-negative SRBF fitting algorithm which we use to transform the measured BRDF data into scattered SRBFs, and discuss the initial guess applied for BRDF data. In Section 4.2, we first briefly introduce the multiple importance sampling technique. Then, we describe how to determine the number of samples of each SRBF kernel and how to generate samples from each SRBF kernel.
4.1 Off-line Process
In this section, we would show how to represent the measured BRDF data by using scattering SRBFs. We first introduce the fitting algorithm, and then we describe the initial guess for BRDF data in the fitting process.
4.1.1 Non-uniform and non-negative SRBF fitting algorithm
We now describe how we transform the original BRDF measurements into the scattered SRBFs representation. The measured BRDF data are usually represented in tabular form, row is the outgoing (viewing) direction, and column is the incident (lighting) direction. When we fix one outgoing direction, we can get a vector which describes the ratio of energy from all incident directions (Figure 4.1). The order of each element contained in the vector is the same as the pixel order defined in the cubic texture, i.e. first unfold each face’s pixels into a vector, then concatenate the unfolded vectors of each face in the order +x, -x, +y, -y, +z, -z.
Figure 4.1: The format of measured BRDF data.
In our fitting process, we fit the scattered SRBFs to the data in a vector of each outgoing direction respectively, and there are three sets of parameters we want to optimize: the set of coefficients L, the set of centers Ξ , and the set of bandwidth 2 parameters Λ2. Our objective is to minimize the square error with the original data:
2
Tsai and Shih [20] have presented a novel method to solve this objective function when modeling the lighting environment. In their fitting approach, the coefficients of each kernel basis are calculated by ordinary least-squares (OLS) or regularized least-squares (RLS). However, fitting by OLS or RLS would make the coefficients become negative in some situations. Because we need to use the coefficients to estimate the probability distribution, we must constrain all coefficients to be positive. Therefore, we also use L-BFGS-B solver [23] to fit the coefficients instead.
The fitting process composes of three main steps:
1. Using L-BFGS-B solver to optimize the set of centers with a given initial guess or the temporary results of the previous iteration.
2. Then, using L-BFGS-B to solve the set of bandwidth parameters and the set of coefficients in turn.
3. If the squared difference errors between current and previous iteration are less than a threshold, the process is terminated. If the count of iterations exceeds a user-defined threshold, break the process as well.
This fitting process is one kind of non-linear optimization. The fitting results are highly dependent on the initial guess. In other words, the initial guess would dominate the accuracy of the representation. In next subsection, we would introduce the initial guess for BRDF data in more details.
4.1.2 Initial Guess for Non-linear Optimization of BRDF data
As mentioned in last subsection, the initial guess plays an important role in non-linear optimization. Tsai and Shih [20] have studied this problem recently, and they present a reasonable guess based on some heuristics to decrease approximation errors. The principle of their initial guess is to make the initial centers of each kernel basis locate at all local peaks of the original data as close as possible, and to estimate a reasonable bandwidth according to the coverage of each initial center. Based on this principle, they generate lots of candidate guesses sorted by the coverage-weighted intensity, and then attempt to choose guesses which avoid the initial centers gathering in local area. This guessing scheme works well for the data with high dynamic range (HDR) property. While for most BRDFs data, the magnitude of the gradient of each data element is much smaller than the HDR image. In order to apply this guessing schema for BRDF data to achieve the same performance in HDR image, we only need to do a minor change – emphasizing the weighting of intensity. We do this straight by squaring the influence of intensity. This modification makes this initial guess schema much tend to select the candidates with high intensity in the entire region of BRDF.
Figure 4.2: Overview of initial guess generation.
(a) Visualization of measured BRDF data on unit sphere.
(b) Estimate the coverage of each unit direction on the sphere.
(c) The candidates selected so far and the summation of their coverage (i.e. the area colored with light orange).
(d) The candidate falls within the range (c) would be rejected (magenta one), yet would be accepted (yellow one).
Figure 4.2 gives the overview of the initial guess. The original measurements of BRDF are discrete. After mapping the measured data onto the sphere, we can get the visualization result (a). Then we estimate the coverage of each direction in a dense set
of unit directions, i.e. each data element on the sphere, as illustrated in (b). With each unit direction η, we want to estimate the coverage by calculating the squared difference of the intensity between η and its neighbors, i.e. I( )η −I( )ηi 2. The coverage of η would spread out from η until the squared difference exceeds a user-defined threshold. Next, we generate a priority queue by sorting the coverage-weighted square of the intensity of each unit direction. Then we choose the appropriate candidate directions from this priority queue. The choosing rule here is intuitive. If the direction of the considered candidate falls within the range covered by the candidates that we previously select, we reject this candidate, and take next candidate from the priority queue into consideration. On the contrary, we accept the candidate to be one of the set of our initial guess. Figure 5.2 (c) illustrates the circumstance we have chosen two candidates, and (d) describes the accepted case (the yellow one) and rejected case (the magenta one) by our choosing rule. This process would continue until the number of the selected candidates exceeds a user-defined threshold. If the priority queue becomes empty before the process terminate, we would narrow the coverage threshold previously defined by user, and restart the process from step (b).
4.1.3 BRDF Representation using Scattered SRBFs
After the fitting process, we have represented the BRDF data with scattered SRBFs for each fixed outgoing direction:
K
kernel function, ξk is the center of basis on unit sphere, λk is the bandwidth of the basis, Fk is the basis coefficient, and K is the number of the bases.
Given an arbitrary viewing direction in 3D space, we need to compute the SRBFs composed of the SRBFs of neighbor outgoing directions. Then, the BRDF is given as: where N is the number of neighbor directions. N would vary according to the given viewing direction. Its value would be 3 or 4 in usual. Wn is the weight for each SRBF expansion. The weight Wn is inversely proportion to the geodesic distance between the given direction and the neighbor outgoing direction. Using this representation we can easily generate samples for rendering with the multiple importance sampling technique. This will be discussed further in next section.
4.2 Run-Time Rendering Process
Importance sampling is a technique to reduce the image variance. Now we have represented the BRDF with the scattered SRBFs. Traditional importance sampling is to find a probability density function (PDF) which describes the density distribution of the entire data region, i.e. the hemisphere for BRDF. However, in this case, it’s hard to combine the PDFs into a new one that describes the distribution of the sum of all related scattered SRBFs. On the other hand, we can’t directly use the distribution function derived from each SRBF neither. It’s because that each SRBF may overlap each other. It’s only correct when all SRBFs are disjoint, and this will become another sampling technique named stratified sampling. Therefore, we apply an alternative
approach presented by Veach and Guibas [21] - multiple importance sampling, which allows all sampling strata to overlap each other. Since each SRBF kernel can be taken as some kind of the density distribution function. Thus we can easily generate samples from each SRBF kernel. With the aid of multiple importance sampling, we can make the rendering process more efficient and compact when we use scattered SRBFs to represent BRDF.
Figure 4.3: Run-Time Rendering Process
Figure 4.3 illustrates the flow of our rendering process. When the view ray hits the object in a scene, we first weight sum the SRBFs of three or four near fixed view direction that we construct in off-line process. Next, we determine the number of
samples of each SRBF according to its integral. Then we generate samples from each SRBF and combine the results by multiple importance sampling technique.
In this section, we start by introducing the multiple importance sampling. We then describe how we determine the number of samples that should be taken from each SRBF, and how we generate sampling direction from each SRBF.
4.2.1 Multiple Importance Sampling
When tracing a path through the scene by BRDF importance sampling, it is desirable to generate rays distributed according to the density of BRDF. We are interested in evaluating the integral of the incident illumination for a fixed outgoing direction ωo located at x with normal n,
where the first line is the reflection equation, and the second line is a Monte Carlo estimator for importance sampling. However, it’s hard to construct a single PDF
( s | o)
γ ω ω that follows the shape of the complex BRDF. There is an alternative technique for importance sampling presented by Veach and Guibas [21] called multiple importance sampling. This technique makes Monte Carlo integration more robust by combining several potentially good estimators. These estimators computed by different PDFs may have different qualities in different regions of the integration domain. Veach and Guibas make a weighted-average of all estimators where the weights depend on the sampling positions. If we want to evaluate the integral of f(x)
( ) f x dx
Ω
∫
,and we have n different estimators, the combined estimator is then given by independent. wi is the weighting function and satisfy the following two conditions:
1
Then, the expected value of the combined estimator F would be equal to the integral of f(x) which we want to evaluate.
So far, we have represented the BRDF using scattered SRBFs. Now we apply this technique to BRDF importance sampling, and the Equation (10) would become:
Ω where n is the number of SRBFs, denotes the number of samples from each SRBF, and
ni
p is the PDF calculated from each SRBF kernel, for i i=1,...,n. are the samples from each SRBF kernel, for
,
Xi j
1,..., i j= n .
Now, there are two issues: one is how to determine , the number of samples should be taken from each SRBF kernel, the other is how to generate , the sampling directions distributed according to each SRBF kernel. These will be
ni
,
Xi j
discussed in the following subsection.
4.2.2 Sampling Algorithm
We now describe how to use our representation for multiple importance sampling. Intuitively, each scattered SRBF covers a part of the entire BRDF region.
Although there would be two or more SRBFs overlap in some regions, we can still gather the intensity of multiple sampling directions within a pixel using multiple importance sampling approach.
As mentioned before, we fit scattered SRBFs to the data in the vector of each discrete outgoing direction respectively. At run time, for a given viewing direction, we need to combine the scattered SRBFs of three or four nearest viewing directions. Then, we have multiple SRBF kernel functions, now we should decide how many samples should be taken from each SRBF kernel. Let’s recall that SRBF is defined on the unit sphere, and its integral is easy to be calculated by using the spherical singular integral property of SRBF. The integral of each SRBF can be taken as its total energy gathered from all directions. Therefore, it’s straightforward to allocate samples according to the ratio of the integral of each SRBF to the sum of all SRBFs’ integrals.
The probability of choosing the SRBF l for each sample is given by:
1
We calculate a 1D CDF over l from these probabilities. In the dispatching process, we generate a uniform random variable in [0, 1] for each sample, and traverse the CDF to determine where the sample should be taken from.
Next, we generate a sample direction in each SRBF kernel by sequentially selecting the elevation angle θ and the azimuth angle φ defined on the plane whose normal is parallel with the SRBF kernel’s center (Figure 4.4).
Figure 4.4: The elevation angle θ and the azimuth angle φ defined against SRBF.
Having chosen the SRBF basis to sample, we can consider its PDF proportion to its kernel function. With the kernels easy for integration and inversion, such as Gaussian kernel, we can use the inversion method to generate random samples directly. We first normalize the kernel function to get a PDF ( )p x , since the integral of PDFs over their domain should equal to 1. After integrating ( )p x , we can get the cumulative distribution function (CDF) . Then, we calculate the inverse of CDF
, and generate a uniformly distributed random number ( )
P x
1( )
P− x ξ in [0, 1] to evaluate
the P−1( )x . We can easily get the elevation angle θ of the sampling direction.
On the other hand, for the kernel which is hard to apply inversion method, we can use metropolis random walk algorithm to generate samples with a desired density instead. It should be noted that the SRBF kernel itself only describes 1D density, while the desired density for sampling is the density over the sphere. Therefore, when we select the elevation angle θ by this approach, we should consider the influence of the circumference around the center of SRBF kernel especially.
The azimuth angle φ defined in each SRBF is easy to be calculated. Because each SRBF kernel is symmetric against the vector that defined by its center and the origin of the unit sphere. Therefore, the distribution of azimuth angle φ is uniform, and we can simply generate a random number in [0, 2π] to evaluate φ. With SRBF
The azimuth angle φ defined in each SRBF is easy to be calculated. Because each SRBF kernel is symmetric against the vector that defined by its center and the origin of the unit sphere. Therefore, the distribution of azimuth angle φ is uniform, and we can simply generate a random number in [0, 2π] to evaluate φ. With SRBF