Analytically projecting lights into SH

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

Projection into SH
...

Given a function on the unit sphere, , we can calculate the coefficients for the projection of this function into SH basis like so:

Where is the unit sphere, is the SH basis function at level , index , and are the SH coefficients for each basis function. This integral is in Solid angle domain.

Directional light
...

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 is the emitted radiance of the light (an authoring input), is the direction of directional light (also an authoring input), and is an input direction. This distribution nicely describes how all the contribution from a directional light comes from a single direction. If we plug it into the formula for SH projection, we get:

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.

Point light
...

The derivation for a point light is quite similar to that for directional light. The distribution of radiance now looks like so:

Where is the distance from the shading point to the light, and is the direction pointing from the shading point to the light. Thus, the SH coefficients for irradiance due to the point light are:

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 light
...

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 is a falloff function describing the cone of the spot light. It's a bit unintuitive to think of it as a function , so we can reparametrize it:

Where is the direction that the spot light is pointing in. The term is the angle between the spot light's direction, and the direction from the light to the shading point. If it exceeds the half-angle defining the light cone, the shading point is outside of the cone, and should retrieve no contribution from the spot light.
Pasted image 20250531193854.png
There's no single 'correct' choice for the falloff function, but one simple choice could be a step function:

Where is the half-angle of the light cone. This results in a very sharp falloff at the edge of the cone. We could make it smoother by defining 2 half-angles - an outer and inner half-angle - and returning a value between 0 and 1 in the region between the angles. We could for example use smoothstep with the outer and inner half-angles being the lower and upper interpolation limits respectively.

Regardless of the choice of , the derivation of the SH coefficients follows the same pattern as the with point and directional lights. We get:

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
...

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 , the center column of the SH pyramid). The projection of a spherical light onto the unit sphere is a spherical cap. If we align the coordinate system such that the Z axis is the direction from the shading point to the center of the light (I'll call this direction ), the spherical cap will be centered on the Z axis, and thus project to only the zonal harmonics. We can then rotate the resulting coefficients back into global space to obtain the final SH coefficients.

To project a spherical light into SH, we want to calculate this integral:

Where is again emitted radiance, and is an indicator function that is 1 for values of that lie within the projection of the light onto the unit sphere centered on the shading point, and 0 elsewhere. Since only directions for which will contribute anything, we can rewrite the integral a bit. Let be the projection of the light onto the unit sphere, we get:

As mentioned above, we will assume that this projection is being done in a coordinate space where the Z axis is , so we can disregard cases where :

Note that I use to indicate that we are talking about SH coefficients given the specific coordinate space where is the Z axis. We can rewrite the integral from solid angle domain to a double integral over angles:

Where is the half-angle of the projected spherical cap as shown in the figure below (see Spherical integrals if you are confused about the ).
Pasted image 20250530235650.png

The zonal harmonics can be defined like this:

Where

And are the associated legendre polynomials with order 0, degree . The first 4 of these are:

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 . We first differentiate both sides with respect to :

Multiplying both sides by , we get:

When , and when , so integral becomes:

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 for all , so we can simplify the integral further:

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 , into global coordinate space. In other words, we want to apply a transformation that takes to . For this, we use a technique from Local, Deformable Precomputed Radiance Transfer (section 3.1), which shows that the rotation via convolution (Spherical Harmonics > Convolution), like so:

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 is representing a function with rotational symmetry around the Z-axis. Rotation of spherical harmonics in the general case is more involved.

are the coefficients for radiance arriving at a point due to the spherical light. As before, we need to use the trick from An Efficient Representation for Irradiance Environment Maps to get coefficients for irradiance:

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;  
}

Polygonal light
...

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, are the coefficients describing incoming radiance at a shading point due to contribution from a (polygonal) light source. is the projection of the polygon onto the unit sphere centered on the shading point, and are the SH basis functions.

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 are the axes of rotation (+Z) for each zonal harmonic, and are weights for each zonal harmonic. is bounded by , but can in practice be lower. These and are pre-computed with an offline optimizer. This result lets the authors consider only a simpler integral for radiance, giving just zonal harmonic coefficients for some specific direction :

If we can evaluate this integral for each rotated zonal harmonic in the pre-computed set, we can evaluate the full integral for by combining the results of each sub-integral. As mentioned in the previous section, the zonal harmonics are defined by:

With being the associated legendre polynomials. This lets us simplify the integral a bit:

The paper refers to as the 'surface integral'. The primary contribution of the paper is deriving a recurrence relation for :

Note that explicit dependence of on () is omitted to simplify notation. The only unaccounted-for term here is . This is a separate integral referred to as the 'boundary integral'. Its domain are the contours (edge) of the polygons projection:
Pasted image 20250601200711.png
can be written as a weighted sum of line integrals, one for each edge. These are denoted with being the index of the edge. The paper gives another recurrence relation for these integrals:

Where is yet another line integral, for which the paper gives yet another recurrence relation:

And is another (!) line integral, with another recurrence relation:

The only symbols left unaccounted for are and . is the angle subtended by each edge of the projected polygon, and are the coordinate of (the axis of rotation of the specific zonal harmonic we are considering) in a coordinate space where the +Z axis is aligned with the direction pointing from the shading point to the first vertex of the edge . In more precise terms: Let be an edge of the projection of the polygon onto the unit sphere, with vertices and , then:

So () forms a coordinate space centered on the vertex , wherein () are the coordinates of .

Phew, that's a lot of recurrence relations. These relations afford an iterative algorithm for computing the original integral for which goes something like this:

  • For each pre-computed zonal harmonic with axis of rotation :
    • Compute a few base cases for the various integrals. The paper computes , , , , , and using relatively simple closed-form solutions that rely on these integrals having simple integrands.
    • Repeatedly apply the 3 recurrence relations to compute for higher and higher values of .
  • For each SH basis function , compute the corresponding coefficient calculating the linear combination of the corresponding pre-computed zonal harmonics with axes of rotation and weights using the result from Nowrouzezahrai et al.

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.

Unity specifics
...

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:

  • Contribution from area lights and emissive meshes is off by . That is, light probes baked in a scene lit only with area light and emissive meshes will contain rather than just irradiance. This is not the case for point, spot, and directional lights. To make the functions for spherical and polygon light match Unity's baked output, one must divide the SH coefficients by .
  • Unity uses a somewhat intuitive format for the shader uniforms storing interpolated light probe data (unity_SHAr..unity_SHC). In particular, the constant part of each SH basis function is baked into the coefficient, the L1 terms are swizzled, and the L0 term contains part of the L2 term due to an optimization. I've provided a helper function SHL2ToUnityUniforms which takes care of conversion from 'proper' SH coefficients to Unity's uniforms.
  • The light color / intensity is specified in the light inspector is in gamma-space, but the probe baker uses linear space. So the intensity obtained from 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).
  • Some scriptable render pipelines (URP) premultiply light intensities by while others (HDRP) don't. This is another source of-by-pi to keep in mind.
  • Unity's light probes are interpolated on CPU side, and one set of coefficients is used for the entire mesh. If you do the analytical SH projection per-pixel, the result will look different (but nicer).

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.

Further reading
...

Referenced papers:
...