Source-linked AI summary
A novel sampling theorem on the sphere
J. D. McEwen, Y. Wiaux
TL;DR
Exact spherical sampling requires O(L^2) measurements, while practical transforms must also control computation and precomputation. The paper uses a periodic extension connecting the sphere to the torus and FFT-based algorithms to construct an exact sampling theorem. It reports fewer samples than competing exact schemes, fast O(L^3) transforms without precomputation, and applicability to both scalar and spin functions.
Problem
Exact representation of band-limited spherical signals requires O(L^2) samples, while spherical harmonic transforms must also be computationally efficient and avoid costly precomputation.
Method
The paper associates the sphere with the torus through a periodic extension and exploits FFTs to construct exact forward and inverse spherical harmonic transforms.
Results
The sampling theorem uses fewer samples than other equiangular theorems, has O(L^3) transform complexity with a reduced constant prefactor, requires no precomputation, and applies to scalar and spin functions.
Takeaways & Limitations
The theorem may benefit high-band-limit CMB analysis, low-band-limit diffusion MRI acquisition, and compressive-sampling reconstruction through reduced sampling or dimensionality.
Abstract
from arXiv · showhide
We develop a novel sampling theorem on the sphere and corresponding fast algorithms by associating the sphere with the torus through a periodic extension. The fundamental property of any sampling theorem is the number of samples required to represent a band-limited signal. To represent exactly a signal on the sphere band-limited at L, all sampling theorems on the sphere require O(L^2) samples. However, our sampling theorem requires less than half the number of samples of other equiangular sampling theorems on the sphere and an asymptotically identical, but smaller, number of samples than the Gauss-Legendre sampling theorem. The complexity of our algorithms scale as O(L^3), however, the continual use of fast Fourier transforms reduces the constant prefactor associated with the asymptotic scaling considerably, resulting in algorithms that are fast. Furthermore, we do not require any precomputation and our algorithms apply to both scalar and spin functions on the sphere without any change in computational complexity or computation time. We make our implementation of these algorithms available publicly and perform numerical experiments demonstrating their speed and accuracy up to very high band-limits. Finally, we highlight the advantages of our sampling theorem in the context of potential applications, notably in the field of compressive sampling.
I. INTRODUCTION
The paper motivates exact and efficient spherical harmonic analysis for data on the sphere, then introduces a sampling theorem and fast algorithms evaluated by sample count, complexity, speed, and precomputation.
- Spherical data arise in computer graphics, planetary science, geophysics, quantum chemistry, and astrophysics, where harmonic analysis can provide insight.
- Scalar and spin spherical harmonic transforms are important for large datasets, including temperature and polarisation observations of the cosmic microwave background.
- Approximate pixelisation-based transforms can offer practical criteria such as equal-area pixels, but do not correspond to exact sampling theorems.
- Exact sampling theorems are assessed by sample count, computational complexity, speed, and precomputation requirements.
- The paper presents a sampling theorem and algorithms applying to scalar and spin functions without changes in asymptotic complexity or computation time.
- Existing exact approaches include Gauss-Legendre sampling with approximately 2L^2 samples and equiangular schemes requiring approximately 4L^2 samples.
- Methods below O(L^3) generally require O(L^3) precomputation and storage, reaching 1.2GB at L = 1024 and approximately 77GB at L = 4096.
III. HARMONIC ANALYSIS ON THE SPHERE
This section establishes scalar spherical harmonics as the orthogonal basis for square-integrable functions on the sphere and expresses functions through their harmonic coefficients.
- Spherical coordinates use colatitude θ ∈ [0, π] and longitude ϕ ∈ [0, 2π), with invariant measure dΩ(θ, ϕ) = sin θ dθ dϕ.
- Scalar spherical harmonics form the canonical orthogonal basis for square-integrable scalar functions on the sphere.
- The adopted Condon-Shortley phase convention ensures the stated conjugate symmetry relation for scalar spherical harmonics.
- The basis satisfies orthogonality and completeness relations involving Kronecker delta symbols and the Dirac delta function.
- Any square-integrable scalar function can be represented by a spherical harmonic expansion with coefficients obtained by projection onto each basis function.
B. Spin spherical harmonics
Spin spherical harmonics generalize scalar harmonics to functions transforming under local tangent-plane rotations, while retaining an orthogonal basis and harmonic expansion.
- A spin function transforms under a rotation by χ in the tangent plane at a sphere point, rather than under a global rotation of the sphere.
- Spin spherical harmonics form an orthogonal basis for spin-s functions when |s| ≤ ℓ.
- Square-integrable spin functions admit harmonic expansions whose coefficients are projections onto the corresponding spin spherical harmonics.
- Spin harmonics and transforms reduce to the scalar case when s = 0.
- Spin raising and lowering operators generate spin s ± 1 functions from spin-s functions and can relate spin functions to scalar functions.
C. Wigner functions
The paper represents spin spherical harmonics through Wigner functions and uses their Fourier-series structure as the foundation for extending spherical transforms to the torus.
- Wigner functions form an orthogonal basis for square-integrable functions on SO(3) and are parameterized by Euler angles.
- The paper adopts the zyz Euler convention for rotations about the z, y, and z axes by γ, β, and α, respectively.
- The Wigner d-functions use real polar d-functions and admit recursion formulae for rapid computation.
- Spin spherical harmonics are defined through Wigner functions, allowing standard Wigner-function decompositions to be applied to spherical harmonics.
- The Fourier-series representation of the d-functions enables the spherical harmonic expansion to be written as a Fourier series on an appropriately extended two-torus.
IV. FAST SPHERICAL HARMONIC TRANSFORM
The paper reformulates scalar and spin spherical harmonic transforms through a periodic extension from the sphere to the torus, enabling FFT-based exact algorithms and a corresponding sampling theorem.
- IV. FAST SPHERICAL HARMONIC TRANSFORM: The resulting construction provides exact forward and inverse transforms and a sampling theorem on the sphere.
- IV. FAST SPHERICAL HARMONIC TRANSFORM: The harmonic formulation covers spin functions band-limited at L, with summations truncated at L−1 and coefficients constrained by |m| ≤ ℓ.
- IV. FAST SPHERICAL HARMONIC TRANSFORM: The forward and inverse transforms are derived in forms that expose similarities with Fourier series representations.
- IV. FAST SPHERICAL HARMONIC TRANSFORM: The approach extends the colatitude domain periodically so Fourier transforms can be applied to spherical harmonic expressions.The extension preserves the representation’s required symmetries.
B. Periodic extension
The periodic extension maps the colatitude domain to [0, 2π) using symmetry-dependent reflection, while sampling choices preserve continuity and align harmonic indices with Fourier modes.
- B. Periodic extension: The colatitude domain is extended to [0, 2π) by reflecting samples with a factor (−1)^(m+s) on the second half.The imposed symmetry is chosen to match the symmetry of the Fourier-transformed spherical function.
- B. Periodic extension: The extension is applied to sG_m(θ), differing from earlier torus-based approaches that used alternative even, odd, or longitude-shifted extensions.
- B. Periodic extension: Odd sample counts in θ and ϕ associate the extended-domain samples directly with harmonic indices m and m′.
- B. Periodic extension: At the poles, sF_m(θ⋆) vanishes when m + s is odd, preventing discontinuities in the periodic extension.
- B. Periodic extension: The sampling theorem uses N_MW = (L − 1)(2L − 1) + 1 ∼ 2L^2 samples on the sphere.
D. Forward transform
The forward transform uses implicit quadrature and FFTs to compute spherical harmonic coefficients efficiently, retaining O(L^3) overall complexity while avoiding precomputation.
- D. Forward transform: The forward algorithm uses implicit quadrature rules for the exact evaluation of the transform integrals.
- D. Forward transform: FFT-based convolution reduces the computation of sG_mm′ from O(L^3) directly to O(L^2 log^2 L).Zero-padding places aliasing outside the coefficient range of interest.
- D. Forward transform: The overall forward algorithm remains O(L^3), dominated by computing the spin spherical harmonic coefficients.
- D. Forward transform: Conceptually, the algorithm Fourier-transforms in ϕ, extends and upsamples in θ, applies reflected weights, and Fourier-transforms in θ.
- D. Forward transform: Continual FFT use lowers the constant prefactor, producing rapid harmonic-coefficient computation without precomputation.
- D. Forward transform: The algorithm applies to scalar and spin functions with unchanged computational complexity and computation time.
E. Inverse transform
The inverse transform computes intermediate Fourier coefficients, reconstructs the function on an extended domain using FFTs, and discards samples outside the original θ domain.
- E. Inverse transform: The inverse algorithm computes sF_mm′ directly from the harmonic coefficients and exploits symmetry to halve the number of computations.
- E. Inverse transform: Its dominant computation has O(L^3) complexity, while subsequent FFT evaluations have O(L^2 log^2 L) complexity.
- E. Inverse transform: The inverse procedure uses an inverse Fourier transform on the extended domain and returns the function after discarding samples in θ ∈ (π, 2π).
- E. Inverse transform: Because the forward and inverse transforms are theoretically exact, the construction constitutes a novel sampling theorem on the sphere.It requires fewer samples than other equiangular schemes and an asymptotically identical but smaller number than Gauss–Legendre sampling.
F. Quadrature
The sampling theorem yields an explicit quadrature rule whose weights are constructed from periodically extended samples and folding, reducing the sampling required for selected coefficients.
- Quadrature construction: The explicit quadrature rule integrates a band-limited function exactly using weights q(θt).The integration corresponds to computing a spherical harmonic coefficient.
- Quadrature construction: Aliasing tolerance reduces the required azimuthal samples from 2L−1 to L while retaining L samples in θ.Only the m′ = 0 coefficient is required in the relevant calculation.
- Weight construction: The weights v(θt) sample sin(θ) on [0, π) and zero on [π, 2π), while q(θt) folds the latter contributions onto [0, π].Both weight sets are compared with sin(θ) in Fig. 1.
G. Real signals
Reality and conjugate symmetries reduce the computational cost of the spherical harmonic transforms, while the evaluation compares sampling requirements and implementation choices across methods.
- Real signals: Reality conditions impose conjugate symmetry on spin spherical harmonic coefficients and reduce to the standard scalar reality condition when s = 0.
- Real signals: Symmetry exploitation reduces the computational cost of both forward and inverse algorithms by approximately two for real spin s = 0 signals.
- Evaluation setup: The implementation uses double-precision C code with FFTW and provides a MATLAB interface, while the Gauss-Legendre comparison algorithm has complexity O(L^3).
- Evaluation setup: The evaluation examines exact sample counts, precomputation, Wigner d-function recursions, numerical accuracy, computation time, and potential applications.
- Evaluation setup: Gauss-Legendre quadrature with L samples in θ computes exact spherical harmonic transforms for scalar and spin functions band-limited at L.The polynomial integrand has maximum degree 2L − 2.
A. Sampling
The proposed sampling theorem uses fewer samples than competing sphere sampling schemes while retaining O(L^3) transform complexity and avoiding required precomputation. Its algorithms are accurate and fast, support scalar and spin transforms, and remain stable to higher band-limits than the compared alternatives.
- Sampling requirements: N_MW ∼ 2L^2 samples are required, versus N_DH ∼ 4L^2 for Driscoll & Healy; the proposed theorem therefore halves the canonical equiangular sampling count.It also requires 3(L − 1) fewer samples than Gauss-Legendre sampling.
- Experimental setup: The numerical experiments use random band-limited signals and measure maximum absolute coefficient error after inverse-then-forward transform round trips.Experiments compare the proposed algorithms with optimized Gauss-Legendre and semi-naive Driscoll & Healy implementations without precomputation.
- Numerical accuracy: At moderate band-limits, all methods achieve errors near machine precision, increasing approximately linearly with band-limit.The corresponding error is measured over repeated random test signals.
- Numerical accuracy: The Gauss-Legendre and Driscoll & Healy implementations become unstable between L = 1024 and L = 2048 because of the pointwise three-term Wigner recursion.The proposed algorithms have superior numerical accuracy because they permit more flexible Wigner-function recursions.
- Computation time: All implementations scale as O(L^3); the semi-naive algorithm is slightly faster, while the proposed algorithm is slightly faster than optimized Gauss-Legendre.Using Trapani & Navaza recursion makes the proposed implementation approximately 20% faster than using Risbo’s recursion.
E. Applications
The sampling theorem is positioned for cosmology, neuroscience, and compressive sampling. Its lower sample count and fast transforms support high-resolution CMB analysis, lower-cost diffusion MRI acquisition, and potentially improved spherical compressive-sampling reconstruction.
- Cosmology: For high-resolution CMB observations, the theorem provides exact fast spherical harmonic transforms for both temperature and polarisation at identical cost.Temperature corresponds to scalar functions, while polarisation gives rise to spin ±2 functions.
- Neuroscience: Fewer required samples can reduce the cost of exactly measuring band-limited functions when acquiring each sample is expensive, including diffusion MRI.The paper identifies diffusion MRI as an application where acquisition cost is important.
- Compressive sampling: Compressive sampling requires fewer measurements for sparse or compressible signals, with the measurement-to-dimensionality ratio scaling linearly with sparsity.The paper expects its reduced spatial-domain dimensionality to enhance compressive-sampling reconstruction on the sphere.
- Applications: The theorem uses fewer samples than other equiangular schemes and a slightly smaller asymptotic count than Gauss-Legendre sampling, with O(L^3) fast algorithms.The conclusion connects these properties to potential application benefits.
- Applications: At L ∼ 4096, the theorem is described as stable and very accurate for CMB transforms, while at L ∼ 10 its lower sample count may reduce diffusion MRI acquisition cost.These are proposed practical benefits at high and low band-limits, respectively.