BSSRDF Importance Sampling 6 - params tweaking, test render and future wishlist

This post actually has nothing to do with importance sampling at all already :P To get some relatively more interesting render images, I tried to tweak the reflection albedo for some different colors, but how does that get mapped to corresponding $\sigma_{a}$ and ${\sigma_{s}}'$ ?(what BSSRDF takes as constructor input)

In Jensen 01 there is a BRDF approximation of BSSRDF based on the assumption that the incident illumination is uniform, and they can integrate the dipole diffusion approximation to get a total diffuse reflectance $R_{d}$ at a surface intersection point:

$$R_{d} = 2\pi\int_{0}^{\infty}rR_{d}(r)dr = \frac{{\alpha}'}{2}(1 + e^{-\frac{4}{3}A\sqrt{3(1-{\alpha}')}}) e^{-\sqrt{3(1-{\alpha}')}}$$
Jensen 02 further purposed that using this approximated $R_{d}$(the diffusion reflection color) and diffuse mean free path $l_{d}$ (the average distance light travels in medium before scattering, my understanding is how translucent the medium is) as input parameters to derived out $\sigma_{a}$ and ${\sigma_{s}}'$: 

1. first get reduced albedo ${\alpha}'$ from above diffuse reflectance equation (the function is hard to invert but it's monotonic so Jensen suggests just using numerical approximation)

2. get effective transport coefficient: $\sigma_{tr} \approx 1/l_{d}$

3. get reduced extinction coefficient: ${\sigma_t}'=\frac{\sigma_{tr}}{\sqrt{3(1-{\alpha}')}}$

4. get reduced scatter coefficient: ${\sigma_{s}}'={\alpha}'{\sigma_{t}}'$

5, get absorb coefficient: $\sigma_{a} = {\sigma_{t}}' - {\sigma_{s}}'$

and we got the ingredient to feed in to BSSRDF. The code looks like this (the numerical approximation part I just modify it from pbrt, it actually looks kinda like some google style interview question :P)

	void BSSRDF::convertFromDiffuse(const Color& Kd,
	const Color& diffuseMeanFreePath, float A,
	Color* absorb, Color* scatterPrime) {
	float diffuse[3] = {Kd.r, Kd.g, Kd.b};
	float sigmaTr[3] = {
	1.0f / diffuseMeanFreePath.r,
	1.0f / diffuseMeanFreePath.g,
	1.0f / diffuseMeanFreePath.b};
	float sigmaSPrime[3];
	float sigmaA[3];
	// RGB channel
	for(int i = 0; i < 3; ++i) {
	// use few binary search iterations to approximate alpha prime
	// since there is no easy way to inverse diffuseReflectance method
	// directly...
	float alphaLow = 0.0f;
	float alphaHigh = 1.0f;
	float rdLow = diffuseReflectance(alphaLow, A);
	float rdHigh = diffuseReflectance(alphaHigh, A);
	for(int j = 0; j < 16; ++j) {
	float alphaMid = 0.5f * (alphaLow + alphaHigh);
	float rdMid = diffuseReflectance(alphaMid, A);
	if(rdMid > diffuse[i]) {
	alphaHigh = alphaMid;
	rdHigh = rdMid;
	} else {
	alphaLow = alphaMid;
	rdLow = rdMid;
	}
	}
	float alphaPrime = 0.5f * (alphaLow + alphaHigh);
	float sigmaTPrime = sigmaTr[i] /
	sqrt(3.0f * (1.0f - alphaPrime));
	sigmaSPrime[i] = alphaPrime * sigmaTPrime;
	sigmaA[i] = sigmaTPrime - sigmaSPrime[i];
	}
	*scatterPrime = Color(sigmaSPrime[0], sigmaSPrime[1], sigmaSPrime[2]);
	*absorb = Color(sigmaA[0], sigmaA[1], sigmaA[2]);
	}
	float BSSRDF::diffuseReflectance(float alphaPrime, float A) {
	float sqrtTerm = sqrt(3.0f * (1.0f - alphaPrime));
	return 0.5f * alphaPrime *
	(1.0f + exp(-(4.0f / 3.0f) * A * sqrtTerm)) *
	exp(-sqrtTerm);
	}

view raw bssrdf_from_diffuse.cpp hosted with ❤ by GitHub

with the above tweaking method in hand, we got some test rendering images, all use a same random $R_{d}$ color (0.478431, 0.513725, 0.521569) from color picker and different mean free path parameters (1-8 from up to down, you can see it gets more and more translucent, while I feel the last few is kinda way too over the top already :P ):

and a diffuse Lambert render reference:

And it's postmortem time, what's not there and should get addressed later:
1. more dynamic splitting: at this moment each BSSRDF probe ray sample only pair with one light sample. For simple lighting setup maybe works fine, but it definitely going to be a pain when the lighting getting more complex. Should make it available to specify how many lighting sample can be pair with each probe ray sample. Also the single scattering evaluation share the same lighting sample at this moment, they can probably get split up too or let user have the option to skip single scattering(which is also quite expensive to evaluate but relatively subtle effect for BSSRDF material like skin)

2. look up surface mesh in scene hierarchy with material id: we are still shooting probe ray to the whole scene graph for intersection test, which is not necessarily and should only do the intersection test against the surface mesh with same material. If the probe ray only intersect with a limited portion of surface, it would definitely speed up the integration process.

3. indirect lighting. Yes I confess...we only calculate the direct lighting for each light sample we cast out for BSSRDF so all our GI information is only appeared on those specular reflection area. The naive implementation on this shouldn't be hard (just like regular path tracing shoot out bounce path for each probe intersection to gather bounce lighting information) but I'm kinda afraid that it's gonna be really expensive. Thinking about solution at this moment......

4. use different max probe radius $R_{max}$ for different spectrum channel. There is only one probe radius for all spectrum channel now, but the diffusion profile shows clearly that red channel last much longer than the other two, it doesn't make sense to use the same probe radius for all of them. Kulla 13 mentioned about they also random picking the probe radius with different diffusion profile and combine the result with MIS for their skin rendering, while I haven't moved forward to this part yet.

5. the ${\sigma_{s}}'$ and $\sigma_{a}$ coefficient are measured in per mm unit, which I believe the scale of geometry will cause look difference at this moment, it would be good if we can come up some normalization mechanics for these coefficient tweaking so the artistic control can be more intuitive.

There should be more, and I'll append them when they pop up in my mind :)

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!

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: 

	// skipRatio = pdf(r) / pdf(0, 0) = exp(-sigmaTr * r^2)
	float skipRatio = 0.01f;
	float Rmax = sqrt(log(skipRatio) / -sigmaTr);

view raw get_rmax.cpp hosted with ❤ by GitHub

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:

	enum BSSRDFSampleAxis {
	UAxis,
	VAxis,
	NAxis
	};
	BSSRDFSampleAxis BSSRDF::sampleProbeRay(const Fragment& fragment,
	const BSSRDFSample& sample, float sigmaTr, float Rmax,
	Ray* probeRay, float* pdf) const {
	Matrix3 shadeToWorld = fragment.getWorldToShade().transpose();
	const Vector3& pwo = fragment.getPosition();
	// sample disk with gaussian falloff
	Vector2 pSample = gaussianSample2D(
	sample.uDisc[0], sample.uDisc[1], sigmaTr, Rmax);
	float halfProbeLength = sqrt(Rmax * Rmax - pSample.squaredLength());
	// figure out we should sample alone U or V or N
	// The chance to pick up U, V, N is 1 : 1 : 2
	BSSRDFSampleAxis axis;
	if(sample.uPickAxis <= 0.5f) {
	probeRay->o = pwo + shadeToWorld *
	Vector3(pSample.x, pSample.y, -halfProbeLength);
	probeRay->d = fragment.getNormal();
	axis = NAxis;
	*pdf = 0.5f;
	} else if(sample.uPickAxis <= 0.75f) {
	probeRay->o = pwo + shadeToWorld *
	Vector3(-halfProbeLength, pSample.x, pSample.y);
	probeRay->d = normalize(fragment.getDPDU());
	axis = UAxis;
	*pdf = 0.25f;
	} else {
	probeRay->o = pwo + shadeToWorld *
	Vector3(pSample.y, -halfProbeLength, pSample.x);
	probeRay->d = normalize(fragment.getDPDV());
	axis = VAxis;
	*pdf = 0.25f;
	}
	probeRay->mint = 0.0f;
	probeRay->maxt = 2.0f * halfProbeLength;
	*pdf *= gaussianSample2DPdf(pSample.x, pSample.y, sigmaTr, Rmax);
	return axis;
	}

view raw sample_probe_ray.cpp hosted with ❤ by GitHub

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:

	float BSSRDF::MISWeight(const Fragment& fo, const Fragment& fi,
	BSSRDFSampleAxis mainAxis, float pdf, float sigmaTr,
	float Rmax) const {
	float weight = 0.0f;
	// The U, V, N ratio is 1 : 1 : 2, and we do MIS with
	// power heuristic, so 1: 1: 4
	const Vector3& pwo = fo.getPosition();
	const Vector3& pwi = fi.getPosition();
	const Vector3& ni = fi.getNormal();
	if(mainAxis == NAxis) {
	const Vector3& u = normalize(fo.getDPDU());
	const Vector3& v = normalize(fo.getDPDV());
	float UPdf = 0.25f *
	gaussianSample2DPdf(pwo, pwi, u, sigmaTr, Rmax) *
	absdot(u, ni);
	float VPdf = 0.25f *
	gaussianSample2DPdf(pwo, pwi, v, sigmaTr, Rmax) *
	absdot(v, ni);
	float numerator = 4 * pdf * pdf;
	weight = numerator / (numerator + UPdf* UPdf + VPdf * VPdf);
	} else if(mainAxis == UAxis) {
	const Vector3& n = fo.getNormal();
	const Vector3& v = normalize(fo.getDPDV());
	float NPdf = 0.5f *
	gaussianSample2DPdf(pwo, pwi, n, sigmaTr, Rmax) *
	absdot(n, ni);
	float VPdf = 0.25f *
	gaussianSample2DPdf(pwo, pwi, v, sigmaTr, Rmax) *
	absdot(v, ni);
	float numerator = pdf * pdf;
	weight = numerator / (4 * NPdf * NPdf + numerator + VPdf * VPdf);
	} else if(mainAxis == VAxis) {
	const Vector3& n = fo.getNormal();
	const Vector3& u = normalize(fo.getDPDU());
	float NPdf = 0.5f *
	gaussianSample2DPdf(pwo, pwi, n, sigmaTr, Rmax) *
	absdot(n, ni);
	float UPdf = 0.25f *
	gaussianSample2DPdf(pwo, pwi, u, sigmaTr, Rmax) *
	absdot(u, ni);
	float numerator = pdf * pdf;
	weight = numerator / (4 * NPdf * NPdf + UPdf * UPdf + numerator);
	}
	return weight;
	}

view raw mult_axis_MIS.cpp hosted with ❤ by GitHub

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. 

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.

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

BSSRDF Importance Sampling 1 - Kickoff

Recently I start to put in subsurface scattering into my personal side project, and there are some quite fun thoughts and learning experience during the research/implementation time. I decide to write it down as memo and documentation. The following few posts are most likely surround on this topic :)

Why subsurface scattering? Simply feel it looks really cool :P the warm and rich color tone rendered from light bouncing beneath the surface is really different from the regular BRDF material. Theoretically, all material in world is actually BSSRDF instead of BRDF, but BRDF do a good job on approximating those surface light didn't penetrate "deep enough". It would be cool, if I can start to render out some non metal/plastic material, something organic! that's pretty much the original motivation.

After some homework research (mainly the legendary Jensen 01, 02 and Donner's 06 thesis), question pop up: which way should I use? Both Pixar and DreamWorks seems using Jensen 02's point cloud approach from interview and paper release (and my personal working experience= =||), but thinking about managing point cloud and a pre-pass before rendering kinda give me a uncomfortable feeling. Part of my daily job is troubleshooting point cloud artifact and managing its gigantic space usage, and thinking about dealing with it at home after work just feel... not really good....Phar 00's Monte Carlo approach sounds like the most "righteous" method, but it also sounds brutally hardcore :| Jensen 01's Monte Carlo disc sampling sounds like the most straightforward method compare to the previous two...

How about Arnold? The most bad ass MCRT player in the industry at this moment, which method are they using? I googled "Arnold subsurface scattering", then this super awesome slide immediately pop into my eyes!@@ After reading the slides, I felt pretty excited: Cool! They are also using the Jensen 01 approach. Furthermore, the multiple importance sampling method they propose looks brilliant and doesn't seem to be really hard to implement. Alright! I am gonna do thing this way!

nothing relate to SSS in this video actually, but Arnold is kicking some serious ass with EPIC meal time and the awesome sandwich!


Now it's TODO list:

1. I need to have a BSSRDF class pack in relative texture info(absorb, scatter coefficient, refraction index and phase coefficient for single scattering), it should also encapsulate the BSSRDF dipole approximation evaluation.

2. I need one separate BSSRDF integration since the current BRDF integration is apparently doesn't fit with BSSRDF (For direct lighting, BSSRDF needs four dimensional integration on surface area and half sphere solid angle for irradiance, while BRDF only evaluate the later part since all things happen in one intersection point)

3. I need to do importance sampling in 2D with exponential falloff, which is what Jensen01 mentioned and sounds right (it does look like the dominant factor in that dipole equation) and 1D exponential falloff (for single scattering part, this should be pretty straightforward)

4. Figure out the Rmax probe radius that Solid Angle mentioned in above SIGGRAPH13 talk, randomly pick up the probe axis with N:U:V = 2:1:1 ratio they suggested.

5. Implement the multiple scattering part with dipole diffusion approximation (ingredients are listed in the first 4 steps), then combine the sample with multiple importance sampling

6. Implement single scattering part based on Jensen01, Hery03 offered some pseudo code for this, should be straightforward.

And here's the paper references:

BSSRDF Importance Sampling : really brilliant MIS approach for variance noise on complex geometry

A Practical Model For Subsurface Light Transport : the basic implementation I did is a step by step go through on this paper

Implement a Skin BSSRDF (or several) : very throughout tutorial with pseudo code, clarify lots of my unclear thoughts on Jensen01

Towards Realistic Image Synthesis of Scattering Materials : this one got a full derivation of dipole diffusion approximation, and ya....it has some really intense math...

A Rapid Hierarchical Rendering For Translucent Material : the classic point cloud approach that Pixar and DreamWorks use (I believe)

Physically Based Rendering : the ultimate offline rendering strategy guide, though the method in the codebase is based on Jensen 02 approach

Monte Carlo Evaluation on Non-Linear Scattering Equations for Subsurface Reflection : just a really brief and quick pass, but this should be the most physically accurate method (and the slowest...)