Introduction
Spherical harmonics are a set of functions that form an orthonormal basis of the space of squareintegrable functions on the threedimensional unit sphere. They play an important role in many scientific fields like chemistry, physics, and computer graphics. I will assume the reader is familiar with the mathematics and focus on describing a small C++ class that puts the theory into practice. The articles [1] and [2] provide a brief introduction into the theory behind spherical harmonics and illustrate their use with exemplary applications. Also note that I will only consider realvalued spherical harmonics. Whenever I refer to spherical harmonics in the following text, I implicitly mean real spherical harmonics.
The intention of this article is to present a class that users new to the topic can use to play around with and explore the field. (This is what I used it for.) Therefore, the class' interface is rather simple and only contains a few basic operations. It is also not directly suitable to be used in highaccuracy / highperformance numerical calculations.
The next section will briefly lay down the requirements for the class. The next section presents the most important parts of the code. The article then closes with a short summary.
Requirements
A class representing a spherical harmonic projection should fulfill the following requirements:

Construction:
Instances of the class should be constructable from function pointers or functors. The functions must take two real numbers as arguments (spherical coordinates) and return a real number (the function value). Upon construction, the function should be projected into SH coefficients by numerical integration. The user should be able to choose the integration method, numerical accuracy, and the order of approximation (i.e. the number of SH bands used). After construction, the function should not be used by the class.
For functions that have already been projected (e.g. because it can be done analytically for them), a constructor that takes a set of SH coefficients should be provided.
As the class should have value semantics, a copy constructor and a copy assignment operator should be supplied.
Implementation
This section will only explain the most important parts of the code. The full source code is available in the accompanying archive.
The class stub to start with is very simple. It is parameterizable and takes the floating point type as a template parameter. It has the number of bands and the SH coefficients as attributes:
template< typename TReal >
class SHProjection
{
public:
typedef std::vector<TReal> SHCoefficients;
private:
int m_numBands; SHCoefficients m_coefficients; }
The number of bands used is not a template parameter since it may be changed by some operations on the class. Then we have the basic accessors with read/write access to single coefficients but readonly access to the number of bands. The code is trivial but for completeness, here are the declarations:
const SHCoefficients& GetProjection() const;
const TReal GetCoefficient( int l, int m ) const;
void SetCoefficient( int l, int m, TReal value );
const TReal GetCoefficient( int i ) const;
void SetCoefficient( int i, TReal value );
const TReal& operator [] ( int i ) const;
TReal& operator [] ( int i );
const int GetNumBands() const;
The class has two constructors besides the trivial copy constructor. The first one constructs a SHProjection
from a vector of coefficients. If the vector has not enough coefficients to fill the highest band completely, the remaining coefficients are assumed to be zero:
template< typename TReal >
inline
SHProjection< TReal >::SHProjection( const SHCoefficients& coeffs )
{
TReal temp = ceil( sqrt( static_cast< TReal >( coeffs.size() ) ) );
m_numBands = static_cast<int>( temp );
m_coefficients.resize( m_numBands * m_numBands, 0.0 );
std::copy( coeffs.begin(), coeffs.end(), m_coefficients.begin() );
}
The second constructor takes a function to project and some parameters that control how this is to be done. It is itself a template parameterizable by the function's type and simply forwards all parameters to the private
member function project
:
template< typename TReal >
template< class TFunction >
inline
SHProjection< TReal >::SHProjection(
const TFunction& function,
int numBands,
int numSamples,
IntegrationMethods integrationMethod )
:
m_numBands( numBands )
{
project( function, numBands, numSamples, integrationMethod );
}
project
then performs numerical integration to calculate the coefficients. Currently, there are three methods to choose from:
enum IntegrationMethods
{
IntegrationMethod_MC, IntegrationMethod_JitteredMC, IntegrationMethod_Simpson,
IntegrationMethod_NumOf
};
The plain Monte Carlo integrator is generally less efficient than the other two, but it has a simple implementation and therefore is shown here for illustrative purposes:
template< typename TReal >
template< class TFunction >
void SHProjection< TReal >::project_MC( const TFunction& function, int numBands,
int numSamples )
{
const int numCoeffs = numBands*numBands;
m_coefficients.clear();
m_coefficients.resize( numCoeffs, 0 );
for( int i=0; i<numSamples; ++i )
{
TReal xi_1 = frand(); TReal xi_2 = frand();
TReal theta = 2 * acos( sqrt( xi_1 ) ); TReal phi = 2 * PI * xi_2;
for( int l=0; l<numBands; ++l )
{
for( int m=l; m<=l; ++m )
{
int index = l*(l+1) + m;
m_coefficients[index] += Y(l,m,theta,phi)
* function(theta, phi); }
}
}
const TReal normalizationFactor = 4*PI / numSamples;
for( int k=0; k<numCoeffs; ++k )
m_coefficients[k] = m_coefficients[k] * normalizationFactor;
}
Y(l,m,theta,phi)
is a static
member function that evaluates the mth sphercial harmonic of band l at (theta, phi). (For code, see the accompanying archive.)
Jittered sampling is (asymptotically) always at least as good as plain Monte Carlo integration. In many cases, it can decrease variance of the estimator and improve efficiency. For smooth functions with high degrees of continuity (C2 or, better yet, C4) the Simpson's rule should be used, because it has the best convergence. The code is not complicated but quite lengthy, so it is not shown here.
After a function has been projected, function values can be queried via two methods. One takes spherical coordinates and the other Cartesian coordinates. The SH basis functions are then evaluated at these coordinates to reconstruct the function in the usual way:
template< typename TReal >
TReal SHProjection< TReal >::Evaluate( TReal theta, TReal phi ) const
{
TReal result = 0.0;
for( int l=0; l<m_numBands; ++l )
for( int m=l; m<=l; ++m )
result += GetCoefficient(l,m) * Y(l,m,theta,phi);
return result;
}
template< typename TReal >
TReal SHProjection< TReal >::Evaluate( TReal x, TReal y, TReal z ) const
{
TReal theta = acos(z);
TReal t = x / sqrt(x*x + y*y);
TReal phi = y >= 0 ? acos(t) : 2*PI  acos(t);
return Evaluate( theta, phi );
}
Besides a simple evaluation, there are some componentwise operations available. They all have trivial implementations:
const SHProjection& operator += ( const SHProjection& other );
const SHProjection& operator *= ( const SHProjection& other );
const SHProjection& operator *= ( TReal factor );
friend const SHProjection operator + <> ( const SHProjection& lhs,
const SHProjection& rhs );
friend const SHProjection operator * <> ( const SHProjection& lhs,
const SHProjection& rhs );
template < typename TReal, class TFunction, typename TScalar >
friend const SHProjection operator * <> ( TScalar lhs, const SHProjection& rhs );
template < typename TReal, class TFunction, typename TScalar >
friend const SHProjection operator * <> ( const SHProjection& lhs, TReal rhs );
Two more interesting operations are the inner product of two vectors of coefficients and, related to that, the convolution of the two projected functions. Because both operations are quite expensive when carried out on the original functions, they are often the main reason to use spherical harmonic expansion. The inner product calculates the integral of the product of the two projected functions (so it returns a scalar). The convolution operation transforms a SHProjection
into its convolution with a kernel function, so it can be queried with multiple coordinates.
template< typename TReal>
inline
TReal SHProjection< TReal >::InnerProduct( const SHProjection< TReal >& lhs,
const SHProjection< TReal >& rhs )
{
const std::size_t n = std::min( lhs.m_coefficients.size(),
rhs.m_coefficients.size() );
return std::inner_product(
lhs.m_coefficients.begin(), lhs.m_coefficients.begin() + n,
rhs.m_coefficients.begin(),
TReal() );
}
template< typename TReal>
inline
void SHProjection< TReal >::Convolve
( const SHProjection< TReal >& kernelProjection )
{
const int numBands = std::min( m_numBands, kernelProjection.m_numBands );
for( int l=0; l<numBands; ++l )
{
const TReal alpha = sqrt( 4*PI / (2*l + 1) );
for( int m=l; m<=l; ++m )
{
m_coefficients[ l*(l+1) + m ] *= alpha *
kernelProjection.GetCoefficient(l,
0);
}
}
}
Summary
A small class that allows for a compact representation of a spherical harmonic expansion was presented. It is not productionready, but intended to illustrate concepts and serve as a starting point for people who want to dive deeper into the topic. Nonetheless, the most common operations are supported and performance should be sufficient for most applications.
References
[1] Robin Green  Spherical Harmonic Lighting: The Gritty Details. http://www.research.scea.com/gdc2003/sphericalharmoniclighting.html. Archives of the Game Developers Conference. 2003.
[2] Volker Schönefeld. Spherical Harmonics. http://heim.cotto.de/~volker/prosem_paper.pdf. Seminal paper. 2005.
History
 17^{th} March, 2008: Initial post