healpix_resample.bicubic (radial cubic-convolution HEALPix regridding)#

healpix_resample.bicubic.BicubicResampler maps unstructured (lon_deg, lat_deg, val) samples onto a HEALPix grid using a radial generalization of Keys’ cubic convolution kernel — the same interpolation kernel behind cv2.INTER_CUBIC / PIL.Image.BICUBIC on regular pixel grids, adapted to scattered data.

It sits between NearestResampler / BilinearResampler (fixed, non-iterative, few neighbours) and PSFResampler (iterative CG deconvolution): a fixed, non-iterative local interpolation using more neighbours than bilinear, for users who want smoother results than bilinear without paying for a full CG solve.

Note on geodesy: the package manages the HEALPix authalic definition and the Earth ellipsoid with WGS84 through its geometry helper (healpix_geo).

What “bicubic” means here#

This package operates on unstructured samples found via a distance-based Gaussian-neighbourhood search (healpix_weighted_nearest in knn.py), not a structured pixel grid with clean (fx, fy) fractional offsets — so a textbook 2D bicubic convolution kernel doesn’t directly apply. What does generalize cleanly to scattered data is the radial profile of Keys’ cubic convolution kernel: replace Euclidean pixel distance with the geodesic distance self.d_m and length scale self.sigma_m already used by every other resampler in this package, and evaluate the standard piecewise-cubic weight (a = -0.5, the common default matching PIL/cv2):

u = d / sigma
w(u) = (a+2)|u|^3 - (a+3)|u|^2 + 1        for |u| <= 1
     = a|u|^3 - 5a|u|^2 + 8a|u| - 4a       for 1 < |u| < 2
     = 0                                   for |u| >= 2

Npt = 16 is the default neighbourhood size — the natural analogue of the classic 4x4 bicubic stencil on a structured grid, since Keys’ kernel has support |u| < 2 (roughly twice the reach of bilinear’s |u| < 1-ish support).

Unlike the nonnegative Gaussian/inverse-distance weights used elsewhere in this package, Keys’ kernel is signed — it goes negative for 1 < |u| < 2, which is what gives cubic convolution its sharpening property relative to bilinear. Two consequences:

  • The per-cell/per-sample weight sums used to normalize the sparse operators can land close to zero from cancellation between the positive and negative lobes. BicubicResampler.comp_matrix() floors these sums (relative to the unsigned accumulated weight) before dividing, rather than dropping affected cells.

  • invert() can genuinely overshoot/ring outside the local sample-value range — this is expected cubic-convolution behaviour, not a bug.

Class hierarchy and construction#

BicubicResampler is a KNeighborsResampler subclass, exactly like BilinearResampler: it overrides __init__ (to fix Npt=16, auto-correct ring_search_max, and accept area=), comp_matrix() (to build the sparse operators from the cubic-convolution weight instead of the Gaussian/IDW weight, plus a conservative-mode operator, see below), and resample() (to add conservative=True, see below). invert() is inherited unchanged from KNeighborsResampler.

from healpix_resample import BicubicResampler

op = BicubicResampler(
    lon_deg=lon,       # (N,) sample longitudes, degrees
    lat_deg=lat,       # (N,) sample latitudes, degrees
    level=level,       # HEALPix level, nside = 2**level
    device="cuda",
    dtype=torch.float32,
    nest=True,
)

Key parameters#

  • lon_deg, lat_deg: unstructured sample coordinates in degrees, shape (N,).

  • level: HEALPix level (nside = 2**level).

  • Npt: number of HEALPix neighbours per sample (default 16).

  • nest: HEALPix indexing scheme.

  • device, dtype: PyTorch placement and numerical type.

  • threshold: minimum accumulated (unsigned, Gaussian-weighted) support for a HEALPix cell to be kept — same meaning as for every other resampler in this package.

  • area: per-sample pixel area/weight, shape (N,). Only used by resample(conservative=True) (see below); ignored otherwise. Defaults to 1.0 for every sample.

Stored attributes#

After construction, BicubicResampler exposes the same attributes as every KNeighborsResampler subclass:

  • N: number of samples.

  • K: number of kept HEALPix cells (cell_ids).

  • cell_ids: (K,) HEALPix cell ids retained.

  • hi: (N, Npt) cell indices per sample (into cell_ids).

  • d_m: (N, Npt) geodesic distances (metres) to those cells.

  • M: sparse CSR (N, K) operator (per-cell-normalized), hval = y @ M.

  • MT: sparse CSR (K, N) operator, val_hat = hval @ MT.

  • M_cons: sparse CSR (N, K) operator (per-sample-normalized — a partition of unity per sample, up to the signed-kernel floor caveat below), used by resample(conservative=True).

  • area: (N,) per-sample area/weight (see above).

Methods#

resample(val, conservative=False)#

Sample space → HEALPix cell space.

  • conservative=False (default): hval = val @ self.M (per-cell-normalized weights) — smooth, but not exactly mass-conserving.

  • conservative=True: hval = (val * self.area) @ self.M_cons — each sample’s own (area-weighted) value is redistributed across its 16 nearest cells using weights normalized so each sample’s own weights sum to 1, guaranteeing sum_k hval[k] == sum_i (valid i) val[i] * area[i] exactly — see “Conservative mode” below.

Accepts val of shape (N,) or (B, N), np.ndarray or torch.Tensor; returns the same type/batch shape it was given, wrapped in a ResampleResults(cell_data, cell_ids).

invert(hval)#

HEALPix cell space → sample space: val_hat = hval @ self.MT (inherited, no CG; unaffected by conservative/area). Same (K,)/(B, K) and NumPy/Torch symmetry as resample.

Conservative mode (area=, resample(conservative=True))#

Added for issue #44 (“conservative bi-linear is missing”), and applied to BicubicResampler for consistency with BilinearResampler. Same construction as the bilinear case (see docs/user-guide/regrid_to_healpix_bilinear.md for the full derivation): M_cons reuses the exact weights already computed for MT (per-sample-normalized) under M’s (sample, cell) index layout, so each sample’s own contribution sums to 1 across the cells it links to and no value is invented or lost globally.

Caveat specific to this class’s signed kernel: unlike BilinearResampler’s non-negative inverse-distance weights, Keys’ cubic kernel is signed, so M_cons’s per-sample rows only sum to exactly 1 for samples whose norm_row wasn’t floored by _floor_signed in comp_matrix() — true for the vast majority of well-conditioned samples. For the rare, pathologically-cancelled sample that does hit the floor, that sample’s own contribution to the conservation identity is only approximate, not bit-exact — the same accepted trade-off _floor_signed already makes for ordinary (non-conservative) interpolation.

NaN handling under conservative=True is the same as BilinearResampler’s: a NaN sample’s value and area are both excluded, so the conservation identity holds over exactly the valid samples; a batch row where every sample is NaN comes back entirely nan.

Notes and practical tips#

  • Bicubic gives no benefit over bilinear on a perfectly linear field (both interpolate a plane exactly); its advantage shows up on fields with curvature.

  • Because the kernel is signed, invert() can produce values slightly outside the local sample range (overshoot/ringing) — expected, not a bug.

  • If your data is very sparse relative to the HEALPix level (few samples within 2*sigma of a cell), prefer BilinearResampler or NearestResampler, which need fewer nearby samples to produce a well-conditioned result.