healpix_analyse.kernel_pyramid#

Compact per-band kernels for pyramidal HEALPix convolution.

Theoretical status (read before using)#

Let W be the analysis operator of a HealPixDecomp (the stacked Down/Up chain) and S its synthesis operator, so that S W = I exactly (this is what makes invert() an exact inverse of compute()). Convolving a map x with a target operator K and then decomposing, or decomposing and then acting on the pyramid coefficients, are not interchangeable in general. Write the “acting on coefficients” operator as B_K such that the intended pipeline is

y = S · B_K · W · x

which reduces to plain convolution (y = K x) only for the specific choice B_K = W K S – and that operator is, in general, dense across bands: detail at one scale can leak into neighbouring scales through B_K.

HealPixKernelPyramid builds a block-diagonal approximation of B_K: one small, compact HealPixConv kernel acting independently on each pyramid band. This is a deliberate simplification, not a mathematical identity, and it is not the same claim as “the pyramid reconstructs exactly” (that property belongs to S W = I alone and holds regardless of which kernel, if any, is inserted). The approximation error of using a block-diagonal B_K in place of the true (block, inter-band-coupled) operator must be measured against a direct, non-pyramidal reference convolution – see healpix_analyse.validation – and is reported, band by band and kernel-size by kernel-size, in the accompanying validation notes rather than assumed to be small.

Note that the literal B_K = W K S choice above is not merely expensive to compute – for the masked/hole-filling use this pyramid exists for, it is useless, because y = S B_K W x = S W K S W x = K x (using S W = I twice): the multiscale structure cancels out completely and the pipeline degenerates to a single, ordinary convolution with K at the finest resolution alone, with only K’s own compact support to reach across a hole. A block-diagonal B_K is not a shortcut around computing the dense operator; it is the thing that actually keeps per-band information distinct so that filling can draw on coarser, wider-reaching bands.

Where does each band’s kernel come from? from_kernel takes a continuous profile kernel(rho_pix, phi_rad) and realizes it, once, on the exact discrete stencil geometry of the finest band (J=0, see weights_from_finest_band below); every coarser band’s HealPixConv reuses that identical tap vector, unchanged. The pyramid is built from a single kernel defined at the finest resolution, not from n_bands independent re-evaluations of the same formula on each band’s own, slightly different real HEALPix geometry (nearest-neighbour distances on the sphere are not perfectly self-similar across nside the way an idealized flat grid’s would be). Re-evaluating per band was the previous default and is still available (weights_from_finest_band=False) for comparison, but it means two bands nominally sharing “the same” kernel can carry slightly different discrete taps for reasons that have nothing to do with the kernel itself.

``from_kernel``/``calibrate`` vs. ``calibrate_joint``: two different problems. Both from_kernel (an analytic profile, shared or resampled per band) and calibrate (a least-squares fit, one band at a time) choose each band’s kernel independently – there is no single common target they are jointly trying to reproduce, and neither can represent a target many times wider than one band’s own compact support (measured directly in test_calibrate_on_a_much_wider_target_does_not_silently_claim_success). HealPixKernelPyramid.calibrate_joint() solves the actually different problem of decomposing one wide, slowly-decaying target kernel – defined once, at the finest resolution – into small per-band kernels that reconstruct it together, through the real S B W pipeline (not the degenerate literal B=W K S case above, which collapses to y=Kx identically). This lets a Lorentzian- or power-law-tailed kernel, far wider than any single compact_kernel_sz stencil, be represented well using only small per-band kernels – the pyramid’s own geometric coarsening supplies most of the reach, and each band’s fitted kernel only has to correct the local shape at its own scale (measured: ~3% relative RMS error for a Lorentzian(scale=6 pixels) target on a 4-band pyramid using 5x5 per-band kernels throughout, vs. ~87% for the same target fit to a single band alone – see test_calibrate_joint_decomposes_a_wide_kernel_across_bands).

This module does not claim to reproduce the specific published construction of Farbman, Fattal & Lischinski, “Convolution Pyramids” (ACM TOG 2011) – their boundary-matched, closed-form per-family derivations are not implemented here. calibrate_joint solves the same class of problem (representing a wide kernel as a small-kernel pyramid) by direct joint least-squares against a brute-force reference instead, which is more general (any kernel profile the reference oracle can evaluate) but purely numerical, with no closed-form guarantees – its accuracy must be measured per kernel family and per Jmax, the same as every other calibration path in this module.

Attributes#

KernelFn

A kernel profile fn(rho_pix, phi_rad) -> weight.

Classes#

HealPixKernelPyramid

One fixed, non-learnable HealPixConv per pyramid band.

Functions#

`kernel_gaussian`(→ KernelFn)	Isotropic Gaussian, `exp(-rho^2 / (2 sigma^2))`. Control/reference case.
`kernel_exponential`(→ KernelFn)	Isotropic exponential ("K_exp"), `exp(-rho / scale)`. Heavier tail than Gaussian.
`kernel_lorentzian`(→ KernelFn)	Isotropic Lorentzian ("K_lor"), `1 / (1 + (rho/scale)^2)`. Slowly-decaying tail.
`kernel_beta`(→ KernelFn)	Isotropic generalized-Lorentzian/"Moffat" profile ("K_beta"),
`kernel_anisotropic_gaussian`(→ KernelFn)	Anisotropic scalar-field Gaussian `K(rho, phi)`, elongated along

Module Contents#

healpix_analyse.kernel_pyramid.KernelFn#

A kernel profile fn(rho_pix, phi_rad) -> weight.

rho_pix is the radial distance from the stencil centre, in pixel units (i.e. sqrt(du**2 + dv**2) on the kernel_sz x kernel_sz index grid, matching healpix_analyse.convol._local_kernel_grid()’s own angular spacing alpha_pix); phi_rad is the azimuth in the fixed North-Pole stencil frame, in radians. Using pixel units (rather than a physical angle) keeps a kernel’s shape in pixels identical at every pyramid band; its physical footprint then grows automatically from band to band because the pixel size itself grows.

healpix_analyse.kernel_pyramid.kernel_gaussian(sigma_pix: float) → KernelFn[source]#: Isotropic Gaussian, exp(-rho^2 / (2 sigma^2)). Control/reference case.

healpix_analyse.kernel_pyramid.kernel_exponential(scale_pix: float) → KernelFn[source]#: Isotropic exponential (“K_exp”), exp(-rho / scale). Heavier tail than Gaussian.

healpix_analyse.kernel_pyramid.kernel_lorentzian(scale_pix: float) → KernelFn[source]#: Isotropic Lorentzian (“K_lor”), 1 / (1 + (rho/scale)^2). Slowly-decaying tail.

healpix_analyse.kernel_pyramid.kernel_beta(scale_pix: float, beta: float = 2.0) → KernelFn[source]#: Isotropic generalized-Lorentzian/”Moffat” profile (“K_beta”), (1 + (rho/scale)^2)^(-beta/2). beta=2 reduces to kernel_lorentzian(); larger beta decays faster (approaches Gaussian-like compactness).

healpix_analyse.kernel_pyramid.kernel_anisotropic_gaussian(sigma_major_pix: float, sigma_minor_pix: float, orientation_rad: float = 0.0) → KernelFn[source]#

Anisotropic scalar-field Gaussian K(rho, phi), elongated along orientation_rad (measured from the stencil’s local x-axis) with principal widths sigma_major_pix/sigma_minor_pix.

Evaluated in the fixed North-Pole stencil frame, this is combined with HealPixConv’s per-pixel gauge rotation, so the anisotropy is expressed consistently in every output pixel’s local tangent frame (the same convention HealPixConv itself uses) – not in a single global frame that would be meaningless on a sphere.

class healpix_analyse.kernel_pyramid.HealPixKernelPyramid[source]#

One fixed, non-learnable HealPixConv per pyramid band.

Built from from_kernel() (or calibrate(), for an optional data-driven refinement). Shares the exact stencil geometry, gauge convention, and cell-id domain of a given HealPixDecomp, band for band, so that apply() can be composed directly with compute()/compute_weighted() and invert().

This is a block-diagonal approximation – see the module docstring (block-diagonal across pyramid bands; see also channels below for the unrelated, additional sense in which multi-channel kernels here are block-diagonal across channels).

Multiple channels (channels > 1, e.g. 3 for RGB) are supported as an explicit second, independent block-diagonal structure: one in_channels=channels, out_channels=channels HealPixConv per band, whose [channels, channels, P] kernel is zero off the diagonal – i.e. every channel is filtered independently with the same profile, with no cross-channel mixing. This is deliberately not full HealPixConv generality (a dense [C_out, C_in, P] kernel, or per-channel gauge counts) – the point is only to filter several co-registered scalar fields (e.g. R, G, B) through the same fixed kernel pyramid in one call, cheaply and without the shape pitfall documented on apply().

decomp: healpix_analyse.decomp.HealPixDecomp#

convs: tuple[healpix_analyse.convol.HealPixConv, ...]#

kernel_sz: int#

gauge_type: str#

n_gauges: int#

channels: int = 1#

classmethod from_kernel(decomp: healpix_analyse.decomp.HealPixDecomp, kernel: KernelFn, *, compact_kernel_sz: int = 5, gauge_type: str = 'phi', n_gauges: int = 1, singularity_lonlat=None, ref_direction=None, bands: Sequence[int] | None = None, channels: int = 1, weights_from_finest_band: bool = True, ellipsoid: str = 'sphere', dtype: torch.dtype | None = None, device=None) → HealPixKernelPyramid[source]#

Build one fixed HealPixConv per band from an analytic kernel profile.

Parameters:

decomp (HealPixDecomp) – The pyramid whose bands this kernel pyramid will act on.
kernel (KernelFn) – fn(rho_pix, phi_rad) -> weight. By default (see weights_from_finest_band) evaluated exactly once, on band 0’s (the finest band’s) own stencil geometry (see _stencil_pixel_polar()), and that single discrete tap vector is then reused unchanged, tap for tap, to build every band’s HealPixConv – the pyramid is constructed from one kernel given at the finest resolution, not from independent per-band reconstructions of it. The kernel’s compact support (in pixel units) is therefore identical band to band by construction, not merely by the profile’s own scale-invariance; its physical footprint still grows with the coarsened pixel size.
compact_kernel_sz (int, default 5) – Odd stencil size, e.g. 5 for a 5x5 stencil (P = 25 taps).
bands (sequence of int, optional) – Restrict construction to these band indices (0 = finest detail, decomp.n_scales = coarse residual). Defaults to every band. The finest band’s geometry (decomp.levels[0]) is still used to realize the shared kernel weights even if band 0 itself is excluded here.
weights_from_finest_band (bool, default True) – If True (default, see above), kernel is sampled once, on band 0’s stencil, and every band shares that identical tap vector. If False, restores the previous behaviour: kernel is re-evaluated independently at each band’s own _stencil_pixel_polar() geometry. The two differ only because real HEALPix pixel geometry is not perfectly self-similar across resolutions – kept for comparison, not because it is the recommended default.
channels (int, default 1) – Number of independent scalar channels each band’s kernel will act on (e.g. 3 for RGB), built as an in_channels=out_channels= channels HealPixConv per band whose kernel is block-diagonal across channels – the same profile applied to every channel with no cross-channel mixing (see the class docstring). Requires n_gauges=1 (not checked against n_gauges passed here unless channels>1, in which case a value other than 1 raises, since the general [G, C, C, P] case with G>1 cannot be reshaped back to a plain [channels, N] band – a currently unaddressed limitation, not a silent approximation). Pass data shaped [..., channels, N] to forward() (see apply()).
ellipsoid (str, default "sphere") – Geometry passed to every band’s HealPixConv. This module defaults to the true sphere (see the module docstring, section A.5), which is not the same convention some real data stores use – e.g. EOPF HEALPix products declare "wgs84". Pass the same ellipsoid the data (and decomp, if it uses one other than its own default) actually use, or geometry is silently mislabeled. Accepted case-insensitively (canonicalized against healpix_geo’s own registry).

classmethod calibrate(decomp: healpix_analyse.decomp.HealPixDecomp, reference_fn: Callable[[numpy.ndarray, int], numpy.ndarray], *, compact_kernel_sz: int = 5, gauge_type: str = 'phi', n_probes: int = 128, n_excitations: int = 4, seed: int = 0, ridge: float = 1e-06, channels: int = 1, ellipsoid: str = 'sphere', dtype: torch.dtype | None = None, device=None) → HealPixKernelPyramid[source]#

Fit one kernel per band by least squares against a reference operator.

This is an explicitly approximate, single-band-at-a-time calibration: it does not account for inter-band coupling (the same block-diagonal caveat as from_kernel() applies), and its accuracy should be checked the same way, against the same reference.

For each band j, this solves for the P kernel taps w that minimise sum_i (sum_p w_p * x_interp[i, p] - target[i])^2 over a set of probe output pixels i and n_excitations independent random excitation fields, where x_interp[i, p] is the value HealPixConv itself bilinear-interpolates for tap p of pixel i (obtained by running the same HealPixConv with a one-hot kernel, so the fit is against the exact runtime stencil geometry) and target = reference_fn(...) applied to the same excitation field. Using a random excitation (rather than probing with per-pixel delta inputs placed at each probe) is required for correctness: for an off-centre tap, a delta exactly at the probe pixel does not generally coincide with that tap’s rotated, bilinear-bound sample location, so per-probe deltas would recover the wrong design matrix.

reference_fn(cell_ids, level) -> operator must return a callable usable as operator(x) -> y mapping a [len(cell_ids)] map to a [len(cell_ids)] map – i.e. it should already be restricted/consistent with this decomposition’s own geometry (masks, units, gauge) for a fair comparison.

channels, as in from_kernel(): the single scalar profile fitted per band is replicated block-diagonally across channels independent channels (no cross-channel mixing); fitting itself is always done against a scalar reference (one profile per channel, not one fit per channel), so reference_fn need not know about channels at all.

classmethod calibrate_joint(decomp: healpix_analyse.decomp.HealPixDecomp, target_kernel: KernelFn, *, target_kernel_sz: int = 41, compact_kernel_sz: int = 5, gauge_type: str = 'phi', n_probes: int = 128, n_excitations: int = 6, seed: int = 0, ridge: float = 1e-06, ellipsoid: str = 'sphere', dtype: torch.dtype | None = None, device=None) → HealPixKernelPyramid[source]#

Decompose a single, wide target kernel (defined once, at the finest resolution) into a genuine pyramid of small per-band kernels.

This is a different problem from from_kernel() and calibrate(): both of those pick, independently for each band, a kernel that acts well at that band’s own resolution – there is no single common target they are jointly trying to reproduce, and (per the module docstring, and demonstrated by test_calibrate_on_a_much_wider_target_does_not_silently_claim_success) neither can represent a target many times wider than one band’s own compact support.

calibrate_joint instead fits all bands’ kernels together against one fixed target target_kernel, evaluated only at the finest band’s own resolution (decomp.levels[0]) via a direct, wide-support brute-force reference (direct_spherical_convolution() with kernel_sz=target_kernel_sz – make this comfortably larger than the target’s actual decay, e.g. target_kernel_sz~=6*scale_pix for a Lorentzian/exponential tail). The fitted operator is the real composed pipeline y = S(B(W(x))) (compute() then each band’s own small HealPixConv then invert()) – not band-by-band in isolation – so each band’s fitted kernel only has to supply whatever the pyramid’s own geometric coarsening does not already give it for free. This is why a wide, slowly-decaying target (Lorentzian, power-law) can be represented well by small (compact_kernel_sz) per-band kernels: most of the target’s reach comes from the pyramid’s coarsening itself, and each band only has to correct the local shape at its own scale. See docs/pyramid_convolution.md for the derivation and honest limits (reach is still bounded by decomp.Jmax; a target wider than the coarsest band’s own reach cannot be represented regardless of how the taps are fit).

Parameters:

decomp (HealPixDecomp) – Must be built with cell_ids=None (full sphere) – fitting needs a well-defined “external” pixel order to compare against the direct reference on the same footprint, and the full sphere keeps that unambiguous. A partial-domain decomposition can still use the resulting HealPixKernelPyramid afterwards (the fitted per-band taps are footprint-independent, like from_kernel()’s), just not for the fit itself.
target_kernel (IsoKernelFn) – fn(rho_pix, phi) -> weight, evaluated in finest-band pixel units (same convention as healpix_analyse.validation). Isotropic only (no phi dependence), matching the independent oracle’s own documented scope.
target_kernel_sz (int, default 41) – Odd stencil size for the brute-force reference – how far the direct oracle searches for neighbours, in finest-band pixel units. This is independent of, and normally much larger than, compact_kernel_sz (that one bounds each band’s own fitted kernel, not the target being approximated).
compact_kernel_sz (int, default 5) – Odd stencil size for every band’s own fitted kernel.
n_probes, n_excitations, seed, ridge (as in :meth:`calibrate`, but) – probing now runs the full analysis/synthesis chain, so cost scales with n_excitations * n_bands * compact_kernel_sz**2.
Only ``channels=1`` is supported here (unlike :meth:`from_kernel`/
:meth:`calibrate`) – fitting a joint multi-channel operator is not
implemented, a currently unaddressed limitation, not a silent
approximation.

apply(bands: Sequence) → tuple[source]#

Apply each band’s fixed kernel to the matching pyramid band.

With channels=1 (the default), each band is passed straight to its HealPixConv – see that class’s own forward for its exact input/output shape convention, including a known asymmetry: a 2-D [B, N] input always returns a 3-D [B, 1, N] output, never squeezed back to 2-D. This method does not paper over that, so a band shaped [C, N] used as a poor-man’s channel batch (e.g. RGB packed as [3, N]) comes back [3, 1, N], not [3, N] – build the kernel pyramid with channels=3 instead (see the class docstring) to get real, correctly shaped multi-channel support.

With channels>1, each band must be shaped [..., channels, N]; this method reshapes it to the [B, channels, N] its HealPixConv (built with matching in_channels=out_channels= channels) actually expects, applies it, and reshapes the result back to [..., channels, N] – any leading batch dims are otherwise untouched.