Stupid Spherical Harmonics (SH)
TricksPeter-Pike Sloan
Microsoft Corporation

Abstract
This paper is a companion to a GDC 2008 Lecture with the same title. It provides a brief
overview of spherical harmonics (SH) and discusses several ways they can be used in interactive
graphics and problems that might arise. In particular it focuses on the following issues: How toevaluate lighting models efficiently using SH, what “ringing” is and what you can do about it,
efficient evaluation of SH products and where they might be used. The most up to date version
is available on the web at http://www.ppsloan.org/publications

Introduction
Harmonic functions [2], the solutions to Laplace’s equation, are used extensively in various
fields. Spherical Harmonics are the solutions when restricted to the sphere 1. They have beenused to solve potential problems in physics, such as the heat equation (modeling the variation oftemperature over time [5][25]), and the gravitational and electric fields[9]. They have also beenused in quantum chemistry and physics to model the electron configuration in atoms and model
quantum angular momentum [16][51]. Closer to graphics they have been used to model
scattering phenomena [7][17]. In computer graphics they have been extensively used, earlyuses were in modeling volumetric scattering effects [18], environmental reflections for microfacet BRDF’s without global shadows[6], non-diffuse off-line light transport simulations[40],
BRDF representations [53], image relighting[28], image based rendering with controllable
lighting [54][55], and modeling light source emission[8]. More recent examples include morework in atmospheric scattering [50] and computer vision [3].
The focus of this article is on techniques related to interactive rendering. The first paper thathas been used extensively in games deals with using Spherical Harmonics to represent
irradiance environment maps efficiently, allowing for interactive rendering of diffuse objects
under distant illumination [35]. This was extended to handle a limited class of BRDF’s with thesame constraints [36]. Precomputed radiance transfer (PRT) [41][20][24] models a static
object/scenes response to a lighting environment, often represented using SH, including
complex global illumination effects like soft-shadows and inter-reflections with diffuse andsimple glossy materials. This was extended to handle more general BRDF models [20][23] [42],1

They are the angular portion of the solution in spherical coordinates.

incorporate sub-surface scattering [42], significantly improve rendering efficiency throughcompression techniques from machine learning [42], and various techniques to model “local”
texture like surface details [43][44][45]. SH have been used to model single scattering fromdistant lighting environments [49]. Other uses have been in using gradients to lift the distantlighting assumption [1], several techniques addressing dynamic objects [56][37] including
support for inter-reflections [46][33], as a representation of visibility to model shadowing ofobjects with general BRDF models [12], using scaling operators to model shadows fromdeformable objects [52], as a parameterization of refraction [11] and as a technique to address
the level of detail problem with normal maps [15].
More practical papers include covering implementation details for PRT [13], how to integratethese techniques into an engine [30], how to use SH+gradients for irradiance volumes [31][32],practical issues around projection and how to efficiently quantize SH coefficients [21] and a nicepaper that projects an analytic skylight model [34] into SH and uses a global polynomial fit toevaluate the SH lightprobe as a function of the models parameters [14]. Numerical techniques
for more robustly projection functions defined over the hemisphere using SH have also been
investigated [22].
Many of the uses in real-time graphics are just as a convenient way for representing sphericalfunctions – visibility, lighting and reflectance. While there are many other basis functions thatcan be used, wavelets [39], wavelets on cube maps [27], spherical radial basis functions [9], andothers [26], spherical harmonics have some nice properties that will be described in thisdocument. It is important to stress that there are scenarios where these other basis functions
are more appropriate.
While spherical harmonics may seem somewhat daunting, they are actually straightforward.
They are the spherical analog to the Fourier basis on the unit circle, and are easy to evaluate
numerically. Like the Fourier basis used in signal processing, care has to be taken whentruncating the series (which will always be done in video games), to minimize the “ringing”
artifacts that can occur. This article will describe how to evaluate and represent lights efficiently
using spherical harmonics, how to pull conventional lights out of a SH representation, describe“ringing” and mitigation techniques to minimize its impact, and go over products of functions
using spherical harmonics, describing where they are useful and special cases that are worth
optimizing for.

BackgroundDefinition Spherical Harmonics define an orthonormal 2 basis over the sphere, S. Using theparameterization
𝑠𝑠 = (𝑥𝑥, 𝑦𝑦, 𝑧𝑧) = (sin 𝜃𝜃 cos 𝜑𝜑 , sin 𝜃𝜃 sin 𝜑𝜑 , cos 𝜃𝜃)2

An orthonormal basis is one that has the property ∫ 𝑓𝑓𝑖𝑖 𝑓𝑓𝑗𝑗 = 𝛿𝛿𝑖𝑖𝑖𝑖 , that is one if i==j, and zero otherwise.

Where s are simply locations on the unit sphere. The basis functions are defined as|𝑚𝑚 |𝑌𝑌𝑙𝑙𝑚𝑚 (𝜃𝜃, 𝜑𝜑) = 𝐾𝐾𝑙𝑙𝑚𝑚 𝑒𝑒 𝑖𝑖𝑖𝑖𝑖𝑖 𝑃𝑃𝑙𝑙(cos 𝜃𝜃), 𝑙𝑙 ∈ 𝐍𝐍, −𝑙𝑙 ≤ 𝑚𝑚 ≤ 𝑙𝑙

Where 𝑃𝑃𝑙𝑙𝑚𝑚 are the associated Legendre polynomials and 𝐾𝐾𝑙𝑙𝑚𝑚 are the normalization constants
𝐾𝐾𝑙𝑙𝑚𝑚 = �(2𝑙𝑙 + 1)(𝑙𝑙 − |𝑚𝑚|)!
4𝜋𝜋(𝑙𝑙 + |𝑚𝑚|)!

The above definition is for the complex form (most commonly used in the non-graphics
literature), a real valued basis is given by the transformation√2𝐾𝐾𝑙𝑙𝑚𝑚 cos 𝑚𝑚𝑚𝑚 𝑃𝑃𝑙𝑙𝑚𝑚 (cos 𝜃𝜃) 𝑚𝑚 > 0√2Re(𝑌𝑌𝑙𝑙𝑚𝑚 ) 𝑚𝑚 > 0𝑦𝑦𝑙𝑙𝑚𝑚 = �√2Im(𝑌𝑌𝑙𝑙𝑚𝑚 ) 𝑚𝑚 < 0 = �√2𝐾𝐾𝑙𝑙𝑚𝑚 sin|𝑚𝑚|𝜑𝜑 𝑃𝑃𝑙𝑙|𝑚𝑚 | (cos 𝜃𝜃) 𝑚𝑚 < 0
𝑚𝑚 = 0𝑚𝑚 = 0
𝑌𝑌𝑙𝑙0
𝐾𝐾𝑙𝑙0 𝑃𝑃𝑙𝑙0 (cos 𝜃𝜃)The index l represents the “band”. Each band is equivalent to polynomials of that degree (sozero is just a constant function, 1 is linear, etc.) and there are 2l+1 functions in a given band.
While spherical coordinates are convenient when computing integrals, they can also be
represented using polynomials, as is commonly done when evaluating them (see Appendix A1
Recursive Rules for Evaluating SH Basis Functions and Appendix A2 Polynomial Forms of SH Basisfor details.) An order n SH expansion uses all of the basis functions through degree n-1.

Spherical Harmonics can be visualized in a couple of ways. One standard way is to distort a unit
sphere, by scaling each point radialy by the absolute value of the function and coloring it basedon the sign (red for positive, blue for negative.) Above are images of the first three bands using
this technique.
The functions in the central column (l=0) are called zonal harmonics (ZH) and will be discussed
later, these functions have rotational symmetry around the z axis and the zeros (locations where
the function is zero) are contours on the sphere parallel to the XY plane. The functions where(l=|m|) are called sectorial harmonics and the zeros define regions like apple slices.
An alternative visualization is to draw them using the parameterization of a cube map unfolded
onto the plane. The unfolding of the cube map is as follows:+Y-X -Z +x-Y
+Z
Here magnitude is encoded with color (red positive, blue negative, zero green) and iso-intensity
contours have been evenly placed (white lines) to give more intuition for the gradient of thefunction (when they bunch together the function is changing faster, etc.)Projection and Reconstruction Because the SH basis is orthonormal the least squares
projection of a scalar function f defined over S is done by simply integrating the function you
want to project, 𝑓𝑓(𝑠𝑠), against the basis functions (proof in Appendix A6 Least Squares
Projection)
𝑓𝑓𝑙𝑙𝑚𝑚 = � 𝑓𝑓(𝑠𝑠)𝑦𝑦𝑙𝑙𝑚𝑚 (𝑠𝑠) 𝑑𝑑𝑑𝑑

These coefficients can be used to reconstruct an approximation of the function f𝑛𝑛𝑙𝑙𝑓𝑓̃(𝑠𝑠) = � � 𝑓𝑓𝑙𝑙𝑚𝑚 𝑦𝑦𝑙𝑙𝑚𝑚 (𝑠𝑠)𝑙𝑙=0 𝑚𝑚 =−𝑙𝑙

Which is increasingly accurate as the number of bands n increases. This paper concentrates on
low-frequency approximations to f, for higher frequency representations other bases tend to doa better job. Projection to n-th order generates n2 coefficients. It is often convenient to use asingle index for both the projection coefficients and the basis function, via
𝑛𝑛 2𝑓𝑓̃(𝑠𝑠) = � 𝑓𝑓𝑖𝑖 𝑦𝑦𝑖𝑖 (𝑠𝑠)𝑖𝑖=0Where i=l(l+1)+m . This formulation makes it clear that evaluating at direction s of theapproximate function is simply a dot product between the n2 coefficient vector fi and the vectorof evaluated basis functions yi(s). The first coefficient (𝑓𝑓00 or 𝑓𝑓0 depending on indexing)
represents the average value of the function over the sphere and will sometimes be referred toas the DC term.
Basic Properties One important property of SH is how projections interact with rotations.
Given a function g(s), which represents a function f(s) rotated by a rotation matrix Q, so g(s) =f(Q(s)) the projection of g is identical to rotating 𝑓𝑓̃ and re-projecting it. This rotational
invariance is similar to the translational invariance in the Fourier transform. This means that, forexample, lighting will be stable under rotations, so there won’t be any aliasing artifacts or“wobbling” of the light sources. Below are images of a sphere illuminated by a directional lightsource, the top row is using SH, the bottom row is using the Ambient Cube basis 3 from Valve[26]. The first column is a best case orientation, and the second is near a worst case one. The
image is invariant using SH. This basis is discussed in more details in Appendix A9 Ambient Cube
Basis. This will happen to some extent with any other basis defined over the sphere.

3

This basis is more efficient to evaluate compared to SH, having 6 basis functions with only 3 being nonzero at any point on the sphere.

SH

HL2

“good”“bad”

Due to the orthonormality of the SH basis, given any two SH functions a and b, the integral of2
the product is simply the dot product of the coefficient vectors: ∫ 𝑎𝑎�(𝑠𝑠)𝑏𝑏�(𝑠𝑠) = ∑𝑛𝑛𝑖𝑖=0 𝑎𝑎𝑖𝑖 𝑏𝑏𝑖𝑖

Convolution Given a kernel function h(z) that has circular symmetry, you can generate a new SH
function that is the result of convolving the kernel with an original function f. h must havecircular symmetry for the result of the convolution to also be represented on the sphere S,
instead of the rotation group SO(3). The convolution can be done directly in the frequency
domain using the following equation:
4𝜋𝜋 0 𝑚𝑚(ℎ ∗ 𝑓𝑓)𝑚𝑚ℎ 𝑓𝑓
𝑙𝑙 = �
2𝑙𝑙 + 1 𝑙𝑙 𝑙𝑙This amounts to simply scaling each band of f by the corresponding m=0 term from h.
Rotation As mentioned before, the SH are closed under rotation. SH rotation matrices are in a
block structure, where each band is rotationally independent and has a dense (2l+1)x(2l+1) submatrix. There are several ways to compute these rotation matrices, for very small orders(quadratic and less) doing so symbolically is most efficient, but for higher orders it seems to be
more efficient to decompose the rotation matrix into zyz Euler angles [19].
Zonal Harmonics Spherical harmonic projections of functions that have rotational symmetryaround an axis are called Zonal Harmonics (ZH.) If they are oriented so that this axis is Z, thezeros of the function will form lines of constant latitude and the functions only depend on 𝜃𝜃.
The coefficient vector in this orientation only has one non-zero per band, so a n-th orderfunction has n instead of n2 coefficients. Zonal Harmonics have been used to approximatetransfer [44], and are the common representation of phase functions in scattering theory[7][17], they will be used extensively in this paper when modeling light sources. Rotation of
Zonal Harmonics is simpler than general SH, it can be done with what is effectively a diagonal
matrix and only requires evaluating the SH basis functions in the new direction d. Given the ZH
coefficients of a function (only the m=0 terms from an SH projection) zl it can be rotated to a
new direction d using this equation:
4𝜋𝜋𝑓𝑓(𝑠𝑠) = � 𝑧𝑧𝑙𝑙 �� 𝑦𝑦𝑙𝑙𝑚𝑚 (𝑑𝑑)𝑦𝑦𝑙𝑙𝑚𝑚 (𝑠𝑠)2𝑙𝑙 + 1𝑙𝑙

So the resulting SH coefficients are:𝑚𝑚

4𝜋𝜋
𝑓𝑓𝑙𝑙𝑚𝑚 = �
𝑧𝑧 𝑦𝑦 𝑚𝑚 (𝑑𝑑)2𝑙𝑙 + 1 𝑙𝑙 𝑙𝑙

SH Products The kth coefficient of the product of two functions f and g represented using nthorder SH projected into SH has the following form:
𝑛𝑛 2

𝑛𝑛 2

𝑖𝑖=0

𝑗𝑗 =0

𝑝𝑝𝑘𝑘 = � 𝑦𝑦𝑘𝑘 (𝑠𝑠) �� 𝑓𝑓𝑖𝑖 𝑦𝑦𝑖𝑖 (𝑠𝑠)� �� 𝑔𝑔𝑗𝑗 𝑦𝑦𝑗𝑗 (𝑠𝑠)� 𝑑𝑑𝑑𝑑 = � Γ𝑖𝑖𝑖𝑖𝑖𝑖 𝑓𝑓𝑖𝑖 𝑔𝑔𝑗𝑗
Where Γ is the triple product tensor:𝑖𝑖𝑖𝑖

Γ𝑖𝑖𝑖𝑖𝑖𝑖 = � 𝑦𝑦𝑖𝑖 (𝑠𝑠) 𝑦𝑦𝑗𝑗 (𝑠𝑠)𝑦𝑦𝑘𝑘 (𝑠𝑠)𝑑𝑑𝑑𝑑

an order 3 sparse symmetric tensor. Since SH are polynomials, the polynomial product will havemaximal degree 2n-2, which means it will have non-zero coefficients through order 2n-1. This
becomes unwieldy as the number of functions being multiplied grows, so it is common totruncate the product early [56][37].The number of non-zero coefficients as a function of n isquite large [47][37] so care has to be taken when generating efficient code. One special casethat is useful to point out is that if the function f is fixed (ie: distant lighting) you can compute a“product matrix” Mf which will significantly reduce the cost. This matrix is symmetric and built
using the following equation:(Mf )ij = � Γ𝑖𝑖𝑖𝑖𝑖𝑖 𝑓𝑓𝑘𝑘kComputing the product with a function g in this case is simply a matrix vector product.

Irradiance Environment Maps
An irradiance environment map is created by convolving a light probe with a clamped cosine
function; this should be normalized by dividing by 𝜋𝜋 to display radiance. This convolution can bedone efficiently using SH [35], and is accurate enough to be efficiently rendered directly from SH
as well. Order 3 SH do a good job approximating this kernel, but if HDR light sources are going
to be used you might want to consider using Order 5 (the order 4 ZH coefficient is zero so thatband can be skipped.) Below are images of the clamped cosine kernel and the order 3 SH
approximation, the red curve is the SH approximation, the figure on the left is a plot as afunction of theta, on the right a polar plot scaled by the absolute value of the function:

Below are plots that also includes the Order 5 projection (blue):

The Order 3 SH approximation over estimates by 1/16th at theta=0 (north pole) and has aspurious lobe at the south pole with magnitude of 1/16th. A directional light source that would
reflect a value of 17 when a normal points at it would reflect a value of 1 pointing in the
opposite direction (should be reflecting 0.) The order 5 approximation has a negative lobe,
which would reflect -1 with a directional light that would reflect 31 with a normal pointing rightat it. While these approximations are accurate, the approximation can cause error, particularly
if very bright light sources are being used.
Appendix A10 Shader/CPU code for Irradiance Environment Maps contains shader and CPU codefor efficiently evaluating irradiance environment maps.

Lighting Models
There are various ways of representing lighting in SH. The simplest is to just project from a cube
map, but there are also analytic models that are inexpensive to evaluate and potentially useful
to expose to artists. A recent paper [14] does a nice job of projecting a practical skylight model[34] into SH and fits a global polynomial of the SH coefficients over the parameter space of themodel.

Projection from Cube Maps
To project from a cube map you simply need to integrate the SH basis functions against the cube
map. This can be done numerically by evaluate the SH basis functions in the direction of eachtexel center, weight it by the differential solid angle for that texel and normalize the results. In
pseudo-code:
float f[],s[];float fWtSum=0;
Foreach(cube map face)
Foreach(texel)
float fTmp = 1 + u^2+v^2;
float fWt = 4/(sqrt(fTmp)*fTmp);
EvalSHBasis(texel,s);
f += t(texel)*fWt*s; // vector
fWtSum += fWt;
f *= 4*Pi/fWtSum; // area of sphere

In the code above u and v represent the two coordinates on the given face that are not +/- 1 4,t(texel) is the texel color. EvalSHBasis needs to normalize the input coordinate (on a cube face),and simply evaluate the SH basis functions in that direction. The last normalization can beomitted (instead you can just divide by the number of samples), since the normalized sum of the4

Ie: on the +X face, it would be y/z, etc.

differential solid angles should be 4*Pi, but it tends to be a little off (particularly when using low
resolution cube maps.)
Below are images of the reconstruction of a HDR light probe into order 1 to 6 Spherical
Harmonics. The final image is the light probe that was projected.

…

Analytic Models
Directional lights are trivial to compute, you simply evaluate the SH basis functions in the given
direction and scale appropriately (see Normalization section.) Spherical Light sources can be
efficiently evaluated using zonal harmonics. Below is a diagram showing an example scene, wewant to compute the incident radiance, a spherical function,at the receiver point P. Given a spherical light source withr
center C, radius r, what is the radiance arriving at a point P dCunits away? The sin of the half-angle subtended by the lightsource is r/d, so you just need to compute a light source that
subtends an appropriate part of the sphere. The ZH
coefficients can be computed in closed form as a function of
d
𝑎𝑎

2𝜋𝜋

this angle: 𝑧𝑧𝑙𝑙 = ∫𝜃𝜃=0 ∫𝜑𝜑=0 𝑦𝑦𝑙𝑙0 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 where a is the half-angle
subtended. See Appendix A3 ZH Coefficients for Spherical
Light Source for the expressions through order 6.

P
Figure 1 Spherical Light Source

The spherical technique can also be used to model a cone (think of it as a disc at infinity) withconstant emission. A softer cone can be modeled that has a smooth fall-off over the visible
portion5 – see Appendix A4 ZH Coefficients for Smooth Cone for the equations.
The columns are for orders 4 and 6, the top row is the projection of a cone and the bottom row
is a cone that has a smooth fall off. The angles are 90 (green), 45 (red), 30 (blue), 12.5 (black)
degrees. The dashed lines represent the actual function (slightly displaced for the cone so they
don’t overlap), the solid lines are the SH approximation. In general the soft cone is better
behaved. How to address the artifacts that arise in projection is the subject of the section on
Ringing below.

Order 4
5

Order 6An “ease spline” is used, this is a cubic polynomial f(x) that has the constraints
32
f(0)=1,f’(0)=0,f(a)=0,f’(a)=0, the unique cubic polynomial is: 𝑓𝑓(𝑥𝑥) = 2𝑥𝑥 �𝑎𝑎3 − 3𝑥𝑥 �𝑎𝑎2 + 1

Normalization
If [0,1] lighting is being used, it is convenient to normalize the radiance vector so that thereflected radiance for an unshadowed receiver with a normal pointing directly at the lightsource would be 1.0. Mathematically you want to compute a scale factor c, that when
multiplied by the lighting vector L will result in unit reflected radiance when integrated against a
vector T that represents an unoccluded clamped cosine lobe (projection of normalized clamped
cosine into SH.) So you have:
𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦

1 = � 𝑐𝑐𝑐𝑐(𝑠𝑠) 𝑇𝑇(𝑠𝑠)𝑑𝑑𝑑𝑑 �⎯⎯⎯� 𝑐𝑐 =1
𝑑𝑑𝑑𝑑𝑑𝑑(𝐿𝐿, 𝑉𝑉)

Only the bands that are going to be used for rendering should be used when computing thenormalization factor. Aligning T with +Z gives a simple analytic formula; here are thecoefficients of T for the first 6 bands:
1

,

√3

,

√5

2√𝜋𝜋 3√𝜋𝜋 8√𝜋𝜋

, 0,

−1

16√𝜋𝜋

,0

For analytic lights you can use an analytic normalization term, for the cone light of angle a thatwould be:1
sin2 𝑎𝑎

However the resulting lights will not reflect unit radiance, because the projection error of boththe clamped cosine function and the light source will not be taken into account. For directional
lights 6 the normalization factor is

16𝜋𝜋
and for “ambient” lights it is 2√𝜋𝜋.
17Extracting Conventional Lights from SH

Given an SH lighting vector, it is possible to approximate it as a single directional light source
and an ambient light source. This has been used on hardware that does not support vertex
shaders. Mathematically we want to compute the intensity of a directional light (c) and theintensity of an ambient light source (a), so that squared error of reflected radiance is minimized
for any surface normal (N). Assuming a fixed direction (d) for the light source, the error function
we want to minimize is:2

𝐸𝐸(𝑐𝑐, 𝑎𝑎) = � �𝑐𝑐𝐻𝐻𝑁𝑁 (𝑑𝑑) + 𝑎𝑎 − � 𝐿𝐿𝑒𝑒 (𝑠𝑠)𝐻𝐻𝑁𝑁 (𝑠𝑠) 𝑑𝑑𝑑𝑑� 𝑑𝑑𝑑𝑑

Where 𝐻𝐻𝑁𝑁(𝑥𝑥) = 𝑚𝑚𝑚𝑚𝑚𝑚(𝑑𝑑𝑑𝑑𝑑𝑑(𝑁𝑁, 𝑥𝑥)⁄𝜋𝜋 , 0) is the normalized clamped cosine. If the lighting 𝐿𝐿𝑒𝑒 (𝑠𝑠) isexpressed in SH, this has a simple solution. The irradiance environment map represented by the
new lighting environment should be as close to the input irradiance environment map –

6

th

th

th

This is assuming no higher than 4 order lighting, for 5 or 6 order the normalization factor is

32𝜋𝜋
31

minimizing squared error of these two environment maps is equivalent to minimizing reflectedradiance over all normals:
2

𝐸𝐸(𝑐𝑐, 𝑎𝑎) = ��𝑐𝑐𝐿𝐿𝑑𝑑 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠) + 𝑎𝑎𝐿𝐿𝑎𝑎 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠) − 𝐿𝐿𝑒𝑒 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠)� 𝑑𝑑𝑑𝑑

Here 𝐿𝐿𝑑𝑑 (𝑠𝑠) is a normalized SH directional light in direction d, and 𝐿𝐿𝑎𝑎 (𝑠𝑠) is a normalized SH
constant light (just depends on DC.) The optimal values a and c are:867
𝑑𝑑𝑑𝑑𝑑𝑑0 �𝐿𝐿�𝑑𝑑 , 𝐿𝐿�𝑒𝑒 �316𝜋𝜋
8√𝜋𝜋 √𝜋𝜋𝑎𝑎 = �𝐿𝐿�𝑒𝑒 [0] − 𝑐𝑐
�
17 2
𝑐𝑐 =

The lighting vectors above are all turned into irradiance environment maps through convolving
with the normalized clamped cosine kernel. The dot product above ignores the DC term and𝐿𝐿�𝑒𝑒 [0] is the DC term of the lighting environment. The above assumes a known direction, a goodcandidate direction is the “optimal linear” direction [44], formed by normalizing the vector
� 𝑒𝑒 [3], −𝐿𝐿�𝑒𝑒 [1], 𝐿𝐿�𝑒𝑒 [2]), which are the linear coefficients of the lighting environment.
(−𝐿𝐿

Extracting Multiple Lights

It is also possible to extract multiple lights from a SH light probe. This might be done to combat
ringing (analytic lights won’t have negative lobes), to model glossy reflection (just from the lightsources pulled out in this manner) or to use a small number of shadow zbuffers (pull the lights
out for both diffuse and glossy parts of the BRDF.) The optimal linear direction works well whena light probe is dominated by a single light source, when this isn’t the case the directions andintensities of the lights should be optimized for in some fashion. One thing to do is to climb “uphill” to find local maxima of the function. Given how smooth SH are, for a given order there is a
finite distance where any point that distance from a distinct peak is guaranteed to reach thepoint using gradient ascent. A set of points can be generated with the property that the pointon the sphere farthest from the center of its Voronoi 7 cell is less then this distance. If you start
your search from each of these points you should find all of the point local maxima. These
distances can be found by looking at the projection of a delta function and computing theangular distance between peaks and the zeros. Using a conservative estimate of two thirds thisradius the number of points needed per order are {1,3,6,10,15,22} for the first 6 orders.
This point set can be computed by computing a simulation: All points have a force with a certain
fall-off (say 1/d2 where d is the Euclidean distance between them) that acts on every otherpoint. Generate the net force acting on each point, and subtract off the normal component (soit is just a tangential force), move the points by this force with a weight w. If the net sum of

7

Given a set of points X on the sphere, the Voronoi cell of a given point xi consists of all the points on thesphere that are closer to xi than any other point in the set.

forces increases, decrease w by half and try again, if it decreases double w and try again. This issolving for a set of electrons that minimize the sum of the electrostatic charges..
Given a SH lighting environment L(s), you want to find a point on the sphere that is a local
maximum, which is the same as finding a local minima of –L(s). This is a non-linear optimization
problem, a small number of BFGS iterations [29] converge on the peaks in my experience. The
gradients of the basis functions need to be computed when using optimization methods. For apoint on the sphere differentiating the polynomials is trivial, but you need to allow the point to
go off the sphere when doing a line search, so you want to use symbolic inputs that normalize
the coefficients (ie: 𝑓𝑓 �

𝑥𝑥

,

𝑦𝑦

,

𝑧𝑧

�𝑥𝑥 2 +𝑦𝑦 2 +𝑧𝑧 2 �𝑥𝑥 2 +𝑦𝑦 2 +𝑧𝑧 2 �𝑥𝑥 2 +𝑦𝑦 2 +𝑧𝑧 2

� in place of 𝑓𝑓(𝑥𝑥, 𝑦𝑦, 𝑧𝑧)). This is done bysimply computing the gradient at the normalized location, and multiplying it by the Jacobian of
the normalization function:

Where L is �𝑥𝑥 2 + 𝑦𝑦 2 + 𝑧𝑧 2 .

Original

𝐿𝐿2 − 𝑥𝑥 2⎡
⎢ 𝐿𝐿3
⎢ −𝑥𝑥𝑥𝑥⎢⎢ 𝐿𝐿⎢ −𝑥𝑥𝑥𝑥⎣ 𝐿𝐿

3 direction

−𝑥𝑥𝑥𝑥
𝐿𝐿
𝐿𝐿2 − 𝑦𝑦 2
𝐿𝐿3
−𝑦𝑦𝑦𝑦
𝐿𝐿

−𝑥𝑥𝑥𝑥
⎤
𝐿𝐿 ⎥
−𝑦𝑦𝑦𝑦 ⎥⎥𝐿𝐿 ⎥2
2
𝐿𝐿 − 𝑧𝑧 ⎥𝐿𝐿3 ⎦2 directions

1 directions

These images compare radiance (top row) and irradiance (bottom row) lighting environments
when approximating with 3, 2 or 1 directional lights and a ambient light. Given the N mostsignificant peaks (in terms of magnitude) this solves for directional and ambient intensities in aleast squares sense.
The equations to compute the light intensities for two lobes are:
𝐴𝐴
𝑐𝑐0
2
2
�𝑐𝑐 � = �𝐴𝐴 − 𝐵𝐵−𝐵𝐵
1
𝐴𝐴2 − 𝐵𝐵2

Where:−𝐵𝐵

𝐴𝐴2 − 𝐵𝐵2𝐴𝐴

𝐴𝐴2 − 𝐵𝐵2

𝐿𝐿 °𝐿𝐿� � 𝑒𝑒 𝑑𝑑0 �
𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑0

𝐴𝐴 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑0 , 𝐵𝐵 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑1

The ambient term is then:𝑎𝑎 = �𝐿𝐿𝑒𝑒 [0] − (𝑐𝑐0 𝐿𝐿𝑑𝑑0 [0] + 𝑐𝑐1 𝐿𝐿𝑑𝑑1 [0])�

For three lights the intensities are:

Where:

𝐴𝐴2 − 𝐷𝐷 2
⎡
𝐸𝐸
⎢
𝑐𝑐0
⎢
𝐶𝐶𝐶𝐶
−
𝐴𝐴𝐴𝐴
�𝑐𝑐1 � = ⎢
𝐸𝐸
𝑐𝑐2
⎢
𝐵𝐵𝐵𝐵
− 𝐴𝐴𝐴𝐴⎢⎣
𝐸𝐸

𝐶𝐶𝐶𝐶 − 𝐴𝐴𝐴𝐴
𝐸𝐸
𝐴𝐴2 − 𝐶𝐶 2
𝐸𝐸
𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴
𝐸𝐸

√𝜋𝜋
2

𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴⎤𝐸𝐸
⎥ 𝐿𝐿 °𝐿𝐿𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴 ⎥ 𝑒𝑒 𝑑𝑑0⎥ �𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑1 �
𝐸𝐸⎥ 𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑2
𝐴𝐴2 − 𝐵𝐵2 ⎥⎦𝐸𝐸

𝐶𝐶 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑2 , 𝐷𝐷 = 𝐿𝐿𝑑𝑑1 ∘ 𝐿𝐿𝑑𝑑2, 𝐸𝐸 = 2𝐵𝐵𝐵𝐵𝐵𝐵 + 𝐴𝐴(𝐴𝐴2 − 𝐵𝐵2 −𝐶𝐶 2 − 𝐷𝐷 2 )This matrix is symmetric, and the ambient coefficient is:𝑎𝑎 = �𝐿𝐿𝑒𝑒 [0] − (𝑐𝑐0 𝐿𝐿𝑑𝑑0 [0] + 𝑐𝑐1 𝐿𝐿𝑑𝑑1 [0] + 𝑐𝑐2 𝐿𝐿𝑑𝑑2 [0])�√𝜋𝜋2

The derivations of the equations are in Appendix A5 Solving for Coefficients to Approximate SH
Environment Map with Directional and Ambient Lights.
Adding all of the degrees of freedom (directions and intensities) and using a non-linear solver
would generate higher quality results, possibly using the technique above as an initial guess.

Ringing
Ringing, also called Gibbs Phenomenon, is a common problem in signal processing. When asignal with a discontinuity is projected into a finite Fourier basis (which can only represent
continuous functions) overshoot and undershoot will happen around the discontinuity.
Functions that don’t have discontinuities can exhibit similar behavior if the projection istruncated. We have already seen these problems when looking at lighting models, and in
representing irradiance environment maps (projection of clamped cosine function.) A similar
problem occurs in surface design, where when trying to satisfy a set of geometric constraints
unwanted oscillations can occur. There are two general solutions to these problems:
1)

2)

Windowing the truncated projection coefficients using sigma factors. This is the mostcommon solution in signal processing, and can be trivially used with Spherical
Harmonics [8][41][37].
Minimizing some form of variational function (minimizing a measure of curvature forexample), instead of just the standard least-squares error. This is commonly done incomputer aided geometric design, but also can be efficiently done using Spherical
Harmonics [38].

Windowing 8
One way to minimize the ringing artifacts is to multiply in the frequency domain (which is a
convolution in the spatial domain) by a kernel with projection coefficients that taper to zero as
you approach the cut-off frequency. If this function is a sinc 9 that is stretched out so that itreaches zero at the truncated frequency band, it is called using Lanczos sigma factors 10.
Intuitively what this is doing (in 1D) is convolving in the spatial domain with a tight box function,
making the function smooth enough to be
represented without excessive ringing. There are
more sophisticated ways to attack Gibbs phenomena[10][4] but they use the SH coefficients to generate
piecewise analytic functions that would not be as
convenient for games.
In our experience the choice of windowing function is
not as important as having flexibility to trade off
between ringing and blurring. The image to the left
shows the two windowing functions (red is sinc, blue
is raised cosine lobe – called a Hanning window)
8

The more common use of the term “windowing” in signal processing is used in the spatial domain –when filtering images the spatial extent of the filter is “windowed”, and when taking a FFT of a 1D signal
the signal might be scaled to make it periodic. In this paper it is being done in the frequency domain.
sin 𝑥𝑥9sinc(𝑥𝑥) =where the limit is 1.0 when x=0.
𝑥𝑥
10Sometimes the sigma factors will be raised to a power to more aggressively reduce the Gibbs
phenomena – which is the same as repeatedly convolving the signal.

scaled for 6th order SH (so they reach zero for the 7th order band, the last value used would be
evaluated at 5 and the functions would only be evaluated at the integers.) The Hanning function
decays faster than the Lancosz function, which makes it blur more aggressively 11.
Equation

W=6

W=10

𝜋𝜋𝜋𝜋
sin 𝑤𝑤
𝜋𝜋𝜋𝜋
𝑤𝑤

𝜋𝜋𝜋𝜋
�1 + cos 𝑤𝑤 �
2

11

The Lancosz squared is close to the Hanning function, but decays a bit faster. For the reconstruction of
the delta function to be a height of 2.25, the Hanning function needs a window of size 12.0105 and theLancosz needs one of size 9.8725. The results look visually indistinguishable.

Above we show the results using Lanczos and Hanning sigma factors. The signal being projected
is an order 6 delta function, which is the “peakiest” signal you can project into SH and exhibits
ringing artifacts. The projection of a delta function is a ZH, so we are showing a cross section of
the sphere, phi is fixed. Radial magnitude is plotted, the sign of the lobes alternate.
Looking at all the graphs together (red is the raw projection of the delta function) you can see
how windowing blurs the signal while eliminating the rings (visible near the origin in the figures.)
Hanning

Lancosz

Minimizing a Functional
An alternative approach is to try and minimize some function besides squared approximation
error. One way to do this is to satisfy a set of constraints (for example exact reconstruction at asmall number of points) and then use the “slack” variables left over (assuming there are enoughdegrees of freedom) to minimize some error functional[38]. Given the low SH order often usedin games/graphics, this approach does not seem that practical so I will not spend any more timeon it. An alternative is to attempt to minimize a norm that penalizes large oscillations. This canbe done with Spherical Harmonics in a straightforward fashion. The Laplace operator, or
Laplacian, is the divergence of the gradient of a scalar function; equivalently it is the sum of
unmixed partial derivatives

Δ𝑓𝑓 =

𝑑𝑑𝑑𝑑 2 𝑑𝑑𝑑𝑑 2 𝑑𝑑𝑑𝑑 2++
𝑑𝑑𝑑𝑑 2 𝑑𝑑𝑑𝑑 2 𝑑𝑑𝑑𝑑 2

In spherical coordinates on the unit sphere this is
Δ𝑓𝑓 =

1 𝑑𝑑
𝑑𝑑𝑑𝑑
1 𝑑𝑑2 𝑓𝑓
�sin 𝜃𝜃 � + 2
sin 𝜃𝜃 𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑
sin 𝜃𝜃 𝑑𝑑𝜙𝜙 2

The integral of the squared Laplacian is a curvature measure used on the sphere [38]. The
function we are going to minimize is:2

2

𝐸𝐸(𝑐𝑐) = � �𝑓𝑓̃(𝑠𝑠) − 𝑓𝑓(𝑠𝑠)� 𝑑𝑑𝑑𝑑 + 𝜆𝜆 � �Δ𝑓𝑓̃(𝑠𝑠)� 𝑑𝑑𝑑𝑑

We will tweak this a little bit, we already know the raw projection coefficients 𝑓𝑓𝑙𝑙𝑚𝑚 , so we wantto find new coefficients 𝑐𝑐𝑙𝑙𝑚𝑚 that are as close to the least squares result as possible, while also
minimizing the weighted squared Laplacian. This can be done in closed form [see Appendix A6
Least Squares Projection], which given a 𝜆𝜆 results in the following coefficients.
𝑐𝑐𝑙𝑙𝑚𝑚 =𝑓𝑓𝑙𝑙𝑚𝑚(1 + 𝜆𝜆𝑙𝑙 2 (𝑙𝑙 + 1)2 )Note that this amounts to a windowing function that is dependent on 𝜆𝜆 . When 𝜆𝜆 is zero, you
get the least squares coefficients, when 𝜆𝜆 is infinity you get just the DC term, which has zerocurvature. One approach to choosing 𝜆𝜆 is to solve that it reduces the squared Lalpacian by afixed amount, say one half. This can be done using any standard root finding technique. See[Appendix A7 Solving for Lamba to Reduce the Squared Laplacian] for an explanation of how to
do this using Newtons method. Below are the results with a 6ht order delta function solving forthe squared Laplacian to be 10% (𝜆𝜆 = .004209 green) and 50% (𝜆𝜆 = .000632 blue) of theoriginal. The final plot has them along with the delta function itself:
Δ2 =0.1

Δ2 =0.5

Composite

Here are images using an actual lighting environment. The first column of images are contour
plots (blue is negative) of the second column of images - windowing with a broad (smoother)window. The forth column of images are contour plots of the third column – smoothing with anarrower (less blurry) filter. The top row is the original image – all results are shown at 6thorder:
Original

Cos

Sinc

Δ2

Windowing should be used judiciously. When using irradiance environment maps, theconvolution with the clamped cosine function aggressively attenuates high frequencies andwindowing is rarely necessary. Also scenes where there is a lot of normal variation tend to notshow off the ringing artifacts as much as ones where the normals of the receiver are smoother.
Below are images of a simple scene (“door” and a ground plane) that show how windowing
affects the shaded result. An area where ringing occurs is highlighted in the second row.

No Windowing

Hanning 6

Hanning 10

When using HDR, or just moderately bright lights, ringing is a more severe problem and evenirradiance environment maps might need some amount of windowing
Another problem that can happen with ringing is color artifacts. In the images below we have asphere illuminated with a bright directional yellow light (almost rim lit from the upper right,
order 6) and a moderate ambient white light. The top row uses irradiance environment maps,
the bottom row uses order 6 PRT (so a more accurate approximation of the cosine kernel isused.) The first column is without windowing, the second column uses a Hanning window oforder 4 on the top row and order 5 on the bottom. The unwindowed version shows both thepositive lobe (order 3, negative order 6) and the blue bands are from the negative rings of thedirectional light (remove red/green from the ambient light, but leave the blue alone.)

Order 3

Order 6

No Windowing

Hanning

Content Sensitive Windowing
Lighting can be globally windowed, but you can also window depending on the influence therings will have on the final shaded image. An example is illustrated below. We are rendering asomewhat matte (but not diffuse – phong with power 10) ball using a bright directional lightsource (the highlight is saturating.) If no windowing is applied, the principle highlight is sharper,
but artifacts from the rings are clear. If the lighting is windowed, the rings disappear but the
principle highlight is blurred as well. Instead you can look at the angle between the reflection
vector and the dominant light direction, if this is small, no windowing is necessary, as itincreases you can blend between the windowed and non-windowed light source, preserving the
sharp highlight while eliminating the ringing artifacts. This sequence of images is shown below:

6th order reflection

Windowed (Hanning 7)

Content sensitive windowing

The equation used to control the blending is:�(𝑛𝑛 ∘ 𝑟𝑟) − 𝑐𝑐𝑡𝑡 �
𝑤𝑤𝑎𝑎 = �𝑚𝑚𝑚𝑚𝑚𝑚 �0,��
1 − 𝑐𝑐𝑡𝑡𝑝𝑝

Where 𝑤𝑤𝑎𝑎 is the weight of the un-windowed light source (windowed is 1 − 𝑤𝑤𝑎𝑎 ), 𝑐𝑐𝑡𝑡 is a thresholdthat determines when you completely use the windowed light source and 𝑝𝑝 controls thetransition region of the blend. In this figure 𝑐𝑐𝑡𝑡 is 0.07 and 𝑝𝑝 is 0.8. You can experiment withthese parameters, the thresholds will largely be based on the amount of windowing and thematerial properties.
If the light source is not a directional one, you can compute how well it can be approximated bya directional light source and take that into account when deciding how much to blend (thesimplest way would be to compute a tensor product between how close the lighting
environment is approximated by a directional light and how close the normal/reflection vector isto the dominant light vector.) For more complex shading, say for example PRT, the dominanttransfer direction can be used, and how well it approximates the transfer function can be yet
another factor in the tensor product (so you window when any one of the terms isn’t
approximated well.)
If the lighting is very dynamic, you will need to watch out for temporal artifacts that could occur
as the lights change (ie: shadows/reflections would become sharper, etc.) For things like static
SH light probes this technique should work well though.

SH Products
Computing the SH representation of the product of two functions represented using SH is often
useful. Example scenarios are:3)

Punch a hole in a skylight model based on a large flying object (visibility times light), or asimple visibility model of the scene (large buildings, etc.)

4)5)Multiply visibility functions. This happens when doing dynamic approximate global
illumination.
Scale or modify a SH light probe. Multiplying by a some constant between zero and 1
can be used to approximate clouds for example.

Computing a product in the frequency domain is quite complex, it boils down to applying a“triple product tensor” times the two SH vectors. Code for this can be generated efficiently [47]and will not be described in this paper. There are several special cases that are worth
mentioning though.

Products with a Constant Function
If one of the SH functions is going to be used a lot you can build a dense matrix called the
product matrix, this makes the triple product a simple matrix times vector product, which issignificantly less expensive. An order 6 product would have 1296 multiplies instead of the 2527
in the code generated by [47].

Products with varying Orders
This is particularly common when the output order is lower, for example quadratic, so that you
can represent a local radiance environment. Special casing the code for these cases cansignificantly decrease the complexity of the code. For example the product of two 6th order SH
has 2527/1995 multiplies/adds when computing 6th order results, but only 933/676 12 when
computing 3rd order results. Another example is simple ambient occlusion, in this case one of
the terms is just DC, and so you simply have to scale the other vector by the DC value. Finally,
one of the two functions might be lower order (ie: just multiply by linear visibility) which canalso reduce the cost.

Product with Zonal Harmonics
If one of the functions is a Zonal Harmonic, you can rotate the other function into that same
frame (less expensive rotation, only need 2 Euler angles due to the symmetry), compute theproduct and rotate back. The sparsity in the ZH frame eliminates a significant amount of work,
which can increase performance. The product of two 6th order functions, one which is a ZH
oriented in z, requires only 380/249 multiplies/adds. For arbitrary ZH the time to compute one
million products (where one of the functions is always a in an arbitrary orientation) is around 1.2
seconds while the time to compute general 6th order products is just over 3 seconds, so thistechnique is almost 3 times faster.

Product with Analytic Function
If one of the functions has an analytic form, it is more accurate to analytically compute whatamounts to the product matrix. An example is zeroing everything under the horizon (useful if
there is a ground plane) or taking a product with a clamped cosine function. Doing this12

More efficient code could be generated, these results are using a naïve algorithm applied to the output
of [47].

analytically is equivalent to having an infinite order expansion of the analytic function, and ingeneral will be much faster than using a SH expansion in these cases. The code for these two
examples is in an appendix [Appendix A8 Code for Multiplying SH by Analytic Functions].

Conclusions
Spherical Harmonics are extremely useful tool, particularly for lighting, in games. Hopefully thisarticle shed some light on how they can be used and how to mitigate some challenges that canarrive when using them. The are several ways to extend the ideas discussed in this paper,
windowing coefficients (using any window function) can be solved for that minimize themagnitude of a negative (or positive) lobe, or possibly the magnitude of reflected radiance whenpointing in a direction that should be zero. When extracting lights, more rigorous techniques
based on non-linear optimization could be used [29][44] and more general lighting models (forexample including the cone angles from the light sources discussed) could be extracted. It also
might be worth investigating tying windowing into the fitting process. Initially fit with asmoothed version of the function, and then dial back the amount of windowing – effectively
steering to better local minima. The content sensitive windowing ideas need to be fleshed out,
particularly when integrating with techniques like PRT.
One tool that I’ve found to be invaluable when playing around with Spherical Harmonics is asymbolic math program. I’ve used Maple (http://www.maplesoft.com), but other programs like
Mathematica (http://www.wolfram.com) would work equally well. The DirectX SDK(http://msdn.microsoft.com/directx) has functions for evaluation, rotation, products and several
analytic lighting models, along with samples using both PRT and irradiance environment maps.

Acknowledgements
Much of the work has been done with other collaborators, in particular a long and fruitful seriesof papers with John Snyder, starting when we were colleagues in the graphics group at
Microsoft Research. I’ve also had the pleasure of discussing Spherical Harmonics with several
game developers: Hao Chen (Bungie), Arn Arndt , James Grieves, Clint Hanson (EA at the time),
Loren McQuade (while at Blizzard), Dan Baker (MS and Firaxis), Alex Evans (Lionhead at thetime), Tom Forsyth (Muckyfoot), Willem de Boer (Muckyfoot), Manchor Ko (Naughty Dog), Naty
Hoffman (SCEA) and Chris Oat (ATI). Jason Sandlin, Ben Luna and Jon Steed worked on thesamples/test code at Microsoft. Finally I’ve enjoyed participating in discussions on the
GDAlgorithms mailing list on SH and other topics.

References
[1]
[2]
[3]

ANEN, T, KAUTZ, J, DURAND, F, SEIDEL, H. Spherical Harmonic Gradients for Mid-Range Illumination,
Eurographics Symposium on Rendering 2004.
AXLER, S, BOURDON, P, RAMEY, W. Harmonic Function Theory, Springer-Verlag, 2000.
BASRI, R, AND JACOBS, D. Lambertian Reflectance and Linear Subspaces, ICCV 2001.

[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]

BLAKELY, C, GELB, A, NAVARRA, A. An Automated Method for Recovering Piecewise Smooth Functions on
Spheres Free from Gibbs Oscillations, Sampling Theory in Signal and Image Processing Vol 6, No. 3, Sep 2007.
BYERLY, W. An Elementary Treatise on Fourier’s Series and Spherical, Cylindrical and Ellipsoidal Harmonics,
with Applications to Problems in Mathematical Phyisics, Dover, 1893.
CABRAL, B, MAX, N, SPRINGMEYER, R. Bidirectional Reflection Functions from Surface Bump Maps, SIGGRAPH
1987.
CHANDRASEKHAR, S. Radiative Transfer, Dover, 1960.
DOBASHI, Y, KANEDA, K, NAKATANI, H, YAMASHITA, H. A Quick Rendering Method Using Basis Functions forInteractive Lighting Design, EUROGRAPHICS 1995.
FREEDEN, W, GERVENS, T, SCHREINER, M. Construction Approximation on the Sphere, Clarendon Press, 1998.
GELB, A. The Resolution of the Gibbs Phenomenon for Spherical Harmonics, Mathematics of Computation,
Volume 66, Number 218, April 1997.
GENEVAUX, O, LARUE, F, DISCHLER, J. Interactive Refraction on Complex Static Geometry using Spherical
Harmonics, Symposium on Interactive 3D Graphics and Games 2006.
GREEN, P, KAUTZ, J, DURAND, F. Efficient Reflectance and Visibility Approximation for Environment Map
Rendering, Eurographics 2007.
GREEN, R. Spherical Harmonic Lighting: The Gritty Details, GDC 2003.
HABEL, R, MUSTATA, B, WIMMER, M. Efficient Spherical Harmonics Lighting with the Preetham Skylight
Model. Eurographics 2008.
HAN, C, SUN, B, RAMAMOORTHI, R, GRINSPUN, E. Frequency Domain Normal Map Filtering, SIGGRAPH 2007.
INUI, T, TANABE, Y, ONODERA, Y. Group Theory and Its Applications in Physics, Springer-Verlag, 1996.
ISHIMARU, A. Wave Propagation and Scattering in Random Media, IEEE Press, 1978.
KAJIYA, J, von HERZEN, B. Ray Tracing Volume Densities, SIGGRAPH 1984.
KAUTZ, J, SLOAN, P, SNYDER, J. Fast, Arbitrary BRDF Shading for Low-Frequency Lighting Using Spherical
Harmonics, EUROGRAPHICS Workshop On Rendering, 2002.
KAUTZ, J, LEHTINEN, J, SLOAN, P. Precomputed Radiance Transfer: Theory and Practice. SIGGRAPH 2005
course. http://www.cs.ucl.ac.uk/staff/j.kautz/PRTCourse/
KO, J, KO, M, ZWICKER, M. Practical Methods for a PRT-based Shader Using Spherical Harmonics, Shader X6:
Advanced Rendering Techniques, 2008.
LANN, P, LEUNG, C, WONG, T. Noise-Resistant Fitting for Spherical Harmonics, IEEE Transactions on
Visualization and Computer Graphics. 12(2) 2006.
LEHTINEN, J, KAUTZ, J. Matrix Radiance Transfer, Symposium on Interactive 3D Graphics, 2003.
LEHTINEN, J. A Framework for Precomputed and Captured Light Transport. ACM Transactions on Graphics
26(4) 2007.
MACROBERTS, T. Spherical Harmonics, Dover, 1948.
MCTAGGART, G. Half-Life 2 Source Shading, GDC 2004.
NG, R, RAMAMOORTHI, R, HANRAHAN, P. All-Frequency Shadows Using Non-Linear Wavelet Lighting
Approximation, SIGGRAPH 2003.
NIMEROFF,J, SIMONCELLI, E, DORSEY, J. Efficient Re-rendering of Natural Environments, EUGROGRAPHICS
Workshop on Rendering, 1994.
NOCEDAL, J, WRIGHT, S. Numerical Optimization, Springer-Verlag 1999.
OAT, C. Adding Spherical Harmonic Lighting to the Sushi Engine. GDC 2004.
http://ati.amd.com/developer/gdc/Oat-GDC04-SphericalHarmonicLighting.pdf
OAT, C. Irradiance Volumes for Games. GDC 2005.
http://ati.amd.com/developer/gdc/GDC2005_PracticalPRT.pdf
OAT, C. Real-Time Irradiance Volumes, ShaderX 5: Advanced Rendering Techniques, 2006.
PAN, M, WANG, R, LIU, X, PENG, Q, HUJUN, B. Precomputed Radiance Transfer Field for Rendering
Interreflections in Dynamic Scenes, Eurographics 2007.
PREETHAM, A, SHIRLEY, P, SMITS, B. A Practical Analytic Model for Daylight, SIGGRAPH 1999.
RAMAMOORTHI, R, AND HANRAHAN, H. An Efficient Representation for Irradiance Environment Maps,
SIGGRAPH 2001.

[36]
[37]
[38]
[39]
[40]
[41]
[42]
[43]
[44]
[45]
[46]
[47]
[48]
[49]
[50]
[51]
[52]
[53]
[54]
[55]
[56]

RAMAMOORTHI, R, AND HANRAHAN, H. Frequency Space Environment Map Rendering, SIGGRAPH 2002.
REN, Z, WANG, R, SNYDER, J, ZHOU, K, SUN, B, SLOAN, P, BAO, H, PENG, Q, GUO, B. Real-Time Soft Shadows in
Dynamic Scenes Using Spherical Harmonic Exponentiation, SIGGRAPH 2006.
RUFFINI, G, MARCO, J, GRAU, C. Spherical Harmonics Interpolation, Computation of Laplacians and Gauge
Theory, arXiv:physics/0206007, 2002.
SCHROEDER, P, SWELDENS, W. Spherical Wavelets: Efficiently Representing Functions on the Sphere,
SIGGRAPH 1995.
SILLION, F, ARVO, J, WESTIN, S, GREENBERG, D. A Global Illumination Solution for General ReflectanceDistributions, SIGGRAPH 1991.
SLOAN, P, KAUTZ, J, SNYDER, J. Precomputed Radiance Transfer for Real-Time Rendering in Dynamic, LowFrequency Lighting Environments, SIGGRAPH 2002.
SLOAN, P, HALL, J, HART, J, SNYDER, J. Clustered Principal Components for Precomputed Radiance Transfer,
SIGGRAPH 2003.
SLOAN, P, LIU, X, SHUM, H-Y, SNYDER, J. Bi-Scale Radiance Transfer, SIGGRAPH 2003.
SLOAN, P, LUNA, B, SNYDER J, Local, Deformable Precomputed Radiance Transfer, SIGGRAPH 2005.
SLOAN, P. Normal Mapping for Precomputed Radiance Transfer, Symposium on Interactive 3D Graphics andGames 2006.
SLOAN, P, GOVINDARAJU, N, NOWROUZEZAHRAI, D, SNYDER, J. Image-Based Proxy Accumulation for RealTime Soft Global Illumination, Pacific Graphics 2007.
SNYDER, J. Code Generation and Factoring for Fast Evaluation of Low-Order Spherical Harmonic Products andSquares, Microsoft Research Technical Report, MSR-TR-2006-53, 2006.
SNYDER, J. Personal Communication, 2008.
SUN, B, RAMAMOORTHI, R, NARASIMHAN, S, NAYAR, S. A Practical Analytic Single Scattering Model for RealTime Rendering, SIGGRAPH 2005.
STAM, J. Multiple Scattering as a Diffusion Process, EUROGRAPHICS Workshop on Rendering, 1995.
VARSHALOVICH, D, MOSKALEV, A, KHERSONSKII, V. Quantum Theory of Angular Momentum, World Scientific
Publishing, 1988.
WANG, J, XU, K, ZHOU, K, LIN, S, HU, S, GUO, B. Spherical Harmonic Scaling, The Visual Computer, 22(9-11),
2006.
WESTIN, S, ARVO, J, TORRANCE, K. Predicting Reflectance Functions from Complex Surfaces”, SIGGRAPH 1992.
WONG, T, HENG, P, OR, S, NG, W. Image-Based Rendering with Controllable Illumination, Eurographics
Workshop on Rendering, 1997.
WONT, T, FU, C, HENG, P, LEUNG, C. The Plenoptic Illumination Function, IEEE Transactions on Multimedia 4(3)2002.
ZHOU, K, HU, Y, LIN, S, GUO, B, SHUM, H-Y. Precomputed Shadow Fields for Dynamic Scenes, SIGGRAPH 2005.

Appendix A1 Recursive Rules for Evaluating SH Basis Functions
Recurrence relations [48] can be used to efficiently evaluate the polynomial forms of the SHbasis functions. Recall the formulas for the basis functions:√2𝐾𝐾𝑙𝑙𝑚𝑚 cos 𝑚𝑚𝑚𝑚 𝑃𝑃𝑙𝑙𝑚𝑚 (cos 𝜃𝜃) 𝑚𝑚 > 0√2Re(𝑌𝑌𝑙𝑙𝑚𝑚 ) 𝑚𝑚 > 0𝑦𝑦𝑙𝑙𝑚𝑚 = �√2Im(𝑌𝑌𝑙𝑙𝑚𝑚 ) 𝑚𝑚 < 0 = �√2𝐾𝐾𝑙𝑙𝑚𝑚 sin|𝑚𝑚|𝜑𝜑 𝑃𝑃𝑙𝑙|𝑚𝑚 | (cos 𝜃𝜃) 𝑚𝑚 < 0
𝑚𝑚 = 0𝑚𝑚 = 0
𝑌𝑌𝑙𝑙0
𝐾𝐾𝑙𝑙0 𝑃𝑃𝑙𝑙0 (cos 𝜃𝜃)Given a point (x,y,z) on the unit sphere, the Associated Legendre polynomials (only depend on Z,
divided by sin 𝜃𝜃 𝑚𝑚 ) can be evaluated using these recurrences (𝑃𝑃00 = 1):

𝑃𝑃𝑚𝑚𝑚𝑚
𝑃𝑃𝑚𝑚𝑚𝑚+1
𝑃𝑃𝑙𝑙𝑚𝑚−1(1 − 2𝑚𝑚)𝑃𝑃𝑚𝑚𝑚𝑚−1(2𝑚𝑚 + 1)𝑧𝑧𝑃𝑃𝑚𝑚𝑚𝑚
𝑚𝑚
𝑚𝑚(2𝑙𝑙 − 1)𝑧𝑧𝑃𝑃𝑙𝑙−1− (𝑙𝑙 + 𝑚𝑚 − 1)𝑃𝑃𝑙𝑙−2=
𝑙𝑙 − 𝑚𝑚==

Where you increment l in the inner loop, m in the outer.

The trigonometric addition formula can be used to evaluate the 𝜑𝜑 dependent terms (multipliedby sin 𝜃𝜃 𝑚𝑚 so that they can be computed from (x,y) coordinates on the sphere), where𝑆𝑆(0) = 0, 𝐶𝐶(0) = 1:
𝑆𝑆(𝑚𝑚) = 𝑥𝑥𝑥𝑥(𝑚𝑚 − 1) + 𝑦𝑦𝑦𝑦(𝑚𝑚 − 1)
𝐶𝐶(𝑚𝑚) = 𝑥𝑥𝑥𝑥(𝑚𝑚 − 1) − 𝑦𝑦𝑦𝑦(𝑚𝑚 − 1)

In the original equation replace sin|𝑚𝑚|𝜑𝜑 with 𝑆𝑆(𝑚𝑚) and cos 𝑚𝑚𝑚𝑚 with 𝐶𝐶(𝑚𝑚). It is generally
more efficient to generate code based on these equations, since it naturally factors common
subexpressions.

Appendix A2 Polynomial Forms of SH Basis
The polynomial forms of the SH basis functions are listed below, L is the band index, M is thebasis function. Note that Maple randomly changes the order of L and M…

{ L = 0, M = 0 },1
2 π

{ L = 1, M = -1 }, −{ L = 1, M = 0 },

3y
2 π

3z
2 π

{ L = 1, M = 1 }, −3x2 π

{ L = 2, M = -2 },15 y x2 π{ M = -1, L = 2 }, −{ M = 0, L = 2 },15 y z2 π

5 ( 3 z2 − 1 )4 π

{ M = 1, L = 2 }, −15 x z2 π{ M = 2, L = 2 },

15 ( x2 − y2 )4 π

{ M = -3, L = 3 }, −105 y x z
2 π

{ M = -2, L = 3 },{ M = -1, L = 3 }, −2

{ M = 1, L = 3 }, −21 x ( −1 + 5 z2 )
8 π

105 ( x2 − y2 ) z
4 π{ L = 3, M = 2 },{ M = 3, L = 3 }, −{ L = 4, M = -4 },2

35 x ( x2 − 3 y2 )8 π

3 35 y x ( x2 − y2 )4 π

{ M = -3, L = 4 }, −3 235 y ( 3 x2 − y2 ) z
8 π

3 5 y x ( −1 + 7 z 2 )4 π{ M = -1, L = 4 }, −{ L = 4, M = 0 },21 y ( −1 + 5 z2 )8 π

27 z ( 5 z2 − 3 )4 π{ L = 3, M = 0 },{ M = -2, L = 4 },35 y ( 3 x2 − y2 )
8 π

2

3 2

5 y z ( −3 + 7 z 2 )8 π

3 ( 35 z4 − 30 z2 + 3 )16 π

{ L = 4, M = 2 },

3 5 ( x2 − y2 ) ( −1 + 7 z 2 )8 π

3 35 ( x4 − 6 y2 x2 + y4 )16 π

{ L = 5, M = -5 }, −{ L = 5, M = -4 },35 x ( x2 − 3 y2 ) z
8 π

3 2

{ M = 3, L = 4 }, −{ M = 4, L = 4 },5 x z ( −3 + 7 z 2 )
8 π

3 2

{ M = 1, L = 4 }, −3 2

77 y ( 5 x4 − 10 y2 x2 + y4 )32 π

3 385 y x ( x2 − y2 ) z
4 π{ L = 5, M = -3 }, −2385 y ( 3 x2 − y2 ) ( −1 + 9 z 2 )32 π1155 y x z ( −1 + 3 z 2 )4 π{ M = -2, L = 5 },{ M = -1, L = 5 }, −165 y ( −14 z2 + 21 z4 + 1 )16 π

11 z ( 63 z 4 − 70 z2 + 15 )16 π

{ L = 5, M = 0 },{ L = 5, M = 1 }, −165 x ( −14 z2 + 21 z4 + 1 )16 π{ L = 5, M = 2 },1155 ( x2 − y2 ) z ( −1 + 3 z 2 )8 π{ L = 5, M = 3 }, −{ L = 5, M = 4 },2385 x ( x2 − 3 y2 ) ( −1 + 9 z 2 )32 π

3 385 ( x4 − 6 y2 x2 + y4 ) z
16 π

{ L = 5, M = 5 }, −3 2

77 x ( x4 − 10 y2 x2 + 5 y4 )32 πAppendix A3 ZH Coefficients for Spherical Light Source
Given a light source that subtends an angle a in radians here are the symbolic integrals for the
first 6 bands:
L=0 − π ( −1 + cos( a ) )1
3
2

π sin( a ) 2

L=2 −

1
5
2

π cos( a ) ( −1 + cos( a ) ) ( cos( a ) + 1 )

L=3 −

1
7
8

π ( −1 + cos( a ) ) ( cos( a ) + 1 ) ( 5 cos( a ) 2 − 1 )

L=4 −

3
π cos( a ) ( −1 + cos( a ) ) ( cos( a ) + 1 ) ( 7 cos( a ) 2 − 3 )
8

L=5 −

1
11
16

L=1

π ( −1 + cos( a ) ) ( cos( a ) + 1 ) ( 21 cos( a ) 4 − 14 cos( a ) 2 + 1 )Appendix A4 ZH Coefficients for Smooth Cone
Given a cone that subtends the angle a in radians, the light source has intensity 1 at the northpole, and falls off to zero at angle a. At 6th order this function should not be evaluated with anangle less than about 8 degrees using single precision. The derivative of the smoothing function
is 0 at the north pole and a. The first 6 bands are:( a 3 + 6 a − 12 sin( a ) + 6 cos( a ) a ) π
a3

1 3 ( a 3 − 3 cos( a ) sin( a ) + 3 cos( a ) 2 a ) π
4
a3

1 5 ( −6 a − 2 cos( a ) 2 sin( a ) − 9 cos( a ) a + 14 sin( a ) + 3 cos( a ) 3 a ) π
9
a3

1
7(
256
4 a 3 + 15 a − 108 cos( a ) 2 a − 30 cos( a ) 3 sin( a ) + 63 cos( a ) sin( a ) + 60 cos( a ) 4 a )πa31( −480 a + 742 sin( a ) + 596 cos( a ) 2 sin( a ) + 225 cos( a ) a − 378 sin( a ) cos( a ) 4
1500− 1650 cos( a ) 3 a + 945 cos( a ) 5 a ) π

a3

1
11 ( −63 a + 12 a 3 + 350 cos( a ) 3 sin( a ) − 1260 cos( a ) 4 a − 15 cos( a ) sin( a )3072− 224 cos( a ) 5 sin( a ) + 672 cos( a ) 6 a + 540 cos( a ) 2 a ) π

a3Appendix A5 Solving for Coefficients to Approximate SH
Environment Map with Directional and Ambient Lights
You can compute the intensity s that minimizes the approximation error between a directional
light in direction d and the original lighting environment:
𝐸𝐸(𝑐𝑐) = (𝐿𝐿𝑒𝑒 − 𝑐𝑐𝐿𝐿𝑑𝑑 )2

Where 𝐿𝐿𝑒𝑒 is the SH representation of the lighting environment, 𝐿𝐿𝑑𝑑 is the SH representation of
the lighting model in direction d 13. The solution is:
𝑐𝑐 =𝐿𝐿𝑒𝑒 ∘ 𝐿𝐿𝑑𝑑
𝐿𝐿𝑑𝑑 ∘ 𝐿𝐿𝑑𝑑

If you want to add an ambient light, you need to minimize the error function:
2

𝐸𝐸(𝑐𝑐, 𝑎𝑎) = ��𝑐𝑐𝐿𝐿𝑑𝑑 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠) + 𝑎𝑎𝐿𝐿𝑎𝑎 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠) − 𝐿𝐿𝑒𝑒 (𝑠𝑠) ∗ 𝐻𝐻𝑁𝑁 (𝑠𝑠)� 𝑑𝑑𝑑𝑑

Absorbing the convolution into the lighting and differentiate with respect to each variable:

13

You can optimize using either radiance (leave vectors alone) or irradiance (convolve vectors.)𝑑𝑑𝑑𝑑
= 2 � �𝑐𝑐𝐿𝐿�𝑑𝑑 (𝑠𝑠) + 𝑎𝑎𝐿𝐿�𝑎𝑎 (𝑠𝑠) − 𝐿𝐿�𝑒𝑒 (𝑠𝑠)� 𝐿𝐿�𝑑𝑑 (𝑠𝑠)𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑
= 2 � �𝑐𝑐𝐿𝐿�𝑑𝑑 (𝑠𝑠) + 𝑎𝑎𝐿𝐿�𝑎𝑎 (𝑠𝑠) − 𝐿𝐿�𝑒𝑒 (𝑠𝑠)� 𝐿𝐿�𝑎𝑎 (𝑠𝑠)𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑

Now solve for the two equations equal to zero 14 to find the minimum. 𝐿𝐿�𝑎𝑎 (𝑠𝑠) is a scalarfunction 15, and due to the orthogonality of SH the integrals are all simple dot products of SH
vectors. This leaves you with the following equations:𝑐𝑐𝑐𝑐 + 𝑎𝑎𝑎𝑎 = 𝐷𝐷
𝑐𝑐𝑐𝑐 + 𝑎𝑎𝑎𝑎 = 𝐸𝐸

Where:
𝐴𝐴 =

16
4
508𝜋𝜋
, 𝐵𝐵 =, 𝐶𝐶 = , 𝐷𝐷 = 𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑑𝑑 , 𝐿𝐿�𝑒𝑒 �, 𝐸𝐸 = 𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑎𝑎 , 𝐿𝐿�𝑒𝑒 �
17
𝜋𝜋
867

Solving for c and a results in:

51
867
𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑑𝑑 , 𝐿𝐿�𝑒𝑒 � −𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑎𝑎 , 𝐿𝐿�𝑒𝑒 �
79
316𝜋𝜋
51
127𝜋𝜋
𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑎𝑎 , 𝐿𝐿�𝑒𝑒 � −𝑑𝑑𝑑𝑑𝑑𝑑�𝐿𝐿�𝑑𝑑 , 𝐿𝐿�𝑒𝑒 �
𝑎𝑎 =
79
316
𝑐𝑐 =

It turns out there is a simpler way to arrive at the same results. Solve for the intensity of adirectional light where both the environment and the light don’t include the DC term, thencomputing an ambient term that reconstructs the environment DC term when using the scaled
directional light. This results in:
867
𝑑𝑑𝑑𝑑𝑑𝑑0 �𝐿𝐿�𝑑𝑑 , 𝐿𝐿�𝑒𝑒 �316𝜋𝜋
8√𝜋𝜋 √𝜋𝜋𝑎𝑎 = �𝐿𝐿�𝑒𝑒 [0] − 𝑐𝑐
�
17 2
𝑐𝑐 =

In the above expression the dot product ignores the DC term and 𝐿𝐿�𝑒𝑒 [0] is the DC term of theenvironment light. We will use similar techniques for 2 and 3 lights (where the equations get
much nastier if you handle a at the same time as the light intensities.)
Multiple lights can be done in a similar manner, first by removing DC from all the light vectors,
with two lights giving you this error function:
2

𝐸𝐸(𝑐𝑐0 , 𝑐𝑐1 ) = ��𝑐𝑐0 𝐿𝐿𝑑𝑑0 (𝑠𝑠) + 𝑐𝑐1 𝐿𝐿𝑑𝑑1 (𝑠𝑠) − 𝐿𝐿𝑒𝑒 (𝑠𝑠)� 𝑑𝑑𝑑𝑑
3 1
The Hessian matrix is �� which has positive eigenvalues (2 ± √2), so this is a minimum.
1 1
15
SH vector only has a non-zero in the DC term.
14

Then differentiate with respect to each of the coefficients and solve for zero. The following
expression gives you the intensities:
𝐴𝐴
𝑐𝑐0
2
2
�𝑐𝑐 � = �𝐴𝐴 − 𝐵𝐵−𝐵𝐵
1
𝐴𝐴2 − 𝐵𝐵2

Where:−𝐵𝐵

𝐴𝐴2 − 𝐵𝐵2𝐴𝐴

𝐴𝐴2 − 𝐵𝐵2

𝐿𝐿 °𝐿𝐿� � 𝑒𝑒 𝑑𝑑0 �
𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑0

𝐴𝐴 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑0 , 𝐵𝐵 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑1

It is worth pointing out that A is constant independent of direction (it depends on order and if
the vectors are convolved or not.) The ambient term is then:𝑎𝑎 = �𝐿𝐿𝑒𝑒 [0] − (𝑐𝑐0 𝐿𝐿𝑑𝑑0 [0] + 𝑐𝑐1 𝐿𝐿𝑑𝑑1 [0])�

For three lights the intensities are:

Where:

𝐴𝐴2 − 𝐷𝐷 2
⎡
𝐸𝐸
⎢
𝑐𝑐0
⎢
𝐶𝐶𝐶𝐶
−
𝐴𝐴𝐴𝐴
�𝑐𝑐1 � = ⎢
𝐸𝐸
𝑐𝑐2
⎢
𝐵𝐵𝐵𝐵
− 𝐴𝐴𝐴𝐴⎢⎣
𝐸𝐸

𝐶𝐶𝐶𝐶 − 𝐴𝐴𝐴𝐴
𝐸𝐸
𝐴𝐴2 − 𝐶𝐶 2
𝐸𝐸
𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴
𝐸𝐸

√𝜋𝜋
2

𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴⎤𝐸𝐸
⎥ 𝐿𝐿 °𝐿𝐿𝐵𝐵𝐵𝐵 − 𝐴𝐴𝐴𝐴 ⎥ 𝑒𝑒 𝑑𝑑0⎥ �𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑1 �
𝐸𝐸⎥ 𝐿𝐿𝑒𝑒 °𝐿𝐿𝑑𝑑2
𝐴𝐴2 − 𝐵𝐵2 ⎥⎦𝐸𝐸

𝐶𝐶 = 𝐿𝐿𝑑𝑑0 ∘ 𝐿𝐿𝑑𝑑2 , 𝐷𝐷 = 𝐿𝐿𝑑𝑑1 ∘ 𝐿𝐿𝑑𝑑2, 𝐸𝐸 = 2𝐵𝐵𝐵𝐵𝐵𝐵 + 𝐴𝐴(𝐴𝐴2 − 𝐵𝐵2 −𝐶𝐶 2 − 𝐷𝐷 2 )This matrix is symmetric, and the ambient coefficient is:𝑎𝑎 = �𝐿𝐿𝑒𝑒 [0] − (𝑐𝑐0 𝐿𝐿𝑑𝑑0 [0] + 𝑐𝑐1 𝐿𝐿𝑑𝑑1 [0] + 𝑐𝑐2 𝐿𝐿𝑑𝑑2 [0])�

Appendix A6 Least Squares Projection√𝜋𝜋2First we will show that for orthonormal basis functions, least squares projection actually isachieved by simply integrating against the basis functions. We want to find the coefficient
vector c that minimizes this expression:
2

𝐸𝐸(𝑐𝑐) = � �� 𝑐𝑐𝑖𝑖 𝑦𝑦𝑖𝑖 (𝑠𝑠) − 𝑓𝑓(𝑠𝑠)� 𝑑𝑑𝑑𝑑𝑖𝑖

This can be done by differentiating with respect to each coefficient, and solving the first
derivatives for zero 16
𝑑𝑑𝑑𝑑= 2 � �� 𝑐𝑐𝑖𝑖 𝑦𝑦𝑖𝑖 (𝑠𝑠) − 𝑓𝑓(𝑠𝑠)� 𝑦𝑦𝑘𝑘 (𝑠𝑠) 𝑑𝑑𝑑𝑑
𝑑𝑑𝑐𝑐𝑘𝑘𝑖𝑖

Because the basis functions are orthonormal we know that ∫ 𝑦𝑦𝑖𝑖 (𝑠𝑠)𝑦𝑦𝑗𝑗 (𝑠𝑠)𝑑𝑑𝑑𝑑 = 𝛿𝛿𝑖𝑖𝑖𝑖 exploiting thisfact and solving for zero we get𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦𝑦
𝑑𝑑𝑑𝑑= 0 �⎯⎯⎯� 𝑐𝑐𝑘𝑘 = � 𝑦𝑦𝑘𝑘 (𝑠𝑠)𝑓𝑓(𝑠𝑠) 𝑑𝑑𝑑𝑑
𝑑𝑑𝑐𝑐𝑘𝑘So that direct integration gives you the least squares result.

Now we derive the coefficient vector g that minimizes an error functional that includes apenalty based on the squared Laplacian integrated over the sphere. We know from above thecoefficient vector c that minimizes pure squared error. We introduce a function h that is used tosimplify indexing
2

𝑛𝑛

𝑛𝑛 2

𝑙𝑙

��Δ𝑓𝑓̃� 𝑑𝑑𝑑𝑑 = � � 𝑙𝑙 2 (𝑙𝑙 + 1)2 (𝑓𝑓𝑙𝑙𝑚𝑚 )2 = � ℎ𝑖𝑖 (𝑓𝑓𝑖𝑖 )2The error function then isDifferentiating you get

Solving for zero you get

𝑖𝑖=0

𝑙𝑙=1 𝑚𝑚 =−1

𝐸𝐸(𝑔𝑔) = ��(𝑔𝑔𝑖𝑖 − 𝑐𝑐𝑖𝑖 )2 � + 𝜆𝜆 � ℎ𝑖𝑖 𝑔𝑔𝑖𝑖 2
𝑖𝑖

𝑖𝑖

𝑑𝑑𝑑𝑑
= 2(𝑔𝑔𝑘𝑘 − 𝑐𝑐𝑘𝑘 ) + 2𝜆𝜆ℎ𝑘𝑘 𝑔𝑔𝑘𝑘
𝑑𝑑𝑔𝑔𝑘𝑘𝑔𝑔𝑘𝑘 =

𝑐𝑐𝑘𝑘(1 + 𝜆𝜆ℎ𝑘𝑘 )Appendix A7 Solving for Lamba to Reduce the Squared LaplacianThis can be done quite easily using Newtons method. The squared Laplacian is:

16

Second derivative mixed partials are all zero, unmixed partials are all 1, so the hessian is an identity
matrix and we have a minima

2

𝑛𝑛

𝑙𝑙

Δ = � � 𝑙𝑙
𝑙𝑙=1 𝑚𝑚 =−1

2 (𝑙𝑙

𝑛𝑛

𝑙𝑙

𝑛𝑛

+ 1)2 (𝑓𝑓𝑙𝑙𝑚𝑚 )2 = � 𝑙𝑙 2 (𝑙𝑙 + 1)2 � (𝑓𝑓𝑙𝑙𝑚𝑚 )2 = � 𝐿𝐿𝑙𝑙 𝐵𝐵𝑙𝑙
𝑙𝑙=1
𝑚𝑚 =−1𝑙𝑙=1

The array of L values (𝐿𝐿𝑙𝑙 = 𝑙𝑙 2 (𝑙𝑙 + 1)2 ) is static, and the array of B values can be computed
𝐵𝐵𝑙𝑙 = ∑𝑙𝑙𝑚𝑚 =−1(𝑓𝑓𝑙𝑙𝑚𝑚 )2 . Newtons method takes an initial guess (0 works well for this problem) andrefines it using the following recurrence:
𝜆𝜆𝑛𝑛+1 = 𝜆𝜆𝑛𝑛 −𝑓𝑓(𝜆𝜆𝑛𝑛 )𝑓𝑓′(𝜆𝜆𝑛𝑛 )Where f is the function we are searching for a root and f’ is the first derivative. Replacing 𝑓𝑓𝑙𝑙𝑚𝑚
with 𝑐𝑐𝑙𝑙𝑚𝑚 which includes 𝜆𝜆 and factoring you have:2𝑛𝑛𝑓𝑓(𝜆𝜆) = Δ − �
𝑓𝑓

′(𝜆𝜆)

𝑙𝑙=1

𝑛𝑛

= 2�
𝑙𝑙=1

𝐿𝐿𝑙𝑙 𝐵𝐵𝑙𝑙(1 + 𝜆𝜆𝐿𝐿𝑙𝑙 )2

𝐿𝐿𝑙𝑙 2 𝐵𝐵𝑙𝑙
(1 + 𝜆𝜆𝐿𝐿𝑙𝑙 )3

Where Δ2 is the original squared Laplacian. Iterate until a maximum number of iteration has
occurred, or the absolute value of successive approximations is below a threshold (1e-6 seems
to work well in practice.) This can be represented as a polynomial (multiplying through by the
product of all the denominators) of degree (𝑛𝑛 − 1)2 , where n is the order. This would only be
useful for quadratic and lower orders, where a closed form solution of the roots can becomputed and even then it is not clear if it would be much faster than an iterative solution.

Appendix A8 Code for Multiplying SH by Analytic Functions
Below are two functions, both defined over the hemisphere. The first is simply constant, which
can be useful in scenes with a dominant ground plane. The second is a clamped cosine in Z.
// generates 6th order SH coefficients from multiplying 6th order lighting by hemisphere in Z
void HemiMult(float *R, float *L){
R[0] = 0.433012702f * L[2] - 0.1653594569f * L[12] + 0.1036445247f * L[30] + 0.5f * L[0];R[1] = 0.5f * L[1] + 0.4192627457f * L[5] - 0.1711632992f * L[19];R[2] = -0.05412658775f * L[20] + 0.5f * L[2] + 0.2420614591f * L[6] + 0.433012702f * L[0];R[3] = 0.5f * L[3] - 0.1711632992f * L[21] + 0.4192627457f * L[7];R[4] = 0.5f * L[4] - 0.1713860232f * L[28] + 0.4133986423f * L[10];R[5] = -0.06477782793f * L[29] + 0.4192627457f * L[1] + 0.5f * L[5] + 0.2614562582f * L[11];R[6] = 0.5f * L[6] - 0.1448476267f * L[30] + 0.2420614591f * L[2] + 0.3697549864f * L[12];R[7] = -0.06477782793f * L[31] + 0.4192627457f * L[3] + 0.5f * L[7] + 0.2614562582f * L[13];R[8] = 0.5f * L[8] - 0.1713860232f * L[32] + 0.4133986423f * L[14];R[9] = 0.5f * L[9] + 0.41015625f * L[17];R[10] = 0.2685102947f * L[18] + 0.5f * L[10] + 0.4133986423f * L[4];R[11] = 0.360244363f * L[19] + 0.5f * L[11] + 0.2614562582f * L[5];R[12] = 0.2790440836f * L[20] + 0.5f * L[12] - 0.1653594569f * L[0] + 0.3697549864f * L[6];R[13] = 0.5f * L[13] + 0.360244363f * L[21] + 0.2614562582f * L[7];R[14] = 0.2685102947f * L[22] + 0.4133986423f * L[8] + 0.5f * L[14];R[15] = 0.41015625f * L[23] + 0.5f * L[15];R[16] = 0.408100316f * L[26] + 0.5f * L[16];R[17] = 0.2720668773f * L[27] + 0.5f * L[17] + 0.41015625f * L[9];R[18] = 0.2685102947f * L[10] + 0.35621916f * L[28] + 0.5f * L[18];R[19] = 0.2856107251f * L[29] - 0.1711632992f * L[1] + 0.5f * L[19] + 0.360244363f * L[11];R[20] = 0.2790440836f * L[12] + 0.5f * L[20] + 0.3498002708f * L[30] - 0.05412658775f * L[2];R[21] = 0.2856107251f * L[31] + 0.360244363f * L[13] + 0.5f * L[21] - 0.1711632992f * L[3];R[22] = 0.5f * L[22] + 0.35621916f * L[32] + 0.2685102947f * L[14];R[23] = 0.2720668773f * L[33] + 0.5f * L[23] + 0.41015625f * L[15];R[24] = 0.5f * L[24] + 0.408100316f * L[34];R[25] = 0.5f * L[25];R[26] = 0.408100316f * L[16] + 0.5f * L[26];R[27] = 0.5f * L[27] + 0.2720668773f * L[17];R[28] = 0.5f * L[28] - 0.1713860232f * L[4] + 0.35621916f * L[18];R[29] = -0.06477782793f * L[5] + 0.5f * L[29] + 0.2856107251f * L[19];R[30] = 0.3498002708f * L[20] + 0.1036445247f * L[0] - 0.1448476267f * L[6] + 0.5f * L[30];R[31] = -0.06477782793f * L[7] + 0.2856107251f * L[21] + 0.5f * L[31];R[32] = 0.5f * L[32] - 0.1713860232f * L[8] + 0.35621916f * L[22];R[33] = 0.5f * L[33] + 0.2720668773f * L[23];R[34] = 0.5f * L[34] + 0.408100316f * L[24];R[35] = 0.5f * L[35];}// generates 6th order SH coefficients from multiplying 6th order lighting by clamped cos in Z
void CosMult(float *R, float *L){// constants to make code fit
const float T1 = 0.09547032698f;
const float T2 = 0.1169267933f;
const float T3 = 0.2581988897f;
const float T4 = 0.2886751347f;
const float T5 = 0.2390457218f;
const float T6 = 0.2535462764f;
const float T7 = 0.2182178903f;
const float T8 = 0.1083940384f;
const float T9 = 0.2519763153f;
const float TA = 0.2439750183f;
const float TB = 0.3115234375f;
const float TC = 0.2512594538f;
const float TD = 0.31640625f;
R[0] = 0.25f * L[0] - 0.03125f * L[20] + T4 * L[2] + 0.1397542486f * L[6];R[1] = 0.2236067977f * L[5] + T2* L[11] - 0.02896952533f * L[29] + 0.1875f * L[1];R[2] = 0.375f * L[2] + T4 * L[0] + T3 * L[6] - 0.01495979856f * L[30] + T1 * L[12];R[3] = 0.2236067977f * L[7] - 0.02896952533f * L[31] + T2* L[13] + 0.1875f * L[3];R[4] = 0.101487352f * L[18] + 0.1889822365f * L[10] + 0.15625f * L[4];R[5] = 0.2236067977f * L[1] + T5 * L[11] + 0.3125f * L[5] + 0.09568319309f * L[19];R[6] = T6 * L[12] + T3 * L[2] + 0.3125f * L[6] + 0.113550327f * L[20] + 0.1397542486f * L[0];R[7] = T5 * L[13] + 0.09568319309f * L[21] + 0.2236067977f * L[3] + 0.3125f * L[7];R[8] = 0.15625f * L[8] + 0.101487352f * L[22] + 0.1889822365f * L[14];R[9] = 0.09068895910f * L[27] + 0.1666666667f * L[17] + 0.13671875f * L[9];R[10] = T7 * L[18] + 0.2734375f * L[10] + 0.09068895910f * L[28] + 0.1889822365f * L[4];R[11] = T8 * L[29] + TA * L[19] + T5 * L[5] + 0.30078125f * L[11] + T2* L[1];R[12] = T9 * L[20] + 0.328125f * L[12] + T6 * L[6] + 0.1028316139f * L[30] + T1 * L[2];R[13] = 0.30078125f * L[13] + T8 * L[31] + TA * L[21] + T5 * L[7] + T2* L[3];R[14] = 0.2734375f * L[14] + 0.09068895910f * L[32] + 0.1889822365f * L[8] + T7 * L[22];R[15] = 0.0906889591f * L[33] + 0.13671875f * L[15] + 0.1666666667f * L[23];R[16] = 0.1507556723f * L[26] + 0.123046875f * L[16];R[17] = 0.24609375f * L[17] + 0.1666666667f * L[9] + 0.201007563f * L[27];R[18] = 0.101487352f * L[4] + T7 * L[10] + 0.2302830933f * L[28] + 0.28125f * L[18];R[19] = TD * L[19] + 0.09568319309f * L[5] + TA * L[11] + 0.2461829819f * L[29];R[20] = T9 * L[12] + TC * L[30] + 0.113550327f * L[6] - 0.03125f * L[0] + TD * L[20];R[21] = 0.2461829819f * L[31] + TA * L[13] + 0.9568319309e-1f * L[7] + TD * L[21];R[22] = 0.101487352f * L[8] + 0.28125f * L[22] + T7 * L[14] + 0.2302830933f * L[32];R[23] = 0.201007563f * L[33] + 0.1666666667f * L[15] + 0.24609375f * L[23];R[24] = 0.1507556723f * L[34] + 0.123046875f * L[24];R[25] = 0.1127929688f * L[25];R[26] = 0.2255859375f * L[26] + 0.1507556723f * L[16];R[27] = 0.9068895910e-1f * L[9] + 0.201007563f * L[17] + 0.2631835938f * L[27];R[28] = 0.2302830933f * L[18] + 0.9068895910e-1f * L[10] + 0.30078125f * L[28];R[29] = 0.1083940385f * L[11] + 0.2461829819f * L[19] - 0.02896952533f * L[1] + TB * L[29];R[30] = 0.322265625f * L[30] + TC * L[20] - 0.01495979856f * L[2] + 0.1028316139f * L[12];R[31] = 0.2461829819f * L[21] - 0.02896952533f * L[3] + TB * L[31] + 0.1083940385f * L[13];R[32] = 0.09068895910f * L[14] + 0.2302830933f * L[22] + 0.30078125f * L[32];R[33] = 0.09068895910f * L[15] + 0.201007563f * L[23] + 0.2631835938f * L[33];R[34] = 0.2255859375f * L[34] + 0.1507556723f * L[24];R[35] = 0.1127929688f * L[35];}Appendix A9 Ambient Cube Basis
The Ambient Cube basis is used by Valve [26]; it consists of 6 basis functions, each defined overa hemisphere:2

𝑉𝑉0 = �𝑥𝑥0𝑥𝑥 > 0 , 𝑉𝑉1 = 0� 2𝑥𝑥
𝑥𝑥 ≤ 0𝑥𝑥 > 0𝑥𝑥 ≤ 02
𝑉𝑉3 = �𝑧𝑧0

𝑧𝑧 > 0 , 𝑉𝑉4 = 0� 2𝑧𝑧𝑧𝑧 ≤ 0𝑧𝑧 > 0𝑧𝑧 ≤ 0𝑦𝑦 2
𝑉𝑉2 = �00𝑦𝑦 > 0, 𝑉𝑉3 = � 2𝑦𝑦 ≤ 0𝑦𝑦𝑦𝑦 > 0𝑦𝑦 ≤ 0This basis is not orthogonal, so coefficients cannot be generated by simply integrating against
the basis functions. To compute optimal projection coefficients you solve a linear least squares
problem (similar to Appendix A6 Least Squares Projection, but without orthonormal basis
functions):
2

𝐸𝐸(𝑐𝑐) = � �� 𝑐𝑐𝑖𝑖 𝑉𝑉𝑖𝑖 (𝑠𝑠) − 𝑓𝑓(𝑠𝑠)� 𝑑𝑑𝑑𝑑𝑖𝑖

Where the 𝑉𝑉𝑖𝑖 (𝑠𝑠) are the Ambient Cube basis functions. Differentiating you get:𝑑𝑑𝑑𝑑
= 2 � �� 𝑐𝑐𝑖𝑖 𝑉𝑉𝑖𝑖 (𝑠𝑠) − 𝑓𝑓(𝑠𝑠)� 𝑉𝑉𝑘𝑘 (𝑠𝑠) 𝑑𝑑𝑑𝑑
𝑑𝑑𝑐𝑐𝑘𝑘𝑖𝑖

Integration is linear so you can re-arange the terms and solve for zero getting:� 𝑐𝑐𝑖𝑖 � 𝑉𝑉𝑖𝑖 (𝑠𝑠)𝑉𝑉𝑘𝑘 (𝑠𝑠)𝑑𝑑𝑑𝑑 = � 𝑓𝑓(𝑠𝑠)𝑉𝑉𝑘𝑘 (𝑠𝑠)𝑑𝑑𝑑𝑑
𝑖𝑖

The left hand side gives you a row of a matrix A, where 𝐴𝐴𝑖𝑖𝑖𝑖 = ∫ 𝑉𝑉𝑖𝑖 (𝑠𝑠)𝑉𝑉𝑗𝑗 (𝑠𝑠)𝑑𝑑𝑑𝑑 and the right handside will give you a vector (integrals of function times this basis function.) This leads to a linear
system: 𝐴𝐴𝐴𝐴 = 𝑏𝑏 the inverse of A is: 11 4 π 1
 4 π 3
− 8 π 3
− 8 π
− 3 8 π 3
− 8 π

1
4π
11
4π
3
−
8π
3
−
8π
3
−
8π
3
−
8π

3
8π
3
−
8π
11
4π
1
4π
3
−
8π
3
−
8π
−

3
8π
3
−
8π
1
4π
11
4π
3
−
8π
3
−
8π
−

3
8π
3
−
8π
3
−
8π
3
−
8π
11
4π
1
4π
−

3 

8 π 
3 
−
8 π 
3 −
8 π 
3 −
8 π 

1 

4 π 
11 
4 π 
−To project a function represented using this basis into SH you just integrate the V’s against theSH basis, for quadratic spherical harmonics that results in the following matrix:

π


3


0



0


− 3 π
4


0

0


π 5
−
15

0

 π 15

 15

π
3

π
3
−

0
0
π

3
4
0
0
−

π 5
15
0
π 15
15

−

−

π
3
π

3

π
3
π

3

π
3

0

4

4

0

0

0

0

0

0
0

0
0

0
0

π 5
15
0

2 π 5
15
0

π 15
15

0

π 5
15
0
π 15
15

−

−

π

3
4





0

3 π 

−
4 


0


0

0

2 π 5 

15 

0


0


The even degrees above quadratic are in the null-space of the Ambient Cube basis, and the odd
degrees above quadratic are in the null space of the clamped cosine function, so if irradianceenvironment maps are being used, no higher order is necessary. The Ambient Cube basis canexactly reconstruct the DC term, approximate the linear term, exactly reconstructs two of thequadratic basis functions (𝑦𝑦20 , 𝑦𝑦22 ) and has three of the quadratic basis functions in the null space(𝑦𝑦2−2 , 𝑦𝑦2−1 , 𝑦𝑦21 ).
To project from SH into the Ambient Cube basis you need to minimize this error function:

2

𝐸𝐸(𝑐𝑐) = � �� 𝑐𝑐𝑖𝑖 𝑉𝑉𝑖𝑖 (𝑠𝑠) − � 𝑙𝑙𝑗𝑗 𝑦𝑦𝑗𝑗 (𝑠𝑠)� 𝑑𝑑𝑑𝑑𝑖𝑖

𝑗𝑗

Where 𝑙𝑙𝑗𝑗 are the SH lighting coefficients (assumed to be convolved.) Once again differentiatewith respect to the unkown variables:

And then solve for zero:

𝑑𝑑𝑑𝑑= 2 � �� 𝑐𝑐𝑖𝑖 𝑉𝑉𝑖𝑖 (𝑠𝑠) − � 𝑙𝑙𝑗𝑗 𝑦𝑦𝑗𝑗 (𝑠𝑠)� 𝑉𝑉𝑘𝑘 (𝑠𝑠) 𝑑𝑑𝑑𝑑
𝑑𝑑𝑐𝑐𝑘𝑘𝑖𝑖

𝑗𝑗� 𝑐𝑐𝑖𝑖 � 𝑉𝑉𝑖𝑖 (𝑠𝑠)𝑉𝑉𝑘𝑘 (𝑠𝑠)𝑑𝑑𝑑𝑑 = � 𝑙𝑙𝑗𝑗 � 𝑦𝑦𝑗𝑗 (𝑠𝑠)𝑉𝑉𝑘𝑘 (𝑠𝑠)𝑑𝑑𝑑𝑑
𝑖𝑖

𝑗𝑗This results in a linear system as well, 𝐴𝐴𝐴𝐴 = 𝐵𝐵𝐵𝐵 , where A is the same as before B is a matrix
where 𝐵𝐵𝑖𝑖𝑖𝑖 = ∫ 𝑉𝑉𝑖𝑖 (𝑠𝑠)𝑦𝑦𝑗𝑗 (𝑠𝑠)𝑑𝑑𝑑𝑑. The matrix 𝐴𝐴−1 𝐵𝐵 can be used to move from SH to the Ambient
Cube basis:

 1
 2 π
 1
 2 π
 1
 2 π
 1
 2 π
 1
 2 π
 1

2 π

0

0

−

5 3
8 π

0 0 −

5
4 π

0

0

0

5 3
8 π

0 0 −

5
4 π

0

5 3
8 π

0

0

0 0 −

5
4 π

0

5 3
8 π

0

0

0 0 −

5
4 π

0

0

5 3
8 π

0

0 0

5
2 π

0

0

−

5 3
8 π

0

0 0

5
2 π

0

−

15 

4 π 

15 

4 π 

15 −

4 π 

15 −

4 π 

0 


0 


Appendix A10 Shader/CPU code for Irradiance Environment Maps
Given a quadratic SH representation of a lighting environment, it is fairly simple to generate
shader code. In [35] a matrix representation is used, however this turns out to require more
instructions 17 (15 vs. 11) and more constants (12 vs. 7) compared to a more direct evaluation of
the quadratic spherical harmonics. You fold leading constants of the polynomials into thelighting coefficients, group all but 𝑦𝑦2−2 by channel (using float4’s) and keep 𝑦𝑦2−2 as a color and17

On scalar GPU’s this gap is even more significant – 60 vs. 42

fold part of 𝑦𝑦20 into DC. Shader code for evaluation and CPU code to setup the constants isshown below. The normal passed in to shader code should be normalized and the 4th channel
should be 1.0. The CPU code takes an array of 3 float pointers to quadratic radiance SH
coefficients and an Effect to bind the constants to.
// constants containing irradiance environmne map
float4 cAr;
float4 cAg;
float4 cAb;
float4 cBr;
float4 cBg;
float4 cBb;
float4 cC;
float3 ShadeIrad(float4 vNormal){
float3 x1, x2, x3;// Linear + constant polynomial terms
x1.r = dot(cAr,vNormal);
x1.g = dot(cAg,vNormal);
x1.b = dot(cAb,vNormal);// 4 of the quadratic polynomials
float4 vB = vNormal.xyzz * vNormal.yzzx;
x2.r = dot(cBr,vB);
x2.g = dot(cBg,vB);
x2.b = dot(cBb,vB);// Final quadratic polynomial
float vC = vNormal.x*vNormal.x - vNormal.y*vNormal.y;x3 = cC.rgb * vC;
return x1+x2+x3;}//----------------------------------------------------------void SetSHEMapConstants( float* fLight[3], ID3DXEffect* pEffect){// Lighting environment coefficients
D3DXVECTOR4 vCoeff[3];static const float s_fSqrtPI = ((float)sqrtf(D3DX_PI));
const float fC0 = 1.0f/(2.0f*s_fSqrtPI);
const float fC1 = (float)sqrt(3.0f)/(3.0f*s_fSqrtPI);
const float fC2 = (float)sqrt(15.0f)/(8.0f*s_fSqrtPI);
const float fC3 = (float)sqrt(5.0f)/(16.0f*s_fSqrtPI);
const float fC4 = 0.5f*fC2;
int iC;
for( iC=0; iC<3; iC++ ){vCoeff[iC].x = -fC1*fLight[iC][3];vCoeff[iC].y = -fC1*fLight[iC][1];vCoeff[iC].z = fC1*fLight[iC][2];vCoeff[iC].w = fC0*fLight[iC][0] - fC3*fLight[iC][6];}pEffect->SetVector( "cAr", &vCoeff[0] );pEffect->SetVector( "cAg", &vCoeff[1] );pEffect->SetVector( "cAb", &vCoeff[2] );for( iC=0; iC<3; iC++ ){vCoeff[iC].x =fC2*fLight[iC][4];vCoeff[iC].y =-fC2*fLight[iC][5];vCoeff[iC].z = 3.0f*fC3*fLight[iC][6];vCoeff[iC].w =-fC2*fLight[iC][7];}pEffect->SetVector( "cBr", &vCoeff[0] );pEffect->SetVector( "cBg", &vCoeff[1] );pEffect->SetVector( "cBb", &vCoeff[2] );vCoeff[0].x = fC4*fLight[0][8];vCoeff[0].y = fC4*fLight[1][8];vCoeff[0].z = fC4*fLight[2][8];vCoeff[0].w = 1.0f;
pEffect->SetVector( "cC", &vCoeff[0] );}Version History
•
•
•
•
•

2/14/2008
4/1/2008
3/17/200911/30/2009
2/10/2010

first public version
Added appendix with recurrence relations when evaluating SH, renumbered appendicesFixed typos in the recurrence relations
Ambient light and optimal linear direction typos fixedTypo in defining product matrix (thanks to Don Holden for pointing it out)