Implementing Needlets, Part 2

In Implementing Needlets we got as far as defining the family of Needlet basis functions, showing what they look like and providing some code for evaluating them. The questions remains, how do I use these to represent a signal on the sphere? Reading the GDC presentation you might get that you need to point Needlets in a bunch of directions and product-integrate them against the signal you want to represent, but which directions and how many points before you get a smooth result?

It’s perfectly possible to take any set of equally spaced directions and use these to approximate the sphere. Take for instance the vertices of tetrahedron.

{{0., 0., 1.}, {-0.471405, -0.816497, -0.333333}, {-0.471405, 
 0.816497, -0.333333}, {0.942809, 0., -0.333333}}


If we point the Needlet N(1,2) along each of these directions and sum the result we get this:


It’s… something, but it doesn’t look like a sphere, not enough points to provide sufficient coverage. If we then take the vertices of an icosahedron normalized onto the unit sphere:

{{0., 0., -1.}, {0., 0., 1.}, {-0.894427, 0., -0.447214}, {0.894427, 
 0., 0.447214}, {0.723607, -0.525731, -0.447214}, {0.723607, 
 0.525731, -0.447214}, {-0.723607, -0.525731, 0.447214}, {-0.723607, 
 0.525731, 0.447214}, {-0.276393, -0.850651, -0.447214}, {-0.276393, 
 0.850651, -0.447214}, {0.276393, -0.850651, 0.447214}, {0.276393, 
 0.850651, 0.447214}}


and do the same operation, something strange happens.


An empty box? Zooming in we find a tiny residual function no larger than 2 \times 10^{-15}, the result of summing the floating point values of each basis:


So given enough density Needlets cancel each other out and sum to zero. What’s going on here?

In general wavelets are designed to represent zero mean functions, i.e. functions that are distributed equally on either side of zero – this is one of the issues with using SH to represent lighting signals where the light probe, visibility and BRDF functions are all positive functions yet we’re using a zero mean system to represent them so we have to carefully avoid reconstruction returning values < 0.

We have successfully projected the unit sphere onto a zero-mean basis and it correctly reconstructs to zero. So somewhere between 4 and 20 points there is a set of points that will accurately reconstruct a band-limited version of our spherical signal using this particular Needlet.

Quadrature and Spherical Designs

Gauss Quadrature is the low cost way to calculate the integral of a function with a known polynomial order and like almost every cool trick in mathematics was invented by Carl Friedrich Gauss [1814]. The idea is that to integrate a function over the range [-1,1], all you need do is sample the polynomial at specific points and make a weighted sum of those samples:

\displaystyle \int_{-1}^{1} f(x) dx = \sum_{i=1}^{n}w_i f(x_i)

More than this, for any n-point Gauss quadrature the samples will correctly integrate any polynomial function f \in \mathcal{P}_{2n+1}. To integrate across any other range you just scale into the [-1,1] domain, calculate the integral and scale the result accordingly. Proofs of this are basic to most undergrad math courses as it’s one of the fundamental relations between Orthogonal bases and function approximation. Gauss quadratures are generated from the zeros of Legendre polynomials, solutions using zeros from Chebychev polynomials are Clenshaw-Curtis quadratures and those using equally-spaced points are Newton-Cotes quadratures (which include the Trapezium and Simpson’s Rules).

Spherical quadrature, sometimes referred to as cubature, has recently gained some interest through it’s use in analyzing directional data and faster integration for Global Illumination.  A weighted quadrature on the \mathbb{S}^2 sphere is known as a spherical design and designs with equal weight on each sample are spherical t-designs. While generating and proving new designs can get involved, using a set of these quadrature points is pretty simple. Lists of designs are available online:

  1. Hardin & Sloane’s page on Spherical t-designs.
  2. Rob Womersley’s page on Cubature on the sphere.

Downloading the data gets you a list of either N pairs of weights and points or, for a t-designs, N points that should be equally weighted \frac{1}{N}. Each design is notated with the number of points and the degree or order of spherical polynomial it will approximate (remembering that a quadrature can accurately integrate degree D and below). Plotting the set of t-designs from Hardin & Sloane’s page we get 112 sets of points here sorted by number of points and colored by the polynomial degree of each set:


You can see that we have runs of designs with the same degree and we’re looking for the minimum number of points that will “cover” the sphere with our Needlet, so which design works the best? Without giving too much away, we can show that N(1,2) is an order-3 function so let’s project it against each of the order-3 designs. Notice how not every set is symmetrical so it’s not just bases opposite each other cancelling out – this cancellation is happening in polynomial frequency space.


The result is that each design can resolve the sphere to within 1 \times 10^{-9}:


So for our purposes we only need to consider the designs with the smallest number of points for each order, lets call it the set of minimum t-designs (shown here for the first 22 orders). There are a few non-intuitive moments in there where order-6 uses less samples than order-5, same with order-8 vs. order-7, so we can pare down the list even further at the low end.


Note that this is a quick result using only the set of t-designs, the designs with uniformly weighted samples. I don’t have any results yet for weighted designs so there could well be point sets out that do the same job for less samples (if you find any tell me and I’ll keep a list updated).

Fractional Needlets

So now we have a way to uncover how many samples each Needlet requires. By matching the Needlet against the list of minimum t-designs we can find the threshold. For the N(j,2) family we find the thresholds are

\begin{tabular}{lll} N(1,2) &= 6 & order 3 \\ N(2,2) &= 24 & order 7 \\ N(3,2) &= 120 & order 18 \\ \\ N(1,3) &= 36 & order 8 \end{tabular}

For fine control Needlets can also be generated with fractional powers, and the function that maps the related bandwidth to polynomial order still needs further investigation (my bet is that it’s related to the highest SH function in the non-zero b() weights). Here’s an animation from N(1,2.0) to N(1,3.0) in steps of +0.05.


Needlets as Frames

So we have Needlet families, we know the sample counts needed to correctly reconstruct a spherical signal, how do we project a spherical function onto Needlets? It works much the same way as SH projection but without the benefits of successive approximation. The level index j in each Needlet family N(j,B) specifies the frequency range that each Needlet is covering and they do so orthogonally like any other wavelet family.


The only difference is that each level requires a different number of points to reliably project the signal. Math libraries that provide Needlet projections are mostly designed to support high-frequency signals and just pick the highest spherical design and project all levels of the Needlet family against that, but for real-time applications we’re looking for the minimal design for each level. Summing the reconstructed signal against successive j indices that were projected using different designs for each level, we can add detail to the approximation in the correct wavelet-like behavior.


Implementing Needlets, Part 2