The four resamplers

This notebook runs all four resamplers on the same dataset so you can compare their behaviour, accuracy, and speed.

Setup

import numpy as np
from healpix_resample import (
    NearestResampler,
    BilinearResampler,
    PSFResampler,
    CellPointResampler,
)

# Shared dataset: a small structured grid near the origin
ndata = 128
lon_grid, lat_grid = np.meshgrid(
    0.3 * np.arange(ndata) / ndata,
    0.3 * np.arange(ndata) / ndata,
)
lon = lon_grid.ravel()
lat = lat_grid.ravel()
val = lon  # simple field: value = longitude

level = 15  # nside = 32768 — high resolution

NearestResampler

Each point is assigned to its single nearest HEALPix cell. Fast and simple, but can produce blocky results.

nr_nearest = NearestResampler(lon_deg=lon, lat_deg=lat, level=level)
res_nearest = nr_nearest.resample(val)

rval_nearest = nr_nearest.invert(res_nearest.cell_data)
mse_nearest = np.mean((rval_nearest - val) ** 2)
print(f"Nearest  — output cells: {res_nearest.cell_data.shape[0]}, MSE: {mse_nearest:.2e}")

Nearest  — output cells: 28268, MSE: 3.09e-37

BilinearResampler

Uses the 4 nearest cells with distance-based weights. Smoother than nearest, good for locally grid-like data.

nr_bili = BilinearResampler(lon_deg=lon, lat_deg=lat, level=level)
res_bili = nr_bili.resample(val, lam=0.0)

rval_bili = nr_bili.invert(res_bili.cell_data)
mse_bili = np.mean((rval_bili - val) ** 2)
print(f"Bilinear — output cells: {res_bili.cell_data.shape[0]}, MSE: {mse_bili:.2e}")

Bilinear — output cells: 28268, MSE: 2.91e-08

/home/runner/work/healpix-resample/healpix-resample/healpix_resample/bilinear.py:74: UserWarning: Sparse CSR tensor support is in beta state. If you miss a functionality in the sparse tensor support, please submit a feature request to https://github.com/pytorch/pytorch/issues. (Triggered internally at /home/conda/feedstock_root/build_artifacts/libtorch_1772252348570/work/aten/src/ATen/SparseCsrTensorImpl.cpp:49.)
  self.M  = M_coo.to_sparse_csr()

PSFResampler

Applies a Gaussian PSF kernel around each sample and solves a damped least-squares problem with Conjugate Gradient. Best reconstruction quality — especially when data is dense or the field has fine structure.

nr_psf = PSFResampler(lon_deg=lon, lat_deg=lat, level=level, threshold=0.5, verbose=False)
res_psf = nr_psf.resample(val, lam=0.0)

rval_psf = nr_psf.invert(res_psf.cell_data)
mse_psf = np.mean((rval_psf - val) ** 2)
print(f"PSF      — output cells: {res_psf.cell_data.shape[0]}, MSE: {mse_psf:.2e}")
print(f"           CG iterations: {res_psf.cg_niters}")

PSF      — output cells: 27892, MSE: 1.44e-16
           CG iterations: 100

CellPointResampler

Special mode: encodes each point as a HEALPix cell ID at level 29. No interpolation — used for exact point indexing.

nr_zuniq = CellPointResampler(lon_deg=lon, lat_deg=lat)
res_zuniq = nr_zuniq.resample(val)

rval_zuniq = nr_zuniq.invert(res_zuniq.cell_data)
max_err = np.max(np.abs(rval_zuniq - val))
print(f"Zuniq    — output cells: {res_zuniq.cell_data.shape[0]}, max error: {max_err:.2e}")

Zuniq    — output cells: 16384, max error: 0.00e+00