BSSRDF Importance Sampling 4 - multiple axis probe ray sampling

First of all, Big big thanks to Chris Kulla for borrowing the slide pictures and the explanation on the way he calculated the max probe radius in paper! This is by far the most fun part in this project spring.

Yet another look on BSSRDF integration equation......
$$L( x_{o},\vec{\omega _{o}} )=\int_{A }\int_{\Omega }L_{i}( x_{i},\vec{\omega _{i}} )S (  x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}} )cos\left ( \theta \right )d\omega dA  $$
Let the inside integration be \(F(x_{i}) = \int_{\Omega }L_{i}( x_{i},\vec{\omega _{i}} )S (  x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}} )cos\left ( \theta \right )d\omega \). This part should be pretty similar to BRDF calculation: we sample the light, calculate the irradiance, apply BSSRDF and \(cos(\theta)\)...Then we integrate this \(F(x_{i})\) over the surface area. It should look like this: \(\int_{A }F(x_{i})dA\) , and its Monte Carlo Estimator will look like this:
$$\frac{1}{N}\frac{F(X_{i})}{pdf_{A}(X_{i})}$$
We have the 2D sample technique and its corresponding \(pdf_{disc}\) in hand now, how we convert it into \(pdf_{A}\) ? If the surface is a flat surface and we project the disc sample in its normal direction (towdown) , then this two pdf should be the same. If we encounter the curvy surface, then \(pdf_{disc} = pdf_{A_{\perp}}\)  . Since \(A_{\perp} = A\parallel\vec{n_{x_{o}}}\cdot\vec{n_{x_{i}}}\parallel\) , and \(pdf_{A_{\perp}} : pdf_{A} = \frac{1}{A_{\perp}} : \frac{1}{A}\) , we have \(pdf_{A} = pdf_{A_{\perp}}\parallel\vec{n_{x_{o}}}\cdot\vec{n_{x_{i}}}\parallel = pdf_{disc}\parallel\vec{n_{x_{o}}}\cdot\vec{n_{x_{i}}}\parallel \) . This can be interpreted in a more intuitive way: we sample by projecting down based on \(n_{x_{o}}\),  and if \(n_{x_{i}}\) is almost perpendicular to  \(n_{x_{o}}\), then chance it got sampled is almost 0 since \(\parallel\vec{n_{x_{o}}}\cdot\vec{n_{x_{i}}}\parallel\) approaches 0 . The following image illustrates the idea:

image taken from Kulla 13, the right most red sample point will have a relative low area pdf since its surface normal is almost perpendicular to the disc sample center \(x_{o}\)

uhoh......In this case, if the nearby surface of \(x_{o}\) all have really different normal direction from \(n_{x_{o}}\), even the top down projection disc sampling make flat surface sampling ideal, we would get bunch of sample with really small \(pdf_{A}\) since lots of the sample may have \(\parallel\vec{n_{x_{o}}}\cdot\vec{n_{x_{i}}}\parallel\) approaches 0, and this can be source of variance noise! It makes sense since highly curvy geometry info means high frequency data in local area, and it does need more sample to converge. Besides, one of dipole diffusion approximation's assumptions is that the surface should be semi-infinite flat.

BSSRDF Importance Sampling (I'll just call it Kulla 13 or BIS later) proposed a simple and elegant solution for the above issue: instead of only doing the disc sampling based on normal direction projection, they also do the sampling alone the orthogonal axis (the slides mentioned that the orthogonal axis can be arbitrary, and their experiment show shading tangents work pretty well). After taking samples from three different axis (normal, dpdu, dpdv), they further weight sum the samples with multiple importance sampling technique so that the sample taken from normal projection but falling on almost vertical neighbors (and also, the sample taken from tangent projection but falling on almost horizontal neighbors) never get high contribution, and the sample that is appropriate for corresponding case can get relatively higher weight. Below two screen shots taken from Kulla 13 illustrate the idea pretty clear:

For further optimization, the talk also suggests that by specifying a max probe radius \(R_{max}\), we can limit the sample probe ray length to just cover nearby region. so the intersection test can be terminated faster. My understanding about this part, is that we sacrifice part of the accuracy for efficiency: since we only casting sample within 0-\({R_{max}}\) range, there will be surface point with non zero irradiance got assigned pdf 0, and it means we introduce bias in this case. The bias probably won't matter that much as long as the pdf integration got ignored is relatively small. After all, the BSSRDF we have is already an approximation model derived from semi-infinite flat surface, while most of the geometry we tried to render don't satisfy this assumption. 

And now we try to get these probe rays, first thing to do is to figure out the max probe radius \(R_{max}\). What I did was using a skip threshold to discard the r distance that its pdf compare to pdf(center) is way to small, something like this: 

There is a magic value Kulla 13 suggest in the paper \(\sqrt{12.46v}\) that I was really not exactly sure how they come up the value. After asking original author Chris about this, he said that the value was calculated by integrate the Gaussian profile from 0 to r to get a value really close to 1, something like 0.9999, and the r value will be the \(R_{max}\). I did a quick experiment with this: \(\int_{0}^{2\pi}\int_{0}^{R_{max}}\frac{1}{2\pi v}e^{-r^2/2v}drd\theta = \int_{0}^{R_{max}}\frac{1}{v}e^{-r^2/2v}dr = 1 - e^{-R_{max}^2/2v} = 0.998\) then we can get result: \(R_{max} = \sqrt{12.429216v}\) That's pretty close to the magic value Chris suggest :) The other fun fact is that the alternative method I used earlier actually come up the same math equation to solve at last. While my threshhold value is a more aggressive 0.01, the idea is pretty similar: discard the sample with too small contribution.

With \(R_{max}\) in hand, now we gonna pick the axis to project sample disc. Kulla 13 suggested the ratio N:U:V = 2:1:1. I pretty much hardcode this dice throwing probability to the code. Few things to notice: 

1. besides the \(pdf_{disc}\), we now need to take into account the probability each axis got pick up(since they are 2:1:1, \(pdf_{N} = 0.5\), \(pdf_{U} = 0.25\), \(pdf_{V} = 0.25\) ). 

2. the probe ray length should be based on its distance to disc center, the longer the distance, the shorter the length.

3. dpdu and dpdv are not necessarily normalized values for their other purpose like anti-alias, we probably need to normalize it for initializing ray's direction (at least in my case)

The code snippet:

And here's a screenshot of testing with debug line visualizing the probe ray on one particular pixel(the red dot), probe rays with normal axis direction are blue, and dpdu for red, dpdv for green, the yellow dots are the probe ray intersection that we gonna do the irradiance calculation. YUP! That pixel is quite expensive!!

 

OK with all those probe rays we had, how we assign the appropriate weight sum for them with multiple importance sampling? According to Veach 97, the multi-importance estimator looks like this:

$$F=\sum_{i=1}^{n}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\omega_{i}(X_{i,j})\frac{f(X_{i,j})}{pdf_{i}(X_{i,j})}$$

where \(n\) is the number of sampling techniques (in this case: \(n=3\) we sampling in three ways, from normal, from dpdu and from dpdv), and \(n_{i}\) is the number of sample for technique \(i\). (in this case, \(n_{0} = 2c\), \(n_{1} = c\), \(n_{2} = c\), \(c\) is a constant since we casting probe ray with ratio N:U:V=2:1:1) I used the power heuristic approach Veach propose:

$$\omega_{i} = (n_{i}pdf_{i}(X_{i,j}))^2/\sum_{k=1}^{n}(n_{k}pdf_{k}(X_{k,j}))^2$$

The constant \(c\) got divided out in above equation so it doesn't matter we assign \(n_{i}\) with 2,1,1 or 4,2,2 or 2c, 1c, 1c. As for \(pdf_{k}(X_{k,j})\), it stands for the pdf evaluation on \(X_{i,j}\) with sampling technique k where k ranges from 1 to \(n_{i}\). (In English: "I got a sample point by picking N axis and Gaussian disc sampling around it, intersect the surface at point p, what's the \(pdf_{A}\) of this sample if I pick U or V axis and Gaussian disc sample around that axis projecting down to point p?") Friendly reminder: this pdf include the pdf picking the axis, the pdf of disc sampling and projecting it down to surface. (I actually forget putting pdf picking the axis into account in the first place and trouble shoot it for a few moments > <)

Put this into code:

This post summarizes the TODO item 4, we have all the fresh ingredients in hand now. The next post will be about putting them together to get the BSSRDF integration.