From 2a5ce761535b2aa3fe7fe2a3c8f5142392a51672 Mon Sep 17 00:00:00 2001
From: Brent <52629076+brentmantooth@users.noreply.github.com>
Date: Fri, 17 Jul 2026 10:24:50 -0400
Subject: [PATCH 01/13] Harden astropy/photutils usage: fix float64 leak, unify
gain parsing, dedupe Moffat fitting
background_subtracted() could silently upcast to float64 if photutils ever
returned a float64 background model; Background2D relied on photutils'
internal SigmaClip default and used the now-deprecated bkgrms_estimator
keyword. Camera gain resolution (EGAIN > GAIN > CCDGAIN > GAINDB) existed
in two places that could disagree for the same file - the GUI metadata
panel and the SNR analyzer - now unified behind one shared helper.
The 2D Moffat PSF fit was duplicated between psf_analyzer.py and
halo_dialog.py with drifted plausibility bounds; extracted into
analysis/moffat_fit.py and migrated off the deprecated LevMarLSQFitter to
TRFLSQFitter with explicit bounds. PSF_FWHM_CLIP_NSIGMA was being applied
to a raw (unscaled) MAD, clipping at ~2sigma despite its name; switched to
mad_std so the multiplier means what it says.
Also: guarded unguarded float() header parses, deduped the pixel-scale
keyword table, removed dead code, and pinned astropy>=6.0/photutils>=3.0
in requirements files.
Co-Authored-By: Claude Sonnet 5 ` tags (e.g. "see 8g for…", "(8b–8f, 8g)"). After adding, removing, or renumbering a subsection, `grep` the function for every old *and* new heading letter — HTML renders a stale cross-reference without error, it just silently misdirects the reader to the wrong subsection. |
| Stale ROI crashes Section 8 with "index -1 is out of bounds for axis 0 with size 0" | `MainWindow._roi` (`gui/main_window.py`) is never cleared when a new image is loaded into either panel. If the user draws an ROI on one image pair, then loads a smaller replacement pair without clearing it, the stale coordinates go out of bounds for the new image. NumPy doesn't raise on an out-of-range slice — `norm_a[ry0:ry1, rx0:rx1]` silently returns a zero-size array — so the crash surfaces much later and far from the real cause: `SpatialDetailAnalyzer._plot_mask_illustration → _stretch_for_display → np.percentile(empty_array, ...)`. The same unguarded `bgsub[y0:y1, x0:x1]` pattern exists in `power_spectrum.py::_extract_roi` and `edge_analyzer.py::analyze`, so a stale ROI can corrupt those sections too (with different, equally misleading errors) if they happen to run. Fixed at the single real boundary — `MainWindow._on_run()`, which is the only path that constructs `AnalysisThread` — by validating `self._roi` against every loaded image's `data.shape` right before `settings["roi"]` is set; an out-of-bounds ROI is cleared (falls back to auto-detect/full-image) with a `QMessageBox` explaining why, rather than patching each analyzer's slice individually. |
+| `binary_dilation` on a loose sigma-threshold mask amplifies noise, not signal | Growing a boolean mask straight from a threshold cut (e.g. Section 8's nebula mask at 1.7σ) dilates *every* True pixel, including scattered single/few-pixel noise-driven false positives — expected in bulk at a loose sigma cut (~4.5% of pixels at 1.7σ one-sided). Each isolated speck balloons into a `~(2·dilation_px+1)²`-px blob, inflating the mask area by 4x+ and diluting any signal-vs-background metric computed over it. Fix: strip small isolated connected components (`scipy.ndimage.label` + `np.bincount` size filter, same size threshold used for hole-filling) *before* calling `binary_dilation` — see `SpatialDetailAnalyzer._remove_small_objects` / `_fill_small_holes` in `image_filters.py`. Caught by comparing mask pixel counts with dilation on vs off on real (noisy) fixture data — a clean synthetic square mask (no noise) will not reveal this bug. |
+| Adding GUI parameter rows clips existing text in the Parameters group | `gui/main_window.py`'s `AnalysisControlPanel.setMaximumHeight(...)` caps the whole control panel's height. Metrics / Parameters / Run are laid out side-by-side (`QHBoxLayout` in `control_panel.py::_build_ui`), so the cap must fit the *tallest* group box's natural content height. Adding new `QFormLayout` rows to any group box (e.g. three new spinboxes) grows that box's required height without growing the cap, clipping/compressing every row's text top-and-bottom. Bump `setMaximumHeight(...)` proportionally (roughly `old_cap * new_row_count / old_row_count`) whenever a group box gains rows. |
---
diff --git a/analysis/image_filters.py b/analysis/image_filters.py
index 2d52d87..36da3a6 100644
--- a/analysis/image_filters.py
+++ b/analysis/image_filters.py
@@ -11,7 +11,7 @@
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
-from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter
+from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter, binary_dilation, binary_fill_holes, label
import pywt
from core.astro_image import AstroImage
@@ -20,7 +20,9 @@
WEBER_KERNEL_SIZES,
XS_LINE_ALPHA, SECTION8_BORDER_CROP_FRACTION, SECTION8_ANALYSIS_CMAP,
XS_SNR_REGION_WIDTH, SECTION8_DIFF_DIST_MAX_SAMPLES,
- SECTION8_LOGRATIO_EPS_PERCENTILE, SECTION8_SCATTER_MAX_SAMPLES)
+ SECTION8_LOGRATIO_EPS_PERCENTILE, SECTION8_SCATTER_MAX_SAMPLES,
+ SECTION8_NEBULA_MASK_SIGMA, SECTION8_NEBULA_MASK_DILATION_PX,
+ SECTION8_NEBULA_MASK_MAX_HOLE_PX)
MAX_DIM_FOR_STD = 2048 # downsample to this before generic_filter (performance)
_DISPLAY_SMOOTH_SIGMA = 1.0 # applied to maps before plotting; does NOT affect metrics
@@ -41,7 +43,10 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None,
weber_kernel_sizes: tuple = WEBER_KERNEL_SIZES,
crosshair: dict | None = None,
roi: tuple | None = None,
- xs_snr_width: int | None = None) -> dict:
+ xs_snr_width: int | None = None,
+ nebula_sigma: float = SECTION8_NEBULA_MASK_SIGMA,
+ nebula_dilation_px: int = SECTION8_NEBULA_MASK_DILATION_PX,
+ nebula_max_hole_px: int = SECTION8_NEBULA_MASK_MAX_HOLE_PX) -> dict:
image_a.estimate_background()
if image_b is not None:
@@ -57,6 +62,9 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None,
"wavelet_snr_b": {},
"sigma_noise_a": None,
"sigma_noise_b": None,
+ "nebula_sigma": nebula_sigma,
+ "nebula_dilation_px": nebula_dilation_px,
+ "nebula_max_hole_px": nebula_max_hole_px,
"weber_contrast_a": {},
"weber_contrast_b": {},
"panels": {},
@@ -80,10 +88,10 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None,
return result
# Nebula / background masks (used for contrast ratio)
- mask_neb_a, mask_bg_a = self._make_masks(image_a)
+ mask_neb_a, mask_bg_a = self._make_masks(image_a, nebula_sigma, nebula_dilation_px, nebula_max_hole_px)
mask_neb_b = mask_bg_b = None
if image_b is not None:
- mask_neb_b, mask_bg_b = self._make_masks(image_b)
+ mask_neb_b, mask_bg_b = self._make_masks(image_b, nebula_sigma, nebula_dilation_px, nebula_max_hole_px)
# When a user ROI is provided, crop the analysis arrays to that region after
# normalisation so the global signal mean is used, not the ROI's local mean.
@@ -124,13 +132,19 @@ def _clip01(v): return max(0.0, min(1.0, v))
result["display_roi"] = display_roi
# Shared nebula/background regions for noise-corrected A/B scoring and diff
- # distributions: pixels BOTH images independently classify the same way.
+ # distributions. Nebula is the two-image UNION: a pixel counts as Nebula if
+ # either image independently classifies it that way, so nebula signal that's
+ # marginal in one image (registration offset, PSF, local noise) still counts.
+ # Background stays the two-image INTERSECTION: a pixel counts as Background
+ # only if both images agree, keeping the noise-floor reference population
+ # clean. Per-image nebula-dominance in _make_masks (bg_mask excludes
+ # nebula_mask) guarantees these two combinations never overlap.
# None in single-image mode.
mask_neb_shared = mask_bg_shared = None
if mask_neb_b is not None:
h_s = min(mask_neb_a.shape[0], mask_neb_b.shape[0])
w_s = min(mask_neb_a.shape[1], mask_neb_b.shape[1])
- mask_neb_shared = mask_neb_a[:h_s, :w_s] & mask_neb_b[:h_s, :w_s]
+ mask_neb_shared = mask_neb_a[:h_s, :w_s] | mask_neb_b[:h_s, :w_s]
mask_bg_shared = mask_bg_a[:h_s, :w_s] & mask_bg_b[:h_s, :w_s]
result["nc_shared_nebula_pixels"] = (
int(np.count_nonzero(mask_neb_shared)) if mask_neb_shared is not None else 0
@@ -313,7 +327,10 @@ def _normalise(self, image: AstroImage) -> np.ndarray | None:
return None
return bgsub / mean_signal
- def _make_masks(self, image: AstroImage) -> tuple[np.ndarray, np.ndarray]:
+ def _make_masks(self, image: AstroImage,
+ nebula_sigma: float = SECTION8_NEBULA_MASK_SIGMA,
+ dilation_px: int = SECTION8_NEBULA_MASK_DILATION_PX,
+ max_hole_px: int = SECTION8_NEBULA_MASK_MAX_HOLE_PX) -> tuple[np.ndarray, np.ndarray]:
rms = image.background_rms
if rms is None:
rms_val = float(np.std(image.background_subtracted()))
@@ -321,7 +338,7 @@ def _make_masks(self, image: AstroImage) -> tuple[np.ndarray, np.ndarray]:
rms_val = float(np.median(rms))
bgsub = image.background_subtracted()
- nebula_mask = bgsub > 2.0 * rms_val
+ nebula_mask = bgsub > nebula_sigma * rms_val
bg_mask = bgsub < 0.5 * rms_val
# Fallback: use top-5% as nebula if no pixels pass threshold
@@ -329,8 +346,57 @@ def _make_masks(self, image: AstroImage) -> tuple[np.ndarray, np.ndarray]:
threshold = np.percentile(bgsub, 95)
nebula_mask = bgsub >= threshold
+ # Fill small enclosed background gaps, strip small isolated noise-driven
+ # specks (same size threshold — a scattered 1-2px false positive at this
+ # sigma level would otherwise balloon into a ~(2*dilation_px+1)^2 blob per
+ # speck once dilated, polluting blank-sky area far from any real nebula),
+ # then grow into adjacent dim/dark transition regions at nebula edges.
+ # All three are applied before the bg_mask exclusion below, since any of
+ # them can pull previously bg-classified pixels into the nebula mask.
+ nebula_mask = self._fill_small_holes(nebula_mask, max_hole_px)
+ nebula_mask = self._remove_small_objects(nebula_mask, max_hole_px)
+ if dilation_px > 0:
+ nebula_mask = binary_dilation(nebula_mask, iterations=dilation_px)
+
+ # Nebula dominates: keep the two masks mutually exclusive after growth.
+ bg_mask = bg_mask & ~nebula_mask
+
return nebula_mask, bg_mask
+ def _fill_small_holes(self, mask: np.ndarray, max_hole_px: int) -> np.ndarray:
+ """Fill enclosed background gaps inside mask up to (max_hole_px)**2 area."""
+ if max_hole_px <= 0:
+ return mask
+ filled = binary_fill_holes(mask)
+ holes = filled & ~mask
+ labeled, n_holes = label(holes)
+ if n_holes == 0:
+ return mask
+ sizes = np.bincount(labeled.ravel())
+ max_area = max_hole_px * max_hole_px
+ keep = np.zeros(sizes.size, dtype=bool)
+ keep[1:] = sizes[1:] <= max_area # label 0 is background, not a hole
+ return mask | keep[labeled]
+
+ def _remove_small_objects(self, mask: np.ndarray, max_size_px: int) -> np.ndarray:
+ """Strip isolated foreground specks up to (max_size_px)**2 area.
+
+ Scattered single/few-pixel noise excursions above the nebula sigma
+ threshold are common at a loose (~1.7 sigma) cut. Left in place, dilation
+ would inflate each one into a much larger blob far from any real nebula
+ structure, so small islands are dropped before growth is applied.
+ """
+ if max_size_px <= 0:
+ return mask
+ labeled, n_objects = label(mask)
+ if n_objects == 0:
+ return mask
+ sizes = np.bincount(labeled.ravel())
+ max_area = max_size_px * max_size_px
+ keep = np.zeros(sizes.size, dtype=bool)
+ keep[1:] = sizes[1:] > max_area # label 0 is background; drop small islands
+ return keep[labeled]
+
def _nebula_bounding_box(self, mask: np.ndarray,
shape: tuple) -> tuple[int, int, int, int] | None:
"""Return (r0, r1, c0, c1) bounding box of the nebula mask with 5% padding.
@@ -506,7 +572,7 @@ def _nc_score(self, detail_map: np.ndarray,
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
+ score = median(|detail|) over the pixels either image classifies 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,
@@ -1256,7 +1322,8 @@ def _plot_mask_illustration(self, base: np.ndarray, mask_neb: np.ndarray,
"Background": "tomato"}) so the two figures read as one visual language.
Unclassified pixels (neither mask) are left plain grayscale. mask_neb/
mask_bg are expected to already be the two-image shared classification
- (mask_neb_a & mask_neb_b) — the exact masks that feed the violin plots
+ (mask_neb_a | mask_neb_b for Nebula, mask_bg_a & mask_bg_b for Background)
+ — the exact masks that feed the violin plots
and correlation scatter plots — so this figure depicts what those plots
are actually gated on.
"""
diff --git a/core/models.py b/core/models.py
index 1289bc9..2447c54 100644
--- a/core/models.py
+++ b/core/models.py
@@ -43,6 +43,9 @@
SECTION8_DIFF_DIST_MAX_SAMPLES = 1000000 # per masked population, per scale — caps violin/KDE cost on full-res diff maps
SECTION8_LOGRATIO_EPS_PERCENTILE = 1.0 # percentile of pooled positive |A|,|B| values used as the epsilon floor in log10(|A|/|B|)
SECTION8_SCATTER_MAX_SAMPLES = 50000 # per masked population, per scale — caps render cost of Section 8g correlation scatter plots
+SECTION8_NEBULA_MASK_SIGMA = 1.7 # ×RMS above background = "Nebula" pixel classification (Section 8 masks); background cut stays fixed at 0.5×RMS
+SECTION8_NEBULA_MASK_DILATION_PX = 3 # px; scipy.ndimage.binary_dilation iterations to grow the nebula mask into adjacent dim/dark nebula regions
+SECTION8_NEBULA_MASK_MAX_HOLE_PX = 5 # px; enclosed background gaps up to this many pixels per side (area ≤ N²) inside the nebula mask are filled before dilation
PSF_SPATIAL_MAP_SIZE = 150 # px; long-axis resolution of FWHM / eccentricity spatial maps
PSF_SPATIAL_MAP_SMOOTH_SIGMA = 5.0 # Gaussian smoothing sigma (px) applied to spatial maps before display
diff --git a/gui/analysis_thread.py b/gui/analysis_thread.py
index c8f9c56..88ea71a 100644
--- a/gui/analysis_thread.py
+++ b/gui/analysis_thread.py
@@ -11,7 +11,8 @@
from PyQt6.QtCore import QThread, pyqtSignal
from core.astro_image import AstroImage
-from core.models import AnalysisResult, LABEL_MAX_LEN
+from core.models import (AnalysisResult, LABEL_MAX_LEN, SECTION8_NEBULA_MASK_SIGMA,
+ SECTION8_NEBULA_MASK_DILATION_PX, SECTION8_NEBULA_MASK_MAX_HOLE_PX)
from analysis.psf_analyzer import PSFAnalyzer
from analysis.halo_analyzer import HaloAnalyzer
from analysis.edge_analyzer import EdgeAnalyzer
@@ -245,13 +246,18 @@ def _power(ps_a=self._starless_a or img_a, ps_b=_ps_b_src):
wavelet_levels = s.get("wavelet_levels", 4)
crosshair = s.get("crosshair")
xs_snr_width = s.get("xs_snr_width")
+ nebula_sigma = s.get("nebula_sigma", SECTION8_NEBULA_MASK_SIGMA)
+ nebula_dilation_px = s.get("nebula_dilation_px", SECTION8_NEBULA_MASK_DILATION_PX)
+ nebula_max_hole_px = s.get("nebula_max_hole_px", SECTION8_NEBULA_MASK_MAX_HOLE_PX)
_sd_b_src = (self._starless_b or img_b) if img_b is not None else None
def _spatial(sd_a=self._starless_a or img_a, sd_b=_sd_b_src,
_ch=crosshair, _roi=roi, _xs_snr_width=xs_snr_width):
sda = SpatialDetailAnalyzer()
spatial = sda.analyze(sd_a, sd_b, levels=wavelet_levels, crosshair=_ch, roi=_roi,
- xs_snr_width=_xs_snr_width)
+ xs_snr_width=_xs_snr_width, nebula_sigma=nebula_sigma,
+ nebula_dilation_px=nebula_dilation_px,
+ nebula_max_hole_px=nebula_max_hole_px)
spatial["used_starless_a"] = self._starless_a is not None
spatial["used_starless_b"] = self._starless_b is not None
result_a.spatial_metrics = spatial
diff --git a/gui/control_panel.py b/gui/control_panel.py
index 6535d9f..34a15a1 100644
--- a/gui/control_panel.py
+++ b/gui/control_panel.py
@@ -14,7 +14,8 @@
from core.models import (
STD_KERNEL_SIZES, LOG_SIGMAS, WAVELET_LEVELS, DEFAULT_PIXEL_SCALE,
MIN_STAR_SNR, SEEING_WARN_FWHM_ARCS, REF_SEEING_ARCSEC, XS_SNR_REGION_WIDTH,
- EPSF_MAX_STARS,
+ EPSF_MAX_STARS, SECTION8_NEBULA_MASK_SIGMA, SECTION8_NEBULA_MASK_DILATION_PX,
+ SECTION8_NEBULA_MASK_MAX_HOLE_PX,
)
@@ -174,6 +175,29 @@ def _build_ui(self) -> None:
self._wavelet_levels.setToolTip("Number of wavelet layers used in spatial detail analysis.")
params_layout.addRow("Wavelet levels:", self._wavelet_levels)
+ self._nebula_sigma = QDoubleSpinBox()
+ self._nebula_sigma.setRange(0.5, 5.0)
+ self._nebula_sigma.setSingleStep(0.1)
+ self._nebula_sigma.setDecimals(2)
+ self._nebula_sigma.setValue(SECTION8_NEBULA_MASK_SIGMA)
+ self._nebula_sigma.setToolTip("Nebula mask threshold for Section 8 spatial detail analysis.\n"
+ "Pixels above this many RMS units over background are classified as Nebula.")
+ params_layout.addRow("Nebula mask threshold (× RMS):", self._nebula_sigma)
+
+ self._nebula_dilation_px = QSpinBox()
+ self._nebula_dilation_px.setRange(0, 20)
+ self._nebula_dilation_px.setValue(SECTION8_NEBULA_MASK_DILATION_PX)
+ self._nebula_dilation_px.setToolTip("Grows the Section 8 nebula mask outward by this many pixels\n"
+ "to capture dim/dark transition regions at nebula edges.")
+ params_layout.addRow("Nebula mask dilation (px):", self._nebula_dilation_px)
+
+ self._nebula_max_hole_px = QSpinBox()
+ self._nebula_max_hole_px.setRange(0, 20)
+ self._nebula_max_hole_px.setValue(SECTION8_NEBULA_MASK_MAX_HOLE_PX)
+ self._nebula_max_hole_px.setToolTip("Fills enclosed background gaps up to this size (px per side)\n"
+ "inside the Section 8 nebula mask, before dilation.")
+ params_layout.addRow("Nebula mask hole-fill (px):", self._nebula_max_hole_px)
+
self._pixel_scale_override = QDoubleSpinBox()
self._pixel_scale_override.setRange(0.0, 20.0)
self._pixel_scale_override.setDecimals(3)
@@ -373,6 +397,9 @@ def settings(self) -> dict:
"pixel_scale_override": pso if pso > 0 else None,
"wavelet_levels": self._wavelet_levels.value(),
"xs_snr_width": self._xs_snr_width.value(),
+ "nebula_sigma": self._nebula_sigma.value(),
+ "nebula_dilation_px": self._nebula_dilation_px.value(),
+ "nebula_max_hole_px": self._nebula_max_hole_px.value(),
"ref_seeing_arcsec": self._ref_seeing_arcsec.value(),
"epsf_max_stars": self._epsf_max_stars.value(),
"roi": self._roi,
diff --git a/gui/main_window.py b/gui/main_window.py
index 5b5cf7f..45b2959 100644
--- a/gui/main_window.py
+++ b/gui/main_window.py
@@ -49,7 +49,7 @@ def _build_ui(self) -> None:
# Control panel below images
self._control = AnalysisControlPanel()
- self._control.setMaximumHeight(240)
+ self._control.setMaximumHeight(340)
main_layout.addWidget(self._control)
# Wire signals
diff --git a/report/report_builder.py b/report/report_builder.py
index b6d9c91..865dad1 100644
--- a/report/report_builder.py
+++ b/report/report_builder.py
@@ -19,7 +19,9 @@
from core.models import (AnalysisResult, HALO_FIT_RADIUS_PX, XS_LINE_ALPHA, GLASS_REFRACTIVE_INDEX,
PSF_SPATIAL_MAP_SIZE, PSF_SPATIAL_MAP_SMOOTH_SIGMA, EDGE_ROI_MAP_INDICATOR_PX,
LABEL_MAX_LEN, REF_SEEING_ARCSEC, REF_SEEING_BETA,
- ABERRATION_MIN_STARS, ABERRATION_OUTER_RADIUS_FRAC)
+ ABERRATION_MIN_STARS, ABERRATION_OUTER_RADIUS_FRAC,
+ SECTION8_NEBULA_MASK_SIGMA, SECTION8_NEBULA_MASK_DILATION_PX,
+ SECTION8_NEBULA_MASK_MAX_HOLE_PX)
from core.astro_image import AstroImage
_TEST_IMAGE_PATH = Path(__file__).parent.parent / "resources" / "ContrastTestImage.png"
@@ -4013,6 +4015,9 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str:
mask_fig = figs.get("mask_illustration")
dist_img, dist_caption = _spatial_diff_distributions_figure(sm.get("diff_dist", {}))
+ nebula_sigma = sm.get("nebula_sigma", SECTION8_NEBULA_MASK_SIGMA)
+ nebula_dilation_px = sm.get("nebula_dilation_px", SECTION8_NEBULA_MASK_DILATION_PX)
+ nebula_max_hole_px = sm.get("nebula_max_hole_px", SECTION8_NEBULA_MASK_MAX_HOLE_PX)
mask_html = (
"
Nebula / Background Mask Regions
" +
_hires_img_tag(mask_fig, "Mask illustration") +
@@ -4020,16 +4025,21 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str:
"How the masks are detected. Each image is independently classified "
"per-pixel from its own background-subtracted flux and noise RMS: a pixel is "
"Nebula if its background-subtracted "
- "value exceeds 2× the image's RMS, and "
+ f"value exceeds {nebula_sigma:g}× the image's RMS, and "
"Background if it's below 0.5× the "
"RMS (a top-5% percentile fallback applies if no pixel clears the nebula threshold). "
"Pixels between the two thresholds are left unclassified (plain grayscale). "
+ f"Enclosed background gaps up to {nebula_max_hole_px}×{nebula_max_hole_px} px "
+ "inside the nebula mask are filled, then the mask is grown outward by "
+ f"{nebula_dilation_px} px to capture dim transition regions at nebula edges, before "
+ "the two images' masks are combined below. "
"How they're used below. The log-ratio distributions and the per-scale "
- "correlation plots embedded in 8b–8f use the two-image intersection "
- "of these masks — "
- "a pixel counts as Nebula only if both images classify it as Nebula (same for "
- "Background). This is deliberately conservative: it excludes pixels where one "
- f"image's registration, PSF, or local noise disagrees with the other's. Shown on "
+ "correlation plots embedded in 8b–8f use the two-image union for Nebula "
+ "— a pixel counts as Nebula if either image classifies it that way, since "
+ "a pixel that's clearly nebula in one image but marginal in the other (e.g. due to "
+ "registration offset, PSF, or local noise) should still count as signal. Background "
+ "stays the two-image intersection — a pixel counts as Background only if "
+ "both images agree, keeping the noise-floor reference population clean. Shown on "
f"{ra.label}, the same array as the "Original" row below."
"
| Metric | Kernel / scale | Primarily measures | Responds to | ||
|---|---|---|---|---|---|
| Local σ map (8b) | 3, 5, 10 px window | Detail (texture / variability) | ' + '|||
| Local σ map (8g) | 3, 5, 10 px window | Detail (texture / variability) | ' "Any local brightness variation — filaments, halos, and noise (can't tell them " ' apart alone) | ||
| Contrast ratio (8b) | same kernel sizes | Contrast, built from the σ map | ' + '|||
| Contrast ratio (8g) | same kernel sizes | Contrast, built from the σ map | ' 'How much more textured the nebula is than blank sky, at this scale | ||
| |LoG| map (8c) | σ = 1.5, 3, 6 px | Detail (edge / curvature strength) | ' + '|||
| |LoG| map (8d) | σ = 1.5, 3, 6 px | Detail (edge / curvature strength) | ' 'Intensity boundaries — filament edges, shell rims — surviving smoothing at scale σ | ||
| Gradient magnitude (8f) | same σ as LoG | Detail (edge sharpness) | ' "How abrupt a boundary is at scale σ; complements Section 6's precise per-edge " ' resolution measurement | ||
| Wavelet decomposition (8d) | levels ≈ 2, 4, 8, 16 px | ' + '||||
| Wavelet decomposition (8e) | levels ≈ 2, 4, 8, 16 px | ' 'Detail by scale band, plus explicit SNR | ' 'Structure whose size matches this scale band; level 1 is noise-only, used to ' ' calibrate the noise floor | ||
| Weber fraction contrast (8e) | 3, 5, 9 px window | " + "||||
| Weber fraction contrast (8h) | 3, 5, 9 px window | " "Contrast (Weber's law, ΔL/L) | " 'Local range vs. local median background — the most literal "Contrast" metric in ' ' this section | ||
| Noise-corrected (NC) score (8b–8f, 8g) | same scale as parent method | ' + '||||
| Noise-corrected (NC) score (8d–8h, 8i) | same scale as parent method | ' 'SNR of Detail | ' "Whether a detail-map response in the nebula is real structure or just this image's " ' own noise floor at that scale |
| Kernel size | {ra.label} | {rb.label} |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
Local σ maps at each kernel size (shared colour scale): Image A (top-left), -Image B (top-right), log-ratio map (middle-left) highlighting where one filter preserves more -local variation, and — when a cross-section line is set — its profile (middle-right), plus a -bottom-row histogram of the log-ratio map's pixel values (same colour scale). Each map -is immediately followed by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{figs_for("nrm_std_")} -Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{figs_for("nrm_log_")} -Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{figs_for("nrm_wavelet_")} +{_family_nrm_figs(_wavelet_rows)} + +| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale per +figure): Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a +cross-section line is set — its profile (middle-right), plus a bottom-row histogram of the +log-ratio map's pixel values (same colour scale). A filter preserving sharper +boundaries shows brighter, more defined gradient response. Each map is immediately followed +by its per-pixel correlation scatter (see methodology above).
+Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
+{_family_nrm_figs(_gradient_rows)} -| Kernel size | {ra.label} | {rb.label} |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
Local σ maps at each kernel size (shared colour scale): Image A (top-left), +Image B (top-right), log-ratio map (middle-left) highlighting where one filter preserves more +local variation, and — when a cross-section line is set — its profile (middle-right), plus a +bottom-row histogram of the log-ratio map's pixel values (same colour scale). Each map +is immediately followed by its per-pixel correlation scatter (see methodology above).
+Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
+{_family_nrm_figs(_std_rows)} +Formula: c = ΔL / L, ' 'where ΔL = Imax − Imin (local range in the K × K kernel) ' @@ -4394,40 +4453,12 @@ def _family_figs_with_corr(rows, map_key_fn) -> str: over dark-sky regions are expected; use a nebula ROI for meaningful filter comparison. Each map is immediately followed by its per-pixel correlation scatter (see methodology above).
Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{figs_for("nrm_weber_")} - -| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale per -figure): Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a -cross-section line is set — its profile (middle-right), plus a bottom-row histogram of the -log-ratio map's pixel values (same colour scale). A filter preserving sharper -boundaries shows brighter, more defined gradient response. Each map is immediately followed -by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{figs_for("nrm_gradient_")} +{_family_nrm_figs(_weber_rows)} -Ratio A/B for every noise-corrected method plotted against its -approximate spatial scale. Scale units differ by method (see 8b–8f methodology +approximate spatial scale. Scale units differ by method (see 8d–8h methodology boxes) — use this chart to spot which spatial-scale regime favours which filter, not to compare absolute ratio values across methods.
""" diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py index 86ef454..081f372 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -33,6 +33,8 @@ def test_original_panel_present_single_image(self, astro_image_a): assert original["a"] is not None assert original["b"] is None assert original["diff"] is None + assert "original" in result["figures"] + assert "corr_original" not in result["figures"] def test_contrast_ratios_are_positive(self, astro_image_a): result = SpatialDetailAnalyzer().analyze(astro_image_a) @@ -218,6 +220,8 @@ def test_original_panel_present_two_image(self, nc_result): assert original["b"] is not None assert original["diff"] is not None assert original["a"].shape == original["b"].shape + assert "original" in nc_result["figures"] + assert "corr_original" in nc_result["figures"] def test_normalized_panel_values_differ_from_raw(self, nc_result): panels = nc_result["panels"] @@ -532,7 +536,8 @@ def test_present_in_two_image_mode(self, nc_result): class TestCorrelationScatterFigures: _CORR_KEYS = ( - [f"corr_std_{ks}px" for ks in STD_KERNEL_SIZES] + ["corr_original"] + + [f"corr_std_{ks}px" for ks in STD_KERNEL_SIZES] + [f"corr_log_{s}" for s in LOG_SIGMAS] + [f"corr_gradient_{s}" for s in LOG_SIGMAS] + ["corr_wavelet_2", "corr_wavelet_3"] @@ -593,3 +598,96 @@ def test_weber_no_crash_without_crosshair(self, nc_result): def test_weber_nrm_figures_present_with_crosshair(self, nc_result_with_crosshair): figs = nc_result_with_crosshair["figures"] assert any(k.startswith("nrm_weber_") for k in figs) + + def test_original_present_with_crosshair(self, nc_result_with_crosshair): + assert "original" in nc_result_with_crosshair["figures"] + + +class TestCrosshairToCroppedPx: + """Unit tests for the helper that overlays the user's cross-section line + directly onto the Image A/B map panels — converts a normalised [0,1] + crosshair into pixel coords matching _crop_border's own offset math.""" + + def test_none_crosshair_returns_none(self): + assert SpatialDetailAnalyzer._crosshair_to_cropped_px(None, (200, 200), 0.05) is None + + def test_full_diagonal_line_cropped(self): + # 200x200 array, 5% crop -> n = 10 px removed from each edge. + crosshair = {"x0": 0.0, "y0": 0.0, "x1": 1.0, "y1": 1.0} + result = SpatialDetailAnalyzer._crosshair_to_cropped_px(crosshair, (200, 200), 0.05) + assert result == pytest.approx((-10.0, -10.0, 190.0, 190.0)) + + def test_center_point_unaffected_by_crop_offset(self): + crosshair = {"x0": 0.5, "y0": 0.5, "x1": 0.5, "y1": 0.5} + x0, y0, x1, y1 = SpatialDetailAnalyzer._crosshair_to_cropped_px( + crosshair, (200, 200), 0.05) + assert x0 == pytest.approx(90.0) + assert y0 == pytest.approx(90.0) + assert (x0, y0) == pytest.approx((x1, y1)) + + def test_small_array_below_crop_threshold_no_offset(self): + # _crop_border only crops when shape > 2n; a tiny array stays uncropped, + # so the returned pixel coords must have zero offset applied. + crosshair = {"x0": 0.0, "y0": 0.0, "x1": 1.0, "y1": 1.0} + result = SpatialDetailAnalyzer._crosshair_to_cropped_px(crosshair, (10, 10), 0.5) + assert result == pytest.approx((0.0, 0.0, 10.0, 10.0)) + + +class TestSectionSpatialReportOrder: + """Integration check on report/report_builder.py::_section_spatial's HTML + output: the reorganized subsection order (8a Background, 8b Original Image, + 8c Log-Ratio Distribution, 8d LoG, 8e Wavelet, 8f Gradient, 8g Local Std, + 8h Weber, 8i NC overview), and the fix for the alphabetic-sort bug that put + e.g. nrm_std_10px before nrm_std_3px/5px.""" + + @pytest.fixture(scope="class") + @classmethod + def section_html(cls, nc_result): + from core.models import AnalysisResult + from report.report_builder import ReportBuilder + + ra = AnalysisResult(label="A", spatial_metrics=nc_result) + rb = AnalysisResult(label="B", spatial_metrics=nc_result) + return ReportBuilder()._section_spatial(ra, rb) + + def test_headings_present_in_order(self, section_html): + import re + expected = ["8a", "8b", "8c", "8d", "8e", "8f", "8g", "8h", "8i"] + found = re.findall(r"' - "Pixel-wise log₁₀(A / B) ratio distributions. " - "Each row is one Section 8 calculation, shown as a " - "violin plot (kernel density estimate of the log-ratio pixel values, " - "randomly subsampled for display) with an IQR box-plot overlay: " + "Masked-region pixel value distributions (Image A vs Image B). " + "Each row is one Section 8j metric/scale, shown as a " + "violin plot (kernel density estimate, randomly subsampled for " + "display) with an IQR box-plot overlay: " "a cyan box spanning Q1–Q3, " "a magenta centre line at the " - "median. Nebula = pixels both " - "images classify as nebula; Background " - "= pixels both images classify as background sky — the same shared masks " - "illustrated above and used for the noise-corrected scores in 8d–8h, here " - "shown as full distributions rather than a single median ratio. " - "A red dashed line marks zero " - "(A = B, equal). Units are log₁₀ — ±0.3 ≈ a " - "2× difference, ±1.0 ≈ a 10× difference. For the " - "Original and Wavelet rows specifically (the only two families that can go " - "negative), sign was discarded before the ratio — these rows compare the " - "magnitude of structure, not signed brightness. " - "How to read it: a Background row centred at zero with a narrow IQR is " - "the expected noise floor; a Nebula row with a similarly narrow, zero-centred " - "distribution means the two filters agree at that scale. A Nebula median shifted " - "away from zero, or an IQR visibly wider than the Background row's, indicates a " - "real structural difference between the filters at that scale rather than noise. " - "Each row's x-axis is independently clipped to its own 1st–99th percentile " - "range so the IQR box stays visible; rows with rare extreme-outlier log-ratios (e.g. " - "Weber contrast near dark-sky pixels, see 8h) may have a small fraction of the " - "violin's tail extend beyond the visible axis — see the per-pixel correlation " - "scatter next to each map figure below (8d–8h) for the full, unclipped " - "upper-tail behaviour." + "median. These are the raw masked-pixel magnitude populations the " + "Mean/Std/Ratio/Significance columns above are computed from — see the " + "table for exact values (computed from the full, unsampled population, not " + "this figure's subsampled copy). " + "Each row's x-axis is independently clipped to its own 1st–99th " + "percentile range so the IQR box stays visible; rows with rare " + "extreme-outlier magnitudes (e.g. Weber contrast near dark-sky pixels) may " + "have a small fraction of the violin's tail extend beyond the visible axis." "
" ) return img_html, caption_html @@ -617,6 +622,38 @@ def _nc_ratio_rows(score_a: dict, score_b: dict, ratio: dict, scale_label, val_f return rows +def _localmax_rows(localmax: dict, rows: list, val_fmt: str = ".3f") -> str: + """BuildSection 8 measures Detail — how much real, resolvable structure ' @@ -1426,12 +1463,6 @@ def _sig(key): vb = [s[key] for s in stars_b if s.get(key) is not None] return _psf_stat_test(va, vb) # (html, p | None) - def _sig_td(html, p): - if p is None: - return f"
Ratio A/B for every noise-corrected method plotted against its approximate spatial scale. Scale units differ by method (see 8d–8h methodology boxes) — use this chart to spot which spatial-scale regime favours which filter, -not to compare absolute ratio values across methods.
""" +not to compare absolute ratio values across methods. Error bars show an +approximate relative uncertainty, propagated from the coefficient of +variation (standard deviation ÷ median) of each image's own nebula-region pixel +population — not a formal confidence interval on the median, since the nebula and +background populations behind each score are not pixel-paired between A and B. + +| Scale | {ra.label} (mean ± SD) | {rb.label} (mean ± SD) | +Ratio A/B (geo. mean) | Significance | N px / % area |
|---|
Local-maxima masked ratio A/B for every metric plotted against its +approximate spatial scale (same scale convention as 8i). Compare within a method's own +line, not numerically across methods. Error bars show ±1 standard deviation of +the per-pixel log-ratio population within each scale's own mask, converted to linear +ratio units — an exact spread measure, since the masked pixels are genuinely paired +between Image A and Image B at each scale.
+{localmax_dist_img} +{localmax_dist_caption} +{_hires_img_tag(figs.get("localmax_mask_illustration"), "Local-maxima mask grid")} +Local-maxima masks for every metric/scale row in the table above — +one panel per combination, rows grouped by metric family, columns ordered +smallest→largest kernel/scale (Wavelet has only 2 display scales, so its third +column is blank). Each panel shows the exact mask used to compute that row's +statistics, overlaid on that metric's own |A| magnitude map.
""" # ── 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 081f372..c469bfe 100644 --- a/tests/test_analysis/test_spatial_detail.py +++ b/tests/test_analysis/test_spatial_detail.py @@ -277,24 +277,24 @@ 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 + score, noise, neb_std = analyzer._nc_score(detail, None, bg_mask) + assert score is None and noise is None and neb_std 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 + score, noise, neb_std = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None and neb_std 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 + score, noise, neb_std = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None and neb_std is None def test_zero_noise_floor_returns_none(self): analyzer = SpatialDetailAnalyzer() @@ -304,8 +304,8 @@ def test_zero_noise_floor_returns_none(self): 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 + score, noise, neb_std = analyzer._nc_score(detail, mask_neb_shared, bg_mask) + assert score is None and noise is None and neb_std is None def test_valid_masks_return_ratio(self): analyzer = SpatialDetailAnalyzer() @@ -316,11 +316,50 @@ def test_valid_masks_return_ratio(self): 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) + score, noise, neb_std = analyzer._nc_score(detail, mask_neb_shared, bg_mask) assert score == pytest.approx(5.0) + assert neb_std == pytest.approx(0.0) # nebula region is a constant 10.0 block assert noise == pytest.approx(2.0) +class TestComputeNcRatioErrors: + """Direct unit tests of _compute_nc_ratio_errors's CV-propagation contract.""" + + def test_known_values_hand_computed(self): + # median_neb_a = score_a * noise_a = 5.0 * 2.0 = 10.0; cv_a = 1.0/10.0 = 0.1 + # median_neb_b = score_b * noise_b = 3.0 * 2.0 = 6.0; cv_b = 0.6/6.0 = 0.1 + ratio = {1: 5.0 / 3.0} + score_a, score_b = {1: 5.0}, {1: 3.0} + noise_a, noise_b = {1: 2.0}, {1: 2.0} + neb_std_a, neb_std_b = {1: 1.0}, {1: 0.6} + out = SpatialDetailAnalyzer._compute_nc_ratio_errors( + ratio, score_a, score_b, noise_a, noise_b, neb_std_a, neb_std_b) + expected = (5.0 / 3.0) * (0.1 ** 2 + 0.1 ** 2) ** 0.5 + assert out[1] == pytest.approx(expected, rel=1e-9) + + def test_missing_input_gives_none(self): + ratio = {1: 1.5} + score_a, score_b = {1: 5.0}, {} # score_b missing for scale 1 + noise_a, noise_b = {1: 2.0}, {1: 2.0} + neb_std_a, neb_std_b = {1: 1.0}, {1: 0.6} + out = SpatialDetailAnalyzer._compute_nc_ratio_errors( + ratio, score_a, score_b, noise_a, noise_b, neb_std_a, neb_std_b) + assert out[1] is None + + def test_zero_median_neb_guard(self): + ratio = {1: 0.0} + score_a, score_b = {1: 0.0}, {1: 3.0} # median_neb_a = 0*2.0 = 0 + noise_a, noise_b = {1: 2.0}, {1: 2.0} + neb_std_a, neb_std_b = {1: 1.0}, {1: 0.6} + out = SpatialDetailAnalyzer._compute_nc_ratio_errors( + ratio, score_a, score_b, noise_a, noise_b, neb_std_a, neb_std_b) + assert out[1] is None + + def test_empty_ratio_dict_returns_empty(self): + out = SpatialDetailAnalyzer._compute_nc_ratio_errors({}, {}, {}, {}, {}, {}, {}) + assert out == {} + + class TestLogRatioHelper: """Direct unit tests of _log_ratio_map's epsilon-floor and sign-discard contract.""" @@ -405,6 +444,161 @@ def test_mixed_sizes_keeps_only_large_object(self): assert result[24, 24] +class TestLocalMaximaMask: + """Direct unit tests of _local_maxima_mask's scale-adaptive peak detection, + region growth, and noise-suppression contract.""" + + def test_finds_single_isolated_peak(self): + data = np.zeros((30, 30), dtype=np.float32) + data[15, 15] = 100.0 + mask = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=5, prominence_percentile=90.0, region_px=0) + assert mask[15, 15] + assert mask.sum() == 1 # no region growth requested + + def test_footprint_scaling_separates_vs_merges_nearby_peaks(self): + data = np.zeros((40, 40), dtype=np.float32) + data[15, 15] = 80.0 + data[15, 19] = 100.0 # 4 px away, slightly taller + small = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=3, prominence_percentile=10.0, region_px=0) + large = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=9, prominence_percentile=10.0, region_px=0) + assert small[15, 15] and small[15, 19] # both survive with a tight footprint + assert not large[15, 15] and large[15, 19] # large footprint suppresses the weaker peak + + def test_prominence_percentile_filters_low_peaks(self): + rng = np.random.default_rng(0) + data = rng.uniform(0, 10, size=(30, 30)).astype(np.float32) + data[15, 15] = 100.0 # one dominant spike + loose = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=3, prominence_percentile=50.0, region_px=0) + strict = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=3, prominence_percentile=99.5, region_px=0) + assert loose.sum() > strict.sum() + assert strict[15, 15] # the dominant spike always survives + + def test_region_px_grows_mask_around_peak(self): + data = np.zeros((30, 30), dtype=np.float32) + data[15, 15] = 100.0 + none = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=5, prominence_percentile=90.0, region_px=0) + grown = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=5, prominence_percentile=90.0, region_px=3) + assert np.count_nonzero(grown) > np.count_nonzero(none) + assert np.all(grown[none]) # superset of the ungrown mask + + def test_presmooth_suppresses_noise_driven_detections(self): + rng = np.random.default_rng(1) + data = rng.normal(0, 5.0, size=(60, 60)).astype(np.float32) + # Broad "real" feature: a Gaussian bump, amplitude well above the noise floor. + yy, xx = np.mgrid[0:60, 0:60] + bump = 60.0 * np.exp(-(((yy - 40) ** 2 + (xx - 40) ** 2) / (2 * 3.0 ** 2))) + data += bump.astype(np.float32) + unsmoothed = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=3, prominence_percentile=90.0, region_px=0, presmooth_sigma=0.0) + smoothed = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=3, prominence_percentile=90.0, region_px=0, presmooth_sigma=2.0) + assert smoothed.sum() < unsmoothed.sum() # fewer spurious noise-driven peaks + assert smoothed[40, 40] # the genuine broad feature still survives + + def test_flat_image_returns_empty_mask(self): + data = np.full((20, 20), 5.0, dtype=np.float32) + mask = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=5, prominence_percentile=90.0, region_px=2) + assert not mask.any() + + def test_empty_array_returns_empty_mask(self): + data = np.zeros((0, 0), dtype=np.float32) + mask = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=5, prominence_percentile=90.0, region_px=2) + assert mask.shape == (0, 0) + + def test_footprint_forced_odd_and_minimum_three(self): + data = np.zeros((20, 20), dtype=np.float32) + data[10, 10] = 50.0 + for fp in (0, 2, 4): + mask = SpatialDetailAnalyzer._local_maxima_mask( + data, footprint_px=fp, prominence_percentile=90.0, region_px=0) + assert mask[10, 10] # no crash, still detects the obvious peak + + +class TestLocalMaxStats: + """Direct unit tests of _localmax_stats's masked-region aggregation contract.""" + + def test_empty_mask_returns_none_stats(self): + abs_a = np.ones((10, 10), dtype=np.float32) + abs_b = np.ones((10, 10), dtype=np.float32) + diff = np.zeros((10, 10), dtype=np.float32) + mask = np.zeros((10, 10), dtype=bool) + rng = np.random.default_rng(0) + stats = SpatialDetailAnalyzer._localmax_stats(abs_a, abs_b, diff, mask, rng) + vals_a, vals_b = stats.pop("vals_a"), stats.pop("vals_b") + assert stats == {"mean_a": None, "mean_b": None, "std_a": None, "std_b": None, + "ratio": None, "log_ratio_std": None, "p_value": None, "cliffs_delta": None, + "n_px": 0, "pct_area": 0.0} + assert vals_a.size == 0 + assert vals_b.size == 0 + + def test_known_values_give_expected_means_and_ratio(self): + abs_a = np.full((10, 10), 4.0, dtype=np.float32) + abs_b = np.full((10, 10), 2.0, dtype=np.float32) + diff = np.full((10, 10), np.log10(2.0), dtype=np.float32) # log10(4/2) + mask = np.zeros((10, 10), dtype=bool) + mask[2:5, 2:5] = True # 9 px + rng = np.random.default_rng(0) + stats = SpatialDetailAnalyzer._localmax_stats(abs_a, abs_b, diff, mask, rng) + assert stats["n_px"] == 9 + assert np.isclose(stats["mean_a"], 4.0) + assert np.isclose(stats["mean_b"], 2.0) + assert stats["std_a"] == 0.0 # constant array within the mask + assert stats["std_b"] == 0.0 + assert np.isclose(stats["ratio"], 2.0, rtol=1e-5) + assert stats["log_ratio_std"] == pytest.approx(0.0) # diff is constant within the mask + assert np.isclose(stats["pct_area"], 9.0) # 9 of 100 px + assert stats["vals_a"].size == 9 + assert stats["vals_b"].size == 9 + # A is uniformly higher than B within the mask -> fully separable. + assert stats["p_value"] is not None and stats["p_value"] < 0.05 + assert stats["cliffs_delta"] is not None and stats["cliffs_delta"] > 0.9 + + def test_log_ratio_std_reflects_varying_diff_within_mask(self): + abs_a = np.full((10, 10), 4.0, dtype=np.float32) + abs_b = np.full((10, 10), 2.0, dtype=np.float32) + diff = np.zeros((10, 10), dtype=np.float32) + mask = np.zeros((10, 10), dtype=bool) + mask[2:5, 2:5] = True # 9 px + # Known, varying log-ratio values within the mask -> hand-computable std. + diff_vals = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8], dtype=np.float32) + diff[2:5, 2:5] = diff_vals.reshape(3, 3) + rng = np.random.default_rng(0) + stats = SpatialDetailAnalyzer._localmax_stats(abs_a, abs_b, diff, mask, rng) + assert stats["log_ratio_std"] == pytest.approx(float(np.std(diff_vals)), rel=1e-5) + assert stats["log_ratio_std"] > 0.0 + + def test_mismatched_shapes_crop_to_common_size(self): + abs_a = np.full((12, 12), 4.0, dtype=np.float32) + abs_b = np.full((10, 10), 2.0, dtype=np.float32) + diff = np.full((10, 10), np.log10(2.0), dtype=np.float32) + mask = np.ones((10, 10), dtype=bool) + rng = np.random.default_rng(0) + stats = SpatialDetailAnalyzer._localmax_stats(abs_a, abs_b, diff, mask, rng) + assert stats["n_px"] == 100 + + def test_vals_capped_at_dist_max_samples(self, monkeypatch): + import analysis.image_filters as image_filters_module + monkeypatch.setattr(image_filters_module, "SECTION8_LOCALMAX_DIST_MAX_SAMPLES", 10) + abs_a = np.full((20, 20), 4.0, dtype=np.float32) + abs_b = np.full((20, 20), 2.0, dtype=np.float32) + diff = np.full((20, 20), np.log10(2.0), dtype=np.float32) + mask = np.ones((20, 20), dtype=bool) # 400 px, well above the patched cap + rng = np.random.default_rng(0) + stats = SpatialDetailAnalyzer._localmax_stats(abs_a, abs_b, diff, mask, rng) + assert stats["n_px"] == 400 # true, uncapped count + assert stats["vals_a"].size == 10 + assert stats["vals_b"].size == 10 + + class _FakeMaskImage: """Minimal duck-typed stand-in for AstroImage, exposing only what _make_masks reads (background_rms, background_subtracted()).""" @@ -556,6 +750,76 @@ def test_absent_in_single_image_mode(self, astro_image_a, key): assert key not in result["figures"] +class TestLocalMaxIntegration: + """Section 8j: result["localmax"] population across all metric/scale + combinations, mirroring the corr_* key set in TestCorrelationScatterFigures.""" + + _LOCALMAX_KEYS = ( + [f"std_{ks}px" for ks in STD_KERNEL_SIZES] + + [f"log_{s}" for s in LOG_SIGMAS] + + [f"gradient_{s}" for s in LOG_SIGMAS] + + ["wavelet_2", "wavelet_3"] + + [f"weber_{ks}px" for ks in WEBER_KERNEL_SIZES] + ) + + @pytest.mark.parametrize("key", _LOCALMAX_KEYS) + def test_entry_present_with_valid_stats(self, nc_result, key): + entry = nc_result["localmax"][key] + assert entry["n_px"] >= 0 + assert 0.0 <= entry["pct_area"] <= 100.0 + assert entry["ratio"] is None or entry["ratio"] > 0 + assert entry["std_a"] is None or entry["std_a"] >= 0 + assert entry["std_b"] is None or entry["std_b"] >= 0 + assert entry["p_value"] is None or 0.0 <= entry["p_value"] <= 1.0 + assert entry["cliffs_delta"] is None or -1.0 <= entry["cliffs_delta"] <= 1.0 + assert entry["vals_a"].size <= entry["n_px"] + assert entry["vals_b"].size <= entry["n_px"] + + def test_empty_in_single_image_mode(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert result["localmax"] == {} + + +class TestLocalMaxFigures: + def test_ratio_overview_present_in_two_image_mode(self, nc_result): + assert "localmax_ratio_overview" in nc_result["figures"] + + def test_ratio_overview_absent_in_single_image_mode(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert "localmax_ratio_overview" not in result["figures"] + + def test_mask_illustration_present_in_two_image_mode(self, nc_result): + assert "localmax_mask_illustration" in nc_result["figures"] + assert isinstance(nc_result["figures"]["localmax_mask_illustration"], str) + assert len(nc_result["figures"]["localmax_mask_illustration"]) > 0 + + def test_mask_illustration_absent_in_single_image_mode(self, astro_image_a): + result = SpatialDetailAnalyzer().analyze(astro_image_a) + assert "localmax_mask_illustration" not in result["figures"] + + def test_mask_illustration_is_a_multi_row_grid(self, nc_result): + # Regression guard: the old single-scale illustration was one ~9-inch-tall + # panel (~1350px at dpi=150); the grid is 5 rows tall (~2700px) -- decode + # the PNG and check its height reflects the multi-row layout, not the old + # single-panel figure. + import base64 + import io + from PIL import Image + png_bytes = base64.b64decode(nc_result["figures"]["localmax_mask_illustration"]) + img = Image.open(io.BytesIO(png_bytes)) + assert img.height > 2000 + + def test_localmax_entries_carry_vals_for_distribution_figure(self, nc_result): + # result["figures"] never stores the A/B distribution violin plot itself + # (that's rendered in report_builder.py from result["localmax"][key]["vals_a"/"vals_b"]), + # so verify those raw-value arrays are actually present in two-image mode. + any_present = any( + entry.get("vals_a") is not None and entry["vals_a"].size > 0 + for entry in nc_result["localmax"].values() + ) + assert any_present + + @pytest.fixture(scope="module") def nc_result_with_crosshair(nc_image_pair) -> dict: img_a, img_b = nc_image_pair @@ -636,9 +900,9 @@ def test_small_array_below_crop_threshold_no_offset(self): class TestSectionSpatialReportOrder: """Integration check on report/report_builder.py::_section_spatial's HTML output: the reorganized subsection order (8a Background, 8b Original Image, - 8c Log-Ratio Distribution, 8d LoG, 8e Wavelet, 8f Gradient, 8g Local Std, - 8h Weber, 8i NC overview), and the fix for the alphabetic-sort bug that put - e.g. nrm_std_10px before nrm_std_3px/5px.""" + 8c Mask Overview, 8d LoG, 8e Wavelet, 8f Gradient, 8g Local Std, + 8h Weber, 8i NC overview, 8j Local-Maxima Masked Metrics), and the fix for + the alphabetic-sort bug that put e.g. nrm_std_10px before nrm_std_3px/5px.""" @pytest.fixture(scope="class") @classmethod @@ -652,8 +916,8 @@ def section_html(cls, nc_result): def test_headings_present_in_order(self, section_html): import re - expected = ["8a", "8b", "8c", "8d", "8e", "8f", "8g", "8h", "8i"] - found = re.findall(r"' + "Masked-region log ratio distributions (log₁₀(A / B)). " + "Each row is one Section 8j metric/scale, shown as a " + "violin plot (kernel density estimate, randomly subsampled for " + "display) of the per-pixel log10(|A|/|B|) population within that row's " + "local-maxima mask, with an IQR box-plot overlay: " + "a cyan box spanning Q1–Q3, " + "a magenta centre line at the " + "median. The dashed vertical line marks 0 (A = B). This is the same " + "masked pixel population the “log ratio A/B (geo. mean ± SD)” " + "table column and the cross-method overview plot's error bars are computed " + "from — a roughly symmetric, unimodal shape here supports summarising " + "it with a mean ± SD. " + "Each row's x-axis is independently clipped to its own 1st–99th " + "percentile range so the IQR box stays visible." + "
" + ) + return img_html, caption_html + + def _epsf_stars_cell(psf_metrics: dict) -> str: """Format the number of stars used to build the ePSF.""" n = psf_metrics.get("epsf_n_stars") @@ -624,12 +711,15 @@ def _nc_ratio_rows(score_a: dict, score_b: dict, ratio: dict, scale_label, val_f def _localmax_rows(localmax: dict, rows: list, val_fmt: str = ".3f") -> str: """BuildA pixel-level difference map (A − B) is computed and displayed with the RdBu_r @@ -2935,6 +3024,18 @@ def panel(arr, title, caption=""): peak-to-valley swing preserved by each PSF.
""" if has_b else "") + conv_a_title = f"Convolved — {sim['label_a']}" + conv_b_title = f"Convolved — {sim['label_b']}" + conv_ref_title = f"Convolved — {sim['label_ref']}" + _testchart_images_html = ( + panel(sim['original'], 'Original test chart') + + panel(sim['conv_a'], conv_a_title) + + (panel(sim['conv_b'], conv_b_title) if has_b else '') + + (panel(sim['conv_ref'], conv_ref_title) if sim.get('conv_ref') is not None else '') + + diff_panel + ) + _testchart_images_box = _info_box(_testchart_images_html, title="Show test chart validation images", open=False) + return f"""@@ -2969,11 +3070,8 @@ def panel(arr, title, caption=""):
{diff_para}Each image is rendered at 1 image-pixel : 1 screen-pixel.
-{panel(sim['original'], 'Original test chart')} -{panel(sim['conv_a'], f"Convolved — {sim['label_a']}")} -{panel(sim['conv_b'], f"Convolved — {sim['label_b']}") if has_b else ''} -{panel(sim['conv_ref'], f"Convolved — {sim['label_ref']}") if sim.get('conv_ref') is not None else ''} -{diff_panel}{xs_block}""" +{_testchart_images_box} +{xs_block}""" # ── Section 4: Halo ─────────────────────────────────────────────────────── @@ -4162,10 +4260,12 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: localmax_rows_html = _localmax_rows(sm.get("localmax", {}), _SPATIAL_CORR_ROWS) localmax_dist_img, localmax_dist_caption = _localmax_distributions_figure(sm.get("localmax", {})) + localmax_log_ratio_dist_img, localmax_log_ratio_dist_caption = _localmax_log_ratio_distribution_figure(sm.get("localmax", {})) lm_footprint_mult = sm.get("localmax_footprint_mult", SECTION8_LOCALMAX_FOOTPRINT_MULT) lm_prominence_pctl = sm.get("localmax_prominence_percentile", SECTION8_LOCALMAX_PROMINENCE_PERCENTILE) lm_region_fraction = sm.get("localmax_region_fraction", SECTION8_LOCALMAX_REGION_FRACTION) lm_presmooth_fraction = sm.get("localmax_presmooth_fraction", SECTION8_LOCALMAX_PRESMOOTH_FRACTION) + lm_top_percent = sm.get("localmax_top_percent", SECTION8_LOCALMAX_TOP_PERCENT) localmax_methodology_box = _info_box( 'For each metric/scale combination in 8d–8h, a local-maxima mask is ' 'built independently: the combined |A|,|B| peak-source array is lightly ' @@ -4174,14 +4274,21 @@ def _section_spatial(self, ra: AnalysisResult, rb: AnalysisResult) -> str: f'their own {lm_footprint_mult:g}×(scale) neighbourhood AND exceed the ' f'{lm_prominence_pctl:g}th percentile of the smoothed array are kept as peaks and grown ' f'by {lm_region_fraction:g}× that same neighbourhood size, so the mask covers the local ' - 'region of pixels around each peak rather than a single pixel. Unlike the Nebula/Background ' - 'masks in 8c, this mask is not restricted to the nebula region — the ' - 'prominence threshold already isolates prominent features on its own, and intersecting with ' - 'the coarser nebula mask could exclude legitimate sharp features (stars, edges) outside it. ' + 'region of pixels around each peak rather than a single pixel. This peak mask is then ' + f'unioned (OR) with a top-{lm_top_percent:g}% brightness mask — pixels in the ' + 'top percentile of Image A\'s or Image B\'s own value distribution — so broad bright plateaus ' + 'that never register as a sharp local maximum are still captured, not just isolated peaks. ' + 'Unlike the Nebula/Background masks in 8c, this mask is not restricted to the ' + 'nebula region — the prominence threshold already isolates prominent features on its own, and ' + 'intersecting with the coarser nebula mask could exclude legitimate sharp features (stars, ' + 'edges) outside it. ' 'Mean A / Mean B are the average metric magnitude ± standard deviation ' - 'over the masked pixels in each image; Ratio (A/B) is the geometric mean of ' - 'the per-pixel A/B ratio at those pixels (10mean(log-ratio)), a plain ' - '×-factor. Significance is a Mann-Whitney U test (two-sided) with ' + 'over the masked pixels in each image. log ratio A/B (geo. mean ± SD) ' + 'is sampled directly from the masked pixels\' per-pixel log10(|A|/|B|) population — not ' + 'derived from Mean A/Mean B — and reported in log10 units: 0 means A = B, positive means A is ' + 'brighter, negative means B is brighter, and the equivalent linear ×-factor is ' + '10value. It is shaded a neutral blue rather than red/green, since it is not an ' + 'A-vs-B comparison. Significance is a Mann-Whitney U test (two-sided) with ' 'Cliff\'s delta effect size, comparing the full masked-pixel populations of A vs B for that ' 'row — the same test and star-rating legend (★★★ large, ★★ ' 'medium, ★ small, n.s. not significant; blue cell = p<0.05) already used for the ' @@ -4336,6 +4443,87 @@ def _family_nrm_figs(rows) -> str: if orig_fig else "" ) + # Section 8d-8h figure-heavy content, collapsed by default (closed|LoG| maps at σ = 1.5, 3, and 6 px (shared colour scale per figure): ' + 'Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a ' + 'cross-section line is set — its profile (middle-right), plus a bottom-row histogram of ' + 'the log-ratio map\'s pixel values (same colour scale). A filter preserving more fine ' + 'detail shows brighter, more defined boundaries at small σ. Each map is immediately ' + 'followed by its per-pixel correlation scatter (see methodology above).
' + 'Noise-normalised (× noise floor) — shared colour scale is a fair ' + 'A/B comparison.
' + + _family_nrm_figs(_log_rows) + ) + _log_images_box = _info_box(_log_images_html, title="Show LoG maps & figures", open=False) + + _wavelet_images_html = ( + _hires_img_tag(figs.get("wavelet_snr"), "Wavelet SNR") + + 'Per-level SNR for both filters. Level 1 SNR < 1 is expected ' + '(noise-dominated). A filter preserving more fine detail shows higher SNR at level 2.
' + + _family_figs_with_corr(_wavelet_rows, lambda k: "wavelet_level" + k.split("_", 1)[1]) + + 'Reconstructed detail images at levels 2 and 3 (shared colour scale, ' + 'diverging colourmap): Image A (top-left), Image B (top-right), log-ratio panel ' + '(middle-left) showing where fine structure differs between the two filters (sign ' + 'discarded — |A|/|B| — since wavelet reconstructions can be negative), and — when a ' + 'cross-section line is set — its profile (middle-right), plus a bottom-row histogram ' + 'of the log-ratio map\'s pixel values (same colour scale). Each map is immediately ' + 'followed by its per-pixel correlation scatter (see methodology above).
' + 'Noise-normalised (× noise floor) — shared colour scale is a fair ' + 'A/B comparison.
' + + _family_nrm_figs(_wavelet_rows) + ) + _wavelet_images_box = _info_box(_wavelet_images_html, title="Show Wavelet maps & figures", open=False) + + _gradient_images_html = ( + _family_figs_with_corr(_gradient_rows, lambda k: k) + + 'Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale ' + 'per figure): Image A (top-left), Image B (top-right), log-ratio map (middle-left), and ' + '— when a cross-section line is set — its profile (middle-right), plus a bottom-row ' + 'histogram of the log-ratio map\'s pixel values (same colour scale). A filter preserving ' + 'sharper boundaries shows brighter, more defined gradient response. Each map is ' + 'immediately followed by its per-pixel correlation scatter (see methodology above).
' + 'Noise-normalised (× noise floor) — shared colour scale is a fair ' + 'A/B comparison.
' + + _family_nrm_figs(_gradient_rows) + ) + _gradient_images_box = _info_box(_gradient_images_html, title="Show Gradient maps & figures", open=False) + + _std_images_html = ( + _family_figs_with_corr(_std_rows, lambda k: k) + + 'Local σ maps at each kernel size (shared colour scale): Image A ' + '(top-left), Image B (top-right), log-ratio map (middle-left) highlighting where one ' + 'filter preserves more local variation, and — when a cross-section line is set — its ' + 'profile (middle-right), plus a bottom-row histogram of the log-ratio map\'s pixel ' + 'values (same colour scale). Each map is immediately followed by its per-pixel ' + 'correlation scatter (see methodology above).
' + 'Noise-normalised (× noise floor) — shared colour scale is a fair ' + 'A/B comparison.
' + + _family_nrm_figs(_std_rows) + ) + _std_images_box = _info_box(_std_images_html, title="Show Local σ maps & figures", open=False) + + _weber_images_html = ( + _family_figs_with_corr(_weber_rows, lambda k: k) + + 'Per-pixel Weber fraction contrast maps (c = ΔL / L, square-root ' + 'colour scale, viridis): Image A (top-left), Image B (top-right), log-ratio panel ' + '(middle-left) showing where one image achieves greater relative contrast, and — when ' + 'a cross-section line is set — its profile (middle-right), plus a bottom-row histogram ' + 'of the log-ratio map\'s pixel values (same colour scale). Brighter regions have higher ' + 'Weber contrast — the local intensity range is large relative to the local median ' + 'luminance. High values over dark-sky regions are expected; use a nebula ROI for ' + 'meaningful filter comparison. Each map is immediately followed by its per-pixel ' + 'correlation scatter (see methodology above).
' + 'Noise-normalised (× noise floor) — shared colour scale is a fair ' + 'A/B comparison.
' + + _family_nrm_figs(_weber_rows) + ) + _weber_images_box = _info_box(_weber_images_html, title="Show Weber contrast maps & figures", open=False) + return f"""|LoG| maps at σ = 1.5, 3, and 6 px (shared colour scale per figure): -Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a -cross-section line is set — its profile (middle-right), plus a bottom-row histogram of -the log-ratio map's pixel values (same colour scale). A filter preserving more fine -detail shows brighter, more defined boundaries at small σ. Each map is immediately -followed by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{_family_nrm_figs(_log_rows)} +{_log_images_box}Per-level SNR for both filters. Level 1 SNR < 1 is expected -(noise-dominated). A filter preserving more fine detail shows higher SNR at level 2.
-| Wavelet level | {ra.label} SNR | {rb.label} SNR | |
|---|---|---|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
Reconstructed detail images at levels 2 and 3 (shared colour scale, -diverging colourmap): Image A (top-left), Image B (top-right), log-ratio panel (middle-left) -showing where fine structure differs between the two filters (sign discarded — |A|/|B| — -since wavelet reconstructions can be negative), and — when a cross-section line is set — -its profile (middle-right), plus a bottom-row histogram of the log-ratio map's pixel values -(same colour scale). Each map is immediately followed by its per-pixel correlation -scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{_family_nrm_figs(_wavelet_rows)} +{_wavelet_images_box}Gradient magnitude maps at σ = 1.5, 3, and 6 px (shared colour scale per -figure): Image A (top-left), Image B (top-right), log-ratio map (middle-left), and — when a -cross-section line is set — its profile (middle-right), plus a bottom-row histogram of the -log-ratio map's pixel values (same colour scale). A filter preserving sharper -boundaries shows brighter, more defined gradient response. Each map is immediately followed -by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{_family_nrm_figs(_gradient_rows)} +{_gradient_images_box}Local σ maps at each kernel size (shared colour scale): Image A (top-left), -Image B (top-right), log-ratio map (middle-left) highlighting where one filter preserves more -local variation, and — when a cross-section line is set — its profile (middle-right), plus a -bottom-row histogram of the log-ratio map's pixel values (same colour scale). Each map -is immediately followed by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{_family_nrm_figs(_std_rows)} +{_std_images_box}Formula: c = ΔL / L, ' @@ -4503,17 +4654,7 @@ def _family_nrm_figs(rows) -> str:
Per-pixel Weber fraction contrast maps (c = ΔL / L, square-root colour -scale, viridis): Image A (top-left), Image B (top-right), log-ratio panel (middle-left) -showing where one image achieves greater relative contrast, and — when a cross-section -line is set — its profile (middle-right), plus a bottom-row histogram of the log-ratio -map's pixel values (same colour scale). Brighter regions have higher Weber contrast — -the local intensity range is large relative to the local median luminance. High values -over dark-sky regions are expected; use a nebula ROI for meaningful filter comparison. -Each map is immediately followed by its per-pixel correlation scatter (see methodology above).
-Noise-normalised (× noise floor) — shared colour scale is a fair A/B comparison.
-{_family_nrm_figs(_weber_rows)} +{_weber_images_box}| Scale | {ra.label} (mean ± SD) | {rb.label} (mean ± SD) | -Ratio A/B (geo. mean) | Significance | N px / % area | log ratio A/B (geo. mean ± SD) | Significance | N px / % area | {localmax_rows_html}
|---|
Local-maxima masked ratio A/B for every metric plotted against its -approximate spatial scale (same scale convention as 8i). Compare within a method's own -line, not numerically across methods. Error bars show ±1 standard deviation of -the per-pixel log-ratio population within each scale's own mask, converted to linear -ratio units — an exact spread measure, since the masked pixels are genuinely paired -between Image A and Image B at each scale.
+Local-maxima masked log₁₀(A/B) for every metric plotted +against its approximate spatial scale (same scale convention as 8i). Compare within a +method's own line, not numerically across methods. Error bars show ±1 standard +deviation of the per-pixel log10(A/B) population within each scale's own mask, plotted +directly with no unit conversion — an exact spread measure, since the masked pixels are +genuinely paired between Image A and Image B at each scale.
{localmax_dist_img} {localmax_dist_caption} +{localmax_log_ratio_dist_img} +{localmax_log_ratio_dist_caption} {_hires_img_tag(figs.get("localmax_mask_illustration"), "Local-maxima mask grid")}Local-maxima masks for every metric/scale row in the table above — one panel per combination, rows grouped by metric family, columns ordered smallest→largest kernel/scale (Wavelet has only 2 display scales, so its third column is blank). Each panel shows the exact mask used to compute that row's -statistics, overlaid on that metric's own |A| magnitude map.
""" +statistics (isolated peaks unioned with the top-brightness mask — see the methodology +box above), overlaid on that metric's own |A| magnitude map.""" # ── Section 9: Signal-to-Noise Ratio ───────────────────────────────────── @@ -5043,13 +5187,6 @@ def row_pm(metric, val_a, val_b, spread_a, spread_b, fmt=".3f", '⚠ = interpret with bandwidth ' 'context (filters had different bandwidths)') - retention_block = getattr(self, "_cached_retention_html", "") - retention_section = ( - "| naming the metric + family (LoG/Wavelet/Gradient/etc.) -- used by the Section 8j combined + cross-method NC table.""" rows = "" + method_td = f" | {method_label} | " if method_label is not None else "" 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"||
| {scale_label(scale)} | " + rows += (f"|||
| {scale_label(scale)} | " f"{_val(va, val_fmt)} | " f"{_val(vb, val_fmt)} | " f"{_val(vr, val_fmt)} |
| {ks} px | " - f"{_val(va)} | " - f"{_val(vb)} | |
| {ks} px | " - f"{_val(va, '.4f')} | " - f"{_val(vb, '.4f')} | |
| Level {lvl} (~{scale_approx}px scale) | " - f"{_val(va)} | " - f"{_val(vb)} |
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
| Wavelet level | {ra.label} SNR | {rb.label} SNR |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
| Kernel size | {ra.label} | {rb.label} |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
| Kernel size | {ra.label} (99th pct c) | {rb.label} (99th pct c) |
|---|
| Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
| Method | Scale | {ra.label} (NC score) | {rb.label} (NC score) | Ratio A/B |
|---|
Local-maxima masked log₁₀(A/B) for every metric plotted
against its approximate spatial scale (same scale convention as 8i). Compare within a
diff --git a/tests/test_analysis/test_edge_analyzer.py b/tests/test_analysis/test_edge_analyzer.py
index bf3636a..2d2ff38 100644
--- a/tests/test_analysis/test_edge_analyzer.py
+++ b/tests/test_analysis/test_edge_analyzer.py
@@ -65,14 +65,14 @@ 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)
+ _, 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)
+ _, esf, _ = ea._extract_esf(roi, edge_info)
assert ea._esf_quality(esf) < EDGE_ESF_MIN_MONOTONICITY
def test_perfectly_flat_scores_zero(self):
@@ -93,7 +93,7 @@ def test_clean_edge_scores_high_at_various_angles(self, angle_deg):
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)
+ _, esf, _ = ea._extract_esf(roi, edge_info)
assert ea._esf_quality(esf) > 0.8
@@ -104,7 +104,7 @@ def test_no_nan_in_returned_esf(self):
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)
+ 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))
@@ -113,14 +113,14 @@ 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)
+ 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)
+ _, esf, _ = ea._extract_esf(roi, edge_info)
assert esf.min() >= -1e-9
assert esf.max() <= 1.0 + 1e-9
From aab18da6983ba4401827ffaf1db11d20dba59ab7 Mon Sep 17 00:00:00 2001
From: Brent <52629076+brentmantooth@users.noreply.github.com>
Date: Sat, 18 Jul 2026 15:10:47 -0400
Subject: [PATCH 09/13] claude setting update
---
.claude/settings.json | 14 +++++++++++++-
1 file changed, 13 insertions(+), 1 deletion(-)
diff --git a/.claude/settings.json b/.claude/settings.json
index ff6952b..e155565 100644
--- a/.claude/settings.json
+++ b/.claude/settings.json
@@ -7,7 +7,19 @@
"WebSearch",
"PowerShell(conda activate astrolab)",
"PowerShell(python -m pytest tests/ -m \"not slow\" -q 2>&1)",
- "PowerShell(Get-ChildItem -Path \"$env:USERPROFILE\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue -Recurse -Depth 2 | Select-Object -First 5 FullName; Get-ChildItem -Path \"C:\\\\ProgramData\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue | Select-Object FullName)"
+ "PowerShell(Get-ChildItem -Path \"$env:USERPROFILE\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue -Recurse -Depth 2 | Select-Object -First 5 FullName; Get-ChildItem -Path \"C:\\\\ProgramData\" -Filter \"condabin\" -Directory -ErrorAction SilentlyContinue | Select-Object FullName)",
+ "PowerShell(Start-Process -FilePath \"python\" -ArgumentList \"AstroImageLab.py\" -WorkingDirectory \"d:\\\\GitHub\\\\AstroImageLab\" -PassThru)",
+ "PowerShell(Stop-Process -Id 25432 -Force -ErrorAction SilentlyContinue)",
+ "Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/measure_control_panel.py\")",
+ "Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/render_control_panel.py\")",
+ "Read(//c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/**)",
+ "Bash(python -c ' *)",
+ "Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/crop_image.py\")",
+ "Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/render_main_window.py\")",
+ "Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_roi_reset.py\")",
+ "Bash(PYTHONIOENCODING=utf-8 \"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_roi_reset.py\")",
+ "Bash(PYTHONIOENCODING=utf-8 \"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_toolbar_sync.py\")",
+ "Bash(python -m py_compile gui/control_panel.py gui/main_window.py gui/image_panel.py)"
]
}
}
From 7c915f55d86814976273ffce6610a8c1c23549e7 Mon Sep 17 00:00:00 2001
From: Brent <52629076+brentmantooth@users.noreply.github.com>
Date: Sat, 18 Jul 2026 21:46:10 -0400
Subject: [PATCH 10/13] Replace Weber contrast with Local Entropy, fix mathtext
race, fix Section 6 edge width bug
Section 8 (Spatial Detail): remove the Weber fraction contrast metric and
add a Local Entropy family (map, contrast ratio, noise-corrected score)
reusing the existing shared NC/contrast helpers. Reuses Weber's former 8h
letter slot so no cross-reference renumbering is needed.
Fix a matplotlib mathtext ParseException race condition: pyparsing's
packrat cache (enabled globally by matplotlib's mathtext grammar) is not
thread-safe, and this app renders figures from multiple analyzer threads
concurrently. Serialize all savefig() calls through a lock in
core/fig_utils.py and remove report_builder.py's duplicate, unprotected
copy of fig_to_b64.
Vectorize the new local entropy map (per-bin uniform_filter box sums
instead of a per-pixel generic_filter callback) after it was found to
regress the test suite runtime ~2x.
Section 6 (Edge Analysis): fix two bugs in EdgeAnalyzer._extract_esf.
The "Scan start" marker was drawn through the Sobel-detected gradient
peak instead of the ROI's own array center, which is the actual pivot
rotate() uses -- the marker could land off the scan-direction line.
Separately, the rotation angle formula aligned edges horizontally
instead of vertically, so the ESF's column-wise averaging integrated
across the transition instead of along it, inflating measured edge
widths by 7-14x since the feature's introduction. Add a ground-truth
width-accuracy test (checking measured width against the analytically
known erf-profile width of a Gaussian-blurred synthetic edge), since the
existing tests only checked monotonicity/shape and never caught this.
Co-Authored-By: Claude Sonnet 5 Local entropy is best read as a secondary texture-complexity signal, '
+ 'not a replacement for σ/LoG/gradient/wavelet — it answers a different question '
+ '("how unpredictable is the local tonal distribution?") rather than "how much does '
+ 'brightness vary?" or "how sharp are the edges?". It is most informative when two images '
+ 'score similarly on σ or wavelet power but still look visibly different in local texture '
+ '(fine nebulosity, dust lanes, mottled galaxy structure) — entropy can surface that '
+ 'difference where amplitude- and edge-based metrics do not. Because entropy alone cannot '
+ 'distinguish real tonal richness from shot noise, always read it alongside its '
+ 'noise-corrected score (8h, 8i), not in isolation. Per-pixel Weber fraction contrast maps (c = ΔL / L, square-root '
- 'colour scale, viridis): Image A (top-left), Image B (top-right), log-ratio panel '
- '(middle-left) showing where one image achieves greater relative contrast, and — when '
- 'a cross-section line is set — its profile (middle-right), plus a bottom-row histogram '
- 'of the log-ratio map\'s pixel values (same colour scale). Brighter regions have higher '
- 'Weber contrast — the local intensity range is large relative to the local median '
- 'luminance. High values over dark-sky regions are expected; use a nebula ROI for '
- 'meaningful filter comparison. Each map is immediately followed by its per-pixel '
- 'correlation scatter (see methodology above). Local Shannon entropy maps (bits, log₂, viridis colour '
+ 'scale): Image A (top-left), Image B (top-right), log-ratio panel (middle-left) '
+ 'showing where one image has richer local tonal structure, and — when a '
+ 'cross-section line is set — its profile (middle-right), plus a bottom-row '
+ 'histogram of the log-ratio map\'s pixel values (same colour scale). Brighter '
+ 'regions have higher local entropy — a richer, less predictable gray-level '
+ 'distribution within the window (nebulosity, mottled dust, unresolved stars — or '
+ 'noise). Each map is immediately followed by its per-pixel correlation scatter '
+ '(see methodology above). Noise-normalised (× noise floor) — shared colour scale is a fair '
'A/B comparison. 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. Formula: H = −Σ pi log₂(pi) '
+ '(Shannon entropy, in bits), computed from the local gray-level histogram over '
+ f'{SECTION8_ENTROPY_N_BINS} bins within the K × K kernel. Output is bounded '
+ f'0 ≤ H ≤ log₂({SECTION8_ENTROPY_N_BINS}) ≈ {math.log2(SECTION8_ENTROPY_N_BINS):.1f} '
+ 'bits; 0 means every pixel in the window is identical, the maximum means every gray '
+ 'level is equally represented. Why deliberate binning: computing entropy directly on continuous '
+ 'float32 data would return a value close to the maximum almost everywhere, since every '
+ 'pixel value in a small window is nearly unique — that measures floating-point '
+ 'precision, not real tonal diversity. The map is instead computed on each image\'s own '
+ f'data, independently quantized into {SECTION8_ENTROPY_N_BINS} levels from its own '
+ f'{SECTION8_ENTROPY_CLIP_PERCENTILE:g}–{100 - SECTION8_ENTROPY_CLIP_PERCENTILE:g} '
+ 'percentile range (not a joint A/B range, and not the display stretch — this is '
+ 'done on linear, mean-signal-normalised data). Interpretation: high entropy means rich, unpredictable local tonal '
+ 'structure — tangled nebulosity, mottled dust, unresolved star fields — '
+ 'but also shot noise; low entropy means a smooth or flat local region '
+ '(uniform sky, a saturated core, a smooth gradient). Entropy alone cannot distinguish '
+ 'real structure from noise — that is what the noise-corrected score (table below, '
+ 'and 8i) is for. How this differs from σ / LoG / gradient / wavelet: those '
+ 'metrics respond to amplitude or edge sharpness; entropy responds to distributional '
+ '"unpredictability" instead, and ignores spatial arrangement entirely — a smooth '
+ 'gradient, scattered noise, and a small filament network can all share a similar local '
+ 'histogram, and therefore similar entropy, even though they look nothing alike. It is '
+ 'most useful as a secondary signal when two images score similarly on σ or wavelet '
+ 'power but still show visibly different internal tonal texture (see the glossary '
+ 'comparison table above). Kernel sizes: 5, 9, and 17 px — deliberately larger than Local '
+ 'σ\'s (8g), since a small window gives a statistically unstable histogram estimate '
+ '(a 3 × 3 window has only 9 samples to build a histogram from). Entropy contrast ratio (table below): median(nebula entropy) / '
+ 'median(background entropy) — the same formula as Local σ\'s contrast ratio '
+ '(8g), applied to the entropy map instead. f9roZ~4)M*$=8m1gyN(1;
z9!mFT?^Sh&wW5mkxV&@jokl3T7?7Z80r>iN;!ze$*#oNps6x-}9O>W3tmAI)ZS&Hf
zZ&;(QIX^dxpzT~!-`5<<4hqft`n`WPutF$=a_(40rocsy?KOV(XMeW++|T{oFPQbW
zA@aH)0z-0D*EWp?Ltto794*E+Oo2R$0yJ oX5-hdx}y
z{ac@ppZJ4+EPnX;d-3qehcZ(U*u*^En%#` wrapper for any Mann-Whitney significance column — used by both Section 4's PSF table and Section 8j's local-maxima table |
| `SpatialDetailAnalyzer._ratio_series_with_errors(ratios_by_method, errors_by_method=None)` | `analysis/image_filters.py` | `{method: {scale: value}}` (+ optional matching errors) → sorted `{method: [(x_px, value, error_or_None), ...]}` point lists for a cross-method overview line plot; shared by `_plot_nc_ratio_overview` (8i) and `_plot_localmax_ratio_overview` (8j) |
-| `_nc_ratio_rows(score_a, score_b, ratio, scale_label, val_fmt=".3f", method_label=None)` | `report_builder.py` | ` ` rows for a noise-corrected score table (Scale \| A \| B \| Ratio A/B). Pass `method_label` to prepend a Method-name ` ` when consolidating several per-family tables that share this schema into one combined table — precedent: Section 8j's cross-method NC table, which concatenates the row-strings from all five `_nc_ratio_rows` calls (LoG/Wavelet/Gradient/Std/Weber) into a single ` ` |
+| `_nc_ratio_rows(score_a, score_b, ratio, scale_label, val_fmt=".3f", method_label=None)` | `report_builder.py` | `
` rows for a noise-corrected score table (Scale \| A \| B \| Ratio A/B). Pass `method_label` to prepend a Method-name ` ` when consolidating several per-family tables that share this schema into one combined table — precedent: Section 8j's cross-method NC table, which concatenates the row-strings from all five `_nc_ratio_rows` calls (LoG/Wavelet/Gradient/Std/Entropy) into a single ` ` |
---
@@ -489,7 +489,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` / `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. |
+| `contrast_ratios_b` / `entropy_contrast_ratio_b` always present | `SpatialDetailAnalyzer.analyze()` always includes `contrast_ratios_b: {}` and `entropy_contrast_ratio_b: {}` even in single-image mode. Neither is ever `None` or absent — check `not b_ratios` / `not ecr_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. |
---
@@ -518,7 +518,7 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html
| Mixed float32/float64 arithmetic silently widens to float64 | NumPy upcasts when operands differ (e.g. `float32_array - float64_scalar`). If photutils ever returns a float64 background model, `background_subtracted()` will silently return float64. Guard by adding `.astype(np.float32)` at the end of `background_subtracted()` in `astro_image.py` if this is observed. |
| `_section_snr` crashed when SNR metric is unchecked | `_plot_snr_pair` (`report_builder.py`) built a `panels` list filtered to non-`None` entries but never checked whether it was empty before calling `plt.subplots(1, len(panels), ...)` — 0 columns raised `ValueError: Number of columns must be a positive integer, not 0`. Hit whenever SNR is unchecked while another metric (e.g. Power Spectrum) is run. Fixed with an early `if not panels: return None` guard, matching `_plot_radial_overlay`/`_plot_radial_ratio_db`; both call sites already pipe the result through `_img_tag`, which turns `None` into `""`. |
| New Section 8 panel key doesn't need Report Inspector code changes | `gui/report_inspector.py` is fully generic — driven entirely by a companion `
',
title="Which concept each metric primarily measures")
+ + _info_box(
+ '` tags (e.g. "see 8i for…", "(8d–8h, 8i)"). After adding, removing, or renumbering a subsection, `grep` the function (and `_SPATIAL_GLOSSARY_HTML`) for every old *and* new heading letter — HTML renders a stale cross-reference without error, it just silently misdirects the reader to the wrong subsection. A self-reference inside a family's own info box (e.g. Gradient's "same framework as the other N families") must enumerate the other letters explicitly rather than use a dash range — under a non-contiguous lettering (Gradient kept letter `8f` while Std/Weber moved past it into the Contrast group), a `8d–8h` range would wrongly include Gradient's own letter. |
+| Renumbering a Section 8 subsection misses caption cross-references | Section 8's sub-heading letters (currently 8a–8j: 8a Background/Key Terms, 8b Original Image, 8c Mask Overview, 8d–8f detail-based families [LoG, Wavelet, Gradient], 8g–8h contrast/texture-based families [Local σ, Local Entropy], 8i Noise-Corrected Cross-Method Overview, 8j Local-Maxima Masked Metrics) are referenced by literal string in caption/info-box text scattered throughout `_section_spatial` — not just in the `
` tags (e.g. "see 8i for…", "(8d–8h, 8i)"). After adding, removing, or renumbering a subsection, `grep` the function (and `_SPATIAL_GLOSSARY_HTML`) for every old *and* new heading letter — HTML renders a stale cross-reference without error, it just silently misdirects the reader to the wrong subsection. A self-reference inside a family's own info box (e.g. Gradient's "same framework as the other N families") must enumerate the other letters explicitly rather than use a dash range — under a non-contiguous lettering (Gradient kept letter `8f` while Std/Entropy moved past it into the Contrast group), a `8d–8h` range would wrongly include Gradient's own letter. Precedent: when Weber contrast (formerly 8h) was replaced by Local Entropy, keeping the same letter for the new family avoided a renumbering pass entirely — every existing `8h`/`8d–8h` cross-reference stayed numerically valid, only the prose describing 8h's content changed. |
| Stale ROI crashes Section 8 with "index -1 is out of bounds for axis 0 with size 0" | `MainWindow._on_image_loaded()` (`gui/main_window.py`) now resets `self._roi`/`self._crosshair` to `None` (plus both panels' visual overlays, via `ImagePanel.clear_roi_overlay()`/`clear_line_overlay()`) on every new main-image load — the choke point is `ImagePanel.image_loaded`, emitted only from `_open_file`/`load_path`, never from `set_starless_path`, so attaching a starless companion correctly does *not* wipe an existing ROI/line. Before this proactive reset existed, a stale ROI drawn against a previous, larger image pair silently went out of bounds for a smaller replacement: NumPy doesn't raise on an out-of-range slice — `norm_a[ry0:ry1, rx0:rx1]` silently returns a zero-size array — so the crash surfaced much later and far from the real cause: `SpatialDetailAnalyzer._plot_mask_illustration → _stretch_for_display → np.percentile(empty_array, ...)`. The same unguarded `bgsub[y0:y1, x0:x1]` pattern exists in `power_spectrum.py::_extract_roi` and `edge_analyzer.py::analyze`, so a stale ROI could corrupt those sections too. `MainWindow._on_run()`'s validation (checking `self._roi` against every loaded image's `data.shape` right before `settings["roi"]` is set, clearing it with a `QMessageBox` if it no longer fits) is kept as defense-in-depth for any future code path that changes loaded-image dimensions without going through `_on_image_loaded`, but the normal load→run flow now clears stale state at the source instead of catching it reactively at Run time. |
| `binary_dilation` on a loose sigma-threshold mask amplifies noise, not signal | Growing a boolean mask straight from a threshold cut (e.g. Section 8's nebula mask at 1.7σ) dilates *every* True pixel, including scattered single/few-pixel noise-driven false positives — expected in bulk at a loose sigma cut (~4.5% of pixels at 1.7σ one-sided). Each isolated speck balloons into a `~(2·dilation_px+1)²`-px blob, inflating the mask area by 4x+ and diluting any signal-vs-background metric computed over it. Fix: strip small isolated connected components (`scipy.ndimage.label` + `np.bincount` size filter, same size threshold used for hole-filling) *before* calling `binary_dilation` — see `SpatialDetailAnalyzer._remove_small_objects` / `_fill_small_holes` in `image_filters.py`. Caught by comparing mask pixel counts with dilation on vs off on real (noisy) fixture data — a clean synthetic square mask (no noise) will not reveal this bug. |
| Adding GUI parameter rows clips existing text in the Parameters group | `gui/main_window.py`'s `AnalysisControlPanel.setMaximumHeight(...)` caps the whole control panel's height. Metrics / Parameters / Region & Run are laid out side-by-side (`QHBoxLayout` in `control_panel.py::_build_ui`), so the cap must fit the *tallest* group box's natural content height. Parameters (`control_panel.py`, "2. Parameters") is itself split into two side-by-side `QFormLayout`s — "General / PSF" (`form1`) and "Nebula & Local-Maxima" (`form2`) — so its own height is driven by `max(form1_rows, form2_rows)`, not the total row count across both. When adding a new parameter row, add it to whichever column keeps the two roughly balanced, then re-measure: construct `AnalysisControlPanel` headlessly and read `QGroupBox.sizeHint()` for all three boxes (see the measurement approach used when this split was introduced — a small script that imports the widget, calls `.adjustSize()`, and prints each `findChildren(QGroupBox)` entry's `sizeHint()` — is far more precise than eyeballing a screenshot, and sidesteps OS/DPI screenshot-scaling inconsistencies entirely) rather than assuming a fixed per-row pixel cost. Set `setMaximumHeight(...)` to comfortably cover the tallest of the three measured heights. |
@@ -531,6 +531,7 @@ pytest tests/ --cov=analysis,core,synthetic,report --cov-report=html
| `QToolBar.addWidget(spacer)` with an `Expanding` size policy can make every action after it vanish | Adding a bare `QWidget` spacer (`QSizePolicy.Policy.Expanding` horizontal, tried both `Preferred` and `Expanding` vertical) to right-align a trailing toolbar action made that action disappear from the rendered toolbar entirely — not shifted, not overflowed into a `»` chevron, just absent — confirmed by re-rendering with the spacer removed (the action reappeared immediately, left-aligned after the preceding separator). Root cause not fully isolated; treat any `addWidget(spacer)`-for-right-alignment idiom in a `QToolBar` in this codebase as unverified until the rendered result is actually checked (see the next pitfall for how to check it without a real display). |
| OS-level screenshots of the PyQt6 GUI are unreliable for verifying a layout change | `Get-WindowRect`/`SetWindowPos`/`Screen.Bounds` from a non-DPI-aware PowerShell process and the actual rendered window disagreed with each other by inconsistent ratios (not a single uniform scale factor) on a scaled Windows display, making a window that should fit on-screen appear clipped, and vice versa — even maximizing the window didn't reliably show all of it. Prefer `QWidget.grab()` (or `QMainWindow.grab()`) from a small headless script that constructs the widget/window directly and saves the returned `QPixmap` to a PNG — this renders through Qt's own coordinate system with a consistent `devicePixelRatio`, sidestepping OS/DPI virtualization entirely, and doesn't require a visible window at all. For precise sizing decisions (e.g. tuning a `setMaximumHeight`), read `QWidget.sizeHint()`/`minimumSizeHint()` directly instead of eyeballing a rendered image — see the "Adding GUI parameter rows..." pitfall above for the exact approach. |
| Overview figure's box/marker count silently drifted from its own caption | `EdgeAnalyzer.analyze()`'s gradient-magnitude overview map was given `rois_used` — every *searched* candidate ROI (`N_CANDIDATE_EDGES = EDGE_N_TOP_EDGES * 3`, e.g. 9) — instead of only the edges actually *accepted* into `edges` (capped at `EDGE_N_TOP_EDGES`, e.g. 3), so the map drew 9 cyan boxes while its own caption said "three selected." Root cause: `rois_used` was assigned once, early, from the full candidate list, and never reassigned after low-quality candidates got filtered out of `edges`. Whenever a display figure loops over a list to draw one marker/box per entry, verify that list is the actually-used subset the caption describes, not the broader search/candidate pool that produced it — fixed by reassigning `rois_used = [e["roi_used"] for e in edges]` right after `edges` is finalized, before it's read by `_plot_gradient_map`/stored in the result dict. |
+| A "plausible-looking" derived metric can still be measuring the wrong thing entirely | `EdgeAnalyzer._extract_esf`'s `rotation_angle = -(90.0 - angle_deg)` (present since the file's first commit) looked like a reasonable 90°-complement but actually rotated the edge **horizontally**, not vertically as its own design requires (`esf_raw = nanmean(rotated, axis=0)` averages *down columns*, so the edge must run vertically for that average to stay on one side of the transition). The bug was invisible to `tests/test_analysis/test_edge_analyzer.py` because every test there checked `_esf_quality` (a monotonicity ratio) or structural shape (no NaN, normalized range) — never the actual measured width against a *known* ground truth — and the wrong-orientation artifact (a disc-boundary/interpolation trend) happened to also be smooth and monotonic, so it passed every existing gate while over-measuring width by 7–14x on a synthetic edge with a known Gaussian blur sigma. Diagnosed by rendering `rotate()`'s output for several known angles and inspecting it directly (ASCII-art / value dump, not just numbers) — the same "don't hand-derive `rotate()`'s convention" principle documented above ("Locating a point across `scipy.ndimage.rotate()`..."), applied to the rotation *angle formula* itself rather than just a post-hoc point lookup. When a metric's test suite only checks shape/monotonicity/range properties, add at least one test with an analytically-known true value (`tests/test_analysis/test_edge_analyzer.py::TestEdgeWidthAccuracy`, using the erf-profile width of a Gaussian-blurred step edge) — shape-only checks can pass on a metric that's confidently, monotonically, consistently wrong. |
---
diff --git a/analysis/edge_analyzer.py b/analysis/edge_analyzer.py
index 399bea7..92bfa75 100644
--- a/analysis/edge_analyzer.py
+++ b/analysis/edge_analyzer.py
@@ -194,9 +194,21 @@ def _build_edge_entry(self, image: AstroImage, bgsub: np.ndarray,
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)
+
+ # The ESF scan/edge-orientation guide lines must pivot on the ROI's own
+ # geometric center, not the Sobel-detected gradient peak (edge_info's
+ # center_x/center_y) -- scipy.ndimage.rotate() in _extract_esf always
+ # rotates roi_data about its own array center ((h-1)/2, (w-1)/2), and
+ # start_xy (the "Scan start" marker) is derived by inverse-rotating a
+ # point on that same pivot. Drawing the guide lines through the
+ # gradient-peak point instead left the red marker floating off the
+ # cyan line whenever the peak wasn't exactly at the ROI's center.
+ h_roi, w_roi = roi_data.shape
+ roi_cx_full = x0 + (w_roi - 1) / 2.0
+ roi_cy_full = y0 + (h_roi - 1) / 2.0
edge_info_display = dict(edge_info)
- edge_info_display["center_x"] = xc_full - dx0
- edge_info_display["center_y"] = yc_full - dy0
+ edge_info_display["center_x"] = roi_cx_full - dx0
+ edge_info_display["center_y"] = roi_cy_full - dy0
# ESF start point (position index 0), in the same display-frame
# coordinates as edge_info_display, for the directional marker.
@@ -298,7 +310,21 @@ def _detect_strongest_edge(self, roi_data: np.ndarray) -> dict | None:
def _extract_esf(self, roi_data: np.ndarray,
edge_info: dict) -> tuple[np.ndarray, np.ndarray | None, tuple | None]:
angle_deg = np.degrees(edge_info["angle_rad"])
- rotation_angle = -(90.0 - angle_deg)
+ # Empirically verified (not hand-derived -- see the impulse-trick
+ # comment below for why that matters with rotate()'s sign/handedness):
+ # rendering rotate(roi_data, angle, ...) for known synthetic edge
+ # angles and inspecting the result directly shows rotation_angle =
+ # angle_deg aligns the edge vertically, exactly as this function
+ # needs (esf_raw below averages DOWN COLUMNS, so the edge must run
+ # vertically for that average to stay on one side of the transition
+ # per column). A previous formula, -(90.0 - angle_deg), looked like
+ # a more "natural" complement but actually aligned the edge
+ # HORIZONTALLY instead -- averaging down columns then integrated
+ # straight across the transition, canceling out nearly all of the
+ # real signal and leaving only a disc-boundary/interpolation
+ # artifact (still smooth and monotonic, so _esf_quality's gate never
+ # caught it).
+ rotation_angle = angle_deg
# 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.
diff --git a/analysis/image_filters.py b/analysis/image_filters.py
index 098d3fe..4de8323 100644
--- a/analysis/image_filters.py
+++ b/analysis/image_filters.py
@@ -11,13 +11,13 @@
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
-from scipy.ndimage import generic_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, minimum_filter, median_filter, binary_dilation, binary_fill_holes, label
+from scipy.ndimage import generic_filter, uniform_filter, gaussian_filter, gaussian_laplace, gaussian_gradient_magnitude, map_coordinates, zoom, maximum_filter, binary_dilation, binary_fill_holes, label
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,
- WEBER_KERNEL_SIZES,
+ ENTROPY_KERNEL_SIZES,
XS_LINE_ALPHA, SECTION8_BORDER_CROP_FRACTION, SECTION8_ANALYSIS_CMAP,
XS_SNR_REGION_WIDTH,
SECTION8_LOGRATIO_EPS_PERCENTILE, SECTION8_SCATTER_MAX_SAMPLES,
@@ -25,7 +25,8 @@
SECTION8_NEBULA_MASK_MAX_HOLE_PX,
SECTION8_LOCALMAX_FOOTPRINT_MULT, SECTION8_LOCALMAX_PROMINENCE_PERCENTILE,
SECTION8_LOCALMAX_PRESMOOTH_FRACTION, SECTION8_LOCALMAX_REGION_FRACTION,
- SECTION8_LOCALMAX_DIST_MAX_SAMPLES, SECTION8_LOCALMAX_TOP_PERCENT)
+ SECTION8_LOCALMAX_DIST_MAX_SAMPLES, SECTION8_LOCALMAX_TOP_PERCENT,
+ SECTION8_ENTROPY_N_BINS, SECTION8_ENTROPY_CLIP_PERCENTILE)
MAX_DIM_FOR_STD = 2048 # downsample to this before generic_filter (performance)
_DISPLAY_SMOOTH_SIGMA = 1.0 # applied to maps before plotting; does NOT affect metrics
@@ -43,7 +44,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,
- weber_kernel_sizes: tuple = WEBER_KERNEL_SIZES,
+ entropy_kernel_sizes: tuple = ENTROPY_KERNEL_SIZES,
crosshair: dict | None = None,
roi: tuple | None = None,
xs_snr_width: int | None = None,
@@ -78,8 +79,8 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None,
"localmax_region_fraction": localmax_region_fraction,
"localmax_presmooth_fraction": localmax_presmooth_fraction,
"localmax_top_percent": localmax_top_percent,
- "weber_contrast_a": {},
- "weber_contrast_b": {},
+ "entropy_contrast_ratio_a": {},
+ "entropy_contrast_ratio_b": {},
"panels": {},
"localmax": {},
"nc_shared_nebula_pixels": 0,
@@ -92,9 +93,9 @@ def analyze(self, image_a: AstroImage, image_b: AstroImage | None = None,
"wavelet_nc_score_a": {}, "wavelet_nc_score_b": {},
"wavelet_nc_noise_a": {}, "wavelet_nc_noise_b": {}, "wavelet_nc_ratio": {},
"wavelet_nc_neb_std_a": {}, "wavelet_nc_neb_std_b": {}, "wavelet_nc_ratio_err": {},
- "weber_nc_score_a": {}, "weber_nc_score_b": {},
- "weber_nc_noise_a": {}, "weber_nc_noise_b": {}, "weber_nc_ratio": {},
- "weber_nc_neb_std_a": {}, "weber_nc_neb_std_b": {}, "weber_nc_ratio_err": {},
+ "entropy_nc_score_a": {}, "entropy_nc_score_b": {},
+ "entropy_nc_noise_a": {}, "entropy_nc_noise_b": {}, "entropy_nc_ratio": {},
+ "entropy_nc_neb_std_a": {}, "entropy_nc_neb_std_b": {}, "entropy_nc_ratio_err": {},
"gm_nc_score_a": {}, "gm_nc_score_b": {},
"gm_nc_noise_a": {}, "gm_nc_noise_b": {}, "gm_nc_ratio": {},
"gm_nc_neb_std_a": {}, "gm_nc_neb_std_b": {}, "gm_nc_ratio_err": {},
@@ -173,7 +174,7 @@ def _clip01(v): return max(0.0, min(1.0, v))
# same input images.
rng = np.random.default_rng(42)
- # Export the exact preprocessed array every std/LoG/wavelet/Weber/gradient
+ # Export the exact preprocessed array every std/LoG/wavelet/entropy/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.
@@ -231,7 +232,7 @@ def _clip01(v): return max(0.0, min(1.0, v))
if orig_corr_fig is not None:
figures["corr_original"] = fig_to_b64(orig_corr_fig, dpi=150)
- # 1-5. Local std, LoG, wavelet, Weber, gradient — all read norm_a/norm_b with no
+ # 1-5. Local std, LoG, wavelet, entropy, 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:
@@ -277,13 +278,15 @@ def _clip01(v): return max(0.0, min(1.0, v))
localmax_presmooth_fraction=localmax_presmooth_fraction,
localmax_top_percent=localmax_top_percent,
)
- _f_web = _ex.submit(self._weber_analysis,
+ _f_ent = _ex.submit(self._entropy_analysis,
analysis_a, analysis_b,
- weber_kernel_sizes,
+ mask_neb_a, mask_bg_a,
+ mask_neb_b, mask_bg_b,
+ entropy_kernel_sizes,
image_a.label, _label_b,
display_roi=display_roi,
crosshair=crosshair_roi,
- mask_neb_shared=mask_neb_shared, mask_bg_a=mask_bg_a, mask_bg_b=mask_bg_b,
+ mask_neb_shared=mask_neb_shared,
mask_bg_shared=mask_bg_shared, rng=rng,
localmax_footprint_mult=localmax_footprint_mult,
localmax_prominence_percentile=localmax_prominence_percentile,
@@ -307,13 +310,13 @@ def _clip01(v): return max(0.0, min(1.0, v))
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()
+ ent_b64, ent_partial = _f_ent.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(ent_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"])
@@ -321,22 +324,22 @@ def _clip01(v): return max(0.0, min(1.0, v))
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["weber_contrast_a"].update(web_partial["weber_contrast_a"])
- result["weber_contrast_b"].update(web_partial["weber_contrast_b"])
+ result["entropy_contrast_ratio_a"].update(ent_partial["entropy_contrast_ratio_a"])
+ result["entropy_contrast_ratio_b"].update(ent_partial["entropy_contrast_ratio_b"])
result["panels"].update(std_partial["panels"])
result["panels"].update(log_partial["panels"])
result["panels"].update(wav_partial["panels"])
- result["panels"].update(web_partial["panels"])
+ result["panels"].update(ent_partial["panels"])
result["panels"].update(grad_partial["panels"])
result["localmax"].update(std_partial["localmax"])
result["localmax"].update(log_partial["localmax"])
result["localmax"].update(wav_partial["localmax"])
- result["localmax"].update(web_partial["localmax"])
+ result["localmax"].update(ent_partial["localmax"])
result["localmax"].update(grad_partial["localmax"])
# 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),
+ ("wavelet", wav_partial), ("entropy", ent_partial),
("gm", grad_partial)):
for suffix in ("nc_score_a", "nc_score_b", "nc_noise_a", "nc_noise_b",
"nc_neb_std_a", "nc_neb_std_b"):
@@ -351,14 +354,14 @@ def _clip01(v): return max(0.0, min(1.0, v))
if image_b is not None:
nc_errors_by_method = {
"std": result["std_nc_ratio_err"], "log": result["log_nc_ratio_err"],
- "wavelet": result["wavelet_nc_ratio_err"], "weber": result["weber_nc_ratio_err"],
+ "wavelet": result["wavelet_nc_ratio_err"], "entropy": result["entropy_nc_ratio_err"],
"gradient": result["gm_nc_ratio_err"],
}
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"],
+ "entropy": result["entropy_nc_ratio"],
"gradient": result["gm_nc_ratio"],
}, nc_errors_by_method)
if nc_fig is not None:
@@ -368,14 +371,14 @@ def _clip01(v): return max(0.0, min(1.0, v))
"std": std_partial["localmax_log_ratio"],
"log": log_partial["localmax_log_ratio"],
"wavelet": wav_partial["localmax_log_ratio"],
- "weber": web_partial["localmax_log_ratio"],
+ "entropy": ent_partial["localmax_log_ratio"],
"gradient": grad_partial["localmax_log_ratio"],
}
localmax_log_ratio_errors_by_method = {
"std": std_partial["localmax_log_ratio_err"],
"log": log_partial["localmax_log_ratio_err"],
"wavelet": wav_partial["localmax_log_ratio_err"],
- "weber": web_partial["localmax_log_ratio_err"],
+ "entropy": ent_partial["localmax_log_ratio_err"],
"gradient": grad_partial["localmax_log_ratio_err"],
}
lm_ratio_fig = self._plot_localmax_ratio_overview(
@@ -394,7 +397,7 @@ def _clip01(v): return max(0.0, min(1.0, v))
("|LoG|", [(f"log_{s}", float(s), f"|LoG| — σ={s} px") for s in log_sigmas]),
("Gradient |G|", [(f"gradient_{s}", float(s), f"Gradient |G| — σ={s} px") for s in log_sigmas]),
("Wavelet", [(f"wavelet_{lvl}", float(2 ** lvl), f"Wavelet — level {lvl}") for lvl in (2, 3)]),
- ("Weber contrast", [(f"weber_{ks}px", float(ks), f"Weber — {ks} px") for ks in weber_kernel_sizes]),
+ ("Local entropy", [(f"entropy_{ks}px", float(ks), f"Entropy — {ks} px") for ks in entropy_kernel_sizes]),
]
grid_rows = []
for family_label, entries in _grid_families:
@@ -920,7 +923,7 @@ def _log_ratio_map(a: np.ndarray, b: np.ndarray) -> np.ndarray:
background-subtracted flux, wavelet band-pass reconstructions) — this
turns the comparison into "which image shows more structure/contrast at
this scale", not "which image is signed-brighter". Non-negative families
- (std, |LoG|, gradient magnitude, Weber contrast) are unaffected by the
+ (std, |LoG|, gradient magnitude, local entropy) are unaffected by the
abs() since they're already >= 0.
Epsilon-floors both operands using a low percentile (not a raw minimum —
@@ -1422,14 +1425,16 @@ def _reconstruct_level(self, coeffs, target_coeff_idx: int,
# ------------------------------------------------------------------
# ------------------------------------------------------------------
- # Weber fraction contrast maps
+ # Local entropy maps
# ------------------------------------------------------------------
- def _weber_analysis(self, norm_a, norm_b, kernel_sizes,
- label_a, label_b,
+ def _entropy_analysis(self, norm_a, norm_b,
+ mask_neb_a, mask_bg_a,
+ mask_neb_b, mask_bg_b,
+ kernel_sizes, label_a, label_b,
display_roi=None,
crosshair=None,
- mask_neb_shared=None, mask_bg_a=None, mask_bg_b=None,
+ mask_neb_shared=None,
mask_bg_shared=None, rng=None,
localmax_footprint_mult=SECTION8_LOCALMAX_FOOTPRINT_MULT,
localmax_prominence_percentile=SECTION8_LOCALMAX_PROMINENCE_PERCENTILE,
@@ -1438,11 +1443,10 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes,
localmax_top_percent=SECTION8_LOCALMAX_TOP_PERCENT) -> 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": {},
- "weber_nc_neb_std_a": {}, "weber_nc_neb_std_b": {},
+ "entropy_contrast_ratio_a": {}, "entropy_contrast_ratio_b": {},
+ "entropy_nc_score_a": {}, "entropy_nc_score_b": {},
+ "entropy_nc_noise_a": {}, "entropy_nc_noise_b": {},
+ "entropy_nc_neb_std_a": {}, "entropy_nc_neb_std_b": {},
"panels": {},
"localmax": {},
"localmax_log_ratio": {},
@@ -1450,96 +1454,98 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes,
}
single = norm_b is None
for ks in kernel_sizes:
- wc_a = self._compute_weber_map(norm_a, ks)
- wc_b = self._compute_weber_map(norm_b, ks) if not single else None
+ ent_a = self._compute_entropy_map(norm_a, ks)
+ ent_b = self._compute_entropy_map(norm_b, ks) if not single else None
- # 99th percentile avoids dark-sky floor driving the scalar metric to extremes
- partial["weber_contrast_a"][ks] = float(np.percentile(wc_a, 99))
+ # Contrast ratios (computed on unsmoothed maps)
+ cr_a = self._contrast_ratio(ent_a, mask_neb_a, mask_bg_a)
+ partial["entropy_contrast_ratio_a"][ks] = cr_a
if not single:
- partial["weber_contrast_b"][ks] = float(np.percentile(wc_b, 99))
+ cr_b = self._contrast_ratio(ent_b, mask_neb_b, mask_bg_b)
+ partial["entropy_contrast_ratio_b"][ks] = cr_b
noise_a = noise_b = None
if not single:
- nc_a, noise_a, neb_std_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
- partial["weber_nc_neb_std_a"][ks] = neb_std_a
- nc_b, noise_b, neb_std_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["weber_nc_neb_std_b"][ks] = neb_std_b
-
- diff = self._log_ratio_map(wc_a, wc_b) if wc_b is not None else None
- partial["panels"][f"weber_{ks}px"] = {
- "a": wc_a.astype(np.float32),
- "b": wc_b.astype(np.float32) if wc_b is not None else None,
+ nc_a, noise_a, neb_std_a = self._nc_score(ent_a, mask_neb_shared, mask_bg_a)
+ partial["entropy_nc_score_a"][ks] = nc_a
+ partial["entropy_nc_noise_a"][ks] = noise_a
+ partial["entropy_nc_neb_std_a"][ks] = neb_std_a
+ nc_b, noise_b, neb_std_b = self._nc_score(ent_b, mask_neb_shared, mask_bg_b)
+ partial["entropy_nc_score_b"][ks] = nc_b
+ partial["entropy_nc_noise_b"][ks] = noise_b
+ partial["entropy_nc_neb_std_b"][ks] = neb_std_b
+
+ diff = self._log_ratio_map(ent_a, ent_b) if not single else None
+ partial["panels"][f"entropy_{ks}px"] = {
+ "a": ent_a.astype(np.float32),
+ "b": ent_b.astype(np.float32) if ent_b is not None else None,
"diff": diff,
}
if diff is not None:
lm_entry = self._localmax_entry(
- wc_a, wc_b, diff, ks,
+ ent_a, ent_b, diff, ks,
localmax_footprint_mult, localmax_prominence_percentile,
localmax_region_fraction, localmax_presmooth_fraction,
localmax_top_percent, rng)
- partial["localmax"][f"weber_{ks}px"] = lm_entry
+ partial["localmax"][f"entropy_{ks}px"] = lm_entry
partial["localmax_log_ratio"][ks] = lm_entry["log_ratio_mean"]
partial["localmax_log_ratio_err"][ks] = lm_entry["log_ratio_std"]
if not single:
corr_fig = self._plot_metric_correlation(
- wc_a, wc_b, diff, mask_neb_shared, mask_bg_shared,
- label_a, label_b, f"Weber contrast (kernel {ks}px)", rng)
+ ent_a, ent_b, diff, mask_neb_shared, mask_bg_shared,
+ label_a, label_b, f"Local entropy (kernel {ks}px)", rng)
if corr_fig is not None:
- figures[f"corr_weber_{ks}px"] = corr_fig
+ figures[f"corr_entropy_{ks}px"] = corr_fig
if not single and noise_a and noise_b:
- partial["panels"][f"nrm_weber_{ks}px"] = {
- "a": (wc_a / noise_a).astype(np.float32),
- "b": (wc_b / noise_b).astype(np.float32),
+ partial["panels"][f"nrm_entropy_{ks}px"] = {
+ "a": (ent_a / noise_a).astype(np.float32),
+ "b": (ent_b / noise_b).astype(np.float32),
"diff": None,
}
xs_raw = None
xs_line = None
if crosshair is not None and not single:
- pos, pa = self._sample_line(wc_a, **crosshair)
- _, pb = self._sample_line(wc_b, **crosshair)
+ pos, pa = self._sample_line(ent_a, **crosshair)
+ _, pb = self._sample_line(ent_b, **crosshair)
xs_raw = (pos, pa, pb, label_a, label_b,
- f"Cross-section — Weber contrast, kernel {ks}px")
- xs_line = self._crosshair_to_cropped_px(crosshair, wc_a.shape, SECTION8_BORDER_CROP_FRACTION)
+ f"Cross-section — Local entropy, kernel {ks}px")
+ xs_line = self._crosshair_to_cropped_px(crosshair, ent_a.shape, SECTION8_BORDER_CROP_FRACTION)
if not single:
fig = self._plot_side_by_side(
- self._crop_border(wc_a, SECTION8_BORDER_CROP_FRACTION),
- self._crop_border(wc_b, SECTION8_BORDER_CROP_FRACTION),
- f"Weber contrast — kernel {ks}px — {label_a}",
- f"Weber contrast — kernel {ks}px — {label_b}",
- diff_title=f"Weber log-ratio (A/B), kernel {ks}px",
+ self._crop_border(ent_a, SECTION8_BORDER_CROP_FRACTION),
+ self._crop_border(ent_b, SECTION8_BORDER_CROP_FRACTION),
+ f"Local entropy — kernel {ks}px — {label_a}",
+ f"Local entropy — kernel {ks}px — {label_b}",
+ diff_title=f"Log ratio (A/B), kernel {ks}px",
cmap=SECTION8_ANALYSIS_CMAP,
- nonlinear_norm=True, # unbounded output; PowerNorm compresses dynamic range
+ nonlinear_norm=True,
display_roi=None, # _crop_border already applied
xs_data=xs_raw,
xs_line=xs_line,
)
else:
fig = self._plot_single(
- self._crop_border(wc_a, SECTION8_BORDER_CROP_FRACTION),
- f"Weber contrast — kernel {ks}px — {label_a}",
+ self._crop_border(ent_a, SECTION8_BORDER_CROP_FRACTION),
+ f"Local entropy — kernel {ks}px — {label_a}",
cmap=SECTION8_ANALYSIS_CMAP,
nonlinear_norm=True,
)
- figures[f"weber_{ks}px"] = fig
+ figures[f"entropy_{ks}px"] = fig
if not single and noise_a and noise_b:
xs_nrm = None
if crosshair is not None:
- pos_n, pa_n = self._sample_line(wc_a / noise_a, **crosshair)
- _, pb_n = self._sample_line(wc_b / noise_b, **crosshair)
+ pos_n, pa_n = self._sample_line(ent_a / noise_a, **crosshair)
+ _, pb_n = self._sample_line(ent_b / noise_b, **crosshair)
xs_nrm = (pos_n, pa_n, pb_n, label_a, label_b,
- f"Cross-section — Weber contrast (× noise floor), kernel {ks}px")
- figures[f"nrm_weber_{ks}px"] = self._plot_side_by_side(
- self._crop_border(wc_a / noise_a, SECTION8_BORDER_CROP_FRACTION),
- self._crop_border(wc_b / noise_b, SECTION8_BORDER_CROP_FRACTION),
- f"Weber (× noise floor) — kernel {ks}px — {label_a}",
- f"Weber (× noise floor) — kernel {ks}px — {label_b}",
+ f"Cross-section — Local entropy (× noise floor), kernel {ks}px")
+ figures[f"nrm_entropy_{ks}px"] = self._plot_side_by_side(
+ self._crop_border(ent_a / noise_a, SECTION8_BORDER_CROP_FRACTION),
+ self._crop_border(ent_b / noise_b, SECTION8_BORDER_CROP_FRACTION),
+ f"Local entropy (× noise floor) — kernel {ks}px — {label_a}",
+ f"Local entropy (× noise floor) — kernel {ks}px — {label_b}",
diff_title=f"Log ratio (A/B), noise-normalised, kernel {ks}px",
cmap=SECTION8_ANALYSIS_CMAP,
display_roi=None,
@@ -1549,12 +1555,27 @@ def _weber_analysis(self, norm_a, norm_b, kernel_sizes,
return figs_to_b64(figures, dpi=150), partial
- 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.
+ def _compute_entropy_map(self, norm: np.ndarray, kernel_size: int,
+ n_bins: int = SECTION8_ENTROPY_N_BINS) -> np.ndarray:
+ """Local Shannon entropy (bits, log2) of the gray-level histogram within a
+ square window. norm is deliberately quantized into n_bins integer levels
+ from a percentile-clipped range *before* filtering -- computing entropy
+ directly on continuous float32 data returns ~log2(window_area) almost
+ everywhere (every pixel value in a window is unique), measuring float
+ precision rather than genuine tonal diversity. The clip/bin range is
+ computed independently per image (matches the existing per-image
+ convention used by every other Section 8 map).
+
+ Vectorized via n_bins uniform_filter passes (one per gray level) rather
+ than a per-pixel generic_filter callback: uniform_filter(mask_b) is a
+ box filter of a 0/1 mask, which *is* exactly the local proportion p_b
+ of bin b within each window -- summing -p_b*log2(p_b) across bins then
+ gives the entropy map with no per-pixel Python-level histogram ever
+ built. A per-pixel generic_filter callback (the initial approach here,
+ mirroring _compute_std_map) measured ~10-50x slower than Weber's fully
+ vectorized maximum/minimum/median_filter maps at the same array size;
+ this reuses that same vectorization principle for entropy.
"""
- _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:
@@ -1564,20 +1585,24 @@ def _compute_weber_map(self, norm: np.ndarray, kernel_size: int) -> np.ndarray:
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)
- 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)
+ lo, hi = np.percentile(data, [SECTION8_ENTROPY_CLIP_PERCENTILE,
+ 100 - SECTION8_ENTROPY_CLIP_PERCENTILE])
+ if hi <= lo:
+ entropy_map = np.zeros_like(data, dtype=np.float64)
+ else:
+ binned = np.clip(((data - lo) / (hi - lo) * n_bins), 0, n_bins - 1).astype(np.intp)
+ entropy_map = np.zeros(data.shape, dtype=np.float64)
+ for b in range(n_bins):
+ p_b = uniform_filter((binned == b).astype(np.float64), size=kernel_size, mode="reflect")
+ with np.errstate(divide="ignore", invalid="ignore"):
+ entropy_map -= np.where(p_b > 0, p_b * np.log2(p_b), 0.0)
if factor < 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)
+ entropy_map = zoom(entropy_map,
+ (norm.shape[0] / entropy_map.shape[0],
+ norm.shape[1] / entropy_map.shape[1]),
+ order=1)
+ return entropy_map.astype(np.float32)
def _plot_side_by_side(self, arr_a: np.ndarray, arr_b: np.ndarray,
title_a: str, title_b: str,
@@ -1865,12 +1890,12 @@ def _plot_nc_ratio_overview(self, ratios_by_method: dict,
approximate symmetric uncertainty (see _compute_nc_ratio_errors),
rendered as error bars when present for a given point."""
_SCALE_LABEL = {
- "std": "px", "weber": "px", "log": "σ px",
+ "std": "px", "entropy": "px", "log": "σ px",
"gradient": "σ px", "wavelet": "level (≈px)",
}
_COLORS = {
"std": "steelblue", "log": "tomato", "wavelet": "mediumpurple",
- "weber": "seagreen", "gradient": "goldenrod",
+ "entropy": "seagreen", "gradient": "goldenrod",
}
series = self._ratio_series_with_errors(ratios_by_method, errors_by_method)
if not series:
@@ -1913,12 +1938,12 @@ def _plot_localmax_ratio_overview(self, log_ratios_by_method: dict,
measure, not an approximation), rendered as error bars when present for
a given point."""
_SCALE_LABEL = {
- "std": "px", "weber": "px", "log": "σ px",
+ "std": "px", "entropy": "px", "log": "σ px",
"gradient": "σ px", "wavelet": "level (≈px)",
}
_COLORS = {
"std": "steelblue", "log": "tomato", "wavelet": "mediumpurple",
- "weber": "seagreen", "gradient": "goldenrod",
+ "entropy": "seagreen", "gradient": "goldenrod",
}
series = self._ratio_series_with_errors(log_ratios_by_method, errors_by_method)
if not series:
diff --git a/core/fig_utils.py b/core/fig_utils.py
index 520d646..b92208c 100644
--- a/core/fig_utils.py
+++ b/core/fig_utils.py
@@ -2,9 +2,24 @@
import base64
import io
+import threading
import matplotlib.pyplot as plt
+# matplotlib's mathtext grammar (matplotlib/_mathtext.py) calls pyparsing's
+# ParserElement.enable_packrat() at import time, turning on a process-wide,
+# non-thread-safe memoization cache used by every mathtext parse (tick labels,
+# titles, legends -- anything rendered through savefig()/draw()). Multiple
+# analyzers render figures concurrently (SpatialDetailAnalyzer's own 5-way
+# ThreadPoolExecutor, PowerSpectrumAnalyzer's per-image ThreadPoolExecutor, and
+# analysis_thread.py's cross-analyzer parallel mode), so unsynchronized
+# savefig() calls can corrupt that shared cache and raise a spurious
+# mathtext ParseException on completely valid text (reproduced directly
+# against matplotlib.mathtext.MathTextParser.parse() under thread
+# concurrency). Serializing every savefig() through one process-wide lock
+# eliminates the race.
+_SAVEFIG_LOCK = threading.Lock()
+
def fig_to_b64(fig: plt.Figure, dpi: int = 120) -> str:
"""Render a matplotlib figure to a base64 PNG string and immediately close it.
@@ -13,7 +28,8 @@ def fig_to_b64(fig: plt.Figure, dpi: int = 120) -> str:
are released as soon as their pixels are captured.
"""
buf = io.BytesIO()
- fig.savefig(buf, format="png", dpi=dpi, bbox_inches="tight")
+ with _SAVEFIG_LOCK:
+ fig.savefig(buf, format="png", dpi=dpi, bbox_inches="tight")
buf.seek(0)
data = base64.b64encode(buf.read()).decode()
plt.close(fig)
diff --git a/core/models.py b/core/models.py
index 0bdd00a..6c4aa72 100644
--- a/core/models.py
+++ b/core/models.py
@@ -32,7 +32,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)
-WEBER_KERNEL_SIZES = (3, 5, 9) # px; local kernel for Weber fraction contrast c = ΔL/L (odd values required)
+ENTROPY_KERNEL_SIZES = (5, 9, 17) # px; local window for Shannon entropy — deliberately larger than STD_KERNEL_SIZES since small windows (<25 samples) give unstable histogram estimates
WAVELET_NAME = "db4"
WAVELET_LEVELS = 4
@@ -43,6 +43,8 @@
SECTION8_ANALYSIS_CMAP = "viridis" # colormap for Section 8 A/B analysis map panels (std, LoG, wavelet)
SECTION8_LOGRATIO_EPS_PERCENTILE = 1.0 # percentile of pooled positive |A|,|B| values used as the epsilon floor in log10(|A|/|B|)
SECTION8_SCATTER_MAX_SAMPLES = 50000 # per masked population, per scale — caps render cost of Section 8g correlation scatter plots
+SECTION8_ENTROPY_N_BINS = 32 # gray-level bins for local entropy histograms; max possible entropy = log2(32) = 5 bits
+SECTION8_ENTROPY_CLIP_PERCENTILE = 0.5 # symmetric percentile clip (0.5-99.5) applied to each image's own normalised data before binning, so a few outlier pixels don't blow out the bin range
SECTION8_NEBULA_MASK_SIGMA = 1.7 # ×RMS above background = "Nebula" pixel classification (Section 8 masks); background cut stays fixed at 0.5×RMS
SECTION8_NEBULA_MASK_DILATION_PX = 3 # px; scipy.ndimage.binary_dilation iterations to grow the nebula mask into adjacent dim/dark nebula regions
SECTION8_NEBULA_MASK_MAX_HOLE_PX = 5 # px; enclosed background gaps up to this many pixels per side (area ≤ N²) inside the nebula mask are filled before dilation
diff --git a/report/report_builder.py b/report/report_builder.py
index bcafe84..4ff3cc8 100644
--- a/report/report_builder.py
+++ b/report/report_builder.py
@@ -2,6 +2,7 @@
import base64
import io
+import math
import webbrowser
from datetime import datetime
from pathlib import Path
@@ -16,6 +17,7 @@
from scipy.interpolate import griddata as _griddata
from PIL import Image as _PILImage
+from core.fig_utils import fig_to_b64 as _fig_to_b64
from core.models import (AnalysisResult, HALO_FIT_RADIUS_PX, XS_LINE_ALPHA, GLASS_REFRACTIVE_INDEX,
PSF_SPATIAL_MAP_SIZE, PSF_SPATIAL_MAP_SMOOTH_SIGMA, EDGE_ROI_MAP_INDICATOR_PX,
EDGE_N_TOP_EDGES,
@@ -25,7 +27,8 @@
SECTION8_NEBULA_MASK_MAX_HOLE_PX,
SECTION8_LOCALMAX_FOOTPRINT_MULT, SECTION8_LOCALMAX_PROMINENCE_PERCENTILE,
SECTION8_LOCALMAX_PRESMOOTH_FRACTION, SECTION8_LOCALMAX_REGION_FRACTION,
- SECTION8_LOCALMAX_TOP_PERCENT)
+ SECTION8_LOCALMAX_TOP_PERCENT,
+ SECTION8_ENTROPY_N_BINS, SECTION8_ENTROPY_CLIP_PERCENTILE)
from core.astro_image import AstroImage
_TEST_IMAGE_PATH = Path(__file__).parent.parent / "resources" / "ContrastTestImage.png"
@@ -159,15 +162,6 @@ def _psf_make_map(pts: list, img_h: int, img_w: int) -> "np.ndarray | None":
# ── Helpers ───────────────────────────────────────────────────────────────────
-def _fig_to_b64(fig: plt.Figure, dpi: int = 120) -> str:
- buf = io.BytesIO()
- fig.savefig(buf, format="png", dpi=dpi, bbox_inches="tight")
- buf.seek(0)
- data = base64.b64encode(buf.read()).decode()
- plt.close(fig)
- return data
-
-
def _img_tag(fig: "plt.Figure | str | None", alt: str = "") -> str:
if fig is None:
return ""
@@ -463,9 +457,9 @@ def _draw_boxwhisker(ax, vals_list):
("gradient_6.0", "Gradient |G| — σ=6 px"),
("wavelet_2", "Wavelet — level 2"),
("wavelet_3", "Wavelet — level 3"),
- ("weber_3px", "Weber contrast — 3 px"),
- ("weber_5px", "Weber contrast — 5 px"),
- ("weber_9px", "Weber contrast — 9 px"),
+ ("entropy_5px", "Local entropy — 5 px"),
+ ("entropy_9px", "Local entropy — 9 px"),
+ ("entropy_17px", "Local entropy — 17 px"),
]
# Same order/labels as above minus "original" (no kernel scale) — used to order the
@@ -535,8 +529,7 @@ def _draw_boxwhisker(ax, vals_list):
palette=palette, inner=None, linewidth=0.8, ax=ax)
_draw_boxwhisker(ax, [va, vb])
- # Some detail maps (e.g. Weber contrast, unbounded near dark-sky pixels —
- # see 8h methodology) have rare extreme-outlier magnitudes that stretch the
+ # Some detail maps have rare extreme-outlier magnitudes that stretch the
# axis so far the IQR box becomes an invisible sliver. Clip the *view* to the
# 1st-99th percentile of this row's own pooled data; the boxplot/violin
# values themselves are unaffected, only what's visible is cropped.
@@ -568,8 +561,8 @@ def _draw_boxwhisker(ax, vals_list):
"this figure's subsampled copy). "
"Each row's x-axis is independently clipped to its own 1st–99th "
"percentile range so the IQR box stays visible; rows with rare "
- "extreme-outlier magnitudes (e.g. Weber contrast near dark-sky pixels) may "
- "have a small fraction of the violin's tail extend beyond the visible axis."
+ "extreme-outlier magnitudes may have a small fraction of the violin's tail "
+ "extend beyond the visible axis."
""
)
return img_html, caption_html
@@ -810,16 +803,43 @@ def _localmax_rows(localmax: dict, rows: list, val_fmt: str = ".3f") -> str:
'
Detail by scale band, plus explicit SNR '
' Structure whose size matches this scale band; level 1 is noise-only, used to '
' calibrate the noise floor '
- " '
+ " Weber fraction contrast (8h) 3, 5, 9 px window "
- " Contrast (Weber's law, ΔL/L) "
- ' Local range vs. local median background — the most literal "Contrast" metric in '
- ' this section '
+ ' Local entropy map (8h) 5, 9, 17 px window "
+ " Detail (texture complexity, bits, log₂) "
+ ' Rich/unpredictable local tonal structure — tangled nebulosity, mottled dust, '
+ ' unresolved star fields — but noise raises this too, even more readily than σ '
' Entropy contrast ratio (8h) same kernel sizes '
+ ' Contrast, built from the entropy map '
+ ' Whether the nebula has richer tonal/textural complexity than blank sky, at this '
+ ' scale '
'Noise-corrected (NC) score (8d–8h, 8i) same scale as parent method '
' SNR of Detail '
" Whether a detail-map response in the nebula is real structure or just this image's "
' own noise floor at that scale '
+ '
'
+ ' '
+ ' Metric High when Weakness '
+ ' Local σ (8g) Pixels differ strongly from the local mean '
+ ' Responds strongly to noise, halos, gradients, bright stars '
+ ' |LoG| / Gradient (8d, 8f) There are edges, curvature, transitions '
+ ' More shape/edge-biased than texture-complexity-biased '
+ ' Wavelets (8e) Structure exists in a specific spatial-frequency band '
+ ' Needs noise calibration (provided here via the explicit per-level SNR) '
+ 'Local entropy (8h) The local intensity distribution is rich / unpredictable '
+ ' Ignores spatial arrangement entirely, and responds to noise even more readily '
+ ' than σ 8. Spatial Detail Comparison ✓ bandwidth-normalised
@@ -4546,30 +4568,44 @@ def _family_nrm_figs(rows) -> str:
'each map figure below, showing how local detail amplitude varies along the selected line.',
title="Local standard deviation")}
{_std_images_box}
-8h. Weber Fraction Contrast Maps
+8h. Local Entropy Maps
{_info_box(
- '8i. Noise-Corrected Contrast — Cross-Method Overview
{_hires_img_tag(figs.get("nc_ratio_overview"), "NC ratio overview")}
@@ -5014,8 +5050,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 {}
- wc_a_s = sm_a.get("weber_contrast_a", {})
- wc_b_s = sm_b.get("weber_contrast_b", {}) if sm_b else {}
+ ecr_a_s = sm_a.get("entropy_contrast_ratio_a", {})
+ ecr_b_s = sm_b.get("entropy_contrast_ratio_b", {}) if sm_b else {}
snr_ma = ra.snr_metrics or {}
snr_mb = rb.snr_metrics or {}
@@ -5073,7 +5109,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("Weber contrast 99th pct (5px)", wc_a_s.get(5), wc_b_s.get(5), fmt=".4f"),
+ row("Entropy contrast ratio (9px)", ecr_a_s.get(9), ecr_b_s.get(9)),
*([row(
"Global SNR — starless (σ) ★",
(snr_ma.get("starless") or {}).get("snr_global"),
diff --git a/tests/test_analysis/test_edge_analyzer.py b/tests/test_analysis/test_edge_analyzer.py
index 2d2ff38..2c6e15c 100644
--- a/tests/test_analysis/test_edge_analyzer.py
+++ b/tests/test_analysis/test_edge_analyzer.py
@@ -4,6 +4,7 @@
import numpy as np
import pytest
from scipy.ndimage import gaussian_filter
+from scipy.special import erfinv
from analysis.edge_analyzer import EdgeAnalyzer
from core.models import EDGE_ESF_MIN_MONOTONICITY
@@ -11,7 +12,8 @@
_RESULT_KEYS = {"edges", "n_edges", "rois_used"}
-def _make_clean_edge_roi(angle_deg: float = 30.0, size: int = 60) -> np.ndarray:
+def _make_clean_edge_roi(angle_deg: float = 30.0, size: int = 60,
+ sigma: float = 1.5) -> 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]
@@ -19,7 +21,7 @@ def _make_clean_edge_roi(angle_deg: float = 30.0, size: int = 60) -> np.ndarray:
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)
+ return gaussian_filter(roi, sigma=sigma)
def _make_double_edge_roi(size: int = 60) -> np.ndarray:
@@ -125,6 +127,34 @@ def test_esf_normalised_to_unit_range(self):
assert esf.max() <= 1.0 + 1e-9
+class TestEdgeWidthAccuracy:
+ """Regression guard for the rotation_angle sign bug in _extract_esf: a
+ step edge blurred by a known Gaussian sigma has an analytically known
+ 10-90% width (the edge spread function of a Gaussian-blurred step is an
+ erf profile), so the measured width can be checked against ground truth
+ instead of only the shape/monotonicity properties the rest of this file
+ tests. A prior formula (-(90.0 - angle_deg), which looked plausible but
+ aligned the edge horizontally instead of vertically) passed every
+ existing test in this file while over-measuring width by 7-14x, because
+ none of them compared against a known true width."""
+
+ @staticmethod
+ def _expected_width(sigma: float) -> float:
+ # ESF(x) = 0.5*(1+erf(x/(sigma*sqrt2))); solve for the 10%/90% crossings.
+ return 2.0 * sigma * np.sqrt(2.0) * erfinv(0.8)
+
+ @pytest.mark.parametrize("angle_deg", [15.0, 30.0, 60.0, 75.0])
+ @pytest.mark.parametrize("sigma", [1.5, 3.0])
+ def test_measured_width_matches_known_sigma(self, angle_deg, sigma):
+ ea = EdgeAnalyzer()
+ roi = _make_clean_edge_roi(angle_deg=angle_deg, sigma=sigma)
+ edge_info = ea._detect_strongest_edge(roi)
+ positions, esf, _ = ea._extract_esf(roi, edge_info)
+ width = ea._measure_edge_width(positions, esf)
+ expected = self._expected_width(sigma)
+ assert width == pytest.approx(expected, rel=0.3)
+
+
class TestQualityGateAutoDetect:
"""Directly control which candidate ROIs _auto_detect_top_rois returns so
the skip/fallback control flow can be tested deterministically, without
diff --git a/tests/test_analysis/test_spatial_detail.py b/tests/test_analysis/test_spatial_detail.py
index 6afc603..616f386 100644
--- a/tests/test_analysis/test_spatial_detail.py
+++ b/tests/test_analysis/test_spatial_detail.py
@@ -7,7 +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
+from core.models import STD_KERNEL_SIZES, LOG_SIGMAS, ENTROPY_KERNEL_SIZES, WAVELET_LEVELS
class TestAnalyze:
@@ -65,20 +65,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):
+ def test_entropy_contrast_ratio_a_present(self, astro_image_a):
result = SpatialDetailAnalyzer().analyze(astro_image_a)
- assert "weber_contrast_a" in result
+ assert "entropy_contrast_ratio_a" in result
- def test_single_image_weber_contrast_b_empty(self, astro_image_a):
+ def test_single_image_entropy_contrast_ratio_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
+ # Single-image mode: entropy_contrast_ratio_b is present but empty (same pattern as contrast_ratios_b)
+ ecr_b = result.get("entropy_contrast_ratio_b")
+ assert ecr_b is None or not ecr_b
- def test_weber_contrast_a_positive(self, astro_image_a):
+ def test_entropy_contrast_ratio_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
+ for v in result.get("entropy_contrast_ratio_a", {}).values():
+ if v is not None:
+ assert v >= 0.0
def test_with_roi(self, astro_image_a):
result = SpatialDetailAnalyzer().analyze(astro_image_a,
@@ -93,7 +94,7 @@ def test_with_roi(self, astro_image_a):
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/
+ amplitude 8x sky noise) is added inside the blob so std/LoG/wavelet/entropy/
gradient all detect meaningfully more local structure than the plain blob."""
rng = np.random.default_rng(seed)
h, w = 256, 256
@@ -159,7 +160,7 @@ def test_shared_nebula_pixels_positive(self, nc_result):
@pytest.mark.parametrize("prefix,scales", [
("std", STD_KERNEL_SIZES),
("log", LOG_SIGMAS),
- ("weber", WEBER_KERNEL_SIZES),
+ ("entropy", ENTROPY_KERNEL_SIZES),
("gm", LOG_SIGMAS),
])
def test_nc_score_dict_keys_match_scales(self, nc_result, prefix, scales):
@@ -171,7 +172,7 @@ def test_wavelet_nc_score_keys_match_levels(self, nc_result):
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"])
+ @pytest.mark.parametrize("prefix", ["std", "log", "wavelet", "entropy", "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():
@@ -190,9 +191,9 @@ 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_entropy_nc_ratio_captures_finer_detail_in_a(self, nc_result):
+ ratio = nc_result["entropy_nc_ratio"][min(ENTROPY_KERNEL_SIZES)]
+ assert ratio is not None and ratio > 1.0
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
@@ -209,8 +210,8 @@ def test_normalized_panels_present_two_image(self, nc_result):
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 ks in ENTROPY_KERNEL_SIZES:
+ assert f"nrm_entropy_{ks}px" in panels
for lvl in (2, 3):
assert f"nrm_wavelet_{lvl}" in panels
@@ -243,7 +244,7 @@ def test_with_roi_two_image_no_crash(self, nc_image_pair):
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
+ contrast_ratios_b / entropy_contrast_ratio_b invariant (never None/absent, never
populated with placeholder values on the A side either)."""
@pytest.fixture(scope="class")
@@ -255,7 +256,7 @@ def single_result(cls, astro_image_a):
"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",
+ "entropy_nc_score_a", "entropy_nc_score_b", "entropy_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):
@@ -826,7 +827,7 @@ class TestCorrelationScatterFigures:
+ [f"corr_log_{s}" for s in LOG_SIGMAS]
+ [f"corr_gradient_{s}" for s in LOG_SIGMAS]
+ ["corr_wavelet_2", "corr_wavelet_3"]
- + [f"corr_weber_{ks}px" for ks in WEBER_KERNEL_SIZES]
+ + [f"corr_entropy_{ks}px" for ks in ENTROPY_KERNEL_SIZES]
)
@pytest.mark.parametrize("key", _CORR_KEYS)
@@ -850,7 +851,7 @@ class TestLocalMaxIntegration:
+ [f"log_{s}" for s in LOG_SIGMAS]
+ [f"gradient_{s}" for s in LOG_SIGMAS]
+ ["wavelet_2", "wavelet_3"]
- + [f"weber_{ks}px" for ks in WEBER_KERNEL_SIZES]
+ + [f"entropy_{ks}px" for ks in ENTROPY_KERNEL_SIZES]
)
@pytest.mark.parametrize("key", _LOCALMAX_KEYS)
@@ -943,8 +944,8 @@ def nc_result_with_crosshair(nc_image_pair) -> dict:
class TestCrosshairEmbeddedCrossSections:
"""Section 8 cross-sections are embedded panels inside the 2x2 map figures,
not separate figures — no standalone xs_* keys should ever appear, and all
- 5 families (including Weber, newly crosshair-capable) must still produce
- their usual figure keys with or without a crosshair, without crashing."""
+ 5 families (including Local entropy) must still produce their usual figure
+ keys with or without a crosshair, without crashing."""
def test_no_crash_with_crosshair(self, nc_result_with_crosshair):
assert isinstance(nc_result_with_crosshair, dict)
@@ -952,15 +953,15 @@ def test_no_crash_with_crosshair(self, nc_result_with_crosshair):
def test_no_standalone_xs_keys_with_crosshair(self, nc_result_with_crosshair):
figs = nc_result_with_crosshair["figures"]
assert not any(k.startswith(("xs_std_", "xs_log_", "xs_wavelet_",
- "xs_gradient_", "xs_weber_")) for k in figs)
+ "xs_gradient_", "xs_entropy_")) for k in figs)
def test_no_standalone_xs_keys_without_crosshair(self, nc_result):
figs = nc_result["figures"]
assert not any(k.startswith(("xs_std_", "xs_log_", "xs_wavelet_",
- "xs_gradient_", "xs_weber_")) for k in figs)
+ "xs_gradient_", "xs_entropy_")) for k in figs)
@pytest.mark.parametrize("key_prefix,scales", [
- ("std_", STD_KERNEL_SIZES), ("weber_", WEBER_KERNEL_SIZES),
+ ("std_", STD_KERNEL_SIZES), ("entropy_", ENTROPY_KERNEL_SIZES),
])
def test_family_figures_present_with_crosshair(self, nc_result_with_crosshair,
key_prefix, scales):
@@ -968,13 +969,14 @@ def test_family_figures_present_with_crosshair(self, nc_result_with_crosshair,
for scale in scales:
assert f"{key_prefix}{scale}px" in figs
- def test_weber_no_crash_without_crosshair(self, nc_result):
- # Regression: weber previously had no crosshair param at all.
- assert all(f"weber_{ks}px" in nc_result["figures"] for ks in WEBER_KERNEL_SIZES)
+ def test_entropy_no_crash_without_crosshair(self, nc_result):
+ # Regression: weber (entropy's predecessor family) previously had no
+ # crosshair param at all.
+ assert all(f"entropy_{ks}px" in nc_result["figures"] for ks in ENTROPY_KERNEL_SIZES)
- def test_weber_nrm_figures_present_with_crosshair(self, nc_result_with_crosshair):
+ def test_entropy_nrm_figures_present_with_crosshair(self, nc_result_with_crosshair):
figs = nc_result_with_crosshair["figures"]
- assert any(k.startswith("nrm_weber_") for k in figs)
+ assert any(k.startswith("nrm_entropy_") for k in figs)
def test_original_present_with_crosshair(self, nc_result_with_crosshair):
assert "original" in nc_result_with_crosshair["figures"]
@@ -1014,8 +1016,8 @@ class TestSectionSpatialReportOrder:
"""Integration check on report/report_builder.py::_section_spatial's HTML
output: the reorganized subsection order (8a Background, 8b Original Image,
8c Mask Overview, 8d LoG, 8e Wavelet, 8f Gradient, 8g Local Std,
- 8h Weber, 8i NC overview, 8j Local-Maxima Masked Metrics), and the fix for
- the alphabetic-sort bug that put e.g. nrm_std_10px before nrm_std_3px/5px."""
+ 8h Local Entropy, 8i NC overview, 8j Local-Maxima Masked Metrics), and the
+ fix for the alphabetic-sort bug that put e.g. nrm_std_10px before nrm_std_3px/5px."""
@pytest.fixture(scope="class")
@classmethod
@@ -1041,9 +1043,9 @@ def test_local_std_now_in_contrast_group_after_gradient(self, section_html):
assert section_html.index("8f. Gradient Magnitude") < section_html.index(
"8g. Local Standard Deviation Maps")
- def test_weber_after_local_std(self, section_html):
+ def test_entropy_after_local_std(self, section_html):
assert section_html.index("8g. Local Standard Deviation Maps") < section_html.index(
- "8h. Weber Fraction Contrast Maps")
+ "8h. Local Entropy Maps")
def test_nrm_std_figures_in_ascending_kernel_order(self, section_html):
# Regression test: figs_for()'s old lexicographic sorted(figs) rendered
@@ -1104,7 +1106,7 @@ def test_old_mask_grid_illustrative_caption_removed(self, section_html):
"Show Wavelet maps & figures",
"Show Gradient maps & figures",
"Show Local σ maps & figures",
- "Show Weber contrast maps & figures",
+ "Show Local entropy maps & figures",
])
def test_family_images_are_collapsed_by_default(self, section_html, title):
# Each of the 5 metric families' map/correlation/noise-normalised figures
diff --git a/tests/test_report/test_inspector_catalog.py b/tests/test_report/test_inspector_catalog.py
index 87f574b..1fd240d 100644
--- a/tests/test_report/test_inspector_catalog.py
+++ b/tests/test_report/test_inspector_catalog.py
@@ -1,9 +1,10 @@
"""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.
+map (std_15px/std_31px from an old STD_KERNEL_SIZES, no Weber/entropy 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
@@ -68,11 +69,13 @@ def test_every_computed_panel_is_cataloged(self, inspector_npz_path):
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.
+ def test_entropy_panels_present(self, inspector_npz_path):
+ # Regression: the old hardcoded _PANEL_IMAGE_SETS never included this
+ # family at all (Weber, entropy's predecessor family, was the original
+ # motivating example).
_, 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"
+ entropy_keys = [k for k in spatial["panels"] if k.startswith("entropy_")]
+ assert entropy_keys, "fixture should have produced Local entropy panels"
def test_all_std_kernel_sizes_present(self, inspector_npz_path):
# Regression: the old map referenced std_15px/std_31px from a stale
diff --git a/tests/test_report/test_report_helpers.py b/tests/test_report/test_report_helpers.py
index dd2fed2..57cc9e7 100644
--- a/tests/test_report/test_report_helpers.py
+++ b/tests/test_report/test_report_helpers.py
@@ -322,7 +322,7 @@ class TestPanelDisplayName:
@pytest.mark.parametrize("pkey,expected", [
("std_3px", "Std Dev 3 px"),
("std_10px", "Std Dev 10 px"),
- ("weber_9px", "Weber 9 px"),
+ ("entropy_9px", "Entropy 9 px"),
("log_1.5", "LoG σ 1.5 px"),
("wavelet_2", "Wavelet level 2"),
("gradient_3.0", "Gradient σ 3.0 px"),
From d65c64cb52316f5b714d4cca5ffacb25514369bc Mon Sep 17 00:00:00 2001
From: Brent <52629076+brentmantooth@users.noreply.github.com>
Date: Sat, 18 Jul 2026 21:46:34 -0400
Subject: [PATCH 11/13] update spash screen graphic
---
.claude/settings.json | 3 ++-
resources/AstroImageLabSplash.png | Bin 0 -> 2371764 bytes
2 files changed, 2 insertions(+), 1 deletion(-)
create mode 100644 resources/AstroImageLabSplash.png
diff --git a/.claude/settings.json b/.claude/settings.json
index e155565..a2ca02c 100644
--- a/.claude/settings.json
+++ b/.claude/settings.json
@@ -19,7 +19,8 @@
"Bash(\"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_roi_reset.py\")",
"Bash(PYTHONIOENCODING=utf-8 \"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_roi_reset.py\")",
"Bash(PYTHONIOENCODING=utf-8 \"/c/Users/bmant/anaconda3/envs/astrolab/python.exe\" \"/c/Users/bmant/AppData/Local/Temp/claude/d--GitHub-AstroImageLab/689ff469-8361-4066-b1c5-da15c208080d/scratchpad/test_toolbar_sync.py\")",
- "Bash(python -m py_compile gui/control_panel.py gui/main_window.py gui/image_panel.py)"
+ "Bash(python -m py_compile gui/control_panel.py gui/main_window.py gui/image_panel.py)",
+ "Bash(/c/Users/bmant/anaconda3/envs/astrolab/python.exe -c ' *)"
]
}
}
diff --git a/resources/AstroImageLabSplash.png b/resources/AstroImageLabSplash.png
new file mode 100644
index 0000000000000000000000000000000000000000..ff6b3b20ae7a2fd4a674d9b0906f97cc49ff1615
GIT binary patch
literal 2371764
zcmV(`K-0g8P)*Fo?uE1NW*^gaicidBdvgreRjlszyjh&C#~29Z!!G7=j1}jsVI7}29O6g*
z=ugFTcrU*5m;c}K^7JC6(-^1usQIPy;mmtILBw=ud~)!73Yu|K1b&vj;_N2(Y7Qo=
zYs|qffW`!u1ax$XXd)Z}9KX-OmqUY&+8mtfTipTJKvjY#b?_#z1VM)f#zs#vv{E2r
zVzU`2_Fq|e7X>0slJzC7lMl2{fyFw$*VxpRaJ<+55)>2wng9GWX{ev?6QyNL0pBu_
zXrHO$eFrf-O-z>0ankO8y$ $nU)XAImi?hhhU4#XecY(>07@b}ff!
z>nd3f?MUBN9wTU@RZ>IcVY^OocN630@o07g;Ewi-O