From ff01c34989d0266e42dd279bfae1fa5bb8887582 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Thu, 2 Jul 2026 15:32:04 -0400 Subject: [PATCH 01/14] Inspector: fix slider alignment timing and add Subtract Image Median checkbox Before/after slider margins now synced via QTimer.singleShot(0) so the deferred call runs after draw_idle() and constrained_layout have settled, preventing the slider from spanning the full window width on mode switch. Adds a Subtract image median checkbox alongside Normalize to peak. When checked, each image's global median is subtracted from its cross-section profile before any normalization, removing the STF background pedestal and shifting profiles to a near-zero baseline. Y-axis label updates to reflect the active combination of both options. Co-Authored-By: Claude Sonnet 4.6 --- gui/report_inspector.py | 56 +++++++++++++++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/gui/report_inspector.py b/gui/report_inspector.py index 876733d..3adf2f4 100644 --- a/gui/report_inspector.py +++ b/gui/report_inspector.py @@ -30,10 +30,10 @@ import matplotlib.pyplot as plt from matplotlib.backends.backend_qtagg import FigureCanvasQTAgg -from PyQt6.QtCore import Qt, pyqtSignal +from PyQt6.QtCore import Qt, QTimer, pyqtSignal from PyQt6.QtWidgets import ( QMainWindow, QWidget, QVBoxLayout, QHBoxLayout, - QLabel, QComboBox, QSlider, QFrame, QSizePolicy, + QLabel, QComboBox, QSlider, QFrame, QSizePolicy, QCheckBox, ) @@ -536,13 +536,22 @@ def __init__(self, parent: QWidget | None = None): def update_profiles(self, dist_a: np.ndarray, prof_a: np.ndarray, dist_b: np.ndarray | None, prof_b: np.ndarray | None, label_a: str, label_b: str, - normalized: bool = True) -> None: + normalized: bool = True, + subtract_median: bool = False) -> None: self._ax.cla() self._ax.plot(dist_a, prof_a, color="steelblue", lw=1.2, label=label_a) if dist_b is not None and prof_b is not None: self._ax.plot(dist_b, prof_b, color="tomato", lw=1.2, label=label_b) self._ax.set_xlabel("Distance (px)", fontsize=8) - self._ax.set_ylabel("Normalized intensity" if normalized else "Value", fontsize=8) + if normalized and subtract_median: + ylabel = "Normalized, bg-subtracted" + elif normalized: + ylabel = "Normalized intensity" + elif subtract_median: + ylabel = "Background-subtracted" + else: + ylabel = "Value" + self._ax.set_ylabel(ylabel, fontsize=8) self._ax.legend(fontsize=7, loc="upper right") self._ax.grid(True, alpha=0.25) self._ax.tick_params(labelsize=7) @@ -686,7 +695,6 @@ def _build_ui(self) -> None: root.addWidget(self._profile_canvas) # ── Profile options row ───────────────────────────────────────── - from PyQt6.QtWidgets import QCheckBox prof_opts = QHBoxLayout() self._normalize_check = QCheckBox("Normalize to peak") self._normalize_check.setChecked(False) # default: absolute values @@ -695,6 +703,13 @@ def _build_ui(self) -> None: "When unchecked, raw array values are shown on a shared scale,\n" "preserving absolute differences between the two images.") prof_opts.addWidget(self._normalize_check) + self._subtract_median_check = QCheckBox("Subtract image median") + self._subtract_median_check.setChecked(False) + self._subtract_median_check.setToolTip( + "Subtract each image's global median from its cross-section profile.\n" + "Removes the background pedestal so profiles start near zero.\n" + "Applied before 'Normalize to peak' when both are checked.") + prof_opts.addWidget(self._subtract_median_check) prof_opts.addStretch() root.addLayout(prof_opts) @@ -711,6 +726,7 @@ def _build_ui(self) -> None: # with the axes bounds (which shift slightly on resize or content change). self._image_canvas._canvas.mpl_connect("draw_event", self._on_canvas_drawn) self._normalize_check.stateChanged.connect(self._on_normalize_changed) + self._subtract_median_check.stateChanged.connect(self._on_normalize_changed) self._reset_zoom_btn.clicked.connect(self._image_canvas.reset_zoom) # ------------------------------------------------------------------ @@ -735,7 +751,7 @@ def _on_mode_changed(self, idx: int) -> None: else: self._slider_row_widget.setVisible(True) self._image_canvas.set_mode("slider") - self._sync_slider_margins() + QTimer.singleShot(0, self._sync_slider_margins) def _on_canvas_drawn(self, event) -> None: if self._mode_combo.currentIndex() == 1: @@ -828,14 +844,24 @@ def _on_line_updated(self, p0, p1) -> None: if p0 is None or self._current_arr_a is None: self._profile_canvas.clear() return - normalized = self._normalize_check.isChecked() - dist_a, prof_a = _sample_line(self._current_arr_a, p0, p1, normalize=normalized) + normalized = self._normalize_check.isChecked() + subtract_median = self._subtract_median_check.isChecked() + + dist_a, prof_a = _sample_line(self._current_arr_a, p0, p1, normalize=False) + if subtract_median: + prof_a = prof_a - float(np.median(self._current_arr_a)) + if normalized: + peak = float(prof_a.max()) + if peak > 0: + prof_a = prof_a / peak + if self._current_arr_b is None: self._profile_canvas.update_profiles( dist_a, prof_a, None, None, self._left_combo.currentText(), "", - normalized=normalized) + normalized=normalized, subtract_median=subtract_median) return + # Scale line coordinates to the right panel's pixel space if sizes differ. # p0/p1 are in the left panel's (arr_a) pixel coordinate system. ha, wa = self._current_arr_a.shape[:2] @@ -845,12 +871,20 @@ def _on_line_updated(self, p0, p1) -> None: p1b = (p1[0] * wb / wa, p1[1] * hb / ha) else: p0b, p1b = p0, p1 - dist_b, prof_b = _sample_line(self._current_arr_b, p0b, p1b, normalize=normalized) + + dist_b, prof_b = _sample_line(self._current_arr_b, p0b, p1b, normalize=False) + if subtract_median: + prof_b = prof_b - float(np.median(self._current_arr_b)) + if normalized: + peak = float(prof_b.max()) + if peak > 0: + prof_b = prof_b / peak + self._profile_canvas.update_profiles( dist_a, prof_a, dist_b, prof_b, self._left_combo.currentText(), self._right_combo.currentText(), - normalized=normalized) + normalized=normalized, subtract_median=subtract_median) def _on_normalize_changed(self, _state: int) -> None: if self._image_canvas._p0 is not None and self._image_canvas._p1 is not None: From 5f37d54c90e1383e30fed2de4128f47e6e771ab8 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Thu, 2 Jul 2026 16:03:54 -0400 Subject: [PATCH 02/14] Inspector: store linear float32 data in NPZ for unbiased cross-sections MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit STF stretch parameters are computed per-image from each image's own sky median and MAD, so cross-sections of the display image can reflect the nonlinear transform rather than real optical differences between filters. Adds _inspector_linear() to report_builder.py that stores the pre-stretch science data as float32 [0, 1] alongside the existing display images. Integer data (uint16 ADUs) is divided by the dtype max; PixInsight float32 data is already in [0, 1]. The 2048 px max-dimension cap from _inspector_display is reused to control NPZ size growth (~2-4x per image vs uint8). Two new entries appear in the inspector section combo — "Original (linear)" and "Starless (linear)" — giving unbiased PSF and filter transmission cross-sections on a common linear scale. Y-axis label updates to reflect the data type and active checkbox state. Old NPZ files without linear_* keys open cleanly with no change in behaviour. Co-Authored-By: Claude Sonnet 4.6 --- gui/report_inspector.py | 21 +++++++++++++------- report/report_builder.py | 41 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 7 deletions(-) diff --git a/gui/report_inspector.py b/gui/report_inspector.py index 3adf2f4..63f7f89 100644 --- a/gui/report_inspector.py +++ b/gui/report_inspector.py @@ -537,20 +537,22 @@ def update_profiles(self, dist_a: np.ndarray, prof_a: np.ndarray, dist_b: np.ndarray | None, prof_b: np.ndarray | None, label_a: str, label_b: str, normalized: bool = True, - subtract_median: bool = False) -> None: + subtract_median: bool = False, + is_linear: bool = False) -> None: self._ax.cla() self._ax.plot(dist_a, prof_a, color="steelblue", lw=1.2, label=label_a) if dist_b is not None and prof_b is not None: self._ax.plot(dist_b, prof_b, color="tomato", lw=1.2, label=label_b) self._ax.set_xlabel("Distance (px)", fontsize=8) + lin = " (linear)" if is_linear else "" if normalized and subtract_median: - ylabel = "Normalized, bg-subtracted" + ylabel = f"Normalized, bg-subtracted{lin}" elif normalized: - ylabel = "Normalized intensity" + ylabel = f"Normalized{lin}" elif subtract_median: - ylabel = "Background-subtracted" + ylabel = f"Background-subtracted{lin}" else: - ylabel = "Value" + ylabel = f"Linear value [0–1]" if is_linear else "Value" self._ax.set_ylabel(ylabel, fontsize=8) self._ax.legend(fontsize=7, loc="upper right") self._ax.grid(True, alpha=0.25) @@ -584,6 +586,7 @@ def __init__(self, npz_path: str | Path, parent=None): self._current_arr_a: np.ndarray | None = None self._current_arr_b: np.ndarray | None = None self._current_shape: tuple | None = None + self._current_is_linear: bool = False self._build_ui() self._populate_section_combo() @@ -846,6 +849,7 @@ def _on_line_updated(self, p0, p1) -> None: return normalized = self._normalize_check.isChecked() subtract_median = self._subtract_median_check.isChecked() + is_linear = self._current_is_linear dist_a, prof_a = _sample_line(self._current_arr_a, p0, p1, normalize=False) if subtract_median: @@ -859,7 +863,8 @@ def _on_line_updated(self, p0, p1) -> None: self._profile_canvas.update_profiles( dist_a, prof_a, None, None, self._left_combo.currentText(), "", - normalized=normalized, subtract_median=subtract_median) + normalized=normalized, subtract_median=subtract_median, + is_linear=is_linear) return # Scale line coordinates to the right panel's pixel space if sizes differ. @@ -884,7 +889,8 @@ def _on_line_updated(self, p0, p1) -> None: dist_a, prof_a, dist_b, prof_b, self._left_combo.currentText(), self._right_combo.currentText(), - normalized=normalized, subtract_median=subtract_median) + normalized=normalized, subtract_median=subtract_median, + is_linear=is_linear) def _on_normalize_changed(self, _state: int) -> None: if self._image_canvas._p0 is not None and self._image_canvas._p1 is not None: @@ -913,6 +919,7 @@ def _load_current_panels(self) -> None: if key_left is None or key_right is None: return + self._current_is_linear = key_left.startswith("linear_") raw_left = _get_array(self._npz, key_left) if raw_left is None: return diff --git a/report/report_builder.py b/report/report_builder.py index c22d319..276ba05 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -36,6 +36,29 @@ def _inspector_display(img: AstroImage, max_dim: int = 2048) -> np.ndarray: return arr +def _inspector_linear(img: AstroImage, max_dim: int = 2048) -> np.ndarray: + """Return pre-stretch linear data as float32 [0, 1], downsampled so max dimension ≤ max_dim. + + Integer data (uint16 ADUs) is divided by the dtype's max representable value. + PixInsight float32 output is already in [0, 1] and is clipped but not rescaled. + Stored as [0, 1] so the inspector's existing >1.5 guard never fires. + """ + arr = img.data.astype(np.float32) + if img.original_dtype is not None and np.issubdtype(img.original_dtype, np.integer): + arr = arr / float(np.iinfo(img.original_dtype).max) + else: + hi = float(arr.max()) + if hi > 1.0: + arr = arr / hi + arr = np.clip(arr, 0.0, 1.0) + h, w = arr.shape[:2] + scale = min(1.0, max_dim / max(h, w)) + if scale < 1.0: + pil = _PILImage.fromarray(arr) # PIL infers mode='F' for float32 2-D arrays + arr = np.asarray(pil.resize((int(w * scale), int(h * scale)), _PILImage.LANCZOS)) + return arr + + def _make_moffat_kernel(fwhm_px: float, beta: float = REF_SEEING_BETA, size: int | None = None) -> np.ndarray: @@ -594,6 +617,14 @@ def _add_options_entry(section: str, name: str, options: dict) -> None: orig_opts["Image B"] = "display_b" _add_options_entry("Input Images", "Original", orig_opts) + # Linear (pre-stretch) counterparts — float32 [0, 1] + _add("linear_a", _inspector_linear(image_a)) + lin_opts: dict[str, str] = {"Image A": "linear_a"} + if image_b is not None: + _add("linear_b", _inspector_linear(image_b)) + lin_opts["Image B"] = "linear_b" + _add_options_entry("Input Images", "Original (linear)", lin_opts) + sl_opts: dict[str, str] = {} if image_a.starless_image is not None: _add("display_sl_a", _inspector_display(image_a.starless_image)) @@ -604,6 +635,16 @@ def _add_options_entry(section: str, name: str, options: dict) -> None: if sl_opts: _add_options_entry("Input Images", "Starless", sl_opts) + lin_sl_opts: dict[str, str] = {} + if image_a.starless_image is not None: + _add("linear_sl_a", _inspector_linear(image_a.starless_image)) + lin_sl_opts["Image A"] = "linear_sl_a" + if image_b is not None and image_b.starless_image is not None: + _add("linear_sl_b", _inspector_linear(image_b.starless_image)) + lin_sl_opts["Image B"] = "linear_sl_b" + if lin_sl_opts: + _add_options_entry("Input Images", "Starless (linear)", lin_sl_opts) + # ── PSF / MTF ───────────────────────────────────────────────────────── pa = result_a.psf_metrics or {} pb = result_b.psf_metrics or {} From 3f4cfe53e6d04c7b3515658aab98f9e40470ab17 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Fri, 3 Jul 2026 11:53:15 -0400 Subject: [PATCH 03/14] perf: accelerate large-array median/nanmean/nanstd with bottleneck Add bottleneck as a dependency and route NaN-aware and median operations on full-image or near-full-image arrays through bn.* instead of np.*, with a transparent ImportError fallback so environments without bottleneck are unaffected. Hot paths benefiting from this change: - core/stretch.py: two bn.median calls in stf_stretch() and stf_stretch_matched() run on the full flattened image array on every display refresh. - analysis/snr_analyzer.py: bn.median on the 2D background model array. - analysis/halo_analyzer.py: bn.median on stacked radial profiles; bn.nanmean / bn.nanstd on the NaN-padded per-star RDF matrix. analysis/image_filters.py has the same one-line bn.median change in _estimate_noise() but is left out of this commit because it also carries unrelated in-progress work. Co-Authored-By: Claude Sonnet 4.6 --- analysis/halo_analyzer.py | 12 ++++++++---- analysis/snr_analyzer.py | 6 +++++- core/stretch.py | 12 ++++++++---- environment.yml | 1 + 4 files changed, 22 insertions(+), 9 deletions(-) diff --git a/analysis/halo_analyzer.py b/analysis/halo_analyzer.py index f219ed3..48995dc 100644 --- a/analysis/halo_analyzer.py +++ b/analysis/halo_analyzer.py @@ -3,6 +3,10 @@ import warnings import numpy as np +try: + import bottleneck as bn +except ImportError: + bn = np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -92,7 +96,7 @@ def analyze(self, image: AstroImage) -> dict: result["star_data"] = per_star common_r = profiles[0][0] - stacked = np.median( + stacked = bn.median( np.array([np.interp(common_r, p[0], p[1]) for p in profiles]), axis=0 ) @@ -327,9 +331,9 @@ def _compute_aggregate_rdf(self, rdf_list: list) -> dict | None: all_stds = np.array(stacked_stds) return { "rdf_radii": r_ref, - "rdf_mean": np.nanmean(all_means, axis=0), - "rdf_std": np.nanmean(all_stds, axis=0), # avg per-bin σ - "rdf_star_std": np.nanstd(all_means, axis=0), # star-to-star σ + "rdf_mean": bn.nanmean(all_means, axis=0), + "rdf_std": bn.nanmean(all_stds, axis=0), # avg per-bin σ + "rdf_star_std": bn.nanstd(all_means, axis=0), # star-to-star σ "rdf_n_stars": len(stacked_means), } diff --git a/analysis/snr_analyzer.py b/analysis/snr_analyzer.py index 7e5dc37..a55ab57 100644 --- a/analysis/snr_analyzer.py +++ b/analysis/snr_analyzer.py @@ -5,6 +5,10 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np +try: + import bottleneck as bn +except ImportError: + bn = np from core.astro_image import AstroImage from core.fig_utils import fig_to_b64 @@ -79,7 +83,7 @@ def analyze(self, image: AstroImage) -> dict: # --- Sky background level ---------------------------------------- # Median of the 2D background model in ADU; lower = darker sky. - background_median = (float(np.median(image.background.background)) + background_median = (float(bn.median(image.background.background)) if image.background is not None else None) # --- Camera gain from FITS header -------------------------------- diff --git a/core/stretch.py b/core/stretch.py index e14992f..dd17b48 100644 --- a/core/stretch.py +++ b/core/stretch.py @@ -1,6 +1,10 @@ from __future__ import annotations import numpy as np +try: + import bottleneck as bn +except ImportError: + bn = np def normalize_unit_interval(data: np.ndarray) -> np.ndarray: @@ -44,8 +48,8 @@ def stf_stretch(data: np.ndarray, finite = norm[np.isfinite(norm)] if finite.size == 0: return np.zeros_like(norm, dtype=np.float32) - med = float(np.median(finite)) - mad = float(np.median(np.abs(finite - med))) + med = float(bn.median(finite)) + mad = float(bn.median(np.abs(finite - med))) c = med + shadow_clip * 1.4826 * mad c = float(np.clip(c, 0.0, 1.0)) if not np.isfinite(c): @@ -87,8 +91,8 @@ def stf_stretch_matched(data: np.ndarray, ref: np.ndarray) -> np.ndarray: finite_ref = norm_ref[np.isfinite(norm_ref)] if finite_ref.size == 0: return norm_data - med = float(np.median(finite_ref)) - mad = float(np.median(np.abs(finite_ref - med))) + med = float(bn.median(finite_ref)) + mad = float(bn.median(np.abs(finite_ref - med))) c = float(np.clip(med - 2.8 * 1.4826 * mad, 0.0, 1.0 - 1e-6)) denom = 1.0 - c if denom <= 0.0: diff --git a/environment.yml b/environment.yml index f551f5e..08b7408 100644 --- a/environment.yml +++ b/environment.yml @@ -12,6 +12,7 @@ dependencies: - matplotlib - astropy - photutils + - bottleneck - PyWavelets - astroalign - Pillow From dc1173e51bb891bc4572cac874c5fcc2c63ea140 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Fri, 3 Jul 2026 11:55:31 -0400 Subject: [PATCH 04/14] finalizing bottleneck updates# On branch dev_NebSynth --- CLAUDE.md | 102 ++++++ README.md | 5 +- analysis/image_filters.py | 104 ++++++- analysis/power_spectrum.py | 26 +- core/models.py | 5 +- gui/main_window.py | 23 ++ gui/spatial_target_dialog.py | 566 ++++++++++++++++++++++++++++++++++ report/report_builder.py | 51 ++- synthetic/target_generator.py | 282 +++++++++++++++++ 9 files changed, 1147 insertions(+), 17 deletions(-) create mode 100644 gui/spatial_target_dialog.py create mode 100644 synthetic/target_generator.py diff --git a/CLAUDE.md b/CLAUDE.md index ff8d272..4d2c752 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -156,6 +156,40 @@ from scipy.signal import fftconvolve convolved = fftconvolve(patch, psf_kernel, mode="same").astype(np.float64) ``` +### Large-array reductions — prefer bottleneck + +`bottleneck` (conda-forge) provides drop-in replacements for the numpy NaN-aware +and median functions that are substantially faster on large arrays (full-image or +background-map sized). Use it for any reduction that operates on arrays with +`size > ~10 000` elements. Always import with a transparent fallback: + +```python +try: + import bottleneck as bn +except ImportError: + bn = np # transparent fallback; bn = np must come after import numpy as np +``` + +**Use `bn.*` instead of `np.*` for these functions on large arrays:** + +| numpy | bottleneck | Notes | +| --- | --- | --- | +| `np.median(a)` | `bn.median(a)` | Supports `axis=` parameter | +| `np.nanmedian(a)` | `bn.nanmedian(a)` | Supports `axis=` parameter | +| `np.nanmean(a, axis=)` | `bn.nanmean(a, axis=)` | NaN-aware row/col aggregation | +| `np.nanstd(a, axis=)` | `bn.nanstd(a, axis=)` | Same default `ddof=0` as numpy | +| `np.nansum(a)` | `bn.nansum(a)` | Only worthwhile when NaNs are actually present | + +**Do not replace:** + +- `np.nanpercentile` / `np.percentile` — bottleneck has no equivalent. +- Any reduction on arrays with fewer than ~1 000 elements — call overhead dominates. + +**Currently in use:** `core/stretch.py` (stf_stretch, stf_stretch_matched), +`analysis/snr_analyzer.py` (background model median), +`analysis/halo_analyzer.py` (stacked radial profiles, RDF nanmean/nanstd), +`analysis/image_filters.py` (wavelet MAD noise estimate). + --- ## Collaboration Rules @@ -248,6 +282,8 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | Closure capture in `secondary_xaxis` lambdas | `lambda x: x * ps` inside a loop captures `ps` by reference. Use default-arg capture: `lambda x, p=ps: x * p` to freeze the value at definition time. | | macOS binary blocked by Gatekeeper | CI-built binaries are unsigned. Users must right-click → Open, or run `xattr -dr com.apple.quarantine AstroImageLab` in Terminal. Code signing requires an Apple Developer certificate ($99/year). | | Linux build needs system Qt libraries | PyInstaller must be able to import PyQt6 during analysis. On `ubuntu-latest` run `sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0` before `pip install -r requirements-build.txt`. | +| `PowerSpectrumAnalyzer` crashes on images smaller than 1024 px | `POWER_SPECTRUM_NPIX = 1024`. The auto-select loop is empty when `min(h, w) < 1024`; the fallback produces negative slice indices → non-square region → `_apply_window` shape mismatch. Fix: `N = min(N, h, w)` before the loop, add `+1` to loop upper bounds, clamp fallback with `max(0, ...)`. | +| `sigma_clip` mask is scalar `False` when nothing is clipped | `clipped.mask` is `np.ma.nomask` (== `False`) when no values are clipped. `region[False]` silently writes only the first row. Use `np.ma.getmaskarray(clipped)` to get a full bool array, then guard with `.any()`. | --- @@ -396,6 +432,72 @@ reject valid matches when residual alignment offset exists. --- +## Spatial Target Generator — Key Patterns + +### Purpose and workflow + +`gui/spatial_target_dialog.py` + `synthetic/target_generator.py` + +Generates a 4-column × 3-row grid of calibrated test zones at known spatial frequencies, +always as a clean/degraded FITS pair. Load clean → Image A and degraded → Image B to +calibrate the spatial-detail metrics against known inputs. + +### Target signal chain + +```text +SpatialTargetDialog.targets_generated = pyqtSignal(str, str, str) # clean_path, degraded_path, mode + → MainWindow._on_target_generated(clean_path, degraded_path, mode) + # mode: "clean_a_deg_b" | "deg_a_clean_b" | "deg_a" | "deg_b" | "" +``` + +`_TargetThread.finished = pyqtSignal(str, str)` (clean, degraded) feeds `_on_gen_done` +which then emits the three-arg `targets_generated` signal. + +### Target return types + +`SpatialTargetGenerator.generate(params, preview=False)`: + +- `preview=True` → `np.ndarray` (float32, degraded image at reduced resolution) +- `preview=False` → `tuple[str, str]` (clean_path, degraded_path) + +### Zone layout + +```text +Row 0: Sine H f=0.04 | Sine H f=0.08 | Sine H f=0.16 | Sine H f=0.32 (c/px) +Row 1: Square H 0.04 | Square H 0.08 | Square H 0.16 | Square H 0.32 +Row 2: Sine V 0.08 | Sine 45° 0.08 | Siemens star | Slant edge ~5° +``` + +Column frequencies align with wavelet levels: 0.04→L4, 0.08→L3, 0.16→L2, 0.32→L1. + +### Contrast ramp + +Each zone ramps Michelson contrast linearly from `contrast_min` (top edge) to +`contrast_max` (bottom edge). At any horizontal strip, all four columns share +the same contrast — enabling direct cross-frequency comparison. Params: +`contrast_min`, `contrast_max` (both 0–1, default 0.02 / 0.50). + +### FITS keywords + +`INSTRUME="SpatialTarget"`, `EGAIN=1.0`, `GAIN=1.0`, `TGT_TYPE`, `TGT_ROWS`, +`TGT_COLS`, `TGT_CMIN`, `TGT_CMAX`, `TGT_SKY`, `TGT_CLEN` (bool: clean flag), +per-zone `TGT_{r}{c}F` / `TGT_{r}{c}W`. Optional: `TGT_FWHM`, `TGT_BETA`, `TGT_RN`. +No `FOCALLEN`, `APTDIA`, `FOCRATIO`, `XPIXSZ`, `EXPTIME` — these are set from +`AstroImage` defaults (pixel scale = `DEFAULT_PIXEL_SCALE`). + +### Power spectrum on spatial target images + +The power spectrum auto-selects one square ROI — not the whole zone grid. The result +reflects whichever zone(s) fall inside that square. Use the explicit crosshair ROI +(user-drawn in the image panel) to target a specific zone for a focused power spectrum. +The frequency axis is always cycles/pixel regardless of ROI size. + +### QSettings key + +`"target_output_dir"` (separate from the synthetic dialog's `"synth_output_dir"`). + +--- + ## Working Effectively with Claude Code ### The most useful problem statement format diff --git a/README.md b/README.md index 4778eb4..74d9573 100644 --- a/README.md +++ b/README.md @@ -73,7 +73,7 @@ This creates a conda environment named `astrolab` with all required packages. ** Install the scientific stack via conda-forge, then add PyQt6 and XISF support via pip: ```bash -conda install -c conda-forge numpy scipy matplotlib astropy photutils pywavelets astroalign pillow lz4 zstandard +conda install -c conda-forge numpy scipy matplotlib astropy photutils bottleneck pywavelets astroalign pillow lz4 zstandard pip install pyqt6 xisf ``` @@ -92,7 +92,7 @@ conda env create -f environment.yml conda activate astrolab # Option B — manual install into an existing environment -conda install -c conda-forge numpy scipy matplotlib astropy photutils pywavelets astroalign pillow lz4 zstandard +conda install -c conda-forge numpy scipy matplotlib astropy photutils bottleneck pywavelets astroalign pillow lz4 zstandard pip install pyqt6 xisf ``` @@ -230,6 +230,7 @@ gui/ |---------|---------| | [astropy](https://www.astropy.org/) | FITS I/O, Moffat2D model, Background2D | | [photutils](https://photutils.readthedocs.io/) | DAOStarFinder, EPSFBuilder, morphology | +| [bottleneck](https://bottleneck.readthedocs.io/) | Fast NaN-aware and median reductions on large arrays | | [scipy](https://scipy.org/) | Optimisation, FFT, image filters | | [PyWavelets](https://pywavelets.readthedocs.io/) | Daubechies-4 wavelet decomposition | | [astroalign](https://astroalign.quatrope.org/) | Image registration | diff --git a/analysis/image_filters.py b/analysis/image_filters.py index 08d609a..2162b4d 100644 --- a/analysis/image_filters.py +++ b/analysis/image_filters.py @@ -2,16 +2,21 @@ import concurrent.futures import numpy as np +try: + import bottleneck as bn +except ImportError: + bn = np import matplotlib import matplotlib.colors as mcolors matplotlib.use("Agg") import matplotlib.pyplot as plt -from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, map_coordinates, zoom +from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, map_coordinates, zoom, maximum_filter, minimum_filter import pywt from core.astro_image import AstroImage from core.fig_utils import fig_to_b64, figs_to_b64 from core.models import (STD_KERNEL_SIZES, LOG_SIGMAS, WAVELET_NAME, WAVELET_LEVELS, + MICHELSON_KERNEL_SIZES, XS_LINE_ALPHA, SECTION8_BORDER_CROP_FRACTION, SECTION8_ANALYSIS_CMAP, XS_SNR_REGION_WIDTH) @@ -31,6 +36,7 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, log_sigmas: tuple = LOG_SIGMAS, wavelet: str = WAVELET_NAME, levels: int = WAVELET_LEVELS, + michelson_kernel_sizes: tuple = MICHELSON_KERNEL_SIZES, crosshair: dict | None = None, roi: tuple | None = None, xs_snr_width: int | None = None) -> dict: @@ -49,6 +55,8 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, "wavelet_snr_b": {}, "sigma_noise_a": None, "sigma_noise_b": None, + "michelson_contrast_a": {}, + "michelson_contrast_b": {}, "panels": {}, } figures: dict = {} @@ -105,7 +113,7 @@ def _clip01(v): return max(0.0, min(1.0, v)) # 1-3. Local std, LoG, wavelet — all read norm_a/norm_b with no shared mutable state, # so they run concurrently. Each method returns (b64_figs, partial_result). - with concurrent.futures.ThreadPoolExecutor(max_workers=3) as _ex: + with concurrent.futures.ThreadPoolExecutor(max_workers=4) as _ex: _f_std = _ex.submit(self._std_analysis, analysis_a, analysis_b, mask_neb_a, mask_bg_a, @@ -127,22 +135,33 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, ) + _f_mic = _ex.submit(self._michelson_analysis, + analysis_a, analysis_b, + michelson_kernel_sizes, + image_a.label, _label_b, + display_roi=display_roi, + ) std_b64, std_partial = _f_std.result() log_b64, log_partial = _f_log.result() wav_b64, wav_partial = _f_wav.result() + mic_b64, mic_partial = _f_mic.result() figures.update(std_b64) figures.update(log_b64) figures.update(wav_b64) + figures.update(mic_b64) result["contrast_ratios_a"].update(std_partial["contrast_ratios_a"]) result["contrast_ratios_b"].update(std_partial["contrast_ratios_b"]) result["sigma_noise_a"] = wav_partial["sigma_noise_a"] result["sigma_noise_b"] = wav_partial["sigma_noise_b"] result["wavelet_snr_a"].update(wav_partial["wavelet_snr_a"]) result["wavelet_snr_b"].update(wav_partial["wavelet_snr_b"]) + result["michelson_contrast_a"].update(mic_partial["michelson_contrast_a"]) + result["michelson_contrast_b"].update(mic_partial["michelson_contrast_b"]) result["panels"].update(std_partial["panels"]) result["panels"].update(log_partial["panels"]) result["panels"].update(wav_partial["panels"]) + result["panels"].update(mic_partial["panels"]) if crosshair is not None: pos_a, prof_a = self._sample_line(norm_a, **crosshair) @@ -457,7 +476,7 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, def _estimate_noise(self, coeffs) -> float: # Finest-level horizontal detail (last element, first sub-band) lh1 = coeffs[-1][0] - return float(np.median(np.abs(lh1))) / 0.6745 + return float(bn.median(np.abs(lh1))) / 0.6745 def _level_snr(self, coeffs, coeff_idx: int, sigma_noise: float, human_level: int) -> float | None: @@ -486,6 +505,85 @@ def _reconstruct_level(self, coeffs, target_coeff_idx: int, # Shared plotting helper # ------------------------------------------------------------------ + # ------------------------------------------------------------------ + # Michelson contrast maps + # ------------------------------------------------------------------ + + def _michelson_analysis(self, norm_a, norm_b, kernel_sizes, + label_a, label_b, + display_roi=None) -> tuple[dict, dict]: + figures = {} + partial: dict = { + "michelson_contrast_a": {}, + "michelson_contrast_b": {}, + "panels": {}, + } + single = norm_b is None + for ks in kernel_sizes: + mc_a = self._compute_michelson_map(norm_a, ks) + mc_b = self._compute_michelson_map(norm_b, ks) if not single else None + + partial["michelson_contrast_a"][ks] = float(np.max(mc_a)) + if not single: + partial["michelson_contrast_b"][ks] = float(np.max(mc_b)) + + partial["panels"][f"michelson_{ks}px"] = { + "a": mc_a.astype(np.float32), + "b": mc_b.astype(np.float32) if mc_b is not None else None, + "diff": (mc_a - mc_b).astype(np.float32) if mc_b is not None else None, + } + + if not single: + fig = self._plot_side_by_side( + self._crop_border(mc_a, SECTION8_BORDER_CROP_FRACTION), + self._crop_border(mc_b, SECTION8_BORDER_CROP_FRACTION), + f"Michelson contrast — kernel {ks}px — {label_a}", + f"Michelson contrast — kernel {ks}px — {label_b}", + diff_title=f"Michelson diff (A−B), kernel {ks}px", + cmap=SECTION8_ANALYSIS_CMAP, + nonlinear_norm=False, # values bounded [0,1]; linear scale is appropriate + display_roi=None, # _crop_border already applied; matches _log_analysis pattern + ) + else: + fig = self._plot_single( + self._crop_border(mc_a, SECTION8_BORDER_CROP_FRACTION), + f"Michelson contrast — kernel {ks}px — {label_a}", + cmap=SECTION8_ANALYSIS_CMAP, + nonlinear_norm=False, + ) + figures[f"michelson_{ks}px"] = fig + + return figs_to_b64(figures, dpi=150), partial + + def _compute_michelson_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray: + """Per-pixel Michelson contrast C = (I_max - I_min) / (I_max + I_min + eps). + Downsamples large images identically to _compute_std_map for performance. + """ + _EPS = 1e-9 + factor = 1.0 + data = norm + if max(norm.shape) > MAX_DIM_FOR_STD: + factor = MAX_DIM_FOR_STD / max(norm.shape) + new_h = int(norm.shape[0] * factor) + new_w = int(norm.shape[1] * factor) + data = zoom(norm, (new_h / norm.shape[0], new_w / norm.shape[1]), order=1) + kernel_size = max(3, int(kernel_size * factor) | 1) + + i_max = maximum_filter(data, size=kernel_size) + i_min = minimum_filter(data, size=kernel_size) + # Background-subtracted images can have negative pixels after mean normalisation. + # Clip to 0 so i_max + i_min is never negative, preventing spurious near-1 contrast. + i_max = np.maximum(i_max, 0.0) + i_min = np.maximum(i_min, 0.0) + mc_map = (i_max - i_min) / (i_max + i_min + _EPS) + + if factor < 1.0: + mc_map = zoom(mc_map, + (norm.shape[0] / mc_map.shape[0], + norm.shape[1] / mc_map.shape[1]), + order=1) + return np.clip(mc_map, 0.0, 1.0) + def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray, title_a: str, title_b: str, diff_title: str = "", diff --git a/analysis/power_spectrum.py b/analysis/power_spectrum.py index f50ed68..d563350 100644 --- a/analysis/power_spectrum.py +++ b/analysis/power_spectrum.py @@ -78,6 +78,10 @@ def _extract_roi(self, bgsub: np.ndarray, sub = zoom(sub, (zy, zx), order=1) return sub[:N, :N] + # For images smaller than POWER_SPECTRUM_NPIX, use the largest square + # that fits — the c/px frequency axis interpretation is unchanged. + N = min(N, h, w) + # Auto-select: find a star-free NxN region catalog = getattr(image, "catalog", None) if catalog is not None and len(catalog) > 0: @@ -86,12 +90,13 @@ def _extract_roi(self, bgsub: np.ndarray, else: star_xs, star_ys = np.array([]), np.array([]) - # Try candidate positions on a grid - step = N // 2 + # Try candidate positions on a grid; +1 so images exactly N px wide/tall + # still yield one candidate at their centre. + step = max(1, N // 2) best_pos = None best_score = -np.inf - for yc in range(N // 2, h - N // 2, step): - for xc in range(N // 2, w - N // 2, step): + for yc in range(N // 2, h - N // 2 + 1, step): + for xc in range(N // 2, w - N // 2 + 1, step): x0 = xc - N // 2 y0 = yc - N // 2 if len(star_xs) > 0: @@ -105,17 +110,20 @@ def _extract_roi(self, bgsub: np.ndarray, best_pos = (x0, y0) if best_pos is None: - # Fallback: centre of image - x0 = w // 2 - N // 2 - y0 = h // 2 - N // 2 + # Fallback: centre of image — clamp so coordinates are never negative + x0 = max(0, w // 2 - N // 2) + y0 = max(0, h // 2 - N // 2) best_pos = (x0, y0) x0, y0 = best_pos region = bgsub[y0:y0 + N, x0:x0 + N].copy() - # Sigma-clip to remove unmasked faint stars + # Sigma-clip to remove unmasked faint stars; getmaskarray always returns + # a full bool array (avoids scalar-False indexing when nothing is clipped). clipped = sigma_clip(region, sigma=3.0, maxiters=3) - region[clipped.mask] = float(np.ma.median(clipped)) + mask = np.ma.getmaskarray(clipped) + if mask.any(): + region[mask] = float(np.ma.median(clipped)) return region # ------------------------------------------------------------------ diff --git a/core/models.py b/core/models.py index b337110..dce7a52 100644 --- a/core/models.py +++ b/core/models.py @@ -26,10 +26,11 @@ FILTER_THICKNESS_MM = 1.0 # narrowband filter substrate thickness (mm); default for UI GLASS_REFRACTIVE_INDEX = 1.9 # dichroic filter substrate refractive index RDF_BIN_WIDTH = 1.0 # px; annular bin width for RDF mean/std computation -POWER_SPECTRUM_NPIX = 1024 +POWER_SPECTRUM_NPIX = 2048 # px; size of square power spectrum maps (must be 2^n) -STD_KERNEL_SIZES = (5, 15, 31) +STD_KERNEL_SIZES = (3, 5, 10) #originally (5, 10, 15) # px; Gaussian kernel sizes for std dev maps LOG_SIGMAS = (1.5, 3.0, 6.0) +MICHELSON_KERNEL_SIZES = (3, 5, 9) # px; local max/min neighbourhood for Michelson contrast (odd values required) WAVELET_NAME = "db4" WAVELET_LEVELS = 4 diff --git a/gui/main_window.py b/gui/main_window.py index ee7f0fe..e57fd6d 100644 --- a/gui/main_window.py +++ b/gui/main_window.py @@ -112,6 +112,10 @@ def _build_menu(self) -> None: act_halo.triggered.connect(self._open_halo_dialog) tools_menu.addAction(act_halo) + act_target = QAction("Spatial &Target…", self) + act_target.triggered.connect(self._open_target_dialog) + tools_menu.addAction(act_target) + help_menu = mb.addMenu("&Help") act_quickstart = QAction("&Quick Start Guide", self) act_quickstart.triggered.connect(self._open_quickstart) @@ -149,6 +153,25 @@ def _on_synthetic_generated(self, main_path: str, starless_path: str, if starless_path: target.set_starless_path(starless_path) + def _open_target_dialog(self) -> None: + from gui.spatial_target_dialog import SpatialTargetDialog + self._target_dialog = SpatialTargetDialog(parent=self) + self._target_dialog.targets_generated.connect(self._on_target_generated) + self._target_dialog.show() + + def _on_target_generated(self, clean_path: str, degraded_path: str, + mode: str) -> None: + if mode == "clean_a_deg_b": + self._panel_a.load_path(clean_path) + self._panel_b.load_path(degraded_path) + elif mode == "deg_a_clean_b": + self._panel_a.load_path(degraded_path) + self._panel_b.load_path(clean_path) + elif mode == "deg_a": + self._panel_a.load_path(degraded_path) + elif mode == "deg_b": + self._panel_b.load_path(degraded_path) + def _on_image_loaded(self, img) -> None: either_loaded = (self._panel_a.image is not None or self._panel_b.image is not None) diff --git a/gui/spatial_target_dialog.py b/gui/spatial_target_dialog.py new file mode 100644 index 0000000..b6f1d38 --- /dev/null +++ b/gui/spatial_target_dialog.py @@ -0,0 +1,566 @@ +"""Spatial Detail Test Target Generator dialog. + +Generates a 4×3 grid of calibrated spatial test patterns with known frequency +content, enabling calibrated interpretation of the app's spatial-detail metrics. + +Workflow: + 1. Configure contrast, sky level, and optional PSF / noise. + 2. Click Generate → produces a clean reference (no PSF/noise) and a + degraded companion (with configured blur and noise). + 3. Auto-load maps clean → Image A, degraded → Image B (configurable). + 4. Run the analysis; compare wavelet SNR, local-std ratios, and LoG maps + against the known zone frequencies to calibrate your metric intuition. +""" + +from __future__ import annotations + +import os +from pathlib import Path + +import numpy as np +from PyQt6.QtCore import Qt, QThread, QTimer, pyqtSignal, pyqtSlot +from PyQt6.QtGui import QImage, QPixmap +from PyQt6.QtWidgets import ( + QCheckBox, QComboBox, QDialog, QDialogButtonBox, QDoubleSpinBox, + QFileDialog, QFormLayout, QGroupBox, QHBoxLayout, QLabel, + QMainWindow, QPushButton, QScrollArea, QSizePolicy, + QTextBrowser, QVBoxLayout, QWidget, +) + +from core.stretch import normalize_for_display +from synthetic.target_generator import SpatialTargetGenerator, ZONE_LABELS +from gui.synthetic_dialog import _SliderRow + + +# --------------------------------------------------------------------------- +# Zone layout info HTML (shown in the ⓘ dialog) +# --------------------------------------------------------------------------- + +_ZONE_INFO_HTML = """ +

Spatial Detail Test Target — Zone Layout

+ +

The target is a 4-column × 3-row grid of calibrated test zones. +Each zone has a known waveform type and spatial frequency. +All zones share the same Michelson contrast and DC sky level.

+ +

Zone Grid

+ + + + + + + + + + + + + + + + + +
RowCol 0Col 1Col 2Col 3
0Sine H
f=0.04 c/px		Sine H
f=0.08 c/px		Sine H
f=0.16 c/px		Sine H
f=0.32 c/px
1Square H
f=0.04 c/px		Square H
f=0.08 c/px		Square H
f=0.16 c/px		Square H
f=0.32 c/px
2Sine V
f=0.08 c/px		Sine 45°
f=0.08 c/px		Siemens star
(radial chirp)		Slant edge
(~5° from vertical)
+ +

Frequency–Wavelet Level Correspondence

+

The four column frequencies are chosen to align with the app's 4-level +wavelet decomposition band centres:

+ + + + + + + +
ColumnFrequency (c/px)Wavelength (px)Wavelet levelApprox scale
00.0425 px4 (coarsest)~16 px
10.0812 px3~8 px
20.166 px2~4 px
30.323 px1 (finest / near Nyquist)~2 px
+ +

Contrast Ramp

+

Each zone ramps Michelson contrast linearly from top to bottom: +the top edge has Contrast min and the bottom edge has Contrast max. +Within any horizontal strip across a row, all four columns share the same contrast — +so you can directly compare how different frequencies or waveform types respond at +identical contrast levels. The ramp restarts at each zone row, so row 0 (sine), +row 1 (square), and row 2 (special) all independently cover the full min→max range.

+ +

How to Use the Target

+
    +
  1. Clean reference (Image A): All zones have perfect contrast at their + design frequency. Wavelet SNR for zone (0,0) should be high at level 4 and + near zero at levels 1–3. Local std at 31 px should respond strongly to the + f=0.04 zone. Use this to establish the ceiling for each metric.
    2. +
  2. Degraded companion (Image B): After PSF convolution, high-frequency + zones lose contrast first. The f=0.32 zone (near Nyquist) should show + near-zero contrast even for moderate FWHM. Wavelet level 1 SNR drops to + near 1.0 (noise floor) when the PSF FWHM exceeds ~3 px.
    3. +
  3. Michelson contrast vs frequency: To measure MTF at each frequency, + load the clean and degraded images, measure local std or wavelet SNR in the + sine zone for each column. Ratio (degraded / clean) gives the MTF at that + spatial frequency.
    4. +
+ +

Waveform Types

+
+
Sine waves
+
Pure single-frequency content. After PSF convolution the Michelson contrast +in the zone equals the MTF at that frequency. Use these for the cleanest +frequency-response measurements.
+
Square waves
+
Contain all odd harmonics (f, 3f, 5f, …). PSF blur first kills the high +harmonics, making the square look trapezoidal then sinusoidal. Gibbs ringing +artifacts from narrow-bandpass filters appear here but not in the sine zones.
+
Siemens star
+
Radial sinusoidal pattern that sweeps from low frequency at the edge to +very high frequency at the centre (where the pattern aliases). Shows all +spatial frequencies and all orientations in one zone. The radius at which +contrast fades to noise is the resolution limit.
+
Slant edge (~5° from vertical)
+
ISO 12233-style step function. The slight tilt oversamples the edge +spread function (ESF) by ~1/sin(5°) ≈ 11×, enabling sub-pixel MTF measurement. +After PSF convolution, the app's edge analyzer can extract the MTF from this +zone directly.
+
+""" + + +# --------------------------------------------------------------------------- +# Background threads +# --------------------------------------------------------------------------- + +class _TargetThread(QThread): + finished = pyqtSignal(str, str) # clean_path, degraded_path + failed = pyqtSignal(str) + + def __init__(self, params: dict, parent=None): + super().__init__(parent) + self._params = params + + def run(self): + try: + gen = SpatialTargetGenerator() + clean_path, degraded_path = gen.generate(self._params, preview=False) + self.finished.emit(clean_path, degraded_path) + except Exception as exc: + self.failed.emit(str(exc)) + + +class _PreviewThread(QThread): + finished = pyqtSignal(object) # np.ndarray + + def __init__(self, params: dict, parent=None): + super().__init__(parent) + self._params = params + + def run(self): + try: + gen = SpatialTargetGenerator() + arr = gen.generate(self._params, preview=True) + self.finished.emit(arr) + except Exception: + pass + + +# --------------------------------------------------------------------------- +# Info dialog +# --------------------------------------------------------------------------- + +class _ZoneInfoDialog(QDialog): + def __init__(self, parent=None): + super().__init__(parent) + self.setWindowTitle("Spatial Target — Zone Layout & Interpretation") + self.setMinimumSize(680, 540) + lay = QVBoxLayout(self) + browser = QTextBrowser() + browser.setHtml(_ZONE_INFO_HTML) + browser.setOpenExternalLinks(False) + lay.addWidget(browser) + btns = QDialogButtonBox(QDialogButtonBox.StandardButton.Close) + btns.rejected.connect(self.accept) + lay.addWidget(btns) + + +# --------------------------------------------------------------------------- +# Constants +# --------------------------------------------------------------------------- + +_SIZE_PRESETS = [ + ("512 × 384 (fast)", 512, 384), + ("1024 × 768 (default)", 1024, 768), + ("2048 × 1536 (slow)", 2048, 1536), +] + +_LOAD_MODES = [ + ("Clean → A, Degraded → B", "clean_a_deg_b"), + ("Degraded → A, Clean → B", "deg_a_clean_b"), + ("Degraded → A only", "deg_a"), + ("Degraded → B only", "deg_b"), + ("Don't load", ""), +] + + +# --------------------------------------------------------------------------- +# Main dialog window +# --------------------------------------------------------------------------- + +class SpatialTargetDialog(QMainWindow): + """Standalone window for generating spatial detail test targets.""" + + # clean_path, degraded_path, load_mode string (see _LOAD_MODES) + targets_generated = pyqtSignal(str, str, str) + + def __init__(self, parent=None): + super().__init__(parent) + self.setWindowTitle("AstroImageLab — Spatial Detail Test Target") + self.setMinimumSize(820, 520) + + self._gen_thread: _TargetThread | None = None + self._prev_thread: _PreviewThread | None = None + self._preview_timer = QTimer(self) + self._preview_timer.setSingleShot(True) + self._preview_timer.setInterval(400) + self._preview_timer.timeout.connect(self._run_preview) + + self._build_ui() + self._preview_cb.setChecked(True) + + # Restore persisted output directory + from PyQt6.QtCore import QSettings + saved = QSettings("FilterImageComparator", + "FilterImageComparator").value("target_output_dir", "") + if saved and os.path.isdir(saved): + self._outdir_lbl.setText(saved) + + # ------------------------------------------------------------------ + # UI construction + # ------------------------------------------------------------------ + + def _build_ui(self) -> None: + central = QWidget() + self.setCentralWidget(central) + root = QVBoxLayout(central) + root.setContentsMargins(8, 8, 8, 8) + root.setSpacing(6) + + top = QHBoxLayout() + top.setSpacing(8) + top.addWidget(self._build_pattern_setup(), stretch=1) + top.addWidget(self._build_degradation(), stretch=1) + top.addWidget(self._build_output(), stretch=1) + root.addLayout(top) + root.addWidget(self._build_preview_section()) + + def _build_pattern_setup(self) -> QGroupBox: + box = QGroupBox("Pattern Setup") + form = QFormLayout(box) + form.setVerticalSpacing(5) + + # Image size preset combo + self._size_combo = QComboBox() + for label, *_ in _SIZE_PRESETS: + self._size_combo.addItem(label) + self._size_combo.setCurrentIndex(1) # 1024×768 default + self._size_combo.setToolTip( + "Image dimensions for the generated target.\n" + "Larger images show more cycles at low frequencies but take longer to generate.") + self._size_combo.currentIndexChanged.connect(self._schedule_preview) + form.addRow("Image size", self._size_combo) + + # Contrast ramp — min (top of zone) and max (bottom of zone) + self._contrast_min_sl = _SliderRow(0.00, 1.0, 0.02, 2) + self._contrast_min_sl.setToolTip( + "Michelson contrast at the TOP edge of each zone.\n\n" + "C = (I_max − I_min) / (I_max + I_min)\n\n" + "Setting this near 0 means the top rows of each zone are nearly flat\n" + "(pattern invisible) — this shows the detection floor for each metric.\n" + "Contrast increases linearly toward the bottom (Contrast max).") + self._contrast_min_sl.valueChanged.connect(self._schedule_preview) + form.addRow("Contrast min (top)", self._contrast_min_sl) + + self._contrast_max_sl = _SliderRow(0.00, 1.0, 0.50, 2) + self._contrast_max_sl.setToolTip( + "Michelson contrast at the BOTTOM edge of each zone.\n\n" + "C = (I_max − I_min) / (I_max + I_min)\n\n" + "Within any horizontal strip across a row, all four columns share\n" + "the same contrast — enabling direct comparison across frequencies\n" + "or waveform types at identical contrast levels.\n\n" + "Typical range: min=0.02, max=0.50 covers the full range from just\n" + "above noise floor to clearly-visible nebula structure.") + self._contrast_max_sl.valueChanged.connect(self._schedule_preview) + form.addRow("Contrast max (bottom)", self._contrast_max_sl) + + # Sky background ADU + self._sky_spin = QDoubleSpinBox() + self._sky_spin.setRange(10, 50_000) + self._sky_spin.setValue(500) + self._sky_spin.setDecimals(0) + self._sky_spin.setSuffix(" ADU") + self._sky_spin.setToolTip( + "DC background level in ADU. Zone patterns are ±amplitude offsets\n" + "around this value where amplitude = sky × Michelson contrast.\n\n" + "Use a value comparable to your actual sky background so the wavelet\n" + "SNR and local-std metrics operate in a similar numerical regime to\n" + "real images. Typical values: 200–2000 ADU.") + self._sky_spin.valueChanged.connect(self._schedule_preview) + form.addRow("Sky DC level", self._sky_spin) + + # Info button + info_btn = QPushButton("ⓘ Zone layout & interpretation…") + info_btn.setToolTip( + "Open a description of what each zone measures and how to\n" + "interpret the spatial-detail metrics against these test patterns.") + info_btn.clicked.connect(self._show_info) + form.addRow(info_btn) + + return box + + def _build_degradation(self) -> QGroupBox: + box = QGroupBox("Degradation (applies to degraded image only)") + form = QFormLayout(box) + form.setVerticalSpacing(5) + + # PSF section + self._psf_cb = QCheckBox("Apply PSF blur") + self._psf_cb.setToolTip( + "Convolve the target with an isotropic Moffat PSF.\n" + "Models atmospheric seeing without field-dependent aberrations.\n\n" + "After blur, high-frequency zones lose contrast first:\n" + " f=0.32 zone loses ~50 % contrast when FWHM ≈ 1 px\n" + " f=0.08 zone loses ~50 % contrast when FWHM ≈ 4 px\n" + "The ratio (blurred contrast / clean contrast) = MTF at that frequency.") + self._psf_cb.toggled.connect(self._on_psf_toggled) + self._psf_cb.toggled.connect(self._schedule_preview) + form.addRow(self._psf_cb) + + self._fwhm_spin = QDoubleSpinBox() + self._fwhm_spin.setRange(0.5, 50.0) + self._fwhm_spin.setValue(4.0) + self._fwhm_spin.setSuffix(" px") + self._fwhm_spin.setSingleStep(0.5) + self._fwhm_spin.setDecimals(1) + self._fwhm_spin.setEnabled(False) + self._fwhm_spin.setToolTip( + "PSF full-width at half-maximum in pixels.\n" + "Match this to the star FWHM measured in your actual images\n" + "so the simulated blur matches a realistic observation.") + self._fwhm_spin.valueChanged.connect(self._schedule_preview) + form.addRow("FWHM", self._fwhm_spin) + + self._beta_spin = QDoubleSpinBox() + self._beta_spin.setRange(1.5, 10.0) + self._beta_spin.setValue(4.77) + self._beta_spin.setSingleStep(0.1) + self._beta_spin.setDecimals(2) + self._beta_spin.setEnabled(False) + self._beta_spin.setToolTip( + "Moffat β — controls the power-law slope of the PSF wings.\n" + "4.77 is the empirical mean for ground-based seeing (Trujillo et al. 2001).\n" + "Higher β → tighter wings (more Gaussian); lower → broader halos.") + self._beta_spin.valueChanged.connect(self._schedule_preview) + form.addRow("Moffat β", self._beta_spin) + + form.addRow(QLabel("")) # visual spacer + + # Noise section + self._noise_cb = QCheckBox("Apply noise") + self._noise_cb.setToolTip( + "Add Poisson sky shot noise and Gaussian read noise.\n" + "The clean reference is always noise-free so the effect of\n" + "noise on each metric can be measured independently of PSF blur.\n\n" + "When both PSF blur and noise are enabled, the degraded image\n" + "matches a realistic single-sub-frame observation.") + self._noise_cb.toggled.connect(self._on_noise_toggled) + self._noise_cb.toggled.connect(self._schedule_preview) + form.addRow(self._noise_cb) + + self._rn_spin = QDoubleSpinBox() + self._rn_spin.setRange(0.0, 500.0) + self._rn_spin.setValue(10.0) + self._rn_spin.setSuffix(" ADU") + self._rn_spin.setSingleStep(1.0) + self._rn_spin.setDecimals(1) + self._rn_spin.setEnabled(False) + self._rn_spin.setToolTip( + "Camera read noise σ in ADU (Gaussian component).\n" + "Sky Poisson noise is added automatically from the sky DC level.\n" + "Typical values: 2–20 ADU. At 500 ADU sky and 10 ADU read noise,\n" + "sky shot noise (√500 ≈ 22 ADU) dominates.") + self._rn_spin.valueChanged.connect(self._schedule_preview) + form.addRow("Read noise", self._rn_spin) + + return box + + def _on_psf_toggled(self, checked: bool) -> None: + self._fwhm_spin.setEnabled(checked) + self._beta_spin.setEnabled(checked) + + def _on_noise_toggled(self, checked: bool) -> None: + self._rn_spin.setEnabled(checked) + + def _build_output(self) -> QGroupBox: + box = QGroupBox("Output") + lay = QVBoxLayout(box) + + dir_row = QHBoxLayout() + self._outdir_lbl = QLabel(str(Path.home())) + self._outdir_lbl.setWordWrap(True) + self._outdir_lbl.setStyleSheet("font-size: 10px; color: #444;") + dir_btn = QPushButton("Browse…") + dir_btn.setToolTip("Choose the folder where generated FITS files will be saved.") + dir_btn.clicked.connect(self._browse_output) + dir_row.addWidget(self._outdir_lbl, stretch=1) + dir_row.addWidget(dir_btn) + lay.addLayout(dir_row) + + # Auto-load combo + self._load_combo = QComboBox() + for label, mode in _LOAD_MODES: + self._load_combo.addItem(label, mode) + self._load_combo.setCurrentIndex(0) # default: clean→A, degraded→B + self._load_combo.setToolTip( + "Which panel(s) to auto-load after generation.\n\n" + "The recommended workflow is Clean → A, Degraded → B so you can\n" + "run the analysis and see exactly how each metric responds to the\n" + "known degradation you configured.") + lay.addWidget(self._load_combo) + + lay.addStretch() + + self._gen_btn = QPushButton("Generate") + self._gen_btn.setDefault(True) + self._gen_btn.setFixedHeight(34) + self._gen_btn.setToolTip( + "Generate the clean reference and degraded companion FITS files.\n\n" + "Always produces two files:\n" + " • target_…_clean.fits — pure pattern, no PSF/noise\n" + " • target_…_degraded.fits — same pattern with configured degradation\n\n" + "Load both into Image A and B, then run the analysis to calibrate\n" + "your intuition for the spatial-detail metrics.") + self._gen_btn.clicked.connect(self._on_generate) + lay.addWidget(self._gen_btn) + + self._status_lbl = QLabel("Ready") + self._status_lbl.setAlignment(Qt.AlignmentFlag.AlignCenter) + self._status_lbl.setStyleSheet("font-size: 10px; color: #555;") + lay.addWidget(self._status_lbl) + + return box + + def _build_preview_section(self) -> QGroupBox: + box = QGroupBox("Live Preview (degraded image)") + lay = QHBoxLayout(box) + + left = QVBoxLayout() + self._preview_cb = QCheckBox("Enable live preview") + self._preview_cb.setToolTip( + "Regenerates a reduced-resolution preview of the degraded target\n" + "~400 ms after any parameter change.") + self._preview_cb.toggled.connect(self._on_preview_toggled) + left.addWidget(self._preview_cb) + left.addStretch() + lay.addLayout(left) + + self._preview_lbl = QLabel("Enable preview to see a live update") + self._preview_lbl.setAlignment(Qt.AlignmentFlag.AlignCenter) + self._preview_lbl.setStyleSheet("background: #1a1a1a; color: #666;") + self._preview_lbl.setMinimumSize(320, 160) + self._preview_lbl.setSizePolicy( + QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Preferred) + lay.addWidget(self._preview_lbl, stretch=1) + + return box + + # ------------------------------------------------------------------ + # Slots + # ------------------------------------------------------------------ + + def _show_info(self) -> None: + dlg = _ZoneInfoDialog(self) + dlg.exec() + + def _browse_output(self) -> None: + d = QFileDialog.getExistingDirectory( + self, "Select output directory", self._outdir_lbl.text()) + if d: + self._outdir_lbl.setText(d) + from PyQt6.QtCore import QSettings + QSettings("FilterImageComparator", "FilterImageComparator").setValue( + "target_output_dir", d) + + def _on_generate(self) -> None: + if self._gen_thread and self._gen_thread.isRunning(): + return + params = self._collect_params() + os.makedirs(params["output_dir"], exist_ok=True) + from PyQt6.QtCore import QSettings + QSettings("FilterImageComparator", "FilterImageComparator").setValue( + "target_output_dir", params["output_dir"]) + self._gen_btn.setEnabled(False) + self._status_lbl.setText("Generating…") + self._gen_thread = _TargetThread(params, self) + self._gen_thread.finished.connect(self._on_gen_done) + self._gen_thread.failed.connect(self._on_gen_failed) + self._gen_thread.start() + + @pyqtSlot(str, str) + def _on_gen_done(self, clean_path: str, degraded_path: str) -> None: + self._gen_btn.setEnabled(True) + self._status_lbl.setText(f"Saved: {Path(clean_path).name}") + mode = self._load_combo.currentData() + self.targets_generated.emit(clean_path, degraded_path, mode or "") + + @pyqtSlot(str) + def _on_gen_failed(self, msg: str) -> None: + self._gen_btn.setEnabled(True) + self._status_lbl.setText(f"Error: {msg}") + + def _on_preview_toggled(self, checked: bool) -> None: + if checked: + self._schedule_preview() + else: + self._preview_lbl.setText("Enable preview to see a live update") + self._preview_lbl.setPixmap(QPixmap()) + + def _schedule_preview(self) -> None: + if self._preview_cb.isChecked(): + self._preview_timer.start() + + def _run_preview(self) -> None: + if self._prev_thread and self._prev_thread.isRunning(): + return + params = self._collect_params() + self._prev_thread = _PreviewThread(params, self) + self._prev_thread.finished.connect(self._on_preview_done) + self._prev_thread.start() + + @pyqtSlot(object) + def _on_preview_done(self, arr: np.ndarray) -> None: + stretched = normalize_for_display(arr) + h, w = stretched.shape + qimg = QImage(stretched.data, w, h, w, QImage.Format.Format_Grayscale8) + pixmap = QPixmap.fromImage(qimg) + self._preview_lbl.setPixmap( + pixmap.scaled(self._preview_lbl.width(), self._preview_lbl.height(), + Qt.AspectRatioMode.KeepAspectRatio, + Qt.TransformationMode.SmoothTransformation)) + + # ------------------------------------------------------------------ + # Parameter collection + # ------------------------------------------------------------------ + + def _collect_params(self) -> dict: + idx = self._size_combo.currentIndex() + _, W, H = _SIZE_PRESETS[min(idx, len(_SIZE_PRESETS) - 1)] + return { + "width": W, + "height": H, + "contrast_min": self._contrast_min_sl.value(), + "contrast_max": self._contrast_max_sl.value(), + "sky_adu": self._sky_spin.value(), + "apply_psf": self._psf_cb.isChecked(), + "fwhm_px": self._fwhm_spin.value(), + "moffat_beta": self._beta_spin.value(), + "apply_noise": self._noise_cb.isChecked(), + "read_noise_adu": self._rn_spin.value(), + "seed": 42, + "output_dir": self._outdir_lbl.text(), + } diff --git a/report/report_builder.py b/report/report_builder.py index 276ba05..9101244 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -3624,6 +3624,18 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: f"{_val(va)}" f"{_val(vb)}") + # Michelson contrast table + mc_a = sm.get("michelson_contrast_a", {}) + mc_b = sm.get("michelson_contrast_b", {}) + mc_rows = "" + for ks in sorted(set(list(mc_a.keys()) + list(mc_b.keys()))): + va = mc_a.get(ks) + vb = mc_b.get(ks) + ca, cb = _better_worse_class(va, vb) + mc_rows += (f"{ks} px" + f"{_val(va, '.4f')}" + f"{_val(vb, '.4f')}") + # Wavelet SNR table snr_a = sm.get("wavelet_snr_a", {}) snr_b = sm.get("wavelet_snr_b", {}) @@ -3790,7 +3802,41 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: {paired_figs_for("wavelet_level", "xs_wavelet_level")}

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.

""" +The difference panel (right) shows where fine structure differs between the two filters.

+ +

8e. Michelson Contrast Maps

+{_info_box( + '

Formula: ' + 'C(x, y) = (Imax − Imin)' + ' / (Imax + Imin + ε), ' + 'where Imax and Imin are the maximum and minimum pixel values ' + 'within a K × K local neighbourhood. ' + 'Output is bounded in [0, 1]: 0 = flat region (no local contrast), ' + '1 = maximum local contrast.

' + '

Scalar metric (table below): The maximum Michelson value ' + 'within the analysis region. This captures the finest local contrast the image achieves ' + '— a higher maximum indicates finer structure survives at that scale.

' + '

Kernel sizes: Small kernels (3 px) are sensitive to ' + 'sub-pixel-scale contrast and noise. Medium kernels (5 px) reflect fine nebula ' + 'filaments and thread-like structure. Large kernels (9 px) capture coarser ' + 'structural contrast such as knots and shell edges.

' + '

Interpreting filter comparisons: ' + 'A narrowband filter preserving fine emission filaments will show higher Michelson ' + 'contrast in nebula regions at small kernel sizes. If both images show max values near ' + '1.0 at all scales, the ROI likely contains stars — select a star-free nebula region ' + 'for a more discriminating comparison. Maps use a starless image when one is available. ' + 'Colour maps are displayed on a linear [0, 1] scale; no gamma compression is applied ' + 'because the values are already bounded.

', + title="Michelson contrast")} + + + {mc_rows} +
Kernel size{ra.label} (max C){rb.label} (max C)
+{figs_for("michelson_")} +

Per-pixel Michelson contrast maps at each kernel size (linear colour +scale, viridis, shared across both panels). Brighter regions have higher local contrast at +that neighbourhood scale. The difference panel (A−B) highlights where one image +achieves greater local contrast. Values near 1.0 indicate strong brightness transitions.

""" # ── Section 9: Signal-to-Noise Ratio ───────────────────────────────────── @@ -4189,6 +4235,8 @@ def row(metric, val_a, val_b, fmt=".3f", cr_b = sm_b.get("contrast_ratios_b", {}) if sm_b else {} snr_wav_a = sm_a.get("wavelet_snr_a", {}) snr_wav_b = sm_b.get("wavelet_snr_b", {}) if sm_b else {} + mc_a_s = sm_a.get("michelson_contrast_a", {}) + mc_b_s = sm_b.get("michelson_contrast_b", {}) if sm_b else {} snr_ma = ra.snr_metrics or {} snr_mb = rb.snr_metrics or {} @@ -4246,6 +4294,7 @@ def row_pm(metric, val_a, val_b, spread_a, spread_b, fmt=".3f", row("Std contrast ratio (15px)", cr_a.get(15), cr_b.get(15)), row("Wavelet SNR level 2", snr_wav_a.get(2), snr_wav_b.get(2)), row("Wavelet SNR level 3", snr_wav_a.get(3), snr_wav_b.get(3)), + row("Michelson contrast max (5px)", mc_a_s.get(5), mc_b_s.get(5), fmt=".4f"), *([row( "Global SNR — starless (σ) ★", (snr_ma.get("starless") or {}).get("snr_global"), diff --git a/synthetic/target_generator.py b/synthetic/target_generator.py new file mode 100644 index 0000000..acb5777 --- /dev/null +++ b/synthetic/target_generator.py @@ -0,0 +1,282 @@ +"""Spatial test target generator for spatial-detail metric calibration. + +Generates a 4×3 grid of calibrated test patterns with known spatial frequencies +and waveform types. Load the clean reference into Image A and the degraded +companion into Image B to calibrate the app's spatial-detail metrics against +known inputs. + +The four column frequencies are deliberately chosen to align with the four +wavelet decomposition band centres (levels 4→1, coarse→fine): + f=0.04 c/px → λ=25 px → wavelet level 4 (~16 px scale) + f=0.08 c/px → λ=12 px → wavelet level 3 (~8 px scale) + f=0.16 c/px → λ=6 px → wavelet level 2 (~4 px scale) + f=0.32 c/px → λ=3 px → wavelet level 1 (~2 px / near Nyquist) +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from scipy.signal import fftconvolve + +from synthetic.generator import _moffat_2d, _fwhm_to_alpha, _siemens_star + + +# (row, col) → (waveform_type, frequency_cycles_per_pixel) +# Row 0: horizontal sine sweep — cleanest MTF measurement +# Row 1: horizontal square-wave sweep — tests harmonic/ringing behavior +# Row 2: direction variants + special targets +_ZONE_SPEC: dict[tuple[int, int], tuple[str, float]] = { + (0, 0): ("sine_h", 0.04), + (0, 1): ("sine_h", 0.08), + (0, 2): ("sine_h", 0.16), + (0, 3): ("sine_h", 0.32), + (1, 0): ("square_h", 0.04), + (1, 1): ("square_h", 0.08), + (1, 2): ("square_h", 0.16), + (1, 3): ("square_h", 0.32), + (2, 0): ("sine_v", 0.08), + (2, 1): ("sine_45", 0.08), + (2, 2): ("siemens", 0.00), + (2, 3): ("edge", 0.00), +} + +N_ROWS = 3 +N_COLS = 4 + +# Human-readable zone descriptions exported for use in the dialog info panel +ZONE_LABELS: dict[tuple[int, int], str] = { + (0, 0): "Sine H f=0.04 c/px (λ≈25 px, wavelet L4)", + (0, 1): "Sine H f=0.08 c/px (λ≈12 px, wavelet L3)", + (0, 2): "Sine H f=0.16 c/px (λ≈6 px, wavelet L2)", + (0, 3): "Sine H f=0.32 c/px (λ≈3 px, wavelet L1 / near Nyquist)", + (1, 0): "Square H f=0.04 c/px", + (1, 1): "Square H f=0.08 c/px", + (1, 2): "Square H f=0.16 c/px", + (1, 3): "Square H f=0.32 c/px", + (2, 0): "Sine V f=0.08 c/px (vertical direction)", + (2, 1): "Sine 45° f=0.08 c/px (diagonal direction)", + (2, 2): "Siemens star (all-angle radial chirp)", + (2, 3): "Slant edge ~5° (ISO 12233-style ESF input)", +} + + +class SpatialTargetGenerator: + """Generate a calibrated spatial test target FITS image pair. + + Always produces two FITS files: + • clean reference — pure pattern, no PSF, no noise + • degraded companion — same pattern with optional PSF blur and/or noise + Load clean → Image A and degraded → Image B, then run the analysis to see + how each metric responds to known spatial-frequency content. + """ + + def generate(self, params: dict, + preview: bool = False) -> tuple[str, str] | np.ndarray: + """Generate the target pair. + + preview=True → returns degraded ndarray (float32) for live display + preview=False → returns (clean_path, degraded_path) as FITS paths + """ + width = int(params["width"]) + height = int(params["height"]) + contrast_min = float(params.get("contrast_min", params.get("contrast", 0.02))) + contrast_max = float(params.get("contrast_max", params.get("contrast", 0.50))) + sky_adu = float(params["sky_adu"]) + fwhm_px = float(params.get("fwhm_px", 0.0)) + beta = float(params.get("moffat_beta", 4.77)) + rn_adu = float(params.get("read_noise_adu", 10.0)) + + if preview: + scale = min(640 / width, 480 / height, 1.0) + width = max(64, int(width * scale) & ~1) + height = max(64, int(height * scale) & ~1) + + clean = self._build_pattern(width, height, contrast_min, contrast_max, sky_adu) + degraded = clean.copy() + + if params.get("apply_psf") and fwhm_px > 0: + degraded = self._apply_psf(degraded, fwhm_px, beta) + if params.get("apply_noise"): + degraded = self._add_noise(degraded, sky_adu, rn_adu, + seed=int(params.get("seed", 42))) + + clean = np.clip(clean, 0.0, 65535.0).astype(np.float32) + degraded = np.clip(degraded, 0.0, 65535.0).astype(np.float32) + + if preview: + return degraded + + clean_path = self._write_fits(clean, params, is_clean=True) + degraded_path = self._write_fits(degraded, params, is_clean=False) + return (clean_path, degraded_path) + + # ------------------------------------------------------------------ + # Pattern construction + # ------------------------------------------------------------------ + + def _build_pattern(self, W: int, H: int, + contrast_min: float, contrast_max: float, + dc: float) -> np.ndarray: + """Assemble the N_ROWS × N_COLS zone grid as a float64 array.""" + image = np.full((H, W), dc, dtype=np.float64) + zw = W // N_COLS + zh = H // N_ROWS + amplitude_min = dc * contrast_min + amplitude_max = dc * contrast_max + + for (row, col), (waveform, freq) in _ZONE_SPEC.items(): + x0 = col * zw + y0 = row * zh + x1 = x0 + zw if col < N_COLS - 1 else W + y1 = y0 + zh if row < N_ROWS - 1 else H + zW, zH = x1 - x0, y1 - y0 + image[y0:y1, x0:x1] = self._make_zone( + waveform, freq, zW, zH, amplitude_min, amplitude_max, dc) + return image + + def _make_zone(self, waveform: str, freq: float, + W: int, H: int, + amplitude_min: float, amplitude_max: float, + dc: float) -> np.ndarray: + """Return an (H × W) float64 zone with contrast ramping top→bottom. + + Contrast (Michelson) increases linearly from amplitude_min/dc at the + top edge to amplitude_max/dc at the bottom edge of the zone. At any + horizontal row across all N_COLS zones the contrast is identical, + enabling direct cross-frequency comparison at the same contrast level. + """ + y_idx, x_idx = np.mgrid[0:H, 0:W].astype(np.float64) + # Y-ramp: 0 at top of zone, 1 at bottom — independent of absolute image position + y_rel = y_idx / max(H - 1, 1) + amplitude_map = amplitude_min + (amplitude_max - amplitude_min) * y_rel + + if waveform == "sine_h": + pattern = np.sin(2.0 * np.pi * freq * x_idx) + + elif waveform == "square_h": + # Hard square wave — preserve edge sharpness to test Gibbs behavior + pattern = np.where(np.sin(2.0 * np.pi * freq * x_idx) >= 0, 1.0, -1.0) + + elif waveform == "sine_v": + pattern = np.sin(2.0 * np.pi * freq * y_idx) + + elif waveform == "sine_45": + # Diagonal at 45°: spatial frequency along (1,1) direction + d = (x_idx + y_idx) / np.sqrt(2.0) + pattern = np.sin(2.0 * np.pi * freq * d) + + elif waveform == "siemens": + size = min(W, H) + if size % 2 == 0: + size -= 1 # odd so the star centre falls on a pixel + ss = _siemens_star(size) # values 0–~0.5, tapered to 0 at centre/edge + # Centre at 0 and normalise to [-1, 1] so the zone DC = dc + ss_c = ss - ss.mean() + ss_max = np.abs(ss_c).max() + if ss_max > 0: + ss_c /= ss_max + pat = np.zeros((H, W), dtype=np.float64) + sx = (W - size) // 2; sy = (H - size) // 2 + tx = max(0, sx); ty = max(0, sy) + ox = max(0, -sx); oy = max(0, -sy) + cw = min(size - ox, W - tx) + ch = min(size - oy, H - ty) + if cw > 0 and ch > 0: + pat[ty:ty+ch, tx:tx+cw] = ss_c[oy:oy+ch, ox:ox+cw] + pattern = pat + + elif waveform == "edge": + # Slant edge at ~5° from vertical — hard step for ISO 12233-style ESF + angle_rad = np.radians(5.0) + cx, cy = W / 2.0, H / 2.0 + # Signed distance from the slant line (positive = right side) + dist = ((x_idx - cx) * np.cos(angle_rad) + - (y_idx - cy) * np.sin(angle_rad)) + pattern = np.where(dist >= 0, 1.0, -1.0) + + else: + pattern = np.zeros((H, W), dtype=np.float64) + + return dc + amplitude_map * pattern + + # ------------------------------------------------------------------ + # Degradation + # ------------------------------------------------------------------ + + def _apply_psf(self, image: np.ndarray, + fwhm_px: float, beta: float) -> np.ndarray: + """Convolve with an isotropic Moffat PSF (field-centre, no aberrations).""" + alpha = _fwhm_to_alpha(fwhm_px, beta) + stamp = max(5, int(fwhm_px * 6)) + if stamp % 2 == 0: + stamp += 1 + kern = _moffat_2d(stamp, alpha, beta) + kern /= kern.sum() + result = fftconvolve(image, kern, mode="same") + return np.maximum(result, 0.0) + + def _add_noise(self, image: np.ndarray, sky_adu: float, + read_noise_adu: float, seed: int = 42) -> np.ndarray: + """Add Poisson sky shot noise + Gaussian read noise (gain=1 assumed).""" + rng = np.random.default_rng(seed) + sky_e = max(0.1, sky_adu) + poisson = rng.poisson(sky_e, size=image.shape).astype(np.float64) + rdnoise = rng.normal(0.0, read_noise_adu, size=image.shape) + return image + poisson + rdnoise + + # ------------------------------------------------------------------ + # FITS output + # ------------------------------------------------------------------ + + def _write_fits(self, image: np.ndarray, + params: dict, is_clean: bool) -> str: + from astropy.io import fits + + W = image.shape[1] + H = image.shape[0] + contrast_min = float(params.get("contrast_min", params.get("contrast", 0.02))) + contrast_max = float(params.get("contrast_max", params.get("contrast", 0.50))) + sky_adu = float(params["sky_adu"]) + + hdr = fits.Header() + hdr["SIMPLE"] = True + hdr["BITPIX"] = -32 + hdr["NAXIS"] = 2 + hdr["NAXIS1"] = W + hdr["NAXIS2"] = H + hdr["INSTRUME"] = "SpatialTarget" + hdr["EGAIN"] = (1.0, "e-/ADU gain=1 assumed") + hdr["GAIN"] = (1.0, "e-/ADU gain=1 assumed") + # Grid metadata — allows future auto-detection of zone layout in analysis + hdr["TGT_TYPE"] = ("spatial_grid", "generator: spatial test target") + hdr["TGT_ROWS"] = (N_ROWS, "grid rows") + hdr["TGT_COLS"] = (N_COLS, "grid columns") + hdr["TGT_CMIN"] = (contrast_min, "Michelson contrast at zone top (min)") + hdr["TGT_CMAX"] = (contrast_max, "Michelson contrast at zone bottom (max)") + hdr["TGT_SKY"] = (sky_adu, "DC sky level ADU") + hdr["TGT_CLEN"] = (is_clean, "True=clean reference image") + if params.get("apply_psf") and not is_clean: + hdr["TGT_FWHM"] = (float(params.get("fwhm_px", 0.0)), "PSF FWHM px") + hdr["TGT_BETA"] = (float(params.get("moffat_beta", 4.77)), "Moffat beta") + if params.get("apply_noise") and not is_clean: + hdr["TGT_RN"] = (float(params.get("read_noise_adu", 10.0)), "read noise ADU") + hdr["TGT_NSED"] = (int(params.get("seed", 42)), "noise RNG seed") + # Per-zone frequency and waveform type for future per-zone analysis + for (row, col), (waveform, freq) in _ZONE_SPEC.items(): + hdr[f"TGT_{row}{col}F"] = (freq, f"r{row}c{col} freq c/px") + hdr[f"TGT_{row}{col}W"] = (waveform[:8], f"r{row}c{col} waveform") + + hdu = fits.PrimaryHDU(image, header=hdr) + hdul = fits.HDUList([hdu]) + + tag = "clean" if is_clean else "degraded" + stem = f"target_Cmin{contrast_min:.2f}_Cmax{contrast_max:.2f}_SKY{sky_adu:.0f}" + if params.get("apply_psf") and not is_clean: + stem += f"_PSF{params.get('fwhm_px', 0.0):.1f}px" + if params.get("apply_noise") and not is_clean: + stem += f"_RN{params.get('read_noise_adu', 10.0):.0f}" + out_path = Path(params["output_dir"]) / f"{stem}_{tag}.fits" + hdul.writeto(str(out_path), overwrite=True) + return str(out_path) From 297f9582e3dc77f391d7de0b2d1960b0980ff6db Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 20:13:13 -0400 Subject: [PATCH 05/14] feat: replace Michelson contrast with Weber fraction contrast in Section 8e MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Michelson contrast (C = ΔL/(I_max+I_min)) proved less useful in practice. Weber fraction contrast (c = ΔL/L where L = kernel median) better represents how much brighter-than-background a feature appears, matching Weber's Law. The median denominator is robust to bright filaments and residual star halos that would inflate the mean. Scalar metric uses the 99th percentile of the Weber map (not max) to avoid dark-sky pixels with near-zero median driving the comparison value to extremes. Maps use PowerNorm for display since the output is unbounded. Co-Authored-By: Claude Sonnet 4.6 --- analysis/image_filters.py | 104 +++++++++++++++++++------------------- core/models.py | 2 +- report/report_builder.py | 82 +++++++++++++++--------------- 3 files changed, 95 insertions(+), 93 deletions(-) diff --git a/analysis/image_filters.py b/analysis/image_filters.py index 2162b4d..0c438bc 100644 --- a/analysis/image_filters.py +++ b/analysis/image_filters.py @@ -10,13 +10,13 @@ import matplotlib.colors as mcolors matplotlib.use("Agg") import matplotlib.pyplot as plt -from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, map_coordinates, zoom, maximum_filter, minimum_filter +from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter import pywt from core.astro_image import AstroImage from core.fig_utils import fig_to_b64, figs_to_b64 from core.models import (STD_KERNEL_SIZES, LOG_SIGMAS, WAVELET_NAME, WAVELET_LEVELS, - MICHELSON_KERNEL_SIZES, + WEBER_KERNEL_SIZES, XS_LINE_ALPHA, SECTION8_BORDER_CROP_FRACTION, SECTION8_ANALYSIS_CMAP, XS_SNR_REGION_WIDTH) @@ -36,7 +36,7 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, log_sigmas: tuple = LOG_SIGMAS, wavelet: str = WAVELET_NAME, levels: int = WAVELET_LEVELS, - michelson_kernel_sizes: tuple = MICHELSON_KERNEL_SIZES, + weber_kernel_sizes: tuple = WEBER_KERNEL_SIZES, crosshair: dict | None = None, roi: tuple | None = None, xs_snr_width: int | None = None) -> dict: @@ -55,8 +55,8 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, "wavelet_snr_b": {}, "sigma_noise_a": None, "sigma_noise_b": None, - "michelson_contrast_a": {}, - "michelson_contrast_b": {}, + "weber_contrast_a": {}, + "weber_contrast_b": {}, "panels": {}, } figures: dict = {} @@ -135,33 +135,33 @@ def _clip01(v): return max(0.0, min(1.0, v)) display_roi=display_roi, crosshair=crosshair_roi, ) - _f_mic = _ex.submit(self._michelson_analysis, + _f_web = _ex.submit(self._weber_analysis, analysis_a, analysis_b, - michelson_kernel_sizes, + weber_kernel_sizes, image_a.label, _label_b, display_roi=display_roi, ) std_b64, std_partial = _f_std.result() log_b64, log_partial = _f_log.result() wav_b64, wav_partial = _f_wav.result() - mic_b64, mic_partial = _f_mic.result() + web_b64, web_partial = _f_web.result() figures.update(std_b64) figures.update(log_b64) figures.update(wav_b64) - figures.update(mic_b64) + figures.update(web_b64) result["contrast_ratios_a"].update(std_partial["contrast_ratios_a"]) result["contrast_ratios_b"].update(std_partial["contrast_ratios_b"]) result["sigma_noise_a"] = wav_partial["sigma_noise_a"] result["sigma_noise_b"] = wav_partial["sigma_noise_b"] result["wavelet_snr_a"].update(wav_partial["wavelet_snr_a"]) result["wavelet_snr_b"].update(wav_partial["wavelet_snr_b"]) - result["michelson_contrast_a"].update(mic_partial["michelson_contrast_a"]) - result["michelson_contrast_b"].update(mic_partial["michelson_contrast_b"]) + result["weber_contrast_a"].update(web_partial["weber_contrast_a"]) + result["weber_contrast_b"].update(web_partial["weber_contrast_b"]) result["panels"].update(std_partial["panels"]) result["panels"].update(log_partial["panels"]) result["panels"].update(wav_partial["panels"]) - result["panels"].update(mic_partial["panels"]) + result["panels"].update(web_partial["panels"]) if crosshair is not None: pos_a, prof_a = self._sample_line(norm_a, **crosshair) @@ -506,60 +506,62 @@ def _reconstruct_level(self, coeffs, target_coeff_idx: int, # ------------------------------------------------------------------ # ------------------------------------------------------------------ - # Michelson contrast maps + # Weber fraction contrast maps # ------------------------------------------------------------------ - def _michelson_analysis(self, norm_a, norm_b, kernel_sizes, - label_a, label_b, - display_roi=None) -> tuple[dict, dict]: + def _weber_analysis(self, norm_a, norm_b, kernel_sizes, + label_a, label_b, + display_roi=None) -> tuple[dict, dict]: figures = {} partial: dict = { - "michelson_contrast_a": {}, - "michelson_contrast_b": {}, + "weber_contrast_a": {}, + "weber_contrast_b": {}, "panels": {}, } single = norm_b is None for ks in kernel_sizes: - mc_a = self._compute_michelson_map(norm_a, ks) - mc_b = self._compute_michelson_map(norm_b, ks) if not single else None + wc_a = self._compute_weber_map(norm_a, ks) + wc_b = self._compute_weber_map(norm_b, ks) if not single else None - partial["michelson_contrast_a"][ks] = float(np.max(mc_a)) + # 99th percentile avoids dark-sky floor driving the scalar metric to extremes + partial["weber_contrast_a"][ks] = float(np.percentile(wc_a, 99)) if not single: - partial["michelson_contrast_b"][ks] = float(np.max(mc_b)) + partial["weber_contrast_b"][ks] = float(np.percentile(wc_b, 99)) - partial["panels"][f"michelson_{ks}px"] = { - "a": mc_a.astype(np.float32), - "b": mc_b.astype(np.float32) if mc_b is not None else None, - "diff": (mc_a - mc_b).astype(np.float32) if mc_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, } if not single: fig = self._plot_side_by_side( - self._crop_border(mc_a, SECTION8_BORDER_CROP_FRACTION), - self._crop_border(mc_b, SECTION8_BORDER_CROP_FRACTION), - f"Michelson contrast — kernel {ks}px — {label_a}", - f"Michelson contrast — kernel {ks}px — {label_b}", - diff_title=f"Michelson diff (A−B), kernel {ks}px", + 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", cmap=SECTION8_ANALYSIS_CMAP, - nonlinear_norm=False, # values bounded [0,1]; linear scale is appropriate - display_roi=None, # _crop_border already applied; matches _log_analysis pattern + nonlinear_norm=True, # unbounded output; PowerNorm compresses dynamic range + display_roi=None, # _crop_border already applied ) else: fig = self._plot_single( - self._crop_border(mc_a, SECTION8_BORDER_CROP_FRACTION), - f"Michelson contrast — kernel {ks}px — {label_a}", + self._crop_border(wc_a, SECTION8_BORDER_CROP_FRACTION), + f"Weber contrast — kernel {ks}px — {label_a}", cmap=SECTION8_ANALYSIS_CMAP, - nonlinear_norm=False, + nonlinear_norm=True, ) - figures[f"michelson_{ks}px"] = fig + figures[f"weber_{ks}px"] = fig return figs_to_b64(figures, dpi=150), partial - def _compute_michelson_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray: - """Per-pixel Michelson contrast C = (I_max - I_min) / (I_max + I_min + eps). - Downsamples large images identically to _compute_std_map for performance. + def _compute_weber_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray: + """Weber fraction contrast c = ΔL / L where ΔL = max−min and L = median of kernel. + L uses median (not mean) — robust background luminance unaffected by bright filaments. + Output is unbounded >= 0; nonlinear_norm=True is used for display. """ - _EPS = 1e-9 + _EPS = 1e-6 # L (median) can be very small over dark sky; larger EPS prevents extremes factor = 1.0 data = norm if max(norm.shape) > MAX_DIM_FOR_STD: @@ -571,18 +573,18 @@ def _compute_michelson_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarr i_max = maximum_filter(data, size=kernel_size) i_min = minimum_filter(data, size=kernel_size) - # Background-subtracted images can have negative pixels after mean normalisation. - # Clip to 0 so i_max + i_min is never negative, preventing spurious near-1 contrast. - i_max = np.maximum(i_max, 0.0) - i_min = np.maximum(i_min, 0.0) - mc_map = (i_max - i_min) / (i_max + i_min + _EPS) + i_med = median_filter(data, size=kernel_size) + + delta_L = i_max - i_min # always >= 0 + L = np.maximum(i_med, 0.0) # bg-subtracted images can have negative medians + weber = delta_L / (L + _EPS) if factor < 1.0: - mc_map = zoom(mc_map, - (norm.shape[0] / mc_map.shape[0], - norm.shape[1] / mc_map.shape[1]), - order=1) - return np.clip(mc_map, 0.0, 1.0) + weber = zoom(weber, + (norm.shape[0] / weber.shape[0], + norm.shape[1] / weber.shape[1]), + order=1) + return np.maximum(weber, 0.0) def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray, title_a: str, title_b: str, diff --git a/core/models.py b/core/models.py index dce7a52..cece8df 100644 --- a/core/models.py +++ b/core/models.py @@ -30,7 +30,7 @@ STD_KERNEL_SIZES = (3, 5, 10) #originally (5, 10, 15) # px; Gaussian kernel sizes for std dev maps LOG_SIGMAS = (1.5, 3.0, 6.0) -MICHELSON_KERNEL_SIZES = (3, 5, 9) # px; local max/min neighbourhood for Michelson contrast (odd values required) +WEBER_KERNEL_SIZES = (3, 5, 9) # px; local kernel for Weber fraction contrast c = ΔL/L (odd values required) WAVELET_NAME = "db4" WAVELET_LEVELS = 4 diff --git a/report/report_builder.py b/report/report_builder.py index 9101244..9fb5636 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -3624,15 +3624,15 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: f"{_val(va)}" f"{_val(vb)}") - # Michelson contrast table - mc_a = sm.get("michelson_contrast_a", {}) - mc_b = sm.get("michelson_contrast_b", {}) - mc_rows = "" - for ks in sorted(set(list(mc_a.keys()) + list(mc_b.keys()))): - va = mc_a.get(ks) - vb = mc_b.get(ks) + # Weber fraction contrast table + wc_a = sm.get("weber_contrast_a", {}) + wc_b = sm.get("weber_contrast_b", {}) + wc_rows = "" + for ks in sorted(set(list(wc_a.keys()) + list(wc_b.keys()))): + va = wc_a.get(ks) + vb = wc_b.get(ks) ca, cb = _better_worse_class(va, vb) - mc_rows += (f"{ks} px" + wc_rows += (f"{ks} px" f"{_val(va, '.4f')}" f"{_val(vb, '.4f')}") @@ -3804,39 +3804,39 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: diverging colourmap), each followed by its cross-section profile. The difference panel (right) shows where fine structure differs between the two filters.

-

8e. Michelson Contrast Maps

+

8e. Weber Fraction Contrast Maps

{_info_box( - '

Formula: ' - 'C(x, y) = (Imax − Imin)' - ' / (Imax + Imin + ε), ' - 'where Imax and Imin are the maximum and minimum pixel values ' - 'within a K × K local neighbourhood. ' - 'Output is bounded in [0, 1]: 0 = flat region (no local contrast), ' - '1 = maximum local contrast.

' - '

Scalar metric (table below): The maximum Michelson value ' - 'within the analysis region. This captures the finest local contrast the image achieves ' - '— a higher maximum indicates finer structure survives at that scale.

' - '

Kernel sizes: Small kernels (3 px) are sensitive to ' - 'sub-pixel-scale contrast and noise. Medium kernels (5 px) reflect fine nebula ' - 'filaments and thread-like structure. Large kernels (9 px) capture coarser ' - 'structural contrast such as knots and shell edges.

' - '

Interpreting filter comparisons: ' - 'A narrowband filter preserving fine emission filaments will show higher Michelson ' - 'contrast in nebula regions at small kernel sizes. If both images show max values near ' - '1.0 at all scales, the ROI likely contains stars — select a star-free nebula region ' - 'for a more discriminating comparison. Maps use a starless image when one is available. ' - 'Colour maps are displayed on a linear [0, 1] scale; no gamma compression is applied ' - 'because the values are already bounded.

', - title="Michelson contrast")} + '

Formula: c = ΔL / L, ' + 'where ΔL = Imax − Imin (local range in the K × K kernel) ' + 'and L = median(kernel) (local background luminance). ' + 'Output is unbounded ≥ 0; a value of 1.0 means the local range equals the background ' + 'luminance, 2.0 means twice, and so on.

' + '

Why median for L: The median represents the background luminance ' + 'the feature is seen against, matching Weber\'s Law. Using the mean would inflate L ' + 'toward bright filaments within the kernel, artificially suppressing contrast values. ' + 'Median is also robust to hot pixels and residual star halos in starless images.

' + '

Scalar metric (table below): 99th percentile of the Weber map ' + 'within the analysis region. Near-zero median pixels over dark sky produce very large ' + 'Weber values; the 99th percentile captures peak structural contrast while ignoring ' + 'isolated dark-floor artefacts.

' + '

Kernel sizes: Small kernels (3 px) respond to sub-pixel-scale ' + 'transitions. Medium kernels (5 px) capture fine filaments. ' + 'Large kernels (9 px) reflect coarser structural contrast such as knots and shell edges.

' + '

Wide dynamic range: Weber contrast is intentionally unbounded. ' + 'Maps are displayed with a square-root colour scale (PowerNorm γ = 0.5) to compress ' + 'the bright end. Selecting a star-free nebula ROI avoids dark-sky pixels that drive ' + 'Weber values very high. Maps use a starless image when one is available.

', + title="Weber fraction contrast")} - - {mc_rows} + + {wc_rows}
Kernel size{ra.label} (max C){rb.label} (max C)
Kernel size{ra.label} (99th pct c){rb.label} (99th pct c)
-{figs_for("michelson_")} -

Per-pixel Michelson contrast maps at each kernel size (linear colour -scale, viridis, shared across both panels). Brighter regions have higher local contrast at -that neighbourhood scale. The difference panel (A−B) highlights where one image -achieves greater local contrast. Values near 1.0 indicate strong brightness transitions.

""" +{figs_for("weber_")} +

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.

""" # ── Section 9: Signal-to-Noise Ratio ───────────────────────────────────── @@ -4235,8 +4235,8 @@ def row(metric, val_a, val_b, fmt=".3f", cr_b = sm_b.get("contrast_ratios_b", {}) if sm_b else {} snr_wav_a = sm_a.get("wavelet_snr_a", {}) snr_wav_b = sm_b.get("wavelet_snr_b", {}) if sm_b else {} - mc_a_s = sm_a.get("michelson_contrast_a", {}) - mc_b_s = sm_b.get("michelson_contrast_b", {}) if sm_b else {} + wc_a_s = sm_a.get("weber_contrast_a", {}) + wc_b_s = sm_b.get("weber_contrast_b", {}) if sm_b else {} snr_ma = ra.snr_metrics or {} snr_mb = rb.snr_metrics or {} @@ -4294,7 +4294,7 @@ def row_pm(metric, val_a, val_b, spread_a, spread_b, fmt=".3f", row("Std contrast ratio (15px)", cr_a.get(15), cr_b.get(15)), row("Wavelet SNR level 2", snr_wav_a.get(2), snr_wav_b.get(2)), row("Wavelet SNR level 3", snr_wav_a.get(3), snr_wav_b.get(3)), - row("Michelson contrast max (5px)", mc_a_s.get(5), mc_b_s.get(5), fmt=".4f"), + row("Weber contrast 99th pct (5px)", wc_a_s.get(5), wc_b_s.get(5), fmt=".4f"), *([row( "Global SNR — starless (σ) ★", (snr_ma.get("starless") or {}).get("snr_global"), From 83c12c90d304b0324dfada510ec44c55d7538573 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 20:21:48 -0400 Subject: [PATCH 06/14] unit test update for weber contrast --- CLAUDE.md | 2 +- tests/test_analysis/test_spatial_detail.py | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/CLAUDE.md b/CLAUDE.md index 4d2c752..dbbf70b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -259,7 +259,7 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | Inline FITS too small to load | `_load_fits` skips HDUs where `max(shape) <= 100`. Test FITS must be at least 101×101; use 128×128 for safety. | | `float(mtf(array, m))` raises TypeError | `mtf()` returns a same-shape array, not a scalar. Index with `[0]` or pass a scalar input. | | `PSFAnalyzer.analyze()["figures"]` KeyError | `figures` is only added when `n_stars_used > 0`. Guard with `if result["n_stars_used"] > 0`. | -| `contrast_ratios_b` always present | `SpatialDetailAnalyzer.analyze()` always includes `contrast_ratios_b: {}` even in single-image mode. It is never `None` or absent — check `not b_ratios` instead. | +| `contrast_ratios_b` / `weber_contrast_b` always present | `SpatialDetailAnalyzer.analyze()` always includes `contrast_ratios_b: {}` and `weber_contrast_b: {}` even in single-image mode. Neither is ever `None` or absent — check `not b_ratios` / `not wc_b` instead. | | Background2D fails on tiny images | `estimate_background()` with default `box_size=64` needs the image to be larger than the box. Any image used in analysis tests should be at least 128×128. | --- diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py index 5986f8f..062844f 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -55,6 +55,21 @@ def test_minimal_image_no_crash(self, tmp_path): result = SpatialDetailAnalyzer().analyze(img) assert isinstance(result, dict) + def test_weber_contrast_a_present(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert "weber_contrast_a" in result + + def test_single_image_weber_contrast_b_empty(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + # Single-image mode: weber_contrast_b is present but empty (same pattern as contrast_ratios_b) + wc_b = result.get("weber_contrast_b") + assert wc_b is None or not wc_b + + def test_weber_contrast_a_positive(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + for v in result.get("weber_contrast_a", {}).values(): + assert v >= 0.0 + def test_with_roi(self, astro_image_a): result = SpatialDetailAnalyzer().analyze(astro_image_a, roi=(50, 50, 450, 450)) From 1c900f37251b134ad755ac3d9ef893e52589c26d Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 20:43:31 -0400 Subject: [PATCH 07/14] perf: switch working dtype from float64 to float32 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Change the central load-time conversion in core/astro_image.py from float64 to float32. All FITS/XISF/TIFF input data is at most 16-bit integer before stacking; float32 (24-bit mantissa) represents every possible value exactly while halving memory footprint and yielding 1.5-2x faster element-wise operations through better cache utilisation and wider SIMD lanes. Byte-order normalisation (big-endian FITS >f4 from astropy) remains a free side-effect of the astype() call. Cascading cleanup removes or narrows nine redundant float64 casts in snr_analyzer.py, edge_analyzer.py, and report_builder.py that were no-ops when self.data was already float64. The astroalign registration call in gui/analysis_thread.py retains its explicit float64 cast — astroalign requires it and the cast now does real conversion work. All 205 fast tests pass after the change. Co-Authored-By: Claude Sonnet 4.6 --- analysis/edge_analyzer.py | 6 +++--- analysis/snr_analyzer.py | 18 +++++++++--------- core/astro_image.py | 6 +++--- report/report_builder.py | 8 ++++---- tests/test_core/test_astro_image.py | 6 +++--- 5 files changed, 22 insertions(+), 22 deletions(-) diff --git a/analysis/edge_analyzer.py b/analysis/edge_analyzer.py index af94c0c..8e6b0ec 100644 --- a/analysis/edge_analyzer.py +++ b/analysis/edge_analyzer.py @@ -133,7 +133,7 @@ def analyze(self, image: AstroImage, from core.stretch import stf_stretch sigma = _gradient_sigma(image.pixel_scale) gm_full = gaussian_gradient_magnitude( - stf_stretch(bgsub).astype(np.float64), sigma=sigma + stf_stretch(bgsub), sigma=sigma ) gradient_fig = self._plot_gradient_map(gm_full, image.label, rois_used) @@ -178,7 +178,7 @@ def _auto_detect_top_rois(self, bgsub: np.ndarray, image: AstroImage, ) -> list[tuple[np.ndarray, tuple]]: """Find n well-separated patches centred on the strongest gradients.""" from core.stretch import stf_stretch - stretched = stf_stretch(bgsub).astype(np.float64) + stretched = stf_stretch(bgsub) sigma = _gradient_sigma(image.pixel_scale) gm = gaussian_gradient_magnitude(stretched, sigma=sigma) @@ -431,7 +431,7 @@ def analyze_crosshair(self, image: AstroImage, crosshair: dict) -> dict | None: n = max(20, int(length)) cols = np.linspace(c0, c1, n) rows = np.linspace(r0, r1, n) - profile = map_coordinates(bgsub.astype(np.float64), [rows, cols], + profile = map_coordinates(bgsub.astype(np.float32), [rows, cols], order=1, mode="nearest") positions = np.linspace(0.0, length, n) lo, hi = profile.min(), profile.max() diff --git a/analysis/snr_analyzer.py b/analysis/snr_analyzer.py index a55ab57..1a2d79a 100644 --- a/analysis/snr_analyzer.py +++ b/analysis/snr_analyzer.py @@ -19,7 +19,7 @@ class SNRAnalyzer: def analyze(self, image: AstroImage) -> dict: image.estimate_background() - bgsub = image.background_subtracted().astype(np.float64) + bgsub = image.background_subtracted() rms = image.background_rms # 2D array from Background2D noise = float(np.median(rms)) if rms is not None else 1.0 @@ -45,7 +45,7 @@ def analyze(self, image: AstroImage) -> dict: star_snr_median: float | None = None star_snr_iqr: float | None = None if cat is not None and len(cat) > 0 and noise > 0: - peaks = np.asarray(cat["peak"], dtype=np.float64) + peaks = np.asarray(cat["peak"], dtype=np.float32) star_snrs = peaks / noise star_snrs = star_snrs[np.isfinite(star_snrs) & (star_snrs > 0)] if star_snrs.size > 0: @@ -56,7 +56,7 @@ def analyze(self, image: AstroImage) -> dict: # --- Local SNR map ----------------------------------------------- if rms is not None: - denom = np.where(rms > 0, rms.astype(np.float64), np.nan) + denom = np.where(rms > 0, rms.astype(np.float32), np.nan) else: denom = noise snr_map = bgsub / denom @@ -129,12 +129,12 @@ def analyze(self, image: AstroImage) -> dict: "snr_global_db": snr_global_db, "noise_median": noise, # median sky RMS (σ_sky) in ADU "background_median": background_median, # median sky background (μ_sky) in ADU - "star_snr_median": star_snr_median, - "star_snr_iqr": star_snr_iqr, - "pct_above_3": pcts[3], - "pct_above_5": pcts[5], - "pct_above_10": pcts[10], - "pct_above_20": pcts[20], + "star_snr_median": star_snr_median, + "star_snr_iqr": star_snr_iqr, + "pct_above_3": pcts[3], + "pct_above_5": pcts[5], + "pct_above_10": pcts[10], + "pct_above_20": pcts[20], "snr_display": snr_display, "snr_p2": snr_p2, "snr_p98": snr_p98, diff --git a/core/astro_image.py b/core/astro_image.py index 3b3b342..add122c 100644 --- a/core/astro_image.py +++ b/core/astro_image.py @@ -42,7 +42,7 @@ def __init__(self, path: str, label: str = ""): self.pixel_scale_is_estimated: bool = False self.bandwidth_nm: float | None = None self.filter_thickness_mm: float = FILTER_THICKNESS_MM - self.original_dtype: np.dtype | None = None # dtype before float64 conversion + self.original_dtype: np.dtype | None = None # dtype before float32 conversion self.background: Background2D | None = None self.background_rms: np.ndarray | None = None self._load_error: str | None = None @@ -68,8 +68,8 @@ def load(self) -> None: raise ValueError(f"Unsupported file format: {suffix}") if self.data is not None: - self.original_dtype = self.data.dtype # capture before float64 conversion - self.data = self.data.astype(np.float64) + self.original_dtype = self.data.dtype # capture before float32 conversion + self.data = self.data.astype(np.float32) self.pixel_scale = self._extract_pixel_scale() self.bandwidth_nm = self._extract_bandwidth() self._extract_metadata() diff --git a/report/report_builder.py b/report/report_builder.py index 9fb5636..ee25b25 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -649,8 +649,8 @@ def _add_options_entry(section: str, name: str, options: dict) -> None: pa = result_a.psf_metrics or {} pb = result_b.psf_metrics or {} if pa.get("epsf_data") is not None and pb.get("epsf_data") is not None: - ea = pa["epsf_data"].astype(np.float64) - eb = pb["epsf_data"].astype(np.float64) + ea = pa["epsf_data"].astype(np.float32) + eb = pb["epsf_data"].astype(np.float32) _add("epsf_a", np.log1p(ea - ea.min()).astype(np.float32)) _add("epsf_b", np.log1p(eb - eb.min()).astype(np.float32)) ref_fwhm = _ref_fwhm_px(pa, pb, self._ref_seeing_arcsec) @@ -2786,7 +2786,7 @@ def _extract_cutout(self, data: np.ndarray, xc: float, yc: float, half: int) -> np.ndarray: h, w = data.shape size = 2 * half + 1 - cut = np.zeros((size, size), dtype=np.float64) + cut = np.zeros((size, size), dtype=np.float32) x0_src = max(0, int(xc) - half) x1_src = min(w, int(xc) + half + 1) y0_src = max(0, int(yc) - half) @@ -2849,7 +2849,7 @@ def _plot_halo_star_map(self, matched: list, img_a: AstroImage, bgsub = img_a.background_subtracted() if img_a.background is not None else img_a.data from core.stretch import stf_stretch - display = stf_stretch(bgsub).astype(np.float64) + display = stf_stretch(bgsub) h, w = display.shape if max(h, w) > 1200: diff --git a/tests/test_core/test_astro_image.py b/tests/test_core/test_astro_image.py index 33bf633..867f67d 100644 --- a/tests/test_core/test_astro_image.py +++ b/tests/test_core/test_astro_image.py @@ -10,8 +10,8 @@ class TestLoad: - def test_data_float64(self, astro_image_a): - assert astro_image_a.data.dtype == np.float64 + def test_data_float32(self, astro_image_a): + assert astro_image_a.data.dtype == np.float32 def test_correct_shape(self, astro_image_a): assert astro_image_a.data.shape == (512, 512) @@ -109,7 +109,7 @@ def test_background_subtracted_shape(self, astro_image_a): def test_background_subtracted_dtype(self, astro_image_a): bs = astro_image_a.background_subtracted() - assert bs.dtype == np.float64 + assert bs.dtype == np.float32 def test_background_subtracted_without_background(self, test_fits_path): img = AstroImage(str(test_fits_path), label="NoBkg") From f21e77555cf8fa3ecfa7d6692a5bbf64f2e09977 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 20:45:49 -0400 Subject: [PATCH 08/14] claude update --- CLAUDE.md | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index dbbf70b..ec30747 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -153,9 +153,33 @@ O(N² log N) and handles any kernel size without performance degradation: ```python from scipy.signal import fftconvolve -convolved = fftconvolve(patch, psf_kernel, mode="same").astype(np.float64) +convolved = fftconvolve(patch, psf_kernel, mode="same").astype(np.float32) ``` +### Working dtype — always float32 + +All image data is converted to `np.float32` at load time (`core/astro_image.py:72`). +This is the single working dtype throughout the pipeline — `self.data`, `background_subtracted()`, +background maps from photutils, and all intermediate analysis arrays. + +**float32 is sufficient:** source data is at most 16-bit integer before stacking; float32 +(24-bit mantissa, ~7 significant digits) represents every possible value exactly. The +switch from float64 halves memory footprint and yields ~1.5–2× faster element-wise +operations through better cache utilisation and wider SIMD lanes. + +**Byte order is handled automatically.** FITS `BITPIX=-32` images arrive from astropy as +big-endian `>f4`; `astype(np.float32)` always produces native-endian output, so there is +no need for `.byteswap()` or `.newbyteorder()`. + +**Do not add float64 casts in analysis code.** The only legitimate exception in the entire +codebase is the `astroalign` registration call in `gui/analysis_thread.py`, which requires +float64 internally. That explicit cast is already in place and must stay. + +**Synthetic generator internals stay float64.** `synthetic/generator.py` and +`synthetic/target_generator.py` accumulate many PSF stamps with `+=` across hundreds of +operations; float64 prevents rounding drift during synthesis. Both generators cast their +output to float32 before writing to FITS. + ### Large-array reductions — prefer bottleneck `bottleneck` (conda-forge) provides drop-in replacements for the numpy NaN-aware @@ -233,7 +257,7 @@ sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0 ```bash conda activate astrolab pip install pytest pytest-cov pytest-timeout # one-time setup; not in environment.yml -pytest tests/ -m "not slow" # fast suite (~90 s, 202 tests) +pytest tests/ -m "not slow" # fast suite (~120 s, 205 tests) pytest tests/ -m slow # slow/integration tests (full FITS generation) pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html ``` @@ -284,6 +308,8 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | Linux build needs system Qt libraries | PyInstaller must be able to import PyQt6 during analysis. On `ubuntu-latest` run `sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0` before `pip install -r requirements-build.txt`. | | `PowerSpectrumAnalyzer` crashes on images smaller than 1024 px | `POWER_SPECTRUM_NPIX = 1024`. The auto-select loop is empty when `min(h, w) < 1024`; the fallback produces negative slice indices → non-square region → `_apply_window` shape mismatch. Fix: `N = min(N, h, w)` before the loop, add `+1` to loop upper bounds, clamp fallback with `max(0, ...)`. | | `sigma_clip` mask is scalar `False` when nothing is clipped | `clipped.mask` is `np.ma.nomask` (== `False`) when no values are clipped. `region[False]` silently writes only the first row. Use `np.ma.getmaskarray(clipped)` to get a full bool array, then guard with `.any()`. | +| 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. | --- From 252a46dcdca20b791b8eac5247dd533495286def Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 21:05:32 -0400 Subject: [PATCH 09/14] synthetic data menu name updates --- gui/main_window.py | 11 ++++++----- gui/spatial_target_dialog.py | 2 +- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/gui/main_window.py b/gui/main_window.py index e57fd6d..7f653f9 100644 --- a/gui/main_window.py +++ b/gui/main_window.py @@ -104,18 +104,19 @@ def _build_menu(self) -> None: analysis_menu.addAction(act_run) tools_menu = mb.addMenu("&Tools") - act_synth = QAction("Synthetic &Data…", self) + act_synth = QAction("Synthetic Star &Data…", self) act_synth.triggered.connect(self._open_synthetic_dialog) tools_menu.addAction(act_synth) + + act_target = QAction("Synthetic Spatial Detail &Target…", self) + act_target.triggered.connect(self._open_target_dialog) + tools_menu.addAction(act_target) + act_halo = QAction("&Halo Analyzer…", self) act_halo.triggered.connect(self._open_halo_dialog) tools_menu.addAction(act_halo) - act_target = QAction("Spatial &Target…", self) - act_target.triggered.connect(self._open_target_dialog) - tools_menu.addAction(act_target) - help_menu = mb.addMenu("&Help") act_quickstart = QAction("&Quick Start Guide", self) act_quickstart.triggered.connect(self._open_quickstart) diff --git a/gui/spatial_target_dialog.py b/gui/spatial_target_dialog.py index b6f1d38..0c6103f 100644 --- a/gui/spatial_target_dialog.py +++ b/gui/spatial_target_dialog.py @@ -167,7 +167,7 @@ def run(self): class _ZoneInfoDialog(QDialog): def __init__(self, parent=None): super().__init__(parent) - self.setWindowTitle("Spatial Target — Zone Layout & Interpretation") + self.setWindowTitle("Spatial Detail Target — Zone Layout & Interpretation") self.setMinimumSize(680, 540) lay = QVBoxLayout(self) browser = QTextBrowser() From 5e9f21ad86eccb665830f3bc40b62e1f11f04b36 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 21:46:17 -0400 Subject: [PATCH 10/14] feat: add dB ratio curve to Section 7 power spectrum comparison The log-Y radial power overlay hides small, consistent A-vs-B differences in the noise floor. Add a 10*log10(P_A/P_B) ratio curve in a separate panel (not twinx -- see CLAUDE.md) for both the primary/starless and with-stars pairs, readable directly against a 0 dB reference line. Also documents the new pattern and a pre-existing, unrelated bug found during verification (_section_snr crashes when SNR is unchecked) in CLAUDE.md. Co-Authored-By: Claude Sonnet 5 --- CLAUDE.md | 29 +++++++- report/report_builder.py | 86 +++++++++++++++++++++++- tests/test_report/test_report_helpers.py | 46 +++++++++++++ 3 files changed, 159 insertions(+), 2 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index ec30747..bc74742 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -58,6 +58,7 @@ synthetic/ | `set_starless_path(path)` | `gui/image_panel.py` | Attach a pre-generated starless FITS to the loaded main image | | `_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 | --- @@ -214,6 +215,31 @@ except ImportError: `analysis/halo_analyzer.py` (stacked radial profiles, RDF nanmean/nanstd), `analysis/image_filters.py` (wavelet MAD noise estimate). +### Ratio/comparison curves in report figures — dB convention, avoid twinx() + +When adding a new A-vs-B ratio curve to a report figure (precedent: `_power_ratio_db` / +`_plot_radial_ratio_db` in `report_builder.py`, Section 7's power-spectrum ratio): + +- **dB convention depends on quantity type.** Power quantities (e.g. `radial_power` in + `analysis/power_spectrum.py`, `= abs(fft2d)**2 / N**2`) use `10 * np.log10(ratio)`. + Amplitude-like quantities (e.g. SNR in `analysis/snr_analyzer.py`) use + `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. +- **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 + np.allclose(freq_a, freq_b)` — check shape first, since `np.allclose` raises + `ValueError` on mismatched shapes rather than returning `False`). This is not + 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. + --- ## Collaboration Rules @@ -306,10 +332,11 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | Closure capture in `secondary_xaxis` lambdas | `lambda x: x * ps` inside a loop captures `ps` by reference. Use default-arg capture: `lambda x, p=ps: x * p` to freeze the value at definition time. | | macOS binary blocked by Gatekeeper | CI-built binaries are unsigned. Users must right-click → Open, or run `xattr -dr com.apple.quarantine AstroImageLab` in Terminal. Code signing requires an Apple Developer certificate ($99/year). | | Linux build needs system Qt libraries | PyInstaller must be able to import PyQt6 during analysis. On `ubuntu-latest` run `sudo apt-get install -y libgl1 libegl1 libxcb-cursor0 libxkbcommon-x11-0` before `pip install -r requirements-build.txt`. | -| `PowerSpectrumAnalyzer` crashes on images smaller than 1024 px | `POWER_SPECTRUM_NPIX = 1024`. The auto-select loop is empty when `min(h, w) < 1024`; the fallback produces negative slice indices → non-square region → `_apply_window` shape mismatch. Fix: `N = min(N, h, w)` before the loop, add `+1` to loop upper bounds, clamp fallback with `max(0, ...)`. | +| `PowerSpectrumAnalyzer` crashes on images smaller than 2048 px | `POWER_SPECTRUM_NPIX = 2048`. The auto-select loop is empty when `min(h, w) < 2048`; the fallback produces negative slice indices → non-square region → `_apply_window` shape mismatch. Fix: `N = min(N, h, w)` before the loop, add `+1` to loop upper bounds, clamp fallback with `max(0, ...)`. | | `sigma_clip` mask is scalar `False` when nothing is clipped | `clipped.mask` is `np.ma.nomask` (== `False`) when no values are clipped. `region[False]` silently writes only the first row. Use `np.ma.getmaskarray(clipped)` to get a full bool array, then guard with `.any()`. | | 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` crashes when SNR metric is unchecked | `_plot_snr_pair` (`report_builder.py`) calls `plt.subplots(1, len(panels), ...)`. If `snr_metrics` is `None` (SNR unchecked while another metric, e.g. Power Spectrum, is run) `panels` is empty and matplotlib raises `ValueError: Number of columns must be a positive integer, not 0`. All metric checkboxes default to checked, so this only bites a deliberately partial run. Discovered during Section 7 ratio-curve verification, not yet fixed — needs the same `metrics or {}` guard pattern used everywhere else in this file. | --- diff --git a/report/report_builder.py b/report/report_builder.py index ee25b25..e28134a 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -461,6 +461,27 @@ def _better_worse_class(val_a, val_b, higher_is_better: bool = True) -> tuple[st return ("better", "worse") if val_a <= val_b else ("worse", "better") +def _power_ratio_db(freq_a, rp_a, freq_b, rp_b) -> tuple[np.ndarray, np.ndarray] | None: + """10*log10(P_A/P_B) per frequency bin. Returns None if data is missing or the + two frequency axes are not bin-for-bin aligned.""" + if freq_a is None or rp_a is None or freq_b is None or rp_b is None: + return None + freq_a = np.asarray(freq_a, dtype=float) + freq_b = np.asarray(freq_b, dtype=float) + # Shape check MUST precede allclose: np.allclose raises ValueError on + # incompatible shapes rather than returning False. + if freq_a.shape != freq_b.shape or not np.allclose(freq_a, freq_b): + return None + rp_a = np.asarray(rp_a, dtype=float) + rp_b = np.asarray(rp_b, dtype=float) + positive = np.concatenate([rp_a[rp_a > 0], rp_b[rp_b > 0]]) + # Relative epsilon, mirroring the existing precedent at + # power_spectrum.py:199 (`ps2d[ps2d > 0].min() * 0.01`). + eps = float(positive.min()) * 0.01 if positive.size > 0 else 1e-12 + ratio_db = 10.0 * np.log10(np.clip(rp_a, eps, None) / np.clip(rp_b, eps, None)) + return freq_a, ratio_db + + def _focal_ratio(img: AstroImage) -> float | None: hdr = img.header if hdr is None: @@ -3518,6 +3539,53 @@ def _plot_radial_overlay(self, ra: AnalysisResult, rb: AnalysisResult, fig.tight_layout() return fig + def _plot_radial_ratio_db(self, ra: AnalysisResult, rb: AnalysisResult, + star_pa: dict | None = None, + star_pb: dict | None = None) -> plt.Figure | None: + """dB ratio of A's to B's radial power spectrum. Solid = primary/starless pair; + dashed = with-stars pair (only when both images have star_power data).""" + pa = ra.power_metrics or {} + pb = rb.power_metrics or {} + primary = _power_ratio_db(pa.get("freq_axis"), pa.get("radial_power"), + pb.get("freq_axis"), pb.get("radial_power")) + star = None + if star_pa and star_pb: + star = _power_ratio_db(star_pa.get("freq_axis"), star_pa.get("radial_power"), + star_pb.get("freq_axis"), star_pb.get("radial_power")) + if primary is None and star is None: + return None # neither pair usable — degrade gracefully, no image, no crash + + import matplotlib + _is_dark = matplotlib.rcParams.get("figure.facecolor", "white") not in ("white", "#ffffff", 1.0) + orig_color = "white" if _is_dark else "black" + + fig, ax = plt.subplots(figsize=(7, 4)) + all_vals = [] + if primary is not None: + freq, ratio_db = primary + ax.plot(freq, ratio_db, color="mediumpurple", linewidth=2, label="Ratio (starless)") + all_vals.append(ratio_db) + if star is not None: + freq, ratio_db = star + ax.plot(freq, ratio_db, color="mediumpurple", linewidth=1.5, linestyle="--", + alpha=0.6, label="Ratio (with stars)") + all_vals.append(ratio_db) + + ax.axhline(0.0, color=orig_color, linestyle="--", linewidth=0.8, label="0 dB (A = B)") + ax.axvline(0.10, color="gray", linestyle="--", linewidth=0.8, + label="Low / mid boundary (0.10 cyc/px)") + ax.set_xlabel("Spatial frequency (cycles/pixel)") + ax.set_ylabel("Ratio (dB) = 10·log10(P_A / P_B)") + ax.set_title(f"Radial power ratio (dB): {ra.label} / {rb.label}") + ax.set_xlim(0, 0.5) + peak = float(np.max(np.abs(np.concatenate(all_vals)))) if all_vals else 3.0 + ylim = max(3.0, peak * 1.1) + ax.set_ylim(-ylim, ylim) + ax.legend(fontsize=9) + ax.grid(True, alpha=0.3) + fig.tight_layout() + return fig + # ── Section 7: Power spectrum ────────────────────────────────────────────── def _section_power(self, ra: AnalysisResult, rb: AnalysisResult) -> str: @@ -3536,6 +3604,12 @@ def _section_power(self, ra: AnalysisResult, rb: AnalysisResult) -> str: star_pb=star_pb or None), "Radial power overlay", ) + img_ratio = _img_tag( + self._plot_radial_ratio_db(ra, rb, + star_pa=star_pa or None, + star_pb=star_pb or None), + "Radial power ratio (dB)", + ) sl_note = "" used_a = pa.get("used_starless", False) @@ -3585,7 +3659,9 @@ def _section_power(self, ra: AnalysisResult, rb: AnalysisResult) -> str: 'rather than a DC artifact. The mid/high-frequency ratio (0.1–0.5 cyc/px vs 0–0.1 cyc/px) ' 'measures fine detail content relative to coarse structure.
' 'Note: This comparison is only meaningful when both images cover ' - 'the same target region.', + 'the same target region. The ratio curve below expresses this directly in decibels ' + '(10·log10(PA / PB)): positive values mean Image A has more ' + 'power (finer detail) at that frequency, negative values mean Image B does.', title="About the power spectrum")} @@ -3598,6 +3674,14 @@ def _section_power(self, ra: AnalysisResult, rb: AnalysisResult) -> str: high frequencies indicate one filter preserves more fine-scale spatial detail. The dashed vertical line marks the boundary between low (coarse structure) and mid/high frequencies.

+{img_ratio} +

Ratio of radial power spectra in decibels. Solid = starless-pair +ratio; dashed = with-stars-pair ratio (only when both images used a starless primary +input). Shown only where both images' frequency bins are exactly aligned bin-for-bin — +always true when an analysis ROI is set, true in auto-ROI mode only when both images' +analysed regions share the same pixel dimensions. If only one pair aligns, only that +curve is shown.

+
{img_a}
{img_b}
diff --git a/tests/test_report/test_report_helpers.py b/tests/test_report/test_report_helpers.py index fe554d9..ab3e7bb 100644 --- a/tests/test_report/test_report_helpers.py +++ b/tests/test_report/test_report_helpers.py @@ -16,6 +16,7 @@ _arr_to_b64_png, _fig_to_b64, _psf_stat_test, + _power_ratio_db, ) EM_DASH = "—" @@ -182,3 +183,48 @@ def test_custom_dpi_changes_size(self): lo = _fig_to_b64(self._make_fig(), dpi=72) hi = _fig_to_b64(self._make_fig(), dpi=200) assert len(base64.b64decode(hi)) > len(base64.b64decode(lo)) + + +class TestPowerRatioDb: + def test_known_ratio_in_db(self): + freq = np.linspace(0.0, 0.5, 10) + rp_b = np.full(10, 1.0) + rp_a = rp_b * 2.0 + result = _power_ratio_db(freq, rp_a, freq, rp_b) + assert result is not None + _, ratio_db = result + assert ratio_db == pytest.approx(10.0 * np.log10(2.0), abs=1e-6) + + def test_identical_curves_zero_db(self): + freq = np.linspace(0.0, 0.5, 10) + rp = np.linspace(1.0, 5.0, 10) + _, ratio_db = _power_ratio_db(freq, rp, freq, rp) + assert ratio_db == pytest.approx(0.0, abs=1e-6) + + def test_zero_bins_stay_finite(self): + freq = np.linspace(0.0, 0.5, 5) + rp_a = np.array([0.0, 1.0, 0.0, 2.0, 0.0]) + rp_b = np.array([1.0, 0.0, 0.0, 1.0, 1.0]) + _, ratio_db = _power_ratio_db(freq, rp_a, freq, rp_b) + assert np.all(np.isfinite(ratio_db)) + + def test_none_inputs_return_none(self): + freq = np.linspace(0.0, 0.5, 5) + rp = np.ones(5) + assert _power_ratio_db(None, rp, freq, rp) is None + assert _power_ratio_db(freq, None, freq, rp) is None + assert _power_ratio_db(freq, rp, None, rp) is None + assert _power_ratio_db(freq, rp, freq, None) is None + + def test_mismatched_length_returns_none(self): + freq_a = np.linspace(0.0, 0.5, 5) + freq_b = np.linspace(0.0, 0.5, 8) + rp_a = np.ones(5) + rp_b = np.ones(8) + assert _power_ratio_db(freq_a, rp_a, freq_b, rp_b) is None + + def test_same_length_different_values_returns_none(self): + freq_a = np.linspace(0.0, 0.5, 5) + freq_b = np.linspace(0.0, 0.4, 5) + rp = np.ones(5) + assert _power_ratio_db(freq_a, rp, freq_b, rp) is None From 7f701281556f63b3791074354b1d2fa6e5fb27fe Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Tue, 7 Jul 2026 21:59:21 -0400 Subject: [PATCH 11/14] fix: guard _plot_snr_pair against zero SNR panels _plot_snr_pair built a panels list filtered to non-None entries but never checked whether it ended up empty before calling plt.subplots(1, len(panels)), which raised ValueError: Number of columns must be a positive integer, not 0. Reachable any time SNR is unchecked while another metric (e.g. Power Spectrum) is run. Return None instead, matching the guard pattern already used by _plot_radial_overlay/_plot_radial_ratio_db -- both call sites already pipe the result through _img_tag, which turns None into "". Verified end-to-end via the real ReportBuilder.generate() pipeline: the original crash scenario now succeeds, and a real SNR-populated run still renders both SNR map panels correctly (no regression). Updates the CLAUDE.md pitfall entry documented in the prior commit to reflect the fix. Co-Authored-By: Claude Sonnet 5 --- CLAUDE.md | 4 ++-- report/report_builder.py | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index bc74742..8f74ecf 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -14,7 +14,7 @@ images. It produces a self-contained HTML report with embedded matplotlib figure ## Architecture -``` +```text AstroImageLab.py PyQt6 app + animated splash screen analysis/ Metric engines — each returns a plain dict psf_analyzer.py Moffat/ePSF fitting, MTF via FFT @@ -336,7 +336,7 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html | `sigma_clip` mask is scalar `False` when nothing is clipped | `clipped.mask` is `np.ma.nomask` (== `False`) when no values are clipped. `region[False]` silently writes only the first row. Use `np.ma.getmaskarray(clipped)` to get a full bool array, then guard with `.any()`. | | 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` crashes when SNR metric is unchecked | `_plot_snr_pair` (`report_builder.py`) calls `plt.subplots(1, len(panels), ...)`. If `snr_metrics` is `None` (SNR unchecked while another metric, e.g. Power Spectrum, is run) `panels` is empty and matplotlib raises `ValueError: Number of columns must be a positive integer, not 0`. All metric checkboxes default to checked, so this only bites a deliberately partial run. Discovered during Section 7 ratio-curve verification, not yet fixed — needs the same `metrics or {}` guard pattern used everywhere else in this file. | +| `_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 `""`. | --- diff --git a/report/report_builder.py b/report/report_builder.py index e28134a..d818f2f 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -3925,10 +3925,12 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: # ── Section 9: Signal-to-Noise Ratio ───────────────────────────────────── def _plot_snr_pair(self, disp_a, label_a, disp_b, label_b, - vmin: float, vmax: float) -> plt.Figure: + vmin: float, vmax: float) -> plt.Figure | None: """Render both SNR maps side-by-side with a shared plasma color scale.""" panels = [(d, lbl) for d, lbl in [(disp_a, label_a), (disp_b, label_b)] if d is not None] + if not panels: + return None fig, axes = plt.subplots(1, len(panels), figsize=(7 * len(panels), 5)) if len(panels) == 1: axes = [axes] From 3502d0585a15facf0ab3fc0c04df53c23daffc8b Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Wed, 8 Jul 2026 21:04:16 -0400 Subject: [PATCH 12/14] feat: add noise-corrected multi-scale local contrast to Section 8 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a shared-nebula-ROI, per-scale A/B ratio metric across std, LoG, wavelet, and Weber contrast (score = median(|detail|) over the region both images agree is nebula, divided by each image's own background noise floor at that scale), plus a new Gradient Magnitude / Edge Sharpness method and noise-normalised display maps so filter A vs B comparisons aren't confounded by differing absolute noise levels. Also fixes the Report Inspector's panel catalog, which used a hardcoded map that had gone stale against STD_KERNEL_SIZES and never included Weber maps at all — it's now derived dynamically from the computed panels so nothing can silently disappear from the Inspector again. Co-Authored-By: Claude Sonnet 5 --- analysis/image_filters.py | 382 +++++++++++++++++++- report/report_builder.py | 150 +++++++- tests/test_analysis/test_spatial_detail.py | 227 ++++++++++++ tests/test_report/test_inspector_catalog.py | 93 +++++ tests/test_report/test_report_helpers.py | 58 +++ 5 files changed, 883 insertions(+), 27 deletions(-) create mode 100644 tests/test_report/test_inspector_catalog.py diff --git a/analysis/image_filters.py b/analysis/image_filters.py index 0c438bc..c44dc81 100644 --- a/analysis/image_filters.py +++ b/analysis/image_filters.py @@ -10,7 +10,7 @@ import matplotlib.colors as mcolors matplotlib.use("Agg") import matplotlib.pyplot as plt -from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter +from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter import pywt from core.astro_image import AstroImage @@ -58,6 +58,17 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None, "weber_contrast_a": {}, "weber_contrast_b": {}, "panels": {}, + "nc_shared_nebula_pixels": 0, + "std_nc_score_a": {}, "std_nc_score_b": {}, + "std_nc_noise_a": {}, "std_nc_noise_b": {}, "std_nc_ratio": {}, + "log_nc_score_a": {}, "log_nc_score_b": {}, + "log_nc_noise_a": {}, "log_nc_noise_b": {}, "log_nc_ratio": {}, + "wavelet_nc_score_a": {}, "wavelet_nc_score_b": {}, + "wavelet_nc_noise_a": {}, "wavelet_nc_noise_b": {}, "wavelet_nc_ratio": {}, + "weber_nc_score_a": {}, "weber_nc_score_b": {}, + "weber_nc_noise_a": {}, "weber_nc_noise_b": {}, "weber_nc_ratio": {}, + "gm_nc_score_a": {}, "gm_nc_score_b": {}, + "gm_nc_noise_a": {}, "gm_nc_noise_b": {}, "gm_nc_ratio": {}, } figures: dict = {} @@ -109,11 +120,23 @@ def _clip01(v): return max(0.0, min(1.0, v)) result["display_roi"] = display_roi + # 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 + ) + _label_b = image_b.label if image_b is not None else None - # 1-3. Local std, LoG, wavelet — all read norm_a/norm_b with no shared mutable state, - # so they run concurrently. Each method returns (b64_figs, partial_result). - with concurrent.futures.ThreadPoolExecutor(max_workers=4) as _ex: + # 1-5. Local std, LoG, wavelet, Weber, gradient — all read norm_a/norm_b with no + # shared mutable state, so they run concurrently. Each method returns + # (b64_figs, partial_result). + with concurrent.futures.ThreadPoolExecutor(max_workers=5) as _ex: _f_std = _ex.submit(self._std_analysis, analysis_a, analysis_b, mask_neb_a, mask_bg_a, @@ -122,34 +145,47 @@ def _clip01(v): return max(0.0, min(1.0, v)) image_a.label, _label_b, display_roi=display_roi, crosshair=crosshair_roi, + mask_neb_shared=mask_neb_shared, ) _f_log = _ex.submit(self._log_analysis, analysis_a, analysis_b, log_sigmas, 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, ) _f_wav = _ex.submit(self._wavelet_analysis, analysis_a, analysis_b, wavelet, levels, 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, ) _f_web = _ex.submit(self._weber_analysis, analysis_a, analysis_b, weber_kernel_sizes, image_a.label, _label_b, display_roi=display_roi, + mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b, + ) + _f_grad = _ex.submit(self._gradient_analysis, + analysis_a, analysis_b, log_sigmas, + 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, ) std_b64, std_partial = _f_std.result() log_b64, log_partial = _f_log.result() wav_b64, wav_partial = _f_wav.result() web_b64, web_partial = _f_web.result() + grad_b64, grad_partial = _f_grad.result() figures.update(std_b64) figures.update(log_b64) figures.update(wav_b64) figures.update(web_b64) + figures.update(grad_b64) result["contrast_ratios_a"].update(std_partial["contrast_ratios_a"]) result["contrast_ratios_b"].update(std_partial["contrast_ratios_b"]) result["sigma_noise_a"] = wav_partial["sigma_noise_a"] @@ -162,6 +198,27 @@ def _clip01(v): return max(0.0, min(1.0, v)) result["panels"].update(log_partial["panels"]) result["panels"].update(wav_partial["panels"]) result["panels"].update(web_partial["panels"]) + result["panels"].update(grad_partial["panels"]) + + # Merge noise-corrected scores/noise-floors and compute A/B ratios centrally. + for prefix, partial in (("std", std_partial), ("log", log_partial), + ("wavelet", wav_partial), ("weber", web_partial), + ("gm", grad_partial)): + for suffix in ("nc_score_a", "nc_score_b", "nc_noise_a", "nc_noise_b"): + result[f"{prefix}_{suffix}"].update(partial[f"{prefix}_{suffix}"]) + result[f"{prefix}_nc_ratio"] = self._compute_nc_ratios( + result[f"{prefix}_nc_score_a"], result[f"{prefix}_nc_score_b"]) + + if image_b is not None: + nc_fig = self._plot_nc_ratio_overview({ + "std": result["std_nc_ratio"], + "log": result["log_nc_ratio"], + "wavelet": result["wavelet_nc_ratio"], + "weber": result["weber_nc_ratio"], + "gradient": result["gm_nc_ratio"], + }) + if nc_fig is not None: + figures["nc_ratio_overview"] = fig_to_b64(nc_fig, dpi=150) if crosshair is not None: pos_a, prof_a = self._sample_line(norm_a, **crosshair) @@ -269,9 +326,15 @@ def _std_analysis(self, norm_a, norm_b, mask_neb_b, mask_bg_b, kernel_sizes, label_a, label_b, display_roi=None, - crosshair=None) -> tuple[dict, dict]: + crosshair=None, + mask_neb_shared=None) -> tuple[dict, dict]: figures = {} - partial: dict = {"contrast_ratios_a": {}, "contrast_ratios_b": {}, "panels": {}} + 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": {}, + } single = norm_b is None for ks in kernel_sizes: std_a = self._compute_std_map(norm_a, ks) @@ -284,11 +347,26 @@ def _std_analysis(self, norm_a, norm_b, cr_b = self._contrast_ratio(std_b, mask_neb_b, mask_bg_b) partial["contrast_ratios_b"][ks] = cr_b + noise_a = noise_b = None + if not single: + nc_a, noise_a = self._nc_score(std_a, mask_neb_shared, mask_bg_a) + partial["std_nc_score_a"][ks] = nc_a + partial["std_nc_noise_a"][ks] = noise_a + nc_b, noise_b = self._nc_score(std_b, mask_neb_shared, mask_bg_b) + partial["std_nc_score_b"][ks] = nc_b + partial["std_nc_noise_b"][ks] = noise_b + 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, } + if not single and noise_a and noise_b: + partial["panels"][f"nrm_std_{ks}px"] = { + "a": (std_a / noise_a).astype(np.float32), + "b": (std_b / noise_b).astype(np.float32), + "diff": None, + } if not single: fig = self._plot_side_by_side( @@ -310,6 +388,17 @@ def _std_analysis(self, norm_a, norm_b, ) figures[f"std_{ks}px"] = fig + if not single and noise_a and noise_b: + 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", + cmap=SECTION8_ANALYSIS_CMAP, + display_roi=None, + ) + if crosshair is not None: pos, pa = self._sample_line(std_a, **crosshair) if not single: @@ -352,6 +441,43 @@ def _contrast_ratio(self, std_map: np.ndarray, return None return float(np.median(neb_vals)) / bg_med + def _nc_score(self, detail_map: np.ndarray, + mask_neb_shared: np.ndarray | None, + mask_bg: np.ndarray) -> tuple[float | None, float | None]: + """Noise-corrected local-contrast score for one detail map at one scale. + + score = median(|detail|) over the pixels BOTH images classify as nebula + (mask_neb_shared), divided by median(|detail|) over THIS image's own + background mask — its empirical per-scale noise floor for this operator. + Returns (score, noise_floor); either is None if a mask selects zero pixels, + mask_neb_shared is unavailable (single-image mode), or noise_floor <= 0. + """ + if mask_neb_shared is None: + return None, None + h = min(detail_map.shape[0], mask_neb_shared.shape[0], mask_bg.shape[0]) + w = min(detail_map.shape[1], mask_neb_shared.shape[1], mask_bg.shape[1]) + absmap = np.abs(detail_map[:h, :w]) + neb_vals = absmap[mask_neb_shared[:h, :w]] + bg_vals = absmap[mask_bg[:h, :w]] + if neb_vals.size == 0 or bg_vals.size == 0: + return None, None + noise_floor = float(bn.median(bg_vals)) + if noise_floor <= 0: + return None, None + return float(bn.median(neb_vals)) / noise_floor, noise_floor + + @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 + (single-image mode).""" + if not score_a or not score_b: + return {} + out = {} + for scale, va in score_a.items(): + vb = score_b.get(scale) + out[scale] = None if (va is None or vb is None or vb == 0) else va / vb + return out + # ------------------------------------------------------------------ # Laplacian of Gaussian maps # ------------------------------------------------------------------ @@ -359,18 +485,40 @@ def _contrast_ratio(self, std_map: np.ndarray, def _log_analysis(self, norm_a, norm_b, sigmas, label_a, label_b, display_roi=None, - crosshair=None) -> tuple[dict, dict]: + crosshair=None, + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: figures = {} - partial: dict = {"panels": {}} + partial: dict = { + "log_nc_score_a": {}, "log_nc_score_b": {}, + "log_nc_noise_a": {}, "log_nc_noise_b": {}, + "panels": {}, + } single = norm_b is None for sigma in sigmas: log_a = np.abs(gaussian_laplace(norm_a, sigma=sigma)) log_b = np.abs(gaussian_laplace(norm_b, sigma=sigma)) if not single else None + + noise_a = noise_b = None + if not single: + nc_a, noise_a = self._nc_score(log_a, mask_neb_shared, mask_bg_a) + partial["log_nc_score_a"][sigma] = nc_a + partial["log_nc_noise_a"][sigma] = noise_a + nc_b, noise_b = self._nc_score(log_b, mask_neb_shared, mask_bg_b) + partial["log_nc_score_b"][sigma] = nc_b + partial["log_nc_noise_b"][sigma] = noise_b + 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, } + if not single and noise_a and noise_b: + partial["panels"][f"nrm_log_{sigma}"] = { + "a": (log_a / noise_a).astype(np.float32), + "b": (log_b / noise_b).astype(np.float32), + "diff": None, + } + if not single: fig = self._plot_side_by_side( self._crop_border(log_a, SECTION8_BORDER_CROP_FRACTION), @@ -390,6 +538,18 @@ def _log_analysis(self, norm_a, norm_b, sigmas, nonlinear_norm=True, ) figures[f"log_sigma{sigma}"] = fig + + if not single and noise_a and noise_b: + 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", + cmap=SECTION8_ANALYSIS_CMAP, + display_roi=None, + ) + if crosshair is not None and not single: pos, pa = self._sample_line(log_a, **crosshair) _, pb = self._sample_line(log_b, **crosshair) @@ -398,6 +558,89 @@ def _log_analysis(self, norm_a, norm_b, sigmas, f"Cross-section — |LoG|, σ={sigma}px") return figs_to_b64(figures, dpi=150), partial + # ------------------------------------------------------------------ + # Gradient magnitude (edge sharpness) + # ------------------------------------------------------------------ + + 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]: + """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.""" + figures = {} + partial: dict = { + "gm_nc_score_a": {}, "gm_nc_score_b": {}, + "gm_nc_noise_a": {}, "gm_nc_noise_b": {}, + "panels": {}, + } + single = norm_b is None + for sigma in sigmas: + gm_a = gaussian_gradient_magnitude(norm_a, sigma=sigma) + gm_b = gaussian_gradient_magnitude(norm_b, sigma=sigma) if not single else None + + noise_a = noise_b = None + if not single: + nc_a, noise_a = self._nc_score(gm_a, mask_neb_shared, mask_bg_a) + partial["gm_nc_score_a"][sigma] = nc_a + partial["gm_nc_noise_a"][sigma] = noise_a + nc_b, noise_b = self._nc_score(gm_b, mask_neb_shared, mask_bg_b) + partial["gm_nc_score_b"][sigma] = nc_b + partial["gm_nc_noise_b"][sigma] = noise_b + + 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, + } + if not single and noise_a and noise_b: + partial["panels"][f"nrm_gradient_{sigma}"] = { + "a": (gm_a / noise_a).astype(np.float32), + "b": (gm_b / noise_b).astype(np.float32), + "diff": None, + } + + 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", + cmap=SECTION8_ANALYSIS_CMAP, + nonlinear_norm=True, + display_roi=None, + ) + else: + fig = self._plot_single( + self._crop_border(gm_a, SECTION8_BORDER_CROP_FRACTION), + f"Gradient |G| σ={sigma}px — {label_a}", + cmap=SECTION8_ANALYSIS_CMAP, + nonlinear_norm=True, + ) + figures[f"gradient_{sigma}"] = fig + + if not single and noise_a and noise_b: + 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", + cmap=SECTION8_ANALYSIS_CMAP, + display_roi=None, + ) + + 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 + # ------------------------------------------------------------------ # Wavelet decomposition # ------------------------------------------------------------------ @@ -405,11 +648,14 @@ def _log_analysis(self, norm_a, norm_b, sigmas, def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, label_a, label_b, display_roi=None, - crosshair=None) -> tuple[dict, dict]: + crosshair=None, + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: figures = {} partial: dict = { "sigma_noise_a": None, "sigma_noise_b": None, "wavelet_snr_a": {}, "wavelet_snr_b": {}, + "wavelet_nc_score_a": {}, "wavelet_nc_score_b": {}, + "wavelet_nc_noise_a": {}, "wavelet_nc_noise_b": {}, "panels": {}, } single = norm_b is None @@ -434,18 +680,39 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, figures["wavelet_snr"] = self._plot_snr_bars( partial["wavelet_snr_a"], partial["wavelet_snr_b"], label_a, label_b, levels) - # Reconstruct and display levels 2 and 3 (best signal content) - for display_level in [2, 3]: - if display_level > levels: - continue - coeff_idx = levels + 1 - display_level + # Reconstruct every level for noise-corrected scoring (raw subband coefficients + # live at reduced spatial resolution and don't align pixel-for-pixel with the + # full-resolution nebula mask — only the inverse-DWT output does). Display + # figures/panels stay restricted to levels 2-3 ("best signal content"). + for human_level in range(1, levels + 1): + coeff_idx = levels + 1 - human_level rec_a = self._reconstruct_level(coeffs_a, coeff_idx, wavelet, levels) rec_b = self._reconstruct_level(coeffs_b, coeff_idx, wavelet, levels) if coeffs_b is not None else None + + noise_a = noise_b = None + if not single: + nc_a, noise_a = self._nc_score(rec_a, mask_neb_shared, mask_bg_a) + partial["wavelet_nc_score_a"][human_level] = nc_a + partial["wavelet_nc_noise_a"][human_level] = noise_a + nc_b, noise_b = self._nc_score(rec_b, mask_neb_shared, mask_bg_b) + partial["wavelet_nc_score_b"][human_level] = nc_b + partial["wavelet_nc_noise_b"][human_level] = noise_b + + if human_level not in (2, 3): + continue # display/panels only for levels 2-3, unchanged from prior behaviour + + display_level = human_level 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, } + 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, + } if not single: fig = self._plot_side_by_side( self._crop_border(rec_a, SECTION8_BORDER_CROP_FRACTION), @@ -464,6 +731,19 @@ def _wavelet_analysis(self, norm_a, norm_b, wavelet, levels, cmap=SECTION8_ANALYSIS_CMAP, ) figures[f"wavelet_level{display_level}"] = fig + + if not single and noise_a and noise_b: + 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", + cmap=SECTION8_ANALYSIS_CMAP, + symmetric_diff=True, + display_roi=None, + ) + if crosshair is not None and not single: pos, pa = self._sample_line(rec_a, **crosshair) _, pb = self._sample_line(rec_b, **crosshair) @@ -511,11 +791,14 @@ 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) -> tuple[dict, dict]: + display_roi=None, + mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None) -> tuple[dict, dict]: figures = {} partial: dict = { "weber_contrast_a": {}, "weber_contrast_b": {}, + "weber_nc_score_a": {}, "weber_nc_score_b": {}, + "weber_nc_noise_a": {}, "weber_nc_noise_b": {}, "panels": {}, } single = norm_b is None @@ -528,11 +811,26 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes, if not single: partial["weber_contrast_b"][ks] = float(np.percentile(wc_b, 99)) + noise_a = noise_b = None + if not single: + nc_a, noise_a = self._nc_score(wc_a, mask_neb_shared, mask_bg_a) + partial["weber_nc_score_a"][ks] = nc_a + partial["weber_nc_noise_a"][ks] = noise_a + nc_b, noise_b = self._nc_score(wc_b, mask_neb_shared, mask_bg_b) + partial["weber_nc_score_b"][ks] = nc_b + partial["weber_nc_noise_b"][ks] = noise_b + 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, } + if not single and noise_a and noise_b: + partial["panels"][f"nrm_weber_{ks}px"] = { + "a": (wc_a / noise_a).astype(np.float32), + "b": (wc_b / noise_b).astype(np.float32), + "diff": None, + } if not single: fig = self._plot_side_by_side( @@ -554,6 +852,17 @@ 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: + 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", + cmap=SECTION8_ANALYSIS_CMAP, + display_roi=None, + ) + return figs_to_b64(figures, dpi=150), partial def _compute_weber_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray: @@ -708,6 +1017,49 @@ def _plot_snr_bars(self, snr_a: dict, snr_b: dict, fig.tight_layout() return fig + def _plot_nc_ratio_overview(self, ratios_by_method: dict) -> plt.Figure | None: + """One line per method: noise-corrected A/B ratio vs. approximate spatial + scale (px, log-x). None if no method has any usable (non-None) value.""" + _SCALE_LABEL = { + "std": "px", "weber": "px", "log": "σ px", + "gradient": "σ px", "wavelet": "level (≈px)", + } + _COLORS = { + "std": "steelblue", "log": "tomato", "wavelet": "mediumpurple", + "weber": "seagreen", "gradient": "goldenrod", + } + series = {} + for method, ratios in ratios_by_method.items(): + if not ratios: + continue + if method == "wavelet": + pts = [(2 ** scale, v) for scale, v in ratios.items() if v is not None] + else: + pts = [(float(scale), v) for scale, v in ratios.items() if v is not None] + if pts: + pts.sort(key=lambda p: p[0]) + series[method] = pts + + if not series: + return None + + fig, ax = plt.subplots(figsize=(7, 4.5)) + for method, pts in series.items(): + xs = [p[0] for p in pts] + ys = [p[1] for p in pts] + ax.plot(xs, ys, marker="o", label=f"{method} ({_SCALE_LABEL.get(method, 'px')})", + color=_COLORS.get(method)) + ax.axhline(1.0, color="black", linestyle="--", linewidth=0.8, + label="Ratio = 1 (A = B)") + ax.set_xscale("log") + ax.set_xlabel("Approximate spatial scale (px)") + ax.set_ylabel("Noise-corrected score ratio (A / B)") + ax.set_title("Noise-corrected local contrast — cross-method overview") + ax.legend(fontsize=8) + ax.grid(True, alpha=0.3) + fig.tight_layout() + return fig + @staticmethod def _sample_line(arr: np.ndarray, x0: float, y0: float, x1: float, y1: float) -> tuple[np.ndarray, np.ndarray]: diff --git a/report/report_builder.py b/report/report_builder.py index d818f2f..76fcfbb 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -461,6 +461,43 @@ def _better_worse_class(val_a, val_b, higher_is_better: bool = True) -> tuple[st return ("better", "worse") if val_a <= val_b else ("worse", "better") +def _nc_ratio_rows(score_a: dict, score_b: dict, ratio: dict, scale_label, val_fmt: str = ".3f") -> str: + """Build
rows for a noise-corrected score table: scale | A | B | Ratio A/B. + scale_label(scale) -> row label string. Ratio cell is colored relative to + 1.0 (parity), independently of the A/B columns' own coloring.""" + rows = "" + for scale in sorted(set(list(score_a.keys()) + list(score_b.keys()))): + va, vb, vr = score_a.get(scale), score_b.get(scale), ratio.get(scale) + ca, cb = _better_worse_class(va, vb) + cr, _ = _better_worse_class(vr, 1.0) + rows += (f"" + f"" + f"" + f"") + return rows + + +def _panel_display_name(pkey: str) -> str: + """Derive a human-readable label from a SpatialDetailAnalyzer panels dict key. + Dynamic (not a hardcoded map) so every scale it computes — including future + additions — automatically appears in the Report Inspector catalog.""" + nrm = pkey.startswith("nrm_") + base = pkey[4:] if nrm else pkey + if base.startswith("std_") and base.endswith("px"): + name = f"Std Dev {base[4:-2]} px" + elif base.startswith("log_"): + name = f"LoG σ {base[4:]} px" + elif base.startswith("wavelet_"): + name = f"Wavelet level {base[8:]}" + elif base.startswith("weber_") and base.endswith("px"): + name = f"Weber {base[6:-2]} px" + elif base.startswith("gradient_"): + name = f"Gradient σ {base[9:]} px" + else: + name = base + return f"{name} (noise-normalized)" if nrm else name + + def _power_ratio_db(freq_a, rp_a, freq_b, rp_b) -> tuple[np.ndarray, np.ndarray] | None: """10*log10(P_A/P_B) per frequency bin. Returns None if data is missing or the two frequency axes are not bin-for-bin aligned.""" @@ -787,19 +824,10 @@ def _add_options_entry(section: str, name: str, options: dict) -> None: {"Image A": "gm_display_a", "Image B": "gm_display_b"}) # ── Spatial Detail subsections ──────────────────────────────────────── - _PANEL_IMAGE_SETS = { - "std_5px": "Std Dev 5 px", - "std_15px": "Std Dev 15 px", - "std_31px": "Std Dev 31 px", - "log_1.5": "LoG σ 1.5 px", - "log_3.0": "LoG σ 3.0 px", - "log_6.0": "LoG σ 6.0 px", - "wavelet_2": "Wavelet level 2", - "wavelet_3": "Wavelet level 3", - } panels_a = (result_a.spatial_metrics or {}).get("panels", {}) panels_b = (result_b.spatial_metrics or {}).get("panels", {}) - for pkey, img_set_name in _PANEL_IMAGE_SETS.items(): + for pkey in sorted(panels_a.keys()): + img_set_name = _panel_display_name(pkey) pa_panel = panels_a.get(pkey) if pa_panel is None: continue @@ -3736,6 +3764,46 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: sigma_a = _val(sm.get("sigma_noise_a"), ".5f") sigma_b = _val(sm.get("sigma_noise_b"), ".5f") + # Noise-corrected local contrast ratio tables (one per method) + std_nc_rows = _nc_ratio_rows( + sm.get("std_nc_score_a", {}), sm.get("std_nc_score_b", {}), + sm.get("std_nc_ratio", {}), lambda ks: f"{ks} px") + log_nc_rows = _nc_ratio_rows( + sm.get("log_nc_score_a", {}), sm.get("log_nc_score_b", {}), + sm.get("log_nc_ratio", {}), lambda s: f"σ = {s} px") + wavelet_nc_rows = _nc_ratio_rows( + sm.get("wavelet_nc_score_a", {}), sm.get("wavelet_nc_score_b", {}), + sm.get("wavelet_nc_ratio", {}), lambda lvl: f"Level {lvl} (~{2 ** lvl}px scale)") + weber_nc_rows = _nc_ratio_rows( + sm.get("weber_nc_score_a", {}), sm.get("weber_nc_score_b", {}), + sm.get("weber_nc_ratio", {}), lambda ks: f"{ks} px") + gm_nc_rows = _nc_ratio_rows( + sm.get("gm_nc_score_a", {}), sm.get("gm_nc_score_b", {}), + sm.get("gm_nc_ratio", {}), lambda s: f"σ = {s} px") + + nc_methodology_box = _info_box( + 'Each detail map above additionally yields a noise-corrected local ' + 'contrast score at every scale: score = median(|detail|) over the ' + 'region BOTH images classify as nebula (intersection of each image\'s own ' + 'nebula mask), divided by median(|detail|) over that image\'s OWN background ' + 'region (its empirical per-scale noise floor). Ratio A/B ' + 'divides Image A\'s score by Image B\'s score at each scale — greater than 1 ' + 'means Image A shows relatively stronger detail than Image B at that scale, ' + '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 ' + '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.', + title="Noise-corrected local contrast (methodology)") + nc_empty_note = "" + if not getattr(self, "_single_image", False) and sm.get("nc_shared_nebula_pixels", 0) == 0: + nc_empty_note = _info_box( + '⚠ Images A and B share no common nebula-mask pixels — noise-corrected ' + 'scores below are unavailable (—) for every scale and method.', + title="No shared nebula region", open=True) + def figs_for(prefix): out = "" for key in sorted(figs): @@ -3834,6 +3902,8 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: 'across different filter bandwidths. Images are shown side-by-side with a shared ' 'colour scale; the third panel shows the difference A−B.', title="Spatial detail maps overview")} +{nc_methodology_box} +{nc_empty_note}

8b. Local Standard Deviation Maps

{_info_box('Measures how much pixel values vary within a neighbourhood. ' @@ -3851,10 +3921,16 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: {cr_rows}
{scale_label(scale)}{_val(va, val_fmt)}{_val(vb, val_fmt)}{_val(vr, val_fmt)}
Kernel size{ra.label}{rb.label}
+ + + {std_nc_rows} +
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.

+

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

+{figs_for("nrm_std_")}

8c. Laplacian of Gaussian (LoG) Maps

{_info_box('The Laplacian of Gaussian highlights regions of rapid intensity ' 'change at a specific spatial scale (controlled by σ). Brighter regions in |LoG| maps ' @@ -3867,10 +3943,16 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: 'brighter LoG response at small σ values. Cross-section profiles reveal 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 σ.

+

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

+{figs_for("nrm_log_")}

8d. Wavelet Decomposition

{_wavelet_box} @@ -3882,11 +3964,17 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: Wavelet level{ra.label} SNR{rb.label} SNR {snr_rows} + + + {wavelet_nc_rows} +
Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B
{paired_figs_for("wavelet_level", "xs_wavelet_level")}

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.

+

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

+{figs_for("nrm_wavelet_")}

8e. Weber Fraction Contrast Maps

{_info_box( @@ -3915,12 +4003,50 @@ def paired_figs_for(img_prefix: str, xs_prefix: str) -> str: Kernel size{ra.label} (99th pct c){rb.label} (99th pct c) {wc_rows} + + + {weber_nc_rows} +
Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B
{figs_for("weber_")}

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.

""" +expected; use a nebula ROI for meaningful filter comparison.

+

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

+{figs_for("nrm_weber_")} + +

8f. Gradient Magnitude (Edge Sharpness)

+{_info_box('The gradient magnitude G = |∇I| = sqrt((∂I/∂x)² + (∂I/∂y)²) highlights regions ' + 'of rapid intensity change at a specific spatial scale (controlled by σ), computed ' + 'as the first spatial derivative rather than the second derivative (curvature) LoG uses. ' + 'Smaller σ highlights finer features; larger σ highlights broader structures. ' + 'Reuses the same σ scales as 8c so gradient and |LoG| are directly comparable at ' + 'identical spatial frequencies. This measures how abrupt structure boundaries are — ' + 'relevant when one filter renders stronger filament boundaries or shock fronts. ' + 'Distinct from Section 6 (Edge Detection): Section 6 measures ' + 'precise sub-pixel edge width and contrast at the 2–3 strongest individual detected ' + 'edges. This section is the opposite granularity — an aggregate, whole-nebula, ' + 'multi-scale sharpness score across the shared ROI, following the same ' + 'noise-corrected framework as 8b–8e. They are complementary, not redundant.', + title="Gradient magnitude / edge sharpness")} + + + {gm_nc_rows} +
Scale{ra.label} (NC score){rb.label} (NC score)Ratio A/B
+{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.

+

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

+{figs_for("nrm_gradient_")} + +

8h. 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 +boxes) — use this chart to spot which spatial-scale regime favours which filter, +not to compare absolute ratio values across methods.

""" # ── Section 9: Signal-to-Noise Ratio ───────────────────────────────────── diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py index 062844f..0f36851 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -7,6 +7,7 @@ from analysis.image_filters import SpatialDetailAnalyzer from core.astro_image import AstroImage +from core.models import STD_KERNEL_SIZES, LOG_SIGMAS, WEBER_KERNEL_SIZES, WAVELET_LEVELS class TestAnalyze: @@ -74,3 +75,229 @@ def test_with_roi(self, astro_image_a): result = SpatialDetailAnalyzer().analyze(astro_image_a, roi=(50, 50, 450, 450)) assert isinstance(result, dict) + + +# --------------------------------------------------------------------------- +# Noise-corrected multi-scale local contrast (two-image A/B comparison) +# --------------------------------------------------------------------------- + +def _make_nc_test_fits(path, add_texture: bool, seed: int) -> None: + """256x256 FITS with a smooth nebula blob (sigma=25, well above 2*rms after + background subtraction). When add_texture, a fine checkerboard (period 6px, + amplitude 8x sky noise) is added inside the blob so std/LoG/wavelet/Weber/ + gradient all detect meaningfully more local structure than the plain blob.""" + rng = np.random.default_rng(seed) + h, w = 256, 256 + sky_level = 1000.0 + sky_noise = 20.0 + data = rng.normal(sky_level, sky_noise, (h, w)).astype(np.float64) + + yy, xx = np.mgrid[0:h, 0:w] + cx, cy = w / 2, h / 2 + blob_sigma = 25.0 + blob_amp = 600.0 + blob = blob_amp * np.exp(-0.5 * (((xx - cx) ** 2 + (yy - cy) ** 2) / blob_sigma ** 2)) + data += blob + + if add_texture: + texture_amp = 8.0 * sky_noise + checker = (((xx // 3).astype(int) + (yy // 3).astype(int)) % 2 == 0) + blob_mask = blob > (0.3 * blob_amp) + data += np.where(blob_mask & checker, texture_amp, 0.0) + + data = np.clip(data, 0, 65535).astype(np.float32) + hdr = ap_fits.Header() + hdr["EGAIN"] = 1.0 + hdr["GAIN"] = 1.0 + hdr["FOCALLEN"] = 500.0 + hdr["XPIXSZ"] = 3.76 + hdr["EXPTIME"] = 300.0 + hdr["INSTRUME"] = "TestCam" + ap_fits.writeto(str(path), data, hdr, overwrite=True) + + +@pytest.fixture(scope="module") +def nc_image_pair(tmp_path_factory) -> tuple[AstroImage, AstroImage]: + """Image A = blob + fine checkerboard texture; Image B = blob only. + Both share the same blob location/extent, so their nebula masks overlap.""" + out = tmp_path_factory.mktemp("nc_fits") + path_a = out / "nc_a.fits" + path_b = out / "nc_b.fits" + _make_nc_test_fits(path_a, add_texture=True, seed=1) + _make_nc_test_fits(path_b, add_texture=False, seed=2) + img_a = AstroImage(str(path_a), label="A") + img_a.load() + img_a.estimate_background() + img_b = AstroImage(str(path_b), label="B") + img_b.load() + img_b.estimate_background() + return img_a, img_b + + +@pytest.fixture(scope="module") +def nc_result(nc_image_pair) -> dict: + img_a, img_b = nc_image_pair + return SpatialDetailAnalyzer().analyze(img_a, img_b) + + +class TestNoiseCorrectedContrast: + """A/B noise-corrected multi-scale local contrast: shared-nebula-ROI scoring, + per-scale A/B ratios, and noise-normalised display panels.""" + + def test_shared_nebula_pixels_positive(self, nc_result): + assert nc_result["nc_shared_nebula_pixels"] > 0 + + @pytest.mark.parametrize("prefix,scales", [ + ("std", STD_KERNEL_SIZES), + ("log", LOG_SIGMAS), + ("weber", WEBER_KERNEL_SIZES), + ("gm", LOG_SIGMAS), + ]) + def test_nc_score_dict_keys_match_scales(self, nc_result, prefix, scales): + assert set(nc_result[f"{prefix}_nc_score_a"].keys()) == set(scales) + assert set(nc_result[f"{prefix}_nc_score_b"].keys()) == set(scales) + + def test_wavelet_nc_score_keys_match_levels(self, nc_result): + expected = set(range(1, WAVELET_LEVELS + 1)) + assert set(nc_result["wavelet_nc_score_a"].keys()) == expected + assert set(nc_result["wavelet_nc_score_b"].keys()) == expected + + @pytest.mark.parametrize("prefix", ["std", "log", "wavelet", "weber", "gm"]) + def test_nc_noise_floor_positive_or_none(self, nc_result, prefix): + for side in ("a", "b"): + for v in nc_result[f"{prefix}_nc_noise_{side}"].values(): + if v is not None: + assert v > 0.0 + + def test_std_nc_ratio_captures_finer_detail_in_a(self, nc_result): + ratio = nc_result["std_nc_ratio"][min(STD_KERNEL_SIZES)] + assert ratio is not None and ratio > 1.05 + + def test_log_nc_ratio_captures_finer_detail_in_a(self, nc_result): + ratio = nc_result["log_nc_ratio"][min(LOG_SIGMAS)] + assert ratio is not None and ratio > 1.05 + + def test_wavelet_nc_ratio_captures_finer_detail_in_a(self, nc_result): + ratio = nc_result["wavelet_nc_ratio"][2] + assert ratio is not None and ratio > 1.05 + + def test_weber_nc_ratio_captures_finer_detail_in_a(self, nc_result): + ratio = nc_result["weber_nc_ratio"][min(WEBER_KERNEL_SIZES)] + assert ratio is not None and ratio > 1.05 + + def test_gradient_nc_ratio_captures_finer_detail_in_a(self, nc_result): + # Gradient magnitude's peak response scale for a given texture need not be + # the finest sigma (Gaussian smoothing at small sigma can attenuate a very + # fine checkerboard before the derivative is taken) — check the strongest + # response across scales rather than pinning to one bin. + values = [v for v in nc_result["gm_nc_ratio"].values() if v is not None] + assert values and max(values) > 1.05 + + def test_normalized_panels_present_two_image(self, nc_result): + panels = nc_result["panels"] + for ks in STD_KERNEL_SIZES: + assert f"nrm_std_{ks}px" in panels + for sigma in LOG_SIGMAS: + assert f"nrm_log_{sigma}" in panels + assert f"nrm_gradient_{sigma}" in panels + for ks in WEBER_KERNEL_SIZES: + assert f"nrm_weber_{ks}px" in panels + for lvl in (2, 3): + assert f"nrm_wavelet_{lvl}" in panels + + def test_normalized_panel_values_differ_from_raw(self, nc_result): + panels = nc_result["panels"] + ks = min(STD_KERNEL_SIZES) + raw_a = panels[f"std_{ks}px"]["a"] + nrm_a = panels[f"nrm_std_{ks}px"]["a"] + assert not np.allclose(raw_a, nrm_a) + + def test_nc_ratio_overview_figure_present(self, nc_result): + assert "nc_ratio_overview" in nc_result["figures"] + + def test_with_roi_two_image_no_crash(self, nc_image_pair): + img_a, img_b = nc_image_pair + result = SpatialDetailAnalyzer().analyze( + img_a, img_b, roi=(20, 20, 236, 236)) + assert isinstance(result, dict) + assert "std_nc_ratio" in result + + +class TestNoiseCorrectedContrastSingleImage: + """Single-image mode: every new NC key must be empty, matching the existing + contrast_ratios_b / weber_contrast_b invariant (never None/absent, never + populated with placeholder values on the A side either).""" + + @pytest.fixture(scope="class") + @classmethod + def single_result(cls, astro_image_a): + return SpatialDetailAnalyzer().analyze(astro_image_a) + + @pytest.mark.parametrize("key", [ + "std_nc_score_a", "std_nc_score_b", "std_nc_ratio", + "log_nc_score_a", "log_nc_score_b", "log_nc_ratio", + "wavelet_nc_score_a", "wavelet_nc_score_b", "wavelet_nc_ratio", + "weber_nc_score_a", "weber_nc_score_b", "weber_nc_ratio", + "gm_nc_score_a", "gm_nc_score_b", "gm_nc_ratio", + ]) + def test_nc_key_empty_in_single_image_mode(self, single_result, key): + assert not single_result[key] + + def test_shared_nebula_pixels_zero(self, single_result): + assert single_result["nc_shared_nebula_pixels"] == 0 + + def test_no_normalized_panels_in_single_image_mode(self, single_result): + assert not any(k.startswith("nrm_") for k in single_result["panels"]) + + +class TestNcScoreHelper: + """Direct unit tests of _nc_score's mask-emptiness contract — deterministic, + unlike relying on two real images happening to produce non-overlapping + nebula masks (background-threshold noise can create incidental overlap).""" + + def test_none_mask_neb_shared_returns_none(self): + analyzer = SpatialDetailAnalyzer() + detail = np.ones((20, 20), dtype=np.float32) + bg_mask = np.ones((20, 20), dtype=bool) + score, noise = analyzer._nc_score(detail, None, bg_mask) + assert score is None and noise is None + + def test_empty_shared_nebula_mask_returns_none(self): + analyzer = SpatialDetailAnalyzer() + detail = np.ones((20, 20), dtype=np.float32) + mask_neb_shared = np.zeros((20, 20), dtype=bool) # no shared nebula pixels + bg_mask = np.ones((20, 20), dtype=bool) + score, noise = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None + + def test_empty_bg_mask_returns_none(self): + analyzer = SpatialDetailAnalyzer() + detail = np.ones((20, 20), dtype=np.float32) + mask_neb_shared = np.ones((20, 20), dtype=bool) + bg_mask = np.zeros((20, 20), dtype=bool) # no background pixels + score, noise = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None + + def test_zero_noise_floor_returns_none(self): + analyzer = SpatialDetailAnalyzer() + detail = np.zeros((20, 20), dtype=np.float32) + detail[:10, :] = 5.0 # nebula half has signal, background half is exactly 0 + mask_neb_shared = np.zeros((20, 20), dtype=bool) + mask_neb_shared[:10, :] = True + bg_mask = np.zeros((20, 20), dtype=bool) + bg_mask[10:, :] = True + score, noise = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None + + def test_valid_masks_return_ratio(self): + analyzer = SpatialDetailAnalyzer() + detail = np.zeros((20, 20), dtype=np.float32) + detail[:10, :] = 10.0 + detail[10:, :] = 2.0 + mask_neb_shared = np.zeros((20, 20), dtype=bool) + mask_neb_shared[:10, :] = True + bg_mask = np.zeros((20, 20), dtype=bool) + bg_mask[10:, :] = True + score, noise = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score == pytest.approx(5.0) + assert noise == pytest.approx(2.0) diff --git a/tests/test_report/test_inspector_catalog.py b/tests/test_report/test_inspector_catalog.py new file mode 100644 index 0000000..87f574b --- /dev/null +++ b/tests/test_report/test_inspector_catalog.py @@ -0,0 +1,93 @@ +"""Regression tests for the Report Inspector's .npz panel catalog. + +Covers the bug where _write_inspector_file used a hardcoded, stale panel-name +map (std_15px/std_31px from an old STD_KERNEL_SIZES, no Weber entries at all) +that silently omitted panels SpatialDetailAnalyzer actually computes. The fix +derives the catalog dynamically from panels_a.keys() via _panel_display_name. +""" +from __future__ import annotations + +import json + +import numpy as np +import pytest +from astropy.io import fits as ap_fits + +from analysis.image_filters import SpatialDetailAnalyzer +from core.astro_image import AstroImage +from core.models import AnalysisResult +from report.report_builder import ReportBuilder + + +def _make_two_image_pair(tmp_path_factory) -> tuple[AstroImage, AstroImage]: + out = tmp_path_factory.mktemp("inspector_fits") + rng_a = np.random.default_rng(1) + rng_b = np.random.default_rng(2) + h, w = 256, 256 + for path, rng in ((out / "a.fits", rng_a), (out / "b.fits", rng_b)): + data = rng.normal(1000.0, 20.0, (h, w)).astype(np.float32) + yy, xx = np.mgrid[0:h, 0:w] + blob = 500.0 * np.exp(-0.5 * (((xx - w / 2) ** 2 + (yy - h / 2) ** 2) / 30.0 ** 2)) + data += blob.astype(np.float32) + hdr = ap_fits.Header() + hdr["EGAIN"] = 1.0 + hdr["GAIN"] = 1.0 + ap_fits.writeto(str(path), data, hdr, overwrite=True) + img_a = AstroImage(str(out / "a.fits"), label="A") + img_a.load() + img_a.estimate_background() + img_b = AstroImage(str(out / "b.fits"), label="B") + img_b.load() + img_b.estimate_background() + return img_a, img_b + + +@pytest.fixture(scope="module") +def inspector_npz_path(tmp_path_factory): + img_a, img_b = _make_two_image_pair(tmp_path_factory) + spatial = SpatialDetailAnalyzer().analyze(img_a, img_b) + result_a = AnalysisResult(label="A", spatial_metrics=spatial) + result_b = AnalysisResult(label="B", spatial_metrics=spatial) + + out_path = tmp_path_factory.mktemp("inspector_out") / "report_inspector.npz" + ReportBuilder()._write_inspector_file(out_path, img_a, img_b, result_a, result_b) + return out_path, spatial + + +class TestInspectorCatalogDynamicPanels: + def test_every_computed_panel_is_cataloged(self, inspector_npz_path): + path, spatial = inspector_npz_path + npz = np.load(str(path), allow_pickle=False) + catalog = json.loads(npz["catalog_json"].tobytes().decode("utf-8")) + spatial_entries = catalog["sections"].get("Spatial Detail", []) + cataloged_names = {e["name"] for e in spatial_entries} + + panels = spatial["panels"] + assert len(spatial_entries) == len(panels) + for pkey in panels: + from report.report_builder import _panel_display_name + assert _panel_display_name(pkey) in cataloged_names + + def test_weber_panels_present(self, inspector_npz_path): + # Regression: the old hardcoded _PANEL_IMAGE_SETS never included Weber at all. + _, spatial = inspector_npz_path + weber_keys = [k for k in spatial["panels"] if k.startswith("weber_")] + assert weber_keys, "fixture should have produced Weber panels" + + def test_all_std_kernel_sizes_present(self, inspector_npz_path): + # Regression: the old map referenced std_15px/std_31px from a stale + # STD_KERNEL_SIZES and dropped std_3px/std_10px from the current one. + from core.models import STD_KERNEL_SIZES + _, spatial = inspector_npz_path + for ks in STD_KERNEL_SIZES: + assert f"std_{ks}px" in spatial["panels"] + + def test_noise_normalized_panels_cataloged(self, inspector_npz_path): + path, spatial = inspector_npz_path + npz = np.load(str(path), allow_pickle=False) + catalog = json.loads(npz["catalog_json"].tobytes().decode("utf-8")) + spatial_entries = catalog["sections"].get("Spatial Detail", []) + cataloged_names = {e["name"] for e in spatial_entries} + nrm_keys = [k for k in spatial["panels"] if k.startswith("nrm_")] + assert nrm_keys + assert any("(noise-normalized)" in n for n in cataloged_names) diff --git a/tests/test_report/test_report_helpers.py b/tests/test_report/test_report_helpers.py index ab3e7bb..a5f7f16 100644 --- a/tests/test_report/test_report_helpers.py +++ b/tests/test_report/test_report_helpers.py @@ -17,6 +17,8 @@ _fig_to_b64, _psf_stat_test, _power_ratio_db, + _nc_ratio_rows, + _panel_display_name, ) EM_DASH = "—" @@ -228,3 +230,59 @@ def test_same_length_different_values_returns_none(self): freq_b = np.linspace(0.0, 0.4, 5) rp = np.ones(5) assert _power_ratio_db(freq_a, rp, freq_b, rp) is None + + +class TestNcRatioRows: + def _label(self, scale): + return f"{scale} px" + + def test_known_ratio_formatted(self): + rows = _nc_ratio_rows({3: 2.0}, {3: 1.0}, {3: 2.0}, self._label) + assert "2.000" in rows + assert "1.000" in rows + + def test_missing_entry_renders_em_dash(self): + rows = _nc_ratio_rows({3: 1.5}, {}, {}, self._label) + assert EM_DASH in rows + + def test_none_ratio_renders_em_dash(self): + rows = _nc_ratio_rows({3: 1.5}, {3: None}, {3: None}, self._label) + assert EM_DASH in rows + + def test_empty_inputs_produce_no_rows(self): + rows = _nc_ratio_rows({}, {}, {}, self._label) + assert rows == "" + + def test_scale_label_used(self): + rows = _nc_ratio_rows({5: 1.0}, {5: 1.0}, {5: 1.0}, lambda s: f"σ = {s} px") + assert "σ = 5 px" in rows + + def test_ratio_cell_colored_independently_of_ab_columns(self): + # A/B columns: A < B -> "worse"/"better"; ratio itself (0.5) vs parity (1.0) + # is colored via a SEPARATE _better_worse_class(vr, 1.0) call. + rows = _nc_ratio_rows({3: 1.0}, {3: 2.0}, {3: 0.5}, self._label) + assert "worse" in rows # A's own value vs B's + # the ratio 0.5 < 1.0 is also "worse" here, but computed independently + # (not derived from the A/B column classes) - verify it renders at all + assert "0.500" in rows + + def test_union_of_scale_keys(self): + rows = _nc_ratio_rows({3: 1.0}, {5: 1.0}, {}, self._label) + assert "3 px" in rows + assert "5 px" in rows + + +class TestPanelDisplayName: + @pytest.mark.parametrize("pkey,expected", [ + ("std_3px", "Std Dev 3 px"), + ("std_10px", "Std Dev 10 px"), + ("weber_9px", "Weber 9 px"), + ("log_1.5", "LoG σ 1.5 px"), + ("wavelet_2", "Wavelet level 2"), + ("gradient_3.0", "Gradient σ 3.0 px"), + ("nrm_log_1.5", "LoG σ 1.5 px (noise-normalized)"), + ("nrm_std_3px", "Std Dev 3 px (noise-normalized)"), + ("nrm_gradient_6.0", "Gradient σ 6.0 px (noise-normalized)"), + ]) + def test_display_name(self, pkey, expected): + assert _panel_display_name(pkey) == expected From a46631404087a6b35d28b44fa230fb2e751880cb Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Wed, 8 Jul 2026 22:32:10 -0400 Subject: [PATCH 13/14] feat: export original image ROI as a Spatial Detail inspector panel Adds the exact mean-normalised, ROI-cropped array every std/LoG/wavelet/ Weber/gradient map is computed from as its own "Original Image (ROI)" panel, so the Report Inspector can show the source image content alongside a map, pixel-aligned to the same crop. Co-Authored-By: Claude Sonnet 5 --- analysis/image_filters.py | 10 ++++++++++ report/report_builder.py | 4 +++- tests/test_analysis/test_spatial_detail.py | 14 ++++++++++++++ tests/test_report/test_report_helpers.py | 1 + 4 files changed, 28 insertions(+), 1 deletion(-) diff --git a/analysis/image_filters.py b/analysis/image_filters.py index c44dc81..d699cce 100644 --- a/analysis/image_filters.py +++ b/analysis/image_filters.py @@ -120,6 +120,16 @@ def _clip01(v): return max(0.0, min(1.0, v)) result["display_roi"] = display_roi + # 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. + 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, + } + # 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 diff --git a/report/report_builder.py b/report/report_builder.py index 76fcfbb..7d1d0f6 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -483,7 +483,9 @@ def _panel_display_name(pkey: str) -> str: additions — automatically appears in the Report Inspector catalog.""" nrm = pkey.startswith("nrm_") base = pkey[4:] if nrm else pkey - if base.startswith("std_") and base.endswith("px"): + if base == "original": + name = "Original Image (ROI)" + elif base.startswith("std_") and base.endswith("px"): name = f"Std Dev {base[4:-2]} px" elif base.startswith("log_"): name = f"LoG σ {base[4:]} px" diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py index 0f36851..2a59626 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -27,6 +27,13 @@ def test_panels_present(self, astro_image_a): result = SpatialDetailAnalyzer().analyze(astro_image_a) assert "panels" in result + def test_original_panel_present_single_image(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + original = result["panels"]["original"] + assert original["a"] is not None + assert original["b"] is None + assert original["diff"] is None + def test_contrast_ratios_are_positive(self, astro_image_a): result = SpatialDetailAnalyzer().analyze(astro_image_a) ratios = result.get("contrast_ratios_a") or [] @@ -205,6 +212,13 @@ def test_normalized_panels_present_two_image(self, nc_result): for lvl in (2, 3): assert f"nrm_wavelet_{lvl}" in panels + def test_original_panel_present_two_image(self, nc_result): + original = nc_result["panels"]["original"] + assert original["a"] is not None + assert original["b"] is not None + assert original["diff"] is not None + assert original["a"].shape == original["b"].shape + def test_normalized_panel_values_differ_from_raw(self, nc_result): panels = nc_result["panels"] ks = min(STD_KERNEL_SIZES) diff --git a/tests/test_report/test_report_helpers.py b/tests/test_report/test_report_helpers.py index a5f7f16..2e47abd 100644 --- a/tests/test_report/test_report_helpers.py +++ b/tests/test_report/test_report_helpers.py @@ -283,6 +283,7 @@ class TestPanelDisplayName: ("nrm_log_1.5", "LoG σ 1.5 px (noise-normalized)"), ("nrm_std_3px", "Std Dev 3 px (noise-normalized)"), ("nrm_gradient_6.0", "Gradient σ 6.0 px (noise-normalized)"), + ("original", "Original Image (ROI)"), ]) def test_display_name(self, pkey, expected): assert _panel_display_name(pkey) == expected From 3d80afc270b204cad69263bcda7e1189cf389d79 Mon Sep 17 00:00:00 2001 From: Brent <52629076+brentmantooth@users.noreply.github.com> Date: Wed, 8 Jul 2026 22:32:30 -0400 Subject: [PATCH 14/14] fix: correct ESF rotation-clipping artifact in edge analyzer, add quality gate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _extract_esf() rotated the entire ROI and averaged across all rows to build the edge spread function. scipy.ndimage.rotate(reshape=False) clips the box's corners for any oblique angle (worst at 45deg, ~30% of the area), and that zero-padding was getting mixed into the column averages, fabricating a second, spurious transition in the ESF for otherwise-clean edges. Averaging is now restricted to the disc inscribed in the ROI square, the largest region mathematically guaranteed to contain only genuine data at any rotation angle. A new ESF monotonicity check (_esf_quality) catches edges that are still contaminated (e.g. a scan line that genuinely crosses two physical edges, like a thin filament) — in auto-detect mode these are skipped in favour of the next-strongest gradient peak, and if every candidate fails, the least-bad one is shown but clearly flagged "low confidence" in the report rather than silently reporting a meaningless width. Co-Authored-By: Claude Sonnet 5 --- analysis/edge_analyzer.py | 206 ++++++++++++++++------ core/models.py | 1 + report/report_builder.py | 28 ++- tests/test_analysis/test_edge_analyzer.py | 173 ++++++++++++++++++ 4 files changed, 353 insertions(+), 55 deletions(-) diff --git a/analysis/edge_analyzer.py b/analysis/edge_analyzer.py index 8e6b0ec..e24b58b 100644 --- a/analysis/edge_analyzer.py +++ b/analysis/edge_analyzer.py @@ -1,18 +1,24 @@ from __future__ import annotations +import warnings + import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt -from scipy.ndimage import sobel, rotate, gaussian_gradient_magnitude, median_filter +from scipy.ndimage import sobel, rotate, gaussian_gradient_magnitude, median_filter, uniform_filter1d from scipy.interpolate import interp1d from core.astro_image import AstroImage from core.fig_utils import figs_to_b64 -from core.models import EDGE_ROI_HALF_WIDTH, EDGE_ROI_MAP_INDICATOR_PX, SECTION8_BORDER_CROP_FRACTION +from core.models import (EDGE_ROI_HALF_WIDTH, EDGE_ROI_MAP_INDICATOR_PX, + SECTION8_BORDER_CROP_FRACTION, EDGE_ESF_MIN_MONOTONICITY) EDGE_DISPLAY_HALF_WIDTH = 250 # half-side of the context window shown in the report figure N_TOP_EDGES = 3 # number of gradient peaks to auto-detect +N_CANDIDATE_EDGES = N_TOP_EDGES * 3 # extra auto-detect candidates so low-quality ones can be skipped +_ESF_DISC_MARGIN_PX = 2.0 # safety margin (px) subtracted from the inscribed-circle radius +_ESF_QUALITY_SMOOTH_FRAC = 0.20 # smoothing window as a fraction of ESF length, for the quality metric only def _gradient_sigma(pixel_scale: float) -> float: @@ -54,9 +60,15 @@ def analyze(self, image: AstroImage, image.estimate_background() bgsub = image.background_subtracted() - # Build list of (roi_data_array, roi_tuple) pairs + # Build list of (roi_data_array, roi_tuple) pairs. Auto-detect mode + # searches more candidates than N_TOP_EDGES so low-quality ones (ESF + # crossed more than one physical edge) can be skipped in favour of the + # next-strongest gradient peak; user-drawn/A-matched ROIs are fixed + # and can't be swapped for an alternative, so quality is only flagged + # for those, never used to drop them. + allow_skip = roi is None if roi is None: - roi_pairs = self._auto_detect_top_rois(bgsub, image) + roi_pairs = self._auto_detect_top_rois(bgsub, image, n=N_CANDIDATE_EDGES) elif isinstance(roi, list): h, w = bgsub.shape roi_pairs = [] @@ -74,7 +86,10 @@ def analyze(self, image: AstroImage, rois_used = [r for _, r in roi_pairs] edges: list[dict] = [] + rejected: list[dict] = [] # low-quality auto-detect candidates, kept only as a fallback for i, (roi_data, roi_tuple) in enumerate(roi_pairs): + if allow_skip and len(edges) >= N_TOP_EDGES: + break if roi_data is None or roi_data.size == 0: continue @@ -90,44 +105,22 @@ def analyze(self, image: AstroImage, if esf is None or len(esf) < 5: continue - lsf = self._compute_lsf(positions, esf) - width = self._measure_edge_width(positions, esf) - ecr = self._measure_edge_contrast_ratio(roi_data, edge_info) - - # 500×500 display context - x0, y0, x1, y1 = roi_tuple - xc_full = x0 + edge_info["center_x"] - yc_full = y0 + edge_info["center_y"] - dw = EDGE_DISPLAY_HALF_WIDTH - dx0 = max(0, xc_full - dw) - dy0 = max(0, yc_full - dw) - dx1 = min(bgsub.shape[1], xc_full + dw) - dy1 = min(bgsub.shape[0], yc_full + dw) - display_roi = bgsub[dy0:dy1, dx0:dx1] - analysis_rect = (x0 - dx0, y0 - dy0, x1 - dx0, y1 - dy0) - edge_info_display = dict(edge_info) - edge_info_display["center_x"] = xc_full - dx0 - edge_info_display["center_y"] = yc_full - dy0 - - edges.append({ - "roi_used": roi_tuple, - "gradient_magnitude": edge_info["gradient_magnitude"], - "angle_rad": edge_info["angle_rad"], - "edge_width_10_90_px": width, - "edge_width_10_90_arcsec": (width * image.pixel_scale - if width is not None else None), - "edge_contrast_ratio": ecr, - "esf": esf, - "lsf": lsf, - "positions": positions.tolist(), - "figures": figs_to_b64({ - "edge": self._plot_results( - roi_data, display_roi, analysis_rect, - positions, esf, lsf, width, image.label, - edge_info_display, edge_num=i + 1, - ) - }), - }) + quality = self._esf_quality(esf) + low_confidence = quality < EDGE_ESF_MIN_MONOTONICITY + entry = self._build_edge_entry( + image, bgsub, roi_data, roi_tuple, edge_info, + positions, esf, i, quality, low_confidence) + + if allow_skip and low_confidence: + rejected.append(entry) + continue + edges.append(entry) + + if allow_skip and not edges and rejected: + # Every candidate crossed more than one edge -- surface the least + # bad one rather than showing nothing, clearly flagged. + rejected.sort(key=lambda e: e["esf_quality"], reverse=True) + edges.append(rejected[0]) # Full-image gradient map for report visualisation from core.stretch import stf_stretch @@ -162,13 +155,66 @@ def analyze(self, image: AstroImage, "figures": figs_to_b64({"gradient_map": gradient_fig}), } if edges: - best = max(edges, key=lambda e: e.get("gradient_magnitude") or 0) + # Prefer a confident measurement over a merely-stronger-gradient one + # (a corner/knot can have higher raw gradient than a clean edge). + best = max(edges, key=lambda e: (not e.get("low_confidence", False), + e.get("gradient_magnitude") or 0)) for k in ("edge_width_10_90_px", "edge_width_10_90_arcsec", "gradient_magnitude", "edge_contrast_ratio", "esf", "lsf"): result[k] = best[k] return result + # ------------------------------------------------------------------ + # Edge-result assembly + # ------------------------------------------------------------------ + + def _build_edge_entry(self, image: AstroImage, bgsub: np.ndarray, + roi_data: np.ndarray, roi_tuple: tuple, + edge_info: dict, positions: np.ndarray, esf: np.ndarray, + edge_num: int, quality: float, low_confidence: bool) -> dict: + lsf = self._compute_lsf(positions, esf) + width = self._measure_edge_width(positions, esf) + ecr = self._measure_edge_contrast_ratio(roi_data, edge_info) + + # 500×500 display context + x0, y0, x1, y1 = roi_tuple + xc_full = x0 + edge_info["center_x"] + yc_full = y0 + edge_info["center_y"] + dw = EDGE_DISPLAY_HALF_WIDTH + dx0 = max(0, xc_full - dw) + dy0 = max(0, yc_full - dw) + dx1 = min(bgsub.shape[1], xc_full + dw) + dy1 = min(bgsub.shape[0], yc_full + dw) + display_roi = bgsub[dy0:dy1, dx0:dx1] + analysis_rect = (x0 - dx0, y0 - dy0, x1 - dx0, y1 - dy0) + edge_info_display = dict(edge_info) + edge_info_display["center_x"] = xc_full - dx0 + edge_info_display["center_y"] = yc_full - dy0 + + return { + "roi_used": roi_tuple, + "gradient_magnitude": edge_info["gradient_magnitude"], + "angle_rad": edge_info["angle_rad"], + "edge_width_10_90_px": width, + "edge_width_10_90_arcsec": (width * image.pixel_scale + if width is not None else None), + "edge_contrast_ratio": ecr, + "esf": esf, + "lsf": lsf, + "positions": positions.tolist(), + "esf_quality": quality, + "low_confidence": low_confidence, + "figures": figs_to_b64({ + "edge": self._plot_results( + roi_data, display_roi, analysis_rect, + positions, esf, lsf, width, image.label, + edge_info_display, edge_num=edge_num + 1, + low_confidence=low_confidence, + ) + }), + } + # ------------------------------------------------------------------ # ROI auto-detection # ------------------------------------------------------------------ @@ -238,20 +284,72 @@ def _extract_esf(self, roi_data: np.ndarray, edge_info: dict) -> tuple[np.ndarray, np.ndarray | None]: angle_deg = np.degrees(edge_info["angle_rad"]) rotation_angle = -(90.0 - angle_deg) - rotated = rotate(roi_data, rotation_angle, reshape=False, order=3) - - esf_raw = np.mean(rotated, axis=0) + # cval=nan (not the default 0.0) marks pixels that rotate() had to + # invent because the source square doesn't cover that output pixel at + # this angle -- see the disc-mask comment below for why this matters. + rotated = rotate(roi_data, rotation_angle, reshape=False, order=3, cval=np.nan) + + # Rotating a square ROI about its own center clips its corners for any + # angle other than 0/90 deg (reshape=False keeps the output the same + # size as the input, so content that rotates outside that frame is + # lost). Averaging over the full box (the previous behaviour) mixes + # this invented/zero-filled corner content into the profile, which can + # fabricate a second, spurious transition -- worst at 45 deg, where + # ~30% of the box area is affected. The disc inscribed in the ROI + # square (radius = half the side length) is the largest region + # guaranteed to contain only genuine data at ANY rotation angle, since + # a rotation about the center never moves points closer to the center + # outside the original square. Restrict averaging to that disc (minus + # a small margin for cubic-spline interpolation at the boundary). + h, w = rotated.shape + cy, cx = (h - 1) / 2.0, (w - 1) / 2.0 + yy, xx = np.mgrid[0:h, 0:w] + radius = min(h, w) / 2.0 - _ESF_DISC_MARGIN_PX + in_disc = ((yy - cy) ** 2 + (xx - cx) ** 2) <= radius ** 2 + masked = np.where(in_disc & ~np.isnan(rotated), rotated, np.nan) + + # Columns beyond the disc radius are entirely NaN by construction + # (trimmed out below) -- nanmean's "empty slice" warning for those is + # expected, not a sign of a problem. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + esf_raw = np.nanmean(masked, axis=0) positions = np.arange(len(esf_raw), dtype=float) - lo, hi = esf_raw.min(), esf_raw.max() + valid_idx = np.where(~np.isnan(esf_raw))[0] + if len(valid_idx) < 5: + return positions, None + lo = float(np.nanmin(esf_raw)) + hi = float(np.nanmax(esf_raw)) if hi - lo < 1e-12: return positions, None esf = (esf_raw - lo) / (hi - lo) - if esf[0] > esf[-1]: + if esf[valid_idx[0]] > esf[valid_idx[-1]]: esf = 1.0 - esf - return positions, esf + # Trim to the contiguous valid range (only the disc-excluded columns + # near the two far ends can be NaN) so downstream LSF/width code never + # has to handle missing values. + i0, i1 = valid_idx[0], valid_idx[-1] + 1 + return positions[i0:i1] - positions[i0], esf[i0:i1] + + @staticmethod + def _esf_quality(esf: np.ndarray) -> float: + """Monotonicity ratio of a lightly-smoothed ESF: net variation over + total variation. 1.0 = a clean single transition; low values mean the + profile doubles back on itself (median(|detail|)-style contamination + from a second edge -- a corner, filament, or nearby knot crossed by + the same scan line).""" + n = len(esf) + win = max(3, int(round(n * _ESF_QUALITY_SMOOTH_FRAC)) | 1) + smoothed = uniform_filter1d(esf.astype(float), size=win, mode="nearest") + diffs = np.diff(smoothed) + total_variation = float(np.sum(np.abs(diffs))) + if total_variation < 1e-9: + return 0.0 + net_variation = float(abs(smoothed[-1] - smoothed[0])) + return net_variation / total_variation # ------------------------------------------------------------------ # LSF and width @@ -357,7 +455,8 @@ def _plot_results(self, roi_data: np.ndarray, width: float | None, label: str, edge_info: dict | None = None, - edge_num: int = 1) -> plt.Figure: + edge_num: int = 1, + low_confidence: bool = False) -> plt.Figure: from matplotlib.patches import Rectangle fig, ax = plt.subplots(figsize=(5, 5)) @@ -366,7 +465,12 @@ def _plot_results(self, roi_data: np.ndarray, ax.imshow(display_roi, origin="upper", cmap="gray", aspect="equal", interpolation="nearest", extent=[0, w_disp, h_disp, 0]) - ax.set_title(f"Edge #{edge_num} ROI — {label}") + title = f"Edge #{edge_num} ROI — {label}" + if low_confidence: + title += " ⚠ low confidence" + ax.set_title(title, color="#c0392b") + else: + ax.set_title(title) ax.set_xlabel("X (px)") ax.set_ylabel("Y (px)") diff --git a/core/models.py b/core/models.py index cece8df..b6a006b 100644 --- a/core/models.py +++ b/core/models.py @@ -23,6 +23,7 @@ SEEING_WARN_FWHM_ARCS = 3.0 EDGE_ROI_HALF_WIDTH = 30 EDGE_ROI_MAP_INDICATOR_PX = 500 # px; full width of the ROI indicator box drawn on the gradient magnitude map +EDGE_ESF_MIN_MONOTONICITY = 0.3 # min net/total variation ratio; below this the ESF likely crossed >1 edge (corner/filament) FILTER_THICKNESS_MM = 1.0 # narrowband filter substrate thickness (mm); default for UI GLASS_REFRACTIVE_INDEX = 1.9 # dichroic filter substrate refractive index RDF_BIN_WIDTH = 1.0 # px; annular bin width for RDF mean/std computation diff --git a/report/report_builder.py b/report/report_builder.py index 7d1d0f6..a4dc325 100644 --- a/report/report_builder.py +++ b/report/report_builder.py @@ -3292,7 +3292,13 @@ def _plot_esf_lsf_pair(self, w_label = f" (A: {w_a:.2f} px)" elif w_b is not None: w_label = f" (B: {w_b:.2f} px)" - ax_esf.set_title(f"Edge #{edge_num} ESF — 10–90% width{w_label}", fontsize=9) + low_conf = edge_a.get("low_confidence") or edge_b.get("low_confidence") + title = f"Edge #{edge_num} ESF — 10–90% width{w_label}" + if low_conf: + title += " ⚠ low confidence" + ax_esf.set_title(title, fontsize=9, color="#c0392b") + else: + ax_esf.set_title(title, fontsize=9) ax_esf.set_xlabel("Position (px)") ax_esf.set_ylabel("Normalised intensity") ax_esf.legend(fontsize=8) @@ -3375,8 +3381,16 @@ def _section_edge(self, ra: AnalysisResult, rb: AnalysisResult, if ea_e or eb_e: esf_lsf_fig = self._plot_esf_lsf_pair(ea_e, eb_e, ra.label, rb.label, i + 1) esf_lsf_html = _img_tag(esf_lsf_fig, f"Edge #{i+1} ESF + LSF") + low_conf_badge = ( + '  ⚠ low confidence — ' + 'ESF profile is not cleanly monotonic, likely because the scan line ' + 'crossed more than one physical edge (a corner, knot, or nearby ' + 'filament); treat the width/contrast numbers for this edge with ' + 'caution' + if (ea_e.get("low_confidence") or eb_e.get("low_confidence")) else "" + ) edge_figures_html += ( - f'

Edge #{i+1}

' + f'

Edge #{i+1}{low_conf_badge}

' f'
' + "".join(f'
{img}
' for img in [img_a_i, img_b_i] if img) @@ -3470,8 +3484,14 @@ def _section_edge(self, ra: AnalysisResult, rb: AnalysisResult, 'ESF scan direction is taken from whichever image has the stronger overall ' 'gradient, then applied to both, so the two ESF curves always sample the ' 'same cross-section orientation and are directly comparable. ' - 'The table below shows metrics from the strongest of the three edges; individual ' - 'per-edge figures follow.', + 'Each candidate\'s ESF is checked for monotonicity (it should be a ' + 'single smooth transition); a candidate whose profile doubles back on itself — ' + 'usually because the scan line crossed a corner, knot, or a second nearby edge — ' + 'is skipped in favour of the next-strongest gradient peak. If every candidate for an ' + 'image fails this check, the least-bad one is still shown, clearly marked ' + '⚠ low confidence, rather than silently reporting a meaningless width. ' + 'The table below shows metrics from the strongest confident edge (falling back to the ' + 'strongest overall if none are confident); individual per-edge figures follow.', title="How the edge regions were selected")} {gradient_row} diff --git a/tests/test_analysis/test_edge_analyzer.py b/tests/test_analysis/test_edge_analyzer.py index 676b3ae..bf3636a 100644 --- a/tests/test_analysis/test_edge_analyzer.py +++ b/tests/test_analysis/test_edge_analyzer.py @@ -1,13 +1,38 @@ """Unit tests for analysis/edge_analyzer.py.""" from __future__ import annotations +import numpy as np import pytest +from scipy.ndimage import gaussian_filter from analysis.edge_analyzer import EdgeAnalyzer +from core.models import EDGE_ESF_MIN_MONOTONICITY _RESULT_KEYS = {"edges", "n_edges", "rois_used"} +def _make_clean_edge_roi(angle_deg: float = 30.0, size: int = 60) -> np.ndarray: + """Single straight edge through the box center, background-subtracted + semantics (background ~ 0, signal positive) -- matches real bgsub data.""" + yy, xx = np.mgrid[0:size, 0:size] + c = (size - 1) / 2.0 + theta = np.radians(angle_deg) + d = (xx - c) * np.cos(theta) + (yy - c) * np.sin(theta) + roi = np.where(d > 0, 200.0, 0.0).astype(float) + return gaussian_filter(roi, sigma=1.5) + + +def _make_double_edge_roi(size: int = 60) -> np.ndarray: + """Thin bright stripe crossing the box -- genuinely two edges for any + perpendicular scan direction; no rotation/masking fix can turn this into + a single clean transition.""" + yy, xx = np.mgrid[0:size, 0:size] + c = (size - 1) / 2.0 + stripe = np.abs((xx - c) - 0.3 * (yy - c)) < 5 + roi = np.where(stripe, 200.0, 0.0).astype(float) + return gaussian_filter(roi, sigma=1.5) + + class TestAnalyze: def test_returns_dict(self, astro_image_a): result = EdgeAnalyzer().analyze(astro_image_a) @@ -35,6 +60,154 @@ def test_with_roi(self, astro_image_a): assert isinstance(result, dict) +class TestEsfQuality: + def test_clean_edge_scores_high(self): + ea = EdgeAnalyzer() + roi = _make_clean_edge_roi(angle_deg=30.0) + edge_info = ea._detect_strongest_edge(roi) + _, esf = ea._extract_esf(roi, edge_info) + assert ea._esf_quality(esf) > 0.8 + + def test_double_edge_scores_low(self): + ea = EdgeAnalyzer() + roi = _make_double_edge_roi() + edge_info = ea._detect_strongest_edge(roi) + _, esf = ea._extract_esf(roi, edge_info) + assert ea._esf_quality(esf) < EDGE_ESF_MIN_MONOTONICITY + + def test_perfectly_flat_scores_zero(self): + ea = EdgeAnalyzer() + assert ea._esf_quality(np.full(60, 0.5)) == 0.0 + + @pytest.mark.parametrize("angle_deg", [15.0, 30.0, 60.0, 75.0]) + def test_clean_edge_scores_high_at_various_angles(self, angle_deg): + # Regression: before the disc-mask fix, a clean oblique edge scored as + # low as ~0.0-0.3 (rotate() zero-padding the clipped box corners + # fabricated a second transition), indistinguishable from a genuinely + # bad edge. 45deg is intentionally excluded: a boundary passing exactly + # through both box corners gives Sobel a perfect gradient tie along the + # whole diagonal, so argmax picks a corner-adjacent pixel instead of + # the center -- a narrow, non-representative synthetic degeneracy + # (real edges are never exactly 45.000 deg through both corners) that + # the quality gate correctly flags rather than silently mismeasures. + ea = EdgeAnalyzer() + roi = _make_clean_edge_roi(angle_deg=angle_deg) + edge_info = ea._detect_strongest_edge(roi) + _, esf = ea._extract_esf(roi, edge_info) + assert ea._esf_quality(esf) > 0.8 + + +class TestExtractEsfDiscMask: + def test_no_nan_in_returned_esf(self): + # _extract_esf trims the disc-masked NaN columns internally -- + # callers should never see NaN. + ea = EdgeAnalyzer() + roi = _make_clean_edge_roi(angle_deg=45.0) + edge_info = ea._detect_strongest_edge(roi) + positions, esf = ea._extract_esf(roi, edge_info) + assert esf is not None + assert not np.any(np.isnan(esf)) + assert not np.any(np.isnan(positions)) + + def test_positions_start_at_zero(self): + ea = EdgeAnalyzer() + roi = _make_clean_edge_roi(angle_deg=45.0) + edge_info = ea._detect_strongest_edge(roi) + positions, _ = ea._extract_esf(roi, edge_info) + assert positions[0] == 0.0 + + def test_esf_normalised_to_unit_range(self): + ea = EdgeAnalyzer() + roi = _make_clean_edge_roi(angle_deg=30.0) + edge_info = ea._detect_strongest_edge(roi) + _, esf = ea._extract_esf(roi, edge_info) + assert esf.min() >= -1e-9 + assert esf.max() <= 1.0 + 1e-9 + + +class TestQualityGateAutoDetect: + """Directly control which candidate ROIs _auto_detect_top_rois returns so + the skip/fallback control flow can be tested deterministically, without + depending on real Background2D estimation picking particular peaks.""" + + def test_bad_candidates_skipped_in_favour_of_clean_one(self, astro_image_a, monkeypatch): + bad1 = _make_double_edge_roi() + bad2 = _make_double_edge_roi() + clean = _make_clean_edge_roi(angle_deg=30.0) + candidates = [ + (bad1, (0, 0, 60, 60)), + (bad2, (100, 100, 160, 160)), + (clean, (200, 200, 260, 260)), + ] + monkeypatch.setattr( + EdgeAnalyzer, "_auto_detect_top_rois", + lambda self, bgsub, image, n=9: candidates) + + result = EdgeAnalyzer().analyze(astro_image_a) + assert result["n_edges"] == 1 + assert result["edges"][0]["low_confidence"] is False + assert result["edges"][0]["roi_used"] == (200, 200, 260, 260) + + def test_all_bad_falls_back_to_least_bad_flagged(self, astro_image_a, monkeypatch): + bad1 = _make_double_edge_roi() + bad2 = _make_double_edge_roi() + candidates = [ + (bad1, (0, 0, 60, 60)), + (bad2, (100, 100, 160, 160)), + ] + monkeypatch.setattr( + EdgeAnalyzer, "_auto_detect_top_rois", + lambda self, bgsub, image, n=9: candidates) + + result = EdgeAnalyzer().analyze(astro_image_a) + assert result["n_edges"] == 1 + assert result["edges"][0]["low_confidence"] is True + + def test_edge_entries_have_quality_keys(self, astro_image_a, monkeypatch): + clean = _make_clean_edge_roi(angle_deg=30.0) + monkeypatch.setattr( + EdgeAnalyzer, "_auto_detect_top_rois", + lambda self, bgsub, image, n=9: [(clean, (0, 0, 60, 60))]) + + result = EdgeAnalyzer().analyze(astro_image_a) + assert result["n_edges"] == 1 + entry = result["edges"][0] + assert "esf_quality" in entry + assert "low_confidence" in entry + assert isinstance(entry["esf_quality"], float) + + +class TestQualityGateExplicitRoi: + """A user-drawn or A-matched ROI is fixed -- low quality must still be + flagged, but the edge must not be silently dropped (there's no + alternative candidate to fall back to).""" + + def test_bad_edge_kept_but_flagged_with_explicit_roi(self, tmp_path): + from astropy.io import fits as ap_fits + from core.astro_image import AstroImage + + h = w = 256 + rng = np.random.default_rng(0) + data = rng.normal(1000.0, 20.0, (h, w)).astype(np.float64) + yy, xx = np.mgrid[0:h, 0:w] + stripe = np.abs((xx - 128) - 0.3 * (yy - 128)) < 5 + data[stripe] += 200.0 + data = np.clip(data, 0, 65535).astype(np.float32) + hdr = ap_fits.Header() + hdr["EGAIN"] = 1.0 + hdr["GAIN"] = 1.0 + path = tmp_path / "stripe.fits" + ap_fits.writeto(str(path), data, hdr, overwrite=True) + + img = AstroImage(str(path), label="Stripe") + img.load() + img.estimate_background() + + result = EdgeAnalyzer().analyze(img, roi=(98, 98, 158, 158)) + assert result["n_edges"] == 1 + assert result["edges"][0]["low_confidence"] is True + + class TestAnalyzeCrosshair: def test_degenerate_crosshair_returns_none(self, astro_image_a): # Zero-length crosshair (start == end)