Computational resources are valuable, even with that sweet, sweet GPU. Sometimes, GI at 60fps per second is simply not attainable, especially when its 2000 or something, when computational resources are way more scarce. Now, what if we time travel back to the past, and we have to render a complex scene, at 60 frames per second, with one catch that the whole scene is diffuse? This is where spherical harmonics comes in handy.
Proposed by Peter-Pike Sloan et al., Precomputed Radiance Transfer (PRT) is a real time, global illumination rendering technique using spherical harmonics. Now, I know what you are asking. Why learn some age-old technique when the GPUs right now have real time raytracing built right into their core? If not for the sake of learning, maybe the fact that it is still very much used in modern day game engines under the name Light Probes can motivate you. So, let’s dive right in!
The first is easily solvable in raytracing or raymarching, since we are sending rays everywhere and it’s bound to sample the environment map in some point. The same cannot be said for rasterization, which real time rendering uses. To understand what an environment map is, we can look at this video illustration made by Keenan Crane (or his students) (the project website is sadly offline now so I can just link to the course website):
Click to play the video. In short, an environment map is a 2D texture
In rasterization, shading a pixel using an environment map as lighting means calculating how light interacts with every pixel on the environment map, which could take up a lot of time. To understand how this works, or why do we need spherical harmonics to begin with, we need to take a look at the Rendering Equation.
Since this is a recursive algorithm, evaluating the
Putting those thoughts aside, let’s focus on the Rendering Equation now. In this rendering equation,
To start it all off, we are going to get rid of the recursive term. If the light never bounces off other objects, and comes straight from the environment map, then there is no recursive term to speak of:
Now, if we combine the scattering function, visibility test, geometry albedo & the cosine term into one function, which we call transport function
If we assume that the only light emitter in the scene is the environment map, then the
In conclusion, if we want to calculate the scene radiance of diffuse surfaces, ignoring interreflections & light emitting objects in scene, this is the final end simplification result. However, the hemispherical integration is still there, and it slows down performance. We need to find a way to get rid of it as well.
Take a look at basis vectors and its definition: in an N-dimensional space, N non-colinear vectors can represent all coordinates in said space with a weighted sum:
What we are most interested in however, are orthonormal basis vectors, which means those vectors have a length of 1, and they are all perpendicular to each other, i.e.
With a set of orthonormal basis vectors, we can propose a projection
procedure, which projects any point
reconstruct
ed back to the canonical space:
Now, to get an intuition about what we are going to do next, imagine the vector dimension goes up. As it goes up and up, there are more and more components in a vector, and eventually, as the dimensionality approaches infinite, dot operation becomes the sum of the products of infinite elements. In math terms, it becomes this:
Since infinite dimension vector is basically a function, we denote it as such. Infinite vector length (function length) is thus defined as
As a result, vector operations on functions make sense - it also has a name, Vector Calculus, and you can check it out on Wikipedia. In fact, function operators on vectors also make sense, such as the Laplace operator. But we are not covering it today.
Take a look back at the simplified rendering equation. If we can somehow perform a projection
on both reconstruction
. Although it feels like some needlessly complex extra steps, this is when magic happens. Because when we plug them back into our simplified Rendering Equation, this is what we get:
We have successfully gotten rid of the integration symbol, and the Rendering Equation devolves into a simple dot operation between two different vectors. Dot operation is, of course, very trivial and quick and definitely can be evaluated real time when the vector dimension is low; in order to obtain the projected vectors however, we need to perform precalculation on both the environment map and all vertices on the scene (
The only problem left is to find a set of orthonormal functions defined on spherical space.
To understand what is spherical harmonics, we have to take a look at the Associated Legendre Polynomials first. ALP is defined as follows, and can be evaluated incrementally:
!!
stands for a double factorial:
int double_factorial(int n)
{
int ret = 1;
for (int i = n; i >= 1; i -= 2)
{
ret *= i;
}
return ret;
}
It looks crazy, I know. But all of those doesn’t matter. What matters is the properties ALP offers:
This is how it looks like when plotted:
This is very cool, an infinite series of orthonormal functions, defined between
Spherical harmonics are orthonormal functions defined over a unit sphere, which at its core, is ALP with extra steps:
In which
This is how spherical harmonics looks like:
Inigo Quilez wrote a spherical harmonics shader on Shadertoy, which renders the image above using the solved version of spherical harmonics. I have written a SH solver using GLSL as well, but it doesn’t look as fancy. Still though, you can check it out!
With the correct orthonormal function series in hand, we can now perform the projection
to the scene. First, we perform projection
to the environment map, up to the first
The projected result, after reconstructed, is a blurred version of the environment map - that is to be expected, as
The blur won’t affect the end result (that much,) since the whole scene is diffuse to begin with, and information is going to be lost anyway. If a specular object exists in the scene, the reflection will be very blurred, and the loss of information will become very apparent.
Here’s a comparison, with the left image being the original environment map and the right one being the reconstructed version:
Then, we perform projection
vertex-wise to the scene:
What can be projected per vertex? This highly depends on what we want to project. In my implementation (which is given below), I have only projected the cosine term
transport function
. But theoretically, everything within the transport function
can be projected, including visibility test, geometry albedo, etc.
During rendering, a simple dot operation in shader is suffice:
in mat4 shCoeffs;
uniform vec3 scene[16];
out vec4 color;
void main() {
vec3 objColor = vec3(1.0);
vec3 r = vec3(0.0);
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
r += (objColor * shCoeffs[i][j]) * scene[i * 4 + j];
}
}
color = vec4(r, 1.0);
}
mat4
is needed here because OpenGL does not allow custom-length vertex array. As the orthonormal function series has an infinite length, we can choose the first
I have actually made a full-blown application for this. You can check my GitHub repo out here. Both source code and various video demos are included. Feel free to play around! It’s CUDA accelerated, so you will need an NVIDIA GPU with CUDA architecture support and CUDA toolkit installed.
That’s it! We have successfully rendered a PBR scene (without interreflection) in real time, 60fps. However, we have traded off quite a few things during this feat:
There are a lot more to discuss. For example, spherical harmonic-projected coefficient vectors can actually be rotated due to its property. The original paper actually implemented self-occlusion, but I didn’t get that far. There is also neighborhood transfer (so interreflection,) and volumetric transfer (clouds!). Here are some more papers to check out, if you are interested, or if you are still confused after reading this:
Comments