healpix_analyse.wide_conv#

Convolution by a wide kernel, given as an image, via a scale pyramid.

Motivation#

HealPixConv convolves with a compact stencil: a kernel_sz x kernel_sz neighbourhood, so its reach is a couple of pixels. A kernel that is genuinely wide – an exponential with a 10-pixel scale, a Lorentzian, anything with a long tail – cannot be written as one compact stencil at all; the stencil would have to be as wide as the kernel, which is exactly what nobody wants to pay for.

The way out is the “convolution pyramids” idea (Farbman, Fattal & Lischinski, ACM TOG 2011): decompose the signal into a Laplacian pyramid, convolve each band with its own small kernel, and synthesize. A wide kernel is wide only in the finest band’s units; at band j the pixels are 2**j times larger, so the same physical width is 2**j times fewer pixels. Past a few stages, a very wide kernel is a handful of pixels across.

What this module adds#

HealPixWideConv is the single-call front end for that:

kernel_image = … # (2n+1, 2n+1) at level conv = HealPixWideConv(kernel_image, level, Jmax=6) y = conv(x, cell_ids) # x: [N] or […, N] -> same shape

You hand it the kernel as an image at one level – the thing you can plot and reason about – and it works out, by itself, the small per-band kernels that reproduce it through the pyramid. There is no per-band kernel to write by hand and no pyramid bookkeeping in the caller.

How the per-band kernels are obtained#

Not by sampling band_j(K) around the centre: the input to band j is band_j(delta), a small Laplacian bump rather than a Dirac, so convolving it by band_j(K)’s own values does not give band_j(K) back (measured: ~99% wrong). They are fitted: for each band, the compact_kernel_sz**2 taps w_j that best satisfy

band_j(delta) (*) w_j ~= band_j(K)

in the least-squares sense, probed through the real HealPixConv so the stencil geometry used to fit is exactly the one used to apply. Since synthesis is the exact inverse of analysis (S W = I), if every band’s fit were exact the result would be K exactly – so the per-band residuals, available from HealPixWideConv.fit_residuals(), are the whole error budget of the method.

The fit needs a domain, so it happens on the first call, on that call’s own cell_ids (and is cached per domain afterwards). That also keeps it honest about HEALPix geometry: pixel shape depends on latitude, and fitting on the data’s own patch beats fitting on a reference patch somewhere else.

Scope#

Deliberately single-channel and mask-free: one kernel, applied independently to every leading dimension of x. Channel mixing ([C_in, C_out, P] kernels) and NaN/weight-aware filtering live in HealPixKernelPyramid / HealPixPyramidConv; folding them in here would complicate the one thing this class is for.

Classes#

HealPixWideConv

Convolve by a kernel given as a (2n+1, 2n+1) image, via a pyramid.

Module Contents#

class healpix_analyse.wide_conv.HealPixWideConv(kernel_image, level: int, *, Jmax: int = 6, compact_kernel_sz: int = 5, gauge_type: str = 'phi', ellipsoid: str = 'sphere', ridge: float = 1e-12, dtype: torch.dtype = torch.float64, device=None)[source]#

Convolve by a kernel given as a (2n+1, 2n+1) image, via a pyramid.

Parameters:

  • kernel_image (array-like, shape ``(2n+1, 2n+1)``) – The kernel, as an image on the HEALPix base-face integer lattice at level: kernel_image[a, b] is the weight at offset (di, dj) = (b - n, a - n) from the centre, where i/j are the base-cell-local coordinates of healpix_geo.nested.healpix_to_base_cell_coordinates() (so rows are j, columns are i – the same convention as unfolding a NESTED tile into a square image). Must be square with an odd side.

  • level (int) – The HEALPix level kernel_image is sampled at. cell_ids passed to __call__ must be at this level.

  • Jmax (int, default 6) – Number of pyramid stages. The kernel’s reach grows roughly like compact_kernel_sz//2 * 2**Jmax finest-band pixels, so this is the knob that has to be large enough for the kernel you are asking for.

  • compact_kernel_sz (int, default 5) – Odd size of each band’s fitted stencil.

  • gauge_type, ellipsoid – Passed through to each band’s HealPixConv.

  • ridge (float, default 1e-12) – Tikhonov regularization on the per-band least-squares fit.

  • dtype (torch.dtype, default torch.float64) – Working precision. float64 is recommended: the fit is a small normal-equations solve and the coarse bands’ taps are large.

  • device (optional) – Torch device for the pyramid and the convolutions.

Notes

Not mask-aware: NaN in x propagates. Use HealPixPyramidConv in mode="normalized" for hole-filling on real, gappy data.

Examples

>>> import numpy as np, healpix_geo.nested as hgn
>>> from healpix_analyse.wide_conv import HealPixWideConv
>>> level, n = 10, 8
>>> di, dj = np.meshgrid(np.arange(-n, n + 1), np.arange(-n, n + 1))
>>> kernel_image = np.exp(-np.hypot(di, dj) / 3.0)          # 17x17
>>> conv = HealPixWideConv(kernel_image, level, Jmax=3)
>>> y = conv(x, cell_ids)
kernel_image#
n#
level#
Jmax = 6#
compact_kernel_sz = 5#
gauge_type = 'phi'#
ellipsoid = 'sphere'#
ridge#
dtype = Ellipsis#
device = None#
static lattice_offsets_m(level: int, lon: float, lat: float, n: int, *, r_earth: float = 6371008.8)[source]#

Metric offsets of the (2n+1, 2n+1) lattice around (lon, lat).

Returns (x_m, y_m), each (2n+1, 2n+1): the position of every lattice cell in a local east/north frame centred on the cell that contains (lon, lat), in metres, using an azimuthal-equidistant projection (so hypot(x_m, y_m) is the exact great-circle distance).

This is where the HEALPix deformation becomes visible and usable: the (i, j) lattice is not a square metric grid – its axes are neither orthogonal nor equally scaled, and how much they are sheared depends on latitude. Anything defined in metres (a function of x, y, a raster with a pixel size in metres) must be evaluated at these positions rather than indexed by (i, j).

Parameters:

  • level (int)

  • lon, lat (float) – Reference point, degrees. The lattice is centred on the cell containing it.

  • n (int) – Half-size; the returned arrays are (2n+1, 2n+1).

  • r_earth (float, default 6371008.8) – Sphere radius used to turn angles into metres.

classmethod from_function(fn, level: int, n: int, lon: float, lat: float, *, r_earth: float = 6371008.8, **kwargs) → HealPixWideConv[source]#

Build from a kernel defined as fn(x_m, y_m), metres.

fn is evaluated at each lattice cell’s true metric offset from the centre (see lattice_offsets_m()), so a kernel written in metres keeps its intended physical shape even though the (i, j) lattice it lands on is sheared.

Parameters:

  • fn (callable) – fn(x_m, y_m) -> weight, applied element-wise to two (2n+1, 2n+1) arrays (east and north offsets, metres).

  • level, n, lon, lat – As for lattice_offsets_m().

  • **kwargs – Passed to the constructor (Jmax, compact_kernel_sz, …).

Examples

>>> conv = HealPixWideConv.from_function(
...     lambda x, y: np.exp(-np.hypot(x, y) / 500.0),      # 500 m e-folding
...     level=17, n=64, lon=2.3198, lat=48.8704, Jmax=6)
classmethod from_radial(fn, level: int, n: int, lon: float, lat: float, *, r_earth: float = 6371008.8, **kwargs) → HealPixWideConv[source]#

Build from an isotropic kernel given as fn(r_m), metres.

The common case: the kernel depends only on the distance to its own centre. r is the exact great-circle distance in metres from the centre cell to each lattice cell (see lattice_offsets_m()), so the kernel is isotropic on the ground, not in pixel indices – which in the polar caps is a materially different object, since one step along i and one along j are not the same number of metres there.

Parameters:

  • fn (callable) – fn(r_m) -> weight, applied element-wise to a (2n+1, 2n+1) array of distances in metres.

  • level, n, lon, lat – As for lattice_offsets_m().

  • **kwargs – Passed to the constructor (Jmax, compact_kernel_sz, …).

Examples

>>> conv = HealPixWideConv.from_radial(
...     lambda r: np.exp(-r / 500.0),                      # 500 m e-folding
...     level=17, n=64, lon=2.3198, lat=48.8704, Jmax=6)
classmethod from_grid(grid, pixel_size_m, level: int, n: int, lon: float, lat: float, *, fill: float = 0.0, r_earth: float = 6371008.8, **kwargs) → HealPixWideConv[source]#

Build from a kernel given as a square raster in metres.

grid is a regular raster whose pixels are pixel_size_m metres apart, centred on its own middle pixel. It is bilinearly interpolated onto the lattice cells’ true metric positions (see lattice_offsets_m()), which is what accounts for the HEALPix deformation: dropping such a raster straight onto (i, j) indices would silently shear it, because the lattice is not a square metric grid.

Parameters:

  • grid (array-like, 2-D) – The kernel raster, grid[b, a] at (x, y) = ((a - cx) * sx, (b - cy) * sy) metres from its centre. Any shape; odd sides put the centre exactly on a pixel.

  • pixel_size_m (float or (float, float)) – Raster spacing in metres, sx or (sx, sy).

  • level, n, lon, lat – As for lattice_offsets_m(). n sets the size of the HEALPix-side kernel image the raster is resampled onto; make it large enough to cover the raster’s own extent.

  • fill (float, default 0.0) – Value for lattice cells falling outside the raster.

kernel_as_field(cell_ids, centre_cell: int | None = None) → numpy.ndarray[source]#

The kernel image laid onto cell_ids, centred on one cell.

This is the reference a convolution must reproduce: convolving a Dirac placed on centre_cell with this operator should give back this field. Useful for validation – and for getting the centring right, since a one-pixel offset between the Dirac and the reference is easily mistaken for a much larger error than it is.

Parameters:

  • cell_ids (array-like of int, shape ``[N]``)

  • centre_cell (int, optional) – Where to centre the kernel. Defaults to the same cell the fit uses (the domain’s own centroid, reference_centre()).

reference_centre(cell_ids) → int[source]#

The cell the per-band fit is centred on (the domain’s centroid).

band_kernels(cell_ids=None) → tuple[source]#

The fitted per-band stencils, each (compact_kernel_sz,) * 2.

fit_residuals(cell_ids=None) → tuple[source]#

Per-band relative residual of the fit – the method’s error budget.

||band_j(delta) (*) w_j - band_j(K)|| / ||band_j(K)||, fine band first. If these were all zero the convolution would reproduce the kernel image exactly.

property decomp: healpix_analyse.decomp.HealPixDecomp#

The pyramid built for the most recent domain.