The name environment lighting is from "10.2 Environment Lighting" of Real-Time Rendering Fourth Edition. By "10.2 Environment
Lighting" of Real-Time Rendering Fourth Edition,
"14.4 The Light Transport Equation" of PBR
Book V3 and "13.1 The Light Transport Equation" of PBR Book
V4, the most important difference between environment lighting (local illumination) and global
illumination is that the shading algorithm of the environment lighting (local illumination) is independent of
the other positions on the surface except the shading position. Hence, when we are discussing the environment
lighting, the position parameter p of most functions, such as f(p,ωi,ωo), Lo(p,ωo) and Li(p,ωi), is omitted.
Let L(ω) be the distant radiance distribution of
which the domain is in world space.
Let n be the normal in world space and ωi be the incident direction in tangent space (where the normal direction is the Z axis [001]).
Evidently, the incident direction in world space can be calculated as R(n)ωi where R(n) is the rotation matrix depending on the
normal direction n in world space. And we have the incident radiance Li(ωi)=L(R(n)ωi).
Since the Lambert BRDF f(ωi,ωo)=π1ρss is constant, we have Lo(ωo)=∫S2f(ωi,ωo)Li(ωi)(cosθi)+dωi=π1ρss⋅E(n) where E(n)=∫S2Li(ωi)(cosθi)+dωi=∫S2L(R(n)ωi)(cosθi)+dωi.
1-1. SH (Spherical Harmonics)
1-1-1. SH Basis
Let Υlm(ω) be the SH (Spherical
Harmonics) basis function of which l is the degree (or band) and m is the basis function index from
-l to l.
The SH basis function Υl0(ω) of which the m is zero is also called
ZH (zonal harmonics).
By "Appendix A2" of [Sloan 2008], we have the polynomial forms of SH basis
Υlm(ω). These polynomial forms are calculated by
sh_eval_basis_2
in DirectXMath, and SHBasisFunction3
in UE4.
Note that the direction vector ω=[xyz] should be
normalized before using the polynomial forms. By "13.5.3 Spherical Coordinates" of PBR
Book V3 and "Spherical Coordinates" of "3.8.3 Spherical Parameterizations" of PBR
Book V4, we have ω=[xyz]=[sinθcosϕsinθsinϕcosθ] where
ϕ is azimuth and θ is zenith (in physics,
while in mathmatics θ is the azimuth and ϕ is the zenith).
Let R be the rotation matrix. By "A.
Rotation of Spherical Harmonics" of "4. SPHERICAL HARMONIC REPRESENTATION" of [Ramamoorthi 2001
A], for each degree (or band) l, we have Υlm(Rω)=j=−l∑lDmjl(R)Υlj(ω) where Dl(R) is the Wigner D-matrix.
This means that ⎣⎡Υl−l(Rω)⋮Υl0(Rω)⋮Υll(Rω)⎦⎤=Dl(R)⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤.
By "Appendix: SH Rotation" of [Kautz 2002], for each degree (or band) l, each element of the
Wigner D-matrix can be calculated as Dijl(R)=∫S2Υli−l(Rω)Υlj−l(ω)dω.
For l = 0, we have D000=∫S2Υ00(Rω)Υ00(ω)dω=∫S22π12π1dω=4π1∫S21dω=4π14π=1. This means that Dl(R)=[1].
For l = 1, it is too complex to calculate the intergral for each element of the Wigner D-matrix. The
trick by [Hable 2014] can be used to calculate the Wigner D-matrix. Since the equation ⎣⎡Υ1−1(Rω)Υ10(Rω)Υ11(Rω)⎦⎤=Dl(R)⎣⎡Υ1−1(ω)Υ10(ω)Υ11(ω)⎦⎤ holds for arbitrary direction ω, we can apply this equation to three linearly independent vectors ⎣⎡100⎦⎤, ⎣⎡010⎦⎤ and ⎣⎡001⎦⎤ to construct an invertible matrix. By the polynomial forms of SH basis, we have ⎣⎡Υ1−1(Rω)Υ10(Rω)Υ11(Rω)⎦⎤=Dl(R)⎣⎡Υ1−1(ω)Υ10(ω)Υ11(ω)⎦⎤⇒⎣⎡00−2π3−2π30002π30⎦⎤Rω=Dl(R)⎣⎡00−2π3−2π30002π30⎦⎤ω⇒⎣⎡00−2π3−2π30002π30⎦⎤R⎣⎡100010001⎦⎤=Dl(R)⎣⎡00−2π3−2π30002π30⎦⎤⎣⎡100010001⎦⎤⇒⎣⎡00−1−100010⎦⎤R⎣⎡100010001⎦⎤=Dl(R)⎣⎡00−1−100010⎦⎤⎣⎡100010001⎦⎤. Evidently, ⎣⎡00−1−100010⎦⎤⎣⎡100010001⎦⎤ is invertible and we have Dl(R)=⎣⎡00−1−100010⎦⎤R⎣⎡100010001⎦⎤⎝⎛⎣⎡00−1−100010⎦⎤⎣⎡100010001⎦⎤⎠⎞−1=⎣⎡00−1−100010⎦⎤R⎣⎡100010001⎦⎤⎣⎡0−10001−100⎦⎤=⎣⎡00−1−100010⎦⎤R⎣⎡0−10001−100⎦⎤=⎣⎡−R10R20−R00−R11R21−R01−R12R22−R02⎦⎤⎣⎡0−10001−100⎦⎤=⎣⎡R11−R21R01−R12R22−R02R10−R20R00⎦⎤.
Let SH be
the SH (Spherical Harmonics) projection operation. Analogous to the Fourier
transform, we have flm=SH(f(ω))=∫S2f(ω)Υlm(ω)dω, and the original function can be reconstructed as the SH series f(ω)=l=0∑∞m=−l∑lflmΥlm(ω)=l=0∑∞[fl−l⋯fl0⋯fll]⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤.
1-1-4. SH Rotational Invariance
Let R be the rotation matrix. Let f′(ω)=f(Rω) and f′(ω)=l=0∑∞m=−l∑lf′lmΥlm(ω). By "Basic Properties" of
"3. Review of Spherical Harmonics" of [Sloan 2002], we have [f′l−l⋯f′l0⋯f′ll]=[fl−l⋯fl0⋯fll]Dl(R) where Dl(R) is the Wigner D-matrix.
Proof
By "SH Projection", we have f(ω)=l=0∑∞m=−l∑lflmΥlm(ω)=l=0∑∞[fl−l⋯fl0⋯fll]⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤.
By "SH Rotation", we have ⎣⎡Υl−l(Rω)⋮Υl0(Rω)⋮Υll(Rω)⎦⎤=Dl(R)⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤.
This means that f′(ω)=f(Rω)=l=0∑∞[fl−l⋯fl0⋯fll]⎣⎡Υl−l(Rω)⋮Υl0(Rω)⋮Υll(Rω)⎦⎤=l=0∑∞[fl−l⋯fl0⋯fll]Dl(R)⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤ where Dl(R) is the Wigner D-matrix.
Since f′(ω)=l=0∑∞m=−l∑lf′lmΥlm(ω)=l=0∑∞[f′l−l⋯f′l0⋯f′ll]⎣⎡Υl−l(ω)⋮Υl0(ω)⋮Υll(ω)⎦⎤, we have [f′l−l⋯f′l0⋯f′ll]=[fl−l⋯fl0⋯fll]Dl(R).
1-1-5. SH Product Integration
Let f(ω)=l=0∑∞m=−l∑lflmΥlm(ω) and g(ω)=l=0∑∞m=−l∑lglmΥlm(ω). By "Basic Properties" of
"3. Review of Spherical Harmonics" of [Sloan 2002], due to the orthonormality of the SH basis, we have
∫S2f(ω)g(ω)dω=∫S2(l=0∑∞m=−l∑lflmΥlm(ω))(l=0∑∞m=−l∑lglmΥlm(ω))dω=l=0∑∞m=−l∑lflmglm.
1-1-6. SH Product Projection
Actually, "SH Product Projection" is related to the Clebsch–Gordan
coefficients which is too complex to be used in rendering. We only need to know that "SH
Product Projection" should be distinguished from "SH Product Integration".
Let h(ω)=f(ω)g(ω). By "Product Projection" of
"3. Review of Spherical Harmonics" of [Sloan 2002], we have that hlm=∫S2h(ω)Υlm(ω)dω=∫S2f(ω)g(ω)Υlm(ω)dω=∫S2f(ω)g(ω)dω which is totally different from "SH Product Integration".
1-2. SH convolution
Analogous to the convolution theorem, by [Ramamoorthi 2001 A], we have E(n)=∫S2Li(ωi)(cosθi)+dωi=∫S2L(R(n)ωi)(cosθi)+dωi=l=0∑∞m=−l∑l2l+14πLlmAlΥlm(n) where Llm=SH(L(ω))=∫S2L(ω)Υlm(ω)dω and Al=ZH((cosθi)+)=∫S2(cosθi)+Υl0(ωi)dωi.
Proof
By "SH Rotational Invariance", we have L(R(n)ωi)=l=0∑∞[Ll−l⋯Ll0⋯Lll]Dl(R(n))⎣⎡Υl−l(ωi)⋮Υl0(ωi)⋮Υll(ωi)⎦⎤ where Dl(R(n)) is the Wigner D-matrix.
Due to circular symmetry of the clamped cosine, only the coefficients of the ZH(Zonal Harmonics)
are non-zero. By "SH Projection", we have (cosθi)+=l=0∑∞AlΥl0(ωi).
By "SH Product Integration", due to the orthonormality of the SH basis, only the
elements Dm0l(R(n)) of the Wigner D-matrix are multiplied
by the non-zero term, and we have E(n)=∫S2L(R(n)ωi)(cosθi)+dωi=∫S2⎝⎛l=0∑∞[Ll−l⋯Ll0⋯Lll]Dl(R(n))⎣⎡Υl−l(ωi)⋮Υl0(ωi)⋮Υll(ωi)⎦⎤⎠⎞(l=0∑∞AlΥl0(ωi))dωi=l=0∑∞m=−l∑lLlmAlDm0l(R(n)).
By "Equation (23)" of [Ramamoorthi 2001 A], we have Dm0l(R(n))=2l+14πΥlm(n), and we have E(n)=l=0∑∞m=−l∑lLlmAlDm0l(R(n))=l=0∑∞m=−l∑l2l+14πLlmAlΥlm(n).
1-3. Cubemap Texel Solid Angle
We merely use numerical quadrature rather than Monte Carlo
integration to integrate over the cubemap. But it should be noted that the solid angle subtended by
each texel of the cubemap is NOT the same. Let u and v be the texcoord of the cubemap. By "5.5.3 Integrals
over Area" of PBRT
Book V3 and "4.2.3 Integrals over Area" of PBR
Book V4, we have dω=r2dAcosθ=dA⋅cosθ⋅r21=cubesize_u⋅cubesize_v(1−(−1))⋅(1−(−1))⋅12+u2+v21⋅12+u2+v21=cubesize_u⋅cubesize_v1⋅12+u2+v2⋅(12+u2+v2)4. Actually the
pseudo code fWt = 4/(sqrt(fTmp)*fTmp) by "Projection from Cube Maps" of [Sloan 2008] is
exactly the 12+u2+v2⋅(12+u2+v2)4. The common
divisor cubesize_u⋅cubesize_v1 can be reduced,
and thus is NOT calculated by [Sloan 2008]. The solid angle is calculated by SHProjectCubeMap
in DirectXMath and DiffuseIrradianceCopyPS
in UE4.
Here is the MATLAB code which verifies the meaning of "fWt = 4/(sqrt(fTmp)*fTmp)".
% retrieved by "textureSize(GLSL)" or "GetDimensions(HLSL)".
cube_size_u = single(4096.0);
cube_size_v = single(4096.0);
[ index_u, index_v ] = meshgrid((single(0.0) : (cube_size_u - single(1.0))), (single(0.0) : (cube_size_v - single(1.0))));
u = (index_u + single(0.5)) ./ cube_size_u .* single(2.0) - single(1.0);
v = (index_v + single(0.5)) ./ cube_size_u .* single(2.0) - single(1.0);
% the common divisor "1/(cube_size_u*cube_size_v)" can be reduced, and thus is NOT calculated in the "fWt = 4/(sqrt(fTmp)*fTmp)".
d_a = (single(1.0) - single(-1.0)) .* (single(1.0) - single(-1.0)) ./ cube_size_u ./ cube_size_v;
r_2 = single(1.0) .* single(1.0) + single(u) .* single(u) + single(v) .* single(v);
cos_theta = single(1.0) ./ sqrt(r_2);
d_omega = single(d_a) .* single(cos_theta) ./ single(r_2);
% "numerical_omega" and "analytic_omega" are expected to be the same.
numerical_omega = sum(sum(d_omega));
analytical_omega = single(4.0) .* single(pi) .* single(1.0) .* single(1.0) / single(6.0);
% output: "numerical:2.094395 analytic_omega:2.094395".
printf("numerical:%f analytical:%f\n", numerical_omega, analytical_omega);
1-4. Clamped Cosine
By [Ramamoorthi 2001 B], 2l+14πAl is constant and can be pre-integrated. For example, by "5.5.1 Integrals over Projected Solid
Angle" of PBR
Book V3 and "4.2.1 Integrals over Projected Solid Angle" of PBR
Book V4, A0 can pre-integrated as A0=∫S2(cosθ)+Υ00(ω)dω=∫S2(cosθ)+2π1dω=2π1∫S2(cosθ)+dω=2π1dω⊥=2π1π. These constants are
calculated by CalcDiffuseTransferSH3
in UE4.
l
2l+14π
Al
2l+14πAl
0
2π
2π1π
π
1
32π
2π332π
32π
2
52π
2π54π
4π
1-5. Form Factor
"Appendix A10" of [Sloan 2008] proposed a more efficient approach to reconstruct from SH
basis than [Ramamoorthi 2001 B]. This approach is used by SampleSH9
in Unity3D and GetSkySHDiffuse
in UE4.
However, it was form factor rather than irradiance that "Appendix
A10" of [Sloan 2008] calculated. Although the terms irradiance and form
factor may be interchangeably used, technically irradiance should NOT be divided
by π. This means that E(n)=πF(n). This form factor is calculated by PackCoefficients
in Unity3D, and SetupSkyIrradianceEnvironmentMapConstantsFromSkyIrradiance
and ComputeSkyEnvMapDiffuseIrradianceCS
in UE4.
And since the value reconstructed from SH basis is form factor rather than irradiance, it should be multiplied
by albedo rather than Lambert BRDF. This multiplication is calculated by GetDiffuseOrDefaultColor
in Unity3D and EnvBRDFApproxFullyRough
in UE4.
fc
value
fc0
π1Υ00A0=π12π1π=2π1
fc1
π1Υ10A1=π12π332π=3π3
fc2
π1Υ2−2A2=π12π154π=8π15
fc3
π1(−Υ20)A2=π14π54π=16π5
0.3 * fc3
π1Υ20A2=π14π354π=16π35
fc4
π1Υ22A2=π14π154π=16π15
2. Specular Environment Lighting
Let L(ω) be the distant radiance distribution.
Evidently, we have Lo(ωo)=∫S2f(ωi,ωo)Li(ωi)(cosθi)+dωi=∫S2f(ωi,ωo)(cosθi)+dωi∫S2f(ωi,ωo)Li(ωi)(cosθi)+dωi⋅∫S2f(ωi,ωo)(cosθi)+dωi=LD(ωo)⋅DFG(ωo) where LD(ωo)=∫S2f(ωi,ωo)(cosθi)+dωi∫S2f(ωi,ωo)Li(ωi)(cosθi)+dωi and DFG(ωo)=∫S2f(ωi,ωo)(cosθi)+dωi.
2-1. Monte Carlo Integration
2-1-1. PDF - NDF
Let ωh=∥ωi+ωo∥ωi+ωo be the half
vector. By "Equation (5.6)" of PBRT-V3,
"Figure 6" of [Walter 2007] and "Figure 14.4" of PBRT-V3,
we have the relationship between dωh and dωi that dωh=r2cosθdA=∥ωi+ωo∥2∣ωi⋅ωh∣dωi=∥ωi+(2(ωi⋅ωh)ωh−ωi)∥2∣ωi⋅ωh∣dωo=4(ωi⋅ωh)2∥ωh∥2∣ωi⋅ωh∣dωi=4∣ωi⋅ωh∣1dωi. Note that since the subtended area on the sphere surface remains the same, the dωh can be calculated even if the ωh is NOT normalized.
By "Equation (9)" of [Heitz 2014] and "Figure 8.15" of PBRT-V3,
we have ∫ΩD(ωh)(cosθh)+dωh=1. According to the relationship between dωh and dωi, we have ∫ΩD(ωh)(cosθh)+dωh=∫ΩD(ωh)(cosθh)+4∣ωi⋅ωh∣1dωi=1. Evidently, this normalized function can be
used as the PDF p(ωi)=D(ωh)(cosθh)+4∣ωi⋅ωh∣1.
2-1-2. PDF - VNDF
Genuinely, although the NDF is the PDF used by UE4, the VNDF by [Heitz 2018] is the better
alternative.
By "Equation (14)" of [Heitz 2014] and "Equation (8.12)" of PBRT-V3, we
have ∫ΩG1(ωo,ωh)(ωo⋅ωh)+D(ωh)dωh=cosθo.
TODO:
2-1-3. Sampling Trowbridge-Reitz Distribution
- NDF
For Beckmann–Spizzichino ("Equation (8.9)" of PBRT-V3)
distribution, the derivation of the NDF sampling is provided by "Equation (14.1)" of PBRT-V3.
For isotropic Trowbridge-Reitz ("Equation (8.11)" of PBRT-V3)
distribution, namely, GGX distribution, the NDF sampling can be analogously derived. Since the
distribution is isotropic, we have p(ϕ∣θh)=2π1⇒ϕ=2πξ2. Since the PDF is p(ωh)=D(ωh)(cosθh)+, by "13.5.3 Spherical Coordinates" of PBRT-V3,
we have the marginal PDF p(θh)=∫02πp(θh,ϕ)dϕ=∫02πp(ωh)sinθhdϕ=∫02πD(ωh)(cosθh)+sinθhdϕ=2πD(ωh)(cosθh)+sinθh and the CDF P(θh)=∫0θhp(θh′)dθh′=ξ1. Fortunately, the inverse of the CDF is closed-form and we have tan2θh=1−ξ1α2ξ1⇒cosθh=1+tan2θh1=1+(α2−1)ξ11−ξ1.
2-1-4. Sampling Trowbridge-Reitz
Distribution - VNDF
For Beckmann–Spizzichino ("Equation (8.9)" of PBRT-V3)
distribution, the derivation of the VNDF sampling is provided by "Equation (14.2)" of PBRT-V3.
TODO:
2-1-5. Low-Discrepancy Sequence
By "20.3 Quasirandom Low-Discrepancy Sequences" of [Colbert 2007] and "13.8.2 Quasi
Monte Carlo" of PBRT-V3,
the low-discrepancy sequence is the better alternative than pseudo-random
sequence to generate the ξ1 and ξ2.
For example, the Hammersley sequence ("7.4.1 Hammersley and Halton
Sequences" of PBRT-V3)
is a typical low-discrepancy sequence.
By "Equation (9.9)" of Real-Time Rendering
Fourth Edition, the DFG term is the HDR (Hemispherical Directional Reflectance) R(ωo)=∫S2f(ωi,ωo)(cosθi)+dωi=∫S2f(ωi,ωo)((0,0,1)⋅ωi)+dωi. Since this integral can be calculated in the tangent space where the normal direction is (0, 0, 1), this
integral is independent of the normal direction and can be determined by the view direction in tangent space
alone.
And by "Equation (13.3)" of PBRT-V3, we have
R(ωo)=∫S2f(ωi,ωo)(cosθi)+dωi=N1i=1∑Np(ωi,ωo)f(ωi,ωo)(cosθi)+. Evidently, D(ωh) exists in both f(ωi,ωo) and p(ωi,ωo), and thus D(ωh) can be omitted when calculating the
integral.
Actually, by "Equation (9)" of [Karis 2013], the Fresnel term can be treated separately, and
we have R(ωo)=∫S2[F0+(F90−F0)(1.0−ωo⋅ωh)5]DV(ωi,ωo)dωi=F0⋅nR(ωo)+F90⋅nG(ωo) where nR(ωo)=∫S2[1.0−(1.0−ωo⋅ωh)5]DV(ωi,ωo)dωi and nG(ωo)=∫S2(1.0−ωo⋅ωh)5DV(ωi,ωo)dωi. Evidently, nR(ωo) and nG(ωo) can be stored in the LUT.
By "Equation (13.3)" of PBRT-V3, we have
LD(ωo)=∫S2f(ωi,ωo)(cosθi)+dωi∫S2f(ωi,ωo)Li(ωi)(cosθi)+dωi=N1∑i=1Np(ωi,ωo)f(ωi,ωo)(cosθi)+N1∑i=1Np(ωi,ωo)f(ωi,ωo)Li(ωi)(cosθi)+=∑i=1N4⋅(NdotL)+⋅(NdotV)+D⋅F⋅G⋅(NdotL)+D⋅(NdotH)+4⋅(LdotH)+∑i=1N4⋅(NdotL)+⋅(NdotV)+D⋅F⋅G⋅L⋅(NdotL)+D⋅(NdotH)+4⋅(LdotH)+=∑i=1N(NdotV)+⋅(NdotH)+F⋅G⋅(LdotH)+∑i=1NL(NdotV)+⋅(NdotH)+F⋅G⋅(LdotH)+=∑i=1NWeight∑i=1NL⋅Weight where Weight=(NdotV)+⋅(NdotH)+F⋅G⋅(LdotH)+.
By "Equation (53)" of [Lagarde 2014], when Li(ωi) is constant, we have LD(ωo)=∑i=1NWeight∑i=1NL⋅Weight=∑i=1NWeightL⋅∑i=1NWeight=L which does
NOT depend on the Weight term. This means that when Li(ωi) is constant, no matter what Weight term we
use, we will always have an exact match.
By [Karis 2013], assuming that the Weight term equals (cosθi)+, we have the split integral approximationLD(ωo)=∑i=1NWeight∑i=1NL⋅Weight≈∑i=1N(NdotL)+∑i=1NL⋅(NdotL)+.
However, we do NOT have enough remaining dimensions to store both N (world space) and NdotV. By [Karis
2013], we assume "V = N" when pre-integrating the LD term, while the reflection direction R is used to
look up the LD term when rendering. By "Figure 10.37" of Real-Time Rendering Fourth Edition, by [Lagarde 2014], the
reflection direction is biased towards the normal direction to reduce the error generated by the non-clipped
specular lobe.
TODO: If we use the Paraboloid Map, perhaps we can use the 3D texture and provide one
dimension for the NdotV. Thus, we can use both N (world space) and NdotV to look up the LD term.
2-3-3. Variance Reduction
By "20.4 Mipmap Filtered Samples" of [Colbert 2007], we use the mip-level sample_mip_level=max(0,log(dωPdωS)) to sample the filtered L value from the
distant radiance distribution to reduce the variance. The dωS is the solid angle subtended by the sample, and we have dωS=N1⋅PDF1. The dωP is the solid angle subtended by a single texel of the zeroth mip-level of the texture, and by
"Projection from Cube Maps" of [Sloan 2008], we have dωP=cubesize_u⋅cubesize_v1⋅12+u2+v2⋅(12+u2+v2)4.