BSSRDF Importance Sampling 3 - exponential falloff sampling on 1D and 2D space

Take one more look on our 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$$

The solid angle part seems to be easier since we already got methods to sampling lights at specified intersection point for our previous work on BRDF integration. The area integral part is more interesting: how we integrate the BSSRDF surface. Sure the naive way will be uniform sample the BSSRDF surface area, but look at that dipole model we got, that radically symmetric fast falloff: how much variance we gonna suffer before converge the clean render if we do things by uniform sampling? (and how much ray casting we gonna do to figure out the projection point on surface?) There are definitely some better ways......we should definitely try importance sampling.

OK importance sampling, then what is the dominant factor in that Rd equation? In Jensen01, he mentioned "we use standard Monte Carlo techniques to randomly sample the surface with density ($\sigma_{tr}e^{-\sigma_{tr}d}$) at some distance $d$ from $x_{o}$" which means he suggests the dominant factor to be $\sigma_{tr}e^{-\sigma_{tr}d}$, though I feel $\sigma_{tr}e^{-\sigma_{tr}d}/d^2$ is the closest one. Jensen said the word, let's try it~ time to brush up the inverse method sampling


-------------------------------WARNING DEAD END ENTRY!-------------------------------
I actually fail to get this answer derived out and change the route at last, but still think I should write down what I thought at that time...if you know how to get this derivation out right, please definitely let me know! I'm really curious about it

pdf in 2D space: $pdf(x, y) = ce^{-\sigma_{tr}\sqrt{x^2+y^2}}$ where c is the constant we try to figure out to normalize the pdf
transform it into polar coordinate, the Jacobian is:
$$J_{r} = \begin{pmatrix} \frac{\mathrm{d} x}{\mathrm{d} r} & \frac{\mathrm{d} x}{\mathrm{d} \theta}\\ \frac{\mathrm{d} y}{\mathrm{d} r} & \frac{\mathrm{d} y}{\mathrm{d} \theta} \end{pmatrix} = r$$
so $pdf(r, \theta) = rpdf(x, y) = cre^{-\sigma_{tr}r}$
$cdf(r, \theta) = \int_{0}^{2\pi }\int_{0}^{\infty }pdf(r, \theta)drd\theta = 2\pi c\int_{0}^{\infty }re^{-\sigma_{tr}r}dr = 1$

now integration by part:

$\int u{v}'dr = uv -\int {u}'vdr$
$u=r$ and $\frac{\mathrm{d} v}{\mathrm{d} r}=e^{-\sigma_{tr}r}$ =>
$v = e^{-\sigma_{tr}r}/-\sigma_{tr}$ and $\frac{\mathrm{d} v}{\mathrm{d} r}=1$
$\int_{0}^{\infty }re^{-\sigma_{tr}r}dr = \frac{re^{-\sigma_{tr}r}}{-\sigma_{tr}}\bigg\vert_{r=0}^{r=\infty}-\int_{0}^{\infty }\frac{e^{-\sigma_{tr}r}}{-\sigma_{tr}}dr = \frac{-1}{\sigma_{tr}^2}(\sigma_{tr} + 1)e^{-\sigma_{tr}r}\bigg\vert_{r=0}^{r=\infty}= \frac{1}{\sigma_{tr}^2} =\frac{1}{2c\pi}$
$c=\frac{\sigma_{tr}^2}{2\pi}$
$pdf(r,\theta) = \frac{\sigma_{tr}^2}{2\pi}re^{-\sigma_{tr}r}$

almost there, now integrate out the $\theta$ part to get the marginal $pdf(r) = \sigma_{tr}^2re^{-\sigma_{tr}r}$
then $cdf(r) = \int_{0}^{r}\sigma_{tr}^2 re^{-\sigma_{tr}r}dr = 1 - (\sigma_{tr}r + 1)e^{-\sigma_{tr}r}$
and now we just need to get the inverse function of r......and then I stuck........I can't seem to find a easy way to inverse this function......

By doing some google work, it seems that I bump into something called Laplace Transform Inverse , and this stuff kinda out of my college freshman level calculus....Hery 03 seems use a precompute table to solve this sampling part, he also kindly offered the source code for precomuting process, but I'm not sure if I really get what that code does......maybe go learn some advanced calculus can shed some lights for me on this issue, until then, I need to find some workaround.....

-----------------------------------EXIT OF THE DEAD END-----------------------------------

How about Gaussian? Like I mentioned in previous post, we can approximate dipole model with weight sum of Gaussian filters, plus Gaussian's falloff speed is higher than the above exponential falloff, consider there's a $(1+\sigma_{tr}d)/d^3$ factor the above pdf didn't take into account. Maybe Gaussian turned out to be a better fit?

Out of curiosity, I throw in the above 3 normalized equation into some online math drawing tool, the result shows up in above image: green for $e^{-\sigma_{tr}d}$ (the Jensen suggestion), red for $e^{-\sigma_{tr}d^2}$ (Gaussian, which is my guess), blue for $(1+\sigma_{tr}d)e^{-\sigma_{tr}}/d^3$ (which is probably the closest factor to dipole model, but looks really hard to figure out the inverse function) , and ya, Gaussian actually fit closer to the dipole model in this case. Worst of all, it falls off with d increase so it's always better than uniform sample. importance sampling works as long as the pdf is closer to original function, but doesn't need to perfectly fit, right?

So we try to do things Gaussian way now, let's try the inverse method again:

pdf in 2D space: $pdf(x, y) = ce^{-\sigma_{tr}(x^2+y^2)}$ where c is the constant we try to figure out to normalize the pdf
so $pdf(r, \theta) = rpdf(x, y) = cre^{-\sigma_{tr}r^2}$
$cdf(r, \theta) = \int_{0}^{2\pi }\int_{0}^{\infty }pdf(r, \theta)drd\theta = 2\pi c\int_{0}^{\infty }re^{-\sigma_{tr}r^2}dr = 2\pi c(\frac{e^{-\sigma_{tr}r^2}}{-2\sigma_{tr}}\bigg\vert_{r=0}^{r=\infty}) = \frac{c\pi}{\sigma_{tr}}= 1$
$c=\frac{\sigma_{tr}}{\pi}$
$pdf(r, \theta) = \frac{\sigma_{tr}}{\pi}re^{-\sigma_{tr}r^2}$

integrate out the theta factor:
$pdr(r) = 2\sigma_{tr}re^{-\sigma_{tr}r^2}$
$cdf(r) = \int_{0}^{r}pdf(r) = -e^{-\sigma_{tr}r^2}\bigg\vert_{0}^{r} = 1 - e^{-\sigma_{tr}r^2}$
inverse function:
$\xi = 1 - e^{-\sigma_{tr}r^2}$
$r =\sqrt{\frac{\ln{(1- \xi)}}{-\sigma_{tr}}}$
and we can further simplified it to:
$r =\sqrt{\frac{\ln{(\xi)}}{-\sigma_{tr}}}$
since $\xi$ is uniform distribution between 0 and 1
and $\theta = 2\pi\xi$

whew~~that's much easier than the previous one....

The above case works when we assume the sample distribute in the whole 2D space( $-\infty ->\infty$, how if we want to sample on a disc with radius $R_{max}$? You can actually just replace the $\infty$ integration range to $R_{max}$ to derive it out the result. I'll just list the result here directly since it's almost the exact same process:
$r =\sqrt{\frac{\ln{(1- \xi(1-e^{-\sigma_{tr}R_{max}^2}))}}{-\sigma_{tr}}}$
$\theta=2\pi\xi$
$pdf(x,y) = ( \frac{\sigma_{tr}}{\pi}e^{-\sigma_{tr}r^{(x^2+y^2)}})/(1-e^{-\sigma_{tr}R_{max}})$
The reason I brought up this explicit radius version is that we gonna use it later and the original BSSRDF Importance Sampling paper use this as their sampling method. This post basically explain how that equation in paper come out. :)  (it may looks tiny different but if you replace the $\frac{1}{2v}$ in paper with $\sigma_{tr}$ , they are exactly the same thing)

Last bonus one, sample in 1D with distribution proportional to $e^{-ax}$ , this one is used in BSSRDF single scattering and actually the easiest one to derive.

$pdf(x) = ce^{-ax}$
$cdf(x) = \int_{0}^{\infty }ce^{-ax}dx = \frac{c}{-a}e^{-ax}\bigg\vert_{0}^{\infty}= \frac{c}{a} = 1$
$c=a$
$pdf(x) = ae^{-ax}$
$cdf(x) =  \int_{0}^{x}ae^{-ax}dx = 1 - e^{-ax}$
$\xi = 1 - e^{-ax}$
$x = \frac{ln{(1-\xi) }}{-a}$
simplified to
$x = \frac{ln{\xi }}{-a}$

Boiling those stuff down to codes are just the following few lines:

	/*
	* let falloff be a
	* since the pdf is proportional to exp(-a * (x * x + y * y))
	* let pdf(x, y) = c * exp(-a * (x * x + y * y))
	* convert to spherical coordinate:
	* pdf(r, theta) = c * r * exp(-a * r * r)
	* integrate c * r * exp(-a * r * r) over r (from 0 to inf) and theta = 1 ->
	* c = a / pi -> pdf(r, theta) = a / pi * r * exp(-a * r * r) ->
	* pdf(x, y) = a / pi * exp(-a * (x * x + y * y))
	* marginal pdf(r) = integrage pdf((r, theta) over 0 to 2pi
	* = 2 * a * r * exp(-a * r * r)
	* conditional pdf(theta|r) = 1 / 2pi
	* cdf(r) = integrate 2 * a * r * exp(-a * r * r) over 0 to r ->
	* cdf(r) = 1 - exp(-a * r * r)
	* cdf(theta) = integrate 1 / 2pi over 0 to theta = theta / 2pi
	* inverse method:
	* u1 = 1 - exp(-a * r * r) -> r = sqrt(ln(1 - u1) / -a) ->
	* r = sqrt(ln(u1) / -a) since u1 is 0-1 uniform distribution
	* u2 = theta / 2pi -> theta = 2pi * u2
	* transform it back from spherical coordinate:
	* x = r * cos(theta)
	* y = r * sin(theta)
	*/
	Vector2 gaussianSample2D(float u1, float u2, float falloff) {
	float r = sqrtf(log(u1) / -falloff);
	float theta = TWO_PI * u2;
	return Vector2(r * cos(theta), r * sin(theta));
	}
	/*
	* similar to the above gaussianSample2D, but the integrate range of r
	* from 0 to Rmax, which makes us able to sample a disc with radius Rmax in
	* gaussian distribution
	*/
	Vector2 gaussianSample2D(float u1, float u2, float falloff, float Rmax) {
	float r = sqrtf(log(1.0f - u1 * (1.0f - exp(-falloff * Rmax * Rmax))) /
	-falloff);
	float theta = TWO_PI * u2;
	return Vector2(r * cos(theta), r * sin(theta));
	}
	/*
	* integrate c * exp(-falloff * x) from 0 to inf = 1 ->
	* c = falloff -> pdf(x) = falloff * exp(-falloff * x)
	* cdf(x) = 1 - exp(-falloff * x) -> u = 1 - exp(-falloff * x) ->
	* x = -ln(1 - u) / falloff -> simplified to x = -ln(u) / falloff
	* since u are [0, 1) uniform distribution
	*/
	float exponentialSample(float u, float falloff) {
	return -log(u) / falloff;
	}
	float gaussianSample2DPdf(const Vector3& pCenter,
	const Vector3& pSample, const Vector3& N, float falloff) {
	Vector3 d = pSample - pCenter;
	Vector3 projected = d - N * dot(d, N);
	return INV_PI * falloff * exp(-falloff * squaredLength(projected));
	}
	float gaussianSample2DPdf(const Vector3& pCenter,
	const Vector3& pSample, const Vector3& N, float falloff, float Rmax) {
	return gaussianSample2DPdf(pCenter, pSample, N, falloff) /
	(1.0f - exp(-falloff * Rmax * Rmax));
	}
	float gaussianSample2DPdf(float x, float y, float falloff) {
	return INV_PI * falloff * exp(-falloff * (x * x + y * y));
	}
	float gaussianSample2DPdf(float x, float y, float falloff,
	float Rmax) {
	return gaussianSample2DPdf(x, y, falloff) /
	(1.0f - exp(-falloff * Rmax * Rmax));
	}
	// project pSample on the plane form by pCenter and N, then calculate
	// the corresponding gaussian 2D pdf
	float gaussianSample2DPdf(const Vector3& pCenter,
	const Vector3& pSample, const Vector3& N, float falloff);
	float gaussianSample2DPdf(const Vector3& pCenter,
	const Vector3& pSample, const Vector3& N, float falloff, float Rmax);
	float exponentialPdf(float x, float falloff) {
	return falloff * exp(-falloff * x);
	}

view raw exponential_sample.cpp hosted with ❤ by GitHub

And here's some testing result with our freshly cooked sampling strategy, for comparison I put in another set of uniform disc sampling implemented based on Shirly 97

left up: 100 Gaussian sample with max radius 2 and falloff 2

right up: 1024 Gaussian sample with max radius 2 and falloff 2

left down: 100 uniform disc sample with radius 2

right down: 1024 uniform disc sample with radius 2

Now we have the TODO list item 3 in hand. We can start to look into the multiple importance sampling method BSSRDF Importance Sampling proposed, which in my opinion, is the most exciting part for this round.