From aa016c7c6fcc4cc4f2e62284a758d22e56f430bd Mon Sep 17 00:00:00 2001 From: Hans Davenport <35202271+hd152@users.noreply.github.com> Date: Tue, 22 Sep 2026 10:40:09 -0700 Subject: [PATCH 1/4] Fix real Phase 1-4 hotspots found by profiling a full run Profiling a full pipeline run (not guessing) found --flat-from-lights's synthetic-flat build silently dominating wall time -- 126.9s of a 187.8s total, hidden in an unlabeled "Other (I/O)" bucket because its timer was gated on real bias/dark/flat frames existing, which is exactly what this path runs without. - Fuse robust_pca_decompose's per-iteration elementwise arithmetic into two new native kernels (robust_pca_pre_svd_input, robust_pca_iterate) and route the L-update matmul through the existing small_times_wide kernel instead of np.matmul (this machine's numpy has no optimized BLAS). 126.9s -> 45.4s. - Time master-calibration building unconditionally, not just when real cal frames exist, and give it its own line in the timing summary instead of letting it vanish into "Other". - --flat-from-lights block-averages each of the 4 CFA sub-planes independently (never across colours) before decomposition, upsampling the recovered flat back to full res after -- vignetting/dust is smooth well above the pixel scale this removes, unlike a real dark/bias master's per-pixel hot pixels. 45.2s -> 3.6s. - DBE's compact-source dilation used scipy's binary_dilation with a large disk structuring element, pathologically slow on one contiguous bright source (16.2s on a synthetic case matching that shape); replaced with a distance-transform threshold, the same result computed near-linearly. - Post-process's hot-pixel step called scipy's median_filter directly instead of this project's own native per-channel median kernel -- same class of miss already fixed once elsewhere, missed here. 5.1s -> ~0.1s. - gaussian_filter_ds's large-sigma upsample switched cubic -> bilinear interpolation; the input just came out of the function's own huge blur, so cubic's extra curvature term recovers nothing real. ~3x/call. Same real session end to end: 3m 7.8s -> 36.9s (~5x). Tried and reverted: finer-grained rayon parallelism in small_times_wide measured no improvement (memory-bandwidth bound, not core bound) -- kept simple. Every change is bit-exact, measured-equivalent, or covered by a new native-vs-fallback parity test. New tests in test_native.py, test_robust_pca.py, test_dbe_gradient.py; tools/bench_robust_pca_scale.py and expanded tools/bench_native.py entries for the new kernels. astro_native bumped to 0.33.0. Full narrative in CLAUDE.md and CHANGELOG.md. Co-Authored-By: Claude Sonnet 5 --- CHANGELOG.md | 36 ++++++++ CLAUDE.md | 7 +- ext/astro_native/Cargo.lock | 2 +- ext/astro_native/Cargo.toml | 2 +- ext/astro_native/pyproject.toml | 2 +- ext/astro_native/src/lib.rs | 151 +++++++++++++++++++++++++++++++ src/background.py | 27 ++++-- src/cli.py | 20 +++-- src/io_fits.py | 10 ++- src/models.py | 53 +++++++++-- src/pipeline.py | 37 ++++++-- src/postprocess.py | 39 +++++++- src/robust_pca.py | 152 +++++++++++++++++++++++++++++--- tests/test_dbe_gradient.py | 51 +++++++++++ tests/test_native.py | 115 ++++++++++++++++++++++++ tests/test_robust_pca.py | 103 ++++++++++++++++++++++ tools/bench_native.py | 32 +++++++ tools/bench_robust_pca_scale.py | 123 ++++++++++++++++++++++++++ 18 files changed, 914 insertions(+), 48 deletions(-) create mode 100644 tools/bench_robust_pca_scale.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e5e216..b8820cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,42 @@ match the `VERSION` file and `v*` git tags. ## [Unreleased] +### Changed + +- **Several real perf fixes, found by profiling a full run rather than guessing.** A synthetic-flat build + (`--flat-from-lights`, auto-triggered whenever no flat frames exist) was silently the single biggest cost + on a real profiled session -- 126.9s of a 187.8s total, hidden inside an unlabeled "Other (I/O)" bucket in + the timing summary. Fixed in stages, each measured before moving to the next: + - `robust_pca_decompose`'s per-iteration numpy arithmetic (soft-threshold, residual, `Y` update, norm) was + ~12 full-array passes/iteration outside the (already-native) SVD step; fused into two native kernels + (`robust_pca_pre_svd_input`, `robust_pca_iterate`), plus routing the L-update's matrix reconstruction + through the existing `small_times_wide` kernel instead of a `np.matmul` call that hit this machine's + unoptimized reference BLAS. 126.9s -> 45.4s. + - The timing summary's "Other (I/O)" bucket was hiding master-calibration building specifically because its + timer was gated on real bias/dark/flat frames existing -- `--flat-from-lights` runs precisely when they + don't, so it was never timed at all. Now timed unconditionally as its own "Calibration" line. + - `--flat-from-lights` now block-averages each of the mosaic's 4 CFA sub-planes independently (never across + colours) by 4x before decomposition, upsampling the recovered flat back to full resolution after -- + a flat's vignetting/dust signal is smooth well above the pixel scale this removes, unlike a real dark or + bias master's per-pixel hot pixels, which is why this is `--flat-from-lights`-only. 45.2s -> 3.6s. + - DBE's compact-source mask dilation (`_build_emission_mask`) called `scipy.ndimage.binary_dilation` with a + large disk structuring element, pathologically slow on one large contiguous bright source (a real galaxy + or comet core): 16.2s on a synthetic case matching that shape. Replaced with a distance-transform + threshold -- the exact same result (verified bit-exact), computed in near-linear time instead. + - Post-processing's first step (per-channel hot-pixel removal) called `scipy.ndimage.median_filter` directly + instead of this project's own native per-channel median kernel -- the same class of miss already fixed + once elsewhere in the codebase, just never caught here. 5.1s -> ~0.1s (3 native calls, one per channel). + - `gaussian_filter_ds`'s large-sigma upsample (used throughout Phase 4: chroma denoising, local contrast, + DBE, and more) used cubic interpolation on a grid that had just come out of the function's own huge + Gaussian blur, where cubic's extra curvature term recovers nothing real; switched to bilinear, ~3x + faster per call, measured difference negligible. + + End to end, the same profiled real session: 3m 7.8s -> 36.9s (~5x). Every change is either bit-exact, + measured-equivalent (documented tolerance), or covered by a new native-vs-fallback parity test; one + attempted optimisation (finer-grained parallelism in `small_times_wide`) was measured, found to make no + difference, and reverted rather than shipped. Full detail, including what was tried and didn't help, in + [CLAUDE.md](CLAUDE.md). + ### Added - **Self-update check.** The CLI prints one line at the end of a run, and the desktop app shows a small diff --git a/CLAUDE.md b/CLAUDE.md index 33fa9b6..44365e7 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -67,7 +67,7 @@ The pipeline is split across `src/` modules. [originstack.py](originstack.py) is | [src/models.py](src/models.py) | `Config`, `FrameInfo`, `ProcessingStats` | | [src/utils.py](src/utils.py) | Print helpers, `format_time`, `get_memory_usage_mb`, `header_get_first` (first present/castable value among a list of FITS-header key spellings), `parse_timestamp` (ISO timestamp incl. a numeric UTC offset -> naive UTC `datetime`; the one home for the Celestron Origin `-0700` normalisation, shared by `sky_model.julian_date` and `observing_geometry`), `disable_astropy_network` (stops astropy reaching for the IERS tables over HTTP mid-run -- an offline machine otherwise eats a socket timeout per mirror inside functions that then fail soft; called from `cli.main` and `tests/conftest.py`) | | [src/io_fits.py](src/io_fits.py) | FITS load/save, `load_frame` (format dispatcher), `make_master`, `populate_fits_header` | -| [src/robust_pca.py](src/robust_pca.py) | Robust PCA (Principal Component Pursuit) master calibration frames (`--master-method robust_pca`): splits a bias/dark/flat stack into a low-rank component (the true shared pattern -- flat-field vignetting, fixed dark current) and a sparse component (dust motes that shifted between sessions, transient hot pixels), instead of `make_master`'s default per-pixel median treating every outlier independently. Opt-in only (needs >= `Config.ROBUST_PCA_MIN_FRAMES` frames, loads the full stack into memory, slower than median/mean); `--auto` narrowly auto-upgrades median->robust_pca per calibration type when its frame count falls in `[ROBUST_PCA_MIN_FRAMES, ROBUST_PCA_AUTO_MAX_FRAMES]` (gated per type since bias/dark/flat counts often differ). `--flat-from-lights` (`src/cli.py::_build_masters`) reuses this same decomposition on a capped sample of *light* frames when no dedicated flats exist -- the low-rank component approximates vignetting/dust, stars and nebula structure fall into the sparse component since dithering shifts them frame-to-frame while vignetting stays sensor-locked. Approximate; opt-in | +| [src/robust_pca.py](src/robust_pca.py) | Robust PCA (Principal Component Pursuit) master calibration frames (`--master-method robust_pca`): splits a bias/dark/flat stack into a low-rank component (the true shared pattern -- flat-field vignetting, fixed dark current) and a sparse component (dust motes that shifted between sessions, transient hot pixels), instead of `make_master`'s default per-pixel median treating every outlier independently. Opt-in only (needs >= `Config.ROBUST_PCA_MIN_FRAMES` frames, loads the full stack into memory, slower than median/mean); `--auto` narrowly auto-upgrades median->robust_pca per calibration type when its frame count falls in `[ROBUST_PCA_MIN_FRAMES, ROBUST_PCA_AUTO_MAX_FRAMES]` (gated per type since bias/dark/flat counts often differ). `--flat-from-lights` (`src/cli.py::_build_masters`) reuses this same decomposition on a capped sample of *light* frames when no dedicated flats exist -- the low-rank component approximates vignetting/dust, stars and nebula structure fall into the sparse component since dithering shifts them frame-to-frame while vignetting stays sensor-locked. Approximate; opt-in. **`--flat-from-lights` only: `downsample=Config.FLAT_FROM_LIGHTS_DOWNSAMPLE` (4)** (`robust_pca_master`, `bayer_block_downsample`/`bayer_block_upsample`) block-averages each of the mosaic's 4 CFA sub-planes independently (never across colours -- R/G/B have real gain differences, not just vignetting, so averaging the raw mosaic directly would corrupt the result) before decomposition, upsampling the recovered master back to full resolution after. Real bias/dark/flat robust_pca masters never pass this -- their sparse component is per-pixel hot pixels, which a downsample would blur away; a flat's low-rank content (vignetting, dust motes) is smooth well above the pixel scale a modest downsample removes. Cuts P by `downsample**2`: measured on the same real session behind the "profiling a real run" note two rows up, robust_pca_master 45.2s -> 3.6s (~12.6x, close to the 16x pixel-count reduction), full-pipeline wall clock on that session 1m26.3s -> 46.7s. Verified to still recover a real vignetting-plus-CFA-gain pattern at full output resolution, not a blurrier wrong one (`tests/test_robust_pca.py::TestRobustPcaMasterDownsample`, `TestBayerBlockDownsampleUpsample`) -- not yet validated against a real vignetted session end-to-end, only synthetic | | [src/dark_temp_model.py](src/dark_temp_model.py) | Temperature-interpolated dark current model (`--dark-temp-model`): fits a per-pixel low-order polynomial of dark signal vs. sensor temperature across the *whole* dark library in one vectorised `np.linalg.lstsq` call (design matrix in temperature, solved for every pixel at once), then evaluates it at the light frames' own temperature -- an alternative to `select_matching_darks`' nearest-temperature selection, letting a smaller library cover a wider temperature range via interpolation. Assumes the library is already homogeneous in ISO/gain/exptime (doesn't also model those); needs >= 3 distinct temperatures, falls back to nearest-match selection otherwise | | [src/io_raw.py](src/io_raw.py) | Camera RAW load (rawpy) — CR2/CR3/NEF/ARW/DNG/ORF/RW2/RAF/PEF/3FR/MRW/X3F/IIQ | | [src/io_tiff.py](src/io_tiff.py) | TIFF load (tifffile) — 16/32-bit linear TIFF lights | @@ -77,7 +77,7 @@ The pipeline is split across `src/` modules. [originstack.py](originstack.py) is | [src/debayer.py](src/debayer.py) | Debayering, hot pixels, white balance. **Session-constant CFA equalisation** (`cfa_frame_stats`/`combine_cfa_stats`/`set_session_cfa`; `frame_processor._measure_session_cfa`, `--no-session-cfa-eq` opts out): `green_equalize`'s G1/G2 gain and `_equalize_bayer_grid`'s four 2x2 green offsets were re-measured on every frame (six sigma-clipped medians over ~100 MB of strided reads -- ~70% of the Debayer step, and memory-bandwidth bound so it did not scale past the 8 physical cores: 248 ms/frame on 1 worker, 1040 on 16). Both are sensor/kernel properties and each per-frame estimate is itself noisy (~2 ADU against offsets of ~2-3 ADU), so eight frames spread through the session go through the real calibration (`_process_single_frame(cfa_probe=True)`), their medians are the session values, and every frame applies them (`_apply_fixed_grid`, in place, clipped at 0). Refused (per-frame fallback) below `Config.SESSION_CFA_MIN_FRAMES` (12) lights, on GPU, for non-Malvar debayer, with < 5 valid samples, or when samples disagree (grid std > 3 ADU / gain std > 0.01). The value lives in a debayer.py global set by the pool initialiser (workers) and `_prepare_session_cfa` (main process, remembered on `args._session_cfa` so `reload_accepted_frames` applies the same numbers) and is cleared by `@_with_session_cfa` when the phase ends. **Measured on the Omega session (114 frames): Debayer 1260 -> 447 ms/frame under 16 workers, Phase 1 39.3 -> 28.8 s, run 131 -> 114 s; stack differs from the per-frame path by 0.1 ADU robust rms against 95-164 ADU noise, FWHM and noise identical, residual 2x2 pattern 0.1-0.3 ADU in both.** The G1/G2 gain varies ~5x more frame-to-frame than noise explains (std 0.0019), so the per-frame gain is picking up something real; a fixed gain leaves a random +-8 ADU checkerboard per frame that averages out in the stack. Sky-level dependence was only tested over a ~10% range. A strided *subsample* of the medians is NOT a substitute: a period-4 structure in the data gave errors of 4-10 ADU against 4-6 ADU offsets | | [src/quality.py](src/quality.py) | `compute_quality_metrics`, star detection, FWHM, `estimate_bortle` (heuristic sky-glow bucket: calibrated background ADU normalised to exposure/gain, bucketed against log-spaced thresholds; same-equipment-only, not survey-grade -- see its docstring). Measured on the sessions behind the README/website sample images (real backyard/rooftop data, not a dark site): Bortle 7 for the two 500-gain sessions (Omega Nebula, Whirlpool Galaxy) and Bortle 9 for the four 200-gain sessions (Orion Nebula, Sagittarius Star Cloud, Andromeda Galaxy, Horsehead Nebula) -- the gain split lines up with a real difference in measured sky-glow rate (post-gain-normalisation ~60-65 vs. ~290-880 in the function's internal units), not just an artifact of the gain term, though with only 6 sessions and one heuristic that isn't fully verified apart from each other | | [src/psf_deconvolution.py](src/psf_deconvolution.py) | PSF estimation, Richardson-Lucy deconvolution | -| [src/background.py](src/background.py) | Mesh-based sky extraction, DBE, residual removal. **DBE re-admits gradient patches** (`_admit_gradient_patches`): the sampler drops any patch brighter than the *luminance-based* `sky_ref + 2*sky_std`, so in the strongest channel it dropped most patches and all the edge ones (measured: R 165 patches accepted vs ~715 for G/B; R edge excess +66..+114 ADU survived DBE untouched while G/B were removed completely). It now fits a 2nd/3rd-order polynomial to the accepted patches, admits candidates within 1 sky-sigma of it, refits and repeats (<= 8 rounds), so the model grows outward along a smooth gradient and stops at anything steeper; candidates come from the same sampler with the brightness cut off, so the emission mask and variance/entropy filters still apply. After the fix every edge is within +-3 ADU and the galaxy disk is unchanged or slightly brighter over sky (64.8 -> 70.0). This -- not a narrow strip -- was the coloured band along the JPG's bottom edge. `remove_edge_bands` (step 15a) stays as a partial backstop. When a caller supplies an `exclusion_mask` (`--galaxy-mode`/`--comet-mode`), `_dbe_prepare_emission_mask` first runs `_flatten_edge_glow`: components of the above-sky region larger than 2% of the frame whose smoothed peak lies in the 5% border band are a gradient cut off by the sensor edge (light pollution, a lamp, moonlight toward a corner), not an object, so they are set to sky before the emission detection passes and DBE models them as background. Without an exclusion mask the nebula-protecting default is unchanged (a nebula's brightest part can sit at an edge). Found on a real low-altitude Whirlpool session where a corner glow survived into the output | +| [src/background.py](src/background.py) | Mesh-based sky extraction, DBE, residual removal. **DBE re-admits gradient patches** (`_admit_gradient_patches`): the sampler drops any patch brighter than the *luminance-based* `sky_ref + 2*sky_std`, so in the strongest channel it dropped most patches and all the edge ones (measured: R 165 patches accepted vs ~715 for G/B; R edge excess +66..+114 ADU survived DBE untouched while G/B were removed completely). It now fits a 2nd/3rd-order polynomial to the accepted patches, admits candidates within 1 sky-sigma of it, refits and repeats (<= 8 rounds), so the model grows outward along a smooth gradient and stops at anything steeper; candidates come from the same sampler with the brightness cut off, so the emission mask and variance/entropy filters still apply. After the fix every edge is within +-3 ADU and the galaxy disk is unchanged or slightly brighter over sky (64.8 -> 70.0). This -- not a narrow strip -- was the coloured band along the JPG's bottom edge. `remove_edge_bands` (step 15a) stays as a partial backstop. When a caller supplies an `exclusion_mask` (`--galaxy-mode`/`--comet-mode`), `_dbe_prepare_emission_mask` first runs `_flatten_edge_glow`: components of the above-sky region larger than 2% of the frame whose smoothed peak lies in the 5% border band are a gradient cut off by the sensor edge (light pollution, a lamp, moonlight toward a corner), not an object, so they are set to sky before the emission detection passes and DBE models them as background. Without an exclusion mask the nebula-protecting default is unchanged (a nebula's brightest part can sit at an edge). Found on a real low-altitude Whirlpool session where a corner glow survived into the output. **`_build_emission_mask`'s compact-source dilation is a distance-transform threshold, not `scipy.ndimage.binary_dilation`**: dilation by a disk of radius r is exactly "within Euclidean distance r of a True pixel" (`ndimage.distance_transform_edt(1 - src_binary) <= r`, verified bit-exact against `binary_dilation(src_binary, structure=disk)`), and unlike scipy's generic morphology machinery this is near-O(H*W) regardless of mask shape. Found by profiling a real stack: scattered-mask synthetic timings looked fine (sub-second either way), but the actual case -- one large *contiguous* bright source (a real galaxy/comet core, not scattered stars) -- measured 16.2s for `binary_dilation` at r=40 on a full-res frame vs 0.37s for the distance-transform route; the profiled run's 2 calls cost 13.6s combined before the fix. **`gaussian_filter_ds`'s large-sigma upsample is bilinear (`order=1`), not cubic (`order=3`)**: profiling the same run found its `affine_transform` upsample step costing 3.5s across 10 calls; cubic's extra curvature term buys nothing on a grid that just came out of this same function's own huge Gaussian blur (measured order=1 vs order=3 mean/max diff 0.001/0.010 against a field std of ~0.75 on realistically-smoothed input) -- ~3x faster per call (0.336 -> 0.111s at ds=4 on a full-res field) | | [src/denoising.py](src/denoising.py) | Curvelet-style wavelet (`directional_wavelet_denoise`, the default primary denoiser), ACDNR, bilateral, anisotropic diffusion, star reduction, local contrast. `--variance-stabilize` applies a generalized Anscombe transform to the luma plane before wavelet thresholding, inverting it after -- makes BayesShrink's single per-subband noise estimate closer to valid everywhere, not just near the sky background level, since shot noise on bright pixels is Poisson not Gaussian; gain/read-noise are self-estimated from the image's own local mean-variance relationship, not sensor specs. `--denoiser curvelet` (`directional_wavelet_denoise`) is curvelet/shearlet-*inspired* (not an actual ridgelet/shearlet transform): locally reduces the BayesShrink threshold wherever a structure-tensor coherence map (`_structure_tensor_coherence`) detects elongated structure (filaments, galaxy arms), protecting it from thresholding more than an isotropic per-subband threshold would; **`--denoiser wavelet` is the one name** (`curvelet` is a hidden-in-spirit alias, normalised to `wavelet` in `parse_args`); how strongly structure is protected is its own option, `--wavelet-protect 0-1` (dest `directional_protect_strength`, default 0.6, 0 = plain BayesShrink; counts as explicit only when passed). Before 2026-09 `--denoiser wavelet` meant *protection forced to 0* -- old command lines using it now get 0.6; spell the old behaviour `--denoiser wavelet --wavelet-protect 0`. The post-processing step is skippable as `--skip-step wavelet` or `curvelet`. **Removed 2026-09 after `tools/bench_denoise_quality.py`** (ground-truth synthetic scenes, correlated noise matched to real stacks): NLM (halved star peak brightness, 11x the input error inside star cores, 17 s/MP), BM3D (best quality but slow and licence-encumbered; its pure-scipy fallback was worse than the pip package), MMT (erased ~93% of fine structure at its default strength), the non-adaptive `wavelet_denoise` + `--denoise-strength`/`--denoise-strength-calibrate` + the Noise2Self calibration module (only served that path), and `adaptive_wavelet_denoise` (identical to curvelet at protect 0). Old `--denoiser nlm/mmt/bm3d` values now error; `--denoiser wavelet` still works. **`reduce_chroma_noise` must not clip negatives**: it runs right after DBE, which centres the sky on zero, and the sky pedestal that keeps noise positive is applied later. A `np.clip(result, 0, None)` there half-wave-rectified the sky (exactly 50% zeros + positive spikes) and, because the pedestal is sized from the sky MAD, that MAD then read 5 ADU against a true ~100 (pedestal +39 instead of +738) -- the stretch rendered the surviving spikes as a field of white dots on a very low-SNR stack (real Whirlpool session, SNR ~1.7/sub) ([tests/test_corner_glow_and_thin_catalog.py](tests/test_corner_glow_and_thin_catalog.py)) | | [src/wavelet.py](src/wavelet.py) | `wavedec2`/`waverec2`: native 2D wavelet transform. bior1.3 (forward+inverse) used by `denoising.py`'s wavelet denoisers; db4 (forward-only) used by `quality.py`'s `compute_multiscale_entropy`. `soft_threshold`'s `value` accepts a per-pixel array (not just a scalar) -- needed by `directional_wavelet_denoise`'s spatially-adaptive threshold and `sparse_wavelet_deconvolve`'s FISTA proximal step | | [src/registration.py](src/registration.py) | `calculate_shift`, affine/RANSAC, `calc_common_crop`, `run_registration_phase`, `fit_displacement_field`. **`registration_stars`** rescues a thin star catalog: Phase 1's default matched-filter threshold returned 2-6 stars per frame on a noisy/hazy session (real Whirlpool, SNR ~1.7), below the 3 `match_stars_affine` needs, so affine matching declined for every frame and registration **silently** went translation-only -- field rotation (~2.3 deg over that session, ~60 px of smear at the edge) then arced stars around the field centre, and the residual check (needs >= 5 reference stars) never ran, so nothing warned. When a catalog is under `Config.REG_MIN_STARS` (12) it re-detects on a 2x2-binned luminance at k_confirm 12/8/6 and keeps the larger catalog (coordinates mapped back to the input grid); used for the reference, each frame, the consensus-reference swap and the residual check's re-detection (they must share one detector or the residual compares mismatched catalogs). `find_extended_source_ellipse` (`--galaxy-mode`) also rejects blobs whose smoothed *peak* lies in the border band, not just those with an edge-near centroid: a 590k-px corner glow with a centroid well inside the frame beat the real galaxy (5k px) and the exclusion ellipse then protected the glow from background extraction ([tests/test_corner_glow_and_thin_catalog.py](tests/test_corner_glow_and_thin_catalog.py)) | @@ -96,7 +96,7 @@ The pipeline is split across `src/` modules. [originstack.py](originstack.py) is | [src/moving_objects.py](src/moving_objects.py) | `--moving-objects[-stack]`: residual vs the stack -> matched-filter smooth -> candidates outside brightness-scaled star masks -> `find_tracks` (velocity-space voting for a cell collecting detections from many *different frames*; v_min = 2 FWHM / span, default v_max 3 px/min) -> least-squares refit -> `tracked_stack` (window median along the track). Hooked in `run_stacking_phase` before the aligned memmap is deleted. Real session: 1497 candidate detections -> 0 tracks (no mover in that field; no false tracks from clutter) | | [src/lightcurve_analysis.py](src/lightcurve_analysis.py) | `--lightcurve-analysis` on `--photometry-timeseries` CSVs: astropy `LombScargle` (periods 4 cadences to half the baseline, with FAP) and `BoxLeastSquares` to place a dip, then a trapezoid `least_squares` fit with Jacobian errors and a BIC test against flat. Uses `astropy.timeseries` (verify it survives a PyInstaller build with `packaging/verify_build.ps1` before relying on it in the packaged app) | | [src/frame_processor.py](src/frame_processor.py) | Parallel workers, `execute_frame_processing`, `quality_gate` | -| [src/postprocess.py](src/postprocess.py) | Full post-processing chain: `postprocess_stack` | +| [src/postprocess.py](src/postprocess.py) | Full post-processing chain: `postprocess_stack`. **Step 1 (hot pixel removal) routes its 5x5 median filter through native per-channel calls** (`_median_filter_per_channel`, `median_filter_native`): the original single `scipy.ndimage.median_filter(stacked, size=(5,5,1))` call -- the first thing Phase 4 does, on every stack -- measured 5.1s on a real full-res frame; scipy's generic N-D rank-filter machinery has no fast path for a size-1 axis, so it silently bypassed this codebase's own already-existing native 2D median kernel (the same class of miss as `_fix_hot_bayer`'s pre-native-routing bug in `debayer.py`, just never caught here). 3 independent native per-channel calls are equivalent by construction and validated against the combined-axis scipy call in `tests/test_native.py` | | [src/difference_imaging.py](src/difference_imaging.py) | Proper image subtraction and transient detection (`--transient-detect REF.fits`). Answers "did anything change?" rather than "what does my target look like?" -- novae, outbursts, supernovae, asteroids, variables. `zogy()` implements Zackay, Ofek & Gal-Yam (2016): rather than degrading one epoch to match the other (Alard-Lupton), it cross-convolves each image with the **other's** PSF, so both sides acquire the same effective PSF and stellar residuals cancel in closed form even across differing seeing -- a plain subtraction leaves a dipole at every star, scaling with brightness, i.e. worst exactly where transients hide. Returns `D` (difference), `S` (match-filtered score) and `S_corr` (score in units of its own propagated sigma, so a threshold is a real significance). **`S_corr`'s astrometric noise term is not optional**: a sub-pixel registration slip leaves a residual proportional to the local gradient, largest at bright stars, and with `astrometric_sigma=0` every bright star reports as a transient (asserted in [tests/test_difference_imaging.py](tests/test_difference_imaging.py)). `_prepare_psf` zero-pads and **rolls the PSF so its centre lands on index [0,0]** -- the FFT's origin; omitting the roll shifts every output by half the frame. `estimate_background_sigma` uses *symmetric* iterative sigma clipping even though stars are a one-sided contaminant: the obvious "keep pixels below the 80th percentile, take their MAD" truncates the Gaussian core and reads 5.7 against an injected 7.0, a 19% underestimate -- and that sigma is the denominator of every significance, so underestimating it manufactures false positives at exactly the threshold users trust. **Runs on the LINEAR stack** (`fits_stacked` in `pipeline.py`, the `RAWSTACK` product the output FITS contains), never the post-processed array: Phase 4's nonlinear stretches/denoise/local-contrast break photometric linearity, and comparing a post-processed frame against a linear reference mismatches the flux scale by a fraction of a percent -- several sigma on a bright star, reporting a transient at every star in the field. **`--transient-detect` refuses a reference without `RAWSTACK=True`** (same rule as `--merge`). The warped reference's uncovered area -- field rotation leaves empty corners -- is masked via a footprint (`_align_reference` warps an all-ones image too; `_erode` uses `border_value=1` so the frame's own edge is not treated as uncovered, which measured 66% "covered" for a 93%-covered frame): unmasked, stars in those corners came out as confident 'brightenings' (6 false candidates on a 9 deg synthetic rotation, 0 masked). The astrometric sigma is the *measured* matched-star RMS residual with `_ASTROMETRIC_SIGMA_FLOOR_PX = 0.3` as a floor -- it used to be a hard-coded 0.3 reported as a measurement. `zogy` pads to `scipy.fft.next_fast_len` (3.4x on the Origin's 1096 = 8x137 axis, whose prime factor pushes pocketfft onto Bluestein; `S_corr` identical to 2e-5 sigma in the interior). `detect_transients` is a single sorted pass, verified identical to the old per-candidate argmax loop on 300 randomized fields including ties and NaNs. `psf_difference` was dropped from `ZogyResult` (no reader). The catalogue's WCS must be built with `naxis=2` (`pipeline.py`): a bare `WCS(header)` on the `(3, H, W)` cube is 3-axis, still passes `has_celestial`, and blanks every RA/Dec | | [src/sky_model.py](src/sky_model.py) | **EXPERIMENTAL** physics-based sky background model (`--bg-method physical`). Fits `sky = c0 + c_air*airglow(z) + c_moon*moonlight(rho,z) + c_zodi*zodiacal(lambda,beta) + c_lp*skyglow(az,z)`, where every component's *spatial shape* is fixed by geometry and only a scalar amplitude is free, so unlike mesh/DBE/wavelet it has nowhere to put a nebula. **Real-data verdict: it does not work on a typical ~1 deg deep-sky field, and now detects that and declines.** On a real Lagoon session the zenith angle varies by only 0.94 deg across the whole frame and the azimuth by 1.5 deg, so every component map is essentially flat, the fit has nothing to grip, and subtracting it made the corner-to-corner gradient *worse* (67 -> 111 ADU) where DBE removed 68%; sweeping `--light-pollution-azimuth` through all 360 deg moved the residual by under 0.01 ADU. Structural, not a tuning problem: physical components vary on ten-degree scales, so a narrow field's gradient is dominated by *instrumental* effects (vignetting, amp glow, filter gradients) a sky model cannot represent. `remove_physical_sky` therefore measures `_corner_gradient` before and after and returns None unless it improved things, letting `postprocess.py::_apply_physical_sky` fall back to DBE with an accurate reason (it must distinguish 'no GPS' from 'could not help' -- reporting the former on a session that had full GPS sent the reader hunting for metadata that was already present). The guard is empirical, not a field-size rule, and its corner-based proxy is weakest on radially symmetric gradients whose four corners are equal by construction. **Two earlier claims are corrected by this testing**: DBE does *not* eat nebulosity (99.5% retained on the Lagoon -- the damage that motivated this work came from the `sky_residual` residual passes, a different step `--auto` already skips), and the '98% synthetic nebula preserved' figure holds only for a gradient built from the model's own basis. Coefficients are constrained non-negative (`scipy.optimize.nnls`) and that *is* load-bearing: across a real field the maps are nearly collinear (condition ~1e5), so an unbounded fit builds an interior maximum from cancelling coefficients and inverts a nebula (-25% preservation). Clipping is upward-only (stars sit above sky). Ephemerides are closed-form, not astropy: `EarthLocation`/`AltAz` pull in the IERS tables [packaging/originstack.spec](packaging/originstack.spec) excludes -- validated against astropy at **0.009 deg (sun) / 0.05 deg (moon)**. `julian_date` honours timezone offsets: a Celestron Origin `info.json` stamps local time (`2026-08-31T20:40:32-0700`), which the first version failed to parse at all -- silently disabling the model on exactly its target data -- and which loses 7 hours of moon position if the offset is merely stripped. Positions are mean-equinox-of-date while WCS pixels are J2000, so separations carry a ~0.35 deg precession offset (uncorrected, deliberately; and worth knowing, since that offset *looks* like an ephemeris bug -- the real one found this way was 19.6 deg, from Schlyter's lunar elements being epoched at 1999-12-31.0 rather than J2000.0). `describe_fit` gates component attribution on the basis condition number: removal can be valid while the *split* between components is unidentifiable. **Later changes:** the FOV gate now runs *before* `build_geometry` (a telescope field is declined in ~0 s / 0 MB instead of ~1.3 s / 578 MB); `fit_sky_model` fits on a strided sample (`_FIT_MAX_SAMPLES`, ~100k pixels; 65x faster, dominant coefficient within 0.25%) but evaluates the model at full resolution; `remove_physical_sky_with_reason` returns `(result, reason)` with the *real* reason (narrow field / failed fit / no improvement / no geometry) and `remove_physical_sky` is a thin wrapper; `_nnls` has **no fallback** -- it used to catch every exception (incl. scipy's max-iteration `RuntimeError`) and substitute a clamped `lstsq`, i.e. the unbounded fit this module documents as inverting nebulae, and the guard cannot see radially symmetric damage. A failure now declines to DBE with a warning. Azimuth and helio-ecliptic longitude are upsampled via sin/cos (`_up_angle`): interpolating across the 360/0 seam ramped the long way round, measured as a 22.5 deg/px step where the truth is ~0.05. The unreleased multi-frame joint fit and `amp_glow_basis` were removed -- nothing in `src/` called them | | [src/uncertainty.py](src/uncertainty.py) | End-to-end uncertainty propagation (`--uncertainty-propagate`) and confidence mapping. Phase 3's `--uncertainty-map` sigma describes the *linear* stack; Phase 4 then reshapes the noise field, so this carries the error bars through to the delivered image. `propagate_uncertainty` is **Monte Carlo, not analytic**: it draws K noise realizations (`--uncertainty-realizations`, default 8) from the Phase 3 sigma map and pushes each through the *unmodified* `postprocess_stack`, taking the per-pixel spread. Deliberate — nine Phase 4 steps are nonlinear denoisers and four are iterative deconvolvers (RL, FISTA, anisotropic diffusion), several spatially adaptive (BayesShrink thresholds, the structure-tensor coherence map), so no closed-form Jacobian exists for most of the chain; the MC estimator is exact up to `~1/sqrt(2K)` *for a chain that does not adapt to its input's noise level*, and stays correct automatically when a denoiser is added. **It is slightly biased low for adaptive steps:** each realization is `stacked + N(0, sigma)` but `stacked` already carries ~sigma, so Phase 4 sees sqrt(2)x the real noise and steps that estimate parameters from the data (BayesShrink, DBE sky sigma) denoise it harder. Measured against this project's own `wavelet_denoise` over sigma {1,4,12} x threshold {2,3,5}: ratio 0.93-1.03 -- a few percent, inside the ~25% MC error at K=8. A first-principles argument predicted ~2x understatement; it did not reproduce (shrinkage and threshold move together), and a reviewer's toy-chain figure was likewise wrong -- measure against the real denoiser before believing either. `propagate_uncertainty` returns `(sigma_post, mean_post, adaptivity)`, where `adaptivity` is sigma(full amp)/sigma(half amp): 2.0 = scale-invariant, and `pipeline.py` warns below 1.5. Realizations run under `_quiet_args` (file-writing/network steps off — `_QUIET_OFF`: `remove_stars`, `nmf_separate`, `photometric_calibration`, `annotate`, `aberration_report`, `diagnostic`, `export_masks`, `keep_intermediates`, `comet_radial_renorm`, `comet_larson_sekanina`, and `denoise_strength_calibrate`, a 9-point sweep. **Any Phase 4 step that writes a file or calls out belongs in that list**: every realization runs under `redirect_stdout`, so a step left on prints its own "Saved:" into a swallowed buffer and the file left on disk is the *last noise realization*, silently) with stdout swallowed, so K passes don't emit K sets of sidecars or K Gaia queries; everything that shapes the *noise* is left exactly as the real pass ran it. `confidence_map` turns the propagated sigma into per-pixel SNR above sky and returns **`NaN` where the propagated sigma is exactly zero** — those pixels were clamped to a constant by the chain (sky pedestal lift, non-negativity clips), so they carry no measurement, and dividing by ~0 would otherwise rank the pipeline's own floor artifacts as the most confident pixels in the frame (observed at ~25% of a real synthetic frame). `error_aware_black_point` (`--error-aware-stretch SIGMA`) returns the faintest pixel clearing SIGMA confidence, used as the preview black point so sub-threshold content clips to black instead of being stretched into apparent structure. Validated against chains with known variance transformation — identity recovers the input sigma, a x3 scale scales it x3, a 3x3 box blur divides it by exactly 3 ([tests/test_uncertainty_propagation.py](tests/test_uncertainty_propagation.py)) | @@ -233,6 +233,7 @@ with `bayerPattern` to override. - **Blind (unknown-rotation) rigid star-pattern match** (`src/blind_match.py` `match_rigid_unknown_rotation`, used by `--merge`'s cross-night registration and the in-session fallback when near-zero-rotation RANSAC fails): the hypothesis-scoring loop (up to `_MAX_HYPOTHESES`=20000 distance-matched pair-of-pairs, each requiring a nearest-neighbour consensus count) ported as a brute-force O(n·m) nearest-neighbour search — n/m are capped small (`max_stars`, default 40) by the caller, so what this actually removes is the per-hypothesis Python-loop + scipy `cKDTree.query` call overhead, not the underlying query cost, which is trivial at this scale either way. The final all-inliers Umeyama refit stays in Python (`src/affine_fit.py::_umeyama_2d`), called once per match, not per-hypothesis. Numpy mirror (`_match_hypotheses_numpy`) parity-tested against the native kernel on synthetic rotated star fields (`tests/test_native.py`). - **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`. - **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/ext/astro_native/Cargo.lock b/ext/astro_native/Cargo.lock index 8d448c2..298b8d1 100644 --- a/ext/astro_native/Cargo.lock +++ b/ext/astro_native/Cargo.lock @@ -49,7 +49,7 @@ checksum = "fb5dfbc6d8d2675589ccbe4d0fd61df2419075625f8c1a62325e718e2b0049f9" [[package]] name = "astro_native" -version = "0.32.0" +version = "0.33.0" dependencies = [ "numpy", "pyo3", diff --git a/ext/astro_native/Cargo.toml b/ext/astro_native/Cargo.toml index af9ee48..0f04e3c 100644 --- a/ext/astro_native/Cargo.toml +++ b/ext/astro_native/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "astro_native" -version = "0.32.0" +version = "0.33.0" edition = "2021" description = "Native (Rust) hot-path kernels for OriginStack: stacking combine, etc." diff --git a/ext/astro_native/pyproject.toml b/ext/astro_native/pyproject.toml index a865297..96aeb8c 100644 --- a/ext/astro_native/pyproject.toml +++ b/ext/astro_native/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "astro_native" -version = "0.32.0" +version = "0.33.0" description = "Native Rust hot-path kernels for OriginStack" requires-python = ">=3.10" classifiers = ["Programming Language :: Rust"] diff --git a/ext/astro_native/src/lib.rs b/ext/astro_native/src/lib.rs index e25a919..cad77a4 100644 --- a/ext/astro_native/src/lib.rs +++ b/ext/astro_native/src/lib.rs @@ -4552,6 +4552,16 @@ fn small_times_wide<'py>( // naive col-outer/i-inner order would stride by P between consecutive // reads -- exactly the huge-stride-read-stream antipattern this file's // gather-transpose driver (`row_parallel`) exists to avoid elsewhere. + // + // Tried and reverted: nested rayon parallelism (also column-chunking + // within each row) on the theory that N-way row parallelism leaves cores + // idle when N (a frame/calibration-stack count, often ~10-20) is well + // under the machine's core count. Measured no improvement (a real + // robust-PCA call: 0.41-0.47s/call either way) -- this operation is + // memory-bandwidth-bound on this machine, same lesson `gram_matrix_wide` + // above already documents for the sibling GEMM in this same trick + // (`numpy's own @ got zero benefit from this machine's cores at this + // shape`). Kept simple since the added chunking bought nothing real. let mut out = vec![0f64; n * p]; py.detach(|| { out.par_chunks_mut(p).enumerate().for_each(|(k, out_row)| { @@ -4571,6 +4581,145 @@ fn small_times_wide<'py>( Ok(arr2.into_pyarray(py)) } +/// Fused `D - S + Y/mu`: the IALM L-update's per-iteration input to the thin +/// SVD in `_thin_svd_wide` (`src/robust_pca.py::robust_pca_decompose`). The +/// numpy reference builds this as two full-`(N,P)`-array passes (`Y/mu` into +/// a temporary, then `D - S + that` into another); profiling a real +/// `--flat-from-lights` run (N=10, P=6.25M) found `robust_pca_decompose` +/// spending 100s of its 126s total in exactly this kind of elementwise numpy +/// arithmetic around the SVD, not in the SVD itself -- this and +/// `robust_pca_iterate` below fuse that arithmetic into one native pass each, +/// same category of win as `calibrate_frame_inplace`. Same float64 operation +/// order as the numpy reference (divide, then left-to-right subtract/add), so +/// the result is bit-identical. +#[pyfunction] +fn robust_pca_pre_svd_input<'py>( + py: Python<'py>, + d: PyReadonlyArray2<'py, f64>, + s: PyReadonlyArray2<'py, f64>, + y: PyReadonlyArray2<'py, f64>, + mu: f64, +) -> PyResult>> { + let d_arr = d.as_array(); + let shape = d_arr.shape().to_vec(); + if s.as_array().shape() != d_arr.shape() || y.as_array().shape() != d_arr.shape() { + return Err(pyo3::exceptions::PyValueError::new_err( + "d, s, y must all have the same shape", + )); + } + let d_flat: &[f64] = d_arr + .as_slice() + .ok_or_else(|| pyo3::exceptions::PyValueError::new_err("d must be C-contiguous"))?; + let s_arr = s.as_array(); + let s_flat: &[f64] = s_arr + .as_slice() + .ok_or_else(|| pyo3::exceptions::PyValueError::new_err("s must be C-contiguous"))?; + let y_arr = y.as_array(); + let y_flat: &[f64] = y_arr + .as_slice() + .ok_or_else(|| pyo3::exceptions::PyValueError::new_err("y must be C-contiguous"))?; + + const CH: usize = 1 << 16; + let mut out = vec![0f64; d_flat.len()]; + py.detach(|| { + out.par_chunks_mut(CH) + .zip(d_flat.par_chunks(CH)) + .zip(s_flat.par_chunks(CH)) + .zip(y_flat.par_chunks(CH)) + .for_each(|(((oc, dc), sc), yc)| { + for i in 0..oc.len() { + oc[i] = dc[i] - sc[i] + yc[i] / mu; + } + }); + }); + + let arr2 = numpy::ndarray::Array2::from_shape_vec((shape[0], shape[1]), out) + .expect("shape mismatch building robust_pca_pre_svd_input output"); + Ok(arr2.into_pyarray(py)) +} + +/// Fused IALM S-update + residual + Y-update + residual-norm: the second half +/// of `robust_pca_decompose`'s per-iteration elementwise arithmetic (see +/// `robust_pca_pre_svd_input` above for the profiling context). The numpy +/// reference computes `temp = D - L + Y/mu`, `S = sign(temp) * +/// max(|temp| - lam/mu, 0)`, `residual = D - L - S`, `Y += mu * residual`, +/// then `norm(residual, 'fro')` as roughly a dozen separate full-array passes +/// (each its own temporary); this does all of it in one pass, writing `S` and +/// `Y` in place (`np.zeros_like(D)`-allocated once by the caller, reused +/// every iteration -- no repeated `(N,P)` allocation) and returning the +/// residual's Frobenius norm directly. Per-element arithmetic matches numpy's +/// operation order exactly (`np.sign` semantics: 0 at exactly 0, not +/// `f64::signum`'s +-1), so `S`/`Y` are bit-identical to the numpy reference; +/// the norm is a parallel reduction over chunks and so may differ from +/// numpy's sequential sum by a few ULPs, same as this file's other +/// parallel-reduction kernels -- immaterial here since it only feeds a +/// convergence-tolerance comparison. +#[pyfunction] +fn robust_pca_iterate<'py>( + py: Python<'py>, + d: PyReadonlyArray2<'py, f64>, + l: PyReadonlyArray2<'py, f64>, + mut s: numpy::PyReadwriteArray2<'py, f64>, + mut y: numpy::PyReadwriteArray2<'py, f64>, + lam_over_mu: f64, + mu: f64, +) -> PyResult { + let d_arr = d.as_array(); + let l_arr = l.as_array(); + if l_arr.shape() != d_arr.shape() + || s.as_array().shape() != d_arr.shape() + || y.as_array().shape() != d_arr.shape() + { + return Err(pyo3::exceptions::PyValueError::new_err( + "d, l, s, y must all have the same shape", + )); + } + let d_flat: &[f64] = d_arr + .as_slice() + .ok_or_else(|| pyo3::exceptions::PyValueError::new_err("d must be C-contiguous"))?; + let l_flat: &[f64] = l_arr + .as_slice() + .ok_or_else(|| pyo3::exceptions::PyValueError::new_err("l must be C-contiguous"))?; + let s_flat: &mut [f64] = s + .as_slice_mut() + .map_err(|_| pyo3::exceptions::PyValueError::new_err("s must be C-contiguous"))?; + let y_flat: &mut [f64] = y + .as_slice_mut() + .map_err(|_| pyo3::exceptions::PyValueError::new_err("y must be C-contiguous"))?; + + const CH: usize = 1 << 16; + let sq_sum: f64 = py.detach(|| { + d_flat + .par_chunks(CH) + .zip(l_flat.par_chunks(CH)) + .zip(s_flat.par_chunks_mut(CH)) + .zip(y_flat.par_chunks_mut(CH)) + .map(|(((dc, lc), sc), yc)| { + let mut acc = 0.0f64; + for i in 0..dc.len() { + let temp = dc[i] - lc[i] + yc[i] / mu; + let abs_temp = temp.abs(); + let shrunk = (abs_temp - lam_over_mu).max(0.0); + let sign = if temp > 0.0 { + 1.0 + } else if temp < 0.0 { + -1.0 + } else { + 0.0 + }; + let sval = sign * shrunk; + sc[i] = sval; + let resid = dc[i] - lc[i] - sval; + yc[i] += mu * resid; + acc += resid * resid; + } + acc + }) + .sum() + }); + Ok(sq_sum.sqrt()) +} + /// Single-pass central moments of a masked (narrowband, continuum) pixel /// pair, for `optimal_continuum_scale`'s closed-form skewness-vs-scale /// polynomial (`src/channel_combine.py`): the subtraction residual @@ -7548,6 +7697,8 @@ fn astro_native(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(warp_affine_kernel_table, m)?)?; m.add_function(wrap_pyfunction!(gram_matrix_wide, m)?)?; m.add_function(wrap_pyfunction!(small_times_wide, m)?)?; + m.add_function(wrap_pyfunction!(robust_pca_pre_svd_input, m)?)?; + m.add_function(wrap_pyfunction!(robust_pca_iterate, m)?)?; m.add_function(wrap_pyfunction!(continuum_scale_moments, m)?)?; m.add_function(wrap_pyfunction!(fit_moffat_native, m)?)?; m.add_function(wrap_pyfunction!(fit_psf_moffat2d_native, m)?)?; diff --git a/src/background.py b/src/background.py index 8fc9563..860a731 100644 --- a/src/background.py +++ b/src/background.py @@ -7,7 +7,7 @@ import numpy as np from scipy import ndimage from scipy.interpolate import RectBivariateSpline -from scipy.ndimage import binary_dilation, gaussian_filter, zoom +from scipy.ndimage import gaussian_filter, zoom try: import astro_native as _native @@ -81,10 +81,19 @@ def gaussian_filter_ds(arr: np.ndarray, sigma: float, # which multiplied by the field's gradient near bright features is the # dominant error term. affine_transform with the matching offset samples # the coarse grid at the true block centres. + # + # order=1 (bilinear), not order=3 (cubic spline): measured 3x faster + # (0.336 -> 0.111s/call on a real full-res field, ds=4) and, on a coarse + # grid that's already been through this function's own huge Gaussian + # blur (the only kind of input this ever upsamples), order=1 vs order=3 + # differ by a mean 0.001 / max 0.010 against a field std of ~0.75 -- + # negligible, for the same reason the module docstring gives for the + # downsample itself: content this smooth has nothing for cubic's extra + # curvature term to recover that bilinear doesn't already get. off = -(ds - 1) / (2.0 * ds) out = ndimage.affine_transform( sm, np.array([1.0 / ds, 1.0 / ds]), offset=[off, off], - output_shape=(H, W), order=3, mode='nearest') + output_shape=(H, W), order=1, mode='nearest') return out @@ -831,9 +840,6 @@ def _build_emission_mask(lum: np.ndarray, star_mask: Optional[np.ndarray], if frac_bright > 0.005: dil_radius = max(15, int(min(H, W) * 0.02)) r = dil_radius - y_idx, x_idx = np.ogrid[-r:r + 1, -r:r + 1] - disk = (y_idx ** 2 + x_idx ** 2 <= r ** 2) - structure = disk.astype(np.uint8) # Binary dilation expects struct array remaining_lum = lum_smooth.copy() primary_peak = float(np.max(remaining_lum)) @@ -852,8 +858,15 @@ def _build_emission_mask(lum: np.ndarray, star_mask: Optional[np.ndarray], src_thresh = sky_med + 0.5 * (detect_thresh - sky_med) src_binary = (lum_smooth > src_thresh).astype(np.uint8) - # Use binary_dilation on the binary mask - dilated = binary_dilation(src_binary, structure=structure).astype(np.float32) + # Dilation by a disk of radius r == "within Euclidean distance r + # of a True pixel", so a distance transform gives the exact same + # mask as scipy.ndimage.binary_dilation(src_binary, structure=disk) + # (verified bit-exact) -- and does it in ~O(H*W), not O(H*W*r^2). + # binary_dilation's generic morphology path is fine for a + # scattered mask but measured 16s on one real contiguous bright + # source at this frame size/radius (a big galaxy/comet core, not + # a handful of stars); the EDT route was 0.37s on the same mask. + dilated = (ndimage.distance_transform_edt(1 - src_binary) <= r).astype(np.float32) np.clip(emission + dilated, 0.0, 1.0, out=emission) # Blank out processed source diff --git a/src/cli.py b/src/cli.py index e41e90c..5efd928 100644 --- a/src/cli.py +++ b/src/cli.py @@ -366,7 +366,13 @@ def _build_masters(frames: dict, stats: "ProcessingStats | None" = None, """ lights = frames.get('light', []) cal_needed = frames.get('dark') or frames.get('flat') or frames.get('bias') - cal_start = time.time() if cal_needed else None + # Timed unconditionally (not just when cal_needed): --flat-from-lights + # builds a synthetic flat from *light* frames when no real flat/dark/bias + # exist at all, so cal_needed is False on exactly that path -- gating the + # timer on it meant the slowest thing this function can do (a robust_pca + # decomposition of sampled lights) was invisible to stats.calibration_time + # and fell into the caller's opaque "Other" bucket instead. + _build_start = time.time() if cal_needed: safe_print("\nCreating master calibration frames...") @@ -483,11 +489,13 @@ def _method_tag(method: str) -> str: sample = sample[:Config.ROBUST_PCA_AUTO_MAX_FRAMES] safe_print(f" No flat frames found -- deriving a synthetic flat from " f"{len(sample)} light frames (--flat-from-lights)...") - synthetic_flat = make_master(sample, method='robust_pca') + synthetic_flat = make_master(sample, method='robust_pca', + downsample=Config.FLAT_FROM_LIGHTS_DOWNSAMPLE) if synthetic_flat is not None: masters['flat'] = synthetic_flat safe_print(f" ✓ Master flat: {len(sample)} light frames (synthetic) -> " - f"{synthetic_flat.shape[0]}×{synthetic_flat.shape[1]} (robust_pca)") + f"{synthetic_flat.shape[0]}×{synthetic_flat.shape[1]} (robust_pca, " + f"{Config.FLAT_FROM_LIGHTS_DOWNSAMPLE}x downsampled decomposition)") else: safe_print(" Synthetic flat-from-lights failed (too few usable frames) " "-- proceeding without a flat") @@ -525,9 +533,6 @@ def _method_tag(method: str) -> str: except Exception as e: safe_print(f" WARNING: Could not read dark EXPTIME ({e}) — dark scaling disabled") - if cal_needed and cal_start is not None and stats is not None: - stats.calibration_time = time.time() - cal_start - # Photon-transfer gain / read-noise from raw bias+flat pairs, for the # Poisson term in --photometry. Only when photometry is requested and no # gain was given explicitly. @@ -567,6 +572,9 @@ def _method_tag(method: str) -> str: safe_print(f" ✓ Vignette map: {os.path.basename(vignette_path)} " f"({vmap.shape[1]}×{vmap.shape[0]})") + if stats is not None: + stats.calibration_time = time.time() - _build_start + return masters diff --git a/src/io_fits.py b/src/io_fits.py index 60a04d2..a6075a3 100644 --- a/src/io_fits.py +++ b/src/io_fits.py @@ -110,9 +110,13 @@ def load_fits(path: str) -> Tuple[np.ndarray, dict]: return data, hdr -def make_master(frames: List[FrameInfo], method: str = 'median') -> Optional[np.ndarray]: +def make_master(frames: List[FrameInfo], method: str = 'median', + downsample: int = 1) -> Optional[np.ndarray]: """Create master calibration frame using streaming (mean), memmap (median), - or robust PCA (low-rank + sparse decomposition, ``method='robust_pca'``).""" + or robust PCA (low-rank + sparse decomposition, ``method='robust_pca'``). + + ``downsample`` only affects the robust_pca path (see + ``robust_pca_master``'s docstring) -- median/mean ignore it.""" if not frames: return None # Probe first frame for shape @@ -124,7 +128,7 @@ def make_master(frames: List[FrameInfo], method: str = 'median') -> Optional[np. if method == 'robust_pca': from src.robust_pca import robust_pca_master - master = robust_pca_master(frames, shape) + master = robust_pca_master(frames, shape, downsample=downsample) if master is not None: return master method = 'median' # too few frames for RPCA -- fall back diff --git a/src/models.py b/src/models.py index c11ac3f..0572196 100644 --- a/src/models.py +++ b/src/models.py @@ -118,12 +118,23 @@ class Config: # underdetermined -- make_master falls back to median ROBUST_PCA_AUTO_MAX_FRAMES = 10 # --auto only auto-upgrades median->robust_pca for a # calibration type at or below this frame count. - # Cost scales O(N^2 x pixels); measured 1264s/~21min - # at N=20 -- too slow for a silent --auto default at - # that scale (see _build_masters's comment). At - # N<=10 that scales to ~(10/20)^2 * 1264s ~= 316s - # (~5.3min), judged tolerable for --auto; still - # opt-in via --master-method robust_pca above this + # Measured 1264s/~21min at N=20 -- too slow for a + # silent --auto default at that scale (see + # _build_masters's comment). Cost vs. N was assumed + # O(N^2 x pixels) (quadratic) when this threshold was + # first set, giving an estimated ~316s/~5.3min at + # N<=10 via (10/20)^2 * 1264s; tools/bench_robust_pca_ + # scale.py later measured the real exponent as N^1.42 + # (sub-quadratic -- iteration count evidently doesn't + # scale down as fast as the per-iteration Gram-matrix + # cost does), giving ~404s/~6.7min at N=10 on that + # run (cross-checked against this file's own N=20 + # anchor at 0.85x -- same ballpark, not exact, since + # it's a different machine/run than the anchor). + # Still judged tolerable for --auto at N<=10; kept at + # 10 rather than widened after seeing the real N=15 + # cost (~717s/~12min) -- opt-in via --master-method + # robust_pca above this frame count regardless ROBUST_PCA_MAX_ITERS = 50 # IALM iterations (each is one economy SVD of an # (N, H*W*C) matrix, via src.robust_pca's # Gram-matrix-trick + native gram_matrix_wide/ @@ -135,6 +146,36 @@ class Config: # 2910s/~48.5min pre-optimization. Bounded, not # adaptive-early-exit beyond the tolerance below) ROBUST_PCA_TOL = 1e-7 # Relative Frobenius-norm residual convergence tolerance + FLAT_FROM_LIGHTS_DOWNSAMPLE = 4 # --flat-from-lights only (never real bias/dark/flat + # robust_pca): block-averages each of the 4 Bayer + # sub-planes by this factor before decomposition, + # cutting P (and the O(N^2 x P) cost) by ~16x -- the + # low-rank content this path actually wants (flat-field + # vignetting, dust motes) is smooth well above a 4-pixel + # scale, unlike a real dark/bias master's per-pixel hot + # pixels, which is why this knob doesn't exist for those. + # Not yet measured end-to-end against a real vignetted + # session (only against a flat synthetic frame with no + # vignetting to recover) -- see robust_pca_master's + # bayer_block_downsample/_upsample for the mechanism. + + # Cosmic-ray rejection auto-skip (see src/debayer.py Sharpness/noise-vs-Siril doc + # in CLAUDE.md for the underlying measurement) + LACOSMIC_LONG_SUB_EXPTIME_S = 25.0 # Median light EXPTIME (s) at/above which + # pipeline.py's >=20-frame auto-skip of + # per-frame L.A.Cosmic is itself skipped (i.e. + # lacosmic stays on) instead of deferring + # entirely to stack-level sigma-clip. Forcing + # lacosmic on measurably cuts stacked noise + # (11-16% across three real sessions) but softens + # stars, and that softening cost tracks sub + # length: ~0% at 30s subs, +6.8% at 20s, +11.8% + # at 10s. 25s sits between the 20s (still a real + # cost) and 30s (negligible) measurements -- + # picked, not measured at that exact value; only + # three real sessions inform this, so treat as a + # starting heuristic pending more data, not a + # calibrated cutoff. # PSF-kernel drizzle resampling DRIZZLE_PSF_KERNEL_SIZE = 9 # Tap radius for the drizzle resample kernel when diff --git a/src/pipeline.py b/src/pipeline.py index ad6dc1f..30e5cd2 100644 --- a/src/pipeline.py +++ b/src/pipeline.py @@ -29,7 +29,7 @@ from src.debayer import autodetect_bayer_orientation, debayer from src.frame_processor import execute_frame_processing, quality_gate, reload_accepted_frames from src.io_fits import load_fits, load_frame, populate_fits_header, save_preview_rgb -from src.models import FrameInfo, ProcessingStats +from src.models import Config, FrameInfo, ProcessingStats from src.plate_solve import solve_plate from src.postprocess import postprocess_stack from src.registration import run_registration_phase, select_reference_frame @@ -511,16 +511,36 @@ def stack_target(frames: List[FrameInfo], output_path: str, args: argparse.Names # the single most expensive Phase-1 step. Drizzle has no per-pixel # rejection, so it keeps lacosmic. Resolved before the resume # branch so checkpoint reloads see the same setting. + # + # That skip trades away a real noise reduction (11-16% measured on + # three real sessions when lacosmic is forced on, from cosmic-ray/ + # hot-pixel spikes the Malvar debayer smears across more pixels in + # R/B than G -- see CLAUDE.md's Sharpness/noise-vs-Siril entry). + # The cost of taking that noise win -- star softening -- tracked + # sub length on those same sessions (~0% at 30s subs, +6.8% at + # 20s, +11.8% at 10s), so on long subs the skip is pure loss: keep + # lacosmic on there even at high frame count. if getattr(args, 'cosmic_ray_rejection', None) is None: _method = getattr(args, 'stack_method', 'auto') _drizzling = float(getattr(args, 'drizzle_scale', 1.0) or 1.0) > 1.0 - if n >= 20 and _method != 'mean' and not _drizzling: + _exptimes = [float(f.header.get('EXPTIME', 0) or 0) for f in lights] + _exptimes = [e for e in _exptimes if e > 0] + _median_exptime = float(np.median(_exptimes)) if _exptimes else None + _long_subs = (_median_exptime is not None + and _median_exptime >= Config.LACOSMIC_LONG_SUB_EXPTIME_S) + if n >= 20 and _method != 'mean' and not _drizzling and not _long_subs: args.cosmic_ray_rejection = False safe_print(f" NOTE: cosmic-ray rejection skipped: {n} frames with " f"rejection stacking removes cosmic rays per-pixel " f"(force with --cosmic-ray-rejection)") else: args.cosmic_ray_rejection = True + if n >= 20 and _long_subs: + safe_print(f" NOTE: cosmic-ray rejection kept on despite {n} " + f"frames: median sub length {_median_exptime:.0f}s " + f">= {Config.LACOSMIC_LONG_SUB_EXPTIME_S:.0f}s, where " + f"it measurably cuts stacked noise at ~no star-softening " + f"cost (disable with --no-cosmic-ray-rejection)") # ====================================================================== # PHASE 1: Process & Analyse @@ -1448,12 +1468,17 @@ def _psf_fallback_suffix() -> str: print(f" Sky (Bortle est): ~{int(round(np.median(bortles)))} " f"(rough, same-equipment-only estimate from background level)") # Per-phase timing with % of wall-clock so the bottleneck is obvious. The - # four phases rarely sum to the total — "Other" captures frame discovery, - # master-calibration building, plate-solve/WCS, colour calibration, and the - # file writes between phases. A large "Other" means the limiter is outside - # the four phases (usually I/O: master building or the output/memmap writes). + # five phases (Calibration + the four stack phases) rarely sum to the + # total — "Other" is what's left over: frame discovery, plate-solve/WCS, + # colour calibration, and the file writes between phases. Calibration + # (master bias/dark/flat -- including --flat-from-lights's synthetic + # flat, the slowest thing that function can do) used to be folded into + # "Other" too, sight unseen: a real --auto run with no flat frames spent + # 2 of its 3 minutes inside _build_masters and "Other (I/O)" was the only + # place that showed, unlabeled, as 69% of the run with no clue why. _tt = max(stats.total_time(), 1e-9) _phases = [ + ("Calibration", stats.calibration_time), ("Quality+Load", stats.quality_time), ("Registration", stats.registration_time), ("Stacking", stats.stacking_time), diff --git a/src/postprocess.py b/src/postprocess.py index 5d565f3..22456b8 100644 --- a/src/postprocess.py +++ b/src/postprocess.py @@ -51,6 +51,43 @@ except Exception: sigma_clipped_stats = None +try: + import astro_native as _native + _HAS_NATIVE_MEDIAN = True +except Exception: + _native = None + _HAS_NATIVE_MEDIAN = False + + +def _median_filter_per_channel(stacked: np.ndarray, size: int) -> np.ndarray: + """``scipy.ndimage.median_filter(stacked, size=(size, size, 1))`` -- + independent per-channel 2D median, not a 3D one -- routed through the + native ``median_filter_native`` kernel per channel when available. + + Found by profiling a real stack: the direct combined-axis + ``ndimage.median_filter(stacked, size=(5,5,1))`` call (post-process's + hot-pixel-removal step, the very first thing Phase 4 does) was 5.1s on a + 2033x3041x3 frame -- scipy's generic N-D rank-filter machinery has no + fast path for a size-1 axis, so it does real (if wasted) work along it. + 3 independent native 2D calls are both algorithmically correct for this + exact shape and match this project's existing per-channel native median + filter (`_median_filter3` in `src/debayer.py`; validated at size 5 + against scipy in `tests/test_native.py::test_median_filter_native_matches_scipy`). + """ + out = np.empty_like(stacked) + for c in range(stacked.shape[2]): + ch = stacked[:, :, c] + med = None + if _HAS_NATIVE_MEDIAN and ch.dtype == np.float32: + try: + med = np.asarray(_native.median_filter_native(np.ascontiguousarray(ch), size)) + except Exception: + med = None + if med is None: + med = ndimage.median_filter(ch, size=size) + out[:, :, c] = med + return out + def _diag_save(img: np.ndarray, diag_dir: Optional[str], counter: list, slug: str) -> None: """Save a float32 FITS snapshot to diag_dir if diagnostic mode is active. @@ -228,7 +265,7 @@ def postprocess_stack( print("\n Removing residual hot pixels (per-channel)...") _hp_start = time.time() _hp_fixed = 0 - _hp_meds = ndimage.median_filter(stacked, size=(5, 5, 1)) + _hp_meds = _median_filter_per_channel(stacked, 5) _hp_diffs = stacked - _hp_meds _hp_mads = np.median(np.abs(_hp_diffs), axis=(0, 1)) _hp_sigmas = np.maximum(_hp_mads * 1.4826, 1e-6) diff --git a/src/robust_pca.py b/src/robust_pca.py index 50378c1..bd9ab7d 100644 --- a/src/robust_pca.py +++ b/src/robust_pca.py @@ -14,6 +14,7 @@ from typing import List, Optional, Tuple import numpy as np +from scipy import ndimage from src.models import Config, FrameInfo from src.utils import get_logger, safe_print @@ -121,31 +122,132 @@ def robust_pca_decompose(D: np.ndarray, max_iters: int = Config.ROBUST_PCA_MAX_I mu_bar = mu * 1e7 rho = 1.5 + # S/Y allocated once and reused every iteration (native writes them in + # place) instead of the ~dozen fresh (N,P) temporaries/iteration the + # plain-numpy expressions below build -- found by profiling a real + # --flat-from-lights run (N=10, P=6.25M): robust_pca_decompose was + # spending 100s of a 126s total in exactly this elementwise arithmetic + # around the SVD, not the SVD itself. L = np.zeros_like(D) S = np.zeros_like(D) Y = np.zeros_like(D) + D_c = np.ascontiguousarray(D) for _ in range(max_iters): # L-update: singular value shrinkage (proximal operator of the nuclear norm) - U, sigma, Vt = _thin_svd_wide(D - S + Y / mu) + svd_input = None + if _HAS_NATIVE: + try: + svd_input = np.asarray(_native.robust_pca_pre_svd_input(D_c, S, Y, mu)) + except Exception: + svd_input = None + if svd_input is None: + svd_input = D - S + Y / mu + U, sigma, Vt = _thin_svd_wide(svd_input) sigma_shrunk = np.maximum(sigma - 1.0 / mu, 0.0) - L = (U * sigma_shrunk) @ Vt + # (U * sigma_shrunk) @ Vt is the exact same small-(N,N)-times-wide- + # (N,P) shape _thin_svd_wide's own Gram-trick GEMM is (small_times_ + # wide above) -- routing it through the same native kernel instead of + # numpy's plain `@` matters because this machine's numpy has no + # optimized BLAS (numpy.show_config reports `blas: name: auto`, + # measured ~1.4 GFLOPS, reference-BLAS-speed): a lone `np.matmul` + # call here was 32.9s of this function's 62.5s once the elementwise + # arithmetic around it was already fused (profiled after the fusion + # above, not guessed). + weighted_U = np.ascontiguousarray(U * sigma_shrunk) + L = None + if _HAS_NATIVE: + try: + L = np.asarray(_native.small_times_wide(weighted_U, np.ascontiguousarray(Vt))) + except Exception: + L = None + if L is None: + L = weighted_U @ Vt + + # S-update + residual + Y-update + residual norm (proximal operator + # of the L1 norm, fused): native writes S and Y in place. + lam_over_mu = lam / mu + resid_norm = None + if _HAS_NATIVE: + try: + resid_norm = float(_native.robust_pca_iterate(D_c, L, S, Y, lam_over_mu, mu)) + except Exception: + resid_norm = None + if resid_norm is None: + temp = D - L + Y / mu + S[...] = np.sign(temp) * np.maximum(np.abs(temp) - lam_over_mu, 0.0) + residual = D - L - S + Y += mu * residual + resid_norm = float(np.linalg.norm(residual, 'fro')) - # S-update: elementwise soft threshold (proximal operator of the L1 norm) - temp = D - L + Y / mu - S = np.sign(temp) * np.maximum(np.abs(temp) - lam / mu, 0.0) - - residual = D - L - S - Y = Y + mu * residual mu = min(mu * rho, mu_bar) - if float(np.linalg.norm(residual, 'fro')) / norm_fro < tol: + if resid_norm / norm_fro < tol: break return L, S -def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...]) -> Optional[np.ndarray]: +def _bayer_planes(img: np.ndarray) -> List[np.ndarray]: + """The 4 same-colour sub-planes of a 2x2-CFA mosaic, R/G1/G2/B order.""" + return [img[0::2, 0::2], img[0::2, 1::2], img[1::2, 0::2], img[1::2, 1::2]] + + +def bayer_block_downsample(img: np.ndarray, k: int) -> np.ndarray: + """Downsample a raw Bayer mosaic by block-averaging *within* each of its 4 + same-colour sub-planes independently, by a factor of ``k`` each, then + re-interleaving into a smaller mosaic with the same 2x2 CFA pattern. + + Averaging across the raw mosaic directly (the naive block-average) would + mix different sensor colour channels together -- R/G/B have very + different absolute response (real gains, not just vignetting), so that + would corrupt the result, not just blur it. Averaging within one colour + at a time is what a flat's own low-rank content (vignetting, dust motes) + survives fine: both vary smoothly over tens-to-hundreds of pixels, well + above the few-pixel scale a modest ``k`` removes. Any remainder pixels + (sub-plane dimension not divisible by ``k``) are cropped, not padded -- + at most ``k-1`` pixels per edge, negligible against a full sensor plane. + """ + if k <= 1: + return img + planes_small = [] + for plane in _bayer_planes(img): + h, w = plane.shape + h2, w2 = h - h % k, w - w % k + p = plane[:h2, :w2].reshape(h2 // k, k, w2 // k, k) + planes_small.append(p.mean(axis=(1, 3))) + sh, sw = planes_small[0].shape + small = np.empty((sh * 2, sw * 2), dtype=planes_small[0].dtype) + small[0::2, 0::2] = planes_small[0] + small[0::2, 1::2] = planes_small[1] + small[1::2, 0::2] = planes_small[2] + small[1::2, 1::2] = planes_small[3] + return small + + +def bayer_block_upsample(img_small: np.ndarray, target_shape: Tuple[int, int]) -> np.ndarray: + """Inverse of ``bayer_block_downsample``: upsample each of the 4 CFA + sub-planes independently (bilinear -- a flat is smooth by assumption, so + this is the same kind of interpolation the flat is later applied + through) and re-interleave, cropped/edge-padded to ``target_shape`` + exactly (the downsample's remainder crop means the scale factor isn't + always a perfectly round number).""" + H, W = target_shape + sub_h, sub_w = H // 2, W // 2 + out = np.empty((H, W), dtype=img_small.dtype) + planes_small = _bayer_planes(img_small) + for (r_off, c_off), p in zip([(0, 0), (0, 1), (1, 0), (1, 1)], planes_small): + zoom_y = sub_h / p.shape[0] + zoom_x = sub_w / p.shape[1] + up = ndimage.zoom(p, (zoom_y, zoom_x), order=1)[:sub_h, :sub_w] + if up.shape != (sub_h, sub_w): + up = np.pad(up, ((0, sub_h - up.shape[0]), (0, sub_w - up.shape[1])), mode='edge') + out[r_off::2, c_off::2] = up + return out + + +def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...], + downsample: int = 1) -> Optional[np.ndarray]: """Build a master calibration frame via robust PCA. Returns the per-pixel median of the recovered low-rank component ``L`` @@ -154,9 +256,27 @@ def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...]) -> Option loaded successfully, the stack won't fit in available memory, or any loaded frame is non-finite -- caller should fall back to a plain median in every ``None`` case. + + ``downsample`` (> 1): block-averages each loaded frame via + ``bayer_block_downsample`` before decomposition and upsamples the + recovered master back to ``shape`` via ``bayer_block_upsample``. + Deliberately **not** used for real bias/dark/flat robust-PCA masters -- + only ``--flat-from-lights`` passes this, where the sparse component only + needs to isolate stars/nebula structure that dithers between frames (not + per-pixel hot pixels the way a real dark/bias decomposition does), and + the low-rank component it's actually after (vignetting, dust motes) is + smooth well above the pixel scale a modest downsample removes. Cuts P + (and so the O(N^2 x P) decomposition cost) by ``downsample**2``. """ from src.io_fits import load_frame + # A CFA-respecting downsample only makes sense on a real 2x2-mosaic raw + # frame; anything else (mono, already-debayered, odd dimensions) falls + # back to full resolution rather than guessing at a layout. + can_downsample = (downsample > 1 and len(shape) == 2 + and shape[0] % 2 == 0 and shape[1] % 2 == 0) + work_shape = shape + imgs = [] skipped = 0 for f in frames: @@ -171,7 +291,11 @@ def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...]) -> Option # mismatched frame can genuinely reach here. skipped += 1 continue - imgs.append(data.astype(np.float64)) + data = data.astype(np.float64) + if can_downsample: + data = bayer_block_downsample(data, downsample) + work_shape = data.shape + imgs.append(data) if skipped: _log.warning( @@ -188,7 +312,7 @@ def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...]) -> Option return None n = len(imgs) - p = int(np.prod(shape)) + p = int(np.prod(work_shape)) # IALM keeps ~7 live (n, p) float64 arrays alive at once (D, L, S, Y, temp, # residual, plus the SVD's working copies) -- guard against OOM the same # way make_master's median path guards its memmap threshold, since --auto @@ -220,5 +344,7 @@ def robust_pca_master(frames: List[FrameInfo], shape: Tuple[int, ...]) -> Option return None L, _S = robust_pca_decompose(D) - master = np.median(L, axis=0).reshape(shape) + master = np.median(L, axis=0).reshape(work_shape) + if work_shape != shape: + master = bayer_block_upsample(master, shape) return master.astype(np.float32) diff --git a/tests/test_dbe_gradient.py b/tests/test_dbe_gradient.py index 8fb5882..9205e4a 100644 --- a/tests/test_dbe_gradient.py +++ b/tests/test_dbe_gradient.py @@ -59,3 +59,54 @@ def test_a_bright_object_is_not_admitted_as_gradient(): out = bg.dynamic_background_extraction(img, patch_size=32) # the bump is nowhere near the edges; the DBE surface must not have absorbed it assert float(out[210, 330].mean()) - float(np.median(out[100:140, 250:400].mean(axis=2))) > 300.0 + + +def test_emission_mask_disk_dilation_matches_scipy_binary_dilation(): + """`_build_emission_mask`'s compact-source dilation switched from + `scipy.ndimage.binary_dilation` with a disk structuring element to a + distance-transform threshold (same result, much faster on a large + contiguous bright region -- see background.py's comment at the call + site). Pin the two as exactly equivalent directly, not just through the + higher-level DBE tests above, since a single wrong offset-by-one in the + distance-transform version would still likely pass those.""" + from scipy.ndimage import binary_dilation, distance_transform_edt + + H, W = 300, 400 + yy, xx = np.indices((H, W)) + # a contiguous elliptical blob, not scattered noise -- this is exactly + # the shape that made scipy's generic binary_dilation pathologically + # slow (a real bright galaxy/comet core), and is what actually exercises + # the interior of the dilated region, not just its edge. + src_binary = (((yy - 150) / 40.0) ** 2 + ((xx - 200) / 60.0) ** 2 <= 1.0).astype(np.uint8) + r = 25 + y_idx, x_idx = np.ogrid[-r:r + 1, -r:r + 1] + structure = (y_idx ** 2 + x_idx ** 2 <= r ** 2).astype(np.uint8) + + want = binary_dilation(src_binary, structure=structure) + got = distance_transform_edt(1 - src_binary) <= r + + np.testing.assert_array_equal(got, want) + + +def test_gaussian_filter_ds_bilinear_upsample_matches_full_resolution(): + """gaussian_filter_ds's large-sigma path downsamples, blurs, then + upsamples back -- its own upsample step switched from cubic (order=3) to + bilinear (order=1) interpolation, ~3x faster, on the reasoning that a + coarse grid which just came out of a huge Gaussian blur has nothing left + for cubic's extra curvature term to recover. Check that claim directly: + the ds-optimized result should stay close to a real full-resolution + gaussian_filter at the same sigma, which is the property the whole + function exists to preserve regardless of interpolation order.""" + from scipy.ndimage import gaussian_filter + + rng = np.random.default_rng(6) + H, W = 300, 400 + field = rng.normal(1000.0, 20.0, (H, W)) + sigma = 40.0 # > ds_threshold=24.0, exercises the downsampled path + + exact = gaussian_filter(field, sigma=sigma) + fast = bg.gaussian_filter_ds(field, sigma) + + 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 diff --git a/tests/test_native.py b/tests/test_native.py index 3cf9f01..5cf2b9c 100644 --- a/tests/test_native.py +++ b/tests/test_native.py @@ -12,6 +12,7 @@ import src.debayer as _debayer_mod import src.denoising as _denoising_mod import src.local_normalize as _local_normalize_mod +import src.postprocess as _postprocess_mod import src.robust_pca as _robust_pca_mod import src.stacking as _stacking_mod import src.star_removal as _star_removal_mod @@ -519,6 +520,30 @@ def test_median_filter_native_matches_scipy(size): assert float(np.max(np.abs(ref.astype(np.float64) - got.astype(np.float64)))) < 1e-4 +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 + real full-res stack -- scipy's N-D rank filter has no fast path for a + size-1 axis) to 3 independent native 2D calls, one per channel. Confirm + that's actually equivalent, not just faster.""" + from scipy import ndimage + rng = np.random.default_rng(4) + stacked = rng.normal(500, 50, (60, 70, 3)).astype(np.float32) + ref = ndimage.median_filter(stacked, size=(5, 5, 1)) + got = _postprocess_mod._median_filter_per_channel(stacked, 5) + assert float(np.max(np.abs(ref.astype(np.float64) - got.astype(np.float64)))) < 1e-4 + + +def test_median_filter_per_channel_falls_back_without_native(monkeypatch): + from scipy import ndimage + rng = np.random.default_rng(5) + stacked = rng.normal(500, 50, (40, 50, 3)).astype(np.float32) + ref = ndimage.median_filter(stacked, size=(3, 3, 1)) + monkeypatch.setattr(_postprocess_mod, '_HAS_NATIVE_MEDIAN', False) + got = _postprocess_mod._median_filter_per_channel(stacked, 3) + assert float(np.max(np.abs(ref.astype(np.float64) - got.astype(np.float64)))) < 1e-4 + + def test_all_nan_pixel_is_zero(): d = _stack(n=8, h=4, w=4, c=1, outliers=False) d[:, 0, 0, 0] = np.nan @@ -993,6 +1018,96 @@ def test_small_times_wide_rejects_shape_mismatch(): native.small_times_wide(small, data) +# --------------------------------------------------------------------------- +# robust_pca_pre_svd_input / robust_pca_iterate: the fused per-iteration IALM +# elementwise kernels (see robust_pca.py's robust_pca_decompose). Found by +# profiling a real --flat-from-lights run where the plain-numpy elementwise +# arithmetic surrounding the SVD -- not the SVD itself -- was ~80% of the +# function's wall time. +# --------------------------------------------------------------------------- + +def test_robust_pca_pre_svd_input_matches_numpy(): + rng = np.random.default_rng(10) + n, p = 9, 500 + d = np.ascontiguousarray(rng.normal(0.0, 5.0, (n, p))) + s = np.ascontiguousarray(rng.normal(0.0, 1.0, (n, p))) + y = np.ascontiguousarray(rng.normal(0.0, 1.0, (n, p))) + mu = 0.37 + got = np.asarray(native.robust_pca_pre_svd_input(d, s, y, mu)) + want = d - s + y / mu + np.testing.assert_array_equal(got, want) # same f64 op order -- bit-exact + + +def test_robust_pca_pre_svd_input_rejects_shape_mismatch(): + rng = np.random.default_rng(11) + d = np.ascontiguousarray(rng.normal(0.0, 1.0, (5, 100))) + s = np.ascontiguousarray(rng.normal(0.0, 1.0, (5, 100))) + y = np.ascontiguousarray(rng.normal(0.0, 1.0, (4, 100))) # N mismatch + with pytest.raises(ValueError): + native.robust_pca_pre_svd_input(d, s, y, 0.5) + + +def test_robust_pca_iterate_matches_numpy(): + rng = np.random.default_rng(12) + n, p = 9, 500 + d = np.ascontiguousarray(rng.normal(0.0, 5.0, (n, p))) + l = np.ascontiguousarray(rng.normal(0.0, 4.0, (n, p))) + y0 = np.ascontiguousarray(rng.normal(0.0, 1.0, (n, p))) + mu = 0.42 + lam_over_mu = 0.1 + + s_native = np.zeros((n, p)) + y_native = y0.copy() + resid_norm = float(native.robust_pca_iterate(d, l, s_native, y_native, lam_over_mu, mu)) + + temp = d - l + y0 / mu + s_want = np.sign(temp) * np.maximum(np.abs(temp) - lam_over_mu, 0.0) + residual_want = d - l - s_want + y_want = y0 + mu * residual_want + resid_norm_want = float(np.linalg.norm(residual_want, 'fro')) + + np.testing.assert_array_equal(s_native, s_want) # per-element op, bit-exact + np.testing.assert_array_equal(y_native, y_want) + # the norm is a parallel reduction, so only close, not bit-exact -- see + # the kernel's own docstring + np.testing.assert_allclose(resid_norm, resid_norm_want, rtol=1e-10) + + +def test_robust_pca_iterate_rejects_shape_mismatch(): + rng = np.random.default_rng(13) + d = np.ascontiguousarray(rng.normal(0.0, 1.0, (5, 100))) + l = np.ascontiguousarray(rng.normal(0.0, 1.0, (5, 100))) + s = np.zeros((5, 100)) + y = np.zeros((4, 100)) # N mismatch + with pytest.raises(ValueError): + native.robust_pca_iterate(d, l, s, y, 0.1, 0.5) + + +def test_robust_pca_decompose_native_matches_numpy_fallback(): + """End-to-end: robust_pca_decompose's native and numpy-fallback paths + must converge to the same low-rank/sparse split, not just agree on the + individual fused kernels in isolation.""" + rng = np.random.default_rng(14) + n, p = 10, 800 + low_rank = np.outer(rng.uniform(0.8, 1.2, n), rng.normal(1000.0, 50.0, p)) + sparse = np.zeros((n, p)) + idx = rng.integers(0, n * p, n * p // 50) + sparse.flat[idx] = rng.uniform(200.0, 2000.0, len(idx)) + D = low_rank + sparse + rng.normal(0.0, 5.0, (n, p)) + + assert _robust_pca_mod._HAS_NATIVE + L_native, S_native = _robust_pca_mod.robust_pca_decompose(D) + + _robust_pca_mod._HAS_NATIVE = False + try: + L_numpy, S_numpy = _robust_pca_mod.robust_pca_decompose(D) + finally: + _robust_pca_mod._HAS_NATIVE = True + + np.testing.assert_allclose(L_native, L_numpy, atol=1e-6, rtol=1e-6) + np.testing.assert_allclose(S_native, S_numpy, atol=1e-6, rtol=1e-6) + + # --------------------------------------------------------------------------- # continuum_scale_moments: single-pass central moments backing # optimal_continuum_scale's closed-form skewness-vs-scale polynomial diff --git a/tests/test_robust_pca.py b/tests/test_robust_pca.py index 89c7118..5cdc1b7 100644 --- a/tests/test_robust_pca.py +++ b/tests/test_robust_pca.py @@ -146,6 +146,109 @@ def test_falls_back_when_memory_insufficient(self): self.assertIsNone(robust_pca_master(frames, shape)) +class TestBayerBlockDownsampleUpsample(unittest.TestCase): + """CFA-respecting down/upsample used by robust_pca_master's downsample= + parameter (only --flat-from-lights passes it -- see FLAT_FROM_LIGHTS_ + DOWNSAMPLE's Config docstring for why real dark/bias/flat masters never do).""" + + def test_does_not_mix_bayer_channels(self): + # Four constant Bayer sub-planes at very different levels -- a + # cross-channel-mixing bug (averaging the raw mosaic directly instead + # of each sub-plane independently) would blend these together. + from src.robust_pca import bayer_block_downsample + img = np.zeros((16, 16)) + img[0::2, 0::2] = 100.0 # R + img[0::2, 1::2] = 500.0 # G1 + img[1::2, 0::2] = 500.0 # G2 + img[1::2, 1::2] = 900.0 # B + small = bayer_block_downsample(img, 2) + self.assertTrue(np.all(small[0::2, 0::2] == 100.0)) + self.assertTrue(np.all(small[0::2, 1::2] == 500.0)) + self.assertTrue(np.all(small[1::2, 0::2] == 500.0)) + self.assertTrue(np.all(small[1::2, 1::2] == 900.0)) + + def test_downsample_shrinks_by_the_requested_factor(self): + from src.robust_pca import bayer_block_downsample + img = np.zeros((64, 96)) + small = bayer_block_downsample(img, 4) + # each sub-plane is (32,48) -> downsampled by 4 -> (8,12) -> reinterleaved (16,24) + self.assertEqual(small.shape, (16, 24)) + + def test_roundtrip_preserves_a_smooth_pattern(self): + # The point of the downsample is to survive exactly this kind of + # content (smooth vignetting), not preserve it exactly. + from src.robust_pca import bayer_block_downsample, bayer_block_upsample + H, W = 200, 300 + yy, xx = np.mgrid[0:H, 0:W] + pattern = 1000.0 - 0.01 * ((yy - H / 2) ** 2 + (xx - W / 2) ** 2) + small = bayer_block_downsample(pattern, 4) + back = bayer_block_upsample(small, (H, W)) + self.assertEqual(back.shape, (H, W)) + rel_err = np.abs(back - pattern).mean() / np.abs(pattern).mean() + self.assertLess(rel_err, 0.02) + + def test_downsample_is_a_noop_at_factor_one(self): + from src.robust_pca import bayer_block_downsample + img = np.random.default_rng(0).normal(size=(20, 20)) + np.testing.assert_array_equal(bayer_block_downsample(img, 1), img) + + +class TestRobustPcaMasterDownsample(unittest.TestCase): + """End-to-end: robust_pca_master(downsample=N) still recovers a real + vignetting pattern at full output resolution, not just a smaller one.""" + + def _write_frame(self, tmpdir: str, name: str, data: np.ndarray) -> FrameInfo: + path = os.path.join(tmpdir, name) + _write_fits(path, data) + return FrameInfo(path=path, type='light', header={}) + + def test_recovers_vignette_at_full_resolution(self): + rng = np.random.default_rng(2) + shape = (200, 300) # even dims, real-mosaic-shaped + yy, xx = np.mgrid[0:shape[0], 0:shape[1]] + # A smooth vignetting-like pattern, per-Bayer-position gain baked in + # (R/G1/G2/B different absolute levels) -- a channel-mixing bug in + # the downsample would show up as a wrong recovered pattern here, + # not just a blurrier one. + base = 1000.0 - 0.02 * ((yy - 100) ** 2 + (xx - 150) ** 2) + gain = np.ones(shape) + gain[0::2, 0::2] = 1.0 + gain[0::2, 1::2] = 1.9 + gain[1::2, 0::2] = 1.9 + gain[1::2, 1::2] = 1.6 + pattern = base * gain + + with tempfile.TemporaryDirectory() as d: + frames = [] + for i in range(10): + frame = pattern + rng.normal(0, 3.0, shape) + frames.append(self._write_frame(d, f'f{i}.fits', frame.astype(np.float32))) + + master = robust_pca_master(frames, shape, downsample=Config.FLAT_FROM_LIGHTS_DOWNSAMPLE) + self.assertIsNotNone(master) + self.assertEqual(master.shape, shape) # upsampled back to full res + self.assertTrue(np.all(np.isfinite(master))) + rel_err = np.abs(master - pattern).mean() / np.abs(pattern).mean() + self.assertLess(rel_err, 0.15) # looser than the full-res test -- it's lossy by design + + def test_downsample_one_matches_undownsampled_call(self): + # downsample=1 must be exactly the pre-existing code path (no + # upsample step, no behavior change for real bias/dark/flat masters + # that never pass downsample at all). + rng = np.random.default_rng(3) + shape = (24, 24) + yy, xx = np.mgrid[0:shape[0], 0:shape[1]] + pattern = 1000.0 - 0.3 * ((yy - 12) ** 2 + (xx - 12) ** 2) + with tempfile.TemporaryDirectory() as d: + frames = [] + for i in range(8): + frame = pattern + rng.normal(0, 2.0, shape) + frames.append(self._write_frame(d, f'f{i}.fits', frame.astype(np.float32))) + master_plain = robust_pca_master(frames, shape) + master_ds1 = robust_pca_master(frames, shape, downsample=1) + np.testing.assert_array_equal(master_plain, master_ds1) + + class TestMakeMasterRobustPcaFallback(unittest.TestCase): """make_master(method='robust_pca') dispatch and graceful fallback.""" diff --git a/tools/bench_native.py b/tools/bench_native.py index 3f8a584..2e90751 100644 --- a/tools/bench_native.py +++ b/tools/bench_native.py @@ -105,6 +105,38 @@ def bench(name, fn, warm=1, reps=3): lambda: nat.aperture_photometry_batch( apb_img, apb_xs, apb_ys, 6.0, 9.0, 15.0, 4)) +# --- Gram-matrix thin-SVD trick: robust-PCA calibration-stack shape --- +if hasattr(nat, "gram_matrix_wide"): + rpca_n, rpca_p = 10, 2048 * 3056 # this project's own real profiled case (--flat-from-lights) + rpca_M = np.ascontiguousarray(rng.normal(1000.0, 50.0, (rpca_n, rpca_p))) + results["gram_matrix_wide"] = bench( + f"gram_matrix_wide (N={rpca_n})", lambda: nat.gram_matrix_wide(rpca_M), reps=2) + rpca_small = np.ascontiguousarray(rng.normal(0.0, 1.0, (rpca_n, rpca_n))) + results["small_times_wide"] = bench( + f"small_times_wide (N={rpca_n})", lambda: nat.small_times_wide(rpca_small, rpca_M), reps=2) + +# --- robust_pca_pre_svd_input / robust_pca_iterate: fused IALM per-iteration update --- +if hasattr(nat, "robust_pca_iterate"): + rpca_D = np.ascontiguousarray(rng.normal(1000.0, 50.0, (rpca_n, rpca_p))) + rpca_L = np.ascontiguousarray(rng.normal(1000.0, 50.0, (rpca_n, rpca_p))) + results["robust_pca_pre_svd_input"] = bench( + f"robust_pca_pre_svd_input (N={rpca_n})", + lambda: nat.robust_pca_pre_svd_input( + rpca_D, np.zeros((rpca_n, rpca_p)), np.zeros((rpca_n, rpca_p)), 0.37), + reps=2) + results["robust_pca_iterate"] = bench( + f"robust_pca_iterate (N={rpca_n})", + lambda: nat.robust_pca_iterate( + rpca_D, rpca_L, np.zeros((rpca_n, rpca_p)), np.zeros((rpca_n, rpca_p)), 0.1, 0.42), + reps=2) + +# --- median_filter_native: postprocess.py's per-channel hot-pixel step shape --- +if hasattr(nat, "median_filter_native"): + mf_plane = np.ascontiguousarray(rng.normal(500.0, 50.0, (2033, 3041)).astype(np.float32)) + results["median_filter_5x5"] = bench( + "median_filter_native (5x5, full-res plane)", + lambda: nat.median_filter_native(mf_plane, 5)) + out = sys.argv[1] if len(sys.argv) > 1 else "bench_results.json" with open(out, "w") as f: json.dump(results, f, indent=1) diff --git a/tools/bench_robust_pca_scale.py b/tools/bench_robust_pca_scale.py new file mode 100644 index 0000000..675ca24 --- /dev/null +++ b/tools/bench_robust_pca_scale.py @@ -0,0 +1,123 @@ +"""Benchmark src.robust_pca.robust_pca_decompose's wall-clock cost vs. frame count N, +to check whether Config.ROBUST_PCA_AUTO_MAX_FRAMES (currently 10, based on a single +N=20/P=18M anchor documented in CLAUDE.md: 1264s) can be safely widened for --auto. + +Running the full N=20, P=18,000,000 (2000x3000x3) case directly takes ~21 minutes per +the existing measurement -- too slow to redo here per candidate N. Instead this measures +the real decomposition (astro_native's Gram-matrix-trick SVD, same code path as +production) at reduced P across a range of N, confirms the O(N^2 x P) scaling the +kernel's docstring claims actually holds on this machine, then extrapolates to the real +P=18M shape -- cross-checked against the documented N=20 anchor. + +Usage: python tools/bench_robust_pca_scale.py +""" +import os +import sys +import time + +import numpy as np + +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +from src.robust_pca import robust_pca_decompose + +try: + import astro_native # noqa: F401 + print(f"native: available ({astro_native.__file__})") +except Exception as e: + print(f"native: NOT available ({e}) -- this benchmark would be measuring the " + f"numpy fallback, not the production path. Build astro_native first.") + raise SystemExit(1) + +rng = np.random.default_rng(7) + +REAL_P = 2000 * 3000 * 3 # matches the documented N=20 anchor's frame shape +ANCHOR_N, ANCHOR_SECONDS = 20, 1264.0 # from CLAUDE.md / src/robust_pca.py docstring + + +def make_calib_stack(n: int, p: int) -> np.ndarray: + """Synthetic (n, p) calibration stack: a shared rank-1 pattern (vignetting/dark + current) + sparse spikes (dust motes, hot pixels) + read noise -- same qualitative + structure robust_pca_decompose is built to separate, so convergence behavior (and + therefore iteration count / wall time) is representative of real calibration data, + not a degenerate all-zero or pure-noise case that would converge trivially fast. + """ + shared = rng.normal(1000, 50, p) + weights = rng.uniform(0.8, 1.2, n) + D = weights[:, None] * shared[None, :] + n_spikes = int(n * p * 0.01) + idx_n = rng.integers(0, n, n_spikes) + idx_p = rng.integers(0, p, n_spikes) + D[idx_n, idx_p] += rng.choice([-1, 1], n_spikes) * rng.uniform(200, 2000, n_spikes) + D += rng.normal(0, 15, (n, p)) + return D.astype(np.float64) + + +def timed_decompose(n: int, p: int, reps: int = 1) -> float: + best = None + for _ in range(reps): + D = make_calib_stack(n, p) + t0 = time.perf_counter() + robust_pca_decompose(D) + dt = time.perf_counter() - t0 + best = dt if best is None else min(best, dt) + return best + + +print("\n--- Phase 1: N-scaling at fixed small P (confirm ~N^2 growth) ---") +P_SMALL = 258 * 258 * 3 # ~200k px, keeps each run to a few seconds +n_values = [10, 15, 20, 25, 30, 40] +n_results = {} +for n in n_values: + dt = timed_decompose(n, P_SMALL) + n_results[n] = dt + print(f" N={n:3d} P={P_SMALL:>9,d} {dt:7.3f}s") + +print("\n--- Phase 2: P-scaling at fixed N=10 (confirm linear growth) ---") +p_values = [50_000, 100_000, 200_000, 400_000] +p_results = {} +for p in p_values: + dt = timed_decompose(10, p) + p_results[p] = dt + print(f" N=10 P={p:>9,d} {dt:7.3f}s") + +# Fit growth exponents from the measured points (log-log slope) instead of assuming +# the docstring's O(N^2 x P) claim holds exactly on this machine/data. +import math + +log_n = [math.log(n) for n in n_values] +log_t_n = [math.log(max(t, 1e-6)) for t in n_results.values()] +n_slope = np.polyfit(log_n, log_t_n, 1)[0] + +log_p = [math.log(p) for p in p_values] +log_t_p = [math.log(max(t, 1e-6)) for t in p_results.values()] +p_slope = np.polyfit(log_p, log_t_p, 1)[0] + +print(f"\nMeasured scaling exponents: time ~ N^{n_slope:.2f} x P^{p_slope:.2f} " + f"(docstring claims O(N^2 x P), i.e. exponents 2.0 and 1.0)") + +# Extrapolate each candidate N to the real P using the largest small-P measurement as +# the base point (closest to the regime we care about) and the measured P exponent. +base_n, base_p = 40, P_SMALL +base_t = n_results[base_n] +p_scale = (REAL_P / base_p) ** p_slope + +print(f"\n--- Phase 3: extrapolated full-resolution (P={REAL_P:,d}) time by N ---") +print(f" {'N':>4s} {'predicted':>10s} {'vs anchor':>10s}") +candidate_ns = [10, 15, 20, 25, 30, 40] +for n in candidate_ns: + n_scale = (n / base_n) ** n_slope + predicted = base_t * n_scale * p_scale + tag = "" + if n == ANCHOR_N: + tag = f" (documented anchor: {ANCHOR_SECONDS:.0f}s / {ANCHOR_SECONDS/60:.1f}min)" + print(f" {n:4d} {predicted:8.1f}s {predicted/60:8.2f}min{tag}") + +anchor_scale = (ANCHOR_N / base_n) ** n_slope * p_scale +anchor_predicted = base_t * anchor_scale +print(f"\nCross-check: predicted N=20 full-res time {anchor_predicted:.0f}s vs. " + f"documented anchor {ANCHOR_SECONDS:.0f}s " + f"(ratio {anchor_predicted / ANCHOR_SECONDS:.2f}x)") +print("A ratio far from 1.0 means this machine/build differs enough from the original " + "measurement that the extrapolation above should not be trusted for a threshold " + "change -- rerun the real N=20 P=18M case directly instead.") From b842bc2aaa6effce2fc54a8edc7b688475a54b86 Mon Sep 17 00:00:00 2001 From: Hans Davenport <35202271+hd152@users.noreply.github.com> Date: Tue, 22 Sep 2026 10:52:26 -0700 Subject: [PATCH 2/4] Refresh Omega benchmark numbers on README/website after the perf pass Ran tools/bench_vs_siril.py for real against C:\astro\Omega_Nebula (Siril is installed locally) on this branch's code. Found the earlier sub-length-aware cosmic-ray fix now correctly triggers on this session's real 30s subs, which the previous commit's Phase 4 speedups don't touch: total time 113s -> 132s (per-frame cosmic-ray detection now runs, at its documented 39-75% cost), noise 0.99-1.19x Siril's -> 0.87-1.03x (its documented 11-16% quieter benefit). Net: slower wall clock, quieter and equally sharp output, on this specific showcase session. Updated README.md and docs/index.html with the real new Omega numbers throughout (table, star-width chart, Speed/Noise/Memory prose). Sunflower and Sculptor's subs (20s, 10s) fall under the threshold and keep the old default behavior, but their Phase 4 time should still be faster from the same commit's fixes -- not re-measured (no local data for those sessions), so both pages say so explicitly rather than leaving them silently stale. Co-Authored-By: Claude Sonnet 5 --- CHANGELOG.md | 19 +++++++++++++++++++ README.md | 12 +++++++----- docs/index.html | 38 +++++++++++++++++++------------------- 3 files changed, 45 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b8820cf..9906b6c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,25 @@ match the `VERSION` file and `v*` git tags. difference, and reverted rather than shipped. Full detail, including what was tried and didn't help, in [CLAUDE.md](CLAUDE.md). +- **A separate, unrelated fix from the same session moved the Omega benchmark numbers too.** The first change + of the session made per-frame cosmic-ray rejection's >=20-frame auto-skip sub-length-aware (see below) -- + Omega's real 30s subs cross that new threshold, so a `tools/bench_vs_siril.py` re-run on the real session + found the README/website comparison table stale: total time 113s -> 132s (per-frame detection now runs, + same documented 39-75% cost), noise 0.99-1.19x Siril's -> 0.87-1.03x (same documented 11-16% quieter + benefit). Sunflower and Sculptor's subs (20s, 10s) are under the threshold and unaffected by it, but their + Phase 4 time should improve too from the fixes above -- not yet re-measured, flagged as such on both pages + rather than left silently stale. README.md and docs/index.html updated with the real Omega numbers. + +- **Per-frame cosmic-ray rejection's auto-skip is now sub-length-aware, not just frame-count-aware.** + Previously any >=20-frame session with a rejection-based stack method skipped per-frame L.A.Cosmic + entirely, on the theory that stack-level sigma-clip catches cosmic rays just as well. True on average, but + the noise-vs-softening tradeoff of forcing it back on (documented in the 2.2.4 entry below) tracks sub + length: ~0% star softening at 30s subs, +6.8% at 20s, +11.8% at 10s, while the noise win (11-16%) holds + across all three. `Config.LACOSMIC_LONG_SUB_EXPTIME_S` (25s) now keeps per-frame detection on even at + >=20 frames when the session's median sub exposure is long enough that the softening cost is negligible, + instead of always deferring to the frame-count skip. Only three real sessions inform the exact threshold, + so it's a starting heuristic, not a calibrated cutoff. + ### Added - **Self-update check.** The CLI prints one line at the end of a run, and the desktop app shows a small diff --git a/README.md b/README.md index bafb768..62c1f76 100644 --- a/README.md +++ b/README.md @@ -815,16 +815,18 @@ OriginStack first, Siril second. "Stack only" leaves out OriginStack's finishing | Session | Stack only | OriginStack, finished image | Peak memory | Star width | |---------|-----------|-----------------------------|-------------|------------| -| Omega, 114 frames | 57 s / 72 s | 110 s | 16.5 / 8.1 GB | 4.40 / 4.56 px | +| Omega, 114 frames | 102 s / 63 s | 132 s | 18.4 / 8.1 GB | 4.39 / 4.56 px | | Sunflower, 158 frames | 106 s / 81 s | 159 s | 16.9 / 8.7 GB | 3.56 / 3.60 px | | Sculptor, 532 frames | 258 s / 150 s | 270 s | 16.7 / 11.2 GB | 2.78 / 2.95 px | - **Sharpness: OriginStack is as sharp or sharper.** Its stars are 3% narrower than Siril's on Omega, level on Sunflower and 6% narrower on Sculptor. It used every frame (114/114, 158/158 and 532/532); Siril used 111, 157 and 525. -- **Speed depends on the session.** OriginStack stacks Omega faster and the two larger sessions slower (1.3× as long on Sunflower, 1.7× on Sculptor). Timings move by up to 15% between runs of the same code: Siril's Omega ranged from 53 to 72 s across our runs, OriginStack's Sunflower from 83 to 106 s. -- **Noise: Siril is still cleaner.** OriginStack's per-pixel noise (after matching the flux scale per channel) is 0.99–1.19× Siril's on Omega, 1.09–1.20× on Sunflower and 1.05–1.33× on Sculptor: about level in green, up to a third higher in red and blue. Part of that is the demosaicing filter: after calibration and debayering, a single OriginStack frame is about 9–10% noisier than Siril's in green and blue and about 11% quieter in red (Malvar-He-Cutler against Siril's RCD). The rest traces to per-frame cosmic-ray/hot-pixel spikes: Malvar's debayer smears a single-pixel spike over more pixels in red/blue than in green, and the pipeline skips per-frame cosmic-ray detection once a session is large enough to rely on stack-level rejection instead — which runs after debayering, on the already-smeared pixels, and doesn't fully make up for it. Forcing `--cosmic-ray-rejection` back on closes most of the gap (measured 11–16% quieter across all three sessions) but softens stars slightly on shorter sub-exposures and adds 39–75% to the stacking time, so it stays opt-in rather than the default — see the changelog for the numbers. -- **OriginStack's memory does not grow with the session.** It stayed at 16.5–16.9 GB from 114 to 532 frames, set by the worker count (`-j` lowers it), while Siril's grew from 8.1 to 11.2 GB. Siril uses less at these sizes. Peak temporary disk was about the same for both (18–59 GB). +- **Speed depends on the session, and now on your sub length.** Omega's stack-only time is now 102 s against Siril's 63 s (1.6× as long) — slower than an earlier version of this section reported, because Omega's 30 s subs now cross this version's threshold for keeping per-frame cosmic-ray detection on through a big session (previously skipped by default above ~20 frames); that's most of the extra time, and it's also what closes most of the noise gap below. Sunflower and Sculptor's shorter subs (20 s, 10 s) fall under that threshold and keep the old default, but this version's finishing step (Phase 4) is roughly twice as fast for every session regardless of sub length — their figures below are from the previous measurement and have not been re-run since this pass. Timings move by up to 15% between runs of the same code: Siril's Omega ranged from 53 to 72 s across our runs, 63 s the latest; OriginStack's Sunflower from 83 to 106 s. +- **Noise: closer than it was, on sessions long enough to earn it.** On Omega, OriginStack's per-pixel noise is now 0.87–1.03× Siril's (was 0.99–1.19×) — the same 30 s-sub cosmic-ray fix above, not a separate change. Sunflower and Sculptor's subs are too short to trigger it by default, so their figures (1.09–1.20× on Sunflower, 1.05–1.33× on Sculptor) stand from the earlier measurement, about level in green, up to a third higher in red and blue. Part of the underlying gap is the demosaicing filter itself: after calibration and debayering, a single OriginStack frame is about 9–10% noisier than Siril's in green and blue and about 11% quieter in red (Malvar-He-Cutler against Siril's RCD). The rest traces to per-frame cosmic-ray/hot-pixel spikes: Malvar's debayer smears a single-pixel spike over more pixels in red/blue than in green, and stack-level rejection alone runs after debayering, on the already-smeared pixels, so it doesn't fully make up for it — per-frame detection catches it earlier, at the time cost above, which is why it's a threshold (`--cosmic-ray-rejection` / `--no-cosmic-ray-rejection` override it either way) rather than always on. See the changelog for the full numbers. +- **OriginStack's memory does not grow with the session.** It ranged 16.7–18.4 GB across the three sessions (Omega's cosmic-ray pass adds some of that), set by the worker count (`-j` lowers it), while Siril's grew from 8.1 to 11.2 GB. Siril uses less at these sizes. Peak temporary disk was about the same for both (18–59 GB). - **DeepSkyStacker** took 14 min 8 s on Omega (defaults: plain average, no rejection), gave stars softer than both others (4.83 px against 4.56 for Siril on the same stars), and its stack showed horizontal hot-pixel streaks and a few misregistered frames. A tuned run would look better. -- **Where OriginStack fits.** It finishes the image in the same run with no setup, explains the night ([diagnostics](docs/advanced.html)), runs fully offline on request (`--offline`) and does much that the other tools leave to you. On a plain stack it is now competitive on sharpness and behind Siril on noise and, for larger sessions, on time. +- **Where OriginStack fits.** It finishes the image in the same run with no setup, explains the night ([diagnostics](docs/advanced.html)), runs fully offline on request (`--offline`) and does much that the other tools leave to you. On a plain stack it is now competitive on sharpness, close to Siril on noise on longer subs where per-frame cosmic-ray detection kicks in, further behind on shorter subs where it doesn't, and behind Siril on time. + +Omega's figures above were refreshed on 2026-09-22 after a performance pass (Phase 4 roughly twice as fast; a new threshold that keeps per-frame cosmic-ray detection on through long sessions with 25 s+ subs, slower but quieter). Sunflower and Sculptor still reflect the prior measurement and are due for a re-run. How star width is measured matters more than it looks. It is a Gaussian fit to each of up to 400 bright, unsaturated, isolated stars, **at the same positions in both stacks and on each stack's own pixel grid** (no resampling), and the median is reported. An earlier version of this section compared each stack's FWHM over its own detected star list and concluded OriginStack's stars were tighter; that came from which stars were picked and was withdrawn. Measured properly, OriginStack's stars were then genuinely 13% wider than Siril's on Omega, until a hot-pixel filter that was clipping the cores of bright stars was found, by switching Phase 1 steps off one at a time, and fixed (see the changelog). Registration accuracy (about 0.1 px between frames), the warp kernel, rejection, weighting and normalisation were all checked and are not the cause of anything above. diff --git a/docs/index.html b/docs/index.html index 10f79d4..ca5ef3d 100644 --- a/docs/index.html +++ b/docs/index.html @@ -86,22 +86,22 @@

114 frames, 30 seconds each, one image.

-

A night of frames in under two minutes

-

The Omega Nebula, 114 frames, on an 8-core desktop: 113 seconds from the raw files to a finished image. The red block is the part that makes the picture look good; the rest is getting 114 frames to line up.

+

A night of frames in about two minutes

+

The Omega Nebula, 114 frames, on an 8-core desktop: 132 seconds from the raw files to a finished image. The red block is the part that makes the picture look good; the rest is getting 114 frames to line up. These 30-second subs cross this version's threshold for keeping per-frame cosmic-ray detection on through a big session, which is most of why loading takes longer than it used to — it also quiets the stack, see the comparison below.

-
From 206e3c4fedec8c27f909137ad9a8dab9df8bb068 Mon Sep 17 00:00:00 2001 From: Hans Davenport <35202271+hd152@users.noreply.github.com> Date: Tue, 22 Sep 2026 11:22:07 -0700 Subject: [PATCH 3/4] Add native Gaussian blur kernel, wire into highest-traffic callers correlate1d (the C function behind scipy.ndimage.gaussian_filter) was the single largest self-time item in every full-pipeline profile taken this session, spread across ~30 call sites project-wide. gaussian_filter_native is a from-scratch separable reimplementation of scipy's default mode='reflect' behaviour (two rayon-parallel passes, same reflect_idx boundary convention the hot-pixel kernels already use). Not a port, so parity is numerical rather than bit-exact: verified to double-precision rounding (rtol/atol 1e-9) across sigma 0.8-96 and two shapes. ~4x faster at small sigma, ~1.3-1.5x at the large sigmas gaussian_filter_ds's downsampled branch uses, measured in isolation. Wired via a new _gaussian_blur wrapper (2-D float only, transparent scipy fallback otherwise) into gaussian_filter_ds's both branches and into background.py/denoising.py's highest-traffic direct callers (reduce_chroma_noise, multiscale_local_contrast, _structure_tensor_coherence, DBE's mesh-fallback surface fit). Deliberately not a full sweep: registration.py, quality.py, moving_objects.py, noise_validation.py and a few others still call scipy directly, each with its own sigma/shape assumptions worth checking before switching over. A 3-D per-axis-sigma call in reduce_stars (sigma=(s,s,0)) was left alone on purpose -- the kernel only supports scalar sigma on a 2-D plane. Measured on a controlled 15-frame real-session subset, same machine, back to back: 46.0s -> 44.7s. The full 114-frame real benchmark run showed no reliable signal either way across repeated runs (system-level timing noise on this machine outweighed the effect at that scale, confirmed via idle CPU load between runs) -- not used as evidence here. astro_native bumped to 0.34.0. Co-Authored-By: Claude Sonnet 5 --- CHANGELOG.md | 13 ++++ CLAUDE.md | 1 + ext/astro_native/Cargo.lock | 2 +- ext/astro_native/Cargo.toml | 2 +- ext/astro_native/pyproject.toml | 2 +- ext/astro_native/src/lib.rs | 101 ++++++++++++++++++++++++++++++++ src/background.py | 38 +++++++++--- src/denoising.py | 22 +++---- tests/test_native.py | 71 ++++++++++++++++++++++ 9 files changed, 230 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9906b6c..2c40fb9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,19 @@ match the `VERSION` file and `v*` git tags. ### Changed +- **A native Gaussian blur kernel, wired into the highest-traffic callers.** `correlate1d` (the C function + behind `scipy.ndimage.gaussian_filter`) was the single largest self-time item in every full-pipeline + profile taken this session, spread across ~30 call sites project-wide (DBE, chroma denoising, local + contrast, structure-tensor coherence, edge-band correction, and more). Added `gaussian_filter_native` + (from-scratch separable reimplementation, `mode='reflect'`, verified against scipy to double-precision + rounding) and wired it into `background.py`'s `gaussian_filter_ds` plus its own and `denoising.py`'s + highest-traffic direct callers (`reduce_chroma_noise`, `multiscale_local_contrast`, + `_structure_tensor_coherence`, DBE's mesh-fallback surface fit) -- ~4x faster at small sigma, ~1.3-1.5x at + the large sigmas `gaussian_filter_ds`'s downsampled branch uses, measured in isolation. This is not a full + sweep of every direct call site (registration.py, quality.py, moving_objects.py, noise_validation.py and a + few others still call scipy directly) -- each has its own sigma/shape assumptions worth checking before + switching over, left for a follow-up. + - **Several real perf fixes, found by profiling a full run rather than guessing.** A synthetic-flat build (`--flat-from-lights`, auto-triggered whenever no flat frames exist) was silently the single biggest cost on a real profiled session -- 126.9s of a 187.8s total, hidden inside an unlabeled "Other (I/O)" bucket in diff --git a/CLAUDE.md b/CLAUDE.md index 44365e7..5b6c5c0 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -234,6 +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) is wired into `gaussian_filter_ds`'s both branches and into `background.py`/`denoising.py`'s highest-traffic *direct* `gaussian_filter` callers (`reduce_chroma_noise`, `multiscale_local_contrast`, `_structure_tensor_coherence`, DBE's mesh-fallback surface fit) -- **not a full sweep**: registration.py's DOG/smoothed-luminance calls, quality.py's star-mask blur, moving_objects.py, noise_validation.py, merge.py and a handful of others still call scipy directly, each with its own sigma/shape assumptions worth checking individually before switching over (a 3-D per-axis-sigma call in `denoising.py::reduce_stars`, `sigma=(s,s,0)`, was deliberately left alone -- the kernel only supports a scalar sigma on a 2-D plane). - **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/ext/astro_native/Cargo.lock b/ext/astro_native/Cargo.lock index 298b8d1..c6fdce9 100644 --- a/ext/astro_native/Cargo.lock +++ b/ext/astro_native/Cargo.lock @@ -49,7 +49,7 @@ checksum = "fb5dfbc6d8d2675589ccbe4d0fd61df2419075625f8c1a62325e718e2b0049f9" [[package]] name = "astro_native" -version = "0.33.0" +version = "0.34.0" dependencies = [ "numpy", "pyo3", diff --git a/ext/astro_native/Cargo.toml b/ext/astro_native/Cargo.toml index 0f04e3c..5e0f08c 100644 --- a/ext/astro_native/Cargo.toml +++ b/ext/astro_native/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "astro_native" -version = "0.33.0" +version = "0.34.0" edition = "2021" description = "Native (Rust) hot-path kernels for OriginStack: stacking combine, etc." diff --git a/ext/astro_native/pyproject.toml b/ext/astro_native/pyproject.toml index 96aeb8c..d14a42a 100644 --- a/ext/astro_native/pyproject.toml +++ b/ext/astro_native/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "astro_native" -version = "0.33.0" +version = "0.34.0" description = "Native Rust hot-path kernels for OriginStack" requires-python = ">=3.10" classifiers = ["Programming Language :: Rust"] diff --git a/ext/astro_native/src/lib.rs b/ext/astro_native/src/lib.rs index cad77a4..db59706 100644 --- a/ext/astro_native/src/lib.rs +++ b/ext/astro_native/src/lib.rs @@ -2861,6 +2861,106 @@ fn median_filter_native<'py>( .into_pyarray(py)) } +/// `scipy.ndimage.gaussian_filter1d`'s exact 1D kernel: normalized samples of +/// the Gaussian PDF over `[-radius, radius]`, `radius = floor(truncate*sigma +/// + 0.5)` -- same formula scipy uses, order-0 (no derivative). +fn gaussian_kernel1d(sigma: f64, truncate: f64) -> Vec { + let radius = (truncate * sigma + 0.5) as isize; + let sigma2 = sigma * sigma; + let mut w: Vec = (-radius..=radius) + .map(|x| { + let xf = x as f64; + (-0.5 * xf * xf / sigma2).exp() + }) + .collect(); + let sum: f64 = w.iter().sum(); + for v in w.iter_mut() { + *v /= sum; + } + w +} + +/// Separable Gaussian blur matching `scipy.ndimage.gaussian_filter(data, +/// sigma, mode='reflect')` (the default mode, and the only one any call site +/// in this codebase uses) for a scalar `sigma` on a 2D array. Two rayon- +/// parallel correlate1d-style passes (row-wise, then column-wise) instead of +/// scipy's generic N-D machinery -- found by profiling a real session: +/// `correlate1d` (the C function `gaussian_filter1d` calls per axis) was the +/// single largest self-time item in every profile taken this session, spread +/// across ~30 call sites project-wide (DBE, chroma denoising, local +/// contrast, structure-tensor coherence, star-mask generation, registration, +/// edge-band correction, and more). This kernel is wired into +/// `background.py`'s `gaussian_filter_ds` (already the single most-reused +/// choke point among those callers) as a first step, not a full sweep of +/// every direct `scipy.ndimage.gaussian_filter` call site -- each of those +/// has its own sigma/shape assumptions worth checking individually before +/// switching it over. Boundary handling reuses `reflect_idx` (scipy's +/// edge-duplicating `mode='reflect'`, not numpy's non-duplicating one). +#[pyfunction] +#[pyo3(signature = (data, sigma, truncate=4.0))] +fn gaussian_filter_native<'py>( + py: Python<'py>, + data: PyReadonlyArray2<'py, f64>, + sigma: f64, + truncate: f64, +) -> PyResult>> { + let arr = data.as_array(); + let s = arr.shape(); + let (h, w) = (s[0], s[1]); + let owned: Vec; + let flat: &[f64] = match arr.as_slice() { + Some(sl) => sl, + None => { + owned = arr.iter().copied().collect(); + &owned + } + }; + + if sigma <= 0.0 { + let arr2 = numpy::ndarray::Array2::from_shape_vec((h, w), flat.to_vec()) + .expect("shape mismatch building gaussian_filter_native passthrough"); + return Ok(arr2.into_pyarray(py)); + } + + let kernel = gaussian_kernel1d(sigma, truncate); + let radius = (kernel.len() / 2) as isize; + + let out = py.detach(|| { + // Pass 1: blur along axis 1 (each row independently), parallel over rows. + let mut tmp = vec![0f64; h * w]; + tmp.par_chunks_mut(w).enumerate().for_each(|(y, out_row)| { + let row = &flat[y * w..(y + 1) * w]; + for x in 0..w { + let mut acc = 0.0f64; + for (k, &kv) in kernel.iter().enumerate() { + let dx = k as isize - radius; + let xi = reflect_idx(x as isize + dx, w); + acc += kv * row[xi]; + } + out_row[x] = acc; + } + }); + // Pass 2: blur along axis 0 (each column), parallel over output rows. + let mut out = vec![0f64; h * w]; + out.par_chunks_mut(w).enumerate().for_each(|(y, out_row)| { + for x in 0..w { + let mut acc = 0.0f64; + for (k, &kv) in kernel.iter().enumerate() { + let dy = k as isize - radius; + let yi = reflect_idx(y as isize + dy, h); + acc += kv * tmp[yi * w + x]; + } + out_row[x] = acc; + } + }); + out + }); + + let arr2 = numpy::ndarray::Array2::from_shape_vec((h, w), out) + .expect("shape mismatch building gaussian_filter_native output"); + Ok(arr2.into_pyarray(py)) +} + // --------------------------------------------------------------------------- // DBE robust background-surface fit // --------------------------------------------------------------------------- @@ -7687,6 +7787,7 @@ fn astro_native(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(anisotropic_diffusion, m)?)?; m.add_function(wrap_pyfunction!(lacosmic_reject_native, m)?)?; m.add_function(wrap_pyfunction!(median_filter_native, m)?)?; + m.add_function(wrap_pyfunction!(gaussian_filter_native, m)?)?; m.add_function(wrap_pyfunction!(dbe_fit_surface, m)?)?; m.add_function(wrap_pyfunction!(dbe_sample_patches, m)?)?; m.add_function(wrap_pyfunction!(detect_stars_matched_filter, m)?)?; diff --git a/src/background.py b/src/background.py index 860a731..11bbb51 100644 --- a/src/background.py +++ b/src/background.py @@ -51,6 +51,28 @@ def safe_print(*args, **kwargs): _gf = None +_HAS_NATIVE_GAUSSIAN = HAS_NATIVE and hasattr(_native, 'gaussian_filter_native') + + +def _gaussian_blur(a: np.ndarray, sigma: float) -> np.ndarray: + """``gaussian_filter(a, sigma=sigma)`` (mode='reflect', the default and + the only mode any caller here uses) routed through the native + ``gaussian_filter_native`` kernel when available. `correlate1d` (the C + function scipy's gaussian_filter1d calls per axis) was the single + largest self-time item in every profile taken of a real session -- + ~4x faster at small sigma, ~1.3-1.5x at the large sigmas + ``gaussian_filter_ds``'s downsampled branch uses, verified within + double-precision rounding of scipy's own result (`tests/test_native.py`). + 2-D float input only; anything else falls back to scipy untouched.""" + if _HAS_NATIVE_GAUSSIAN and a.ndim == 2: + try: + return np.asarray(_native.gaussian_filter_native(np.ascontiguousarray(a, dtype=np.float64), + float(sigma))) + except Exception: + pass + return gaussian_filter(a, sigma=sigma) + + 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. @@ -67,14 +89,14 @@ def gaussian_filter_ds(arr: np.ndarray, sigma: float, sigma = float(sigma) a = np.asarray(arr, dtype=np.float64) if sigma < ds_threshold or a.ndim != 2: - return gaussian_filter(a, sigma=sigma) + return _gaussian_blur(a, sigma) ds = 4 if sigma < 96.0 else 8 H, W = a.shape h2, w2 = (H // ds) * ds, (W // ds) * ds if h2 < ds or w2 < ds: - return gaussian_filter(a, sigma=sigma) + return _gaussian_blur(a, sigma) coarse = a[:h2, :w2].reshape(h2 // ds, ds, w2 // ds, ds).mean(axis=(1, 3)) - sm = gaussian_filter(coarse, sigma=sigma / ds) + sm = _gaussian_blur(coarse, sigma / ds) # Centre-aligned upsample: coarse pixel j represents fine pixels # [j*ds, (j+1)*ds), whose centre is j*ds + (ds-1)/2. Plain zoom is # corner-aligned and would return the whole field shifted by ~ds/2 px, @@ -322,7 +344,7 @@ def extract_background(img: np.ndarray, mesh_size: int = 256, filter_size: int = # Gaussian smooth if min(ny, nx) >= 4: - bg_grid = ndimage.gaussian_filter(bg_grid.astype(np.float64), sigma=0.8) + bg_grid = _gaussian_blur(bg_grid.astype(np.float64), 0.8) # --- Interpolation to Full Res --- grid_y = (np.arange(ny) + 0.5) * cell_h @@ -372,7 +394,7 @@ def extract_background(img: np.ndarray, mesh_size: int = 256, filter_size: int = # Final Gaussian blur to suppress high-frequency mesh ripple blur_sigma = cell_h * 0.5 if blur_sigma > 0: - background = ndimage.gaussian_filter(background.astype(np.float64), sigma=blur_sigma) + background = _gaussian_blur(background.astype(np.float64), blur_sigma) return np.asarray(background, dtype=np.float32) @@ -647,7 +669,7 @@ def _process_channel(c): bg_grid = ndimage.median_filter(bg_grid, size=filter_size) if min(ny, nx) >= 4: - bg_grid = ndimage.gaussian_filter(bg_grid.astype(np.float64), sigma=0.8) + bg_grid = _gaussian_blur(bg_grid.astype(np.float64), 0.8) grid_y = (np.arange(ny) + 0.5) * (H / ny) grid_x = (np.arange(nx) + 0.5) * (W / nx) @@ -1320,7 +1342,7 @@ def _fit_background_surface(coords: np.ndarray, values: np.ndarray, coords, values, H, W, outlier_sigma, max_iter, sigma_px, Hc, Wc, verbose) surface = zoom(coarse, (H / Hc, W / Wc), order=3)[:H, :W] - surface = ndimage.gaussian_filter(surface, sigma=patch_size * 0.5) + surface = _gaussian_blur(surface, patch_size * 0.5) return np.clip(surface, surf_lo, surf_hi) @@ -1351,7 +1373,7 @@ def poly2(y, x): xx = np.linspace(0.0, 1.0, W) grid_y, grid_x = np.meshgrid(yy, xx, indexing='ij') surface = poly(grid_y.ravel(), grid_x.ravel()).dot(coeffs).reshape(H, W) - return ndimage.gaussian_filter(surface, sigma=patch_size * 0.5) + return _gaussian_blur(surface, patch_size * 0.5) def dynamic_background_extraction( diff --git a/src/denoising.py b/src/denoising.py index 5c60434..5138601 100644 --- a/src/denoising.py +++ b/src/denoising.py @@ -7,7 +7,7 @@ from scipy import ndimage from src import wavelet -from src.background import _estimate_sky_sigma, gaussian_filter_ds +from src.background import _estimate_sky_sigma, _gaussian_blur, gaussian_filter_ds from src.models import Config from src.utils import get_logger, safe_print @@ -60,9 +60,9 @@ def _structure_tensor_coherence(plane: np.ndarray, sigma: float = 1.5) -> np.nda *magnitude* alone, not orientation coherence). """ gy, gx = np.gradient(plane.astype(np.float64)) - jxx = ndimage.gaussian_filter(gx * gx, sigma) - jyy = ndimage.gaussian_filter(gy * gy, sigma) - jxy = ndimage.gaussian_filter(gx * gy, sigma) + jxx = _gaussian_blur(gx * gx, sigma) + jyy = _gaussian_blur(gy * gy, sigma) + jxy = _gaussian_blur(gx * gy, sigma) trace = jxx + jyy disc = np.sqrt(np.maximum((jxx - jyy) ** 2 + 4 * jxy ** 2, 0.0)) lam1 = 0.5 * (trace + disc) @@ -366,7 +366,7 @@ def acdnr_denoise(img: np.ndarray, smoothing_sigma: float = 1.5, return img.copy() # Local contrast map at the smoothing scale - smooth_luma = ndimage.gaussian_filter(luma, sigma=smoothing_sigma) + smooth_luma = _gaussian_blur(luma, smoothing_sigma) contrast = np.abs(luma - smooth_luma) # YCbCr split (same BT.601 coefficients as the other denoisers) @@ -381,9 +381,9 @@ def acdnr_denoise(img: np.ndarray, smoothing_sigma: float = 1.5, chroma_w = np.exp(-0.5 * (contrast / chroma_thr) ** 2) # Smooth each YCbCr plane, then adaptively blend - Y_smooth = ndimage.gaussian_filter(Y, sigma=smoothing_sigma) - Cb_smooth = ndimage.gaussian_filter(Cb, sigma=smoothing_sigma) - Cr_smooth = ndimage.gaussian_filter(Cr, sigma=smoothing_sigma) + Y_smooth = _gaussian_blur(Y, smoothing_sigma) + Cb_smooth = _gaussian_blur(Cb, smoothing_sigma) + Cr_smooth = _gaussian_blur(Cr, smoothing_sigma) Y_d = luma_w * Y_smooth + (1.0 - luma_w) * Y Cb_d = chroma_w * Cb_smooth + (1.0 - chroma_w) * Cb @@ -475,7 +475,7 @@ def reduce_chroma_noise(img: np.ndarray, sigma: float = 2.0, sky_mask = 1.0 - protect # float [0,1] result = np.empty_like(img, dtype=np.float64) - blurred_weight = ndimage.gaussian_filter(sky_mask, sigma=sigma) + blurred_weight = _gaussian_blur(sky_mask, sigma) safe_weight = np.maximum(blurred_weight, 1e-9) # Optional coarse pass — smooths medium-scale colour blotches (walking / @@ -491,7 +491,7 @@ def reduce_chroma_noise(img: np.ndarray, sigma: float = 2.0, for c in range(img.shape[2]): chroma = img[:, :, c].astype(np.float64) - lum # Weighted blur: star pixels contribute 0, background contributes 1 - smooth_chroma = ndimage.gaussian_filter(chroma * sky_mask, sigma=sigma) / safe_weight + smooth_chroma = _gaussian_blur(chroma * sky_mask, sigma) / safe_weight # Stars keep original chroma; background gets smoothed chroma out_chroma = chroma * protect + smooth_chroma * sky_mask if do_large: @@ -688,7 +688,7 @@ def multiscale_local_contrast( for sigma, w in zip(scales, scale_weights): if w <= 0 or strength <= 0: continue - blurred = ndimage.gaussian_filter(detail_src, sigma=float(sigma)) + blurred = _gaussian_blur(detail_src, float(sigma)) detail = detail_src - blurred # high-frequency detail at this scale enhanced_lum += strength * w * detail * mask diff --git a/tests/test_native.py b/tests/test_native.py index 5cf2b9c..96209e3 100644 --- a/tests/test_native.py +++ b/tests/test_native.py @@ -7,6 +7,7 @@ import pytest import originstack as astro +import src.background as _background_mod import src.blind_match as _blind_match_mod import src.channel_combine as _channel_combine_mod import src.debayer as _debayer_mod @@ -520,6 +521,76 @@ def test_median_filter_native_matches_scipy(size): assert float(np.max(np.abs(ref.astype(np.float64) - got.astype(np.float64)))) < 1e-4 +@pytest.mark.parametrize("sigma", [0.8, 2.0, 5.0, 24.0, 32.0]) +@pytest.mark.parametrize("shape", [(80, 96), (301, 257)]) +def test_gaussian_filter_native_matches_scipy(sigma, shape): + """gaussian_filter_native is a from-scratch separable reimplementation of + scipy.ndimage.gaussian_filter's default mode='reflect' -- not a port, so + parity is judged by numerical agreement, not shared code. Real-shape + profiling showed correlate1d (gaussian_filter1d's C function) as the top + self-time item in every profile taken of a full pipeline run this + session, spread across ~30 call sites; this kernel and its wiring into + background.py's _gaussian_blur / gaussian_filter_ds don't sweep all of + them, just the highest-traffic ones.""" + from scipy.ndimage import gaussian_filter + rng = np.random.default_rng(4) + a = rng.normal(1000.0, 50.0, shape) + want = gaussian_filter(a, sigma=sigma) + got = np.asarray(native.gaussian_filter_native(a, sigma)) + # separable-pass floating-point summation order differs from scipy's; + # both are exact renditions of the same closed-form kernel, so the gap + # is double-precision rounding, not an approximation choice + np.testing.assert_allclose(got, want, rtol=1e-9, atol=1e-9) + + +def test_gaussian_filter_native_zero_sigma_is_passthrough(): + rng = np.random.default_rng(5) + a = rng.normal(size=(20, 30)) + got = np.asarray(native.gaussian_filter_native(a, 0.0)) + np.testing.assert_array_equal(got, a) + + +def test_gaussian_filter_native_rejects_non_2d_by_signature(): + # the numpy binding itself enforces 2D; the Python-side _gaussian_blur + # wrapper is what actually gates this in production (see + # test_gaussian_blur_falls_back_for_non_2d below) + rng = np.random.default_rng(6) + a = rng.normal(size=(10, 10, 3)) + with pytest.raises(Exception): + native.gaussian_filter_native(a, 2.0) + + +def test_gaussian_blur_wrapper_matches_scipy_native_and_fallback(): + """background.py's _gaussian_blur (the wrapper wired into + gaussian_filter_ds and swapped into ~10 direct call sites in background.py + and denoising.py) must agree with plain scipy whether or not native is + available.""" + from scipy.ndimage import gaussian_filter + rng = np.random.default_rng(7) + a = rng.normal(1000.0, 50.0, (150, 200)) + want = gaussian_filter(a, sigma=5.0) + + got_native = _background_mod._gaussian_blur(a, 5.0) + np.testing.assert_allclose(got_native, want, rtol=1e-9, atol=1e-9) + + had = _background_mod._HAS_NATIVE_GAUSSIAN + _background_mod._HAS_NATIVE_GAUSSIAN = False + try: + got_fallback = _background_mod._gaussian_blur(a, 5.0) + finally: + _background_mod._HAS_NATIVE_GAUSSIAN = had + np.testing.assert_array_equal(got_fallback, want) + + +def test_gaussian_blur_falls_back_for_non_2d(): + from scipy.ndimage import gaussian_filter + rng = np.random.default_rng(8) + a = rng.normal(size=(20, 30, 3)) + want = gaussian_filter(a, sigma=2.0) + got = _background_mod._gaussian_blur(a, 2.0) + np.testing.assert_array_equal(got, want) + + 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 From 243cdbdb5da1a13b5a5ed7a672486f8f73367301 Mon Sep 17 00:00:00 2001 From: Hans Davenport <35202271+hd152@users.noreply.github.com> Date: Tue, 22 Sep 2026 11:36:35 -0700 Subject: [PATCH 4/4] Sweep the native Gaussian blur kernel across remaining call sites Follow-up to the previous commit, which wired gaussian_filter_native into only the highest-traffic callers and flagged the rest for later. Swept every remaining direct scipy.ndimage.gaussian_filter call that uses a plain scalar sigma on a 2-D array: cli.py (master bias/dark/flat 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 on purpose: 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 with its own tight numeric tolerances against a reference implementation) -- the kernel only supports a scalar sigma on a 2-D plane. These already fall back to scipy safely on their own (a tuple sigma raises inside the wrapper's float(sigma) cast, caught by the existing try/except), so nothing needed changing there. Also fixed a real bug in _gaussian_blur found while auditing the sweep: it didn't preserve the caller's dtype, so a float32 array (several master-calibration and Phase 4 call sites pass float32) silently came back as float64, doubling memory for no reason. The native kernel is always float64 internally like every other kernel in this file; the wrapper now casts the result back to the input's dtype, matching scipy's own contract. New dtype-parity test in test_native.py catches this going forward. Full test suite (1647 tests) passes; smoke-tested --trail-reject and --noise-validate (two of the newly-wired opt-in paths) against a real 15-frame session subset with no crash. Co-Authored-By: Claude Sonnet 5 --- CHANGELOG.md | 28 ++++++++++++++++------------ CLAUDE.md | 2 +- src/background.py | 13 ++++++++++--- src/cli.py | 8 ++++---- src/exposure_fusion.py | 4 +++- src/merge.py | 6 +++--- src/moving_objects.py | 5 ++++- src/noise_validation.py | 7 +++++-- src/postprocess.py | 5 +++-- src/quality.py | 6 ++++-- src/registration.py | 17 ++++++++++------- src/star_removal.py | 7 +++++-- src/trail_reject.py | 7 +++++-- tests/test_native.py | 18 ++++++++++++++++++ 14 files changed, 91 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2c40fb9..8001fac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,18 +8,22 @@ match the `VERSION` file and `v*` git tags. ### Changed -- **A native Gaussian blur kernel, wired into the highest-traffic callers.** `correlate1d` (the C function - behind `scipy.ndimage.gaussian_filter`) was the single largest self-time item in every full-pipeline - profile taken this session, spread across ~30 call sites project-wide (DBE, chroma denoising, local - contrast, structure-tensor coherence, edge-band correction, and more). Added `gaussian_filter_native` - (from-scratch separable reimplementation, `mode='reflect'`, verified against scipy to double-precision - rounding) and wired it into `background.py`'s `gaussian_filter_ds` plus its own and `denoising.py`'s - highest-traffic direct callers (`reduce_chroma_noise`, `multiscale_local_contrast`, - `_structure_tensor_coherence`, DBE's mesh-fallback surface fit) -- ~4x faster at small sigma, ~1.3-1.5x at - the large sigmas `gaussian_filter_ds`'s downsampled branch uses, measured in isolation. This is not a full - sweep of every direct call site (registration.py, quality.py, moving_objects.py, noise_validation.py and a - few others still call scipy directly) -- each has its own sigma/shape assumptions worth checking before - switching over, left for a follow-up. +- **A native Gaussian blur kernel, swept across essentially every 2-D scalar-sigma call site.** `correlate1d` + (the C function behind `scipy.ndimage.gaussian_filter`) was the single largest self-time item in every + full-pipeline profile taken this session, spread across ~30 call sites project-wide. Added + `gaussian_filter_native` (from-scratch separable reimplementation, `mode='reflect'`, verified against scipy + to double-precision rounding) behind a `_gaussian_blur` wrapper (`src/background.py`) that preserves the + caller's dtype (float32 stays float32 -- a first version of the wrapper silently upcast to float64, caught + by a new dtype-parity test before it could inflate memory on every caller) and falls back to scipy for + anything the kernel doesn't cover. ~4x faster at small sigma, ~1.3-1.5x at large. Wired into every direct + `gaussian_filter` call site that uses a plain scalar sigma on a 2-D array: `background.py`, + `denoising.py`, `cli.py` (master bias/dark/flat smoothing), `exposure_fusion.py`, `merge.py`, + `moving_objects.py`, `noise_validation.py`, `postprocess.py`, `quality.py`, `registration.py`, + `star_removal.py`, `trail_reject.py`. Deliberately left alone: the handful of call sites that pass a + per-axis sigma tuple (e.g. `(sigma, sigma, 0)` to blur spatial axes only, or `--fix-atmospheric-dispersion` + -adjacent inference preprocessing with its own tight numeric tolerances against a reference implementation) + -- the kernel only supports a scalar sigma on a 2-D plane, and those calls already fall back to scipy + safely on their own (a tuple sigma raises inside the wrapper's `float(sigma)` cast, caught and handled). - **Several real perf fixes, found by profiling a full run rather than guessing.** A synthetic-flat build (`--flat-from-lights`, auto-triggered whenever no flat frames exist) was silently the single biggest cost diff --git a/CLAUDE.md b/CLAUDE.md index 5b6c5c0..097d5f0 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) is wired into `gaussian_filter_ds`'s both branches and into `background.py`/`denoising.py`'s highest-traffic *direct* `gaussian_filter` callers (`reduce_chroma_noise`, `multiscale_local_contrast`, `_structure_tensor_coherence`, DBE's mesh-fallback surface fit) -- **not a full sweep**: registration.py's DOG/smoothed-luminance calls, quality.py's star-mask blur, moving_objects.py, noise_validation.py, merge.py and a handful of others still call scipy directly, each with its own sigma/shape assumptions worth checking individually before switching over (a 3-D per-axis-sigma call in `denoising.py::reduce_stars`, `sigma=(s,s,0)`, was deliberately left alone -- the kernel only supports a scalar sigma on a 2-D plane). +- **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). - **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 11bbb51..14bbb9c 100644 --- a/src/background.py +++ b/src/background.py @@ -63,11 +63,18 @@ def _gaussian_blur(a: np.ndarray, sigma: float) -> np.ndarray: ~4x faster at small sigma, ~1.3-1.5x at the large sigmas ``gaussian_filter_ds``'s downsampled branch uses, verified within double-precision rounding of scipy's own result (`tests/test_native.py`). - 2-D float input only; anything else falls back to scipy untouched.""" + 2-D float input only; anything else falls back to scipy untouched. + + Output dtype always matches the input's, same as scipy -- the native + kernel always computes in float64 internally (like every other kernel in + this file), so a float32 input is cast back down after the blur rather + than silently widening every caller's array to float64. + """ if _HAS_NATIVE_GAUSSIAN and a.ndim == 2: try: - return np.asarray(_native.gaussian_filter_native(np.ascontiguousarray(a, dtype=np.float64), - float(sigma))) + out = np.asarray(_native.gaussian_filter_native( + np.ascontiguousarray(a, dtype=np.float64), float(sigma))) + return out.astype(a.dtype, copy=False) if out.dtype != a.dtype else out except Exception: pass return gaussian_filter(a, sigma=sigma) diff --git a/src/cli.py b/src/cli.py index 5efd928..62c50b5 100644 --- a/src/cli.py +++ b/src/cli.py @@ -15,8 +15,8 @@ import numpy as np from astropy.io import fits -from scipy import ndimage +from src.background import _gaussian_blur from src.cleanup import deregister as _cleanup_deregister from src.cleanup import register as _cleanup_register from src.debayer import build_hot_pixel_map @@ -511,18 +511,18 @@ def _method_tag(method: str) -> str: if masters.get('bias') is not None: n_bias = len(frames['bias']) sigma_b = max(1, 30 // max(1, int(np.sqrt(n_bias)))) - masters['bias'] = ndimage.gaussian_filter(masters['bias'].astype(np.float32), sigma=sigma_b) + masters['bias'] = _gaussian_blur(masters['bias'].astype(np.float32), sigma_b) if masters.get('dark') is not None: n_dark = len(frames['dark']) sigma_d = max(1, 20 // max(1, int(np.sqrt(n_dark)))) - masters['dark'] = ndimage.gaussian_filter(masters['dark'].astype(np.float32), sigma=sigma_d) + masters['dark'] = _gaussian_blur(masters['dark'].astype(np.float32), sigma_d) if masters.get('flat') is not None: n_flat = len(frames['flat']) sigma_f = max(1, 15 // max(1, int(np.sqrt(n_flat)))) flat_raw = masters['flat'].astype(np.float32) for r_off, c_off in [(0, 0), (0, 1), (1, 0), (1, 1)]: ch = flat_raw[r_off::2, c_off::2] - flat_raw[r_off::2, c_off::2] = ndimage.gaussian_filter(ch, sigma=sigma_f) + flat_raw[r_off::2, c_off::2] = _gaussian_blur(ch, sigma_f) masters['flat'] = flat_raw masters['dark_exptime'] = None diff --git a/src/exposure_fusion.py b/src/exposure_fusion.py index fdbbc9b..3fdaf05 100644 --- a/src/exposure_fusion.py +++ b/src/exposure_fusion.py @@ -23,6 +23,8 @@ import numpy as np from scipy.ndimage import gaussian_filter, zoom +from src.background import _gaussian_blur + def _fit_shape(x: np.ndarray, target_shape) -> np.ndarray: """Pad (edge) or crop x's first two axes to exactly target_shape[:2] -- @@ -93,7 +95,7 @@ def _quality_weights(norm_img: np.ndarray, contrast_w: float, saturation_w: floa gray = norm_img.mean(axis=-1) # Contrast: a difference-of-Gaussians magnitude, a simple standard # stand-in for the paper's own Laplacian-magnitude measure. - contrast = np.abs(gaussian_filter(gray, 1.0) - gaussian_filter(gray, 2.0)) + contrast = np.abs(_gaussian_blur(gray, 1.0) - _gaussian_blur(gray, 2.0)) # Saturation: std across channels -- low for a washed-out/near-grey pixel. saturation = norm_img.std(axis=-1) diff --git a/src/merge.py b/src/merge.py index a292e43..ed17d80 100644 --- a/src/merge.py +++ b/src/merge.py @@ -200,13 +200,13 @@ def _match_flux_scale(ref: np.ndarray, img: np.ndarray, valid: np.ndarray merge is not perturbed by estimator noise. Returns (gain[3], offset[3], pixels used in the worst-measured channel). """ - from scipy.ndimage import gaussian_filter + from src.background import _gaussian_blur gains = np.ones(3) offsets = np.zeros(3) n_min = None for c in range(3): - sr = gaussian_filter(ref[:, :, c].astype(np.float32), _SCALE_SMOOTH_SIGMA) - si = gaussian_filter(img[:, :, c].astype(np.float32), _SCALE_SMOOTH_SIGMA) + sr = _gaussian_blur(ref[:, :, c].astype(np.float32), _SCALE_SMOOTH_SIGMA) + si = _gaussian_blur(img[:, :, c].astype(np.float32), _SCALE_SMOOTH_SIGMA) vr, vi = sr[valid].astype(np.float64), si[valid].astype(np.float64) if vr.size < _SCALE_MIN_PIXELS: n_min = 0 diff --git a/src/moving_objects.py b/src/moving_objects.py index e5a51a9..b0b7df1 100644 --- a/src/moving_objects.py +++ b/src/moving_objects.py @@ -33,9 +33,12 @@ try: from scipy import ndimage + + from src.background import _gaussian_blur _HAS_SCIPY = True except Exception: # pragma: no cover ndimage = None + _gaussian_blur = None _HAS_SCIPY = False _log = logging.getLogger("originstack") @@ -82,7 +85,7 @@ def detect_in_residual(residual: np.ndarray, fwhm: float, threshold: float, The residual is smoothed at the PSF width and divided by its own robust noise (MAD), so the threshold is in sigma of the smoothed image.""" - sm = ndimage.gaussian_filter(residual.astype(np.float32), max(fwhm / 2.355, 0.8)) + sm = _gaussian_blur(residual.astype(np.float32), max(fwhm / 2.355, 0.8)) med = float(np.median(sm)) sig = 1.4826 * float(np.median(np.abs(sm - med))) if not np.isfinite(sig) or sig <= 0: diff --git a/src/noise_validation.py b/src/noise_validation.py index 57ed412..c20fc99 100644 --- a/src/noise_validation.py +++ b/src/noise_validation.py @@ -25,9 +25,12 @@ try: from scipy import ndimage + + from src.background import _gaussian_blur _HAS_SCIPY = True except Exception: # pragma: no cover ndimage = None + _gaussian_blur = None _HAS_SCIPY = False @@ -93,8 +96,8 @@ def consistency_map(a: np.ndarray, b: np.ndarray, smooth: float = 2.0, window: i def lum(x): return 0.299 * x[:, :, 0] + 0.587 * x[:, :, 1] + 0.114 * x[:, :, 2] - la = ndimage.gaussian_filter(lum(a).astype(np.float64), smooth) - lb = ndimage.gaussian_filter(lum(b).astype(np.float64), smooth) + la = _gaussian_blur(lum(a).astype(np.float64), smooth) + lb = _gaussian_blur(lum(b).astype(np.float64), smooth) m = lambda x: ndimage.uniform_filter(x, window, mode='nearest') ma, mb = m(la), m(lb) cov = m(la * lb) - ma * mb diff --git a/src/postprocess.py b/src/postprocess.py index 22456b8..5ef8c22 100644 --- a/src/postprocess.py +++ b/src/postprocess.py @@ -12,6 +12,7 @@ from src.background import ( _border_pixels, _dbe_prepare_emission_mask, + _gaussian_blur, apply_background_extraction, dynamic_background_extraction, gaussian_filter_ds, @@ -805,9 +806,9 @@ def postprocess_stack( # minimum. _ring_r = max(4, int(float(psf.shape[0])), int(np.ceil(psf_fwhm * 5.0))) _core = ndimage.binary_dilation(_pts, iterations=_ring_r) - deconv_mask = ndimage.gaussian_filter( + deconv_mask = _gaussian_blur( _core.astype(np.float32), - sigma=max(3.0, float(_ring_r) * 0.45)) + max(3.0, float(_ring_r) * 0.45)) _dmx = float(deconv_mask.max()) if _dmx > 0: deconv_mask /= _dmx diff --git a/src/quality.py b/src/quality.py index f4ade4a..e4b06d6 100644 --- a/src/quality.py +++ b/src/quality.py @@ -75,9 +75,11 @@ def detect_stars_auto(lum: np.ndarray, noise: float, # attribute resolution on every call to compute_quality_metrics. try: from scipy.ndimage import gaussian_filter, laplace, maximum_filter + + from src.background import _gaussian_blur _SCIPY_AVAILABLE = True except ImportError: - laplace = maximum_filter = gaussian_filter = None + laplace = maximum_filter = gaussian_filter = _gaussian_blur = None _SCIPY_AVAILABLE = False @@ -105,7 +107,7 @@ def generate_star_mask(shape: Tuple[int, int], star_positions: Optional[object], point_mask = np.zeros(shape, dtype=np.float32) # np.maximum.at handles multiple stars landing on the same pixel safely. np.maximum.at(point_mask, (ys[in_bounds], xs[in_bounds]), 1.0) - blurred = gaussian_filter(point_mask, sigma=sigma) + blurred = _gaussian_blur(point_mask, sigma) max_val = blurred.max() if max_val > 0: mask = blurred / max_val # normalise to [0, 1] diff --git a/src/registration.py b/src/registration.py index 4cf088d..f47c51a 100644 --- a/src/registration.py +++ b/src/registration.py @@ -12,6 +12,7 @@ import scipy.fft as sfft from scipy import ndimage +from src.background import _gaussian_blur from src.gpu_context import get_gpu from src.models import Config, FrameInfo, ProcessingStats from src.utils import get_logger, safe_print @@ -1586,9 +1587,11 @@ def run_registration_phase( if ref_stars is None: try: # 3. Local-maxima fallback — pure scipy, always available - from scipy.ndimage import gaussian_filter, maximum_filter + from scipy.ndimage import maximum_filter + + from src.background import _gaussian_blur _redet_tried.append('local-maxima') - smoothed = gaussian_filter(ref_lum.astype(np.float64), sigma=2.0) + smoothed = _gaussian_blur(ref_lum.astype(np.float64), 2.0) bg = float(np.median(smoothed)) thresh = bg + 5.0 * max(noise_val, float(np.std(smoothed)) * 0.5) local_max = maximum_filter(smoothed, size=11) @@ -2196,8 +2199,8 @@ def _find_in_map(filtered: np.ndarray) -> Optional[Tuple[float, float]]: try: sigma_small = Config.COMET_DOG_SIGMA_SMALL sigma_large = Config.COMET_DOG_SIGMA_LARGE - dog = (ndimage.gaussian_filter(lum64, sigma=sigma_large) - - ndimage.gaussian_filter(lum64, sigma=sigma_small)) + dog = (_gaussian_blur(lum64, sigma_large) - + _gaussian_blur(lum64, sigma_small)) dog = np.clip(dog, 0.0, None) # keep only positive (bright blob) response if dog.max() > 0: result = _find_in_map(dog) @@ -2207,7 +2210,7 @@ def _find_in_map(filtered: np.ndarray) -> Optional[Tuple[float, float]]: pass # --- Fallback: original Gaussian-blur + threshold method --- - smoothed = ndimage.gaussian_filter(lum64, sigma=smooth_sigma) + smoothed = _gaussian_blur(lum64, smooth_sigma) result = _find_in_map(smoothed) if result is not None: return result @@ -2255,7 +2258,7 @@ def find_extended_source_ellipse( lum64 = lum.astype(np.float64) if smooth_sigma is None: smooth_sigma = max(20.0, min(H, W) / 50.0) - lum_smooth = ndimage.gaussian_filter(lum64, sigma=smooth_sigma) + lum_smooth = _gaussian_blur(lum64, smooth_sigma) border = max(4, int(min(H, W) * 0.05)) edge = np.concatenate([ @@ -2411,7 +2414,7 @@ def find_comet_tail_pa(lum: np.ndarray, nucleus_y: float, nucleus_x: float, try: H, W = lum.shape # Smooth strongly to get a clean gradient at the coma scale - smoothed = ndimage.gaussian_filter(lum.astype(np.float64), sigma=max(sample_radius * 0.3, 5.0)) + smoothed = _gaussian_blur(lum.astype(np.float64), max(sample_radius * 0.3, 5.0)) # Compute gradient at nucleus location (via sobel or finite differences on the smoothed image) gy = ndimage.sobel(smoothed, axis=0) gx = ndimage.sobel(smoothed, axis=1) diff --git a/src/star_removal.py b/src/star_removal.py index bd89d6f..41ea2c8 100644 --- a/src/star_removal.py +++ b/src/star_removal.py @@ -24,9 +24,12 @@ try: from scipy.ndimage import gaussian_filter + + from src.background import _gaussian_blur _HAS_SCIPY = True except Exception: # pragma: no cover gaussian_filter = None + _gaussian_blur = None _HAS_SCIPY = False try: @@ -133,11 +136,11 @@ def remove_stars(rgb: np.ndarray, sources, fwhm: float, fill_sigma = max(6.0, 1.5 * max_r) out = rgb.astype(np.float32, copy=True) keep = (~mask).astype(np.float32) - denom = gaussian_filter(keep, sigma=fill_sigma) + denom = _gaussian_blur(keep, fill_sigma) denom = np.maximum(denom, 1e-6) for c in range(3): ch = out[:, :, c] - bg = gaussian_filter(ch * keep, sigma=fill_sigma) / denom + bg = _gaussian_blur(ch * keep, fill_sigma) / denom ch[mask] = bg[mask] out[:, :, c] = ch diff --git a/src/trail_reject.py b/src/trail_reject.py index 264e3c2..0f4240b 100644 --- a/src/trail_reject.py +++ b/src/trail_reject.py @@ -34,9 +34,12 @@ try: from scipy.ndimage import binary_dilation, gaussian_filter + + from src.background import _gaussian_blur _HAS_SCIPY = True except Exception: # pragma: no cover binary_dilation = gaussian_filter = None + _gaussian_blur = None _HAS_SCIPY = False try: @@ -244,11 +247,11 @@ def reject_trails(rgb: np.ndarray, lum: Optional[np.ndarray] = None, keep = (~mask).astype(np.float32) # Normalised convolution: smooth background from the non-trail pixels only, # so the fill never smears trail flux back into the gap. - denom = gaussian_filter(keep, sigma=fill_sigma) + denom = _gaussian_blur(keep, fill_sigma) denom = np.maximum(denom, 1e-6) for c in range(3): ch = out[:, :, c] - bg = gaussian_filter(ch * keep, sigma=fill_sigma) / denom + bg = _gaussian_blur(ch * keep, fill_sigma) / denom ch[mask] = bg[mask] out[:, :, c] = ch if verbose: diff --git a/tests/test_native.py b/tests/test_native.py index 96209e3..1153911 100644 --- a/tests/test_native.py +++ b/tests/test_native.py @@ -591,6 +591,24 @@ def test_gaussian_blur_falls_back_for_non_2d(): np.testing.assert_array_equal(got, want) +def test_gaussian_blur_preserves_input_dtype(): + """The native kernel always computes in float64 (like every other kernel + in this file); _gaussian_blur must cast back down so a float32 caller + (several master-calibration and Phase 4 call sites pass float32) doesn't + silently get a float64 array back and double its memory footprint.""" + assert _background_mod._HAS_NATIVE_GAUSSIAN # sanity: this run actually has native available + rng = np.random.default_rng(9) + a32 = rng.normal(1000.0, 50.0, (60, 70)).astype(np.float32) + out32 = _background_mod._gaussian_blur(a32, 3.0) + assert out32.dtype == np.float32 + + a64 = a32.astype(np.float64) + out64 = _background_mod._gaussian_blur(a64, 3.0) + assert out64.dtype == np.float64 + # same blur either way, just float32-rounded + np.testing.assert_allclose(out32.astype(np.float64), out64, 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