healpix_analyse.wide_conv ========================= .. py:module:: healpix_analyse.wide_conv .. autoapi-nested-parse:: Convolution by a *wide* kernel, given as an image, via a scale pyramid. Motivation ---------- :class:`~healpix_analyse.convol.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 --------------------- :class:`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 :class:`~healpix_analyse.convol.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 :meth:`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 :class:`~healpix_analyse.kernel_pyramid.HealPixKernelPyramid` / :class:`~healpix_analyse.pyramid_conv.HealPixPyramidConv`; folding them in here would complicate the one thing this class is for. Classes ------- .. autoapisummary:: healpix_analyse.wide_conv.HealPixWideConv Module Contents --------------- .. py:class:: 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) Convolve by a kernel given as a ``(2n+1, 2n+1)`` image, via a pyramid. :Parameters: * **kernel_image** (:py:class:`array-like`, :py:class:`shape ``(2n+1`, :py:class:`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 :func:`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** (:py:class:`int`) -- The HEALPix level ``kernel_image`` is sampled at. ``cell_ids`` passed to ``__call__`` must be at this level. * **Jmax** (:py:class:`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** (:py:class:`int`, *default* ``5``) -- Odd size of each band's fitted stencil. * **gauge_type, ellipsoid** -- Passed through to each band's :class:`~healpix_analyse.convol.HealPixConv`. * **ridge** (:py:class:`float`, *default* ``1e-12``) -- Tikhonov regularization on the per-band least-squares fit. * **dtype** (:py:class:`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. .. admonition:: Notes Not mask-aware: NaN in ``x`` propagates. Use :class:`~healpix_analyse.pyramid_conv.HealPixPyramidConv` in ``mode="normalized"`` for hole-filling on real, gappy data. .. admonition:: 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) # doctest: +SKIP >>> y = conv(x, cell_ids) # doctest: +SKIP .. py:attribute:: kernel_image .. py:attribute:: n .. py:attribute:: level .. py:attribute:: Jmax :value: 6 .. py:attribute:: compact_kernel_sz :value: 5 .. py:attribute:: gauge_type :value: 'phi' .. py:attribute:: ellipsoid :value: 'sphere' .. py:attribute:: ridge .. py:attribute:: dtype :value: Ellipsis .. py:attribute:: device :value: None .. py:method:: lattice_offsets_m(level: int, lon: float, lat: float, n: int, *, r_earth: float = 6371008.8) :staticmethod: 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** (:py:class:`int`) * **lon, lat** (:py:class:`float`) -- Reference point, degrees. The lattice is centred on the cell containing it. * **n** (:py:class:`int`) -- Half-size; the returned arrays are ``(2n+1, 2n+1)``. * **r_earth** (:py:class:`float`, *default* ``6371008.8``) -- Sphere radius used to turn angles into metres. .. py:method:: from_function(fn, level: int, n: int, lon: float, lat: float, *, r_earth: float = 6371008.8, **kwargs) -> HealPixWideConv :classmethod: 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 :meth:`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** (:py:class:`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 :meth:`lattice_offsets_m`. * **\*\*kwargs** -- Passed to the constructor (``Jmax``, ``compact_kernel_sz``, ...). .. admonition:: Examples >>> conv = HealPixWideConv.from_function( # doctest: +SKIP ... 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) .. py:method:: from_radial(fn, level: int, n: int, lon: float, lat: float, *, r_earth: float = 6371008.8, **kwargs) -> HealPixWideConv :classmethod: 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 :meth:`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** (:py:class:`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 :meth:`lattice_offsets_m`. * **\*\*kwargs** -- Passed to the constructor (``Jmax``, ``compact_kernel_sz``, ...). .. admonition:: Examples >>> conv = HealPixWideConv.from_radial( # doctest: +SKIP ... lambda r: np.exp(-r / 500.0), # 500 m e-folding ... level=17, n=64, lon=2.3198, lat=48.8704, Jmax=6) .. py:method:: 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 :classmethod: 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 :meth:`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** (:py:class:`array-like`, :py:class:`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** (:py:class:`float` or :py:class:`(float`, :py:class:`float)`) -- Raster spacing in metres, ``sx`` or ``(sx, sy)``. * **level, n, lon, lat** -- As for :meth:`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** (:py:class:`float`, *default* ``0.0``) -- Value for lattice cells falling outside the raster. .. py:method:: kernel_as_field(cell_ids, centre_cell: Optional[int] = None) -> numpy.ndarray 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** (:py:class:`array-like` of :py:class:`int`, :py:class:`shape ``[N]```) * **centre_cell** (:py:class:`int`, *optional*) -- Where to centre the kernel. Defaults to the same cell the fit uses (the domain's own centroid, :meth:`reference_centre`). .. py:method:: reference_centre(cell_ids) -> int The cell the per-band fit is centred on (the domain's centroid). .. py:method:: band_kernels(cell_ids=None) -> tuple The fitted per-band stencils, each ``(compact_kernel_sz,) * 2``. .. py:method:: fit_residuals(cell_ids=None) -> tuple 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. .. py:property:: decomp :type: healpix_analyse.decomp.HealPixDecomp The pyramid built for the most recent domain.