Error in spherical harmonics implementation

[thetianshuhuang]

Description

The spherical harmonics implementation has some math errors in components (4,-1) and (4,4).

Component (4,-1): a straightforward missing variable.

components[..., 19] = 0.6690465435572892 * y * (7 * zz - 3)
# should be
components[..., 19] = 0.6690465435572892 * y * z * (7 * zz - 3)

Component (4,4): the coefficient appears to be the imaginary coefficient 3/16 * sqrt(35 / (2pi) instead of the real coefficient 3/16 * sqrt(35 / pi).

components[..., 24] = 0.4425326924449826 * (xx * (xx - 3 * yy) - yy * (3 * xx - yy))
# should be
components[..., 24] = 0.6258357354491761 * (xx * (xx - 3 * yy) - yy * (3 * xx - yy))

To Reproduce

This can be checked through numerical integration; I wrote a quick test harness (pardon the modified code -- I found this bug only when working on a different project based on this code, not actually using it), which can be used to verify that the spherical harmonic coefficients are orthonormal:

from matplotlib import pyplot as plt
from jax import numpy as jnp
import numpy as np
from jax import vmap


def spherical_harmonics(dx, harmonics: int = 25):
    """Compute spherical harmonic coefficients.

    Parameters
    ----------
    dx: offset direction from the point to the sensor (on the unit sphere).
    harmonics: number of coefficients to compute.

    References
    ----------
    https://en.wikipedia.org/wiki/Table_of_spherical_harmonics
    https://github.com/nerfstudio-project/nerfstudio:
        nerfstudio/utils/math.py
    """
    x, y, z = dx
    xx = x**2
    yy = y**2
    zz = z**2

    coef = [0.28209479177387814]
    if harmonics >= 4:
        coef += [
            0.4886025119029199 * y,
            0.4886025119029199 * z,
            0.4886025119029199 * x]
    if harmonics >= 9:
        coef += [
            1.0925484305920792 * x * y,
            1.0925484305920792 * y * z,
            0.9461746957575601 * zz - 0.31539156525251999,
            1.0925484305920792 * x * z,
            0.5462742152960396 * (xx - yy)]
    if harmonics >= 16:
        coef += [
            0.5900435899266435 * y * (3 * xx - yy),
            2.890611442640554 * x * y * z,
            0.4570457994644658 * y * (5 * zz - 1),
            0.3731763325901154 * z * (5 * zz - 3),
            0.4570457994644658 * x * (5 * zz - 1),
            1.445305721320277 * z * (xx - yy),
            0.5900435899266435 * x * (xx - 3 * yy)]
    if harmonics >= 25:
        coef += [
            2.5033429417967046 * x * y * (xx - yy),
            1.7701307697799304 * y * z * (3 * xx - yy),
            0.9461746957575601 * x * y * (7 * zz - 1),
            0.6690465435572892 * y * (7 * zz - 3),
            0.10578554691520431 * (35 * zz * zz - 30 * zz + 3),
            0.6690465435572892 * x * z * (7 * zz - 3),
            0.47308734787878004 * (xx - yy) * (7 * zz - 1),
            1.7701307697799304 * x * z * (xx - 3 * yy),
            0.6258357354491761 * (xx * (xx - 3 * yy) - yy * (3 * xx - yy))]

    return jnp.array(coef)


N = 1000000

dx = np.random.normal(size=(N, 3))
dx = dx / np.linalg.norm(dx, axis=1).reshape(-1, 1)
sh = vmap(spherical_harmonics)(dx)
orthogonality = jnp.matmul(sh.T, sh) / N * 4 * jnp.pi

fig, ax = plt.subplots(1, 1, figsize=(8, 8))
cmap = ax.imshow(orthogonality)
fig.colorbar(cmap)

Expected behavior

The spherical harmonics are orthonormal, with L2 norm of 1 (on the diagonal) and inner products zero (off-diagonal).

Instead, we see that component (4,-1) is neither normal and is not orthogonal to all components, and component (4,4) is not normal.

[maturk]

Hey @thetianshuhuang, looks like you indeed found some errors. Maybe you could consider making a pull request with the same comments and fixes to the spherical harmonics components so that this could get more exposure and could be fixed in the code.

[thetianshuhuang]

Currently busy with a deadline, but I'll try to put one together soon!