diff --git a/.claude/settings.json b/.claude/settings.json index b15e19f..ff6952b 100644 --- a/.claude/settings.json +++ b/.claude/settings.json @@ -4,7 +4,10 @@ "PowerShell(cmd /c \"type C:\\\\Users\\\\bmant\\\\anaconda3\\\\envs\\\\astrolab\\\\qt6.conf\")", "PowerShell(cmd /c \"C:\\\\Users\\\\bmant\\\\anaconda3\\\\Scripts\\\\conda.exe run -n astrolab python -c \"\"import struct, sys; data = open\\(r'C:\\\\Users\\\\bmant\\\\anaconda3\\\\envs\\\\astrolab\\\\Lib\\\\site-packages\\\\PyQt6\\\\Qt6\\\\bin\\\\Qt6Core.dll', 'rb'\\).read\\(\\); print\\('DLL size:', len\\(data\\)\\)\"\" 2>&1\")", "PowerShell(cmd /c \"C:\\\\Users\\\\bmant\\\\anaconda3\\\\Scripts\\\\conda.exe run -n astrolab pip install pyqt6==6.11.0 pyqt6-qt6==6.11.0 pyqt6-sip==13.11.1 --force-reinstall --ignore-installed 2>&1\")", - "WebSearch" + "WebSearch", + "PowerShell(conda activate astrolab)", + "PowerShell(python -m pytest tests/ -m \"not slow\" -q 2>&1)", + "PowerShell(Get-ChildItem -Path \"$env:USERPROFILE\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue -Recurse -Depth 2 | Select-Object -First 5 FullName; Get-ChildItem -Path \"C:\\\\ProgramData\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue | Select-Object FullName)" ] } } diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 34b5ceb..977da0d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,6 +13,7 @@ jobs: os: [windows-latest, ubuntu-latest, macos-latest] runs-on: ${{ matrix.os }} + timeout-minutes: 20 steps: - uses: actions/checkout@v4 @@ -44,6 +45,7 @@ jobs: os: [windows-latest, ubuntu-latest, macos-latest] runs-on: ${{ matrix.os }} + timeout-minutes: 45 steps: - uses: actions/checkout@v4 @@ -57,6 +59,7 @@ jobs: # to import PyQt6 during the analysis phase. - name: Install Linux Qt system libraries if: runner.os == 'Linux' + timeout-minutes: 5 run: | sudo apt-get update -qq sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index be16023..589c54a 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -16,6 +16,7 @@ jobs: os: [windows-latest, ubuntu-latest, macos-latest] runs-on: ${{ matrix.os }} + timeout-minutes: 45 steps: - uses: actions/checkout@v4 @@ -27,6 +28,7 @@ jobs: - name: Install Linux Qt system libraries if: runner.os == 'Linux' + timeout-minutes: 5 run: | sudo apt-get update -qq sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0 diff --git a/AstroImageLab.py b/AstroImageLab.py index a5117b8..85ad46e 100644 --- a/AstroImageLab.py +++ b/AstroImageLab.py @@ -4,8 +4,8 @@ # # PR → merge to main (CI runs tests + build to verify everything works) # Tag the merge commit on main → triggers the release workflow -# git tag v0.0.7 -# git push origin v0.0.7 +# git tag v0.0.8 +# git push origin v0.0.8 import sys import os diff --git a/CLAUDE.md b/CLAUDE.md index 8f74ecf..81743bf 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -59,6 +59,11 @@ synthetic/ | `_extract_cutout(data, xc, yc, radius)` | `gui/halo_dialog.py` | 2r×2r patch centred on star, zero-padded at image edges | | `_annular_rdf(log_data, xc, yc, radius)` | `gui/halo_dialog.py` | 1-px annular mean/std in log10 space; mirrors `HaloAnalyzer._annular_stats` | | `_power_ratio_db(freq_a, rp_a, freq_b, rp_b)` | `report_builder.py` | 10·log10 dB ratio between two radial power curves; returns `None` on missing data or misaligned frequency bins | +| `_log_ratio_map(a, b)` | `analysis/image_filters.py` | Per-pixel `log10(\|A\|/\|B\|)` map with percentile-based epsilon floor and defensive shape crop — the Section 8 replacement for plain `A − B` diff | +| `_log_ratio_color_range(diff)` | `analysis/image_filters.py` | Symmetric `(vmin, vmax)` for the `bwr` log-ratio colormap, shared by the log-ratio map panel, its histogram, and the correlation scatter dot coloring | +| `_plot_mask_illustration(base, mask_neb, mask_bg)` | `analysis/image_filters.py` | Translucent steelblue/tomato mask overlay on a grayscale base image | +| `_plot_metric_correlation(map_a, map_b, log_ratio, mask_neb, mask_bg, ...)` | `analysis/image_filters.py` | 1×2 masked-region scatter (A vs B) with a 1:1 line; each point colored by its pixel's log-ratio value using the same `bwr` scale as the adjacent map figure | +| `_family_figs_with_corr(rows, map_key_fn)` | `report_builder.py` | Emits a Section 8 family's map figure immediately followed by its `corr_*` correlation scatter, one scale at a time, in numeric order (`_SPATIAL_CORR_ROWS`) — the pattern to follow when adding any new per-scale Section 8 figure pair | --- @@ -226,11 +231,13 @@ When adding a new A-vs-B ratio curve to a report figure (precedent: `_power_rati `20 * np.log10(ratio)`. Using the wrong constant is silently off by 2× in dB — no exception, no obviously-wrong output, just a subtly incorrect number. - **Don't add the ratio via `ax.twinx()`** onto the existing absolute-value plot unless - both axes are the same kind of quantity (linear-vs-linear, as in - `analysis/image_filters.py::_plot_cross_section`'s A−B difference line). A linear, - zero-centered ratio next to a log-scale absolute axis has no principled vertical - alignment between the two scales — matplotlib's independent autoscaling invents a - relationship that isn't in the data. Build a separate, dedicated figure/panel instead. + both axes are the same kind of quantity (linear-vs-linear). A linear, zero-centered + ratio next to a log-scale absolute axis has no principled vertical alignment between + the two scales — matplotlib's independent autoscaling invents a relationship that + isn't in the data. Build a separate, dedicated figure/panel instead. (The codebase's + prior linear-vs-linear precedent, `_draw_cross_section`'s A−B difference line, was + removed as unnecessary clutter — there is currently no `ax.twinx()` usage anywhere + in the codebase, so treat this as a rule to apply fresh, not an existing pattern to copy.) - **Guard bin alignment before dividing two arrays from different analyses.** Two per-image radial/frequency arrays are only safely divisible bin-for-bin when they share the same shape *and* values (`freq_a.shape == freq_b.shape and @@ -239,6 +246,37 @@ When adding a new A-vs-B ratio curve to a report figure (precedent: `_power_rati guaranteed whenever an auto-selected ROI is involved (`_extract_roi` in `analysis/power_spectrum.py` computes `N` independently per image when no explicit ROI is set). Degrade gracefully — return `None` / skip the curve — rather than crash. +- **Epsilon-flooring a per-pixel ratio map needs a percentile, not a raw minimum.** + `_power_ratio_db`'s `positive.min() * 0.01` floor is fine for small 1-D arrays + (frequency bins, ~10²–10³ samples) but fragile at per-pixel map scale (10⁵–10⁷ + samples): the minimum order statistic over that many samples can be pathologically + tiny and let one spurious pixel dominate the ratio's dynamic range. `_log_ratio_map` + (`analysis/image_filters.py`) instead floors both operands at a low percentile + (`SECTION8_LOGRATIO_EPS_PERCENTILE`, default 1st) of the pooled positive `|A|,|B|` + values — same tool as the existing display-clipping precedent + (`_plot_side_by_side`'s `np.percentile(arr, 0.5)`), applied to the epsilon floor + instead of just the color scale. + +### Background estimation — compute once via the pre-pass, never redundantly + +`AstroImage.estimate_background()` (`core/astro_image.py`) is idempotent: it returns +immediately if `self.background is not None`, since `self.data` is only ever set once, +during `load()`. Every analyzer (`SNRAnalyzer`, `PSFAnalyzer`, `HaloAnalyzer`, +`EdgeAnalyzer`, `PowerSpectrumAnalyzer`, `SpatialDetailAnalyzer`) still calls +`estimate_background()` unconditionally at the top of its `analyze()` — that's +intentional and does not need to change; the idempotency guard just makes each of +those calls a cheap no-op once the object's background has already been computed. + +`gui/analysis_thread.py::_execute()` runs a pre-pass — after alignment, before task +dispatch — that calls `estimate_background()` once per distinct `AstroImage` object +(`img_a`, `img_b`, `self._starless_a`, `self._starless_b`) via a small +`ThreadPoolExecutor`. This exists because multiple analyzers share the same image +object and can run concurrently under `parallel=True`; without the pre-pass each one +would independently trigger a full `Background2D` computation (expensive) and race to +write `self.background` / `self.background_rms` on the same object. When adding a new +analyzer that needs background stats, just call `image.estimate_background()` as +normal at the top of `analyze()` — do not add another pre-pass call site; the existing +one in `_execute()` already covers every image object the thread constructs. --- @@ -337,6 +375,9 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | Adding a float64 cast in analysis code | Don't. All image data is float32 after `AstroImage.load()`. The only float64 exception is `astroalign` in `gui/analysis_thread.py`. Redundant float64 casts waste memory and defeat the float32 performance gains. | | Mixed float32/float64 arithmetic silently widens to float64 | NumPy upcasts when operands differ (e.g. `float32_array - float64_scalar`). If photutils ever returns a float64 background model, `background_subtracted()` will silently return float64. Guard by adding `.astype(np.float32)` at the end of `background_subtracted()` in `astro_image.py` if this is observed. | | `_section_snr` crashed when SNR metric is unchecked | `_plot_snr_pair` (`report_builder.py`) built a `panels` list filtered to non-`None` entries but never checked whether it was empty before calling `plt.subplots(1, len(panels), ...)` — 0 columns raised `ValueError: Number of columns must be a positive integer, not 0`. Hit whenever SNR is unchecked while another metric (e.g. Power Spectrum) is run. Fixed with an early `if not panels: return None` guard, matching `_plot_radial_overlay`/`_plot_radial_ratio_db`; both call sites already pipe the result through `_img_tag`, which turns `None` into `""`. | +| New Section 8 panel key doesn't need Report Inspector code changes | `gui/report_inspector.py` is fully generic — driven entirely by a companion `_inspector.npz` (raw float32/uint8 arrays) plus an embedded `catalog_json` built in `report_builder.py::_write_inspector_file`. `_panel_display_name`/`_panel_concept` dynamically parse any `panels` dict key prefix, so a new `SpatialDetailAnalyzer` panel family auto-appears in the inspector with zero inspector-side changes. A genuinely new *visual type* is a different story: the inspector only knows how to `imshow` 2D/RGB arrays (side-by-side or slider-reveal), so scatter-style plots (Section 8's `corr_*` correlation figures, interleaved into 8b–8f right after each map figure via `_family_figs_with_corr`) must stay static-HTML-only unless new inspector canvas code is written. | +| Renumbering a Section 8 subsection misses caption cross-references | Section 8's sub-heading letters (8a–8g) are referenced by literal string in caption/info-box text scattered throughout `_section_spatial` — not just in the `

` tags (e.g. "see 8g for…", "(8b–8f, 8g)"). After adding, removing, or renumbering a subsection, `grep` the function for every old *and* new heading letter — HTML renders a stale cross-reference without error, it just silently misdirects the reader to the wrong subsection. | +| Stale ROI crashes Section 8 with "index -1 is out of bounds for axis 0 with size 0" | `MainWindow._roi` (`gui/main_window.py`) is never cleared when a new image is loaded into either panel. If the user draws an ROI on one image pair, then loads a smaller replacement pair without clearing it, the stale coordinates go out of bounds for the new image. NumPy doesn't raise on an out-of-range slice — `norm_a[ry0:ry1, rx0:rx1]` silently returns a zero-size array — so the crash surfaces much later and far from the real cause: `SpatialDetailAnalyzer._plot_mask_illustration → _stretch_for_display → np.percentile(empty_array, ...)`. The same unguarded `bgsub[y0:y1, x0:x1]` pattern exists in `power_spectrum.py::_extract_roi` and `edge_analyzer.py::analyze`, so a stale ROI can corrupt those sections too (with different, equally misleading errors) if they happen to run. Fixed at the single real boundary — `MainWindow._on_run()`, which is the only path that constructs `AnalysisThread` — by validating `self._roi` against every loaded image's `data.shape` right before `settings["roi"]` is set; an out-of-bounds ROI is cleared (falls back to auto-detect/full-image) with a `QMessageBox` explaining why, rather than patching each analyzer's slice individually. | --- diff --git a/analysis/image_filters.py b/analysis/image_filters.py index 9415f5f..2d52d87 100644 --- a/analysis/image_filters.py +++ b/analysis/image_filters.py @@ -10,6 +10,7 @@ import matplotlib.colors as mcolors matplotlib.use("Agg") import matplotlib.pyplot as plt +from matplotlib.patches import Patch from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter import pywt @@ -18,7 +19,8 @@ from core.models import (STD_KERNEL_SIZES, LOG_SIGMAS, WAVELET_NAME, WAVELET_LEVELS, WEBER_KERNEL_SIZES, XS_LINE_ALPHA, SECTION8_BORDER_CROP_FRACTION, SECTION8_ANALYSIS_CMAP, - XS_SNR_REGION_WIDTH) + XS_SNR_REGION_WIDTH, SECTION8_DIFF_DIST_MAX_SAMPLES, + SECTION8_LOGRATIO_EPS_PERCENTILE, SECTION8_SCATTER_MAX_SAMPLES) MAX_DIM_FOR_STD = 2048 # downsample to this before generic_filter (performance) _DISPLAY_SMOOTH_SIGMA = 1.0 # applied to maps before plotting; does NOT affect metrics @@ -58,6 +60,7 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, "weber_contrast_a": {}, "weber_contrast_b": {}, "panels": {}, + "diff_dist": {}, "nc_shared_nebula_pixels": 0, "std_nc_score_a": {}, "std_nc_score_b": {}, "std_nc_noise_a": {}, "std_nc_noise_b": {}, "std_nc_ratio": {}, @@ -120,26 +123,43 @@ def _clip01(v): return max(0.0, min(1.0, v)) result["display_roi"] = display_roi + # Shared nebula/background regions for noise-corrected A/B scoring and diff + # distributions: pixels BOTH images independently classify the same way. + # None in single-image mode. + mask_neb_shared = mask_bg_shared = None + if mask_neb_b is not None: + h_s = min(mask_neb_a.shape[0], mask_neb_b.shape[0]) + w_s = min(mask_neb_a.shape[1], mask_neb_b.shape[1]) + mask_neb_shared = mask_neb_a[:h_s, :w_s] & mask_neb_b[:h_s, :w_s] + mask_bg_shared = mask_bg_a[:h_s, :w_s] & mask_bg_b[:h_s, :w_s] + result["nc_shared_nebula_pixels"] = ( + int(np.count_nonzero(mask_neb_shared)) if mask_neb_shared is not None else 0 + ) + + # Fixed seed so the diff-distribution subsampling below is reproducible + # across report generations for the same input images. + diff_dist_rng = np.random.default_rng(42) + # Export the exact preprocessed array every std/LoG/wavelet/Weber/gradient # map is computed from (mean-normalised, ROI-cropped if a ROI was given) as # its own panel, so the Report Inspector can show the source image content # for a map alongside the map itself, pixel-aligned to the same crop. + original_diff = self._log_ratio_map(analysis_a, analysis_b) if analysis_b is not None else None result["panels"]["original"] = { "a": analysis_a.astype(np.float32), "b": analysis_b.astype(np.float32) if analysis_b is not None else None, - "diff": (analysis_a - analysis_b).astype(np.float32) if analysis_b is not None else None, + "diff": original_diff, } + if original_diff is not None: + result["diff_dist"]["original"] = self._diff_distribution( + original_diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) - # Shared nebula ROI for noise-corrected A/B scoring: pixels BOTH images - # independently classify as nebula. None in single-image mode. - mask_neb_shared = None - if mask_neb_b is not None: - h_s = min(mask_neb_a.shape[0], mask_neb_b.shape[0]) - w_s = min(mask_neb_a.shape[1], mask_neb_b.shape[1]) - mask_neb_shared = mask_neb_a[:h_s, :w_s] & mask_neb_b[:h_s, :w_s] - result["nc_shared_nebula_pixels"] = ( - int(np.count_nonzero(mask_neb_shared)) if mask_neb_shared is not None else 0 - ) + # Mask illustration: only meaningful in two-image mode, mirroring the + # violin plots' own empty-in-single-image-mode behaviour. + if mask_neb_shared is not None: + mask_fig = self._plot_mask_illustration( + result["panels"]["original"]["a"], mask_neb_shared, mask_bg_shared) + figures["mask_illustration"] = fig_to_b64(mask_fig, dpi=150) _label_b = image_b.label if image_b is not None else None @@ -156,6 +176,7 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, mask_neb_shared=mask_neb_shared, + mask_bg_shared=mask_bg_shared, diff_dist_rng=diff_dist_rng, ) _f_log = _ex.submit(self._log_analysis, analysis_a, analysis_b, log_sigmas, @@ -163,6 +184,7 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b, + mask_bg_shared=mask_bg_shared, diff_dist_rng=diff_dist_rng, ) _f_wav = _ex.submit(self._wavelet_analysis, analysis_a, analysis_b, wavelet, levels, @@ -170,13 +192,16 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b, + mask_bg_shared=mask_bg_shared, diff_dist_rng=diff_dist_rng, ) _f_web = _ex.submit(self._weber_analysis, analysis_a, analysis_b, weber_kernel_sizes, image_a.label, _label_b, display_roi=display_roi, + crosshair=crosshair_roi, mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b, + mask_bg_shared=mask_bg_shared, diff_dist_rng=diff_dist_rng, ) _f_grad = _ex.submit(self._gradient_analysis, analysis_a, analysis_b, log_sigmas, @@ -184,6 +209,7 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b, + mask_bg_shared=mask_bg_shared, diff_dist_rng=diff_dist_rng, ) std_b64, std_partial = _f_std.result() log_b64, log_partial = _f_log.result() @@ -209,6 +235,11 @@ def _clip01(v): return max(0.0, min(1.0, v)) result["panels"].update(wav_partial["panels"]) result["panels"].update(web_partial["panels"]) result["panels"].update(grad_partial["panels"]) + result["diff_dist"].update(std_partial["diff_dist"]) + result["diff_dist"].update(log_partial["diff_dist"]) + result["diff_dist"].update(wav_partial["diff_dist"]) + result["diff_dist"].update(web_partial["diff_dist"]) + result["diff_dist"].update(grad_partial["diff_dist"]) # Merge noise-corrected scores/noise-floors and compute A/B ratios centrally. for prefix, partial in (("std", std_partial), ("log", log_partial), @@ -337,13 +368,15 @@ def _std_analysis(self, norm_a, norm_b, kernel_sizes, label_a, label_b, display_roi=None, crosshair=None, - mask_neb_shared=None) -> tuple[dict, dict]: + mask_neb_shared=None, + mask_bg_shared=None, diff_dist_rng=None) -> tuple[dict, dict]: figures = {} partial: dict = { "contrast_ratios_a": {}, "contrast_ratios_b": {}, "std_nc_score_a": {}, "std_nc_score_b": {}, "std_nc_noise_a": {}, "std_nc_noise_b": {}, "panels": {}, + "diff_dist": {}, } single = norm_b is None for ks in kernel_sizes: @@ -366,11 +399,21 @@ def _std_analysis(self, norm_a, norm_b, partial["std_nc_score_b"][ks] = nc_b partial["std_nc_noise_b"][ks] = noise_b + diff = self._log_ratio_map(std_a, std_b) if not single else None partial["panels"][f"std_{ks}px"] = { "a": std_a.astype(np.float32), "b": std_b.astype(np.float32) if std_b is not None else None, - "diff": (std_a - std_b).astype(np.float32) if std_b is not None else None, + "diff": diff, } + if diff is not None: + partial["diff_dist"][f"std_{ks}px"] = self._diff_distribution( + diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) + if not single: + corr_fig = self._plot_metric_correlation( + std_a, std_b, diff, mask_neb_shared, mask_bg_shared, + label_a, label_b, f"Local σ (kernel {ks}px)", diff_dist_rng) + if corr_fig is not None: + figures[f"corr_std_{ks}px"] = corr_fig if not single and noise_a and noise_b: partial["panels"][f"nrm_std_{ks}px"] = { "a": (std_a / noise_a).astype(np.float32), @@ -378,16 +421,24 @@ def _std_analysis(self, norm_a, norm_b, "diff": None, } + xs_raw = None + if crosshair is not None and not single: + pos, pa = self._sample_line(std_a, **crosshair) + _, pb = self._sample_line(std_b, **crosshair) + xs_raw = (pos, pa, pb, label_a, label_b, + f"Cross-section — Local σ, kernel {ks}px") + if not single: fig = self._plot_side_by_side( self._crop_border(std_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(std_b, SECTION8_BORDER_CROP_FRACTION), f"Local σ — kernel {ks}px — {label_a}", f"Local σ — kernel {ks}px — {label_b}", - diff_title=f"Diff (A−B), kernel {ks}px", + diff_title=f"Log ratio (A/B), kernel {ks}px", cmap=SECTION8_ANALYSIS_CMAP, nonlinear_norm=True, display_roi=None, + xs_data=xs_raw, ) else: fig = self._plot_single( @@ -399,24 +450,23 @@ def _std_analysis(self, norm_a, norm_b, figures[f"std_{ks}px"] = fig if not single and noise_a and noise_b: + xs_nrm = None + if crosshair is not None: + pos_n, pa_n = self._sample_line(std_a / noise_a, **crosshair) + _, pb_n = self._sample_line(std_b / noise_b, **crosshair) + xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b, + f"Cross-section — Local σ (× noise floor), kernel {ks}px") figures[f"nrm_std_{ks}px"] = self._plot_side_by_side( self._crop_border(std_a / noise_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(std_b / noise_b, SECTION8_BORDER_CROP_FRACTION), f"Local σ (× noise floor) — kernel {ks}px — {label_a}", f"Local σ (× noise floor) — kernel {ks}px — {label_b}", - diff_title=f"Diff (A−B), noise-normalised, kernel {ks}px", + diff_title=f"Log ratio (A/B), noise-normalised, kernel {ks}px", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_nrm, ) - if crosshair is not None: - pos, pa = self._sample_line(std_a, **crosshair) - if not single: - _, pb = self._sample_line(std_b, **crosshair) - figures[f"xs_std_{ks}px"] = self._plot_cross_section( - pos, pa, pb, label_a, label_b, - f"Cross-section — Local σ, kernel {ks}px") - return figs_to_b64(figures, dpi=150), partial def _compute_std_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray: @@ -476,6 +526,75 @@ def _nc_score(self, detail_map: np.ndarray, return None, None return float(bn.median(neb_vals)) / noise_floor, noise_floor + @staticmethod + def _log_ratio_map(a: np.ndarray, b: np.ndarray) -> np.ndarray: + """Per-pixel log10(|A|/|B|) map, replacing plain A-B difference. + + Sign is discarded (via abs()) before the ratio so the result is always + well-defined even for map families that can go negative (original + background-subtracted flux, wavelet band-pass reconstructions) — this + turns the comparison into "which image shows more structure/contrast at + this scale", not "which image is signed-brighter". Non-negative families + (std, |LoG|, gradient magnitude, Weber contrast) are unaffected by the + abs() since they're already >= 0. + + Epsilon-floors both operands using a low percentile (not a raw minimum — + a raw minimum over a multi-megapixel array can be pathologically tiny and + let a single spurious pixel dominate the log-ratio's dynamic range) of the + pooled positive values from both inputs, mirroring report_builder.py's + _power_ratio_db epsilon pattern but adapted for per-pixel map sizes. + + A and B are defensively cropped to their common shape first — two-image + analysis can reach here with mismatched shapes when astroalign + registration fails and analysis proceeds unaligned. + """ + abs_a, abs_b = np.abs(a), np.abs(b) + h = min(abs_a.shape[0], abs_b.shape[0]) + w = min(abs_a.shape[1], abs_b.shape[1]) + abs_a, abs_b = abs_a[:h, :w], abs_b[:h, :w] + + positive = np.concatenate([abs_a[abs_a > 0].ravel(), abs_b[abs_b > 0].ravel()]) + eps = max(float(np.percentile(positive, SECTION8_LOGRATIO_EPS_PERCENTILE)), 1e-12) \ + if positive.size > 0 else 1e-12 + + ratio = np.clip(abs_a, eps, None) / np.clip(abs_b, eps, None) + return np.log10(ratio).astype(np.float32) + + @staticmethod + def _log_ratio_color_range(diff: np.ndarray) -> tuple[float, float]: + """Symmetric (vmin, vmax) for the bwr log-ratio colormap, shared by the + log-ratio map panel, its histogram, and the correlation scatter dots.""" + d_max = float(np.percentile(np.abs(diff), 99.5)) or 1.0 + return -d_max, d_max + + @staticmethod + def _diff_distribution(diff_map: np.ndarray, + mask_neb_shared: np.ndarray | None, + mask_bg_shared: np.ndarray | None, + rng: np.random.Generator) -> dict: + """Random-subsampled log10(|A|/|B|) ratio pixel populations for nebula vs + background. + + Returns {"nebula": ndarray, "background": ndarray} (float32, signed + log-ratio values — 0 means A=B — up to SECTION8_DIFF_DIST_MAX_SAMPLES + each). Either array is empty if the corresponding mask is unavailable + (single-image mode) or selects zero pixels. + """ + out = {"nebula": np.empty(0, dtype=np.float32), + "background": np.empty(0, dtype=np.float32)} + if mask_neb_shared is None or mask_bg_shared is None: + return out + h = min(diff_map.shape[0], mask_neb_shared.shape[0], mask_bg_shared.shape[0]) + w = min(diff_map.shape[1], mask_neb_shared.shape[1], mask_bg_shared.shape[1]) + cropped = diff_map[:h, :w] + for key, mask in (("nebula", mask_neb_shared), ("background", mask_bg_shared)): + vals = cropped[mask[:h, :w]] + if vals.size > SECTION8_DIFF_DIST_MAX_SAMPLES: + idx = rng.choice(vals.size, SECTION8_DIFF_DIST_MAX_SAMPLES, replace=False) + vals = vals[idx] + out[key] = vals.astype(np.float32) + return out + @staticmethod def _compute_nc_ratios(score_a: dict, score_b: dict) -> dict: """Per-scale A/B ratio of noise-corrected scores; {} if either side is empty @@ -496,12 +615,14 @@ def _log_analysis(self, norm_a, norm_b, sigmas, label_a, label_b, display_roi=None, crosshair=None, - mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None, + mask_bg_shared=None, diff_dist_rng=None) -> tuple[dict, dict]: figures = {} partial: dict = { "log_nc_score_a": {}, "log_nc_score_b": {}, "log_nc_noise_a": {}, "log_nc_noise_b": {}, "panels": {}, + "diff_dist": {}, } single = norm_b is None for sigma in sigmas: @@ -517,11 +638,21 @@ def _log_analysis(self, norm_a, norm_b, sigmas, partial["log_nc_score_b"][sigma] = nc_b partial["log_nc_noise_b"][sigma] = noise_b + diff = self._log_ratio_map(log_a, log_b) if log_b is not None else None partial["panels"][f"log_{sigma}"] = { "a": log_a.astype(np.float32), "b": log_b.astype(np.float32) if log_b is not None else None, - "diff": (log_a - log_b).astype(np.float32) if log_b is not None else None, + "diff": diff, } + if diff is not None: + partial["diff_dist"][f"log_{sigma}"] = self._diff_distribution( + diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) + if not single: + corr_fig = self._plot_metric_correlation( + log_a, log_b, diff, mask_neb_shared, mask_bg_shared, + label_a, label_b, f"|LoG| (σ={sigma}px)", diff_dist_rng) + if corr_fig is not None: + figures[f"corr_log_{sigma}"] = corr_fig if not single and noise_a and noise_b: partial["panels"][f"nrm_log_{sigma}"] = { "a": (log_a / noise_a).astype(np.float32), @@ -529,16 +660,24 @@ def _log_analysis(self, norm_a, norm_b, sigmas, "diff": None, } + xs_raw = None + if crosshair is not None and not single: + pos, pa = self._sample_line(log_a, **crosshair) + _, pb = self._sample_line(log_b, **crosshair) + xs_raw = (pos, pa, pb, label_a, label_b, + f"Cross-section — |LoG|, σ={sigma}px") + if not single: fig = self._plot_side_by_side( self._crop_border(log_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(log_b, SECTION8_BORDER_CROP_FRACTION), f"|LoG| σ={sigma}px — {label_a}", f"|LoG| σ={sigma}px — {label_b}", - diff_title=f"LoG diff (A−B), σ={sigma}px", + diff_title=f"|LoG| log-ratio (A/B), σ={sigma}px", cmap=SECTION8_ANALYSIS_CMAP, nonlinear_norm=True, display_roi=None, + xs_data=xs_raw, ) else: fig = self._plot_single( @@ -550,22 +689,22 @@ def _log_analysis(self, norm_a, norm_b, sigmas, figures[f"log_sigma{sigma}"] = fig if not single and noise_a and noise_b: + xs_nrm = None + if crosshair is not None: + pos_n, pa_n = self._sample_line(log_a / noise_a, **crosshair) + _, pb_n = self._sample_line(log_b / noise_b, **crosshair) + xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b, + f"Cross-section — |LoG| (× noise floor), σ={sigma}px") figures[f"nrm_log_{sigma}"] = self._plot_side_by_side( self._crop_border(log_a / noise_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(log_b / noise_b, SECTION8_BORDER_CROP_FRACTION), f"|LoG| (× noise floor) σ={sigma}px — {label_a}", f"|LoG| (× noise floor) σ={sigma}px — {label_b}", - diff_title=f"Diff (A−B), noise-normalised, σ={sigma}px", + diff_title=f"Log ratio (A/B), noise-normalised, σ={sigma}px", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_nrm, ) - - if crosshair is not None and not single: - pos, pa = self._sample_line(log_a, **crosshair) - _, pb = self._sample_line(log_b, **crosshair) - figures[f"xs_log_sigma{sigma}"] = self._plot_cross_section( - pos, pa, pb, label_a, label_b, - f"Cross-section — |LoG|, σ={sigma}px") return figs_to_b64(figures, dpi=150), partial # ------------------------------------------------------------------ @@ -576,7 +715,8 @@ def _gradient_analysis(self, norm_a, norm_b, sigmas, label_a, label_b, display_roi=None, crosshair=None, - mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None, + mask_bg_shared=None, diff_dist_rng=None) -> tuple[dict, dict]: """G = |gradient| at Gaussian scale sigma (first spatial derivative magnitude). Reuses the LOG_SIGMAS scale set so gradient and |LoG| are directly comparable at identical spatial scales. Structured identically to _log_analysis.""" @@ -585,6 +725,7 @@ def _gradient_analysis(self, norm_a, norm_b, sigmas, "gm_nc_score_a": {}, "gm_nc_score_b": {}, "gm_nc_noise_a": {}, "gm_nc_noise_b": {}, "panels": {}, + "diff_dist": {}, } single = norm_b is None for sigma in sigmas: @@ -600,11 +741,21 @@ def _gradient_analysis(self, norm_a, norm_b, sigmas, partial["gm_nc_score_b"][sigma] = nc_b partial["gm_nc_noise_b"][sigma] = noise_b + diff = self._log_ratio_map(gm_a, gm_b) if gm_b is not None else None partial["panels"][f"gradient_{sigma}"] = { "a": gm_a.astype(np.float32), "b": gm_b.astype(np.float32) if gm_b is not None else None, - "diff": (gm_a - gm_b).astype(np.float32) if gm_b is not None else None, + "diff": diff, } + if diff is not None: + partial["diff_dist"][f"gradient_{sigma}"] = self._diff_distribution( + diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) + if not single: + corr_fig = self._plot_metric_correlation( + gm_a, gm_b, diff, mask_neb_shared, mask_bg_shared, + label_a, label_b, f"Gradient |G| (σ={sigma}px)", diff_dist_rng) + if corr_fig is not None: + figures[f"corr_gradient_{sigma}"] = corr_fig if not single and noise_a and noise_b: partial["panels"][f"nrm_gradient_{sigma}"] = { "a": (gm_a / noise_a).astype(np.float32), @@ -612,16 +763,24 @@ def _gradient_analysis(self, norm_a, norm_b, sigmas, "diff": None, } + xs_raw = None + if crosshair is not None and not single: + pos, pa = self._sample_line(gm_a, **crosshair) + _, pb = self._sample_line(gm_b, **crosshair) + xs_raw = (pos, pa, pb, label_a, label_b, + f"Cross-section — Gradient, σ={sigma}px") + if not single: fig = self._plot_side_by_side( self._crop_border(gm_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(gm_b, SECTION8_BORDER_CROP_FRACTION), f"Gradient |G| σ={sigma}px — {label_a}", f"Gradient |G| σ={sigma}px — {label_b}", - diff_title=f"Gradient diff (A−B), σ={sigma}px", + diff_title=f"Gradient log-ratio (A/B), σ={sigma}px", cmap=SECTION8_ANALYSIS_CMAP, nonlinear_norm=True, display_roi=None, + xs_data=xs_raw, ) else: fig = self._plot_single( @@ -633,22 +792,22 @@ def _gradient_analysis(self, norm_a, norm_b, sigmas, figures[f"gradient_{sigma}"] = fig if not single and noise_a and noise_b: + xs_nrm = None + if crosshair is not None: + pos_n, pa_n = self._sample_line(gm_a / noise_a, **crosshair) + _, pb_n = self._sample_line(gm_b / noise_b, **crosshair) + xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b, + f"Cross-section — Gradient (× noise floor), σ={sigma}px") figures[f"nrm_gradient_{sigma}"] = self._plot_side_by_side( self._crop_border(gm_a / noise_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(gm_b / noise_b, SECTION8_BORDER_CROP_FRACTION), f"Gradient (× noise floor) σ={sigma}px — {label_a}", f"Gradient (× noise floor) σ={sigma}px — {label_b}", - diff_title=f"Diff (A−B), noise-normalised, σ={sigma}px", + diff_title=f"Log ratio (A/B), noise-normalised, σ={sigma}px", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_nrm, ) - - if crosshair is not None and not single: - pos, pa = self._sample_line(gm_a, **crosshair) - _, pb = self._sample_line(gm_b, **crosshair) - figures[f"xs_gradient_{sigma}"] = self._plot_cross_section( - pos, pa, pb, label_a, label_b, - f"Cross-section — Gradient, σ={sigma}px") return figs_to_b64(figures, dpi=150), partial # ------------------------------------------------------------------ @@ -659,7 +818,8 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, label_a, label_b, display_roi=None, crosshair=None, - mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None, + mask_bg_shared=None, diff_dist_rng=None) -> tuple[dict, dict]: figures = {} partial: dict = { "sigma_noise_a": None, "sigma_noise_b": None, @@ -667,6 +827,7 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, "wavelet_nc_score_a": {}, "wavelet_nc_score_b": {}, "wavelet_nc_noise_a": {}, "wavelet_nc_noise_b": {}, "panels": {}, + "diff_dist": {}, } single = norm_b is None @@ -712,26 +873,47 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, continue # display/panels only for levels 2-3, unchanged from prior behaviour display_level = human_level + diff = self._log_ratio_map(rec_a, rec_b) if rec_b is not None else None partial["panels"][f"wavelet_{display_level}"] = { "a": rec_a.astype(np.float32), "b": rec_b.astype(np.float32) if rec_b is not None else None, - "diff": (rec_a - rec_b).astype(np.float32) if rec_b is not None else None, + "diff": diff, } + if diff is not None: + partial["diff_dist"][f"wavelet_{display_level}"] = self._diff_distribution( + diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) + if not single: + # Raw signed reconstructions (not abs()) — complementary to the + # sign-discarding log-ratio map, shows whether band-pass detail + # flips sign between the two filters at a given pixel. + corr_fig = self._plot_metric_correlation( + rec_a, rec_b, diff, mask_neb_shared, mask_bg_shared, + label_a, label_b, f"Wavelet level {display_level}", diff_dist_rng) + if corr_fig is not None: + figures[f"corr_wavelet_{display_level}"] = corr_fig if not single and noise_a and noise_b: partial["panels"][f"nrm_wavelet_{display_level}"] = { "a": (rec_a / noise_a).astype(np.float32), "b": (rec_b / noise_b).astype(np.float32), "diff": None, } + xs_raw = None + if crosshair is not None and not single: + pos, pa = self._sample_line(rec_a, **crosshair) + _, pb = self._sample_line(rec_b, **crosshair) + xs_raw = (pos, pa, pb, label_a, label_b, + f"Cross-section — Wavelet level {display_level}") + if not single: fig = self._plot_side_by_side( self._crop_border(rec_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(rec_b, SECTION8_BORDER_CROP_FRACTION), f"Wavelet level {display_level} — {label_a}", f"Wavelet level {display_level} — {label_b}", - diff_title=f"Level {display_level} diff (A−B)", + diff_title=f"Level {display_level} log-ratio (A/B)", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_raw, ) else: fig = self._plot_single( @@ -742,23 +924,23 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, figures[f"wavelet_level{display_level}"] = fig if not single and noise_a and noise_b: + xs_nrm = None + if crosshair is not None: + pos_n, pa_n = self._sample_line(rec_a / noise_a, **crosshair) + _, pb_n = self._sample_line(rec_b / noise_b, **crosshair) + xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b, + f"Cross-section — Wavelet level {display_level} (× noise floor)") figures[f"nrm_wavelet_{display_level}"] = self._plot_side_by_side( self._crop_border(rec_a / noise_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(rec_b / noise_b, SECTION8_BORDER_CROP_FRACTION), f"Wavelet level {display_level} (× noise floor) — {label_a}", f"Wavelet level {display_level} (× noise floor) — {label_b}", - diff_title=f"Level {display_level} diff (A−B), noise-normalised", + diff_title=f"Level {display_level} log-ratio (A/B), noise-normalised", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_nrm, ) - if crosshair is not None and not single: - pos, pa = self._sample_line(rec_a, **crosshair) - _, pb = self._sample_line(rec_b, **crosshair) - figures[f"xs_wavelet_level{display_level}"] = self._plot_cross_section( - pos, pa, pb, label_a, label_b, - f"Cross-section — Wavelet level {display_level}") - return figs_to_b64(figures, dpi=150), partial def _estimate_noise(self, coeffs) -> float: @@ -800,7 +982,9 @@ def _reconstruct_level(self, coeffs, target_coeff_idx: int, def _weber_analysis(self, norm_a, norm_b, kernel_sizes, label_a, label_b, display_roi=None, - mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: + crosshair=None, + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None, + mask_bg_shared=None, diff_dist_rng=None) -> tuple[dict, dict]: figures = {} partial: dict = { "weber_contrast_a": {}, @@ -808,6 +992,7 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes, "weber_nc_score_a": {}, "weber_nc_score_b": {}, "weber_nc_noise_a": {}, "weber_nc_noise_b": {}, "panels": {}, + "diff_dist": {}, } single = norm_b is None for ks in kernel_sizes: @@ -828,11 +1013,21 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes, partial["weber_nc_score_b"][ks] = nc_b partial["weber_nc_noise_b"][ks] = noise_b + diff = self._log_ratio_map(wc_a, wc_b) if wc_b is not None else None partial["panels"][f"weber_{ks}px"] = { "a": wc_a.astype(np.float32), "b": wc_b.astype(np.float32) if wc_b is not None else None, - "diff": (wc_a - wc_b).astype(np.float32) if wc_b is not None else None, + "diff": diff, } + if diff is not None: + partial["diff_dist"][f"weber_{ks}px"] = self._diff_distribution( + diff, mask_neb_shared, mask_bg_shared, diff_dist_rng) + if not single: + corr_fig = self._plot_metric_correlation( + wc_a, wc_b, diff, mask_neb_shared, mask_bg_shared, + label_a, label_b, f"Weber contrast (kernel {ks}px)", diff_dist_rng) + if corr_fig is not None: + figures[f"corr_weber_{ks}px"] = corr_fig if not single and noise_a and noise_b: partial["panels"][f"nrm_weber_{ks}px"] = { "a": (wc_a / noise_a).astype(np.float32), @@ -840,16 +1035,24 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes, "diff": None, } + xs_raw = None + if crosshair is not None and not single: + pos, pa = self._sample_line(wc_a, **crosshair) + _, pb = self._sample_line(wc_b, **crosshair) + xs_raw = (pos, pa, pb, label_a, label_b, + f"Cross-section — Weber contrast, kernel {ks}px") + if not single: fig = self._plot_side_by_side( self._crop_border(wc_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(wc_b, SECTION8_BORDER_CROP_FRACTION), f"Weber contrast — kernel {ks}px — {label_a}", f"Weber contrast — kernel {ks}px — {label_b}", - diff_title=f"Weber diff (A−B), kernel {ks}px", + diff_title=f"Weber log-ratio (A/B), kernel {ks}px", cmap=SECTION8_ANALYSIS_CMAP, nonlinear_norm=True, # unbounded output; PowerNorm compresses dynamic range display_roi=None, # _crop_border already applied + xs_data=xs_raw, ) else: fig = self._plot_single( @@ -861,14 +1064,21 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes, figures[f"weber_{ks}px"] = fig if not single and noise_a and noise_b: + xs_nrm = None + if crosshair is not None: + pos_n, pa_n = self._sample_line(wc_a / noise_a, **crosshair) + _, pb_n = self._sample_line(wc_b / noise_b, **crosshair) + xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b, + f"Cross-section — Weber contrast (× noise floor), kernel {ks}px") figures[f"nrm_weber_{ks}px"] = self._plot_side_by_side( self._crop_border(wc_a / noise_a, SECTION8_BORDER_CROP_FRACTION), self._crop_border(wc_b / noise_b, SECTION8_BORDER_CROP_FRACTION), f"Weber (× noise floor) — kernel {ks}px — {label_a}", f"Weber (× noise floor) — kernel {ks}px — {label_b}", - diff_title=f"Diff (A−B), noise-normalised, kernel {ks}px", + diff_title=f"Log ratio (A/B), noise-normalised, kernel {ks}px", cmap=SECTION8_ANALYSIS_CMAP, display_roi=None, + xs_data=xs_nrm, ) return figs_to_b64(figures, dpi=150), partial @@ -909,7 +1119,10 @@ def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray, cmap: str = "viridis", nonlinear_norm: bool = False, display_roi=None, - smooth_display: bool = True) -> plt.Figure: + smooth_display: bool = True, + xs_data: tuple | None = None) -> plt.Figure: + """xs_data, if given, is (pos, prof_a, prof_b, label_a, label_b, xs_title) + for the embedded cross-section panel; None leaves that panel blank.""" # Crop to bright-feature ROI if available if display_roi is not None: r0, r1, c0, c1 = display_roi @@ -931,25 +1144,36 @@ def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray, # Sqrt (PowerNorm gamma=0.5) compresses bright stars, reveals faint nebula norm = mcolors.PowerNorm(gamma=0.5, vmin=vmin, vmax=vmax) if nonlinear_norm else None - # Difference panel (computed before any possible shape mismatch) - h_min = min(arr_a.shape[0], arr_b.shape[0]) - w_min = min(arr_a.shape[1], arr_b.shape[1]) - diff = arr_a[:h_min, :w_min] - arr_b[:h_min, :w_min] + # Log-ratio panel (helper handles any shape mismatch defensively) + diff = self._log_ratio_map(arr_a, arr_b) # Symmetric about zero so the "bwr" midpoint (white) always means no difference. - d_max = float(np.percentile(np.abs(diff), 99.5)) or 1.0 - dvmin, dvmax = -d_max, d_max - - # 3×1 column layout — A, B, then diff stacked vertically. - # Use 1:1 pixel aspect; size figure based on the cropped array dimensions. + dvmin, dvmax = self._log_ratio_color_range(diff) + + # Dark-mode-aware reference-line color (project convention). + is_dark = matplotlib.rcParams.get("figure.facecolor", "white") not in ("white", "#ffffff", 1.0) + orig_color = "white" if is_dark else "black" + + # 3-row grid: A|B on top, log-ratio diff | cross-section in the middle, + # a log-ratio pixel-value histogram spanning both columns on the bottom. + # A/B/diff share the source array's pixel aspect ratio (aspect="equal" + # imshow); the cross-section panel is a line plot with no such constraint + # but occupies an equal-size grid cell so the other three panels stay + # geometrically identical whether or not a crosshair (and thus xs_data) + # is set. h, w = arr_a.shape[:2] aspect_ratio = h / max(w, 1) - panel_w = 10.0 + panel_w = 5.0 # half the old single-column width — 2 columns now share it panel_h = panel_w * aspect_ratio - fig_h = panel_h * 3 + 2.0 # 3 panels + headroom for colorbars/titles - fig, axes = plt.subplots(3, 1, figsize=(panel_w, fig_h), - constrained_layout=True) - ax_a, ax_b, ax_diff = axes + hist_h = 2.2 + fig = plt.figure(figsize=(panel_w * 2, panel_h * 2 + hist_h + 1.5), + constrained_layout=True) + gs = fig.add_gridspec(3, 2, height_ratios=[panel_h, panel_h, hist_h]) + ax_a = fig.add_subplot(gs[0, 0]) + ax_b = fig.add_subplot(gs[0, 1]) + ax_diff = fig.add_subplot(gs[1, 0]) + ax_xs = fig.add_subplot(gs[1, 1]) + ax_hist = fig.add_subplot(gs[2, :]) for ax, arr, title in zip([ax_a, ax_b], [arr_a, arr_b], [title_a, title_b]): im = ax.imshow(arr, origin="upper", cmap=cmap, @@ -968,6 +1192,33 @@ def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray, ax_diff.axis("off") fig.colorbar(im_diff, ax=ax_diff, fraction=0.046, pad=0.04) + if xs_data is not None: + pos, prof_a, prof_b, xs_label_a, xs_label_b, xs_title = xs_data + self._draw_cross_section(ax_xs, pos, prof_a, prof_b, + xs_label_a, xs_label_b, xs_title) + else: + ax_xs.axis("off") + + # Histogram of the log-ratio map's pixel distribution, colored to match + # the diff panel above: full data range on the x-axis (no pixels hidden), + # but each bin's fill color is clipped to [dvmin, dvmax] so extreme-tail + # bins saturate to the same end colors imshow already uses for its own + # outliers. + counts, bin_edges, patches = ax_hist.hist(diff.ravel(), bins=60) + hist_norm = mcolors.Normalize(vmin=dvmin, vmax=dvmax, clip=True) + hist_cmap = plt.get_cmap("bwr") + bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + for patch, center in zip(patches, bin_centers): + patch.set_facecolor(hist_cmap(hist_norm(center))) + ax_hist.axvline(0.0, color=orig_color, linestyle="--", linewidth=1.0, label="A = B") + ax_hist.set_yscale("log") + ax_hist.set_xlabel("Log ratio, log10(|A|/|B|)", fontsize=8) + ax_hist.set_ylabel("Pixel count", fontsize=8) + ax_hist.tick_params(labelsize=7) + ax_hist.legend(fontsize=6.5, loc="upper right") + ax_hist.grid(True, alpha=0.3) + ax_hist.set_title("Log-ratio pixel distribution", fontsize=9) + return fig def _plot_single(self, arr_a: np.ndarray, title_a: str, @@ -996,6 +1247,48 @@ def _plot_single(self, arr_a: np.ndarray, title_a: str, fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) return fig + def _plot_mask_illustration(self, base: np.ndarray, mask_neb: np.ndarray, + mask_bg: np.ndarray, alpha: float = 0.45) -> plt.Figure: + """Image A shown with a translucent nebula/background mask overlay. + + Uses the same steelblue/tomato color convention as the Section 8a violin + plots (report_builder.py's palette = {"Nebula": "steelblue", + "Background": "tomato"}) so the two figures read as one visual language. + Unclassified pixels (neither mask) are left plain grayscale. mask_neb/ + mask_bg are expected to already be the two-image shared classification + (mask_neb_a & mask_neb_b) — the exact masks that feed the violin plots + and correlation scatter plots — so this figure depicts what those plots + are actually gated on. + """ + gray = self._stretch_for_display(base) + h = min(gray.shape[0], mask_neb.shape[0], mask_bg.shape[0]) + w = min(gray.shape[1], mask_neb.shape[1], mask_bg.shape[1]) + gray = gray[:h, :w] + mask_neb = mask_neb[:h, :w] + mask_bg = mask_bg[:h, :w] + + rgb = np.stack([gray, gray, gray], axis=-1) + neb_color = np.array(mcolors.to_rgb("steelblue")) + bg_color = np.array(mcolors.to_rgb("tomato")) + rgb[mask_neb] = (1 - alpha) * rgb[mask_neb] + alpha * neb_color + rgb[mask_bg] = (1 - alpha) * rgb[mask_bg] + alpha * bg_color + + aspect_ratio = h / max(w, 1) + panel_w = 8.0 + fig, ax = plt.subplots(figsize=(panel_w, panel_w * aspect_ratio + 1.0)) + ax.imshow(rgb, origin="upper", interpolation="nearest", aspect="equal") + ax.axis("off") + ax.set_title("Nebula / background mask regions (shared A∩B classification), " + "shown on Image A", fontsize=10) + legend_handles = [ + Patch(facecolor="steelblue", edgecolor="none", alpha=0.7, label="Nebula"), + Patch(facecolor="tomato", edgecolor="none", alpha=0.7, label="Background"), + Patch(facecolor="0.5", edgecolor="none", label="Unclassified"), + ] + ax.legend(handles=legend_handles, loc="lower right", fontsize=8, framealpha=0.8) + fig.tight_layout() + return fig + def _plot_snr_bars(self, snr_a: dict, snr_b: dict, label_a: str, label_b: str | None, levels: int) -> plt.Figure: fig, ax = plt.subplots(figsize=(7, 4)) @@ -1081,29 +1374,106 @@ def _sample_line(arr: np.ndarray, x0: float, y0: float, return positions, values @staticmethod - def _plot_cross_section(pos: np.ndarray, prof_a: np.ndarray, prof_b: np.ndarray, - label_a: str, label_b: str, title: str) -> plt.Figure: + def _draw_cross_section(ax, pos: np.ndarray, prof_a: np.ndarray, prof_b: np.ndarray, + label_a: str, label_b: str, title: str) -> None: + """Draw a cross-section profile into an existing Axes: the bottom-right + quadrant of _plot_side_by_side's 2×2 grid. Fonts/linewidths are tuned down + from the metric's old standalone-figure sizes since this panel is now + roughly a quarter the area.""" # Images with slightly different pixel dimensions produce different-length profiles n = min(len(pos), len(prof_a), len(prof_b)) pos, prof_a, prof_b = pos[:n], prof_a[:n], prof_b[:n] - fig, ax1 = plt.subplots(figsize=(9, 4), constrained_layout=True) - ax1.plot(pos, prof_a, color="steelblue", linewidth=1.5, alpha=XS_LINE_ALPHA, label=label_a) - ax1.plot(pos, prof_b, color="tomato", linewidth=1.5, alpha=XS_LINE_ALPHA, label=label_b) - ax1.set_xlabel("Position along line (px)") - ax1.set_ylabel("Map value") - ax1.legend(loc="upper left", fontsize=8) - ax1.grid(True, alpha=0.3) - - ax2 = ax1.twinx() - n = min(len(prof_a), len(prof_b)) - diff = prof_a[:n] - prof_b[:n] - ax2.plot(pos[:n], diff, color="#2ca02c", linewidth=1.2, - linestyle="--", alpha=0.85, label="A−B") - ax2.axhline(0, color="#2ca02c", linewidth=0.8, alpha=0.3) # zero-crossing reference - ax2.set_ylabel("Difference (A−B)", color="#2ca02c") - ax2.tick_params(axis="y", labelcolor="#2ca02c") - ax2.legend(loc="upper right", fontsize=8) - ax1.set_title(title, fontsize=10) + ax.plot(pos, prof_a, color="steelblue", linewidth=1.1, alpha=XS_LINE_ALPHA, label=label_a) + ax.plot(pos, prof_b, color="tomato", linewidth=1.1, alpha=XS_LINE_ALPHA, label=label_b) + ax.set_xlabel("Position along line (px)", fontsize=8) + ax.set_ylabel("Map value", fontsize=8) + ax.tick_params(labelsize=7) + ax.legend(loc="upper left", fontsize=6.5, labelspacing=0.3) + ax.grid(True, alpha=0.3) + ax.set_title(title, fontsize=9) + + @staticmethod + def _plot_metric_correlation(map_a: np.ndarray, map_b: np.ndarray, + log_ratio: np.ndarray, + mask_neb_shared: np.ndarray, mask_bg_shared: np.ndarray, + label_a: str, label_b: str, metric_title: str, + rng: np.random.Generator, + max_samples: int = SECTION8_SCATTER_MAX_SAMPLES) -> plt.Figure | None: + """1x2 correlation scatter (Nebula | Background): metric value in A (y) + vs. metric value in B (x), with a dashed 1:1 reference line. Each point + is colored by that pixel's log-ratio value (log_ratio — the same array + driving the adjacent log-ratio map panel), using the same bwr colormap + and range, so the scatter visually links back to the map. + + Axis limits reflect the FULL pooled population range per subplot (not + percentile-clipped like the Section 8a violin plots) so upper-tail + divergence from the 1:1 line — the signal this plot exists to surface — + stays visible. Point clouds are randomly subsampled (up to max_samples) + for render cost only; the axis range is always computed from the full, + unsampled population. Returns None if both subplots have too few points. + """ + h = min(map_a.shape[0], map_b.shape[0], log_ratio.shape[0], + mask_neb_shared.shape[0], mask_bg_shared.shape[0]) + w = min(map_a.shape[1], map_b.shape[1], log_ratio.shape[1], + mask_neb_shared.shape[1], mask_bg_shared.shape[1]) + map_a = map_a[:h, :w] + map_b = map_b[:h, :w] + log_ratio = log_ratio[:h, :w] + mask_neb_shared = mask_neb_shared[:h, :w] + mask_bg_shared = mask_bg_shared[:h, :w] + + # Dark-mode-aware reference-line color (project convention). + is_dark = matplotlib.rcParams.get("figure.facecolor", "white") not in ("white", "#ffffff", 1.0) + orig_color = "white" if is_dark else "black" + + # Shared across both panels so a given color always means the same + # log-ratio value, matching the adjacent map figure's own scale. + dvmin, dvmax = SpatialDetailAnalyzer._log_ratio_color_range(log_ratio) + + fig, axes = plt.subplots(1, 2, figsize=(9, 4.5)) + any_data = False + for ax, region_name, mask in zip(axes, ("Nebula", "Background"), + (mask_neb_shared, mask_bg_shared)): + a_vals = map_a[mask] + b_vals = map_b[mask] + c_vals = log_ratio[mask] + n = a_vals.size + if n < 3: + ax.set_visible(False) + continue + any_data = True + + lo = float(min(a_vals.min(), b_vals.min())) + hi = float(max(a_vals.max(), b_vals.max())) + pad = 0.1 * (hi - lo) if hi > lo else 1.0 + lo -= pad + hi += pad + + if n > max_samples: + idx = rng.choice(n, max_samples, replace=False) + a_plot, b_plot, c_plot = a_vals[idx], b_vals[idx], c_vals[idx] + else: + a_plot, b_plot, c_plot = a_vals, b_vals, c_vals + + sc = ax.scatter(b_plot, a_plot, c=c_plot, cmap="bwr", vmin=dvmin, vmax=dvmax, + alpha=0.55, s=8, zorder=3, edgecolors="none", rasterized=True) + ax.plot([lo, hi], [lo, hi], color=orig_color, linestyle="--", + linewidth=1.2, zorder=4, label="Slope = 1 (A = B)") + fig.colorbar(sc, ax=ax, fraction=0.046, pad=0.04, label="log10(|A|/|B|)") + ax.set_xlim(lo, hi) + ax.set_ylim(lo, hi) + ax.set_aspect("equal") + ax.set_xlabel(f"{metric_title} — {label_b} (x)", fontsize=8) + ax.set_ylabel(f"{metric_title} — {label_a} (y)", fontsize=8) + ax.set_title(f"{region_name} (n={n})", fontsize=9) + ax.tick_params(labelsize=7) + ax.legend(fontsize=6.5, loc="best", labelspacing=0.3) + ax.grid(True, alpha=0.3) + + if not any_data: + plt.close(fig) + return None + fig.tight_layout() return fig @staticmethod diff --git a/core/astro_image.py b/core/astro_image.py index add122c..e293209 100644 --- a/core/astro_image.py +++ b/core/astro_image.py @@ -272,6 +272,8 @@ def _extract_metadata(self) -> None: def estimate_background(self, box_size: int = 64) -> None: if self.data is None: raise RuntimeError("Image not loaded") + if self.background is not None: + return # already computed for this instance's data; self.data never changes post-load with warnings.catch_warnings(): warnings.simplefilter("ignore") self.background = Background2D( diff --git a/core/models.py b/core/models.py index 8c03695..f413f67 100644 --- a/core/models.py +++ b/core/models.py @@ -40,6 +40,9 @@ EPSF_MAX_STARS = 600 # maximum candidate stars passed to EPSFBuilder; limits computation time SECTION8_BORDER_CROP_FRACTION = 0.05 # fraction of each image dimension cropped from perimeter in Section 8 display maps SECTION8_ANALYSIS_CMAP = "viridis" # colormap for Section 8 A/B analysis map panels (std, LoG, wavelet) +SECTION8_DIFF_DIST_MAX_SAMPLES = 1000000 # per masked population, per scale — caps violin/KDE cost on full-res diff maps +SECTION8_LOGRATIO_EPS_PERCENTILE = 1.0 # percentile of pooled positive |A|,|B| values used as the epsilon floor in log10(|A|/|B|) +SECTION8_SCATTER_MAX_SAMPLES = 50000 # per masked population, per scale — caps render cost of Section 8g correlation scatter plots PSF_SPATIAL_MAP_SIZE = 150 # px; long-axis resolution of FWHM / eccentricity spatial maps PSF_SPATIAL_MAP_SMOOTH_SIGMA = 5.0 # Gaussian smoothing sigma (px) applied to spatial maps before display diff --git a/gui/analysis_thread.py b/gui/analysis_thread.py index 7289ba1..c8f9c56 100644 --- a/gui/analysis_thread.py +++ b/gui/analysis_thread.py @@ -109,6 +109,20 @@ def _execute(self) -> None: else: self.progress.emit(2, "Single-image mode — skipping alignment…") + # Pre-compute background once per distinct image object, ahead of the parallel + # dispatch below. Every analyzer (SNR, PSF, Halo, Edge, PowerSpectrum, Spatial) + # calls estimate_background() unconditionally on whichever image it's given; + # with parallel=True several of them can share the same img_a/img_b/starless + # instance and race to recompute the same Background2D concurrently. Computing + # it here first makes every downstream call a no-op (see the idempotency guard + # in AstroImage.estimate_background) instead of a redundant recomputation. + _bg_images = {id(im): im for im in + (img_a, img_b, self._starless_a, self._starless_b) + if im is not None} + if _bg_images: + with concurrent.futures.ThreadPoolExecutor(max_workers=len(_bg_images)) as ex: + list(ex.map(lambda im: im.estimate_background(), _bg_images.values())) + # Build ordered task list: (metric_key, display_label, callable) # Each callable is a zero-arg function that writes into result_a / result_b. tasks: list[tuple[str, str, Callable[[], None]]] = [] diff --git a/gui/main_window.py b/gui/main_window.py index 7f653f9..5b5cf7f 100644 --- a/gui/main_window.py +++ b/gui/main_window.py @@ -296,6 +296,35 @@ def _on_run(self, settings: dict) -> None: self._control.set_run_enabled(True) return + # A previously-drawn ROI is never auto-cleared when a new image is loaded, so it + # can silently go out of bounds for a smaller replacement image. An out-of-range + # ROI slices to a zero-size array in every ROI-aware analyzer (Section 8, Power + # Spectrum, Edge Detection) with no exception at the slice itself — it fails much + # later with a confusing "index -1 is out of bounds for axis 0 with size 0" (or + # similar) deep inside a percentile/zoom/regression call. Validate once here so + # every analyzer downstream always receives either None or an ROI that fits. + if self._roi is not None: + rx0, ry0, rx1, ry1 = self._roi + roi_fits = True + for img in (img_a, img_b): + if img is None: + continue + h, w = img.data.shape[:2] + if not (0 <= rx0 < rx1 <= w and 0 <= ry0 < ry1 <= h): + roi_fits = False + break + if not roi_fits: + self._roi = None + self._control.set_roi(None) + QMessageBox.warning( + self, "ROI no longer fits loaded images", + "The previously selected ROI no longer fits within the currently " + "loaded image(s) and has been cleared. Analysis will proceed using " + "the full image (or an auto-detected region) instead.\n\n" + "Use \"Select ROI…\" to redraw a region if you want to restrict " + "analysis to a specific area." + ) + # Merge ROI and crosshair from window state into settings settings["roi"] = self._roi settings["crosshair"] = self._crosshair diff --git a/report/report_builder.py b/report/report_builder.py index 96bdd27..b6d9c91 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -428,6 +428,144 @@ def _draw_boxwhisker(ax, vals_list): return img_html, caption_html +# Display order and labels for _spatial_diff_distributions_figure. Keys must match +# analysis/image_filters.py's partial["diff_dist"] keys exactly (same keys used for +# partial["panels"], see SpatialDetailAnalyzer._std_analysis/_log_analysis/etc.). +_SPATIAL_DIFF_DIST_ROWS = [ + ("original", "Original (normalised image)"), + ("std_3px", "Local σ — 3 px"), + ("std_5px", "Local σ — 5 px"), + ("std_10px", "Local σ — 10 px"), + ("log_1.5", "|LoG| — σ=1.5 px"), + ("log_3.0", "|LoG| — σ=3 px"), + ("log_6.0", "|LoG| — σ=6 px"), + ("gradient_1.5", "Gradient |G| — σ=1.5 px"), + ("gradient_3.0", "Gradient |G| — σ=3 px"), + ("gradient_6.0", "Gradient |G| — σ=6 px"), + ("wavelet_2", "Wavelet — level 2"), + ("wavelet_3", "Wavelet — level 3"), + ("weber_3px", "Weber contrast — 3 px"), + ("weber_5px", "Weber contrast — 5 px"), + ("weber_9px", "Weber contrast — 9 px"), +] + +# Same order/labels as above minus "original" (no kernel scale) — used to order the +# per-scale correlation scatter plots interleaved into Section 8b-8f, each one placed +# immediately after its corresponding map figure. +_SPATIAL_CORR_ROWS = [(k, label) for k, label in _SPATIAL_DIFF_DIST_ROWS if k != "original"] + + +def _spatial_diff_distributions_figure(diff_dist: dict) -> tuple[str, str]: + """Combined figure: one row per Section 8 calculation, each row showing a + nebula-region violin and a background-region violin of the A-B diff pixel + values, with an IQR box-plot overlay (median magenta, IQR cyan — same styling + as the ePSF section's _psf_distributions_figure). Always violin+box, never + strip/swarm — diff populations here are always high-N. + + Returns (img_html, caption_html), or ("", "") if no row has enough data + (including single-image mode, where diff_dist is empty). + """ + import seaborn as sns + import pandas as pd + + rows = [(k, label) for k, label in _SPATIAL_DIFF_DIST_ROWS if k in diff_dist] + has_data = any( + diff_dist[k]["nebula"].size >= 3 and diff_dist[k]["background"].size >= 3 + for k, _ in rows + ) + if not rows or not has_data: + return "", "" + + fig, axes = plt.subplots(len(rows), 1, figsize=(7, 1.3 * len(rows) + 1)) + fig.subplots_adjust(hspace=0.65, left=0.22, right=0.97, top=0.97, bottom=0.04) + palette = {"Nebula": "steelblue", "Background": "tomato"} + order = ["Nebula", "Background"] + + def _draw_boxwhisker(ax, vals_list): + for i, vals in enumerate(vals_list): + ax.boxplot( + [vals], positions=[i], vert=False, + widths=0.45, zorder=5, + patch_artist=True, + manage_ticks=False, + boxprops=dict(facecolor="none", edgecolor="#00e5ff", linewidth=1.5, alpha=0.9), + medianprops=dict(color="magenta", linewidth=2.0, alpha=0.9), + whiskerprops=dict(color="#00e5ff", linewidth=1.5, alpha=0.9), + capprops=dict(color="#00e5ff", linewidth=1.5, alpha=0.9), + flierprops=dict(marker="", visible=False), + ) + + for ax, (key, title) in zip(np.atleast_1d(axes), rows): + neb = diff_dist[key]["nebula"] + bg = diff_dist[key]["background"] + if neb.size < 3 or bg.size < 3: + ax.set_visible(False) + continue + + combined = np.concatenate([neb, bg]) + df = pd.DataFrame({ + "value": combined, + "region": (["Nebula"] * neb.size) + (["Background"] * bg.size), + }) + sns.violinplot(data=df, x="value", y="region", order=order, + palette=palette, inner=None, linewidth=0.8, ax=ax) + _draw_boxwhisker(ax, [neb, bg]) + + # Some detail maps (e.g. Weber contrast, unbounded near dark-sky pixels — + # see 8e methodology) have rare extreme-outlier log-ratios that stretch the + # axis so far the IQR box becomes an invisible sliver. Clip the *view* to the + # 1st-99th percentile of this row's own pooled data; the boxplot/violin + # values themselves are unaffected, only what's visible is cropped. + lo, hi = np.percentile(combined, [1.0, 99.0]) + if hi > lo: + pad = 0.05 * (hi - lo) + ax.set_xlim(lo - pad, hi + pad) + + ax.axvline(0.0, color="red", linestyle="--", linewidth=1.0, alpha=0.8, zorder=3) + ax.set_title(title, fontsize=8, loc="left", pad=2) + ax.set_xlabel("", fontsize=7) + ax.set_ylabel("", fontsize=7) + ax.tick_params(axis="x", labelsize=7) + ax.tick_params(axis="y", labelsize=7, pad=1) + ax.spines[["top", "right"]].set_visible(False) + + img_html = _img_tag(fig, "spatial_diff_distributions") + + caption_html = ( + '

' + "Pixel-wise log₁₀(A / B) ratio distributions. " + "Each row is one Section 8 calculation, shown as a " + "violin plot (kernel density estimate of the log-ratio pixel values, " + "randomly subsampled for display) with an IQR box-plot overlay: " + "a cyan box spanning Q1–Q3, " + "a magenta centre line at the " + "median. Nebula = pixels both " + "images classify as nebula; Background " + "= pixels both images classify as background sky — the same shared masks " + "illustrated above and used for the noise-corrected scores in 8b–8f, here " + "shown as full distributions rather than a single median ratio. " + "A red dashed line marks zero " + "(A = B, equal). Units are log₁₀ — ±0.3 ≈ a " + "2× difference, ±1.0 ≈ a 10× difference. For the " + "Original and Wavelet rows specifically (the only two families that can go " + "negative), sign was discarded before the ratio — these rows compare the " + "magnitude of structure, not signed brightness. " + "How to read it: a Background row centred at zero with a narrow IQR is " + "the expected noise floor; a Nebula row with a similarly narrow, zero-centred " + "distribution means the two filters agree at that scale. A Nebula median shifted " + "away from zero, or an IQR visibly wider than the Background row's, indicates a " + "real structural difference between the filters at that scale rather than noise. " + "Each row's x-axis is independently clipped to its own 1st–99th percentile " + "range so the IQR box stays visible; rows with rare extreme-outlier log-ratios (e.g. " + "Weber contrast near dark-sky pixels, see 8e) may have a small fraction of the " + "violin's tail extend beyond the visible axis — see the per-pixel correlation " + "scatter next to each map figure below (8b–8f) for the full, unclipped " + "upper-tail behaviour." + "

" + ) + return img_html, caption_html + + def _epsf_stars_cell(psf_metrics: dict) -> str: """Format the number of stars used to build the ePSF.""" n = psf_metrics.get("epsf_n_stars") @@ -541,7 +679,7 @@ def _nc_ratio_rows(score_a: dict, score_b: dict, ratio: dict, scale_label, val_f " Contrast (Weber's law, ΔL/L)" ' Local range vs. local median background — the most literal "Contrast" metric in ' ' this section' - ' Noise-corrected (NC) score (8b–8f, 8h)same scale as parent method' + ' Noise-corrected (NC) score (8b–8f, 8g)same scale as parent method' ' SNR of Detail' " Whether a detail-map response in the nebula is real structure or just this image's " ' own noise floor at that scale' @@ -972,7 +1110,7 @@ def _add_options_entry(section: str, name: str, options: dict, if pa_panel.get("diff") is not None: npz_diff = f"sp_{pkey}_diff" _add(npz_diff, pa_panel["diff"]) - sp_opts["Diff (A−B)"] = npz_diff + sp_opts["Log Ratio (A/B)"] = npz_diff _add_options_entry("Spatial Detail", img_set_name, sp_opts, concept=_panel_concept(pkey)) @@ -3873,6 +4011,34 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: sm = ra.spatial_metrics or {} figs = sm.get("figures", {}) + mask_fig = figs.get("mask_illustration") + dist_img, dist_caption = _spatial_diff_distributions_figure(sm.get("diff_dist", {})) + mask_html = ( + "

Nebula / Background Mask Regions

" + + _hires_img_tag(mask_fig, "Mask illustration") + + '

' + "How the masks are detected. Each image is independently classified " + "per-pixel from its own background-subtracted flux and noise RMS: a pixel is " + "Nebula if its background-subtracted " + "value exceeds 2× the image's RMS, and " + "Background if it's below 0.5× the " + "RMS (a top-5% percentile fallback applies if no pixel clears the nebula threshold). " + "Pixels between the two thresholds are left unclassified (plain grayscale). " + "How they're used below. The log-ratio distributions and the per-scale " + "correlation plots embedded in 8b–8f use the two-image intersection " + "of these masks — " + "a pixel counts as Nebula only if both images classify it as Nebula (same for " + "Background). This is deliberately conservative: it excludes pixels where one " + f"image's registration, PSF, or local noise disagrees with the other's. Shown on " + f"{ra.label}, the same array as the "Original" row below." + "

" + if mask_fig else "" + ) + dist_html = ( + "

8a. Log-Ratio Distribution & Mask Overview

" + mask_html + dist_img + dist_caption + if dist_img else "" + ) + cr_a = sm.get("contrast_ratios_a", {}) cr_b = sm.get("contrast_ratios_b", {}) @@ -3942,7 +4108,7 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: 'after accounting for each image\'s own noise level. Scale units differ by ' 'method (kernel px for std/Weber, Gaussian σ px for LoG/gradient, ≈2level ' 'px for wavelet) and ratios should not be compared numerically across methods — ' - 'see 8h for a cross-method overview. See also Section 7 for the frequency-domain ' + 'see 8g for a cross-method overview. See also Section 7 for the frequency-domain ' 'view of this same question. Maps below are also shown in noise-normalised form ' '(map ÷ noise floor) so that a shared colour scale is a fair visual comparison ' 'between A and B, even when their absolute noise levels differ.', @@ -3961,29 +4127,57 @@ def figs_for(prefix): out += _hires_img_tag(figs[key], key) + "\n" return out - def xs_figs_for(prefix: str) -> str: + def _family_figs_with_corr(rows, map_key_fn) -> str: + """Emit each row's raw map figure immediately followed by its + per-pixel correlation scatter (when present), in _SPATIAL_CORR_ROWS + (numeric-scale) order — not figs_for's alphabetic key sort, which + would put e.g. std_10px before std_3px.""" out = "" - for key in sorted(figs): - if key.startswith(prefix): - out += _hires_img_tag(figs[key], key) + "\n" + for key, _label in rows: + map_key = map_key_fn(key) + fig = figs.get(map_key) + if fig: + out += _hires_img_tag(fig, map_key) + "\n" + corr_fig = figs.get(f"corr_{key}") + if corr_fig: + out += _hires_img_tag(corr_fig, f"corr_{key}") + "\n" return out - def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: - """Emit each image immediately followed by its matching cross-section.""" - out = "" - for img_key in sorted(k for k in figs if k.startswith(img_prefix)): - suffix = img_key[len(img_prefix):] - out += _hires_img_tag(figs[img_key], img_key) + "\n" - xs_key = xs_prefix + suffix - if xs_key in figs: - out += _hires_img_tag(figs[xs_key], xs_key) + "\n" - return out + _std_rows = [(k, l) for k, l in _SPATIAL_CORR_ROWS if k.startswith("std_")] + _log_rows = [(k, l) for k, l in _SPATIAL_CORR_ROWS if k.startswith("log_")] + _gradient_rows = [(k, l) for k, l in _SPATIAL_CORR_ROWS if k.startswith("gradient_")] + _wavelet_rows = [(k, l) for k, l in _SPATIAL_CORR_ROWS if k.startswith("wavelet_")] + _weber_rows = [(k, l) for k, l in _SPATIAL_CORR_ROWS if k.startswith("weber_")] + + _has_any_corr = any(f"corr_{k}" in figs for k, _l in _SPATIAL_CORR_ROWS) + corr_methodology_box = ( + _info_box('Each per-pixel correlation scatter (embedded next to its map figure ' + 'below) plots the raw (pre-ratio) metric value at every shared pixel: ' + f'{ra.label} on the y-axis vs. {rb.label} on the x-axis, split into a ' + 'Nebula panel and a Background panel using the shared masks explained in 8a. ' + 'The black dashed diagonal is the 1:1 line (perfect agreement, ' + 'A = B). Axis limits span the full data range for that panel — ' + 'unlike the 8a violin plots, they are not clipped to a percentile — so ' + 'the behaviour of the upper tail is always visible. ' + 'Each point is colored by that pixel\'s log-ratio value ' + '(log10(|A|/|B|)), using the same bwr colour scale as the ' + 'adjacent log-ratio map panel and its histogram — red points are pixels where A ' + 'is stronger, blue where B is stronger, white where the two agree — so the ' + 'scatter\'s spatial pattern (visible in the map) links directly to its position ' + 'in this plot. ' + 'Point clouds are randomly subsampled for rendering; the axis range and each ' + 'panel\'s point count (n) are always computed from the full, unsampled ' + 'population.', + title="Per-pixel correlation methodology") + if _has_any_corr else "" + ) has_crosshair = sm.get("crosshair") is not None xs_note = _info_box( - 'ℹ Cross-section profiles below are extracted along ' - 'the line selected in the viewer. Left axis: both images ' - '(steelblue = A, tomato = B). Right axis (green dashed): difference A−B.', + 'ℹ When a cross-section line is set in the viewer, its profile is embedded ' + 'as the middle-right panel of each map-pair figure below (8b–8f): Image A ' + '(top-left), Image B (top-right), log-ratio map (middle-left), cross-section ' + 'profile (middle-right, steelblue = A, tomato = B).', title="Cross-section profiles", open=True, ) if has_crosshair else "" @@ -4037,7 +4231,8 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: f'thin filaments), Level 3 ≈ 8 px (medium structures — emission knots, shell edges), ' f'Level 4 ≈ 16 px (broader features). A higher SNR at Level 2 indicates the filter ' f'preserves sub-arcsecond detail better; Level 3 reflects medium-scale structure. ' - f'Cross-section profiles show how detail amplitude varies spatially along the selected line.', + f'When a cross-section line is set, its profile (embedded in each map figure below) ' + f'shows how detail amplitude varies spatially along the selected line.', title="Wavelet decomposition", ) @@ -4050,11 +4245,23 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: {smooth_note} {_info_box('All maps below are computed on mean-signal-normalised data ' '(each image divided by its own mean signal), making them dimensionless and comparable ' - 'across different filter bandwidths. Images are shown side-by-side with a shared ' - 'colour scale; the third panel shows the difference A−B.', + 'across different filter bandwidths. Each figure has Image A (top-left) and Image B ' + '(top-right) sharing a colour scale; the middle-left panel shows ' + 'the per-pixel log-ratio log10(A / B): red (positive) means A is stronger at ' + 'that pixel, blue (negative) means B is stronger, white means the two are equal. ±0.3 ≈ ' + 'a 2× difference, ±1.0 ≈ a 10× difference. For the Original and Wavelet maps ' + '— the only two families whose values can be negative — the ratio uses |A| and |B|, so ' + 'it compares the magnitude of structure rather than signed brightness. The ' + 'middle-right panel shows the cross-section profile along the selected line when one is ' + 'set, otherwise left blank. A third row spans both columns with a histogram of the ' + 'log-ratio map\'s pixel values, coloured bin-by-bin with the same bwr scale and range as ' + 'the log-ratio panel above it.', title="Spatial detail maps overview")} +{dist_html} {nc_methodology_box} {nc_empty_note} +{corr_methodology_box} +{xs_note}

8b. Local Standard Deviation Maps

{_info_box('Measures how much pixel values vary within a neighbourhood. ' @@ -4065,8 +4272,8 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: 'window. Brighter regions contain more local variation — typically nebula filaments, ' 'star halos, or noise. A filter with higher std values in targeted emission regions ' 'preserves more structure; higher std in blank sky regions indicates more photon noise. ' - 'The cross-section profiles below each map pair show how local detail amplitude varies ' - 'along the selected line.', + 'When a cross-section line is set, its profile is embedded in the middle-right panel of ' + 'each map figure below, showing how local detail amplitude varies along the selected line.', title="Local standard deviation")} @@ -4076,10 +4283,12 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: {std_nc_rows}
Kernel size{ra.label}{rb.label}
Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B
-{xs_note}{paired_figs_for("std_", "xs_std_")} -

Side-by-side local σ maps at each kernel size (shared colour scale), -each followed by its cross-section profile. -The difference map (right) highlights where one filter preserves more local variation.

+{_family_figs_with_corr(_std_rows, lambda k: k)} +

Local σ maps at each kernel size (shared colour scale): Image A (top-left), +Image B (top-right), log-ratio map (middle-left) highlighting where one filter preserves more +local variation, and — when a cross-section line is set — its profile (middle-right), plus a +bottom-row histogram of the log-ratio map's pixel values (same colour scale). Each map +is immediately followed by its per-pixel correlation scatter (see methodology above).

Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.

{figs_for("nrm_std_")}

8c. Laplacian of Gaussian (LoG) Maps

@@ -4091,17 +4300,20 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: 'then computing the Laplacian (second spatial derivative), which peaks at intensity ' 'boundaries. |LoG| is shown so bright-to-dark and dark-to-bright edges are treated ' 'equally. Compare maps at each σ: a sharper or higher-contrast filter will show ' - 'brighter LoG response at small σ values. Cross-section profiles reveal subtle ' - 'differences in edge sharpness along the selected line.', + 'brighter LoG response at small σ values. When a cross-section line is set, its profile ' + 'reveals subtle differences in edge sharpness along the selected line.', title="Laplacian of Gaussian (LoG)")} {log_nc_rows}
Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B
-{paired_figs_for("log_", "xs_log_")} -

|LoG| maps at σ = 1.5, 3, and 6 px (shared colour scale per row), -each followed by its cross-section profile. -A filter preserving more fine detail shows brighter, more defined boundaries at small σ.

+{_family_figs_with_corr(_log_rows, lambda k: "log_sigma" + k.split("_", 1)[1])} +

|LoG| maps at σ = 1.5, 3, and 6 px (shared colour scale per figure): +Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a +cross-section line is set — its profile (middle-right), plus a bottom-row histogram of +the log-ratio map's pixel values (same colour scale). A filter preserving more fine +detail shows brighter, more defined boundaries at small σ. Each map is immediately +followed by its per-pixel correlation scatter (see methodology above).

Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.

{figs_for("nrm_log_")}

8d. Wavelet Decomposition

@@ -4120,10 +4332,14 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: {wavelet_nc_rows} -{paired_figs_for("wavelet_level", "xs_wavelet_level")} +{_family_figs_with_corr(_wavelet_rows, lambda k: "wavelet_level" + k.split("_", 1)[1])}

Reconstructed detail images at levels 2 and 3 (shared colour scale, -diverging colourmap), each followed by its cross-section profile. -The difference panel (right) shows where fine structure differs between the two filters.

+diverging colourmap): Image A (top-left), Image B (top-right), log-ratio panel (middle-left) +showing where fine structure differs between the two filters (sign discarded — |A|/|B| — +since wavelet reconstructions can be negative), and — when a cross-section line is set — +its profile (middle-right), plus a bottom-row histogram of the log-ratio map's pixel values +(same colour scale). Each map is immediately followed by its per-pixel correlation +scatter (see methodology above).

Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.

{figs_for("nrm_wavelet_")} @@ -4158,12 +4374,15 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B {weber_nc_rows} -{figs_for("weber_")} +{_family_figs_with_corr(_weber_rows, lambda k: k)}

Per-pixel Weber fraction contrast maps (c = ΔL / L, square-root colour -scale, viridis). Brighter regions have higher Weber contrast — the local intensity range -is large relative to the local median luminance. The difference panel (A−B) shows where -one image achieves greater relative contrast. High values over dark-sky regions are -expected; use a nebula ROI for meaningful filter comparison.

+scale, viridis): Image A (top-left), Image B (top-right), log-ratio panel (middle-left) +showing where one image achieves greater relative contrast, and — when a cross-section +line is set — its profile (middle-right), plus a bottom-row histogram of the log-ratio +map's pixel values (same colour scale). Brighter regions have higher Weber contrast — +the local intensity range is large relative to the local median luminance. High values +over dark-sky regions are expected; use a nebula ROI for meaningful filter comparison. +Each map is immediately followed by its per-pixel correlation scatter (see methodology above).

Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.

{figs_for("nrm_weber_")} @@ -4185,14 +4404,17 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B {gm_nc_rows} -{paired_figs_for("gradient_", "xs_gradient_")} -

Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale per row), -each followed by its cross-section profile. -A filter preserving sharper boundaries shows brighter, more defined gradient response.

+{_family_figs_with_corr(_gradient_rows, lambda k: k)} +

Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale per +figure): Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a +cross-section line is set — its profile (middle-right), plus a bottom-row histogram of the +log-ratio map's pixel values (same colour scale). A filter preserving sharper +boundaries shows brighter, more defined gradient response. Each map is immediately followed +by its per-pixel correlation scatter (see methodology above).

Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.

{figs_for("nrm_gradient_")} -

8h. Noise-Corrected Contrast — Cross-Method Overview

+

8g. Noise-Corrected Contrast — Cross-Method Overview

{_hires_img_tag(figs.get("nc_ratio_overview"), "NC ratio overview")}

Ratio A/B for every noise-corrected method plotted against its approximate spatial scale. Scale units differ by method (see 8b–8f methodology diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py index 2a59626..b0fde89 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -315,3 +315,120 @@ def test_valid_masks_return_ratio(self): score, noise = analyzer._nc_score(detail, mask_neb_shared, bg_mask) assert score == pytest.approx(5.0) assert noise == pytest.approx(2.0) + + +class TestLogRatioHelper: + """Direct unit tests of _log_ratio_map's epsilon-floor and sign-discard contract.""" + + def test_equal_inputs_give_zero(self): + a = np.full((10, 10), 5.0, dtype=np.float32) + b = np.full((10, 10), 5.0, dtype=np.float32) + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert np.allclose(result, 0.0, atol=1e-5) + + def test_double_ratio_gives_log10_two(self): + b = np.full((10, 10), 3.0, dtype=np.float32) + a = 2.0 * b + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert np.allclose(result, np.log10(2.0), atol=1e-5) + + def test_opposite_sign_equal_magnitude_gives_zero(self): + a = np.full((10, 10), -5.0, dtype=np.float32) + b = np.full((10, 10), 5.0, dtype=np.float32) + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert np.allclose(result, 0.0, atol=1e-5) + + def test_near_zero_pixel_amid_nonzero_data_stays_finite(self): + a = np.full((10, 10), 5.0, dtype=np.float32) + b = np.full((10, 10), 5.0, dtype=np.float32) + b[0, 0] = 0.0 + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert np.all(np.isfinite(result)) + + def test_mismatched_shapes_crop_to_common_size(self): + a = np.full((10, 10), 5.0, dtype=np.float32) + b = np.full((8, 9), 5.0, dtype=np.float32) + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert result.shape == (8, 9) + + def test_all_zero_inputs_stay_finite(self): + a = np.zeros((10, 10), dtype=np.float32) + b = np.zeros((10, 10), dtype=np.float32) + result = SpatialDetailAnalyzer._log_ratio_map(a, b) + assert np.all(np.isfinite(result)) + + +class TestMaskIllustrationFigure: + def test_absent_in_single_image_mode(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert "mask_illustration" not in result["figures"] + + def test_present_in_two_image_mode(self, nc_result): + assert "mask_illustration" in nc_result["figures"] + assert isinstance(nc_result["figures"]["mask_illustration"], str) + assert len(nc_result["figures"]["mask_illustration"]) > 0 + + +class TestCorrelationScatterFigures: + _CORR_KEYS = ( + [f"corr_std_{ks}px" for ks in STD_KERNEL_SIZES] + + [f"corr_log_{s}" for s in LOG_SIGMAS] + + [f"corr_gradient_{s}" for s in LOG_SIGMAS] + + ["corr_wavelet_2", "corr_wavelet_3"] + + [f"corr_weber_{ks}px" for ks in WEBER_KERNEL_SIZES] + ) + + @pytest.mark.parametrize("key", _CORR_KEYS) + def test_present_in_two_image_mode(self, nc_result, key): + assert key in nc_result["figures"] + assert isinstance(nc_result["figures"][key], str) + assert len(nc_result["figures"][key]) > 0 + + @pytest.mark.parametrize("key", _CORR_KEYS) + def test_absent_in_single_image_mode(self, astro_image_a, key): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert key not in result["figures"] + + +@pytest.fixture(scope="module") +def nc_result_with_crosshair(nc_image_pair) -> dict: + img_a, img_b = nc_image_pair + crosshair = {"x0": 0.1, "y0": 0.5, "x1": 0.9, "y1": 0.5} + return SpatialDetailAnalyzer().analyze(img_a, img_b, crosshair=crosshair) + + +class TestCrosshairEmbeddedCrossSections: + """Section 8 cross-sections are embedded panels inside the 2x2 map figures, + not separate figures — no standalone xs_* keys should ever appear, and all + 5 families (including Weber, newly crosshair-capable) must still produce + their usual figure keys with or without a crosshair, without crashing.""" + + def test_no_crash_with_crosshair(self, nc_result_with_crosshair): + assert isinstance(nc_result_with_crosshair, dict) + + def test_no_standalone_xs_keys_with_crosshair(self, nc_result_with_crosshair): + figs = nc_result_with_crosshair["figures"] + assert not any(k.startswith(("xs_std_", "xs_log_", "xs_wavelet_", + "xs_gradient_", "xs_weber_")) for k in figs) + + def test_no_standalone_xs_keys_without_crosshair(self, nc_result): + figs = nc_result["figures"] + assert not any(k.startswith(("xs_std_", "xs_log_", "xs_wavelet_", + "xs_gradient_", "xs_weber_")) for k in figs) + + @pytest.mark.parametrize("key_prefix,scales", [ + ("std_", STD_KERNEL_SIZES), ("weber_", WEBER_KERNEL_SIZES), + ]) + def test_family_figures_present_with_crosshair(self, nc_result_with_crosshair, + key_prefix, scales): + figs = nc_result_with_crosshair["figures"] + for scale in scales: + assert f"{key_prefix}{scale}px" in figs + + def test_weber_no_crash_without_crosshair(self, nc_result): + # Regression: weber previously had no crosshair param at all. + assert all(f"weber_{ks}px" in nc_result["figures"] for ks in WEBER_KERNEL_SIZES) + + def test_weber_nrm_figures_present_with_crosshair(self, nc_result_with_crosshair): + figs = nc_result_with_crosshair["figures"] + assert any(k.startswith("nrm_weber_") for k in figs)