Switch DBE's regression-surface upsample to bilinear

CHANGELOG.md

@@ -8,6 +8,14 @@ match the `VERSION` file and `v*` git tags.
 
 ### Changed
 
+- **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
+  have added -- same reasoning already applied to `gaussian_filter_ds`'s own upsample. Validated directly
+  against real regression-grid output (post-blur mean diff <0.1, max <1.0 ADU), not assumed. A sibling call
+  in `wavelet_background_extraction` (no follow-up blur there) was left alone -- it has no real quality
+  test to validate a change against, only a mocked one, and a synthetic check was inconclusive.
+
 - **DBE's entropy filter runs a native kernel instead of a Python loop of ~4000 `np.histogram` calls.**
   Profiling a real `--auto` run (which sets `entropy_bg=True` for most target types) found this filter --
   documented as "cheap" in the code it lives next to -- costing 4.2s on its own, because it runs on every

CLAUDE.md

@@ -223,7 +223,7 @@ with `bayerPattern` to override.
 - **Online (streaming) sigma-clip** (`src/stacking.py`, used by `--stream`/`src/stream_stack.py`): `online_sigma_clip_seed_burnin` seeds a running Welford `(mean, m2, n_acc)` state from a small MAD-rejected burn-in window, and `online_sigma_clip_fold_frame` folds one further frame in at a time — the per-frame step of a true frame-at-a-time streaming combine, as an alternative to materializing the whole `(N,H,W,C)` aligned stack `sigma_clip_combine` needs. `online_sigma_clip_fold_frame` updates its `(mean, m2, n_acc)` state arrays **in place** (the one exception to this file's alloc-and-return convention) — reallocating and copying three full-resolution float64 arrays on every accepted frame was real allocator churn at full sensor resolution. `online_sigma_clip_combine` (whole-array, benchmark-only) is also available for throughput comparison against the batch kernels.
 - **Per-frame warp** (`src/registration.py` `apply_transform`): `warp_affine_lanczos3`, a 2D Lanczos-3 affine/shift resample (~5×). CPU path only (GPU unchanged). This is a quality-validated change, not numeric parity: FWHM and flux match scipy order-3, and the mild Lanczos ringing is incoherent across dithered frames (averages to zero in the stack).
 - **Anisotropic diffusion** (`src/denoising.py` `anisotropic_diffusion`): native Jacobi iteration with periodic boundary (~37×), exact parity.
-- **DBE surface fit + patch sampler** (`src/background.py`): `dbe_fit_surface`, Gaussian-weighted local-linear regression with Tukey-biweight IRLS (~2.4× over the scipy RBF it replaced), and `dbe_sample_patches`, the per-patch rejection cascade (~31×). Not just a speedup — the RBF (thin-plate spline + hard outlier-rejection loop) was unbounded and could extrapolate wildly into rejected-sample gaps near bright stars; the local fit is bounded near the sample values by construction. Numpy mirror in `_dbe_fit_surface_numpy`, parity-tested. **`patch_entropy_batch`**: `dbe_sample_patches`'s own docstring calls the entropy filter "cheap, operates on the small per-patch result" -- true per-call, but under `--auto` (`entropy_bg=True` for most target types) it runs on every DBE pass regardless of session size, and profiling a real run found `_filter_sampled_patches`'s Python loop over up to `Config.DBE_MAX_SAMPLES` (4000) candidate patches, each calling `np.histogram`, costing 4.2s on its own. Fused into one rayon-parallel native pass. Not bit-exact -- numpy's histogram fast path has a floating-point edge-correction pass (compares each value against the actual `linspace` bin edges and nudges the index where rounding put it a bin off) this doesn't replicate, so entropy values agree to ~3e-4 absolute on a real-shaped case, not bit-for-bit; immaterial for what this feeds (a median+MAD outlier threshold over thousands of patches). Cut a real profiled run's Post-process time by ~2.6s.
+- **DBE surface fit + patch sampler** (`src/background.py`): `dbe_fit_surface`, Gaussian-weighted local-linear regression with Tukey-biweight IRLS (~2.4× over the scipy RBF it replaced), and `dbe_sample_patches`, the per-patch rejection cascade (~31×). Not just a speedup — the RBF (thin-plate spline + hard outlier-rejection loop) was unbounded and could extrapolate wildly into rejected-sample gaps near bright stars; the local fit is bounded near the sample values by construction. Numpy mirror in `_dbe_fit_surface_numpy`, parity-tested. **`patch_entropy_batch`**: `dbe_sample_patches`'s own docstring calls the entropy filter "cheap, operates on the small per-patch result" -- true per-call, but under `--auto` (`entropy_bg=True` for most target types) it runs on every DBE pass regardless of session size, and profiling a real run found `_filter_sampled_patches`'s Python loop over up to `Config.DBE_MAX_SAMPLES` (4000) candidate patches, each calling `np.histogram`, costing 4.2s on its own. Fused into one rayon-parallel native pass. Not bit-exact -- numpy's histogram fast path has a floating-point edge-correction pass (compares each value against the actual `linspace` bin edges and nudges the index where rounding put it a bin off) this doesn't replicate, so entropy values agree to ~3e-4 absolute on a real-shaped case, not bit-for-bit; immaterial for what this feeds (a median+MAD outlier threshold over thousands of patches). Cut a real profiled run's Post-process time by ~2.6s. **`_fit_background_surface`'s coarse-grid upsample is bilinear (`order=1`), not cubic (`order=3`)**: the very next line blurs the surface again (`sigma=patch_size*0.5`, typically >=32px), same "nothing left for cubic's curvature term to recover" reasoning as `gaussian_filter_ds`'s own upsample -- validated directly (not assumed) against real `_dbe_regression_coarse_grid` output, post-blur mean diff <0.1, max <1.0 (`tests/test_dbe_gradient.py::test_fit_background_surface_bilinear_upsample_matches_cubic`). `wavelet_background_extraction`'s sibling `zoom(order=3)` call (no follow-up blur there) was deliberately left alone -- it has no real quality test (only a mocked one in `test_sky_model.py`) and a synthetic check was inconclusive, so there wasn't enough evidence to change it with confidence.
 - **Cosmic-ray rejection** (`src/stacking.py` `lacosmic_reject`): `lacosmic_reject_native`, L.A.Cosmic-style Laplacian spike detection + 5×5 median replacement, f32 internally (~2×+ vs numpy; the memory-traffic halving matters more than raw compute under 16-way ProcessPool contention).
 - **Median filter** (`median_filter_native`): windowed median (reflect boundary) used by lacosmic and the hot-pixel detectors. Interior fast path with contiguous reads; 3×3 uses Paeth's 19-op branchless median network (~13×), 5×5 quickselect (~1.6×).
 - **Malvar-He-Cutler debayer** (`src/debayer.py` `debayer_malvar`, `--debayer-method malvar`, the default): the true published algorithm (Malvar, He & Cutler 2004), not an approximation — sparse per-pixel tap gather (each output pixel needs at most 2 of the 4 kernels; its own channel is the raw sample) rather than 4 whole-image convolutions, with an interior fast path (no per-tap boundary-index modulo) for all but a 2px border (~2×). Replaces the old cv2 `COLOR_BAYER_*_EA` path, which required requantizing to uint16 first (real precision loss) and only worked when cv2 was installed; this runs on the native float32 data with no dependency. Numpy mirror (`_debayer_malvar_numpy`) validated bit-exact against the `colour-demosaicing` package's reference Malvar2004 implementation across all 4 Bayer patterns (validation-only dependency, not required at runtime — see `tests/test_debayer_malvar.py`).

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
@@ -1359,7 +1359,11 @@
     coarse, surf_lo, surf_hi = _dbe_regression_coarse_grid(
         coords, values, H, W, outlier_sigma, max_iter, sigma_px, Hc, Wc, verbose)
 
-    surface = zoom(coarse, (H / Hc, W / Wc), order=3)[:H, :W]
+    # order=1 (bilinear), not order=3 (cubic): the very next line blurs this
+    # surface again (sigma=patch_size*0.5, typically >=32px), which absorbs
+    # whatever cubic's extra curvature term would have added -- same
+    # reasoning as gaussian_filter_ds's own upsample (see its comment).
+    surface = zoom(coarse, (H / Hc, W / Wc), order=1)[:H, :W]
     surface = _gaussian_blur(surface, patch_size * 0.5)
     return np.clip(surface, surf_lo, surf_hi)
 
@@ -1598,7 +1602,6 @@
         safe_print(f"    Wavelet BG: {scales} scales on {Hc}x{Wc} regression grid "
                    f"(cutoff ~{(1 << scales) * stride}px)")
 
-    coarse_approx = np.clip(coarse_approx, surf_lo, surf_hi)
     surface = zoom(coarse_approx, (H / Hc, W / Wc), order=3)[:H, :W]
     return np.clip(surface, surf_lo, surf_hi)
 

tests/test_dbe_gradient.py

@@ -110,3 +110,41 @@ def test_gaussian_filter_ds_bilinear_upsample_matches_full_resolution():
     diff = np.abs(exact - fast)
     assert float(diff.mean()) < 0.5           # field std is ~20 -- this is noise-floor small
     assert float(diff.max()) < 5.0
+
+
+def test_fit_background_surface_bilinear_upsample_matches_cubic():
+    """_fit_background_surface's coarse-grid upsample switched cubic ->
+    bilinear too, on the same reasoning as gaussian_filter_ds's -- the very
+    next line blurs the surface again (sigma=patch_size*0.5), which absorbs
+    whatever cubic's extra curvature term would have added. Check the
+    post-blur surfaces stay close regardless of which interpolation order
+    fed into that blur (a real DBE run through dynamic_background_extraction
+    already exercises this function end-to-end; this pins the specific
+    order=1-vs-3 claim directly, on the same regression-grid machinery
+    _fit_background_surface itself uses, not a synthetic stand-in)."""
+    from scipy.ndimage import zoom
+
+    from src.background import _dbe_regression_coarse_grid, _gaussian_blur
+
+    rng = np.random.default_rng(7)
+    H, W = 420, 640
+    patch_size = 32
+    stride = max(4, patch_size // 4)
+    Hc, Wc = max(4, H // stride), max(4, W // stride)
+
+    n_pts = 500
+    coords = np.column_stack([rng.uniform(0, H, n_pts), rng.uniform(0, W, n_pts)])
+    values = (1000.0 + 40.0 * np.sin(coords[:, 0] / 80.0)
+              + 25.0 * np.cos(coords[:, 1] / 100.0) + rng.normal(0, 8.0, n_pts))
+    sigma_px = 1.25 * patch_size
+
+    coarse, surf_lo, surf_hi = _dbe_regression_coarse_grid(
+        coords, values, H, W, 2.5, 3, sigma_px, Hc, Wc, False)
+
+    blur_sigma = patch_size * 0.5
+    up1 = _gaussian_blur(zoom(coarse, (H / Hc, W / Wc), order=1)[:H, :W], blur_sigma)
+    up3 = _gaussian_blur(zoom(coarse, (H / Hc, W / Wc), order=3)[:H, :W], blur_sigma)
+
+    diff = np.abs(up1 - up3)
+    assert float(diff.mean()) < 0.1
+    assert float(diff.max()) < 1.0