BSSRDF Importance Sampling 2 - Dipole Diffusion Approximation101

Super quick review on the dipole diffusion approximation model.
Compare to the classic BRDF integration:
$$L\left ( x_{o},\vec{\omega _{o}} \right )=L_{e}\left ( x_{o},\vec{\omega _{o}} \right )+\int_{\Omega }L_{i}\left ( x_{i},\vec{\omega _{i}} \right )f_{r}\left ( \vec{\omega _{i}},\vec{\omega _{o}} \right )cos\left ( \theta \right )d\omega$$
This is BSSRDF integration:
$$L\left ( x_{o},\vec{\omega _{o}} \right )=\int_{A}\int_{\Omega  }L_{i}\left ( x_{i},\vec{\omega _{i}} \right )S\left (  x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}} \right )cos\left ( \theta \right )d\omega dA$$

In short words, to evaluate the radiance come out from surface point $x_{o}$ with direction $\omega_{o}$, we need to integrate over the whole surface area (the $\int_{A }dA$ part ), at each surface point $x_{i}$, we integrate over the hemisphere(the $\int_{\Omega }d\omega$ part) that for radiance $L( x_{i},\vec{\omega _{i}})$, how much portion arrive at point $x_{o}$ with direction $\omega_{o}$(the $S( x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}})$ part)

the BSSRDF $S(  x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}})$ is a 4 dimensional function, but most of the time it got simplified to radially symmetric and become a 1D distance funtion $R_{d}(d)$ where d is distance between $x_{o}$ and $x_{i}$.

Dipole Diffusion Approximation made assumptions that:
-light distribution  in highly scattering media tends to become isotropic. (My understanding for this part is that when light enter into medium with really high Albedo, after multiple scattering event happens, the light scatter distribution will become more and more uniform, and approach $1/4\pi$ at last)
-the surface is semi-infinite flat surface since for finite media there is no general way to solve diffusion equation in analytic way.

By putting two point light source above and under the medium(that's why it's called dipole method), we can get the analytic solution for fluence on the semi-infinite flat surface (OK...."we" not including "me" at this moment......I confess...... Donner 06 thesis Chapter 5 contains a detail derivation on this, while I just grab the result equation directly...)

$$R_{d}\left (r\right )=\frac{{\alpha}'}{4\pi}\biggl [z_{r}( \sigma_{tr}d_{r}+1)\frac{e^{-\sigma_{tr}d_{r}}}{ {d_{r}}^{3} } + z_{v}( \sigma_{tr}d_{v}+1)\frac{e^{-\sigma_{tr}d_{v}}}{ {d_{v}}^{3} }\biggr ]$$

After a topological sort on those scary symbols:
$\sigma_{a}$: absorb coefficient, describe how much radiance got absorbed in medium
$\sigma_{s}$: scatter coefficient, describe how much radiance got scattered in medium
$g$: anisotropy parameter (range from -1 to 1, -1 for fully backward scattering, 1 for fully forward scattering, 0 for isotropy)
${\sigma_{s}}' =\sigma_{s}(1-g)$: reduced scattering coefficient
${\sigma_{t}}' ={\sigma_{s}}' + \sigma_{a}$: reduced extinction coefficient
${\alpha}' = {\sigma_{s}}'/{\sigma_{t}}'$: the reduced albedo
$\sigma_{tr} = \sqrt{3\sigma_{a}{\sigma_{t}}'}$: the effective transport coefficient
$F_{r}$: the Fresnel formula for the reflection at a dielectric surface
$F_{dr} = \int_{\Omega }F_{r}(\eta, \vec{n}\cdot \vec{\omega})(\vec{n}\cdot \vec{\omega})d\omega$: the average diffuse Fresnel reflectance
$A=(1+F_{dr})/(1-F_{dr})$: the coefficient takes into account the effect of internal diffuse reflection
$D=\frac{1}{3{\sigma_{t}}'}$: the diffusion constant
$z_{r} = 1/{\sigma_{t}}'$: the distance beneath surface where we put the positive dipole light
$z_{v} = z_{r} + 4AD$: the distance above surface where we put the negative dipole light
$d_{r} =\parallel x-x_{r}\parallel$: the distance from x to the positive dipole light
$d_{v} =\parallel x-x_{v}\parallel$: the distance from x to the negative dipole light

Hopefully this dipole image illustrates the thing a little bit clearer...

At last, we put into account the Fresnel reflection at the boundary for both the incoming light and the outgoing light, we get:

$$S_{d}\left (x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}}\right)=\frac{1}{\pi}F_{r}(\eta, \vec{\omega_{i}})R_{d}(\parallel x_{i}, x_{o} \parallel)F_{r}(\eta, \vec{\omega_{o}})$$

Image taken from GPU Gems3, showing the BSSRDF falloff pretty quick when the distance between $x_{o}$ and $x_{i}$ increases, we can also see that red channel tend to resist more before fully extinct, this sorta explains the idea human skin have that red texture underneath

For Fresnel diffuse reflection $F_{dr}$ , it's an integration that can be approximated with the following code snippet. It was taken from Donner 06 too:

	inline float BSSRDF::Fdr(float eta) {
	// see Donner. C 2006 Chapter 5
	// the internal Fresnel reflectivity
	// approximated with a simple polynomial expansion
	if(eta < 1.0f) {
	return -0.4399f + 0.7099f / eta - 0.3319f / (eta * eta) +
	0.0636f / (eta * eta * eta);
	} else {
	return -1.4399f / (eta * eta) + 0.7099f / eta + 0.6681f +
	0.0636f * eta;
	}
	}

view raw fdr_approximation.hpp hosted with ❤ by GitHub

YUUUP.....we got the diffusion approximation part BSSRDF $S_{d}(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}})$ now, and it looks pretty mathmatically daunting (again, for me....). So, I encapsulate all these stuff into the BSSRDF black box, and hopefully I don't need to see them most of the time :P The constructor takes $\sigma_{a}$, ${\sigma_{s}}'$, refraction index $\eta$, anisotropy parameter $g$ , the rest of the stuff it can cook indoor. Something looks like this:

	BSSRDF::BSSRDF(const ColorTexturePtr& absorb,
	const ColorTexturePtr& scatterPrime, float eta, float g):
	mAbsorb(absorb), mScatterPrime(scatterPrime),
	mEta(eta), mG(g) {
	float fdr = Fdr(eta);
	mA = (1.0f + fdr) / (1.0f - fdr);
	}
	Color BSSRDF::Rd(const Fragment& fragment, float d2) const {
	// see Donner. C 2006 Chapter 5 for the full derivation
	// of the following disffusion dipole approximation equation
	Color sigmaA = mAbsorb->lookup(fragment);
	Color sigmaSPrime = mScatterPrime->lookup(fragment);
	Color sigmaTPrime = sigmaA + sigmaSPrime;
	Color sigmaTr = sqrtColor(3.0f * sigmaA * sigmaTPrime);
	Color one(1.0f);
	Color zr = one / sigmaTPrime;
	// zv = zr + 4AD where D = 1/(3 * sigmaT') = zr / 3
	Color zv = zr * (1.0f + 4.0f / 3.0f * mA);
	Color dr = sqrtColor(zr * zr + Color(d2));
	Color dv = sqrtColor(zv * zv + Color(d2));
	Color alphaPrime = sigmaSPrime / sigmaTPrime;
	Color sTrDr = sigmaTr * dr;
	Color sTrDv = sigmaTr * dv;
	Color rd = 0.25f * INV_PI * alphaPrime * (
	(zr * (one + sTrDr) * expColor(-sTrDr) / (dr * dr * dr)) +
	(zv * (one + sTrDv) * expColor(-sTrDv) / (dv * dv * dv)));
	return clampColor(rd);
	}

view raw BSSRDF_Rd.cpp hosted with ❤ by GitHub

OK! That pretty much sums up the TODO list item1, now moves on to item2, get a separate integration for BSSRDF. This is just open another cut in point for some future implementation, pretty straightforward for small side project scale code base:

	Color Renderer::Lsubsurface(const ScenePtr& scene,
	const Intersection& intersection, const Vector3& wo,
	const Sample& sample,
	const BSSRDFSampleIndex* bssrdfSampleIndex) const {
	const MaterialPtr& material =
	intersection.primitive->getMaterial();
	const BSSRDF* bssrdf = material->getBSSRDF();
	const vector<Light*>& lights = scene->getLights();
	if(bssrdf == NULL || lights.size() == 0) {
	return Color(0.0f);
	}
	const Fragment& fragment = intersection.fragment;
	Color Lsinglescatter = LbssrdfSingle(scene, fragment, bssrdf,
	wo, sample, bssrdfSampleIndex);
	// multiple scattering part with diffusion approximation
	Color Lmultiscatter = LbssrdfDiffusion(scene, fragment, bssrdf,
	wo, sample, bssrdfSampleIndex);
	return Lsinglescatter + Lmultiscatter;
	}

view raw sss_entry.cpp hosted with ❤ by GitHub

Done!

One particular thought pop in my mind while reading Jensen01 was: "For the sake of god!!! To just get something radially symmetrically falloff you really need to be that hardcore? I mean....there should be something so much easier symmetrically falloff model you can use right? like Gaussian distribution maybe?" The fun fact is yes! There are already some research about using multiple weight sum of Gaussian filters to approximate BSSRDF and it can be really close to the dipole model. This GPU gems 3 article mentioned the idea and give you the weight sum parameters directly for skin BSSRDF. If you are interested in the full story under the hood, d'Eon11 "A Quantized-Diffusion Model for Rendering Translucent Materials" is a cool ride on it.

Image taken from GPU Gems3 , we can see the 4 Gaussian weight sum is really close to the Dipole model already