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

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

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.

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:

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

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:

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:

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:

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.