We typically think of elements of vector spaces, vectors, as lists of numbers. However, this is limiting a definition. In fact, if we loosen this definition, we can use functions as our vectors, forming a function space. The addition and scaling operators are defined as such:

Where is a scalar, and and are functions.

The elements of a basis for a function space is then called a basis function. A simple example of a function space is a the space of single-indeterminate polynomials. One possible basis for this vector space is the monomial basis:

This is clearly a valid basis as every polynomial can be written as a linear combination of the elements:

Where are constant coefficients for each basis function.

As mentioned just above, 2 operations of particular importance are projection and reconstruction. They are both fairly straight forward to define. If we want to project a function into the CH basis, we do it like so:

Where is the i'th basis function, and is the coefficient to associate with the i'th basis function. Naturally, to reconstruct our approximate version of , we take a simple inner product:

A pragmatic way to think about Circular Harmonics are as a kind of compression - a way to store a -periodic function numerically using very little space. The more coefficients we store, the better the approximation will be, but we'll quickly face diminishing returns. In other words, the approximation will already be quite good with only a few coefficients for many kinds of functions, and in particular those of low-frequency.

Next I'll show a concrete example of the aforementioned projection and reconstruction operations using some Python code.

First, let's define a function which we will project into CH basis and then reconstruct. I've chosen a saw wave:

def f(x):
    x = np.mod(x, 2*np.pi)
    return x/(2*np.pi)
Copy

Next, let's define a function that gives us the n'th basis function:

import numpy as np

# n is the index of the basis function, x is the input to the function
def fourier_basis(n, x):
    multiplier = np.ceil(n / 2)
    if n == 0:
        return 1 / np.sqrt(2 * np.pi)
    elif n % 2 == 1:
        return np.cos(multiplier * x) / np.sqrt(np.pi)
    else:
        return np.sin(multiplier * x) / np.sqrt(np.pi)
Copy

Next, we need a way to calculate integrals. I've implemented a simple Riemann sum:

# function is the function to integrate
# a and b are the bounds of the integral
# n is the number of discrete steps to take
def riemann(function,a,b,n):
    sumval = 0
    h = (b-a)/n
    for i in range(0,n-1):
        current_x = a+i*h
        sumval = sumval + function(current_x) * h
    return sumval
Copy

With those 2 building blocks, we can define a function that calculates coefficient for the n'th basis function:

# f is the function to project into the CH basis
# n is the index of the basis function
# calculates the integral of the n'th basis function multiplied by f
def fourier_project(n, f):
    return riemann(lambda x: f(x) * fourier_basis(n, x), 0, 2*np.pi, 1000)
Copy

We can then define a function that performs reconstruction. We'll use 30 basis functions for this example:

# f is the function to project into CH basis, then reconstruct
# x is the is the input to the function
def fourier_reconstruct(f, x):
    sum = 0
    for n in range(30):
        sum += fourier_project(n, f) * fourier_basis(n, x)
    return sum
Copy

Finally, we can compare the actual and reconstructed functions:

from matplotlib import pyplot as plt

x = np.linspace(-np.pi*2, np.pi*2, 1000)
fig, ax = plt.subplots()
ax.plot(x, f(x), label='f(x)')
ax.plot(x, fourier_reconstruct(f, x), label='f_reconstructed(x)')
ax.legend()
plt.show()
Copy

Looks pretty good to me!

Projection and reconstruction with the SH basis is completely analogous to case of circular harmonics. For projection:

Where is the coefficient to associate with basis function . If you are confused about this integral, please see Spherical integrals.

Reconstruction is again just an inner product:

Where is the order (amount of levels) of the basis we are using.

Note: The reason why I use Spherical Coordinates in the above description, despite having just explained that we usually prefer to use unit vectors, is primarily to make the integral more tenable. The space of unit vectors doesn't map intuitively to a pleasant domain of integration. Alternatively, we could write the double integral as a single integral in Solid Angle domain, but integrals in solid angle domain are not easy to calculate directly, so we usually translate to spherical coordinates before calculating anyways.

Two quantities of particular importance when rendering scenes with global illumination (GI) are radiance and irradiance. I'll briefly summarize what these quantities mean below:

Radiance is the quantity typically associated with a single ray of light, while irradiance is the quantity associated with a point on a surface. You can think of radiance as roughly being being the energy carried by a ray of light, and irradiance as the total energy incident on a point on a surface from all directions.

We can calculate the irradiance at a given point on a surface by integrating radiance over the hemisphere of directions centered around the normal vector at the point, attenuated by the cosine of the angle between incoming light direction and normal vector:

Here denotes the hemisphere centered around normal vector , is the point on the surface. denotes incoming light direction. The term is the dot of the normal and incoming light direction, ie. the cosine of the angle between the 2. The entire is sometimes aptly called the "clamped cosine term" and is sometimes written as .

Note: The cosine term is necessary due to how projected area varies with the angle between the surface normal and the normal of the area being projected. Roughly speaking, we want to weigh light coming from shallow angles less, as that light is spread over a larger area. For some context, see the images in Cosine weighted sampling > Why.

One technique for lighting scenes with GI is light probes. These are positions in space with which we associate a quantity of light. When it is time to shade a surface point, we can grab the few nearest probes to said point, interpolate between them, and use them to get an estimate of the indirect light contribution. That begs the question - what exactly do we store in each light probe? Spherical Harmonics turn out to be a great a fit, for a few main reasons:

• They are trivial to interpolate between.
• They can approximate low-frequency signals (such as is the case for irradiance) using very little data. We'll typically L2 Spherical Harmonics, which amounts to 27 coefficients per probe assuming we are using RGB (as opposed to a different method of sampling the visible spectrum).
• They are very cheap to evaluate in realtime. Evaluation boils down to basically a few dot products.

Recall that the radiometric quantity associated with light incident on a surface point is irradiance. It should then be no surprise that the quantity we should project into the SH basis and store in each light probe is also irradiance. To compute the coefficients for each probe, we can use path tracing.

Path tracing involves shooting a bunch of light rays from a point, bouncing them around the scene, and attenuating their contribution whenever they hit a surface. Note however that the radiometric quantity associated with a ray of light is radiance, not irradiance. I'll ignore this mismatch in desired quantities for now, and first explain how to project radiance into the L2 SH basis. Below is some pseudo-C#-code to do this. Note that I'm not going to bother with color - we'll just be using monochromatic lighting for simplicity.

// Given a probe position, returns the SH coefficients for radiance
// incoming from all directions to that position, using monte carlo
// integration.
float[] ProjectRadianceIntoSH(Vector3 probePosition)
{
	float[] radianceSH = new float[9];
    for (int i = 0; i < TotalSamplesPerProbe; i++)
    {
	    // Trace a light path going in a random direction to get radiance
	    // associated with the direction
	    Vector3 rayDirection = GetRandomDirection();
	    float radiance = TracePath(probePosition, rayDirection);

		for (int n = 0; n < 9; n++)
		{
			// SHBasis(n, d) evaluates the n'th SH basis function
			// given a direction vector d.
			radianceSH[n] += radiance * SHBasis(n, rayDirection);
		}
    }
    
    // Monte carlo normalization
    float reciprocalSampleCount = 1.0f / TotalSamplesPerProbe;
    float reciprocalUniformSphereDensity = 4.0f * Math.PI;
    return radianceSH * reciprocalSampleCount * reciprocalUniformSphereDensity;
}
Copy

For convenience, I've used a flat index for indexing SH basis functions, rather than an explicit level and index-in-level . Index 0 thus means basis function , index 1 is , index 2 is etc. Using L2 SH means we have to store 9 coefficients for the 9 basis functions. This pseudocode is essentially just implementing the SH projection operator described earlier, with radiance in a given direction, as the function to project.

Now we are faced the issue I outlined earlier: If we want to actually shade a surface with this, we need irradiance rather than radiance, which means that we essentially need to calculate the integral described in A brief primer on radiometry:

Ramamoorthi and Hanrahan describe how to do just this in their paper On the relationship between radiance and irradiance: determining the illumination from images of a convex Lambertian object. I'm not going to go over all the details of their derivations, but focus on the conclusions in broad strokes - you can read the paper for more info.

Ramamoorthi defines an operator for the clamped cosine term called :
Where:

• is the angle between the surface normal and the incident light direction.
• is the coefficient associated with the l'th SH level, of the clamped cosine term projected into the SH basis.
• is the SH basis function at level , index 0.

It's important to note here that we are only indexing the SH coefficients of projected into the SH basis using the level , and omitting the index in said level. The clamped cosine term is only dependent on the polar angle , and not the azimuthal angle . This should make intuitive sense - the total amount of light reflected at a point on a lambertian surface is only dependent on how shallow the angle of incidence is. Due to this property, we can discard any coefficient that isn't associated with a zonal harmonic (remember, these are basis functions with index 0) - and we end up with only as many coefficients for the projection of as we have levels of SH.

Ramamoorthi then goes on to derive an analytical expression for the coefficients :

Let's now revisit the irradiance integral with these definitions. We'll start by rewriting it a tiny bit:

Note that we have dropped the surface position from irradiance function . This is because we assume the illumination is distant, which means changes in position have negligible effect.

Given the nice property I described in the section Spherical Harmonics > Integral of product, one may now be tempted calculate irradiance integral like so:

However, this would be incorrect. The reason is a mismatch in coordinate systems. Incoming radiance, is defined with respect to global coordinates, while the clamped cosine term is with a local coordinates ie. the coordinate system where the surface normal at the point we are shading is pointing straight upwards. We need to introduce a rotation term in order to account for this. I'll omit some details here, as this term is quite complicated to calculate, and just name it somewhat opaquely . With this term, we can now write irradiance as:

The effect that this term has when multiplied onto the summand is to transform the input to the radiance function from local coordinate space, which is the domain of our irradiance integral, into global coordinate space, which is the domain of the radiance function, i.e.:

Ramamoorthi shows that the rotation term can be rewritten to:

Thus the expression for irradiance becomes:

In the follow up paper An Efficient Representation for Irradiance Environment Maps, Ramamoorthi further simplifies this expression by defining a new set of coefficients :

Which lets us finally express irradiance as:

Next, we rewrite irradiance in terms of SH coefficients, using what we know from the sections on projection:

It then becomes evident from the last 2 equations that:

Phrased in plain English, if we have radiance in SH, given by coefficients , and we want irradiance in SH, given by coefficients , we need only multiply the radiance coefficients by a factor which varies only per SH level.

Note: You may notice that the above expression looks quite similar to the expression in the Spherical Harmonics > Convolution section. This is because the transformation from radiance to irradiance just that - a convolution of radiance with the clamped cosine term.

Since we already have an analytical expression for , and differs from only by a constant factor per SH level, we can derive an analytical expression for too, which is exactly what Ramamoorthi does next:

The term for the first 3 levels of SH (which is what we typically care about for light probes), is:

We can put this all together to write some pseudocode that convolves our SH coefficients for radiance incoming from a given direction, to irradiance given the normal vector of a surface point:

// Given 9 SH coefficients (L2 SH) encoding directional radiance,
// return 9 SH coefficients encoding directional irradiance
float[] ConvolveRadianceToIrradianceSH(float[] radianceSH)
{
	float[] irradianceSH = new float[9];
	
	// L0
	float aHat0 = 3.14159;
	irradianceSH[0] = aHat0 * radianceSH[0];
	
	// L1
	float aHat1 = 2.09439;
	irradianceSH[1] = aHat1 * radianceSH[1];
	irradianceSH[2] = aHat1 * radianceSH[2];
	irradianceSH[3] = aHat1 * radianceSH[3];
	
	// L2
	float aHat2 = 0.78539;
	irradianceSH[4] = aHat2 * radianceSH[4];
	irradianceSH[5] = aHat2 * radianceSH[5];
	irradianceSH[6] = aHat2 * radianceSH[6];
	irradianceSH[7] = aHat2 * radianceSH[7];
	irradianceSH[8] = aHat2 * radianceSH[8];

    return irradianceSH;
}
Copy

As well as a function to shade a surface given a surface normal, surface albedo and a set of SH coefficients corresponding to the a light probe:

float ShadeSH(float[] irradianceSH, float surfaceAlbedo, Vector3 surfaceNormal)
{
    float outputColor = 0;
    for (int n = 0; n < 9; n++)
    {
        outputColor += irradianceSH[n] * SHBasis(n, surfaceNormal);
    }
    return outputColor * surfaceAlbedo;
}
Copy

Note that I am still using monochromatic lighting for simplicity, so the surface albedo is just a single scalar.

With this, we can finally write pseudocode to shade a point using a light probe:

// First, we "bake" a probe at a given position:
Vector3 probePosition = ...;
float[] radianceSH = ProjectRadianceIntoSH(probePosition);
float[] irradianceSH = ConvolveRadianceToIrradiance(radianceSH);

// Next, we use it to shade a point with a given normal:
Vector3 surfaceNormal = ...;
float surfaceAlbedo = ...;
float colorAtSurface = ShadeSH(irradianceSH, surfaceAlbedo, surfaceNormal);
Copy

In reality, instead of using a single probe, we would bake multiple probes, yielding several sets of SH coefficients, and then interpolate the few nearest ones to get a set of interpolated SH coefficients for irradiance. This is just a dumbed down example.