We present a new BSSRDF for rendering images of translucent materials. Previous diffusion BSSRDFs are limited by the accuracy of classical diffusion theory. We introduce a modified diffusion theory that is more accurate for highly absorbing materials and near the point of illumination. The new diffusion solution accurately decouples single and multiple scattering. We then derive a novel, analytic, extended-source solution to the multilayer search-light problem by quantizing the diffusion Green's function. This allows the application of the diffusion multipole model to material layers several orders of magnitude thinner than previously possible and creates accurate results under high-frequency illumination. Quantized diffusion provides both a new physical foundation and a variable-accuracy construction method for sum-of-Gaussians BSSRDFs, which have many useful properties for efficient rendering and appearance capture. Our BSSRDF maps directly to previous real-time rendering algorithms. For film production rendering, we propose several improvements to previous hierarchical point cloud algorithms by introducing a new radial-binning data structure and a doubly-adaptive traversal strategy.


Eugene D'Eon;  Geoffrey Irving

Publons users who've claimed - I am an author
Contributors on Publons
  • 1 author
  • 2 reviewers
Followers on Publons
  • Following presentation of the paper at SIGGRAPH 2011 I can only recall one question from Wenzel Jakob: "How were the parameters of the BSSRDF derived when rendering a textured surface, like the face images?" Response: A pre-computed lookup table is used to map a requested diffuse albedo to a single-scattering albedo such that normally incident illumination of a semi-infinite material produces the desired diffuse albedo. For the skin renders, this was done for a fixed index of refraction (1.4). How this is applied during rendering using the radial binning acceleration method is explained further in the paper.

    Reviewed by
    Ongoing discussion
  • On page 4 the footnote incorrectly states that equation (5) is not absolutely convergent. Equation (5) provides the convergent form that separates the ballistic fluence from the scattered fluence and does converge. It is the alternative form $$ \phi(r) = \frac{\mu_t^2}{2 \pi^2 r} \int_0^\infty \frac{\text{arctan}\, u}{\mu_t-\mu_s \, \frac{\text{arctan} \, u}{u}} \, \sin (r \, \mu_t \, u) \, du. $$ that does not separate the ballistic fluence explicitly and only converges in the sense of Cesaro summability.

    Reviewed by
    Ongoing discussion
  • Let $\tau_0:=\frac{1}{s}\tau_1$, then $\tau_i=s^{i-1}\tau_1$ and $v_i=D(\tau_{i+1}+\tau_i)=s^i v_0$ $$\int\limits_{0}^{\infty} G_{3D}(2D\tau,r)e^{-\mu_a\tau}d\tau \approx \sum\limits_{i=0}^{k-1}\int\limits_{\tau_i}^{\tau_{i+1}} G_{3D}(2D\tau,r)e^{-\mu_a\tau}d\tau = \sum\limits_{i=0}^{k-1}(\tau_{i+1}-\tau_i) G_{3D}(D(\tau_{i+1}+\tau_i),r)e^{-\mu_a\frac{\tau_{i+1}+\tau_i}{2}} = \frac{s-1}{s+1}\sum\limits_{i=0}^{k-1}(\tau_{i+1}+\tau_i) G_{3D}(D(\tau_{i+1}+\tau_i),r)e^{-\mu_a\frac{\tau_{i+1}+\tau_i}{2}} = \frac{s-1}{s+1}\sum\limits_{i=0}^{k-1}\frac{s^i v_0}{D} G_{3D}(s^i v_0,r)e^{-\mu_a\frac{s^i v_0}{2D}}$$

    Given $s=\frac{1+\sqrt{5}}{2}$, $ \frac{s-1}{s+1} \approx 0.236068$ rather than $0.240606$ in the paper. Can someone tell me how to get that magic number? And why $v_0$ is specified in terms of $\mu_a$? Thanks for your kindness.

    Reviewed by
    Ongoing discussion (2 comments - click to toggle)
    • Eugene d'Eon | 6 years, 11 months ago

      Excellent questions, G Zhou, thanks for asking them. Second question first---the dependence of $v_0$ on $\mu_a$ arises for two reasons. The range of times we need to consider to describe all the possible photons is $t \in [0,\infty]$ and because we quantize the time variable with an ever-widening set of time intervals (each one $s > 1$ times wider than the last) we, in theory, require an infinite number of time packets to account for all the photons and might consider an expansion $$\phi(r) \approx \sum_{i=-\infty}^\infty w(s^i) G_{3D}(v(s^i),r)$$ where now the exponent on the spread factor $s$ can be negative as well as positive. However, in practice we can allow this sum to become finite because, for very negative $i$, $s^i$ becomes negligible, and because absorption (if non-zero) will make the terms with large $i$ negligible. With a finite range of variances, we then seek a range of $50$ or so that encompasses most of the energy (relative to the total energy reflecting from a semi-infinite half-space). Because of how absorption dampens higher order terms (wider Gaussians), the lower the aborption level, the more energy we will see appearing in very large Gaussian terms, and the relative importance of the very narrow ones goes down. To keep our budget of $50$ Gaussians describing the most important parts of the signal as the absorption level goes down, we shift the variance of the smallest in our set, neglecting some minimial amounts of energy, and leaving more Gaussians for the very wide ones that become increasingly important. The absorption coefficient $\mu_a$ varies with absorption level, so we see it appear in the selection of $v_0$. Also, the absorption coefficient is $a \mu_t$, where $a$ is the single-scattering albedo and $\mu_t$ is the inverse of the mean-free path---the fundamental length scale of the transport problem. If we scale $\mu_t$ and move all the scattering events closer together or farther apart, all the Gaussians in our Green's function need to scale with the mean-free path, and so the variances change with $1/\mu_t^2$. For both of these reasons the optimal selection of $v_0$ depends on $\mu_a$.

      Regarding your first question---the derivation of (21). We described one simple method for producing temporally quantized Greens function expansions, but the explicit form we provided in (21) has an additional optimization step applied that we did not describe in the paper. For the following, without loss of generality, let $\mu_t = 1, \mu_s = a, \mu_a = 1-a$, where $a$ is the single-scattering albedo. We noted that the simple quantization method produced a sum of the form $$ \phi(r) \approx M \sum_{i=-\infty}^{\infty} \frac{s^i e^{-\frac{(1-a) s^i}{2 D}}}{D} G_{3D}(s^i,r) $$ where $M$ is some constant that is independent of $a$. To require that the approximate expansion satisfy a moment relation that the exact Green's function also satisifes, we required that the coefficients of $a^k$ in the Maclaurin series expansion of the zeroth moment $$ \int_0^{\infty} 4 \pi r^2 \phi(r) dr = \frac{a}{1-a} = a + a^2 + a^3 + a^4 + ... $$ be $1$. Given that the spherical Gaussians have zeroth moments of $1$, it is possible to form the series expansion in $a$ of the approximate summation and derive that this moment requirement reduces to setting $M$ to $$ M = \frac{1}{ \sum _{i=-\infty }^{\infty } \frac{3}{2} e^{-\frac{3 s^i}{4}} s^i } $$ which for $s$ equal to the golden ratio gives the magic constant $M = 0.240605907165256938044$, and overall outperforms the more simple expansion derivation.

    • Guo Zhou | 6 years, 11 months ago

      Hi Dr. d'Eon,

      Many thanks for shedding me so much light. Please let me insist to put it in a nutshell.

      So to satisfy the moment relation (equation 10) $$\frac{s-1}{s+1}\sum\limits_{i=-\infty}^{\infty}\frac{s^i v_0}{D} (\int_{0}^{\infty}G_{3D}(s^i v_0,r)dr)e^{-\mu_a\frac{s^i v_0}{2D}}=\frac{s-1}{s+1}\sum\limits_{i=-\infty}^{\infty}\frac{s^i v_0}{D} e^{-\mu_a\frac{s^i v_0}{2D}}$$ should approximate to $\frac{1}{\mu_a}$ to compensate the quantization error, i.e. $$M:=\frac{1}{\sum\limits_{i=-\infty}^{\infty} \mu_a \frac{s^i v_0}{D} e^{-\mu_a\frac{s^i v_0}{2D}}}$$ $$\frac{M}{\frac{s-1}{s+1}}\frac{s-1}{s+1}\sum\limits_{i=-\infty}^{\infty}\frac{s^i v_0}{D}G_{3D}(s^i v_0,r)e^{-\mu_a\frac{s^i v_0}{2D}}=M \sum\limits_{i=-\infty}^{\infty}\frac{s^i v_0}{D}G_{3D}(s^i v_0,r)e^{-\mu_a\frac{s^i v_0}{2D}}$$ whose moment is just $1$. To get that magic constant it turns out $$\mu_a\frac{v_0}{D}=\mu_a\tau_1(\frac{1}{s}+1)=\frac{3}{2}$$ but how does it come out?

      Moreover, since only finite range of variances is considered, the truncation error is likely to break the moment relation. So it should be compensated as well?



  • The equation between (19) and (20) has incorrect indices on the $\tau$s and should read:

    $$ w_i = \int_{\tau_i}^{\tau_{i+1}} e^{-\tau \mu_a} d\tau = \frac{e^{-\tau_i \mu_a} - e^{-\tau_{i+1} \mu_a }}{\mu_a} $$

    (thanks to Toshiya Hachisuka for bringing this to my attention)

    Reviewed by
    Ongoing discussion
All peer review content displayed here is covered by a Creative Commons CC BY 4.0 license.