Linear Spherical Harmonics

In graphics, we often deal with spherical function, i.e. function that is defined on the surface of a sphere. Spherical Harmonics (SH) is a commonly used tool to represent spherical function by writing it as an infinite sum of these spherical harmonics. However, in games, we usually limit the number of basis functions to either 4 or 9 basis functions. They are called linear SH (SH2) and quadratic SH (SH3) respectively.

SH is an orthonormal basis over the sphere $S$ with the following parameterization:

$\quad s = \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} sin\theta \; cos\phi \\ sin\theta \; sin\phi \\ cos\theta \end{bmatrix}$

Least Square Formulation

Since there are only 4 coefficients for linear SH, it’s not that difficult to derive fully expanded equations which make things more readable. Let’s go back to basic and see how things work.

Without loss of generality, given a function $f(s)$ we want to approximate with basis functions $B_0(s)$, $B_1(s)$, $B_2(s)$ and $B_3(s)$. Let’s call this approximation as $g(s)$.

$\quad \begin{aligned} f(s) &\approx g(s) \\ g(s) &= b_0 B_0(s) + b_1 B_1(s) + b_2 B_2(s) + b_3 B_3(s) \end{aligned} $

In order to find the coefficients $b_0, b_1, b_2, b_3$ we simply solve least square equation:

$\quad \begin{aligned} E(s) &= \int (g(s) - f(s))^2 ds \\ &= \int ((b_0 B_0(s) + b_1 B_1(s) + b_2 B_2(s) + b_3 B_3(s)) - f(s))^2 ds \end{aligned} $

To find the minimum we just set the derivative of each coefficient to 0:

$\quad \frac{dE(s)}{db_0} = 0, \frac{dE(s)}{db_1} = 0, \frac{dE(s)}{db_2} = 0, \frac{dE(s)}{db_3} = 0$

Solving those equations and arranging it in a matrix, this will yield:

$\quad \begin{aligned} \begin{bmatrix} \int B_0(s) B_0(s) ds & \int B_0(s) B_1(s) ds & \int B_0(s) B_2(s) ds & \int B_0(s) B_3(s) ds \\ \int B_1(s) B_0(s) ds & \int B_1(s) B_1(s) ds & \int B_1(s) B_2(s) ds & \int B_1(s) B_3(s) ds \\ \int B_2(s) B_0(s) ds & \int B_2(s) B_1(s) ds & \int B_2(s) B_2(s) ds & \int B_2(s) B_3(s) ds \\ \int B_3(s) B_0(s) ds & \int B_3(s) B_1(s) ds & \int B_3(s) B_2(s) ds & \int B_3(s) B_3(s) ds \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{bmatrix} &= \begin{bmatrix} \int f(s) \; B_0(s) ds \\ \int f(s) \; B_1(s) ds \\ \int f(s) \; B_2(s) ds \\ \int f(s) \; B_3(s) ds \end{bmatrix} \\ A \; b &= c \end{aligned} $

We can simply choose our basis functions $B_0(s), B_1(s), B_2(s), B_3(s)$ and to find our coefficients, we can simply do

$\quad b = A^{-1} c = Proj(f)$

This operation of finding coefficients is often called projection.

Standard SH basis

Standard SH basis are orthonormal so the matrix $A$ is simply the identity matrix, the basis functions [1] and projection are given by:

The original function can then be approximated with:

This is all good however translating this directly in code yields many constants, for example $\frac{1}{2 \sqrt{\pi}}, \frac{\sqrt{3}}{2 \sqrt{\pi}}$. It turns out we simplify the equations by removing the constants in the basis functions.

Simplified SH basis

Following Geomerics formulation [2] using orthogonal (but not orthonormal) basis functions:

$\quad B_0(s) = 1, B_1(s) = x, B_2(s) = y, B_3(s) = z $

computing $b = A^{-1} \; c$ (note that A isn’t identity) yields the following coefficients:

$\quad \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{bmatrix} = \begin{bmatrix} \frac{1}{4\pi} \int f(s) \; B_0(s) ds \\ \frac{3}{4\pi} \int f(s) \; B_1(s) ds \\ \frac{3}{4\pi} \int f(s) \; B_2(s) ds \\ \frac{3}{4\pi} \int f(s) \; B_3(s) ds \end{bmatrix} $

Furthermore, we pull the constant factor $3$ from $b_1(s), b_2(s), b_3(s)$ to the reconstruction function. All of these yields the following:

Note that there’s linear relationship between the two coefficients which mean we can use standard SH math and then convert to the simplified version or vice versa.

Common SH operations

We will discuss common SH operations for both simplified and standard SH.

Projection of Constant function

Constant function $f(s) = 1$ defined in a spherical cap $\theta \in [0, \pi]$ and $\phi \in [0, 2\pi]$.

It is clear to see that the reconstructed function $g(s)$ equals exactly to $f(s)$.

Projection of Dirac Delta Function

Diract delta function $d$ is a function that is zero everywhere except at direction $d$ and integrates to $1$.

$\quad \begin{matrix} \delta(s) = \begin{cases} \infty &\text{if } s = d \\ 0 &\text{otherwise} \end{cases} & \; \; \; \displaystyle \int \delta(s) ds = 1 \end{matrix} $

Projection of dirac delta function to SH is simply the evaluation of the basis functions. Representing dirac delta function requires infinite sum so truncating the terms yield inaccurate approximation. This however is the energy conserving projection.

Projection of (Clamped) Cos Function

Since full integration of cos function evaluates to 0, we need to clamp the domain to be hemispherical. Given $f(\theta, \phi) = cos(\theta)$ with $\theta \in [0, \frac{\pi}{2}]$ and $\phi \in [0, 2\pi]$, the projection to SH is as following:

Oriented (clamped) cos function can be obtained by rotating/aligning clamped cos SH to direction $d$.

Convolution with Cos Function

This can be done directly in the frequency domain easily due to circular symmetry (look at Convolution section in [1]). Given a function $f_{SH}(s)$ convolved with Cos Function $h_{SH}(s)$, the resulting function $g_{SH}(s)$ coefficients are given by:

**^{$f_{g0}, f_{g1}, f_{g2}, f_{g3}$ are coefficients for $f_{SH}$ with simplified SH basis}
**^{$f_{b0}, f_{b1}, f_{b2}, f_{b3}$ are coefficients for $f_{SH}$ with standard SH basis}

Application to Rendering

Ok, now we have the knowledge of Spherical Harmonics, let’s now apply it to solve rendering equation:

$\displaystyle \quad L_o(p, \omega_o) = L_e(p,\omega_o) + \int_{S^2} f(p,\omega_o,\omega_i) \; L_i(p,\omega_i) \; (n \cdot \omega_i) \; d\omega_i$

In order to use the equation above for games, we use the following assumptions:

• $c_d$ for the value read from diffuse/albedo texture.
• Light intensity $I$ which can be defined as $I = \frac{L_i}{\pi}$.
• Energy conserving Lambertian BRDF $f(p,\omega_o,\omega_i) = \frac{c_d}{\pi}$.
• We don’t care about emission at the moment, so $L_e(p,\omega_o) = 0$.

The rendering equation can then be rewritten as:

$\quad L_o = c_d \int I_{\omega_i} \; (n \cdot \omega_i) \; d\omega_i = c_d \cdot Irr(n)$

**^{For directional light, the equation above becomes $L_o = c_d \; I \; (n \cdot L)$}
**^{For point light, we will need to add distance attenuation factor $\frac{1}{d^2}$, so it becomes $L_o = c_d \; \frac{I}{d^2} \; (n \cdot \ L)$}

Note that $(n \cdot \omega_i)$ is the clamped cosine function and the irradiance $Irr(n)$ can be computed by convolving $I$ with cos function and evaluate it at $n$.

$\quad \begin{aligned} Irr(n) &= \int I_{\omega_i} \; (n \cdot \omega_i) \; d\omega_i \\ &\approx g_0 G_0(n) + 3 \; ( g_1 G_1(n) + g_2 G_2(n) + g_3 G_3(n) ) \\ &= (\pi \; i_{g0}) \; G_0(n) + 3 \; ((\frac{2}{3}\pi \; i_{g1}) \; G_1(n) + (\frac{2}{3}\pi \; i_{g2}) \; G_2(n) + (\frac{2}{3}\pi \; i_{g3}) \; G_3(n)) \\ &= \pi \; (i_{g0} \; G_0(n) + 2 \; (i_{g1} \; G_1(n) + i_{g2} \; G_2(n) + i_{g3} \; G_3(n))) \\ &= \pi \; (i_{g0} + 2 \; (i_{g1} \; n_x + i_{g2} \; n_y + i_{g3} \; n_z)) \end{aligned} $

Projection of Ambient Light

Ambient light can be modeled as a constant in all direction and constant function can be reconstructed exactly in Spherical Harmonics.

Projection of Directional Light

Directional light can be represented with delta function. So projection of directional light can be done with projection of delta function.

Reconstruction Quality

Unfortunately, linear SH isn’t very accurate even for a simple cos function. There are two main problems:

1. Negative value (ringing)
2. Peak value doesn’t match

Handling Negative Value

Let’s look at the reconstruction again:

$\quad \begin{aligned} g(s) &= g_0 G_0(s) + 3 \; ( g_1 G_1(s) + g_2 G_2(s) + g_3 G_3(s) ) \\ &= g_0 + 3 \; ( g_1 x + g_2 y + g_3 z) \\ &= g_0 + 3 \; \begin{bmatrix} g_1 \\ g_2 \\ g_3 \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ z \end{bmatrix} \\ &= g_0 + 3 \; g \cdot d \end{aligned} $

Since we are projecting a non-negative signal $g_0$ will be non-negative. In order to guarantee non-negative result, we just need to scale the remaining term so that the minimum value is $0$.

Basically we want to find $c$ for the following:

$\quad \begin{aligned} min(g(s)) &= 0 \; \\ g_0 + c \; ( 3 \; min( g \cdot d ) ) &= 0 \\ g_0 + c \; ( 3 \; min( |g| \; |d| \; cos(\theta) ) ) &= 0 \\ g_0 + c \; ( -3 \; |g| ) &= 0 \\ c = \frac{g_0}{3 \; |g|} \end{aligned} $

Note that we will have 3 scaling factor for each RGB channel. In order to avoid color shifting we choose the minimum scaling factor $c = min(c_r, c_g, c_b)$.

Matching the Maximum Value for Directional Light

The projection of delta function requires infinite sum and because we truncate the terms, there will be approximation error. For directional light it might be desirable to match the peak so that irradiance is 1 when normal is exactly the same as light direction.

The simplest way to do this is to find a constant $c$ to scale the SH coefficients of the directional light. We first project the directional light as delta function with direction $(n_x, n_y, n_z)$ and solve the following:

$\quad \begin{aligned} c \cdot Irr(n) &= 1 \\ c \cdot \pi \; (i_{g0} + 2 \; (i_{g1} \; n_x + i_{g2} \; n_y + i_{g3} \; n_z)) &= 1 \\ c \cdot \pi \; (\frac{1}{4 \pi} + 2 \; (\frac{1}{4 \pi} n_x \; n_x + \frac{1}{4 \pi} n_y \; n_y + \frac{1}{4 \pi} n_z \; n_z)) &= 1 \\ c \cdot \pi \; \frac{3}{4 \pi} &= 1 \\ c &= \frac{4}{3} \end{aligned} $

We can then use this to find the coefficients for directional light (we can simply scale the coefficients with intensity $I$):

Non Linear Reconstruction

One benefit of pulling the constant 3 from the simplified SH projection is that we can reason about the illumination based on the ratio of the coefficients, i.e. $r = \frac{|R_1|}{R_0}$, whereby $R_0 = g_0$ and $|R_1| = \sqrt{{g_1}^2 + {g_2}^2 + {g_3}^2}$.

Graham Hazel ([3]) provides alternative way to reconstruct irradiance non-linearly by looking at multiple cases of $r$ and generalizing them. The original reference has excellent analysis and explanation, the following is just my note/understanding.

The irradiance reconstruction we want to look at:

$\quad \begin{aligned} Irr(n) &= \pi \; (g_0 + 2 \; (g_1 \; n_x + g_2 \; n_y + g_3 \; n_z)) \\ &= \pi \; (R_0 + 2 \; R_1 \cdot n) \\ &= \pi \; R_0 \; (1 + 2 \; \frac{R_1}{R_0} \cdot n) \end{aligned} $

Ambient light

For ambient light, it’s easy to see $r = 0$.

Directional light

For directional light (with direction $d$) projected to SH, the ratio will be:

$\quad \begin{aligned} r &= \frac{|R_1|}{R_0} \\ &= \frac{\sqrt{{g_1}^2 + {g_2}^2 + {g_3}^2}}{g_0} \\ &= \frac{\sqrt{(\frac{1}{4 \pi} \; d_x)^2 + (\frac{1}{4 \pi} \; d_y)^2 + (\frac{1}{4 \pi} \; d_z)^2}}{\frac{1}{4 \pi}} \\ &= 1 \end{aligned} $

Given above, the irradiance function can be simplified to $Irr(n) = \pi \; R_0 \; (1 + 2 \; \frac{R_1}{|R_1|} \cdot n)$.

Let’s try to fit $g = (1 + 2 \; \frac{R_1}{|R_1|} \cdot n)$ with quartic polynomial $f(x)$ where $x = \frac{R_1}{|R_1|}$ and $x \in [-1, 1]$.

$\quad \begin{aligned} f(x) &= a x^4 + b x^3 + c x^2 + d x + e \\ g(x) &= (1 + 2 x) \end{aligned} $

We want the following properties:

• Min value: $f(-1) = 0$.
  ^{We want non-negative reconstruction}
• Max value: $f(1) = \frac{4}{3} \; g(n) = 4$.
  ^{look at Matching the Maximum Value section}
• Preserve Energy: $\displaystyle \int_{-1}^{1} f(x) \; dx = \int_{-1}^{1} g(x) \; dx = 2$
• Derivative conditions: $f’(-1) = 0$ and $f”(-1) = 0$

If we put the 5 equations into a form of $Ax = b$, we will have the following:

That yields:

$\quad \begin{aligned} f(x) &= \frac{1}{2} x^3 + \frac{3}{2} x^2 + \frac{3}{2} x + \frac{1}{2} \\ &= \frac{1}{2} (x + 1)^3 \end{aligned} $

References

1. https://www.ppsloan.org/publications/StupidSH36.pdf
2. https://community.arm.com/cfs-file/__key/telligent-evolution-components-attachments/01-2066-00-00-00-01-27-70/Simplifying_2D00_Spherical_2D00_Harmonics_2D00_for_2D00_Lighting.pdf
3. https://grahamhazel.com/blog/2017/12/22/converting-sh-radiance-to-irradiance/
4. https://www.sjbrown.co.uk/posts/ambient-and-directional-lighting-in-spherical-harmonics/