Spherical Harmonic Basis Functions
This page has been updated from its 2004 original.
Code for the function generator and rotations can now be found at https://github.com/sjb3d/sh.
Spherical Harmonic Function Generator
The real-valued spherical harmonics can be defined as:
$$\begin{align*}
Y_l^m(\phi, \theta) &= \Theta_l^{|m|}(\theta) \, \Phi_m(\phi) \\
\Theta_l^m(\theta) &= \sqrt{\frac{2l + 1}{4\pi}\frac{(l - m)!}{(l + m)!}}P^m_l(\cos\theta) \\
\Phi_m(\phi) &= \begin{cases}
\sqrt{2}\cos{m\theta}, &m>0 \\
1, &m=0 \\
\sqrt{2}\sin{|m|\theta}, &m<0
\end{cases}
\end{align*}$$
The associated Legendre polynomials $P^m_l$used above can be defined recursively as:
$$\begin{align*}
P^0_0(\cos\theta) &= 1 \\
P^m_m(\cos\theta) &= (2m - 1)\sin\theta \, P^{m-1}_{m-1}(\cos\theta) \\
P^m_{m+1}(\cos\theta) &= (2m + 1)\cos\theta \, P^m_m(\cos\theta) \\
P^m_l(\cos\theta) &= \frac{2l - 1}{l - m}\cos\theta \, P^m_{l-1}(\cos\theta) + \frac{1 - l - m}{l - m}P^m_{l-2}(\cos\theta)
\end{align*}$$
Out of curiosity I wrote a program to multiply everything out in terms of sines and cosines of phi and theta.
Here are the results for the first 3 orders:
$$\begin{align*}
Y_0^{0} &= \frac{1}{2}\sqrt{\frac{1}{\pi}}\\
Y_1^{-1} &= \frac{1}{2}\sqrt{\frac{3}{\pi}}\sin\phi\,\sin\theta\\
Y_1^{0} &= \frac{1}{2}\sqrt{\frac{3}{\pi}}\cos\theta\\
Y_1^{1} &= \frac{1}{2}\sqrt{\frac{3}{\pi}}\cos\phi\,\sin\theta\\
Y_2^{-2} &= \frac{1}{2}\sqrt{\frac{15}{\pi}}\sin\phi\,\cos\phi\,\sin^2\theta\\
Y_2^{-1} &= \frac{1}{2}\sqrt{\frac{15}{\pi}}\sin\phi\,\sin\theta\,\cos\theta\\
Y_2^{0} &= \frac{1}{4}\sqrt{\frac{5}{\pi}}(3\cos^2\theta-1)\\
Y_2^{1} &= \frac{1}{2}\sqrt{\frac{15}{\pi}}\cos\phi\,\sin\theta\,\cos\theta\\
Y_2^{2} &= \frac{1}{4}\sqrt{\frac{15}{\pi}}(\cos^2\phi-\sin^2\phi)\sin^2\theta\\
\end{align*}$$
And the next 2 orders:
$$\begin{align*}
Y_3^{-3} &= \frac{1}{4}\sqrt{\frac{35}{2\pi}}(3\sin\phi\,\cos^2\phi-\sin^3\phi)\sin^3\theta\\
Y_3^{-2} &= \frac{1}{2}\sqrt{\frac{105}{\pi}}\sin\phi\,\cos\phi\,\sin^2\theta\,\cos\theta\\
Y_3^{-1} &= \frac{1}{4}\sqrt{\frac{21}{2\pi}}\sin\phi(5\sin\theta\,\cos^2\theta-\sin\theta)\\
Y_3^{0} &= \frac{1}{4}\sqrt{\frac{7}{\pi}}(5\cos^3\theta-3\cos\theta)\\
Y_3^{1} &= \frac{1}{4}\sqrt{\frac{21}{2\pi}}\cos\phi(5\sin\theta\,\cos^2\theta-\sin\theta)\\
Y_3^{2} &= \frac{1}{4}\sqrt{\frac{105}{\pi}}(\cos^2\phi-\sin^2\phi)\sin^2\theta\,\cos\theta\\
Y_3^{3} &= \frac{1}{4}\sqrt{\frac{35}{2\pi}}(\cos^3\phi-3\sin^2\phi\,\cos\phi)\sin^3\theta\\
Y_4^{-4} &= \frac{3}{4}\sqrt{\frac{35}{\pi}}(\sin\phi\,\cos^3\phi-\sin^3\phi\,\cos\phi)\sin^4\theta\\
Y_4^{-3} &= \frac{3}{4}\sqrt{\frac{35}{2\pi}}(3\sin\phi\,\cos^2\phi-\sin^3\phi)\sin^3\theta\,\cos\theta\\
Y_4^{-2} &= \frac{3}{4}\sqrt{\frac{5}{\pi}}\sin\phi\,\cos\phi(7\sin^2\theta\,\cos^2\theta-\sin^2\theta)\\
Y_4^{-1} &= \frac{3}{4}\sqrt{\frac{5}{2\pi}}\sin\phi(7\sin\theta\,\cos^3\theta-3\sin\theta\,\cos\theta)\\
Y_4^{0} &= \frac{3}{16}\sqrt{\frac{1}{\pi}}(35\cos^4\theta-30\cos^2\theta+3)\\
Y_4^{1} &= \frac{3}{4}\sqrt{\frac{5}{2\pi}}\cos\phi(7\sin\theta\,\cos^3\theta-3\sin\theta\,\cos\theta)\\
Y_4^{2} &= \frac{3}{8}\sqrt{\frac{5}{\pi}}(\cos^2\phi-\sin^2\phi)(7\sin^2\theta\,\cos^2\theta-\sin^2\theta)\\
Y_4^{3} &= \frac{3}{4}\sqrt{\frac{35}{2\pi}}(\cos^3\phi-3\sin^2\phi\,\cos\phi)\sin^3\theta\,\cos\theta\\
Y_4^{4} &= \frac{3}{16}\sqrt{\frac{35}{\pi}}(\cos^4\phi-6\sin^2\phi\,\cos^2\phi+\sin^4\phi)\sin^4\theta\\
\end{align*}$$
The code on github can generate higher orders if required.
Convolutions
An advantage of working in spherical harmonics is that convolutions can be implemented very efficiently.
For example, suppose you have incoming radiance in spherical harmonics $L^m_l$, and you wish to compute the irradiance of a diffuse patch with some normal $E^m_l$ (or put more simply, the response of a diffuse surface parameterised by normal).
This can be computed by convolving the incoming radiance with a clamped cosine kernel.
Summarising "An Efficient Representation for Irradiance Environment Maps" and "On the Relationship between Radiance and Irradiance":
$$ E^m_l = A_l \, L^m_l$$
Where $A_l$ is defined as:
$$A_n = \begin{cases}
2\pi\frac{(-1)^{\frac{n}{2} - 1}}{(n+2)(n-1)}(\frac{n!}{2^n(\frac{n}{2}!)^2}), &n\text{ even} \\
\frac{2\pi}{3}, &n=1 \\
0, &n>1,\text{ odd}
\end{cases}$$
The first few terms of this multiply out to:
$$\begin{align*}
A_0 &= \pi \\
A_1 &= \frac{2\pi}{3} \\
A_2 &= \frac{\pi}{4} \\
A_3 &= 0 \\
A_4 &= -\frac{\pi}{24}
\end{align*}$$
Since the terms are small beyond $A_2$, representing irradiance with the first 3 orders (i.e. 9 spherical harmonic basis functions) has low error.
Another neat result is for the Henyey-Greenstein phase function. The phase function for some angle $\theta$ between viewer and light has a single parameter $g$ as follows:
$$\rho(\theta) = \frac{1 - g^2}{4\pi(1 + g^2 - 2g\cos\theta)^{1.5}}$$
The spherical harmonics for outgoing radiance (parameterised by viewer direction) can be computed using $A_l$ defined as:
$$A_n = g^n$$
Which is awesome.
Spherical Harmonic Rotations
Rotations of high order spherical harmonics can be implemented using the method in "Rapid and Stable Determination of Rotation Matrices between Spherical Harmonics by Direct Recursion".
To get this to work I had to replace the equations 8.8 through 8.14 with the following:
$$\begin{align*}
R_l^{m,n} &= \begin{cases}
\frac{1}{\sqrt{2}}(S_l^{m,-n}) + (-1)^n S_l^{m,n}), &n>0 \\
S_l^{m,n}, &n=0 \\
\frac{1}{\sqrt{2}}(-T_l^{m,n} + (-1)^n T_l^{m,-n}), &n<0
\end{cases} \\
S_l^{m,n} &= \begin{cases}
\frac{1}{\sqrt{2}}(F_l^{-m,n} + (-1)^m F_l^{m,n}), &m>0 \\
F_l^{m,n}, &m=0 \\
\frac{1}{\sqrt{2}}(G_l^{m,n} - (-1)^m G_l^{-m,n}), &m<0
\end{cases} \\
T_l^{m,n} &= \begin{cases}
\frac{1}{\sqrt{2}}(G_l^{-m,n} + (-1)^m G_l^{m,n}), &m>0 \\
G_l^{m,n}, &m=0 \\
\frac{1}{\sqrt{2}}(-F_l^{m,n} + (-1)^m F_l^{-m,n}), &m<0
\end{cases}
\end{align*}$$
This is a direct evaluation of $R = W D W^*$ (where $W^*$ is the hermitian conjugate of $W$).