Spherical Discrete Sampling Techniques

• Spherical discrete sampling is the method of replacing continuous measures on the sphere with finite point sets and quadrature rules that preserve key invariants.
• These techniques optimize criteria like equal-area coverage, minimal angular separation, and exact harmonic recovery, adapting to various problem-specific requirements.
• They enable applications in diffusion MRI, rendering, and graph-based analyses through innovations such as optimal-dimensionality designs and area-preserving mappings.

Spherical discrete sampling is the construction of finite point sets, quadrature rules, graph discretizations, mixture expansions, or deterministic warps that replace continuous measures on $S^2$ or on spherical subsets by structured discrete objects while preserving task-specific invariants. In the literature, the target may be exact recovery of band-limited spherical harmonic coefficients, nearly equal-area coverage, maximal minimal angular separation, efficient sampling of rotationally invariant operators, low-variance solid-angle integration on spherical patches, or exact simulation of spherical probability laws through countable mixtures (McEwen et al., 2011, Elahi et al., 2018, Perraudin et al., 2018, Guillén et al., 2018, Cheng et al., 2017, Baringhaus et al., 2023).

1. Formal scope and evaluation criteria

A central formal setting treats a scalar field on the sphere as a spherical-harmonic expansion

$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$

with $L^2$ harmonic degrees of freedom when the field is band-limited at $L$ (McEwen et al., 2011). For spin-$s$ fields, the natural basis is ${}_sY_{\ell m}$, and the harmonic dimensionality becomes

$N_o=L^2-s^2,$

while antipodal scalar fields retain only even degrees and therefore have

$N_0=\frac{L(L+1)}{2}$

degrees of freedom (Elahi et al., 2018, Bates et al., 2015). In this regime, spherical discrete sampling means choosing spatial samples and a transform so that these coefficients can be recovered accurately and stably.

A second formal setting is geometric rather than harmonic. For a finite set $\{\mathbf u_i\}_{i=1}^K\subset S^2$, one may optimize the covering radius

$$d\bigl(\{\mathbf u_\ell\}_{\ell=1}^K\bigr)=\min_{i\neq j}\arccos\bigl|\mathbf u_i^\top \mathbf u_j\bigr|,$$

which directly measures minimal angular separation under antipodal symmetry (Cheng et al., 2017). Other geometric criteria appearing in the literature include minimum geodesic distance, mesh norm, and mesh ratio for comparing covering and packing quality of spherical point sets (Elahi et al., 2018). These criteria are distinct from exact quadrature or exact harmonic inversion: a configuration that is optimal for spherical code packing need not be optimal for spherical harmonic transforms, and vice versa.

A third setting discretizes not the whole sphere but a spherical subset or a spherical distribution. In rendering, the domain may be a spherical ellipse or spherical cap, and the objective is area preservation with respect to solid angle rather than global harmonic exactness (Guillén et al., 2018, Dupuy et al., 2023). In probability, the goal may be to express a spherical law as a countable mixture over simpler spherical components, reducing continuous sampling to discrete sampling over an index set followed by conditional sampling on the sphere (Baringhaus et al., 2023). This suggests that “uniformity” on the sphere is not a single notion but a family of task-dependent invariance requirements.

2. Exact sampling theorems and optimal-dimensionality constructions

Classical spherical sampling theorems seek exact reconstruction of band-limited functions from finitely many samples. An influential equiangular construction associates the sphere with the torus through a periodic extension and yields exact scalar and spin transforms with sample count

$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$0

while maintaining $$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$1 transform complexity and avoiding precomputation (McEwen et al., 2011). In the same account, this count is less than half that of Driscoll–Healy-type equiangular schemes and asymptotically identical, but smaller, than Gauss–Legendre sampling (McEwen et al., 2011). The underlying idea is that iso-latitude sampling and a periodic extension in colatitude convert the spherical transform into FFT-compatible Fourier structure plus dense but structured Wigner-$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$2 algebra.

Optimal-dimensionality designs sharpen this principle by matching the number of spatial samples to the harmonic degrees of freedom. For antipodally symmetric diffusion MRI signals, an iso-latitude construction uses exactly

$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$3

samples, equal to the number of nonzero even-degree coefficients, and supports a corresponding spherical harmonic transform with near machine-precision accuracy for practically relevant band-limits (Bates et al., 2015). For spin-$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$4 functions, an iso-latitude, variable-longitude grid

$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$5

uses exactly

$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$6

samples and is accompanied by an $$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$7 spin-SHT whose co-latitudes are chosen by condition-number minimization of Wigner-$$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$8-based linear systems; a multi-pass refinement further reduces residual error (Elahi et al., 2018).

These results correct a common assumption that exact spherical harmonic analysis necessarily requires approximately $$f(\theta,\phi)=\sum_{\ell=0}^{L-1}\sum_{m=-\ell}^{\ell} f_{\ell m}Y_{\ell m}(\theta,\phi),$$9 or more samples. The optimal-dimensionality constructions show that the oversampling of classical schemes is not information-theoretically necessary, although it may simplify quadrature or implementation (Elahi et al., 2018, Bates et al., 2015). Conversely, they also show that exactness is tied to highly structured sample placement and numerically stable inversion, not merely to point count.

3. Grid geometries, equal-area discretizations, and graph representations

When the priority is scalable computation on full or partial sky maps, spherical discrete sampling is often realized as a pixelization. HEALPix partitions the sphere into $L^2$0 base pixels and $L^2$1 equal-area pixels, while preserving iso-latitude rings and an exact four-child hierarchical refinement (Perraudin et al., 2018). These three properties—equal area, iso-latitude, and hierarchy—are emphasized as jointly crucial for multiscale processing and fast transforms, and the same discretization underlies graph-based spherical CNNs, where one node is assigned to each pixel and neighboring pixels define a sparse weighted graph (Perraudin et al., 2018). In DeepSphere, the graph Laplacian

$L^2$2

acts as a discrete Laplace–Beltrami operator; Chebyshev polynomials of $L^2$3 yield localized, radial graph filters with $L^2$4 convolution cost and approximate rotation equivariance (Perraudin et al., 2018). The approximation is not exact because HEALPix is not graph-regular: most pixels have eight neighbors, but $L^2$5 pixels have only seven neighbors, and edge lengths vary slightly (Perraudin et al., 2018).

Equal-area sampling also appears in differentiable rendering. UniTriSplat reformulates 3D Gaussian splatting on the unit sphere by HEALPix discretization, with each pixel covering

$L^2$6

The grid resolution is chosen to match input angular resolution through

$L^2$7

then rounded to the nearest power of two (Zhu et al., 29 Jun 2026). The result is uniform solid-angle weighting across perspective, fisheye, and omnidirectional cameras, together with HEALPix-aware SSIM and spherical gradient propagation in arc-length coordinates (Zhu et al., 29 Jun 2026).

A different design point is SymPix, which keeps Gauss–Legendre ring latitudes but varies the number of pixels per ring by bands of $L^2$8 rings, with neighboring band ratios restricted to

$L^2$9

This sacrifices exact equal area but introduces many repeated local neighborhoods, allowing rotationally invariant local kernels to be sampled far more efficiently than on HEALPix (Seljebotn et al., 2015). For $L$0, the reported average speed-ups are $L$1 for $L$2 and $L$3 for $L$4 relative to HEALPix when constructing a representative preconditioner (Seljebotn et al., 2015). A plausible implication is that equal-area sampling and operator-sampling efficiency are separate optimization axes: SymPix optimizes symmetry reuse, whereas HEALPix optimizes area equality and hierarchy.

Graph constructions can also be made sampling-agnostic. Interpolated SelectionConv represents spherical or surface data as a graph of sampled points with directional adjacency matrices defined in local tangent frames, and explicitly supports equirectangular, Fibonacci spiral, icosphere, layering, and random samplings (Hart et al., 2022). Its ablation study reports that layering clustering is consistently best across samplings, and that angle-based interpolation generally outperforms barycentric interpolation because the latter introduces more smoothing (Hart et al., 2022). This decouples convolution from any single spherical grid and shifts the design problem to neighborhood geometry, interpolation, and multiresolution clustering.

4. Area-preserving warps and sampling of spherical subsets

A major branch of spherical discrete sampling does not seek a global grid at all. Instead, it constructs a deterministic map

$L$5

that transports 2D sample sets—random, stratified, or low-discrepancy—onto a spherical subset $L$6 while preserving area in solid angle. For direct illumination from disk lights, the relevant subset is a spherical ellipse: the image on the unit sphere of the disk as seen from a shading point (Guillén et al., 2018). The paper develops two exact area-preserving parameterizations from the unit square to this spherical ellipse, a parallel map and a radial map, together with a low-distortion radial variant, all based on a generalized Archimedes hat-box theorem and slice-area formulas involving incomplete elliptic integrals (Guillén et al., 2018). The result is uniform solid-angle sampling of disk lights without rejection, lower variance than area sampling, and practical integration into Mitsuba and Arnold through analytic or tabulated inversion (Guillén et al., 2018).

This framework is important because source-area uniformity and solid-angle uniformity are not equivalent. If one samples points uniformly on the emitting disk and converts them to directions, the directional density varies strongly over the spherical ellipse, especially when the disk is near the shading point (Guillén et al., 2018). By contrast, a constant directional density $L$7 over the spherical ellipse reduces variation in $L$8 for the Monte Carlo estimator and therefore reduces variance (Guillén et al., 2018). The misconception that planar sampling of a light source is already “uniform on the sphere” is therefore false in this setting.

A related construction appears in visible-normal sampling for GGX microfacet models. After transforming the anisotropic GGX problem to a unit hemisphere, visible-normal sampling can be reduced to uniform sampling over a spherical cap

$L$9

with cap density

$s$0

A sampled cap direction is then combined with the incident direction by a half-vector construction to produce a visible normal (Dupuy et al., 2023). This spherical-cap formulation yields the same VNDF as the previous method of Heitz while simplifying the implementation and producing systematic speed-ups in CPU and GPU benchmarks (Dupuy et al., 2023). Together, these rendering examples show that local spherical warps are a core part of spherical discrete sampling whenever the integration domain is a structured subset of $s$1 rather than the entire sphere.

5. Sparse, probabilistic, and constrained spherical sampling

Sparse spherical sampling studies signals that are not band-limited but admit a finite-parametric representation. A key example is a sum of $s$2 Diracs on the sphere convolved with a band-limited kernel. By rewriting low-pass spherical harmonic coefficients into a structured data matrix whose columns are exponential sums in $s$3, the paper generalizes annihilating filters to the sphere and proves that $s$4 spikes can be reconstructed from

$s$5

spatial samples of the low-pass observation (Dokmanic et al., 2015). This improves previous spherical FRI requirements by a factor of four for large $s$6, requires no separation condition between Diracs, and supports applications to diffusion source localization, shot-noise removal, and spherical microphone arrays (Dokmanic et al., 2015). Here spherical discrete sampling is organized around finite rate of innovation rather than harmonic exactness.

Probability theory provides another discrete layer: a continuous spherical law may admit a countable mixture representation

$s$7

so that sampling reduces to drawing a discrete index $s$8 and then sampling from a simpler spherical component $s$9 (Baringhaus et al., 2023). The paper derives exact discrete mixtures for the von Mises–Fisher, Watson, angular Gaussian, and spherical Cauchy families, using bases such as ${}_sY_{\ell m}$0, ${}_sY_{\ell m}$1, or an ultraspherical family ${}_sY_{\ell m}$2 built from normalized Gegenbauer polynomials (Baringhaus et al., 2023). The mixing laws include confluent hypergeometric-series, negative-binomial, hypergeometric-series, discrete parabolic-cylinder, and generalized positive Skellam distributions (Baringhaus et al., 2023). This makes “spherical discrete sampling” literal: the continuous spherical law is represented as a discrete probability distribution over spherical component families.

A related geometric strand uses the sphere as an auxiliary sampling manifold for constrained Euclidean domains. Spherical Augmentation maps, for example, the unit ball ${}_sY_{\ell m}$3 to the sphere ${}_sY_{\ell m}$4 by

${}_sY_{\ell m}$5

then runs HMC or LMC on the sphere so that proposals automatically satisfy the original constraints when mapped back (Lan et al., 2015). The same framework treats ${}_sY_{\ell m}$6-balls, boxes, quadratic constraints, and the simplex; for simplex variables in LDA, the Fisher metric becomes, up to scale, the canonical spherical metric after the square-root map (Lan et al., 2015). Although the original domain is not ${}_sY_{\ell m}$7, the method extends the conceptual reach of spherical sampling as a way of removing boundaries through spherical geometry.

6. Angular optimality, discrete operators, and broader analogues

In diffusion MRI, spherical discrete sampling is often judged directly by angular resolution. The spherical code formulation defines the local nearest-neighbor angle by

${}_sY_{\ell m}$8

and seeks to maximize the global minimum over all pairs on a single shell, or a weighted combination of within-shell and across-shell covering radii for multiple shells (Cheng et al., 2017). This differs from Electrostatic Energy Minimization, which minimizes a sum of repulsive pairwise energies and only approaches covering-radius maximization in the limit ${}_sY_{\ell m}$9 (Cheng et al., 2017). The paper introduces incremental SC, IMOC, 1-Opt, MILP, and CNLO algorithms for continuous design, discrete subsampling, and acquisition ordering, and reports that SC methods obtain larger angular separation and better rotational invariance than EEM and GEEM (Cheng et al., 2017). The contrast is instructive: maximizing minimal angle and minimizing electrostatic energy are related but not identical objectives.

Discrete spherical means provide yet another interpretation. For directional derivatives on the circle or sphere, continuous rotationally invariant averages can be replaced by finite weighted sums over sampled directions. In two dimensions, weights $N_o=L^2-s^2,$0 and directions $N_o=L^2-s^2,$1 satisfying

$N_o=L^2-s^2,$2

yield discrete identities for the Laplacian and for $N_o=L^2-s^2,$3, while in arbitrary dimension the construction is generalized via harmonic Veronese maps and Minkowski’s existence theorem (Belyaev et al., 2011). This produces finite-difference Laplacians and quasi-Laplacians with good rotation-invariance properties, including familiar $N_o=L^2-s^2,$4 stencils as special cases (Belyaev et al., 2011). The broader principle is that spherical sampling sets can be designed by forcing discrete annihilation of low-degree harmonic moments.

The topic also has a lattice analogue. On $N_o=L^2-s^2,$5, one may average over discrete spheres

$N_o=L^2-s^2,$6

and study the lacunary maximal operator formed from these spherical averages. For $N_o=L^2-s^2,$7, and for $N_o=L^2-s^2,$8 with radii $N_o=L^2-s^2,$9, the discrete lacunary spherical maximal function is bounded on $N_0=\frac{L(L+1)}{2}$0 for all

$N_0=\frac{L(L+1)}{2}$1

using Magyar’s decomposition and a Kloosterman refinement (Anderson et al., 2020). This is not sampling on $N_0=\frac{L(L+1)}{2}$2 itself, but it shows that the phrase “spherical discrete sampling” can also refer to spheres defined intrinsically inside a discrete ambient space.

Across these strands, a recurring lesson is that no single discretization dominates all others. Exact harmonic inversion favors carefully designed iso-latitude schemes; equal solid-angle coverage favors HEALPix-like constructions; efficient localized operator sampling favors SymPix-like symmetries; rendering on spherical subsets favors area-preserving warps; sparse recovery favors annihilating structures; and distributional sampling may favor discrete mixture bases (McEwen et al., 2011, Seljebotn et al., 2015, Guillén et al., 2018, Dokmanic et al., 2015, Baringhaus et al., 2023). Spherical discrete sampling is therefore best understood as a family of constructions that tailor finite spherical representations to the invariants, transforms, and operators of the problem at hand.