Spherical Feature Interpolation

• Spherical feature interpolation reconstructs scalar and feature fields on the sphere or hypersphere while preserving its intrinsic geometric structure.
• Key methods include chart-based local interpolation, spectral techniques with spherical harmonics, and geodesic SLERP-based approaches.
• Applications span graphics, acoustics, and point-cloud processing, enhancing smoothness, derivative quality, and computational efficiency.

Spherical feature interpolation is the reconstruction of values, derivatives, or learned representations indexed by directions on the sphere, or, in a broader hyperspherical sense, by unit-norm feature vectors. In current literature, the subject spans scalar fields $r:S^2\to\mathbb R$ stored on cubemaps for graphics, spectral interpolation of fields on spherical grids, geodesic interpolation of sampled curves on $\mathbb S^2$ via repeated spherical linear interpolation, kernel and nearest-neighbor schemes for scattered spherical data, and learned interpolation operators for point-cloud and HRTF features (Abouagour et al., 23 Feb 2026, Chiwere et al., 2024, Fong et al., 2022, Chen et al., 2023). Across these settings, the central difficulty is to preserve intrinsic spherical geometry while controlling smoothness, pole and seam behavior, rotational structure, derivative quality, and runtime cost.

1. Scope and problem setting

A basic formulation treats the target as a spherical signal sampled at finitely many directions. In graphics, one representative model is a scalar function

$r:S^2\to\mathbb R$

sampled on a regular cubemap parameterization, with each face carrying local chart coordinates $(u,v)\in[0,1]^2$ (Abouagour et al., 23 Feb 2026). In acoustics, one frequency slice of a head-related transfer function is treated as a real-valued function on the sphere, written as

$HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$

with $\Omega\in S^2$ and frequency index $l$ (Chen et al., 2023). In conditional neural-process formulations, interpolation is posed on unordered context-target sets $\{(x_c,y_c)\}$ with $x_c\in S^2$ and targets queried at unseen spherical locations (Thuillier et al., 2023).

The term also appears in settings where the “feature” is not a scalar field on $S^2$ but a unit-length vector in a high-dimensional embedding space. In that case, the relevant geometry is the hypersphere, and interpolation proceeds along great-circle arcs rather than Euclidean chords (Zhong et al., 2024). This broader usage is technically distinct from interpolation of sampled spherical functions, but it shares the same geometric principle: Euclidean linearization is replaced by an interpolation rule that respects the manifold structure.

Several failure modes recur across applications. Tensor-product grids in spherical coordinates can introduce artificial boundaries at the poles (Chiwere et al., 2024). Bilinear look-up-table sampling is fast but produces faceting, while finite-difference normals require extra samples and introduce noise (Abouagour et al., 23 Feb 2026). In point-cloud processing, traditional grid-shaped 3D convolution operators are difficult to apply directly because of the unstructured property of raw point clouds (Wang et al., 2020). In HRTF interpolation, image-super-resolution formulations can neglect spatial acoustic features (Chen et al., 2023). These observations motivate geometry-aware interpolation operators rather than planar or purely Cartesian substitutes.

2. Geometric and analytic formulations

Three formulations dominate the literature: chart-based local interpolation, global spectral interpolation, and geodesic interpolation.

In chart-based formulations, the sphere is covered by parameterized patches, and interpolation is performed in local coordinates before being mapped back to the tangent geometry of $\mathbb S^2$0. Spherical Hermite Maps use a padded cubemap, with each face stored with a one-texel gutter so that adjacent-face values and derivatives are available at seams (Abouagour et al., 23 Feb 2026). This setup supports local tensor-product cubic reconstruction while retaining access to chart derivatives.

In spectral formulations, the spherical signal is expanded in spherical harmonics or related trigonometric bases. For example, HRTF interpolation uses the truncated expansion

$\mathbb S^2$1

with coefficients estimated by least squares from sparse measurements (Chen et al., 2023). Closely related FFT-based and DFS-based methods treat spherical interpolation through trigonometric structure. The Double Fourier Sphere method transforms the sphere into a doubly periodic domain and thereby avoids artificial boundaries at the poles (Chiwere et al., 2024). Uniform-grid FFT interpolants and DFFT-type constructions likewise exploit separability in spherical-polar variables and achieve $\mathbb S^2$2 complexity on regular grids (Dominguez et al., 2014, Dominguez et al., 2012).

In geodesic formulations, interpolation is performed intrinsically on the sphere through spherical linear interpolation. For unit vectors or unit quaternions, the standard SLERP formula is

$\mathbb S^2$3

which traces the shorter great-circle arc between endpoints (Zhong et al., 2024). Higher-order constructions such as SIDER replace affine Neville blends by recursive SLERP blends on local stencils (Leung, 11 Jun 2026).

This diversity suggests that “spherical feature interpolation” is not a single operator class but a geometric design space. The appropriate formulation depends on whether the data are grid-based or scattered, whether derivatives are required, whether the task is local or global, and whether the target object is a scalar field, a curve on $\mathbb S^2$4, or a hyperspherical feature embedding.

3. Major methodological families

The literature can be organized into a small number of recurrent method families.

Family	Core construction	Representative papers
Spectral and barycentric	Spherical harmonics, DFS, FFT, parity-modified Fourier bases	(Chiwere et al., 2024, Dominguez et al., 2014, Dominguez et al., 2012, Erb, 2018)
Local kernel and scattered-data methods	ZBFs, Shepard partition of unity, zonal kernels, k-nearest neighbors	(Rossi, 2014, Sun et al., 2024, Sun et al., 14 Oct 2025, Trempe, 2019)
Constrained polynomial approximation	Interpolation-regression with Vandermonde and KKT systems	(Dell'Accio et al., 13 May 2026)
Geodesic high-order interpolation	SLERP, SIDER-$\mathbb S^2$5, SENO	(Fong et al., 2022, Leung, 11 Jun 2026, Leung, 11 Jun 2026)
Learned spherical operators	Spherical interpolated convolution, spherical CNNs, spherical neural processes	(Wang et al., 2020, Chen et al., 2023, Thuillier et al., 2023)
Derivative-augmented LUTs	Bicubic-Hermite reconstruction on padded cubemaps	(Abouagour et al., 23 Feb 2026)

Spectral methods emphasize global smoothness, exactness properties, and fast transforms. The DFS barycentric formulas rely on precomputed weights that depend on the grid structure and not the data itself, and for standard tensor-product grids they exhibit exponential convergence when approximating smooth functions (Chiwere et al., 2024). Uniform-grid FFT interpolants provide spectrally accurate convergence in Sobolev norms and constructive $\mathbb S^2$6 algorithms (Dominguez et al., 2014). Lissajous-based schemes build interpolation spaces from a parity-modified double Fourier basis and report logarithmic growth of the numerical condition number (Erb, 2018).

Local kernel and scattered-data methods emphasize flexibility under irregular sampling. Hybrid spherical approximation combines zonal basis functions with spherical harmonics inside a partition-of-unity framework (Rossi, 2014). Scaled zonal-kernel quasi-interpolation constructs approximate identities on the sphere and then discretizes the convolution by quadrature, balancing approximation and sampling errors through the scale parameter (Sun et al., 2024). Monte Carlo and quasi-Monte Carlo quasi-interpolation extend this line to noisy data and multilevel constructions without solving linear systems (Sun et al., 14 Oct 2025). Spherical k-nearest neighbors interpolation instead uses geodesic proximity and neighborhood-distribution debiasing for sparse heterogeneous geospatial data (Trempe, 2019).

Constrained polynomial approaches occupy an intermediate position. Interpolation-regression on the sphere computes a spherical polynomial of degree $\mathbb S^2$7 that interpolates data on a prescribed subset of nodes while fitting the remaining samples in a least-squares sense, with a KKT characterization and symmetry reductions in antipodally symmetric settings (Dell'Accio et al., 13 May 2026).

4. Derivative-aware local interpolation: Spherical Hermite Maps

A recent derivative-aware formulation is the Spherical Hermite Map, introduced as a derivative-augmented LUT representation for spherical functions in computer graphics. At each texel of a padded cubemap, four scalars are stored in RGBA:

$\mathbb S^2$8

where $\mathbb S^2$9 is texel spacing. Bicubic-Hermite reconstruction is then performed over a $r:S^2\to\mathbb R$0 footprint using tensor-product Hermite basis functions, so smooth value reconstruction and smooth chart-space derivatives arise from the same four texture fetches. The chart derivatives $r:S^2\to\mathbb R$1 are converted to the true tangent-plane gradient on the sphere by solving a $r:S^2\to\mathbb R$2 metric system built from $r:S^2\to\mathbb R$3 and $r:S^2\to\mathbb R$4. For a star-shaped surface $r:S^2\to\mathbb R$5 with $r:S^2\to\mathbb R$6, the outward normal is

$r:S^2\to\mathbb R$7

All ingredients needed for this normal come from the same Hermite evaluation, so no extra texture samples are required (Abouagour et al., 23 Feb 2026).

The runtime process is explicit: select cubemap face and chart coordinates from the view direction, compute the integer cell and local fractional coordinates, fetch the four RGBA texels, evaluate the Hermite basis and derivative polynomials, reconstruct $r:S^2\to\mathbb R$8 and its chart partials, convert them through the local metric tensor, and form the spherical gradient and normal. In the reported complexity comparison, bilinear value-only interpolation uses approximately 5 texture ops once finite-difference normals are included; 16-tap bicubic uses 16 texture ops; fast bicubic uses 4 bilinear fetches for values and typically 8 ops once gradients are added; Hermite uses 4 RGBA fetches, 16 scalars, and approximately 60 ALU total, while producing continuous $r:S^2\to\mathbb R$9 values and analytic gradients at no extra fetch cost. The storage overhead is 4 floats per texel versus 1 float for value-only LUTs (Abouagour et al., 23 Feb 2026).

Quantitatively, the method reports substantial gains. Compared to bilinear interpolation, Hermite yields $(u,v)\in[0,1]^2$0 dB PSNR at $(u,v)\in[0,1]^2$1, $(u,v)\in[0,1]^2$2 dB at $(u,v)\in[0,1]^2$3, and $(u,v)\in[0,1]^2$4 dB at $(u,v)\in[0,1]^2$5. Against 16-tap bicubic, it is $(u,v)\in[0,1]^2$6 dB at $(u,v)\in[0,1]^2$7 and $(u,v)\in[0,1]^2$8 dB at $(u,v)\in[0,1]^2$9, at one-quarter the fetch cost. On procedural planet terrain with a $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$0 LUT, mean angular normal error drops from $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$1 for bilinear+FD and $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$2 for bicubic+analytic $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$3 to $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$4 for Hermite, a 9% reduction versus bilinear. On an asteroid impostor with the same LUT size, mean error drops from $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$5 for bilinear+FD and $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$6 for bicubic+analytic to $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$7, a 13% reduction. The paper attributes stable specular highlights and the absence of “swimming” to the continuity and noise-free character of the analytic normals (Abouagour et al., 23 Feb 2026).

The method is demonstrated in spherical-harmonic glyph visualization, radial depth-map impostors for mesh level-of-detail, and procedural planet or asteroid rendering with spherical heightfields. Reported usage includes hundreds of thousands of glyphs at more than 30 FPS, billboard raycasting with only 4 fetches for LOD shading, and high-frequency procedural detail with minimal shader cost. The principal integration constraints are also explicit: storage is quadrupled, gutter texels must encode derivatives as well as values, mipmap generation requires recomputing derivatives per level, alternative chartings would require analogous derivative transport and seam handling, and at extremely low resolutions finite-difference errors in baked derivatives can cause minor seam artifacts (Abouagour et al., 23 Feb 2026).

5. Geodesic high-order interpolation and non-oscillatory reconstruction

For sampled curves on $HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$8, high-order interpolation is built from repeated SLERP rather than from local polynomial fits in Euclidean coordinates. SIDER-$HM(\Omega,l)=20\log_{10}|H(\Omega,l)|,$9 is a recursive construction in which SIDER2 is formed from three points using extrapolated controls, and higher orders are generated by a Neville-type recursion

$\Omega\in S^2$0

The original SENO paper describes SIDER-$\Omega\in S^2$1 as a family of $\Omega\in S^2$2 spherical Bézier-like interpolants and reports $\Omega\in S^2$3 convergence on sufficiently smooth generating paths (Fong et al., 2022).

Later local-consistency analyses make the order structure precise. In geodesic normal coordinates, SLERP is affine to leading order and has no quadratic correction; the first nonlinear term is cubic (Leung, 11 Jun 2026). In the quaternion–Lie algebra formulation, the BCH expansion shows that the first nonlinear correction is cubic and commutator-valued, and that the recursive SIDER step raises the order by one under smoothness and small-stencil assumptions (Leung, 11 Jun 2026). The resulting estimates are

$\Omega\in S^2$4

with explicit low-order cases: SIDER2 is $\Omega\in S^2$5, SIDER3 is $\Omega\in S^2$6, and SIDER4 is $\Omega\in S^2$7 (Leung, 11 Jun 2026).

A distinct issue is oscillation near kinks or discontinuities in higher derivatives. SENO addresses this by constructing multiple local SIDER-$\Omega\in S^2$8 candidates on overlapping stencils and selecting the one with least variation on the target interval. The variation metric is defined by sampling the candidate curve on subdivided subintervals and summing geodesic distances; the chosen stencil is therefore the least-varying local high-order reconstruction (Fong et al., 2022). In the reported experiments, SLERP is second order, SQUAD is third order on smooth data but degrades near the kink, while SENO2 remains $\Omega\in S^2$9 and SENO3 remains $l$0 even in the presence of the kink (Fong et al., 2022).

These results distinguish two notions of spherical interpolation that are sometimes conflated. Spherical field interpolation reconstructs sampled values over $l$1; geodesic curve interpolation reconstructs a trajectory whose samples lie on $l$2. The two settings share the same manifold but use different consistency arguments and different smoothness criteria.

6. Learned interpolation operators and broader feature-space extensions

In learned spherical feature interpolation, the interpolation operator is embedded into a trainable model. For point clouds, the spherical interpolated convolutional network replaces the traditional grid-shaped 3D convolution operator by a spherical interpolated convolution operator defined on densely packed spherical cells around FPS-selected centers. A self-learned distance-feature density modifies inverse-distance interpolation weights by feature correlation. On ScanNet, the method reports 63.6% mIoU, compared with 55.6% for PointConv, 61.0% for SPH3D-GCN, 61.8% for HPEIN, 63.3% for MCCNN, and 63.4% for joint point-based; on Paris-Lille-3D it reports 70.1% versus 68.3% for HDGCN, 66.9% for MS3-DVS, and 56.3% for RF-MSSF. The paper also reports that a spherical kernel uses 9.06 M parameters, 5.10 G FLOPs, and 264 ms per 8192 points, compared with 15.34 M parameters, 7.43 G FLOPs, and 397 ms for a cubic kernel (Wang et al., 2020).

For HRTFs, two related directions appear. A spherical CNN formulation expands HRTF magnitudes in spherical harmonics, performs zonal spectral convolutions, and reconstructs dense HRTF maps from sparse measurements. On 7 HUTUBS test subjects, the reported average distortion is approximately 2.17 dB LSD, versus approximately 6 dB for SH-only interpolation, approximately 5.11 dB for SH+DNN, approximately 4.19 dB for a neural-field method, and approximately 4.36 dB for CNN+GAN with planar projections (Chen et al., 2023). A later spherical neural-process meta-learner combines spherical set convolutions, spherical-harmonic CNN blocks, and a Gaussian decoder, and reports up to 3 dB relative error reduction compared to state-of-the-art interpolation methods. The same work states that the required correcting HRTF data point count can be reduced from 50 to 28 points to reach an average of $l$3 dB relative error per interpolated feature, and that uncertainty calibration achieves approximately 1.0 dB Mean Calibration Distance (Thuillier et al., 2023).

Beyond $l$4, the hyperspherical case uses SLERP directly on unit-length feature embeddings. SlerpFace rotates face templates toward a secret key by spherical linear interpolation, then applies group-wise partitioning and dropout to improve irreversibility. In the reported LFW inversion experiment, the unprotected setting has inversion similarity approximately 0.46 and SRRA approximately 98.8%, whereas SlerpFace reports inversion similarity approximately 0.05 and SRRA approximately 0.10% (Zhong et al., 2024). Although this is not interpolation of a spherical field, it shows that “spherical feature interpolation” can also denote operations on learned feature vectors constrained to a hypersphere.

Across these learned and non-learned formulations, a common technical pattern emerges. Geometry-aware interpolation is used to avoid artifacts that arise when spherical structure is ignored: artificial pole boundaries in tensor-product spherical coordinates, neglect of spatial acoustic features in planarized HRTF models, or noisy derivative estimation in finite-difference shading (Chiwere et al., 2024, Chen et al., 2023, Abouagour et al., 23 Feb 2026). A plausible implication is that the choice of interpolation geometry is not merely a numerical detail; it determines which invariances, continuity properties, and downstream differential quantities remain available after reconstruction.