# SG Series Part 2: Spherical Gaussians 101

You can find an ad-free static site version of this post here: https://therealmjp.github.io/posts/sg-series-part-2-spherical-gaussians-101/

This is part 2 of a series on Spherical Gaussians and their applications for pre-computed lighting. You can find the other articles here:

In the previous article, I gave a quick rundown of some of the available techniques for representing a pre-computed distribution of radiance or irradiance for each lightmap texel or probe location. In this article, I’m going cover the basics of Spherical Gaussians, which are a type of spherical radial basis function (SRBF for short). The concepts introduced here will serve as the core set of tools for working with Spherical Gaussians, and in later articles I’ll demonstrate how you can use those tools to form an alternative for approximating incoming radiance in pre-computed lightmaps or probes.

I should point out that this article is still going to be somewhat high-level, in that it won’t provide full derivations and background details for all formulas and operations. However it is my hope that the material here will be sufficient to gain a basic understanding of SG’s, and also use them in practical scenarios.

### What’s a Spherical Gaussian?

A Spherical Gaussian, or “SG” for short, is essentially a Gaussian function that’s defined on the surface of a sphere. If you’re reading this, then you’re probably already familar with how a Gaussian function works in 1D: you compute the distance from the center of the Gaussian, and use this distance as part of a base-e exponential. This produces the characteristic “hump” that you see when you graph it:

A Gaussian in 1D centered at x=0, with a height of 3

You’re probably also familiar with how it looks in 2D, since it’s very commonly used in image processing as a filter kernel. It ends up looking like what you would get if you took the above graph and revolved it around its axis

A Gaussian filter applied to a 2D image of a white dot, showing that the impulse response is effectively a Gaussian function in 2D

A Spherical Gaussian still works the same way, except that it now lives on the surface of a sphere instead of on a line or a flat plane. If you’re having trouble visualizing that, imagine if you took the above image and wrapped it around a sphere like wrapping paper. It ends up looking like this:

A Spherical Gaussian visualized on the surface of a sphere

Since an SG is defined on a sphere rather than a line or plane, it’s parameterized differently than a normal Gaussian. A 1D Gaussian function always has the following form: $ae^{\frac{-(x - b)^{2}}{2c^{2}}}$

The part that we need to change in order to define the function on a sphere is the “(x – b)” term. This part of the function essentially makes the Gaussian a function of the cartesian distance between a given point and the center of the Gaussian, which can be trivially extended into 2D using the standard distance formula. To make this work on a sphere, we must instead make our Gaussian a function of the angle between two unit direction vectors. In practice we do this by making an SG a function of the cosine of the angle between two vectors, which can be efficiently computed using a dot product like so: $G(\mathbf{v};\mathbf{\mu},\lambda,a) = ae^{\lambda(\mathbf{\mu} \cdot \mathbf{v} - 1)}$

Just like a normal Gaussian, we have a few parameters that control the shape and location of the resulting lobe. First we have μ, which is the axis, or direction of the lobe. It effectively controls where the lobe is located on the sphere, and always points towards the exact center of the lobe. Next we have λ, which is the sharpness of the lobe. As this value increases, the lobe will get “skinnier”, meaning that the result will fall off more quickly as you get further from the lobe axis. Finally we have a, which is the amplitude or intensity of the lobe. If you were to look at a polar graph of an SG, it would correspond to the height of the lobe at its peak. The amplitude can be a scalar value, or for graphics applications we may choose to make it an RGB triplet in order to support varying intensities for different color channels. This all lends itself to a simple HLSL code definition:

struct SG
{
float3 Amplitude;
float3 Axis;
float Sharpness;
};


Evaluating an SG is also easily expressible in HLSL. All we need is a normalized direction vector representing the point on the sphere where we’d like to compute the value of the SG:

float3 EvaluateSG(in SG sg, in float3 dir)
{
float cosAngle = dot(dir, sg.Axis);
return sg.Amplitude * exp(sg.Sharpness * (cosAngle - 1.0f));
}


### Why Spherical Gaussians?

Now that we know what a Spherical Gaussian is, what’s so useful about them anyway? One pontential benefit is that they’re fairly intuitive: it’s not terribly hard to understand how the 3 parameters work, and how each parameter affects the resulting lobe. The other main draw is that they inherit a lot of useful properties of “regular” Gaussians, which makes them useful for graphics and other related applications. These properties have been explored and utilized in several research papers that were primarily aimed at achieving pre-computed radiance transfer (PRT) with both diffuse and specular material response. In particular,  the paper entitled “All-Frequency Rendering of Dynamic, Spatially-Varying Reflectance” by Wang et al. was our main inspiration for pursuing SG’s at RAD.

### Products

So what are these useful Gaussian properties that we can exploit? For starters, taking the product of 2 Gaussians functions produces another Gaussian. For an SG, this is equivalent to visiting every point on the sphere, evaluating 2 different SG’s, and multiplying the two results. Since it’s an operation that takes 2 SG’s and produces another SG, it is sometimes referred to as a “vector” product. It’s defined as the following: $G_{1}(\mathbf{v})G_{2}(\mathbf{v}) = G(\mathbf{v};\frac{\mu_{m}}{||\mu_{m}||},a_{1}a_{2}e^{\lambda_{m}(||\mu_{m}|| - 1)})$ $\lambda_{m} = \lambda_{1} + \lambda_{2}$ $\mu_{m} = \frac{\lambda_{1}\mu_{1} + \lambda_{2}\mu_{2}}{\lambda_{1} + \lambda_{2}}$

In HLSL code, it looks like this:

SG SGProduct(in SG x, in SG y)
{
float3 um = (x.Sharpness * x.Axis + y.Sharpness * y.Axis) /
(x.Sharpness + y.Sharpness);
float umLength = length(um);
float lm = x.Sharpness + y.Sharpness;

SG res;
res.Axis = um * (1.0f / umLength);
res.Sharpness = lm * umLength;
res.Amplitude = x.Amplitude * y.Amplitude *
exp(lm * (umLength - 1.0f));

return res;
}


### Integrals

Gaussians have another really nice property in that their integrals have a closed-form solution, which is known as the error function. The property also extends to SG’s, where we can compute the integral of an SG over the entire sphere: $\int_{\Omega} G(\mathbf{v})d\mathbf{v} = 2\pi\frac{a}{\lambda}(1 - e^{-2\lambda})$

Computing an integral will essentially tell us the total “energy” of an SG, which can be useful for lighting calculations. It can also be useful for normalizing an SG, which produces an SG that integrates to 1. Such a normalized SG is suitable for representing a probability distribution, such as an NDF. In fact, a normalized SG is actually equivalent to a von Mises-Fisher distribution in 3D!

An SG integral is actually very cheap to compute…or at least it would be if we removed the exponential term. It turns out that the $(1 - e^{-2\lambda})$ term actually approaches 1 very quickly as the SG’s sharpness increases, which means we can potentially drop it with little error as long as we know that the sharpness is high enough. Here’s what a graph of $(1 - e^{-2\lambda})$ looks like for increasing sharpness:

A graph of the exponential term in computing the integral of an SG over a sphere, which approaches 1 as the sharpness increases. The X-axis is sharpness, and the Y-axis is the value of $(1 - e^{-2\lambda})$.

This all lends itself naturally to HLSL implementations for accurate and approximate versions of an SG integral:

float3 SGIntegral(in SG sg)
{
float expTerm = 1.0f - exp(-2.0f * sg.Sharpness);
return 2 * Pi * (sg.Amplitude / sg.Sharpness) * expTerm;
}

float3 ApproximateSGIntegral(in SG sg)
{
return 2 * Pi * (sg.Amplitude / sg.Sharpness);
}


### Inner Product

If we were to use our SG integral formula to compute the integral of the product of two SG’s, we can compute what’s known as the inner product, or dot product of those SG’s. The operation is usually defined like this: $\int_{\Omega} G_{1}(\mathbf{v}) G_{2}(\mathbf{v}) d\mathbf{v} = \frac{4 \pi a_{0} a_{1}}{e^{\lambda_{m}}} \frac{sinh(d_{m})}{d_{m}}$ $d_{m} = || \lambda_{1}\mu_{1} + \lambda_{2}\mu_{2} ||$

However we can avoid numerical precision issues by using an alternate arrangement: $\int_{\Omega} G_{1}(\mathbf{v}) G_{2}(\mathbf{v}) d\mathbf{v} = 2 \pi a_{0} a_{1}\frac{e^{d_{m} - \lambda_{m}} - e^{-d_{m} - \lambda_{m}}}{d_{m}}$

…which looks like this in HLSL:

float3 SGInnerProduct(in SG x, in SG y)
{
float dm = length(x.Sharpness * x.Axis + y.Sharpness * y.Axis);
float3 expo = exp(dm - x.Sharpness - y.Sharpness) * x.Amplitude * y.Amplitude;
float other = 1.0f - exp(-2.0f * dm);
return (2.0f * Pi * expo * other) / dm;
}


### Treshold

SG’s have what’s known as “compact-ε” support, which means that it’s possible to determine an angle θ such that all points within θ radians of the SG’s axis will have a value greater than ε. This property is potentially more useful if we flip it around so that we calculate a sharpness λ that results in a given θ for a particular value of ε: $a e^{\lambda(cos\theta - 1)} = \epsilon$ $\lambda = \frac{ln(\epsilon) - ln(a)}{cos\theta - 1}$

float SGSharpnessFromThreshold(in float amplitude,
in float epsilon,
in float cosTheta)
{
return (log(epsilon) - log(amplitude)) / (cosTheta - 1.0f);
}


### Rotation

One last operation I’ll discuss is rotation. Rotating an SG is trivial: all you need to do is apply your rotation transform to the SG’s axis vector and you have a rotated SG! You can apply the transform using a matrix, a quaternion, or any other means you might have for rotating a vector. This is a welcome change from SH, which requires a very complex transform once you go above L1.

 Gaussian Function – https://en.wikipedia.org/wiki/Gaussian_function
 All-Frequency Rendering of Dynamic, Spatially-Varying Reflectance – http://research.microsoft.com/en-us/um/people/johnsny/papers/sg.pdf
 Error Function – https://en.wikipedia.org/wiki/Error_function
 von-Mises Fisher Distribtion – https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution

## 5 thoughts on “SG Series Part 2: Spherical Gaussians 101”

1. ChenA says:

G_{1}(\mathbf{v})G_{2}(\mathbf{v}) = G(\mathbf{v};\frac{\mu_{m}}{||\mu_{m}||},a_{1}a_{2}e^{\lambda_{m}(||\mu_{m}|| – 1)})
is it lost the sharpness λm∥pm∥?

and
SG SGProduct(in SG x, in SG y)
{
float3 um = (x.Sharpness * x.Axis + y.Sharpness + y.Axis) /
(x.Sharpness + y.Sharpness);

is this should be:
float3 um = (x.Sharpness * x.Axis + y.Sharpness * y.Axis) /
(x.Sharpness + y.Sharpness);

2. MJP says:

Yes indeed, that was an error on my part. Thank you for pointing that out!

3. Mike R. says:

Just wanted say thanks for the awesome explanations, you make everything so clear and easy to understand.

4. Nikita Lisitsa (@lisyarus) says:

Saying that the integral of a Gaussian has a closed-form solution in the form of an error function is unfair – the error function is pretty much defined to be that integral anyway. You also probably mean the definite integral over the whole range (the real axis in the case of a 1D Gaussian), which for Gaussians doesn’t involve the error function: https://en.wikipedia.org/wiki/Gaussian_integral

5. Oh okay thanks. Loved the vibrant colour