# SG Series Part 4: Specular Lighting From an SG Light Source

You can find an ad-free static site version of this post here: https://therealmjp.github.io/posts/sg-series-part-4-specular-lighting-from-an-sg-light-source/

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

In the previous article, we explored a few ways to compute the contribution of an SG light source when using a diffuse BRDF. While this is already useful, it would be nice to be able work with more complex view-dependent BRDF’s so that we can also compute a specular contribution. In this article I’ll explain a possible approach for approximating the response of a microfacet specular BRDF when applied to an SG light, and also introduce the concept of Anisotropic Spherical Gaussians.

### Microfacet Specular Recap

You probably recall from the past 5 years of physically based rendering presentations that a standard microfacet BRDF takes the following structure: $f(\mathbf{i}, \mathbf{o}) = \frac{F(\mathbf{o}, \mathbf{h})\,G(\mathbf{i}, \mathbf{o}, \mathbf{h})\,D(\mathbf{h})} {4\,(\mathbf{n} \cdot \mathbf{i})\,(\mathbf{n} \cdot \mathbf{o})}$

Recall that D(h) is the distribution term, which tells us the percentage of active microfacets for a particular combination of view and light vectors. It’s also commonly known as the normal distribution function, or “NDF” for short. It’s generally parameterized on a roughness parameter, which essentially describes the “bumpiness” of the underlying microgeometry. Lower roughness values lead to sharp mirror-like reflections with a very narrow and intense specular lobe, while higher roughness values lead to more broad reflections with a wider specular lobe. Most modern games (including The Order: 1886) use the GGX (AKA Trowbridge-Reitz) distribution for this purpose.  The top image shows a GGX distribution term with a roughness parameter of 0.5. The bottom image shows the same GGX distribution term with a roughness parameter of 0.1. For both graphs, the X axis represents the angle between the surface normal and the half vector. Click on either image for an interactive graph.

Next we have G(i, o, h) which is referred to as the geometry term or the masking-shadow function. Eric Heitz’s paper on the masking-shadow function has a fantastic explanation of how these terms work and why they’re important, so I would strongly recommend reading through it if you haven’t already. As Eric explains, the geometry term actually accounts for two different phenomena. The first is the local occlusion of reflections by other neighboring microfacets. Depending on the angle of the incoming lighting and the bumpiness of the microsurface, a certain amount of the lighting will be occluded by the surface itself, and this term attempts to model that. The other phenomenon handled by this term is visibility of the surface from the viewer. A surface that isn’t visible naturally can’t reflect light towards the viewer, and so the geometry term models this masking effect using the view direction and the roughness of the microsurface. The Smith visibility term for GGX as a function of the angle between the surface normal and the light direction. The roughness used is 0.25, and the angle between the normal and the view direction is 0. Click on the image for an interactive graph.

Finally we have F(o, h), which is the Fresnel term. This familiar term determines the amount of light that is reflected vs. the amount that is refracted or absorbed, which varies depending on the index of refraction for a particular interface as well as the angle of incidence. For a microfacet BRDF we compute the Fresnel term using the angle between the active microfacet direction (the half vector) and the light or view direction (it is equivalent to use either). This generally causes the reflected intensity  to increase when both the viewing direction and the light direction are grazing with respect to the surface normal. In real-time graphics Schlick’s approximation is typically used to represent the Fresnel term, since it’s a bit cheaper than evaluating the actual Fresnel equations. It is also common to parameterize the Fresnel term on the reflectance at zero incidence (referred to as “F0”) instead of directly working with an index of refraction. Schlick’s approximation of the Fresnel term as a function of the angle between the half vector and the light direction. The graph uses a value of 0.04 for F0. Click on the image for an interactive graph.

### Specular From an SG Light

Let’s now return to our example of a Spherical Gaussian light source. In the previous article we explored how to approximate the purely diffuse contribution from such a light source, but how do we handle the specular response? If we go down this road, we ideally want to do it in a way that uses a microfacet specular model so that we can remain consistent with our lighting from punctual light sources and IBL probes. If we start out by just plugging our BRDF and SG light into the rendering equation we get this monstrosity: $L_{o}(\mathbf{o}, \mathbf{x}) = \int_{\Omega} \frac{F(\mathbf{o}, \mathbf{h})\,G(\mathbf{i}, \mathbf{o}, \mathbf{h})\,D(\mathbf{h})G_{l}(\mathbf{i};\mathbf{\mu},\lambda,a)cos(\theta_{i})d\Omega} {4cos(\theta_{i})cos(\theta_{o})}$ $L_{o}(\mathbf{o}, \mathbf{x}) = \int_{\Omega} \frac{F(\mathbf{o}, \mathbf{h})\,G(\mathbf{i}, \mathbf{o}, \mathbf{h})\,D(\mathbf{h})ae^{\lambda(\mathbf{\mu} \cdot \mathbf{i} - 1)}d\Omega} {4cos(\theta_{o})}$

Unlike diffuse we now have multiple terms inside the integral, many of which are view-dependent. This suggests that we will need to combine several aggressive optimizations in order to get anywhere close to our desired result.

#### The Distribution Term

Let’s start out by taking the distribution term on its own, since it’s arguably the most important part of the specular BRDF. The distribution determines the overall shape and intensity of the resulting specular highlight, and has a widely-varying frequency response over the domain of possible roughness values. If we look at the graph of the GGX distribution from the earlier section, the shape does resemble a Gaussian. Unfortunately it’s not an exact match: the GGX distribution has the characteristic narrow highlight and long  tails which we can’t match using a single SG. We could get a closer fit by summing multiple Gaussians, but for this article we’ll keep things simple by sticking with a single lobe. If we go back to Wang et al.‘s paper that we referenced earlier, we can see that they suggest a very simple fit of an SG to a Cook-Torrance distribution term: $D(\mathbf{h})= e^{-(arccos(\mathbf{h} \cdot \mathbf{n})/m)^{2}} \approx G(\mathbf{h};\mathbf{n},\frac{2}{m^2},1)$

If we look closely, we can see that they’re actually fitting for the Gaussian model mentioned in the original Cook-Torrance paper, which is similar to the one used in the Torrance-Sparrow model. Note that this should not be confused the with Beckmann distribution that’s also mentioned in the Cook-Torrance paper, which is actually a 2D Gaussian in the slope domain (AKA parallel plane domain). However the shape isn’t even the biggest problem here, as the variant they’re using isn’t normalized. This means the peak of the Gaussian will always have a height of 1.0, rather than shrinking when the roughness increases and growing when the roughness decreases. Fortunately this is really easy to fix, since we have a simple analytical formula for computing the integral of an SG. Therefore if we set the amplitude to 1 over the integral, we end up with a normalized distribution: $D_{SG}(\mathbf{h})= G(\mathbf{h};\mathbf{n},\frac{2}{m^2}, \frac{1}{\pi m^2})$

SG DistributionTermSG(in float3 direction, in float roughness)
{
SG distribution;
distribution.Axis = direction;
float m2 = roughness * roughness;
distribution.Sharpness = 2 / m2;
distribution.Amplitude = 1.0f / (Pi * m2);

return distribution;
}

Let’s take a look at a graph of our distribution term, and see how far off it is from our target:  Top graph shows a comparison between GGX, Beckmann, normalized Gaussian, and SG distribution terms with a roughness of 0.25. The bottom shows the same comparison with  a roughness of 0.5. Click on either image for an interactive graph.
It should come as no surprise that our SG distribution is almost an exact match for a normalized version of a Gaussian distribution. When the roughness is lower it’s also a fairly close match for Beckmann, but for higher roughness the difference gets to be quite large. In both cases our distribution isn’t a perfect fit for GGX, but it’s a workable approximation.
Warping Between Domains
So we now have a normal distribution function in SG format, but unfortunately we’re not quite ready to use it as-is. The problem is that we’ve defined our distribution in the half-vector domain: the axis of the SG points in the direction of the surface normal, and we use the half-vector as our sampling direction. If we want to use an SG product to compute the result of multiplying our distribution with an SG light source, then we need to ensure that the distribution lobe is in the same domain as our light source. Another way to think about this is that center of our distribution shifts depending on viewing angle, since the half-vector also shifts as the camera moves.
In order to make sure that our distribution lobe is in the correct domain, we need “warp” our distribution so that it lines up with the current BRDF slice for our viewing direction. If you’re wondering what a BRDF slice is, it basically tells you “if I pick a particular view direction. what is the value of my  BRDF for a given light direction?”. So if you had a mirror BRDF, the slice would just be a single ray pointing in the direction of the view ray reflected off the surface normal. For microfacet specular BRDF’s, you get a lobe that’s roughly centered around the reflected view direction. Here’s what a polar graph of a GGX BRDF slice looks like if we only consider the distribution term: Polar graph of the GGX distribution term from two different viewing angles. The blue line is the view direction, the green line is the surface normal, and the pink line is the reflected view direction. The top image shows the resulting slice when the viewing angle is 0 degrees, and the bottom images shows the resulting slice when the viewing angle is 45 degrees.
Wang et al. proposed a simple spherical warp operator that would orient the distribution lobe about the reflected view direction, while also modifying the sharpness to take into account the differential area at the original center of the lobe: $\mu_{w}=2(\mathbf{o} \cdot \mu_{d})\mu_{d} - \mathbf{o}$ $\lambda_{w}=\frac{\lambda_{d}}{4|\mu_{d} \cdot \mathbf{o}|}$ $a_{w} = a_{d}$
SG WarpDistributionSG(in SG ndf, in float3 view)
{
SG warp;

warp.Axis = reflect(-view, ndf.Axis);
warp.Amplitude = ndf.Amplitude;
warp.Sharpness = ndf.Sharpness;
warp.Sharpness /= (4.0f * max(dot(ndf.Axis, view), 0.0001f));

return warp;
}

Let’s now look at the result of that warp, and compare it to what the actual GGX distribution looks like: Result of applying a spherical warp to the SG distribution term (green) compared with the actual GGX distribution (red). The top graph shows a viewing angle of 0 degrees, and the bottom graph shows a viewing angle of 45 degrees.
A quick look at the graph shows us that the shape is a bit off, but our warped lobe is ending up in approximately in the right spot. We’ll revisit the lobe shape later, but for now let’s try combining our distribution with the rest of our BRDF.
Approximating The Remaining Terms
In the previous section we figured out how to obtain an SG approximation of our distribution term, and also warp it so that it’s in the correct domain for multiplication with our SG light source. Using our SG product operator would allow to us to represent the result of multiplying the distribution with the light source as another SG, which we could then multiply with other terms using SG operators…or at least we could if we were to represent the remaining terms as SG’s. Unfortunately this turns out to be a problem: the geometry and Fresnel terms are nothing at all like a Gaussian, which rules out approximating them as an SG. Wang et al. sidestep this issue by making the somewhat-weak assumption that the values of these terms will be constant across the entire BRDF lobe, which allows them to pull the terms out of the integral and evaluate them only for the axis direction of the BRDF lobe. This allows the resulting BRDF to still capture some of the glancing angle effects, with similar performance cost to evaluating those terms for a punctual light source. The downside is that the error of these terms will increase as the BRDF lobe becomes wider (increasing roughness), since the value of the geometry and Fresnel terms will vary more the further they are from the lobe center. Putting it all together gives the following specular BRDF: $f_{sg}(\mathbf{i}, \mathbf{o}) = \frac{F(\mathbf{o}, \mathbf{h_{w}})\,G(\mathbf{\mu_{w}}, \mathbf{o}, \mathbf{h_{w}})\, \frac{1}{\pi m^2}e^{\lambda_{w}(\mathbf{\mu_{w}} \cdot \mathbf{i} - 1)}} {4\,(\mathbf{n} \cdot \mathbf{\mu_{w}})\,(\mathbf{n} \cdot \mathbf{o})}$ $\mathbf{h_{w}} = \frac{\mathbf{o} + \mathbf{\mu_{w}}}{||\mathbf{o} + \mathbf{\mu_{w}}||}$
The last thing we need to account for is the cosine term that needs to be multiplied with the BRDF inside of the hemispherical integral. Wang et al. suggest using an SG product to compute an SG representing the result of multiplying the distribution term SG and their SG approximation of a clamped cosine lobe, which can then be multiplied with the lighting lobe using an SG inner product. In order to avoid another expensive SG operation, we will instead use the same approach that we used for geometry and Fresnel terms and evaluate the cosine term using the BRDF lobe axis direction. Implementing it in shader code gives us the following:
float GGX_V1(in float m2, in float nDotX)
{
return 1.0f / (nDotX + sqrt(m2 + (1 - m2) * nDotX * nDotX));
}

float3 SpecularTermSGWarp(in SG light, in float3 normal,
in float roughness, in float3 view,
in float3 specAlbedo)
{
// Create an SG that approximates the NDF.
SG ndf = DistributionTermSG(normal, roughness);

// Warp the distribution so that our lobe is in the same
// domain as the lighting lobe
SG warpedNDF = WarpDistributionSG(ndf, view);

// Convolve the NDF with the SG light
float3 output = SGInnerProduct(warpedNDF, light);

// Parameters needed for the visibility
float3 warpDir = warpedNDF.Axis;
float m2 = roughness * roughness;
float nDotL = saturate(dot(normal, warpDir));
float nDotV = saturate(dot(normal, view));
float3 h = normalize(warpedNDF.Axis + view);

// Visibility term evaluated at the center of
// our warped BRDF lobe
output *= GGX_V1(m2, nDotL) * GGX_V1(m2, nDotV);

// Fresnel evaluated at the center of our warped BRDF lobe
float powTerm = pow((1.0f - saturate(dot(warpDir, h))), 5);
output *= specAlbedo + (1.0f - specAlbedo) * powTerm;

// Cosine term evaluated at the center of the BRDF lobe
output *= nDotL;

return max(output, 0.0f);
}

Let’s now (finally) take a look at what our specular approximation looks like for a scene being lit by an SG light source: A scene being lit by an SG light source using our diffuse and specular approximations. The scene is using a uniform roughness of 0.128.
Going Anisotropic
So if we look at the specular highlights in the above image, you may notice that while the highlights on the red and green walls look pretty good, the highlight on floor seems a bit off. The highlight is rather wide and rounded, and our intuition tells that a highlight viewed at a grazing angle should appear vertically stretched across the floor. To determine why the look is so off, we need to revisit our warp of the distribution term. Previously when we looked at the polar graph of the distribution I noted that the shape of the resulting lobe was a bit off, even though it was oriented in approximately the right direction to line up with the BRDF slice. To get a better idea of what’s going on, let’s now take a look at a 3D graph of the GGX distribution term: 3D graph of the GGX distribution term. The left image shows the distribution when the viewing angle is very low, while the right image shows the distribution when the viewing angle is very steep.
Looking at the left image where the angle between the view direction and  the surface normal are very low, the distribution is radially symmetrical just like an SG lobe. However as the viewing angle increases the lobe begins to stretch, looking more and more like the non-symmetrical lobe that we see in the right image. The stretching of the lobe is what causes the stretched highlights that occur when applying the BRDF, and our warped SG is unable to properly represent it since it must remain radially symmetric about its axis.
Luckily for us, there is a better way. In 2013 Xu et al.released a paper titled Anisotropic Spherical Gaussians, where they explain how they extended SG’s to support anisotropic lobe width/sharpness. They’re defined like this: $G(\mathbf{v};[\mu_x, \mu_y, \mu_z],[\lambda_x, \lambda_y],a) = a \cdot S(\mathbf{v},\mu_z) e^{-\lambda_x(\mathbf{v} \cdot \mu_x) - \lambda_y(\mathbf{v} \cdot \mu_y)}$
Instead of having a single axis direction, an ASG now has $\mu_x$, $\mu_y$, and $\mu_z$, which are three orthogonal vectors forming a complete basis. It’s very similar to a tangent frame, where the normal, tangent, and bitangent together make up an orthonormal basis. With an ASG you also now have two separate sharpness parameters, $\lambda_x$ and $\lambda_y$, which control the sharpness with respect to $\mu_x$ and $\mu_y$. So for example setting $\lambda_x$ to 16 and $\lambda_y$ to 64 will result in stretched Gaussian lobe that’s skinnier along the $\lambda_y$ direction, and with its center located at $\mu_z$. Visualizing such an ASG on the surface of a sphere gives you this: An Anisotropic Spherical Gaussian visualized on the surface of a sphere. $\lambda_x$ has a value of 16, and $\lambda_y$ has a value of 64.
Like SG’s, the equations lend themselves to simple HLSL implementations:
struct ASG
{
float3 Amplitude;
float3 BasisZ;
float3 BasisX;
float3 BasisY;
float SharpnessX;
float SharpnessY;
};

float3 EvaluateASG(in ASG asg, in float3 dir)
{
float sTerm = saturate(dot(asg.BasisZ, dir));
float lambdaTerm = asg.SharpnessX * dot(dir, asg.BasisX)
* dot(dir, asg.BasisX);
float muTerm = asg.SharpnessY * dot(dir, asg.BasisY)
* dot(dir, asg.BasisY);
return asg.Amplitude * sTerm* exp(-lambdaTerm - muTerm);
}

The ASG paper provides us with formulas for two operations that we can use to improve the quality of our specular approximation for SG light sources. The first is a new warping operator that can take an NDF represented as an isotropic SG, and stretch it along the view direction to produce an ASG that better represents the actual BRDF. The other helpful forumula it provides is for convolving an ASG with an SG, which we can use to convolve a anisotropically warped NDF lobe with an SG lighting lobe. Let’s take a look at how their improved warp looks when graphing the NDF for a large viewing angle: The left image is a 3D graph of the distribution term when using a spherical warp. The middle image is the resulting distribution term when using an anisotropic warp. The right image is the actual GGX distribution term.
The anisotropic distribution looks much closer to the actual GGX NDF, since it now has the vertical stretching that we were missing. Let’s now implement their formulas in HLSL so we can try the new warp in our test scene:
float3 ConvolveASG_SG(in ASG asg, in SG sg) {
// The ASG paper specifes an isotropic SG as
// exp(2 * nu * (dot(v, axis) - 1)),
// so we must divide our SG sharpness by 2 in order
// to get the nup parameter expected by the ASG formula
float nu = sg.Sharpness * 0.5f;

ASG convolveASG;
convolveASG.BasisX = asg.BasisX;
convolveASG.BasisY = asg.BasisY;
convolveASG.BasisZ = asg.BasisZ;

convolveASG.SharpnessX = (nu * asg.SharpnessX) /
(nu + asg.SharpnessX);
convolveASG.SharpnessY = (nu * asg.SharpnessY) /
(nu + asg.SharpnessY);

convolveASG.Amplitude = Pi / sqrt((nu + asg.SharpnessX) *
(nu + asg.SharpnessY));

float3 asgResult = EvaluateASG(convolveASG, sg.Axis);
return asgResult * sg.Amplitude * asg.Amplitude;
}

ASG WarpDistributionASG(in SG ndf, in float3 view)
{
ASG warp;

// Generate any orthonormal basis with Z pointing in the
// direction of the reflected view vector
warp.BasisZ = reflect(-view, ndf.Axis);
warp.BasisX = normalize(cross(ndf.Axis, warp.BasisZ));
warp.BasisY = normalize(cross(warp.BasisZ, warp.BasisX));

float dotDirO = max(dot(view, ndf.Axis), 0.0001f);

// Second derivative of the sharpness with respect to how
// far we are from basis Axis direction
warp.SharpnessX = ndf.Sharpness / (8.0f * dotDirO * dotDirO);
warp.SharpnessY = ndf.Sharpness / 8.0f;

warp.Amplitude = ndf.Amplitude;

return warp;
}

float3 SpecularTermASGWarp(in SG light, in float3 normal,
in float roughness, in float3 view,
in float3 specAlbedo)
{
// Create an SG that approximates the NDF
SG ndf = DistributionTermSG(normal, roughness);

// Apply a warpring operation that will bring the SG from
// the half-angle domain the the the lighting domain.
ASG warpedNDF = WarpDistributionASG(ndf, view);

// Convolve the NDF with the light
float3 output = ConvolveASG_SG(warpedNDF, light);

// Parameters needed for evaluating the visibility term
float3 warpDir = warpedNDF.BasisZ;
float m2 = roughness * roughness;
float nDotL = saturate(dot(normal, warpDir));
float nDotV = saturate(dot(normal, view));
float3 h = normalize(warpDir + view);

// Visibility term
output *= GGX_V1(m2, nDotL) * GGX_V1(m2, nDotV);

// Fresnel
float powTerm = pow((1.0f - saturate(dot(warpDir, h))), 5);
output *= specAlbedo + (1.0f - specAlbedo) * powTerm;

// Cosine term
output *= nDotL;

return max(output, 0.0f);
}

If we swap out our old specular approximation for one that uses an anisotropic warp, our test scene now looks much better! Diffuse and specular approximations applied to an SG light source, using an anisotropic warp to approximate the NDF as an ASG.
Visually Comparing the BRDF
Many of the images that I used for visually comparing the SG BRDF approximations were generated using Disney’s BRDF Explorer. When we were doing our initial research into SG’s and figuring out how to implement them, BRDF Explorer was extremely valuable both for understanding the concepts and for experimenting with different variations. If you’d like to this yourself, there’s a very easy way to do that courtesy of Nick Brancaccio. Nick was kind of enough to create his own awesome WebGL version of BRDF Explorer, and it comes pre-loaded with options for comparing an approximate SG specular BRDF with the GGX BRDF. I would recommend checking it out if you’d like to play around with the BRDF’s and make some pretty 3D graphs!
References
 All-Frequency Rendering of Dynamic, Spatially-Varying Reflectance – http://research.microsoft.com/en-us/um/people/johnsny/papers/sg.pdf
 A Reflectance Model for Computer Graphics – http://inst.cs.berkeley.edu/~cs294-13/fa09/lectures/cookpaper.pdf
 Theory for Off-Specular Reflection From Roughened Surfaces – http://www.graphics.cornell.edu/~westin/pubs/TorranceSparrowJOSA1967.pdf
 Anisotropic Spherical Gaussians – http://cg.cs.tsinghua.edu.cn/people/~kun/asg/
 BRDF Explorer – https://github.com/wdas/brdf
 WebGL BRDF Explorer – https://depot.floored.com/brdf_explorer

Related


## 4 thoughts on “SG Series Part 4: Specular Lighting From an SG Light Source”

1. Evgenii Golubev says:

The first paragraph is missing a link. 🙂

Great series! I especially like the included graphs.

2. MJP says:

Thank you Evgenii! I added the missing link. I’m also glad that you like the graphs, it took a while to generate all of them. 🙂

3. Devsh says:

Hi MJP,

I’ve noticed that we could fit more gaussians to an arbitrary function.

For example I’ve managed to match the Beckmann distribution with 3 gaussians for the roughness value of 1
https://www.desmos.com/calculator/o1lces1o1l

My aim is to fit between 2-5 gaussians to common BRDFs and sample from environment maps pre-convolved with SG, then that way I can get a much better approximation to Importance Sampled references than the current IBL methods.

4. mmerchante says:

Hey MJP,

Just wanted to mention that from the ASG paper (and the code you presented), it seems the LaTeX definition you wrote down has the wrong exponential term, as it’s lacking the pow(dot(v, u_x), 2).

Thanks for such a good blog series btw!