Spherical Harmonics(SH) functions are a set of orthogonal basis functions defined in spherical coordinates using imaginary numbers. In this post, we use the following conversion between spherical and cartesian coordinates:
Since we are dealing with real value functions, we only need to deal with real spherical harmonics functions which in the form of:
The index l of the SH function is called the band index which is an integer >= 0 and index m is an integer with range -l<=m<=l , so there will be (2l + 1) functions in a given band. You may refer to the Appendix A2 of Stupid Spherical Harmonics(SH) Trick to look up the evaluated value of the SH basis function for a pair of (l, m).
The linear combination of SH basis functions with scalar values can be used to approximate a function as below:
With an approximation up to band l = n - 1, which n×n coefficients are needed.
So the remaining problem to approximate a function is to compute the coefficient c which can be solved either analytically or numerically by Monte Carlo Integration.
For shading lambert diffuse surface without shadow, we can simplify the rendering equation into:
To solve this integral, we can project the functions L(x, ω) and max(N.ω, 0) into SH functions using Monte Carlo Integration, then by the property 2 described above, the integral equals to dot product of the SH coefficients of the 2 SH projected functions.
If a SH projected function is rotational symmetric about a fixed axis, it is called Zonal Harmonics(ZH). If this axis is the z-axis, this will make the ZH function only depends on θ, which will result in only one non-zero coefficient in each band with m= 0. Then rotation of the ZH function can be greatly simplified. When the ZH function is rotated to a new axis d, the coefficients of the rotated SH function will equals to:
,which is faster than the general SH rotation. The ZH function is well suit to approximate the function max(N.ω, 0) in the above diffuse surface rendering equation since the SH projected L(x, ω) is usually done in world space while the shading surface can be re-oriented to the same space to perform lighting calculation.
Here is a webGL demo (which need a webGL enabled browser such as Chrome) using the cube map below as light source and projected to SH function using Monte Carlo Integration.
Both the white and the blue color on the model is reflected from the sun and the blue sky using SH coefficients generated from the cube map and the ZH coefficients projected from max(N.ω, 0) which rotated to world space according the surface normal. The approximation is done up to band l=2. You can drag in the viewport to rotate the camera. Below is a screen shot of the demo:
The source code of the webGL can be downloaded here.
SH functions can be used to approximate the rendering equation with only a few coefficients and a simple dot product to evaluate lighting during run time. But it also has its disadvantage while SH can only approximate low frequency function as it needs large number of bands to represent high frequency details.
 Spherical Harmonics Lighting: The Gritty Details: http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
 Stupid Spherical Harmonics(SH) Trick: http://www.ppsloan.org/publications/StupidSH36.pdf
 Physically Based Rendering: http://www.amazon.com/gp/product/0123750792/ref=pd_lpo_k2_dp_sr_1?pf_rd_p=486539851&pf_rd_s=lpo-top-stripe-1&pf_rd_t=201&pf_rd_i=012553180X&pf_rd_m=ATVPDKIKX0DER&pf_rd_r=09AG8FQQWKJHC2AEFPD1
 Sky box texture downloaded from: http://www.codemonsters.de/home/content.php?show=cubemaps