Understanding spherical harmonics code

Question

I am trying to understand the function computeColorFromSH, which is a part of the differentiable gaussian rasterizer used by the paper 3D Gaussian Splatting for Real-Time Radiance Field Rendering.

I am looking for the formula that the code is based on. Here's an excerpt for the degree > 2 part:

            result = result +
                SH_C3[0] * y * (3.0f * xx - yy) * sh[9] +
                SH_C3[1] * xy * z * sh[10] +
                SH_C3[2] * y * (4.0f * zz - xx - yy) * sh[11] +
                SH_C3[3] * z * (2.0f * zz - 3.0f * xx - 3.0f * yy) * sh[12] +
                SH_C3[4] * x * (4.0f * zz - xx - yy) * sh[13] +
                SH_C3[5] * z * (xx - yy) * sh[14] +
                SH_C3[6] * x * (xx - 3.0f * yy) * sh[15];

It seems to be some kind of polynomial. The Wikipedia article on Spherical Harmonics has only integrals, so I'm guessing that this is a special case and maybe even a well-known one.

Is anyone familiar with this? Where I can read more about it?

Answer 1

You can check this out (the appendix of PlenOctree, Section B.1), since the SH implementation of 3D GS is directly borrowed from PlenOctree. I can walk you through it with an example. First, take a look at the real form of SH (I won't use Latex here since, you know it's painful): Where $P_l^{|m|}(\cdot)$ is the associative Legendre Polynomial:

(From Wikipedia)

So where are those SH_C3? They have some thing to do with: $$ K_{lm} = \sqrt{\frac{2l+1}{4\pi}\frac{(l - |m|)!}{(l+|m|)!}} $$ Ok, so basically, we first get $\theta$ and $\phi$, as the direction $(x, y, z)$ can be easily converted to the spherical coordinates $(\theta, \phi)$. Therefore, the corresponding $\cos(\cdot), \sin(\cdot)$ values can be obtained.

Then we use the formula of associative Legendre Polynomial to get, say $m=0, l = 0$ case: $$ P_0^0=1 $$ We next evaluate SH basis function (given in the first image), for example: $Y_{00}, l=0,m=0$ will be $\sqrt{1 / 4\pi} \approx 0.282095$, which equals to SH_C0 in the code. Other coefficients can be derived similarly: substitute $l, m$ with the order you need, then evaluate the whole formula and store the constant part. With all the precomputed coefficients, now you can multiply your learned SH coefficients with these precomputed ones. Hope this will help.

---

The SH itself is computed in the way describe above. Yet the coefficients for the decomposition in 3D Cartesian space are a bit different:

The image is taken from: Spherical Harmonic Lighting: The Gritty Details, page 22.

The code in diff-gaussian-rasterization:

// cuda_rasterizer/auxiliary.cu
__device__ const float SH_C2[] = {
    1.0925484305920792f,
    -1.0925484305920792f,
    0.31539156525252005f,
    -1.0925484305920792f,
    0.5462742152960396f
};
// cuda_rasterizer/forward.cu
if (deg > 1)
{
    float xx = x * x, yy = y * y, zz = z * z;
    float xy = x * y, yz = y * z, xz = x * z;
    result = result +
        SH_C2[0] * xy * sh[4] +
        SH_C2[1] * yz * sh[5] +
        SH_C2[2] * (2.0f * zz - xx - yy) * sh[6] +
        SH_C2[3] * xz * sh[7] +
        SH_C2[4] * (xx - yy) * sh[8];
}

The code implements exactly the same concept.