diff --git a/CHANGELOG.md b/CHANGELOG.md index e2f262e..e821124 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CLAUDE.md b/CLAUDE.md index ba82e88..a1f2068 100644 --- a/CLAUDE.md +++ b/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. diff --git a/src/background.py b/src/background.py index 0d82bfb..bbba936 100644 --- a/src/background.py +++ b/src/background.py @@ -80,6 +80,28 @@ def _gaussian_blur(a: np.ndarray, sigma: float) -> np.ndarray: 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. diff --git a/src/denoising.py b/src/denoising.py index 5138601..97ce3b7 100644 --- a/src/denoising.py +++ b/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] diff --git a/src/exposure_fusion.py b/src/exposure_fusion.py index 3fdaf05..eef474f 100644 --- a/src/exposure_fusion.py +++ b/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]: diff --git a/tests/test_native.py b/tests/test_native.py index 628459d..593ee8c 100644 --- a/tests/test_native.py +++ b/tests/test_native.py @@ -724,6 +724,62 @@ def test_gaussian_blur_preserves_input_dtype(): np.testing.assert_allclose(out32.astype(np.float64), out64, rtol=1e-5, atol=1e-3) +def test_gaussian_blur_spatial_matches_scipy_per_axis_sigma(): + """gaussian_blur_spatial replaces this codebase's sigma=(s, s, 0) call + sites (star reduction's halo blur -- a default Phase 4 step -- and HDR + exposure-fusion's Laplacian pyramid), which _gaussian_blur alone can't + route to native since it only accepts a scalar sigma on a 2-D plane.""" + from scipy.ndimage import gaussian_filter + rng = np.random.default_rng(10) + a = rng.normal(1000.0, 50.0, (60, 70, 3)).astype(np.float64) + want = gaussian_filter(a, sigma=(2.5, 2.5, 0)) + got = _background_mod.gaussian_blur_spatial(a, 2.5) + np.testing.assert_allclose(got, want, rtol=1e-9, atol=1e-9) + + +def test_gaussian_blur_spatial_2d_matches_gaussian_blur(): + rng = np.random.default_rng(11) + a = rng.normal(size=(50, 60)) + np.testing.assert_array_equal( + _background_mod.gaussian_blur_spatial(a, 3.0), + _background_mod._gaussian_blur(a, 3.0)) + + +def test_gaussian_blur_spatial_falls_back_for_other_ndim(): + from scipy.ndimage import gaussian_filter + rng = np.random.default_rng(12) + a = rng.normal(size=(4, 5, 6, 2)) + want = gaussian_filter(a, sigma=1.5) + got = _background_mod.gaussian_blur_spatial(a, 1.5) + np.testing.assert_array_equal(got, want) + + +def test_reduce_stars_matches_original_scipy_three_axis_call(): + """reduce_stars (default Phase 4 star-reduction step, --no-star-reduce + to disable) switched its sigma=(blur_sigma, blur_sigma, 0) blur to + gaussian_blur_spatial -- pin its real output against the original raw + scipy call, not just the isolated gaussian_blur_spatial kernel.""" + from scipy.ndimage import gaussian_filter + + rng = np.random.default_rng(13) + H, W = 80, 100 + img = rng.uniform(500.0, 5000.0, (H, W, 3)).astype(np.float32) + star_mask = np.zeros((H, W), dtype=np.float32) + star_mask[30:36, 40:46] = 1.0 + star_mask[10:14, 70:74] = 0.6 + + blur_sigma = 1.5 + reduction_factor = 0.4 + ref_blurred = gaussian_filter(img.astype(np.float64), sigma=(blur_sigma, blur_sigma, 0)) + blend = (star_mask * reduction_factor).astype(np.float64) + mask3 = blend[:, :, np.newaxis] + ref = np.clip(img.astype(np.float64) * (1.0 - mask3) + ref_blurred * mask3, + 0.0, None).astype(np.float32) + + got = _denoising_mod.reduce_stars(img, star_mask, reduction_factor, blur_sigma) + np.testing.assert_allclose(got, ref, rtol=1e-5, atol=1e-3) + + def test_median_filter_per_channel_matches_combined_axis_scipy_call(): """postprocess.py's hot-pixel step switched from one scipy ndimage.median_filter(stacked, size=(5,5,1)) call (measured 5.1s on a