# Sampling Burley's Normalized Diffusion Profiles

A couple of years ago, I worked on an implementation of Burley's Normalized Diffusion (a.k.a. Disney SSS). The original paper claims that the CDF is not analytically invertible. I have great respect for both authors, Brent Burley and Per Christensen, so I haven't questioned their claim for a second. Turns out, "Question Everything" is probably a better mindset.

I've been recently alerted by @stirners_ghost on Twitter (thank you!) that the CDF is, in fact, analytically invertible. Moreover, the inversion process is almost trivial, as I will demonstrate below.

The diffusion model is radially symmetric, and is defined as

$$\tag{1} R(r) = \frac{A s}{8 \pi r} \Big( e^{-s r} + e^{-s r / 3} \Big).$$

To find the corresponding PDF, it must be normalized:

$$\tag{2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{R(\sqrt{x^2+y^2})}{A} dx dy = \int_{0}^{2 \pi} \int_{0}^{\infty} \frac{R(r)}{A} r dr d\theta = \int_{0}^{2 \pi} \int_{0}^{\infty} p_r(r, \theta) dr d\theta = 1.$$

Therefore, in polar coordinates, the PDF is given as

$$\tag{3} p_r(r, \theta) = \frac{R(r)}{A} r = \frac{s}{8 \pi} \Big( e^{-s r} + e^{-s r / 3} \Big).$$

The marginal CDF is then

$$\tag{4} P(r) = \int_{0}^{2 \pi} \int_{0}^{r} p_r(z, \theta) dz d\theta = 1 - \frac{1}{4} e^{-s r} - \frac{3}{4} e^{-s r / 3}.$$

Let us define $$x = s r$$. If we ask Mathematica to solve $$y = P(x)$$ for $$x$$ in order to obtain the inverse $$x = P^{-1}(y)$$, we get a somewhat terrifying (but, of course, correct) output:

The trick is to solve $$P(x) = 1 - u$$ instead (where $$u = 1 - P(x)$$ is the complementary CDF):

To get rid of the imaginary part, we can set the free constant $$c_{1} = 0$$, which results in the following inverse:

$$\tag{5} x = s r = 3 \log{\Bigg(\frac{1 + G(u)^{-1/3} + G(u)^{1/3}}{4 u} \Bigg)},$$

where

$$\tag{6} G(u) = 1 + 4 u \Big( 2 u + \sqrt{1 + 4 u^2} \Big).$$

We can verify that it works by feeding the complementary CDF with its inverse.

For importance sampling, we can uniformly sample either the complementary or the regular CDF - it makes no difference $$($$except for reversing the order of samples, s.t. $$x(0) = \infty$$ and $$x(1) = 0 )$$.

An optimized implementation is listed below.

// Performs sampling of a Normalized Burley diffusion profile in polar coordinates.
// 'u' is the random number (the value of the CDF): [0, 1).
// rcp(s) = 1 / ShapeParam = ScatteringDistance.
// 'r' is the sampled radial distance, s.t. (u = 0 -> r = 0) and (u = 1 -> r = Inf).
// rcp(Pdf) is the reciprocal of the corresponding PDF value.
void SampleBurleyDiffusionProfile(float u, float rcpS, out float r, out float rcpPdf)
{
u = 1 - u; // Convert CDF to CCDF; the resulting value of (u != 0)

float g = 1 + (4 * u) * (2 * u + sqrt(1 + (4 * u) * u));
float n = exp2(log2(g) * (-1.0/3.0));                    // g^(-1/3)
float p = (g * n) * n;                                   // g^(+1/3)
float c = 1 + p + n;                                     // 1 + g^(+1/3) + g^(-1/3)
float x = (3 / LOG2_E) * log2(c * rcp(4 * u));           // 3 * Log[c / (4 * u)]

// x      = s * r
// exp_13 = Exp[-x/3] = Exp[-1/3 * 3 * Log[c / (4 * u)]]
// exp_13 = Exp[-Log[c / (4 * u)]] = (4 * u) / c
// exp_1  = Exp[-x] = exp_13 * exp_13 * exp_13
// expSum = exp_1 + exp_13 = exp_13 * (1 + exp_13 * exp_13)
// rcpExp = rcp(expSum) = c^3 / ((4 * u) * (c^2 + 16 * u^2))
float rcpExp = ((c * c) * c) * rcp((4 * u) * ((c * c) + (4 * u) * (4 * u)));

r      = x * rcpS;
rcpPdf = (8 * PI * rcpS) * rcpExp; // (8 * Pi) / s / (Exp[-s * r / 3] + Exp[-s * r])
}


## Acknowledgments

I would like to thank @stirners_ghost for informing me that the CDF is invertible, and Brent Burley for spotting a missing minus sign in the derivation.