# HEALPix Support: Regridding and Tracking

This document describes spectral regridding to HEALPix and object detection on a one-dimensional HEALPix grid.

## 1. Overview
HEALPix represents a global field as equal-area, iso-latitude pixels with a neighbor graph. The one-dimensional topology avoids the pole singularity and latitude-dependent cell areas of a regular latitude-longitude mesh. PyStormTracker can regrid a regular or reduced-Gaussian field to HEALPix and run object detection on that graph.

## 2. Spectral Regridding (`SpectralRegridder`)
`SpectralRegridder` uses `ducc0.sht` for spherical harmonic analysis of an input grid and synthesis on a target grid.

- **Supported Inputs**: Clenshaw-Curtis (CC) and Gauss-Legendre (GL).
- **Supported Outputs**: CC, GL, and HEALPix.
- **Spectral Logic**:
    - **Analysis**: Extracts spherical harmonic coefficients ($a_{lm}$) from 2D grids using `ducc0.sht.analysis_2d`.
    - **Synthesis**: Projects coefficients onto the target grid. For HEALPix, it uses `ducc0.sht.synthesis` with `geometry` parameters derived from `ducc0.healpix.Healpix_Base.sht_info()`.
    - **Spectral Truncation**: Supports explicit $L_{max}$ and $M_{max}$ band limits. If omitted, the truncation is inferred from the input longitude count.

## 3. HEALPix Tracking Algorithm (`HealpixTracker`)
`HealpixTracker` implements the `Tracker` protocol for HEALPix pixel arrays. Tracking currently supports the serial backend.

### 3.1. 1D Graph Topology
Unlike 2D meshes where neighbors are found via index offsets, HEALPix neighbors are determined via a precomputed adjacency list.
- During initialization, the `HealpixDetector` generates a **neighbor table** of shape `(8, N_pixels)` using `ducc0.healpix.Healpix_Base.neighbors()`.
- All detection kernels operate on this 1D graph, using the table for topological lookups.

### 3.2. Connected Component Labeling (CCL)
The tracker groups adjacent pixels into objects with the Numba kernel `_numba_healpix_ccl`.
- **Algorithm**: Iterative label propagation over the 1D graph until convergence.
- **Constraints**: Supports `threshold` filtering and `min_points` object-size constraints.

### 3.3. Spherical Subgrid Refinement
Optional `subgrid_refine_healpix` applies these steps:
1.  **Local Projection**: For each detected extremum at pixel $P$, it projects $P$ and its 8 neighbors onto a local **equirectangular plane** centered at $P$.
2.  **Numerical Stability**: Coordinate scaling/normalization is applied to the local projected coordinates to prevent matrix ill-conditioning when solving the least-squares system.
3.  **Surface Fitting**: An unconstrained least-squares fit is performed to find the coefficients of a local quadratic surface: $z = Ax^2 + By^2 + Cxy + Dx + Ey + F$.
4.  **Analytical Stationary Point**: The candidate offset $(dx, dy)$ is found by solving the system where partial derivatives $\frac{\partial z}{\partial x} = 0$ and $\frac{\partial z}{\partial y} = 0$.
5.  **Inverse Projection**: The refined coordinates are projected back to standard Latitude and Longitude.

## 4. Engineering Standards
- **Dependencies**: HEALPix operations use the existing `ducc0` dependency; `healpy` is not required.
- **Numba JIT**: Graph traversal and local quadratic fitting kernels are compiled with `cache=True` and `nogil=True`.
- **Defaults**: HEALPix tracking defaults to T5-42 filtering and enabled subgrid refinement. Both can be disabled explicitly. The `min_lifetime` constructor value is currently stored but not applied during HEALPix linking.
- **Tracker protocol and output**: `HealpixTracker` implements the common `Tracker` protocol and returns the array-backed `Tracks` model. Dask and MPI tracking are not implemented.

## 5. Usage Example
```python
import xarray as xr

from pystormtracker import HealpixTracker, SpectralRegridder

# Regrid one ERA5 frame (CC) to HEALPix (Nside=64)
ds = xr.open_dataset("data.nc")
field = ds["msl"]
frame = field.isel({field.dims[0]: 0})
regridder = SpectralRegridder(lmax=42)
da_hp = regridder.to_healpix(frame, nside=64)

# Track the first eight time steps after automatic HEALPix conversion
tracker = HealpixTracker()
tracks = tracker.track(
    infile=field.isel({field.dims[0]: slice(0, 8)}),
    varname="msl",
    mode="min",
    threshold=1000.0,
    filter=True,  # Apply T5-42 spectral filtering
)
```

The first part demonstrates one-frame regridding; the second passes a regular three-dimensional latitude-longitude field to `HealpixTracker.track()`, which triggers T5-42 filtering and HEALPix conversion by default. For an already regridded time-by-cell field, preprocessing must be performed before tracking.