Spherical harmonics rotation is one of those problems that you will occasionally run into as a graphics programmer. There has been some recent work, most notably Sparse Zonal Harmonic Factorization for Efficient SH Rotation (Project , PDF) which was presented at Siggraph 2012. In games we usually care about low order SH, especially 3rd order. According to Table 2, it takes 112 multiplications to rotate 3rd order SH by the zxzxz method, and 90 by the improved Zonal Harmonics method. I’ll show you a simple way to do it in 57, and source is provided at the end of the post.
As mentioned in the last post, one way that you can improve the lighting quality is with spherical harmonic ambient occlusion. In my case, I have a skinned mesh with baked spherical harmonics on each vertex. But as the mesh moves the spherical harmonics have to move with it. And since the mesh rotates the SH vector has to rotate too.
For this post, I’m assuming that you understand the basics of spherical harmonics. If you don’t, the best place to start is with Peter-Pike Sloan’s excellent presentation Stupid Spherical Harmonics Tricks (Slides , PDF). It contains everything that you really need to know about spherical harmonics.
As a quick review, spherical harmonics are split into bands. Band 0 has 1 coefficient, band 1 has 3, and band 2 has 5. You can see visualizations of the bands in the image below (which is from the PDF). Also, when we say “Order N” we mean the first N bands. So “Order 3” means “The first 3 bands”, which includes band 0, 1, and 2.
The most important property of Spherical Harmonics is rotational invariance. No matter what direction the light comes from the projected image looks the same. The image below (also taken from Stupid SH Tricks) compares SH to the Ambient Cube from Half-Life 2. With the Ambient Cube a light coming directly along the X/Y/Z axis is much brighter than a light coming from an angle whereas it looks the same from any angle with Spherical Harmonics. The formulas for projecting a normal vector into spherical harmonics are in the appendix of Stupid SH Tricks.
Each of the bands is independent. If we want to rotate 3rd order SH then we need to rotate bands 0, 1, and 2 separately.
Finally, each band can be rotated by a linear transformation. In other words, band N will have 2N+1 coefficients, and we can rotate that band with a square matrix of size 2N+1.
In summary, here are the important properties for rotating spherical harmonics:
- A light direction vector can be projected into spherical harmonics with a simple, closed form solution.
- A direction projected into spherical harmonics looks the same regardless of which direction it comes from.
- We can rotate spherical harmonics with a linear transformation.
- Each band is rotated independently.
Getting to the point, if we have an SH vector and a 3x3 rotation matrix M, how can we rotate the vector? There are many options to do it. We could:
- Rotate around Z and rotate 90 degrees with a closed form solution. So we could decompose our matrix into Euler angles, and multiply by 5 sparse matrices. This is the zxzxz solution.
- Use a Taylor series to approximate the rotation function (as in some PRT work). This option has problems with large angles.
- Recent work (mentioned above) involves factorizing into Zonal Harmonics.
As mentioned before, for 3rd Order the sparse matrix zxzxz solution requires 112 multiplications, the sparse ZH solution requires 90, but we can do it in 57. So what’s the trick?
First, for band 0 we don’t have to do anything because band 0 is just a constant and has no direction. Band 1 is a simple matrix multiplication. So we’ll focus on band 2, but in theory this approach should work for any band. The trick is that rotation followed by projection is the same as projection followed by rotation.
Let’s define a few things.
- x: our SH band 2 that we want to rotate. It has 5 components.
- P: A function which projects a normal vector into band 2. So it takes a 3 component normalized vector as input and outputs a 5 component SH vector.
- M: Our 3x3 rotation matrix. It's the rotation that we want to somehow apply to our SH vector.
- R: The 5x5 (unknown) rotation matrix that want to apply to x.
- N: Some 3D normalized vector.
As mentioned before, if we rotate a vector and then project it into SH we would get the same result by projecting it into SH first and then rotating it. We can describe this algebraically as:
When you think about it this way, solving for R is easy. We can do this same operation for 5 vectors and solve for R.
We can clean this up a little bit by setting a matrix A to the left side, as in:
So this becomes:
And as long as we choose our normal vectors so that A is invertible, we can solve this directly. Which turns into:
That’s it. Since our normal vectors are chosen once we can precalculate invA as A
The final algorithm to rotate our SH vector x by the 3x3 rotation matrix M is:
- Multiply x by invA
- Rotate our 5 pre-chosen normal vectors by M
- Project those rotated vectors into SH which creates a dense 5x5 matrix.
- Multiply our the result of invA*x by the new dense matrix.
If you look at the Zonal Harmonics paper you will see that this algorithm is almost identical. But the advantage is that we can get sparser data. The Zonal Harmonics paper is restricted to finding, you know, Zonal Harmonics. We can just choose 5 vectors out of thin air and it works as long as the projections of those 5 vectors are linearly independent. So we’ll choose these vectors:
And here is what our invA looks like. Inverting a sparse matrix will not necessarily preserve sparsity but in this case it does.
But the really nice thing is that most of these terms end up cancelling out. In the optimized version, the sparse matrix calculation just requires one multiplication. We can divide the whole matrix by k0 and multiply by it at the end which turns most elements of the matrix into ones.
We have to rotate our 5 normal vectors by our 3x3 matrix. In the general case, a 3x3 matrix multiplied with a 3 component vector is 9 multiplies and 6 adds. But two of our vectors don’t need any operations and the 1/sqrt(2) cancels out. So multiplying all 5 normal vectors by M is just 9 adds.
During the projection step the constants in front of each term cancel out too. In the general case projecting into band 2 requires 14 multiplications. Due to canceling of terms, we can skip the 5 multiplications by constants. Projecting and accumulating each vector turns out to require 9 multiplications each.
Then we have a few multiplications on the end and we are finished. All in, rotating band 2 requires 48 multiplications, band 1 requires 9, and band 0 is free. So that’s 57 multiplications for full 3rd Order rotation.
The actual performance gained is heavily dependent on the architecture. On hardware that has a fused multiply add some of these multiplications would have paired with an add anyways so nothing is gained. But the algorithm probably has fewer total instructions regardless.
Finally, here is the source code. Please try it out and let me know if there are any issues. I included the JAMA library for the inverse. I haven’t tested with larger orders but it the same approach should work. Although for really high SH bands the zxzxz approach probably wins out since the dense matrix multiplication would start to dominate.
Source: ShRotation.zip
comments powered by Disqus