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
ω
→
⋅
[
∂
L
(
p
→
,
ω
→
)
∂
x
∂
L
(
p
→
,
ω
→
)
∂
y
∂
L
(
p
→
,
ω
→
)
∂
z
]
=
L
e
(
p
→
,
ω
→
)
−
σ
t
(
p
→
,
ω
→
)
L
(
p
→
,
ω
→
)
+
σ
s
(
p
→
,
ω
→
)
∫
S
2
L
(
p
→
,
ω
′
→
)
p
(
p
→
,
ω
→
,
ω
′
→
)
d
ω
′
→
=
L
e
(
p
→
,
ω
→
)
−
(
σ
a
(
p
→
,
ω
→
)
+
σ
s
(
p
→
,
ω
→
)
)
L
(
p
→
,
ω
→
)
+
σ
s
(
p
→
,
ω
→
)
∫
S
2
L
(
p
→
,
ω
′
→
)
p
(
p
→
,
ω
→
,
ω
′
→
)
d
ω
′
→
\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}'}
ω
⋅ [ ∂ x ∂ L ( p
, ω
) ∂ y ∂ L ( p
, ω
) ∂ z ∂ L ( p
, ω
) ] = L e ( p
, ω
) − σ t ( p
, ω
) L ( p
, ω
) + σ s ( p
, ω
) ∫ S 2 L ( p
, ω ′
) p ( p
, ω
, ω ′
) d ω ′
= L e ( p
, ω
) − ( σ a ( p
, ω
) + σ s ( p
, ω
) ) L ( p
, ω
) + σ s ( p
, ω
) ∫ S 2 L ( p
, ω ′
) p ( p
, ω
, ω ′
) d ω ′
.
TODO
By [Jensen 2001] and "Equation
(15.15)" of PBR Book V3 , we have the
reduced scattering coefficient
σ
s
′
(
p
→
)
=
(
1
−
g
)
σ
s
(
p
→
)
\displaystyle
\operatorname{\sigma_s'}(\overrightarrow{p}) = (1 - g)
\operatorname{\sigma_s}(\overrightarrow{p})
σ s ′ ( p
) = ( 1 − g ) σ s ( p
) , reduced extinction
coefficient
σ
t
′
(
p
→
)
=
σ
a
(
p
→
)
+
σ
s
′
(
p
→
)
\displaystyle
\operatorname{\sigma_t'}(\overrightarrow{p}) =
\operatorname{\sigma_a}(\overrightarrow{p}) +
\operatorname{\sigma_s'}(\overrightarrow{p})
σ t ′ ( p
) = σ a ( p
) + σ s ′ ( p
) and isotropic phase
function
p
(
p
→
,
ω
→
,
ω
′
→
)
=
1
4
π
\displaystyle \operatorname{p}(\overrightarrow{p},
\overrightarrow{\omega}, \overrightarrow{{\omega}'}) = \frac{1}{4 \pi}
p ( p
, ω
, ω ′
) = 4 π 1 . This means that
the RTE can be simplified as
ω
→
⋅
[
∂
L
(
p
→
,
ω
→
)
∂
x
∂
L
(
p
→
,
ω
→
)
∂
y
∂
L
(
p
→
,
ω
→
)
∂
z
]
=
L
e
(
p
→
,
ω
→
)
−
(
σ
a
(
p
→
)
+
σ
s
′
(
p
→
)
)
L
(
p
→
,
ω
→
)
+
σ
s
′
(
p
→
)
4
π
∫
S
2
L
(
p
→
,
ω
′
→
)
d
ω
′
→
\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}'}
ω
⋅ [ ∂ x ∂ L ( p
, ω
) ∂ y ∂ L ( p
, ω
) ∂ z ∂ L ( p
, ω
) ] = L e ( p
, ω
) − ( σ a ( p
) + σ s ′ ( p
) ) L ( p
, ω
) + 4 π σ s ′ ( p
) ∫ S 2 L ( p
, ω ′
) d ω ′
.
Diffusion Approximation
Let
ϕ
(
p
→
)
=
∫
S
2
L
(
p
→
,
ω
→
)
d
ω
→
\displaystyle
\operatorname{\phi}(\overrightarrow{p}) = \int_{\mathrm{S}^2}
\operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega}) \, d
\overrightarrow{\omega}
ϕ ( p
) = ∫ S 2 L ( p
, ω
) d ω
be the fluence rate and
E
→
(
p
→
)
=
[
∫
S
2
L
(
p
→
,
ω
→
)
ω
x
d
ω
→
∫
S
2
L
(
p
→
,
ω
→
)
ω
y
d
ω
→
∫
S
2
L
(
p
→
,
ω
→
)
ω
z
d
ω
→
]
\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}
E
( p
) = [ ∫ S 2 L ( p
, ω
) ω x d ω
∫ S 2 L ( p
, ω
) ω y d ω
∫ S 2 L ( p
, ω
) ω z d ω
] be the
vector irradiance .
By [Jensen 2001] and "Equation
(15.16)" of PBR Book V3 , the unknown
function
L
(
p
→
,
ω
→
)
\displaystyle \operatorname{L}(\overrightarrow{p},
\overrightarrow{\omega})
L ( p
, ω
) can be approximated by the expansion
L
(
p
→
,
ω
→
)
≈
1
4
π
ϕ
(
p
→
)
+
3
4
π
ω
→
⋅
E
→
(
p
→
)
\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})
L ( p
, ω
) ≈ 4 π 1 ϕ ( p
) + 4 π 3 ω
⋅ E
( 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
L
(
p
→
,
ω
→
)
≈
L
0
0
Υ
0
0
(
ω
→
)
+
L
1
−
1
Υ
1
−
1
(
ω
→
)
+
L
1
0
Υ
1
0
(
ω
→
)
+
L
1
−
1
Υ
1
−
1
(
ω
→
)
\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})
L ( p
, ω
) ≈ L 0 0 Υ 0 0 ( ω
) + L 1 − 1 Υ 1 − 1 ( ω
) + L 1 0 Υ 1 0 ( ω
) + L 1 − 1 Υ 1 − 1 ( ω
) where
L
l
m
=
∫
S
2
L
(
p
→
,
ω
→
)
Υ
l
m
(
ω
→
)
d
ω
→
\displaystyle L_l^m = \int_{\mathrm{S}^2}
\operatorname{L}(\overrightarrow{p}, \overrightarrow{\omega})
\operatorname{\Upsilon_l^m}(\overrightarrow{\omega}) \, d \overrightarrow{\omega}
L l m = ∫ S 2 L ( p
, ω
) Υ l m ( ω
) d ω
.
By "Appendix A2" of [Sloan 2008], we have the polynomial forms of basis
function
Υ
0
0
(
ω
→
)
=
1
2
π
\displaystyle
\operatorname{\Upsilon_0^0}(\overrightarrow{\omega}) = \frac{1}{2 \sqrt{\pi}}
Υ 0 0 ( ω
) = 2 π
1 ,
Υ
1
−
1
(
ω
→
)
=
−
3
2
π
ω
y
\displaystyle
\operatorname{\Upsilon_1^{-1}}(\overrightarrow{\omega}) = - \frac{\sqrt{3}}{2
\sqrt{\pi}} \omega_y
Υ 1 − 1 ( ω
) = − 2 π
3
ω y
,
Υ
1
0
(
ω
→
)
=
3
2
π
ω
z
\displaystyle
\operatorname{\Upsilon_1^0}(\overrightarrow{\omega}) = \frac{\sqrt{3}}{2 \sqrt{\pi}}
\omega_z
Υ 1 0 ( ω
) = 2 π
3
ω z
and
Υ
1
1
(
ω
→
)
=
−
3
2
π
ω
x
\displaystyle
\operatorname{\Upsilon_1^1}(\overrightarrow{\omega}) = - \frac{\sqrt{3}}{2
\sqrt{\pi}} \omega_x
Υ 1 1 ( ω
) = − 2 π
3
ω x
.
This means that we have
L
0
0
=
∫
S
2
L
(
p
→
,
ω
→
)
1
2
π
d
ω
→
=
1
2
π
∫
S
2
L
(
p
→
,
ω
→
)
d
ω
→
=
1
2
π
ϕ
(
p
→
)
\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})
L 0 0 = ∫ S 2 L ( p
, ω
) 2 π
1 d ω
= 2 π
1 ∫ S 2 L ( p
, ω
) d ω
= 2 π
1 ϕ ( p
) ,
L
1
−
1
=
∫
S
2
L
(
p
→
,
ω
→
)
(
−
3
2
π
ω
y
)
d
ω
→
=
−
3
2
π
∫
S
2
L
(
p
→
,
ω
→
)
ω
y
d
ω
→
=
−
3
2
π
E
y
(
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})
L 1 − 1 = ∫ S 2 L ( p
, ω
) ( − 2 π
3
ω y ) d ω
= − 2 π
3
∫ S 2 L ( p
, ω
) ω y d ω
= − 2 π
3
E y ( p
) ,
L
1
0
=
∫
S
2
L
(
p
→
,
ω
→
)
3
2
π
ω
z
d
ω
→
=
3
2
π
∫
S
2
L
(
p
→
,
ω
→
)
ω
z
d
ω
→
=
3
2
π
E
z
(
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})
L 1 0 = ∫ S 2 L ( p
, ω
) 2 π
3
ω z d ω
= 2 π
3
∫ S 2 L ( p
, ω
) ω z d ω
= 2 π
3
E z ( p
) and
L
1
1
=
∫
S
2
L
(
p
→
,
ω
→
)
(
−
3
2
π
ω
x
)
d
ω
→
=
−
3
2
π
∫
S
2
L
(
p
→
,
ω
→
)
ω
x
d
ω
→
=
−
3
2
π
E
x
(
p
→
)
\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})
L 1 1 = ∫ S 2 L ( p
, ω
) ( − 2 π
3
ω x ) d ω
= − 2 π
3
∫ S 2 L ( p
, ω
) ω x d ω
= − 2 π
3
E x ( p
) .
And the expansion based on the SH can be calculated as
L
(
p
→
,
ω
→
)
≈
L
0
0
Υ
0
0
(
ω
→
)
+
L
1
−
1
Υ
1
−
1
(
ω
→
)
+
L
1
0
Υ
1
0
(
ω
→
)
+
L
1
−
1
Υ
1
−
1
(
ω
→
)
=
1
2
π
ϕ
(
p
→
)
1
2
π
+
(
−
3
2
π
E
y
(
p
→
)
)
(
−
3
2
π
ω
y
)
+
3
2
π
E
z
(
p
→
)
3
2
π
ω
z
+
(
−
3
2
π
E
x
(
p
→
)
)
(
−
3
2
π
ω
x
)
=
1
2
π
ϕ
(
p
→
)
+
3
4
π
(
E
y
(
p
→
)
ω
y
+
E
z
(
p
→
)
ω
z
+
E
x
(
p
→
)
ω
x
)
=
1
2
π
ϕ
(
p
→
)
+
3
4
π
[
ω
x
ω
y
ω
z
]
[
E
x
(
p
→
)
E
y
(
p
→
)
E
z
(
p
→
)
]
=
1
4
π
ϕ
(
p
→
)
+
3
4
π
ω
→
⋅
E
→
(
p
→
)
\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})
L ( p
, ω
) ≈ L 0 0 Υ 0 0 ( ω
) + L 1 − 1 Υ 1 − 1 ( ω
) + L 1 0 Υ 1 0 ( ω
) + L 1 − 1 Υ 1 − 1 ( ω
) = 2 π
1 ϕ ( p
) 2 π
1 + ( − 2 π
3
E y ( p
) ) ( − 2 π
3
ω y ) + 2 π
3
E z ( p
) 2 π
3
ω z + ( − 2 π
3
E x ( p
) ) ( − 2 π
3
ω x ) = 2 π 1 ϕ ( p
) + 4 π 3 ( E y ( p
) ω y + E z ( p
) ω z + E x ( p
) ω x ) = 2 π 1 ϕ ( p
) + 4 π 3 [ ω x ω y ω z ] [ E x ( p
) E y ( p
) E z ( p
) ] = 4 π 1 ϕ ( p
) + 4 π 3 ω
⋅ E
( p
) .
Diffusion Equation
Let
Q
0
(
p
→
)
=
∫
S
2
L
e
(
p
→
,
ω
→
)
d
ω
→
\displaystyle
\operatorname{Q_0}(\overrightarrow{p}) = \int_{\mathrm{S}^2}
\operatorname{L_e}(\overrightarrow{p}, \overrightarrow{\omega}) \, d
\overrightarrow{\omega}
Q 0 ( p
) = ∫ S 2 L e ( p
, ω
) d ω
be the zeroth moment of the emission and
Q
1
→
(
p
→
)
=
[
∫
S
2
L
e
(
p
→
,
ω
→
)
ω
x
d
ω
→
∫
S
2
L
e
(
p
→
,
ω
→
)
ω
y
d
ω
→
∫
S
2
L
e
(
p
→
,
ω
→
)
ω
z
d
ω
→
]
\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}
Q 1
( p
) = [ ∫ S 2 L e ( p
, ω
) ω x d ω
∫ S 2 L e ( p
, ω
) ω y d ω
∫ S 2 L e ( p
, ω
) ω z d ω
] 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
∂
E
→
(
p
→
)
∂
x
+
∂
E
→
(
p
→
)
∂
y
+
∂
E
→
(
p
→
)
∂
z
=
−
σ
a
(
p
→
,
ω
→
)
ϕ
(
p
→
)
+
Q
0
(
p
→
)
\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})
∂ x ∂ E
( p
) + ∂ y ∂ E
( p
) + ∂ z ∂ E
( p
) = − σ a ( p
, ω
) ϕ ( p
) + Q 0 ( p
) .
By "Equation (2)" of [Jensen 2001] and "Equation
(15.18)" of PBR Book V3 , by substituting
the approximatation into the RTE, we have
[
∂
ϕ
(
p
→
)
∂
x
∂
ϕ
(
p
→
)
∂
y
∂
ϕ
(
p
→
)
∂
z
]
=
−
3
(
σ
a
(
p
→
,
ω
→
)
+
σ
s
′
(
p
→
,
ω
→
)
)
E
→
(
p
→
)
+
Q
1
→
\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}}
[ ∂ x ∂ ϕ ( p
) ∂ y ∂ ϕ ( p
) ∂ z ∂ ϕ ( p
) ] = − 3 ( σ a ( p
, ω
) + σ s ′ ( p
, ω
) ) E
( p
) + 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
D
(
p
→
)
(
∂
2
ϕ
(
p
→
)
∂
2
x
+
∂
2
ϕ
(
p
→
)
∂
2
y
+
∂
2
ϕ
(
p
→
)
∂
2
z
)
=
σ
a
(
p
→
)
ϕ
(
p
→
)
−
Q
0
(
p
→
)
+
3
D
(
p
→
)
(
∂
Q
1
→
(
p
→
)
∂
x
+
∂
Q
1
→
(
p
→
)
∂
y
+
∂
Q
1
→
(
p
→
)
∂
z
)
\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
D ( p
) ( ∂ 2 x ∂ 2 ϕ ( p
) + ∂ 2 y ∂ 2 ϕ ( p
) + ∂ 2 z ∂ 2 ϕ ( p
) ) = σ a ( p
) ϕ ( p
) − Q 0 ( p
) + 3 D ( p
) ( ∂ x ∂ Q 1
( p
) + ∂ y ∂ Q 1
( p
) + ∂ z ∂ Q 1
( p
) ) where
D
(
p
→
)
=
1
3
(
σ
a
(
p
→
)
+
σ
s
′
(
p
→
)
)
\displaystyle \operatorname{D}(\overrightarrow{p})
= \frac{1}{3 \left\lparen \operatorname{\sigma_a}(\overrightarrow{p}) +
\operatorname{\sigma_s'}(\overrightarrow{p}) \right\rparen}
D ( p
) = 3 ( σ a ( p
) + σ s ′ ( p
) ) 1 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.