by Leonard Ritter, Duangle GbR
This is a writeup on my derivation work on first order Spherical Harmonics, that is, integration of the SH basis with l=1, exclusively, which is useful for coarsely integrating and simulating highly diffuse, but energy preserving propagation of light.
It should be said first that everything we're doing here does not expand at all to higher order spherical harmonics, which is, at heart, a method to perform fourier transformations in the spherical domain.
It just so happens that under l=1, the Y function of SH has only one constant and three directional components which turn out to be entirely compatible with linear algebra, and so the math can be radically simplified here, in a way that has not been sufficiently documented yet, which this document shall rectify.
The full Y function to build argument dependent coefficients is defined as follows, in Scopes code:
fn Y (l m theta phi)
if (m > 0)
(sqrt 2.0:f64) * (K l m) * (cos (m as f64 * phi)) * (P l m (cos theta))
elseif (m < 0)
(sqrt 2.0:f64) * (K l m) * (sin (-m as f64 * phi)) * (P l -m (cos theta))
else # m == 0
(K l 0) * (P l 0 (cos theta))
As we source our angles from normal vectors as follows
theta := (acos n.z)
phi := (atan2 n.y n.x) # this normalizes n.xy by 1/(sin theta)
We can plug these equations into Y, and arrive at
fn Y (l m n)
if (m > 0)
(sqrt 2.0:f64) * (K l m) * (cos (m as f64 * (atan2 n.y n.x))) * (P l m n.z)
elseif (m < 0)
(sqrt 2.0:f64) * (K l m) * (sin (-m as f64 * (atan2 n.y n.x))) * (P l -m n.z)
else # m == 0
(K l 0) * (P l 0 n.z)
We only need the first four coefficients of the series, so let's plug those in, and also substitute for K:
fn Y_0_0 (n)
(sqrt (1 / (4 * pi))) * (P 0 0 n.z)
fn Y_1_-1 (n)
# (sqrt 2) * (sqrt (3 / (8 * pi))) = (sqrt (3 / (4 * pi)))
(sqrt (3 / (4 * pi))) * n.y / sin(acos(n.z)) * (P 1 1 n.z)
fn Y_1_0 (n)
(sqrt (3 / (4 * pi))) * (P 1 0 n.z)
fn Y_1_1 (n)
(sqrt (3 / (4 * pi))) * n.x / sin(acos(n.z)) * (P 1 1 n.z)
The definition of P is fairly complex, but for l=1, it hardly matters, as P(0,0,x)=1, P(1,0,x)=x and P(1,1,x)=-sqrt((1-x)*(1+x))=-abs(sin(acos(x))); so let's make those substitions as well:
fn Y_0_0 (n)
sqrt (1 / (4 * pi))
fn Y_1_-1 (n)
(sqrt (3 / (4 * pi))) * -n.y
fn Y_1_0 (n)
(sqrt (3 / (4 * pi))) * n.z
fn Y_1_1 (n)
(sqrt (3 / (4 * pi))) * -n.x
Wow, not much left of the mighty Y, eh? We can now build a basis for a normal vector from a simple vector multiplication. The signs can be dropped altogether, as they are an artifact of a mirrored coordinate system.
Y_0 := (sqrt (1 / (4 * pi))
Y_1 := (sqrt (3 / (4 * pi))
fn SH1 (n)
(vec4 (vec3 Y_1) Y_0) * (vec4 n 1)
As to rotating a vector in SH1 basis, we exploit the fact that the first four rows and columns of the SH rotation matrix are simply swizzled, mirrored versions of a regular rotation. We fixed all the swizzling and mirroring, which means we can now use our existing 3D matrix rotation functions. This only applies to l=1! SH vectors of higher order need special treatment.
Alright. When do we need this function? There are practically only three operations we perform on SH basis vectors frequently: construction (scaling of components, also called renormalization), integration (addition) and probing (dot product).
The construction part, you've just seen. Under addition, the coefficients remain constant, which means we can delay the renormalization until after the integration loop, which also helps with maintaining precision:
(SH1 v1) + (SH1 v2) + ... = (SH1 (v1 + v2 + ...))
Probing shows us why these coefficients exist in the first place; it reduces our multidimensional vector to a single irradiance value that can be used to scale a different basis. Here too, we can defer renormalization:
((SH1 v1) ⋅ (SH1 v2)) = (3 * (v1.x * v2.x + v1.y * v2.y + v1.z * v2.z) + v1.w * v2.w) / (4 * pi)
Moving on to analysis. What values do the individual components represent in practice? Our Monte Carlo integration kernel would, for example, now look like this:
let L_surface =
fold (L_surface = (vec4 0)) for s in (range N)
n := (uniformly-random-point-on-sphere-surface rng)
# compute light intensity at this point
L := (f n)
L_surface + (vec4 n 1) * L
# average and scale by surface area
L_surface := L_surface * (4 * pi / N)
# now we can renormalize it, if we like
SH_surface := (SH1 L_surface)
We can see here that the directional light components are averaged per coordinate axis, and the light intensity over the entire surface is stored in w. So when we analytically build integrals, we can factor out these renormalization coefficients entirely.
You might have noticed that 4 * pi, the area of the unit sphere (and also solid angle of its entire surface), comes up a lot, and you're probably wondering if you can normalize it out; The answer is yes, but code will be harder to understand, as your solid angles will have to be in the range 0..1 as well.
As a bonus, here are some popular analytical integrals:
- Ambient term:
SH1 (vec4 0 0 0 (4 * pi)) - A lambert shaded surface oriented towards
n:SH1 ((vec4 (n * 2 / 3) 1) * pi) - A half lambert shaded surface oriented towards
n:SH1 ((vec4 (n / 3) 1) * 2 * pi) - A hemispheric (i.e. directional) light source coming from
n:SH1 ((vec4 n 2) * pi) - A light source coming from
nwith solid angler:
k := (cos (r / 4))
`SH1 ((vec4 (n * (1 - k * k)) (2 * (1 - k))) * pi)`
Probing for irradiance, i.e. getting a normalized lighting term we can shade surfaces with, requires use of this function, which I lifted from the LPV paper:
fn irradiance_probe (v sharpness)
sh_c0 := (2.0 - sharpness) * 1.0
sh_c1 := sharpness * 2.0 / 3.0
vec4 (v.xyz * sh_c1) (v.w * sh_c0)
With a sharpness default value of 1, we get a half lambert shading. sharpness > 1 produces a ringing that destroys the ambient term but emphasizes the directional component. sharpness at sqrt(0.5) is my personal favorite. sharpness == 0 nulls the directional term entirely, giving everything a smooth, buttery, subsurface scattered look.
For probing a SH term with a surface normal, consider alternatively probing it with a reflection vector, as if you were sampling a cubemap. As the SH term is highly diffuse already, a specular probe can reconstruct some of the directional energy. At higher SH orders, the specular probe will look more and more like a mirror of the environment, and so here a diffuse probe is not a bad choice.
In my Light Propagation Volumes shadertoy, I've been going even further and tried to erase superfluous coefficients from intermediate computations as much as possible. This led to these forms:
vec4 sh_project(vec3 n) { return vec4(n, 0.57735026918963); }
This is all that's left of the renormalization function. There's only one coefficient left, which has been reduced to sqrt(1/3).
float sh_dot(vec4 a, vec4 b) {
return max(dot(a,b),0.0);
}
This is the probing function. The changes to sh_project have erased all coefficients required here. After this operation, to compute flux, the result can be scaled by the desired solid angle, which must be rescaled by 3 / (4 * pi) first.
Manufacturing a probe for shading goes as follows:
vec4 sh_irradiance_probe(vec4 v) {
return vec4(v.xyz * 2.0 / 3.0, v.w);
}
float shade_probe(vec4 sh, vec3 n) {
return sh_dot(sh_irradiance_probe(sh), sh_project(n)) * 3 / (4 * pi);
}