Extend the native Gaussian blur kernel to per-channel (H,W,C) blurs

CHANGELOG.md

@@ -8,6 +8,16 @@ match the `VERSION` file and `v*` git tags.
 
 ### Changed
 
+- **The native Gaussian blur kernel now covers per-channel `(H, W, C)` blurs, not just single 2-D planes.**
+  New `gaussian_blur_spatial` (`src/background.py`) blurs each channel independently by looping the
+  existing native kernel per channel -- the same per-channel-native-call pattern this project's
+  `_median_filter_per_channel` already uses. Wired into `denoising.py::reduce_stars` (the default Phase 4
+  star-reduction step, `--no-star-reduce` to disable) and `exposure_fusion.py`'s Laplacian-pyramid blur
+  (`--hdr-blend-mode fusion`), neither of which the earlier 2-D-only kernel could reach. Measured ~2.8x at
+  a real full-resolution shape. `originvision_infer.py`'s resize prefilter uses the same tuple-sigma
+  pattern but was left on scipy -- it has its own tight numeric tolerances already validated against a
+  reference implementation, and swapping its blur would add unvalidated drift for no clear benefit.
+
 - **DBE's regression-surface upsample switched cubic -> bilinear interpolation.** `_fit_background_surface`
   (the standard, non-dense-field DBE path) blurs its coarse-grid upsample again immediately afterward
   (`sigma=patch_size*0.5`, typically >=32px), which absorbs whatever cubic's extra curvature term would

CLAUDE.md

@@ -234,7 +234,7 @@ with `bayerPattern` to override.
 - **Bayer hot-pixel fix, sigma-clipped median, and sparse hot-pixel replacement** (`src/debayer.py`, all three on the default Phase 1 hot path, found via profiling a real session): `_fix_hot_bayer` (`mode='bayer'`) was calling `scipy.ndimage.median_filter` directly on each 2x2 Bayer sub-channel instead of the already-existing native `median_filter_native` dispatcher (`_median_filter3`) — because those sub-channels are non-contiguous strided views (`result[dy::2, dx::2]`) and the dispatcher's fast-path guard requires C-contiguous input. Fixed by routing through a contiguous copy first (`np.ascontiguousarray(sub)`) — no new Rust kernel needed, just a dispatch fix; ~7.5x faster on its own (242ms → 32ms on a 1936x1096 frame, measured in-process with warm-up to avoid cold-start noise). `sigma_clipped_median_native` ports `_sigma_clipped_median` (iterative sigma-clip + median, quickselect-based like `median_inplace` but f64 throughout since it returns a single scalar with no coefficient-chaining concern) — shared by `green_equalize` (G1/G2 balance) and `_equalize_bayer_grid` (per-Bayer-position sky level), both calling it several times per frame on quarter-resolution strided views. `hot_pixel_box_replace_native` replaces `_fix_hot_rgb_impl`/`_fix_hot_mono_impl`'s `scipy.ndimage.uniform_filter(ch, size=3)` replacement step, which computed the box mean over the *entire* channel just to keep a handful of masked (hot) positions — the native kernel computes the 3x3 box mean only at flagged positions, a sparse win on top of the native-vs-scipy win. Together these cut a real sequential (single-threaded) profile of Phase 1 from 173.5s to 93.3s on a 214-frame session (`_fix_hot_bayer` 52.0s→6.7s, `_fix_hot_rgb_impl` 33.6s→13.8s, `_sigma_clipped_median`/`_equalize_bayer_grid` region 31s→16.4s) — the wall-clock improvement under full parallelism is smaller (Quality+Load is already ~4x parallelized across cores, and these three functions are only part of its total cost), but the underlying per-frame work is genuinely ~46% cheaper.
 - **Gram-matrix thin-SVD trick** (`src/robust_pca.py` `_thin_svd_wide`, used by `--master-method robust_pca`'s IALM solver): `gram_matrix_wide` (`D @ D.T` for a wide `(N, P)` matrix, rayon-parallel over the `N*(N+1)/2` upper-triangle pairs, mirrored to the lower triangle) and `small_times_wide` (`small @ data` for a small `(N,N)` left operand and wide `(N,P)` right operand, rayon-parallel over output rows, axpy-style accumulation to keep each inner loop reading a data row sequentially rather than striding by P) — together the two GEMMs of eigendecomposing the small `N x N` Gram matrix instead of calling `np.linalg.svd` directly on the full wide matrix. Measured on the realistic robust-PCA shape (N=20, P=18M): direct `np.linalg.svd` 37.8s, the Gram trick in plain numpy 15.7s (~2.4x — numpy's own SVD isn't specialized for this shape), with these native kernels 7.9s (~4.8x on top of the algorithmic win, ~9x combined) — numpy's own `@` got zero benefit from this machine's cores at this shape (8.1s default-threaded vs 8.9s forced single-threaded), which is the headroom the native kernels close. Numpy mirror is `_thin_svd_wide`'s own fallback path (same function, not a separate mirror), so no separate parity test file — `tests/test_robust_pca.py` validates the whole decomposition end-to-end instead.
 - **IALM per-iteration fusion** (`src/robust_pca.py` `robust_pca_pre_svd_input` + `robust_pca_iterate`, used by `robust_pca_decompose`): profiling a real `--flat-from-lights` run (auto-triggered by `--auto` when no flat frames exist, N=10 sampled lights, P=6.25M mono pixels) found `robust_pca_decompose` spending 100s of a 126s total in plain-numpy elementwise arithmetic around the SVD (`D - S + Y/mu`, the S soft-threshold, `D - L - S`, `Y += mu*residual`, the residual norm) — roughly a dozen full-`(N,P)`-array passes/iteration, each its own temporary, not counted in `_thin_svd_wide`'s own (already-native) time at all. `robust_pca_pre_svd_input` fuses the first into one native pass; `robust_pca_iterate` fuses the rest into one native pass, writing `S`/`Y` in place (allocated once outside the loop, not re-allocated every iteration) and returning the residual's Frobenius norm directly. Per-element arithmetic (including matching `np.sign`'s 0-at-exactly-0 semantics, not `f64::signum`'s +-1) is bit-identical to the numpy reference; the norm is a parallel-chunk reduction, so only close to numpy's sequential sum (immaterial — it's just a convergence-tolerance check). Cut this same real case from 126.9s to 62.5s (~2x). Profiling *that* found the next bottleneck was no longer numpy at all: the L-update's `(U * sigma_shrunk) @ Vt` reconstruction — the exact same small-`(N,N)`-times-wide-`(N,P)` shape `_thin_svd_wide`'s own Gram-trick GEMM is — was calling plain `np.matmul`, and this machine's numpy has no optimized BLAS (`numpy.show_config()` reports `blas: name: auto`; measured ~1.4 GFLOPS, reference-BLAS speed), so that one line was 32.9s of the post-fusion 62.5s on its own. Routing it through the already-existing `small_times_wide` native kernel instead cut the real case to 45.4s (~2.8x combined vs. the 126.9s baseline; full pipeline wall-clock on the same session 3m7.8s -> 1m46.0s). Numpy fallback is the original expressions (`np.sign`/`np.maximum`/`np.matmul`), unchanged and still exercised by `_HAS_NATIVE` toggled off in tests. End-to-end native-vs-numpy-fallback parity (not just the individual kernels) in `tests/test_native.py`.
-- **Gaussian blur** (`gaussian_filter_native`, `src/background.py`'s `_gaussian_blur` wrapper + `gaussian_filter_ds`): `correlate1d` (the C function `scipy.ndimage.gaussian_filter1d` calls per axis) was the single largest self-time item in every full-pipeline profile taken this session -- ~30 call sites project-wide (DBE, chroma denoising, local contrast, structure-tensor coherence, edge-band correction, registration, star masking, and more). From-scratch separable reimplementation (two rayon-parallel correlate1d-style passes, row-wise then column-wise) of scipy's default `mode='reflect'` behaviour: same kernel formula (`radius = floor(truncate*sigma + 0.5)`, normalised Gaussian PDF samples, order-0), same edge-duplicating boundary (`reflect_idx`, shared with the hot-pixel kernels). Not a port -- parity is numerical, verified to double-precision rounding (rtol/atol 1e-9) across sigma 0.8-96 and two shapes, not bit-exact, since the two implementations sum in a different order. ~4x faster at small sigma (0.8-5), ~1.3-1.5x at the large sigmas `gaussian_filter_ds`'s downsampled branch uses (24-96) -- larger sigma means a wider kernel radius, so the per-tap loop cost grows and narrows the win. `_gaussian_blur` (2-D float only, transparent scipy fallback otherwise) preserves the caller's dtype (a first version silently upcast float32 to float64 -- the native kernel is always float64 internally like every other kernel in this file, and the wrapper wasn't casting the result back down; caught by a dtype-parity test before it could double the memory footprint of every float32 caller). Wired into every direct `gaussian_filter` call site project-wide that uses a plain scalar sigma on a 2-D array (`background.py`, `denoising.py`, `cli.py`'s master smoothing, `exposure_fusion.py`, `merge.py`, `moving_objects.py`, `noise_validation.py`, `postprocess.py`, `quality.py`, `registration.py`, `star_removal.py`, `trail_reject.py`) -- left alone: the handful of call sites passing a per-axis sigma *tuple* (`denoising.py::reduce_stars`'s `sigma=(s,s,0)`, `exposure_fusion.py`'s pyramid up/downsample, `originvision_infer.py`'s resize prefilter, the last with its own tight numeric tolerances against a reference implementation already documented in that module) -- the kernel only supports a scalar sigma on a 2-D plane, and a tuple sigma there already falls back to scipy safely on its own (raises inside the wrapper's `float(sigma)` cast, caught).
+- **Gaussian blur** (`gaussian_filter_native`, `src/background.py`'s `_gaussian_blur` wrapper + `gaussian_filter_ds`): `correlate1d` (the C function `scipy.ndimage.gaussian_filter1d` calls per axis) was the single largest self-time item in every full-pipeline profile taken this session -- ~30 call sites project-wide (DBE, chroma denoising, local contrast, structure-tensor coherence, edge-band correction, registration, star masking, and more). From-scratch separable reimplementation (two rayon-parallel correlate1d-style passes, row-wise then column-wise) of scipy's default `mode='reflect'` behaviour: same kernel formula (`radius = floor(truncate*sigma + 0.5)`, normalised Gaussian PDF samples, order-0), same edge-duplicating boundary (`reflect_idx`, shared with the hot-pixel kernels). Not a port -- parity is numerical, verified to double-precision rounding (rtol/atol 1e-9) across sigma 0.8-96 and two shapes, not bit-exact, since the two implementations sum in a different order. ~4x faster at small sigma (0.8-5), ~1.3-1.5x at the large sigmas `gaussian_filter_ds`'s downsampled branch uses (24-96) -- larger sigma means a wider kernel radius, so the per-tap loop cost grows and narrows the win. `_gaussian_blur` (2-D float only, transparent scipy fallback otherwise) preserves the caller's dtype (a first version silently upcast float32 to float64 -- the native kernel is always float64 internally like every other kernel in this file, and the wrapper wasn't casting the result back down; caught by a dtype-parity test before it could double the memory footprint of every float32 caller). Wired into every direct `gaussian_filter` call site project-wide that uses a plain scalar sigma on a 2-D array (`background.py`, `denoising.py`, `cli.py`'s master smoothing, `exposure_fusion.py`, `merge.py`, `moving_objects.py`, `noise_validation.py`, `postprocess.py`, `quality.py`, `registration.py`, `star_removal.py`, `trail_reject.py`). **`gaussian_blur_spatial`** (`src/background.py`) extends coverage to this codebase's `sigma=(s, s, 0)` call sites -- blur each channel of an `(H, W, C)` array independently (never across channels) by looping `_gaussian_blur` per channel, same per-channel-native-call pattern as `postprocess.py`'s `_median_filter_per_channel`. Wired into `denoising.py::reduce_stars` (the default Phase 4 star-reduction step, `--no-star-reduce` to disable) and `exposure_fusion.py`'s Laplacian-pyramid up/downsample (`--hdr-blend-mode fusion`). Measured ~2.8x at a real full-res `(2033,3041,3)` shape, sigma=1.5 (495.6ms -> 177.7ms) -- lower than the 2-D kernel's own ~4x at that sigma since each channel still pays native-call overhead separately, not fused into one 3-D pass. Left alone: `originvision_infer.py`'s resize prefilter, which also uses a `sigma=(s,s,0)` tuple but has its own tight numeric tolerances already documented and validated against a reference implementation (cv2) -- swapping its blur to a different reimplementation would add a new, unvalidated source of drift for no clear need.
 - **PSF-matched drizzle kernel + IBP super-resolution** (`src/stacking.py`, `--drizzle-kernel psf` and `--super-res-iters`): `warp_affine_kernel_table` is `warp_affine_lanczos3`'s sibling for an arbitrary (non-separable) precomputed tap-weight table instead of the fixed Lanczos-3 formula — lets drizzle resample with the session's own estimated PSF (`build_drizzle_psf_table`, a Wiener-regularized inverse filter of the PSF, not the raw PSF shape — using the raw shape measurably broadens stars, convolving an already-blurred profile with itself again) as a matched filter. No fast-path branch like Lanczos-3's separable case: a general Moffat/Gaussian PSF isn't X/Y-separable, so every output pixel gathers the full `(2*halo+1)^2` neighbourhood; numpy mirror is `_warp_affine_kernel_table_numpy`. `iterative_back_projection` (Irani & Peleg 1991, `--super-res-iters`) refines a drizzle output afterward in pure Python/numpy (no native kernel — the per-frame forward-simulate/back-project loop is FFT- and `affine_transform`-bound, not a hot inner-loop shape this file's kernels target): forward-simulates what each original frame should look like given the current estimate (inverse-warp + PSF blur via `scipy.signal.fftconvolve`), compares to what was actually observed, and back-projects the residual correction, reusing `_drizzle_matrix` (shared with the main drizzle resample loop) for the exact same per-frame affine mapping so the forward-simulate step is a literal inverse of the resample that built the initial estimate. Not compatible with `--elastic-registration` yet (IBP's forward model doesn't account for the local displacement field).
 
 - **Continuum-subtraction moments** (`src/channel_combine.py` `continuum_scale_moments`, used by `optimal_continuum_scale`, `--continuum`): the naive approach re-scans the full masked pixel array once per swept scale to compute `scipy.stats.skew` directly — this instead exploits that the subtraction residual `narrowband - s*continuum` is *linear* in `s`, so its skewness at any scale is a closed-form polynomial of 7 scalar central moments computed once, an algorithmic win independent of native vs. numpy (verified exact to float64 precision against `scipy.stats.skew` at several scales before replacing the per-scale loop). The native kernel is the constant-factor win on top of that: two passes (mean, then central moments, deliberately mirroring the numpy fallback's own two-step logic rather than a single-pass raw-moment shift-formula that would need separate algebra to get right) fusing what numpy computes as ~7 separate full-array elementwise-power passes (`ap*bp`, `ap*ap*bp`, ...) into one pass per stage. Not rayon-parallelised — called once per `optimal_continuum_scale` call, not once per swept scale, so there's no outer loop multiplying its cost the way this file's per-frame kernels have.

src/background.py

@@ -37,7 +37,7 @@
 except ImportError:
     # Fallback logging for standalone usage
     def safe_print(*args, **kwargs):
         print(*args, **kwargs)
 
 try:
     from astropy.stats import sigma_clipped_stats
@@ -80,6 +80,28 @@
     return gaussian_filter(a, sigma=sigma)
 
 
+def gaussian_blur_spatial(a: np.ndarray, sigma: float) -> np.ndarray:
+    """``gaussian_filter(a, sigma=(sigma, sigma, 0))`` for an ``(H, W, C)``
+    array -- blur each channel independently, never across channels --
+    routed through ``_gaussian_blur`` per channel instead of scipy's N-D
+    path. ``_gaussian_blur`` only accepts a scalar sigma on a 2-D plane, so
+    every one of this codebase's ``sigma=(s, s, 0)`` call sites (star
+    reduction's halo blur, HDR exposure-fusion's Laplacian pyramid) fell
+    back to scipy untouched even after ``_gaussian_blur`` existed. Same
+    per-channel-native-call pattern as ``postprocess.py``'s
+    ``_median_filter_per_channel``. 2-D input passes straight through to
+    ``_gaussian_blur``; anything else (no trailing channel axis to split)
+    falls back to plain ``gaussian_filter``."""
+    if a.ndim == 2:
+        return _gaussian_blur(a, sigma)
+    if a.ndim != 3:
+        return gaussian_filter(a, sigma=sigma)
+    out = np.empty_like(a)
+    for c in range(a.shape[2]):
+        out[:, :, c] = _gaussian_blur(a[:, :, c], sigma)
+    return out
+
+
 def gaussian_filter_ds(arr: np.ndarray, sigma: float,
                        ds_threshold: float = 24.0) -> np.ndarray:
     """Gaussian blur evaluated on a block-downsampled copy for large sigmas.

src/denoising.py

@@ -7,7 +7,12 @@
 from scipy import ndimage
 
 from src import wavelet
-from src.background import _estimate_sky_sigma, _gaussian_blur, gaussian_filter_ds
+from src.background import (
+    _estimate_sky_sigma,
+    _gaussian_blur,
+    gaussian_blur_spatial,
+    gaussian_filter_ds,
+)
 from src.models import Config
 from src.utils import get_logger, safe_print
 
@@ -737,9 +742,7 @@ def reduce_stars(
         return img
 
     # Blur per spatial dimension only (channel axis excluded)
-    blurred = ndimage.gaussian_filter(
-        img.astype(np.float64),
-        sigma=(blur_sigma, blur_sigma, 0))
+    blurred = gaussian_blur_spatial(img.astype(np.float64), blur_sigma)
 
     blend = (star_mask * reduction_factor).astype(np.float64)
     mask3 = blend[:, :, np.newaxis]

src/exposure_fusion.py

@@ -21,9 +21,9 @@
 from typing import List, Sequence
 
 import numpy as np
-from scipy.ndimage import gaussian_filter, zoom
+from scipy.ndimage import zoom
 
-from src.background import _gaussian_blur
+from src.background import _gaussian_blur, gaussian_blur_spatial
 
 
 def _fit_shape(x: np.ndarray, target_shape) -> np.ndarray:
@@ -43,8 +43,7 @@ def _fit_shape(x: np.ndarray, target_shape) -> np.ndarray:
 def _pyr_down(x: np.ndarray) -> np.ndarray:
     """Blur + downsample by 2 in the first two axes; a trailing channel
     axis, if present, is left untouched."""
-    sigma = (1.0, 1.0, 0.0) if x.ndim == 3 else (1.0, 1.0)
-    blurred = gaussian_filter(x, sigma=sigma, mode='reflect')
+    blurred = gaussian_blur_spatial(x, 1.0)
     return blurred[::2, ::2, ...]
 
 
@@ -57,8 +56,7 @@ def _pyr_up(x: np.ndarray, out_shape) -> np.ndarray:
     zoom_factors = (zy, zx, 1.0) if x.ndim == 3 else (zy, zx)
     up = zoom(x, zoom_factors, order=1)
     up = _fit_shape(up, out_shape)
-    sigma = (1.0, 1.0, 0.0) if x.ndim == 3 else (1.0, 1.0)
-    return gaussian_filter(up, sigma=sigma, mode='reflect')
+    return gaussian_blur_spatial(up, 1.0)
 
 
 def _gaussian_pyramid(x: np.ndarray, levels: int) -> List[np.ndarray]: