I describe how to project various types of lights into Spherical Harmonics (SH) basis analytically. In particular, I want to calculate the irradiance at a given point in space, which I'll call the shading point, due to contribution from a single light source.
I've implemented everything described in these notes for L2 SH as a single-file HLSL library: https://gist.github.com/pema99/226022f4e78173a2c21c56044e86d4a4
Given a function on the unit sphere,
Where
Directional lights are modeled nicely as Dirac delta distributions, in particular:
They have a property that lets us 'pluck' values out of an integral. This is known as the sifting property:
We can define the distribution of radiance for a directional light using the Dirac delta like so:
Where
This derivation gives the SH coefficients for radiance arriving at a point due to the light. To instead calculate the irradiance we can use the trick presented in An Efficient Representation for Irradiance Environment Maps by Ramamoorthi et al (more details in Spherical Harmonics > Radiance to irradiance). The authors show that this can be done by simple pointwise multiplication with a few constants that depend on the degree of SH. For L2 SH, these constants are:
So the final L2 SH coefficients for irradiance become:
Written out in HLSL shader code, we get something like this:
// Constants for radiance -> irradiance convolution
static const float A_HAT_0 = PI;
static const float A_HAT_1 = 2.0 * PI / 3.0;
static const float A_HAT_2 = PI / 4.0;
SphericalHarmonicsL2 ProjectDirectionalLightIrradianceSHL2(
float3 lightDirection,
float3 lightColor)
{
SphericalHarmonicsL2 sh;
float3 toLightNorm = normalize(-lightDirection);
sh.coefficients[0] = lightColor * Y0() * A_HAT_0;
sh.coefficients[1] = lightColor * Y1_1(toLightNorm) * A_HAT_1;
sh.coefficients[2] = lightColor * Y10(toLightNorm) * A_HAT_1;
sh.coefficients[3] = lightColor * Y11(toLightNorm) * A_HAT_1;
sh.coefficients[4] = lightColor * Y2_2(toLightNorm) * A_HAT_2;
sh.coefficients[5] = lightColor * Y2_1(toLightNorm) * A_HAT_2;
sh.coefficients[6] = lightColor * Y20(toLightNorm) * A_HAT_2;
sh.coefficients[7] = lightColor * Y21(toLightNorm) * A_HAT_2;
sh.coefficients[8] = lightColor * Y22(toLightNorm) * A_HAT_2;
return sh;
}
Where Y0
..Y22
are the SH basis functions for up to L2 SH.
The derivation for a point light is quite similar to that for directional light. The distribution of radiance now looks like so:
Where
Written out in HLSL shader code:
SphericalHarmonicsL2 ProjectPointLightIrradianceSHL2(
float3 shadingPosition,
float3 lightPosition,
float3 lightColor)
{
SphericalHarmonicsL2 sh;
float3 toLight = lightPosition - shadingPosition;
float distanceSquared = dot(toLight, toLight);
float3 toLightNorm = normalize(toLight);
sh.coefficients[0] = lightColor / distanceSquared * Y0() * A_HAT_0;
sh.coefficients[1] = lightColor / distanceSquared * Y1_1(toLightNorm) * A_HAT_1;
sh.coefficients[2] = lightColor / distanceSquared * Y10(toLightNorm) * A_HAT_1;
sh.coefficients[3] = lightColor / distanceSquared * Y11(toLightNorm) * A_HAT_1;
sh.coefficients[4] = lightColor / distanceSquared * Y2_2(toLightNorm) * A_HAT_2;
sh.coefficients[5] = lightColor / distanceSquared * Y2_1(toLightNorm) * A_HAT_2;
sh.coefficients[6] = lightColor / distanceSquared * Y20(toLightNorm) * A_HAT_2;
sh.coefficients[7] = lightColor / distanceSquared * Y21(toLightNorm) * A_HAT_2;
sh.coefficients[8] = lightColor / distanceSquared * Y22(toLightNorm) * A_HAT_2;
return sh;
}
Spot lights (aka cone lights) are again quite similar to point lights. The distribution of radiance for a spot light can be modeled as so:
Where
Where
There's no single 'correct' choice for the falloff function, but one simple choice could be a step function:
Where
Regardless of the choice of
Written out in HLSL shader code, assuming the smoothstep as the falloff function:
SphericalHarmonicsL2 ProjectSpotLightIrradianceSHL2(
float3 shadingPosition,
float3 lightPosition,
float3 lightDirection,
float3 lightColor,
float outerHalfAngle,
float innerHalfAngle)
{
SphericalHarmonicsL2 sh;
float3 toLight = lightPosition - shadingPosition;
float distanceSquared = dot(toLight, toLight);
float3 toLightNorm = normalize(toLight);
float cosTheta = dot(toLightNorm, -lightDirection);
float theta = acos(cosTheta);
float falloff = smoothstep(outerHalfAngle, innerHalfAngle, theta);
sh.coefficients[0] = lightColor / distanceSquared * falloff * Y0() * A_HAT_0;
sh.coefficients[1] = lightColor / distanceSquared * falloff * Y1_1(toLightNorm) * A_HAT_1;
sh.coefficients[2] = lightColor / distanceSquared * falloff * Y10(toLightNorm) * A_HAT_1;
sh.coefficients[3] = lightColor / distanceSquared * falloff * Y11(toLightNorm) * A_HAT_1;
sh.coefficients[4] = lightColor / distanceSquared * falloff * Y2_2(toLightNorm) * A_HAT_2;
sh.coefficients[5] = lightColor / distanceSquared * falloff * Y2_1(toLightNorm) * A_HAT_2;
sh.coefficients[6] = lightColor / distanceSquared * falloff * Y20(toLightNorm) * A_HAT_2;
sh.coefficients[7] = lightColor / distanceSquared * falloff * Y21(toLightNorm) * A_HAT_2;
sh.coefficients[8] = lightColor / distanceSquared * falloff * Y22(toLightNorm) * A_HAT_2;
return sh;
}
Spherical light sources are a bit more tricky, since they can't be modeled as delta distributions. There is however still an analytical method for these, described in the paper Harmonics Virtual Lights : fast projection of luminance field on spherical harmonics for efficient rendering by Mézières et al.
The technique uses a neat property of SH to simplify the projection - functions that are rotationally symmetric around the Z axis will project to only the zonal harmonics (basis functions
To project a spherical light into SH, we want to calculate this integral:
Where
As mentioned above, we will assume that this projection is being done in a coordinate space where the Z axis is
Note that I use
Where
The zonal harmonics can be defined like this:
Where
And
They can also be defined for negative degree:
We can use these definitions to rewrite the integral once more:
And integrate over $phi, which isn't actually used in the integrand:
Let's focus on the remaining definite integral for a moment. We can solve it by a change of variables, substituting
Multiplying both sides by
When
Which reduces the integrand to just the associated legendre polynomials with order 0. The antiderivative of these polynomials are:
So we can obtain the solution to the integral by simply applying the fundamental theorem of calculus:
The polynomials have the property that
So the original integral has the solution:
The full expression for the SH coefficients then becomes:
We can expand terms and simplify a bit:
Finally, we need to rotate these SH coefficients, which are assuming a coordinate centered on the light direction
We are essentially just multiplying with the SH basis functions evaluated in the +Z direction in the coordinate space we are rotating from. Note that this is only possible because
Translating all of this to HLSL shader code, we get:
// Associated legendre polynomial for zonal harmonics, up to L3
float AssociatedLegendreZH(uint l, float x)
{
[forcecase] switch (l)
{
case 0: return 1.0;
case 1: return x;
case 2: return 0.5 * (3.0 * x * x - 1.0);
case 3: return 0.5 * x * (5.0 * x * x - 3.0);
default: return 0.0;
}
}
SphericalHarmonicsL2 ProjectSphereLightIrradianceSHL2(
float3 shadingPosition,
float3 lightPosition,
float lightRadius,
float3 lightColor)
{
SphericalHarmonicsL2 sh;
float3 toLight = lightPosition - shadingPosition;
float distanceSquared = dot(toLight, toLight);
float3 toLightNorm = normalize(toLight);
// sin = opposite / hypotenuse = radius / distance
float sinAngle = lightRadius / max(sqrt(distanceSquared), 0.000001f);
// sin^2 + cos^2 = 1 => cos^2 = 1 - sin^2 => cos = sqrt(1 - sin^2)
float cosAngle = sqrt(saturate(1.0f - sinAngle * sinAngle));
// Calculate projection in a space where +Z is the direction towards the light
float l0Hat = sqrt(PI) * (1.0 - cosAngle);
float l1Hat = sqrt(PI / 3.0) * (AssociatedLegendreZH(0, cosAngle) - AssociatedLegendreZH(2, cosAngle));
float l2Hat = sqrt(PI / 5.0) * (AssociatedLegendreZH(1, cosAngle) - AssociatedLegendreZH(3, cosAngle));
// Rotate into global space, multiply on intensity, and convolve to irradiance to get final coefficients
sh.coefficients[0] = lightColor * sqrt(4.0 * PI) * Y0() * l0Hat * A_HAT_0;
sh.coefficients[1] = lightColor * sqrt(4.0 * PI / 3.0) * Y1_1(toLightNorm) * l1Hat * A_HAT_1;
sh.coefficients[2] = lightColor * sqrt(4.0 * PI / 3.0) * Y10(toLightNorm) * l1Hat * A_HAT_1;
sh.coefficients[3] = lightColor * sqrt(4.0 * PI / 3.0) * Y11(toLightNorm) * l1Hat * A_HAT_1;
sh.coefficients[4] = lightColor * sqrt(4.0 * PI / 5.0) * Y2_2(toLightNorm) * l2Hat * A_HAT_2;
sh.coefficients[5] = lightColor * sqrt(4.0 * PI / 5.0) * Y2_1(toLightNorm) * l2Hat * A_HAT_2;
sh.coefficients[6] = lightColor * sqrt(4.0 * PI / 5.0) * Y20(toLightNorm) * l2Hat * A_HAT_2;
sh.coefficients[7] = lightColor * sqrt(4.0 * PI / 5.0) * Y21(toLightNorm) * l2Hat * A_HAT_2;
sh.coefficients[8] = lightColor * sqrt(4.0 * PI / 5.0) * Y22(toLightNorm) * l2Hat * A_HAT_2;
return sh;
}
When researching analytical projection of various light types into SH, I was reminded of the paper Real-Time Polygonal-Light Shading with Linearly Transformed Cosines by Heitz et al, which describes a method for analytically evaluating direct illumination from polygonal light sources. The method is fairly straightforward, building upon a closed-form expression for irradiance of a polygonal light source that has been known since the 80s, and extending it with handling for non-Lambertian BRDFs.
I wondered whether we can build upon this work to enable projecting irradiance from polygonal light sources into SH analytically. I found one 2018 paper of particular interest: Analytic Spherical Harmonic Coefficients for Polygonal Area Lights by Wang et al. It presents a fully analytical method for doing just this. The mathematical derivation of their method in all its gory detail goes over my head, so I'll settle for describing it in broad strokes. Hope I don't get this too wrong 😅
The paper seeks to provide an analytical solution to the integral:
Just as in the section on Spherical lights,
The paper breaks down the problem by building upon a result from another paper, Sparse Zonal Harmonic Factorization for Efficient SH Rotation by Nowrouzezahrai et al, which showed that every SH basis function can be written as a linear combination of rotated versions of just the zonal harmonics:
Where
If we can evaluate this integral for each rotated zonal harmonic in the pre-computed set, we can evaluate the full integral for
With
The paper refers to
Note that explicit dependence of
Where
And
The only symbols left unaccounted for are
So (
Phew, that's a lot of recurrence relations. These relations afford an iterative algorithm for computing the original integral for
The code for this is quite ugly, involving a bunch of hardcoded pre-computed values. It's also rather long, so I won't go over it here. It's a fairly direct port and simplification of the implementation provided as supplemental materials for the paper. If you want to read my implementation, see here.
The HLSL library in which I've implemented these techniques is primarily intended for use with Unity, but doesn't actually depend on it in any way. I've strived for correctness in the implementation, which leads to a few pitfalls when trying to use the library in Unity, due to some quirks and inconsistencies in the engine.
Ideally, the result of analytical and numerical integration (i.e. baked light probes) should match exactly for each light type. A few issues make this tricky in practice:
SHL2ToUnityUniforms
which takes care of conversion from 'proper' SH coefficients to Unity's uniforms.Light.intensity
should be passed through GammaToLinearSpaceExact
before being fed to the functions in the library to match baked probes. This is assuming you are using linear color space (please don't use gamma color space).I managed to make every function in the library match Unity's baked output in the builtin render pipeline, but I had to apply a bunch of fixes due to the aforementioned. They do however all match reference values generated by simple monte carlo integration.