Diffusion Profile

Diffusion Profile

RTE (Radiative Transfer Equation)

By [Christensen 2015] and [Golubev 2018], the popular disney diffusion profile, which is used by RadiusRootFindAnalytic in UE4 and SampleBurleyDiffusionProfile in Unity3D, is empirical. But, by "15.5.2 Diffusion Theory" of PBR Book V3, technically, the diffusion profile is derived from the RTE (Radiative Transfer Equation).

By "15.5.2 Diffusion Theory" and "11.1.3 Out-Scattering and Attenuation" of PBR Book V3, we have the RTE $\displaystyle \overrightarrow{\omega} \cdot \begin{bmatrix} \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial x} & \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial y} & \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial z} \end{bmatrix} = \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) - \operatorname{\sigma_t}(\overrightarrow{p}, \overrightarrow{\omega}) \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) + \operatorname{\sigma_s}(\overrightarrow{p}, \overrightarrow{\omega}) \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{{\omega}'}) \operatorname{p}(\overrightarrow{p}, \overrightarrow{\omega}, \overrightarrow{{\omega}'}) \, d \overrightarrow{{\omega}'} = \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) - \left\lparen \operatorname{\sigma_a}(\overrightarrow{p}, \overrightarrow{\omega}) + \operatorname{\sigma_s}(\overrightarrow{p}, \overrightarrow{\omega}) \right\rparen \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) + \operatorname{\sigma_s}(\overrightarrow{p}, \overrightarrow{\omega}) \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{{\omega}'}) \operatorname{p}(\overrightarrow{p}, \overrightarrow{\omega}, \overrightarrow{{\omega}'}) \, d \overrightarrow{{\omega}'}$ .

TODO

By [Jensen 2001] and "Equation (15.15)" of PBR Book V3, we have the reduced scattering coefficient $\displaystyle \operatorname{\sigma_s'}(\overrightarrow{p}) = (1 - g) \operatorname{\sigma_s}(\overrightarrow{p})$, reduced extinction coefficient $\displaystyle \operatorname{\sigma_t'}(\overrightarrow{p}) = \operatorname{\sigma_a}(\overrightarrow{p}) + \operatorname{\sigma_s'}(\overrightarrow{p})$ and isotropic phase function $\displaystyle \operatorname{p}(\overrightarrow{p}, \overrightarrow{\omega}, \overrightarrow{{\omega}'}) = \frac{1}{4 \pi}$. This means that the RTE can be simplified as $\displaystyle \overrightarrow{\omega} \cdot \begin{bmatrix} \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial x} & \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial y} & \displaystyle \frac{\partial \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})}{\partial z} \end{bmatrix} = \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) - \left\lparen \operatorname{\sigma_a}(\overrightarrow{p}) + \operatorname{\sigma_s'}(\overrightarrow{p}) \right\rparen \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) + \frac{\operatorname{\sigma_s'}(\overrightarrow{p})}{4 \pi} \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{{\omega}'}) \, d \overrightarrow{{\omega}'}$ .

Diffusion Approximation

Let $\displaystyle \operatorname{\phi}(\overrightarrow{p}) = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \, d \overrightarrow{\omega}$ be the fluence rate and $\displaystyle \operatorname{\overrightarrow{E}}(\overrightarrow{p}) = \begin{bmatrix} \displaystyle \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_x \, d \overrightarrow{\omega} & \displaystyle \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_y \, d \overrightarrow{\omega} & \displaystyle \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_z \, d \overrightarrow{\omega} \end{bmatrix}$ be the vector irradiance.

By [Jensen 2001] and "Equation (15.16)" of PBR Book V3, the unknown function $\displaystyle \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})$ can be approximated by the expansion $\displaystyle \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \approx \frac{1}{4 \pi} \operatorname{\phi}(\overrightarrow{p}) + \frac{3}{4 \pi} \overrightarrow{\omega} \cdot \operatorname{\overrightarrow{E}}(\overrightarrow{p})$ .

Proof

Since the SH(Spherical Harmonics) were NOT popular, this expansion was based on the spherical moments by [Jensen 2001].

However, it is much more convenient to prove this expansion based on the SH.

By "Projection and Reconstruction" of [Sloan 2008], we have the expansion based on the SH $\displaystyle \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \approx L_0^0 \operatorname{\Upsilon_0^0}(\overrightarrow{\omega}) + L_1^{-1} \operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega}) + L_1^0 \operatorname{\Upsilon_1^0}(\overrightarrow{\omega}) + L_1^{-1} \operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega})$ where $\displaystyle L_l^m = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \operatorname{\Upsilon_l^m}(\overrightarrow{\omega}) \, d \overrightarrow{\omega}$ .

By "Appendix A2" of [Sloan 2008], we have the polynomial forms of basis function $\displaystyle \operatorname{\Upsilon_0^0}(\overrightarrow{\omega}) = \frac{1}{2 \sqrt{\pi}}$, $\displaystyle \operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega}) = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_y$ , $\displaystyle \operatorname{\Upsilon_1^0}(\overrightarrow{\omega}) = \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_z$ and $\displaystyle \operatorname{\Upsilon_1^1}(\overrightarrow{\omega}) = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_x$ .

This means that we have $\displaystyle L_0^0 = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \frac{1}{2 \sqrt{\pi}} \, d \overrightarrow{\omega} = \frac{1}{2 \sqrt{\pi}} \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \, d \overrightarrow{\omega} = \frac{1}{2 \sqrt{\pi}} \operatorname{\phi}(\overrightarrow{p})$, $\displaystyle L_1^{-1} = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_y \right\rparen \, d \overrightarrow{\omega} = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_y \, d \overrightarrow{\omega} = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_y}(\overrightarrow{p})$, $\displaystyle L_1^0 = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_z \, d \overrightarrow{\omega} = \frac{\sqrt{3}}{2 \sqrt{\pi}} \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_z \, d \overrightarrow{\omega} = \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_z}(\overrightarrow{p})$ and $\displaystyle L_1^1 = \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_x \right\rparen \, d \overrightarrow{\omega} = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \int_{\mathrm{S}^2} \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_x \, d \overrightarrow{\omega} = - \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_x}(\overrightarrow{p})$.

And the expansion based on the SH can be calculated as $\displaystyle \operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \approx L_0^0 \operatorname{\Upsilon_0^0}(\overrightarrow{\omega}) + L_1^{-1} \operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega}) + L_1^0 \operatorname{\Upsilon_1^0}(\overrightarrow{\omega}) + L_1^{-1} \operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega}) = \frac{1}{2 \sqrt{\pi}} \operatorname{\phi}(\overrightarrow{p}) \frac{1}{2 \sqrt{\pi}} + \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_y}(\overrightarrow{p}) \right\rparen \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_y \right\rparen + \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_z}(\overrightarrow{p}) \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_z + \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \operatorname{E_x}(\overrightarrow{p}) \right\rparen \left\lparen - \frac{\sqrt{3}}{2 \sqrt{\pi}} \omega_x \right\rparen = \frac{1}{2 \pi} \operatorname{\phi}(\overrightarrow{p}) + \frac{3}{4 \pi} \left\lparen \operatorname{E_y}(\overrightarrow{p}) \omega_y + \operatorname{E_z}(\overrightarrow{p}) \omega_z + \operatorname{E_x}(\overrightarrow{p}) \omega_x \right\rparen = \frac{1}{2 \pi} \operatorname{\phi}(\overrightarrow{p}) + \frac{3}{4 \pi} \begin{bmatrix} \displaystyle \omega_x & \displaystyle \omega_y & \displaystyle \omega_z \end{bmatrix} \begin{bmatrix} \displaystyle \operatorname{E_x}(\overrightarrow{p}) & \displaystyle \operatorname{E_y}(\overrightarrow{p}) & \displaystyle \operatorname{E_z}(\overrightarrow{p}) \end{bmatrix} = \frac{1}{4 \pi} \operatorname{\phi}(\overrightarrow{p}) + \frac{3}{4 \pi} \overrightarrow{\omega} \cdot \operatorname{\overrightarrow{E}}(\overrightarrow{p})$.

Diffusion Equation

Let $\displaystyle \operatorname{Q_0}(\overrightarrow{p}) = \int_{\mathrm{S}^2} \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) \, d \overrightarrow{\omega}$ be the zeroth moment of the emission and $\displaystyle \operatorname{\overrightarrow{Q_1}}(\overrightarrow{p}) = \begin{bmatrix} \displaystyle \int_{\mathrm{S}^2} \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_x \, d \overrightarrow{\omega} & \displaystyle \int_{\mathrm{S}^2} \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_y \, d \overrightarrow{\omega} & \displaystyle \int_{\mathrm{S}^2} \operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) \omega_z \, d \overrightarrow{\omega} \end{bmatrix}$ be the first moment of the emission.

By "Equation (1)" of [Jensen 2001] and "Equation (15.17)" of PBR Book V3, by integrating the RTE over all directions $\displaystyle \overrightarrow{\omega}$ , we have $\displaystyle \frac{\partial \operatorname{\overrightarrow{E}}(\overrightarrow{p})}{\partial x} + \frac{\partial \operatorname{\overrightarrow{E}}(\overrightarrow{p})}{\partial y} + \frac{\partial \operatorname{\overrightarrow{E}}(\overrightarrow{p})}{\partial z} = -\operatorname{\sigma_a}(\overrightarrow{p}, \overrightarrow{\omega}) \operatorname{\phi}(\overrightarrow{p}) + \operatorname{Q_0}(\overrightarrow{p})$.

By "Equation (2)" of [Jensen 2001] and "Equation (15.18)" of PBR Book V3, by substituting the approximatation into the RTE, we have $\displaystyle \begin{bmatrix} \displaystyle \frac{\partial \operatorname{\phi}(\overrightarrow{p})}{\partial x} & \displaystyle \frac{\partial \operatorname{\phi}(\overrightarrow{p})}{\partial y} & \displaystyle \frac{\partial \operatorname{\phi}(\overrightarrow{p})}{\partial z} \end{bmatrix} = -3 \left\lparen \operatorname{\sigma_a}(\overrightarrow{p}, \overrightarrow{\omega}) + \operatorname{\sigma_s'}(\overrightarrow{p}, \overrightarrow{\omega}) \right\rparen \operatorname{\overrightarrow{E}}(\overrightarrow{p}) + \operatorname{\overrightarrow{Q_1}}$ .

By [Jensen 2001] and "15.5.2 Diffusion Theory" of PBR Book V3, by combining the previous two equations, we have the classic diffusion equation $\displaystyle \operatorname{D}(\overrightarrow{p}) \left\lparen \frac{\partial^2 \operatorname{\phi}(\overrightarrow{p})}{\partial^2 x} + \frac{\partial^2 \operatorname{\phi}(\overrightarrow{p})}{\partial^2 y} + \frac{\partial^2 \operatorname{\phi}(\overrightarrow{p})}{\partial^2 z} \right\rparen = \operatorname{\sigma_a}(\overrightarrow{p}) \operatorname{\phi}(\overrightarrow{p}) - \operatorname{Q_0}(\overrightarrow{p}) + 3 \operatorname{D}(\overrightarrow{p}) \left\lparen \frac{\partial \operatorname{\overrightarrow{Q_1}}(\overrightarrow{p})}{\partial x} + \frac{\partial \operatorname{\overrightarrow{Q_1}}(\overrightarrow{p})}{\partial y} + \frac{\partial \operatorname{\overrightarrow{Q_1}}(\overrightarrow{p})}{\partial z} \right\rparen$ where $\displaystyle \operatorname{D}(\overrightarrow{p}) = \frac{1}{3 \left\lparen \operatorname{\sigma_a}(\overrightarrow{p}) + \operatorname{\sigma_s'}(\overrightarrow{p}) \right\rparen}$ is the classic diffusion coefficient.

Monopole

DOM (Discrete Ordinates Method)

angular discretization: integro-differential RTE - DOM -> simultaneous partial differential equations

spatial discretization: simultaneous partial differential equations - FEM -> simultaneous linear equations

By [Fattal 2009] and "Further Reading" of "15 Light Transport II: Volume Rendering" of PBR Book V3, the DOM (Discrete Ordinates Method) can also be used to solve the RTE.

References

[Jensen 2001] Henrik Jensen, Steve Marschner, Marc Levoy, Pat Hanrahan. "A Practical Model for Subsurface Light Transport." SIGGRAPH 2001.
[Donner 2005] Craig Donner, Henrik Wann Jensen. "Light Diffusion in Multi-Layered Translucent Materials." SIGGRAPH 2005.
[Sloan 2008] Peter-Pike Sloan. "Stupid Spherical Harmonics (SH) Tricks." GDC 2008.
[Fattal 2009] Raanan Fattal. "Participating Media Illumination using Light Propagation Maps." TOG 2009.
[Arbree 2010] Adam Arbree. "Heterogeneous Subsurface Scattering Using the Finite Element Method." TVCG 2010.
[Christensen 2015] Per Christensen, Brent Burley. "Approximate Reflectance Profiles for Efficient Subsurface Scattering." SIGGRAPH 2015.
[Golubev 2018] Evgenii Golubev. "Efficient Screen-Space Subsurface Scattering Using Burley's Normalized Diffusion in Real-Time." SIGGRAPH 2018.