BSSRDF Importance Sampling 5 - Integrating BSSRDF and single scattering evaluation

The BSSRDF model Jensen 01 propose is "a sum of the diffusion approximation and the single scattering term"
$$S(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}}) = S_{d}(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}}) + S^{(1)}(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}})$$
The previous few posts all surrounds on the topic about: "how to get $S_{d}(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}})$ sample?", "How to get the sample distributed well for importance sampling?" Now we can put these parts all together :) Last look (for now) on the integration equation:

$$L( x_{o},\vec{\omega _{o}} )=\int_{A }\int_{\Omega }L_{i}( x_{i},\vec{\omega _{i}} )S_{d} (  x_{i},\vec{\omega _{i}}, x_{o},\vec{\omega _{o}} )cos\left ( \theta \right )d\omega dA$$
The basic flowchart of this part will be:
1. when intersect a BSSRDF material, start figure out the probe radius $R_{max}$ (this can be precomputed actually but considering $\sigma_{a}$ and $\sigma_{tr}$ may varying in texture space, I made the computation on the fly for now)
2. casting out probe rays for scene intersection test (we can actually just do the intersection test with this particular surface with BSSRDF material instead of intersection test the whole scene. This is a BIG TODO in my mind now for efficiency concern)
3. for each probe ray intersection sample point, calculate the irradiance (the current method is picking light based on power distribution, and only casting one lighting sample for each probe ray intersection, we can definitely do splitting here!(multiple lighting sample with one BSSRDF sample) another BIG TODO item)
4. evaluate $S_{d}(x_{i}, \vec{\omega_{i}}, x_{o}, \vec{\omega_{o}})=\frac{1}{\pi}F_{r}(\eta, \vec{\omega_{i}})R_{d}(\parallel x_{i}, x_{o} \parallel)F_{r}(\eta, \vec{\omega_{o}})$, the two Fresnel terms can be calculated with ray direction, intersection normal and the refraction index
5. add up the 3. 4. calculation with appropriate MIS weight we mentioned in the previous post, and divide it with corresponding pdf (pick the light, pick the geometry point on light surface, pick the axis, pick the disc sample, the ray/intersection normal dot product, all need to be taken into account)

and here's the naive implementation with the above flowchart:

	Color Renderer::LbssrdfDiffusion(const ScenePtr& scene,
	const Fragment& fragment, const BSSRDF* bssrdf, const Vector3& wo,
	const Sample& sample,
	const BSSRDFSampleIndex* bssrdfSampleIndex,
	WorldDebugData* debugData) const {
	const Vector3& pwo = fragment.getPosition();
	float coso = absdot(wo, fragment.getNormal());
	float eta = bssrdf->getEta();
	float Ft = 1.0f - Material::fresnelDieletric(coso, 1.0f, eta);
	float sigmaTr = bssrdf->getSigmaTr(fragment).luminance();
	const vector<Light*>& lights = scene->getLights();
	// figure out the sample probe radius, we ignore the integration
	// of area with pdf too small compare to center, yes...this introduces
	// bias...but can limit the probe ray as short as possible
	// skipRatio = pdf(r) / pdf(0, 0) = exp(-sigmaTr * r^2)
	float skipRatio = 0.01f;
	float Rmax = sqrt(log(skipRatio) / -sigmaTr);
	Color Lmultiscatter(0.0f);
	for(uint32_t i = 0; i < bssrdfSampleIndex->samplesNum; ++i) {
	const BSSRDFSample bssrdfSample(sample, *bssrdfSampleIndex, i);
	// sample a probe ray with gaussian falloff pdf from intersection
	Ray probeRay;
	float discPdf;
	BSSRDFSampleAxis axis = bssrdf->sampleProbeRay(fragment,
	bssrdfSample, sigmaTr, Rmax, &probeRay, &discPdf);
	Intersection probeIntersect;
	float epsilon;
	// TODO we can use a material indexed lookup table to query
	// surface with this particular BSSRDF instead of doing a whole
	// scene intersaction test
	if(scene->intersect(probeRay, &epsilon, &probeIntersect)) {
	if(probeIntersect.getMaterial()->getBSSRDF() == bssrdf) {
	const Fragment& probeFragment = probeIntersect.fragment;
	const Vector3& pProbe = probeFragment.getPosition();
	Color Rd = bssrdf->Rd(probeFragment,
	squaredLength(pProbe - pwo));
	// calculate the irradiance on the sample point
	float pickLightPdf;
	int lightIndex = mPowerDistribution->sampleDiscrete(
	bssrdfSample.uPickLight, &pickLightPdf);
	const Light* light = lights[lightIndex];
	Vector3 wi;
	float lightPdf;
	Ray shadowRay;
	const Vector3& ni = probeFragment.getNormal();
	Color L = light->sampleL(pProbe, epsilon, bssrdfSample.ls,
	&wi, &lightPdf, &shadowRay);
	if(L == Color::Black || lightPdf == 0.0f ||
	scene->intersect(shadowRay)) {
	continue;
	}
	float cosi = absdot(ni, wi);
	Color irradiance = L * cosi / (lightPdf * pickLightPdf);
	float Fti = 1.0f -
	Material::fresnelDieletric(cosi, 1.0f, eta);
	// evaluate the MIS weight
	float pdf = discPdf * absdot(probeRay.d, ni);
	float w = bssrdf->MISWeight(fragment, probeFragment, axis,
	pdf, sigmaTr, Rmax);
	Lmultiscatter +=
	(w * INV_PI * Ft * Fti * Rd * irradiance) / pdf;
	}
	}
	}
	Lmultiscatter /= (float)bssrdfSampleIndex->samplesNum;
	return Lmultiscatter;
	}

view raw diffusion_approximation.cpp hosted with ❤ by GitHub

Lots of stuff can be further optimized at this moment, but let's give it a quick taste on the new dishes :P, with just a little bit seasoning pepper: 25 samples

For comparison I put another simplest Lambert shade version. Same camera/lighting settings, same samples number:

we can see the dark edges(the head and tail part) in BRDF shading got soften a lot in BSSRDF since the lighting contribution from nearby area, also the bumpy body become not that hard lined in BSSRDF version.

Now moving to the single scattering part, as the term suggests, this part simulate the case that light specular refract into medium, encounter a scatter event to change its direction, refract out the medium and enter the camera. (or think in the other way: camera shoot eye rays with importance, refract into the medium, during the travel in the medium, it absorbs photon that emit from light and refract into the medium) The following images illustrate the single scattering idea:

'

Image taken from Hery 03, the whole radiance received by this eye ray would be the integration alone the refract ray $\vec{pT_{o}}$ in the medium

Put in Fresnel, scatter event, phase function into accout, Jensen 01 derive the equation for each sample point alone the refracted eye ray:

$$L^{(1)}(x_{o}, \vec{\omega_{o}}) = (\sigma_s(x_{o})Fp(\vec{{\omega_{i}}'},\vec{{\omega_{o}}'})/\sigma_{tc})e^{-{s_{i}}'\sigma_{t}(x_{i})}e^{-{s_{o}}'\sigma_{t}(x_{o})}L_{i}(x_{i}, \vec{\omega_{i}})$$

another round of scary symbol sort:
$F = F_{r}(x_{o}, \eta)F_{r}(x_{i}, \eta)$ : the Fresnel factor
$\sigma_{a}$: the absorb coefficient
$\sigma_{s}$: the scatter coefficient
$p$: Henry-Greenstein phase function, similar to BRDF for surface, this describe the scatter property for volume
$\vec{{\omega_{i}}'}$: photon $L_{i}$ refraction direction
$\vec{{\omega_{o}}'}$: eye ray refraction direction
$G = \parallel\vec{n_{x_{i}}}\cdot\vec{{\omega_{o}}'}\parallel/\parallel\vec{n_{x_{i}}}\cdot\vec{{\omega_{i}}'}\parallel$ : geometry factor
$\sigma_{t}=\sigma_{a} + \sigma_{s}$ : the extinction coefficient
$\sigma_{tc}=\sigma_{t}(x_{o}) + G\sigma_{t}(x_{i})$ : combined extinction coefficient
${s_{i}}'$: the refracted photon traveling distance
${s_{o}}'$: the refracted eye ray traveling distance

For the refract photon distance ${s_{i}}'$, it is actually pretty difficult to sample correctly since there is a specular refraction direction change from the inside medium sample to light sample point. Jensen proposed an approximation for the true refracted distance with Snell's law:
$${s_{i}}'=s_{i}\frac{\parallel\vec{n_{x_{i}}}\cdot\vec{\omega_{i}}\parallel}{\sqrt{1-(\frac{1}{\eta})^2(1-\parallel\vec{n_{x_{i}}}\cdot\vec{\omega_{i}}\parallel^2)}}$$
where $s_{i}$ is the distance between sample point in medium to light sample point minus the part this ray is not in the medium.

Jensen also suggested that we importance sampling alone the eye ray with $pdf(d) = \sigma_{t}e^{-\sigma_{t}d}$ since the exponential attenuation alone the medium. This 1D sample was covered in previous post , so we are pretty much ready to go!

The flowchar for single scattering will be:
1. when intersect a BSSRDF material, change the ray direction with specular refraction
2. alone the refraction ray sample point with 1D exponential falloff
3. for each sample point from 2. , sample light source and casting shadow ray to see if there is other geometry blocking the light besides the current medium (need to be a little bit careful on this part when playing around with the shadow ray length)
4. approximate ${s_{i}}'$ with Snell's law, evaluate Fresnel term and phase function

put them into codes:

	Color Renderer::LbssrdfSingle(const ScenePtr& scene,
	const Fragment& fragment, const BSSRDF* bssrdf, const Vector3& wo,
	const Sample& sample,
	const BSSRDFSampleIndex* bssrdfSampleIndex,
	WorldDebugData* debugData) const {
	const Vector3& pwo = fragment.getPosition();
	const Vector3& no = fragment.getNormal();
	float coso = absdot(wo, fragment.getNormal());
	float eta = bssrdf->getEta();
	float Ft = 1.0f - Material::fresnelDieletric(coso, 1.0f, eta);
	float sigmaTr = bssrdf->getSigmaTr(fragment).luminance();
	Color scatter = bssrdf->getScatter(fragment);
	Color sigmaT = bssrdf->getAttenuation(fragment);
	Vector3 woRefract = Goblin::specularRefract(wo, no, 1.0f, eta);
	const vector<Light*>& lights = scene->getLights();
	Color Lsinglescatter(0.0f);
	// single scattering part
	for(uint32_t i = 0; i < bssrdfSampleIndex->samplesNum; ++i) {
	const BSSRDFSample bssrdfSample(sample, *bssrdfSampleIndex, i);
	// sample a distance with exponential falloff
	float d = exponentialSample(bssrdfSample.uSingleScatter, sigmaTr);
	Vector3 pSample = pwo + d * woRefract;
	float samplePdf = exponentialPdf(d, sigmaTr);
	// sample a light from pSample
	float pickLightPdf;
	int lightIndex = mPowerDistribution->sampleDiscrete(
	bssrdfSample.uPickLight, &pickLightPdf);
	const Light* light = lights[lightIndex];
	Vector3 wi;
	float epsilon, lightPdf;
	Ray shadowRay;
	Color L = light->sampleL(pSample, 1e-5f, bssrdfSample.ls,
	&wi, &lightPdf, &shadowRay);
	if(L == Color::Black || lightPdf == 0.0f) {
	continue;
	}
	// preserve the maxt since the next intersection test will modify
	// it to the closest intersection, but we still need this maxt to
	// cast shadow ray from that intersection to light
	float maxt = shadowRay.maxt;
	Intersection wiIntersect;
	// TODO we can use a material indexed lookup table to query
	// surface with this particular BSSRDF instead of doing a whole
	// scene intersaction test
	if(scene->intersect(shadowRay, &epsilon, &wiIntersect)) {
	// if not, the pSample is out of the BSSRDF geometry already
	if(wiIntersect.getMaterial()->getBSSRDF() == bssrdf) {
	// update shadow ray to start from pwi
	const Fragment& fwi = wiIntersect.fragment;
	const Vector3& pwi = fwi.getPosition();
	const Vector3& ni = fwi.getNormal();
	shadowRay.mint = shadowRay.maxt + epsilon;
	shadowRay.maxt = maxt;
	if(!scene->intersect(shadowRay)) {
	float p = bssrdf->phase(wi, woRefract);
	float cosi = absdot(ni, wi);
	float Fti = 1.0f -
	Material::fresnelDieletric(cosi, 1.0f, eta);
	Color sigmaTi = bssrdf->getAttenuation(fwi);
	float G = absdot(ni, woRefract) / cosi;
	Color sigmaTC = sigmaT + G * sigmaTi;
	float di = length(pwi - pSample);
	float et = 1.0f / eta;
	float diPrime = di * absdot(wi, ni) /
	sqrt(1.0f - et * et * (1.0f - cosi * cosi));
	Lsinglescatter += (Ft * Fti * p * scatter / sigmaTC) *
	expColor(-diPrime * sigmaTi) *
	expColor(-d * sigmaT) * L /
	(lightPdf * pickLightPdf * samplePdf);
	}
	}
	}
	}
	Lsinglescatter /= (float)bssrdfSampleIndex->samplesNum;
	return Lsinglescatter;
	}

view raw bssrdf_single_scattering.cpp hosted with ❤ by GitHub

and we got the single scattering part :) (Though we still face the similar further optimization in TODO list like splitting and direct surface intersection test)

One more! in Jensen 01's final demo image marble bust he added a Fresnel reflection pass. Though it's not necessary, I feel that shiny sparkling highlight is pretty charming. Plus, we already have that in our BRDF, so I just directly plug it into the new Subsurface material:

All the jigsaw pieces are put in place now. Finally! Put them altogether, we have this shiny marble cutie :)

So how's the multiple importance sampling benefit we have? That's actually the most exciting part! For comparison, I render another image with sample alone normal axis only without multiple importance sampling, same 25 samples.

And the side by side comparison (left side no MIS, right side with MIS):

In the highly curvy surface, The MIS version did a pretty good job removing those most noticeable variance noise, while in the relatively flat surface(the dragon right side body part), project all samples alone normal has a better result.  Kulla 13 did mentioned about this part that there will be performance lost in overall flat surface since we would have extra wasted sample projected alone dpdu/dpdv and intersect nothing. The algorithm generally remove the spiky high variance but introduce some relatively more uniform variance, which is easier to converge when the sample amount got cranked up.

Image taken from Kulla 13, explains the issue that MIS doesn't work particular well for overall flat surface

Alright! TODO item 5, 6 done! I probably will post up some more interesting rendering images later, but these few posts summarize what my research/interpretation about Kulla 13 and Jensen 01. Math is definitely not my strongest field so if there's anything wrong in the post. Please let me know. Thanks!