A Simple Class for Spherical Harmonic Expansion

• Download source - 8.82 KB

Introduction

Spherical harmonics are a set of functions that form an orthonormal basis of the space of square-integrable functions on the three-dimensional 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 real-valued 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 high-accuracy / high-performance 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.

• Parameterizable with floating point type:

  The class should be able to use different floating point data types for storing the SH coefficients of the projection.

• Operations:

  The following operations should be provided (where appropriate by operator overloading):

• Component-wise addition (represents an addition of the projected functions)
• Component-wise multiplication (useful as a building block for convolution)
• Multiplication by a scalar value (represents a scaling of the projected function)
• Accessors for the SH coefficients (in terms of band and band index (m,l pairs) or in terms of the canonical one-dimensional index (i=(l+1)l+m))
• Outputting the projection to a stream in a formatted way
• Reconstruction of the projected function for points given in spherical and Cartesian coordinates (function evaluation)
• Inner product (yields (an approximation of) the integral of the product of the two projected functions; a very useful property that is the basis of all SH lighting algorithms)
• Convolution with a kernel function also given in SH (can be used to filter the projected function)

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:
  /// Vector of TReals representing the SH-projection of a function
  typedef std::vector<TReal> SHCoefficients;

private:
  int m_numBands;   //< number of bands used for this projection (order of approximation)
  SHCoefficients m_coefficients;    //< the SH coefficients of this projection
}

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 read-only access to the number of bands. The code is trivial but for completeness, here are the declarations:

/// Returns the SH coefficients of this projection.
const SHCoefficients& GetProjection() const;

/// Returns the SH coefficient belonging to Y(l,m).
const TReal GetCoefficient( int l, int m ) const;

/// Sets the SH coefficient belonging to Y(l,m).
void SetCoefficient( int l, int m, TReal value );

/// Returns the i-th SH coefficient.
const TReal GetCoefficient( int i ) const;

/// Sets the i-th SH coefficient.
void SetCoefficient( int i, TReal value );

/// Const access to the i-th SH coefficient.
const TReal& operator [] ( int i ) const;

/// Non-const access to the i-th SH coefficient.
TReal& operator [] ( int i );

/// Returns the number of bands used for this projection (ie the order of approximation).
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 )
{
    // init m_numBands
    TReal temp = ceil( sqrt( static_cast< TReal >( coeffs.size() ) ) );
    m_numBands = static_cast<int>( temp );

    // copy coefficients over (empty places in the highest band will be set to zero)
    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:

/// Enumeration of available numerical integration methods
enum IntegrationMethods
{
    IntegrationMethod_MC,            //< plain Monte Carlo integration over the
                                     //sphere's surface
    IntegrationMethod_JitteredMC,    //< Monte Carlo integration using jittered sampling
    IntegrationMethod_Simpson,       //< integration using Simpson's rule

    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;

    // prepare the coefficients vector
    m_coefficients.clear();
    m_coefficients.resize( numCoeffs, 0 );

    for( int i=0; i<numSamples; ++i )
    {
        TReal xi_1 = frand();    // get two independent, uniformly distributed
        TReal xi_2 = frand();    // random numbers in [0,1) ...

        TReal theta = 2 * acos( sqrt( xi_1 ) );        // ... and generate points in
                                                       // spherical coordinates
        TReal phi = 2 * PI * xi_2;                     // uniformly distributed over
                                                       // the sphere's surface.

        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); // take
                // a sample and add it to the sum
            }
        }
    }

    // divide by the number of samples and multiply by the area
    // of the integration region to get the final result
    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 m-th 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
{
    // first get the spherical coordinates ...
    TReal theta = acos(z);
    TReal t = x / sqrt(x*x + y*y);
    TReal phi = y >= 0 ? acos(t) : 2*PI - acos(t);

    // ... then call the spherical evaluate
    return Evaluate( theta, phi );
}

Besides a simple evaluation, there are some component-wise operations available. They all have trivial implementations:

/// Adds another SHProjection to this one (component-wise).
/// If other has a higher order this SHProjection's order is increased.
/// This is equivalent to a addition of the represented functions.
const SHProjection& operator += ( const SHProjection& other );

/// Multiplies this SHProjection by another one (component-wise).
/// If other has a higher order this SHProjection's order is increased.
const SHProjection& operator *= ( const SHProjection& other );

/// Multiplies this SHProjection by a constant factor (component-wise).
/// This is equivalent to a scaling of the represented function.
const SHProjection& operator *= ( TReal factor );

/// Performs component-wise addition of two SHProjections.
/// The order of the result is the maximum order of the arguments.
/// This is equivalent to an addition of the represented functions.
friend const SHProjection operator + <> ( const SHProjection& lhs,
    const SHProjection& rhs );

/// Performs component-wise multiplication of two SHProjections.
/// The order of the result is the maximum order of the arguments.
friend const SHProjection operator * <> ( const SHProjection& lhs,
    const SHProjection& rhs );

/// Multiplies a SHProjection by a constant factor.
/// This is equivalent to a scaling of the represented function.
template < typename TReal, class TFunction, typename TScalar >
friend const SHProjection operator * <> ( TScalar lhs, const SHProjection& rhs );

/// Multiplies a SHProjection by a constant factor.
/// This is equivalent to a scaling of the represented function.
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 production-ready, 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/spherical-harmonic-lighting.html. Archives of the Game Developers Conference. 2003.

[2] Volker Schönefeld. Spherical Harmonics. http://heim.c-otto.de/~volker/prosem_paper.pdf. Seminal paper. 2005.

History

• 17^th March, 2008: Initial post

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)