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.