diff --git a/CHANGELOG.md b/CHANGELOG.md index 9fe8fe6..e2f262e 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CLAUDE.md b/CLAUDE.md index f565b35..ba82e88 100644 --- a/CLAUDE.md +++ b/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`). diff --git a/src/background.py b/src/background.py index e313e27..0d82bfb 100644 --- a/src/background.py +++ b/src/background.py @@ -1359,7 +1359,11 @@ def _fit_background_surface(coords: np.ndarray, values: np.ndarray, 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 @@ def _wavelet_background_surface(coords: np.ndarray, values: np.ndarray, 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) diff --git a/tests/test_dbe_gradient.py b/tests/test_dbe_gradient.py index 9205e4a..5a69944 100644 --- a/tests/test_dbe_gradient.py +++ b/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