From b18d2261344e7408322bd208742c164c2a1b9cb3 Mon Sep 17 00:00:00 2001 From: Peter Jacobson Date: Mon, 15 Jun 2026 12:35:39 +1000 Subject: [PATCH 1/2] Add sub-pixel Gaussian-blur sigma with a physical (nm/FWHM/kernel) readout MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The image-viewer "Smooth: Gaussian" control was an integer slider (1-20 px) and reported nothing about the physical extent of the kernel, so a user could not tell what distance was being blurred or do gentle sub-pixel smoothing. The blur itself (scipy gaussian_filter, sigma in px, kernel truncated at +-4 sigma) is sound and ~equivalent to ImageJ for sigma >= 1 px on float data; the gap was the GUI control and units feedback. Backend already accepts float sigma, so this is GUI-only. Where to look: - probeflow/gui/processing.py: * format_gaussian_readout(sigma_px, px_nm) — pure, unit-tested formatter (FWHM = 2.3548*sigma; kernel half-width = int(4*sigma+0.5), matching scipy truncate=4.0). * _sub_slider() gained a `scale` arg; the smooth slider now spans 0.2-20.0 px in 0.1 steps (_SMOOTH_SIGMA_SCALE). state()/set_state() convert accordingly. * set_pixel_size_nm() + _update_smooth_readout() drive a readout label that shows sigma/FWHM/kernel extent in nm (or px-only when uncalibrated). - probeflow/gui/dialogs/image_viewer.py: _refresh_display_array pushes the loaded scan's pixel size into the panel via set_pixel_size_nm. Tests: test_format_gaussian_readout_pure, test_viewer_full_smooth_sigma_is_subpixel_float. Co-Authored-By: Claude Opus 4.8 --- probeflow/gui/dialogs/image_viewer.py | 5 ++ probeflow/gui/processing.py | 87 ++++++++++++++++++++++++--- tests/test_gui_processing_panel.py | 37 ++++++++++++ 3 files changed, 120 insertions(+), 9 deletions(-) diff --git a/probeflow/gui/dialogs/image_viewer.py b/probeflow/gui/dialogs/image_viewer.py index 7704d14..09e6f40 100644 --- a/probeflow/gui/dialogs/image_viewer.py +++ b/probeflow/gui/dialogs/image_viewer.py @@ -365,6 +365,11 @@ def _processing_pixel_sizes_m(self) -> tuple[float | None, float | None]: return w_m / Nx, h_m / Ny def _refresh_display_array(self, reset_zoom_if_shape_changed: bool = False): + # Keep the Gaussian-blur σ readout calibrated to the loaded scan. + panel = getattr(self, "_processing_panel", None) + if panel is not None and hasattr(panel, "set_pixel_size_nm"): + psx, _psy = self._processing_pixel_sizes_m() + panel.set_pixel_size_nm(psx * 1e9 if psx else None) old_shape = self._display_arr.shape if self._display_arr is not None else None # display array: raw with processing applied (no grain overlay — that's visual only) if self._raw_arr is not None and self._processing: diff --git a/probeflow/gui/processing.py b/probeflow/gui/processing.py index 7a09a5a..b0f2d9a 100644 --- a/probeflow/gui/processing.py +++ b/probeflow/gui/processing.py @@ -10,6 +10,38 @@ QSizePolicy, QSlider, QVBoxLayout, QWidget, ) +# The smooth-σ slider is an integer QSlider scaled by this factor so it can +# express sub-pixel σ (e.g. a slider value of 5 → σ = 0.5 px). +_SMOOTH_SIGMA_SCALE = 10 +# FWHM of a Gaussian = 2·sqrt(2·ln 2)·σ. +_GAUSSIAN_FWHM_FACTOR = 2.35482004503 + + +def format_gaussian_readout(sigma_px: float, px_nm: float | None) -> str: + """One-line description of the Gaussian-blur kernel's physical extent. + + ``sigma_px`` is the standard deviation in pixels. The kernel is truncated at + ±4σ to match scipy's ``gaussian_filter(truncate=4.0)`` default, so the + half-width in pixels is ``int(4σ + 0.5)``. When ``px_nm`` (pixel size in nm) + is known the σ, FWHM and kernel extent are also reported in nanometres. + """ + sigma_px = float(sigma_px) + fwhm_px = _GAUSSIAN_FWHM_FACTOR * sigma_px + r_px = int(4.0 * sigma_px + 0.5) + if px_nm and px_nm > 0: + sigma_nm = sigma_px * px_nm + fwhm_nm = fwhm_px * px_nm + r_nm = r_px * px_nm + return ( + f"σ {sigma_px:.1f} px · {sigma_nm:.3g} nm " + f"FWHM {fwhm_nm:.3g} nm " + f"kernel ±{r_px} px (±{r_nm:.3g} nm)" + ) + return ( + f"σ {sigma_px:.1f} px FWHM {fwhm_px:.2f} px kernel ±{r_px} px" + ) + + class ProcessingControlPanel(QWidget): """Internal processing controls shared by Browse and Viewer.""" @@ -26,6 +58,7 @@ def __init__(self, mode: str, parent=None): if mode not in ("browse_quick", "viewer_full"): raise ValueError(f"Unknown processing panel mode: {mode}") self._mode = mode + self._smooth_px_nm: float | None = None self._build() def _build(self): @@ -57,7 +90,10 @@ def _combo_row(label: str, items: list[str], return cb def _sub_slider(label: str, mn: int, mx: int, init: int, - fmt="{v}") -> tuple[QWidget, QSlider, QLabel]: + fmt="{v}", scale: int = 1) -> tuple[QWidget, QSlider, QLabel]: + # ``scale`` > 1 makes the integer slider represent a fractional value: + # the displayed/logical value is ``slider_value / scale`` (e.g. scale=10 + # gives 0.1 steps). scale=1 is the original integer behaviour. w = QWidget() rl = QHBoxLayout(w) rl.setContentsMargins(0, 0, 0, 0) @@ -67,11 +103,15 @@ def _sub_slider(label: str, mn: int, mx: int, init: int, sl = QSlider(Qt.Horizontal) sl.setRange(mn, mx) sl.setValue(init) - val_lbl = QLabel(fmt.format(v=init)) + + def _disp(v: int) -> str: + return fmt.format(v=(v / scale if scale != 1 else v)) + + val_lbl = QLabel(_disp(init)) val_lbl.setFont(ui_font(8)) - val_lbl.setFixedWidth(28) + val_lbl.setFixedWidth(28 if scale == 1 else 44) sl.valueChanged.connect( - lambda v, vl=val_lbl, f=fmt: vl.setText(f.format(v=v))) + lambda v, vl=val_lbl: vl.setText(_disp(v))) rl.addWidget(lbl) rl.addWidget(sl, 1) rl.addWidget(val_lbl) @@ -239,14 +279,27 @@ def _col_lbl(text: str, target): self._smooth_combo = _combo_row("Smooth:", ["None", "Gaussian"], L, 54) self._smooth_combo.setToolTip( "Gaussian blur to suppress noise. Larger sigma (px) smooths more, but " - "also blurs genuine fine features." + "also blurs genuine fine features. The kernel spans ±4σ; σ may be " + "sub-pixel (down to 0.2 px) for gentle denoising." ) self._smooth_sigma_w, self._smooth_sigma_sl, _ = _sub_slider( - "sigma:", 1, 20, 1, "{v}px") + "sigma:", 2, 200, _SMOOTH_SIGMA_SCALE, "{v:.1f}px", + scale=_SMOOTH_SIGMA_SCALE) L.addWidget(self._smooth_sigma_w) self._smooth_sigma_w.setVisible(False) + # Physical readout: σ/FWHM/kernel extent in nm (when calibrated) or px. + self._smooth_readout_lbl = QLabel() + self._smooth_readout_lbl.setFont(ui_font(7)) + self._smooth_readout_lbl.setWordWrap(True) + self._smooth_readout_lbl.setVisible(False) + L.addWidget(self._smooth_readout_lbl) + self._smooth_sigma_sl.valueChanged.connect( + lambda _v: self._update_smooth_readout()) self._smooth_combo.currentIndexChanged.connect( - lambda i: self._smooth_sigma_w.setVisible(i != 0)) + lambda i: (self._smooth_sigma_w.setVisible(i != 0), + self._smooth_readout_lbl.setVisible(i != 0), + self._update_smooth_readout())) + self._update_smooth_readout() self._highpass_combo = _combo_row("Hi-pass:", ["None", "Gaussian"], L, 54) self._highpass_sigma_w, self._highpass_sigma_sl, _ = _sub_slider( @@ -314,7 +367,10 @@ def state(self) -> dict: "remove_bad_lines_max_adjacent_bad_lines": int( self._bad_line_adjacent_spin.value() ), - "smooth_sigma": self._smooth_sigma_sl.value() if smooth_i != 0 else None, + "smooth_sigma": ( + self._smooth_sigma_sl.value() / _SMOOTH_SIGMA_SCALE + if smooth_i != 0 else None + ), "highpass_sigma": self._highpass_sigma_sl.value() if highpass_i != 0 else None, "edge_method": edge_map[edge_i], "edge_sigma": self._edge_sigma_sl.value(), @@ -348,7 +404,8 @@ def set_state(self, state: dict | None) -> None: sigma = state.get("smooth_sigma") if sigma: self._smooth_combo.setCurrentIndex(1) - self._smooth_sigma_sl.setValue(int(sigma)) + self._smooth_sigma_sl.setValue( + int(round(float(sigma) * _SMOOTH_SIGMA_SCALE))) else: self._smooth_combo.setCurrentIndex(0) @@ -365,6 +422,18 @@ def set_state(self, state: dict | None) -> None: "sobel": 4, "scharr": 5}.get(edge, 0)) self._edge_sigma_sl.setValue(int(state.get("edge_sigma", 1))) + def set_pixel_size_nm(self, px_nm: float | None) -> None: + """Tell the panel the loaded scan's pixel size (nm) for the σ readout.""" + self._smooth_px_nm = float(px_nm) if px_nm else None + self._update_smooth_readout() + + def _update_smooth_readout(self) -> None: + lbl = getattr(self, "_smooth_readout_lbl", None) + if lbl is None: # browse_quick panel has no smooth control + return + sigma_px = self._smooth_sigma_sl.value() / _SMOOTH_SIGMA_SCALE + lbl.setText(format_gaussian_readout(sigma_px, self._smooth_px_nm)) + def bad_line_method(self) -> str | None: return self.state().get("remove_bad_lines") diff --git a/tests/test_gui_processing_panel.py b/tests/test_gui_processing_panel.py index 50b1aaf..63f7f29 100644 --- a/tests/test_gui_processing_panel.py +++ b/tests/test_gui_processing_panel.py @@ -188,6 +188,43 @@ def test_viewer_full_panel_round_trips_standard_processing_state(qapp): ) +def test_format_gaussian_readout_pure(): + from probeflow.gui.processing import format_gaussian_readout + + # Calibrated: σ, FWHM and kernel extent reported in nm (σ_nm = σ_px × px_nm, + # FWHM = 2.3548·σ, kernel half-width = int(4σ+0.5) to match scipy truncate=4). + text = format_gaussian_readout(1.0, 0.5) + assert "σ 1.0 px" in text + assert "0.5 nm" in text # σ_nm + assert "FWHM 1.18 nm" in text # 2.3548 * 0.5 + assert "kernel ±4 px" in text + assert "±2 nm" in text # 4 px × 0.5 nm + + # Uncalibrated: px-only, no nm. + px_only = format_gaussian_readout(2.0, None) + assert "nm" not in px_only + assert "kernel ±8 px" in px_only + + +def test_viewer_full_smooth_sigma_is_subpixel_float(qapp): + from probeflow.gui import ProcessingControlPanel + + panel = ProcessingControlPanel("viewer_full") + + # Sub-pixel σ now round-trips through the GUI slider. + panel.set_state({"smooth_sigma": 0.5}) + assert panel.state()["smooth_sigma"] == pytest.approx(0.5) + + # Physical readout reflects the calibration once the pixel size is known. + panel.set_pixel_size_nm(0.5) + assert "nm" in panel._smooth_readout_lbl.text() + assert "σ 0.5 px" in panel._smooth_readout_lbl.text() + + # Without calibration the readout falls back to pixels only. + panel.set_pixel_size_nm(None) + assert "nm" not in panel._smooth_readout_lbl.text() + + def test_viewer_align_rows_applies_immediately(qapp, monkeypatch): from probeflow.gui import ImageViewerDialog, SxmFile, THEMES From 8394da429dfcdd5be43e46b7549bac889f32c656 Mon Sep 17 00:00:00 2001 From: Peter Jacobson Date: Mon, 15 Jun 2026 12:35:56 +1000 Subject: [PATCH 2/2] Fix bad-line MAD detection: column-aware limit + matched-filter for extended lines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The MAD/outlier bad-line method failed on real STM/AFM scans (createc qPlus AFM) in three compounding ways. All three are fixed in probeflow/processing/bad_lines.py (detect_bad_scanline_segments MAD branch + helpers). 1. Correction could not be applied. _split_segments_by_adjacent_limit was column-blind: a line counted as "bad" if it had a segment anywhere, so noise across many lines formed one long consecutive-line run that exceeded "max adjacent" and ALL segments were skipped (nothing repaired). It is now column-overlap-aware (union-find over segments linked only when adjacent rows overlap in columns); only a genuine vertical stack taller than the limit is skipped. The repair itself was already column-aware. 2. Long segments were invisible/shattered. The method subtracted each row's own median before thresholding, so a defect covering >~half the row became the row median and vanished (lowering the threshold only surfaced noise). It now thresholds the neighbour-referenced residual directly, bridges short noise gaps (_close_small_gaps, _NOISE_GAP_TOL_PX) and keeps a run only if it is a defect on its MEDIAN — robust to per-pixel noise; no upper-length cap. 3. Low-contrast lines were undetectable. The real bright lines are only ~1.5 sigma per pixel, so per-pixel thresholding can't separate them from texture. The residual is now smoothed along each row first (_smooth_rows_nanaware, _MAD_SMOOTH_WINDOW_PX = 11) as a matched filter for extended defects: SNR gain ~sqrt(window) lets a faint-but-long line clear a clean threshold while texture averages out. On the sample file, threshold-5 detection went from 0 to ~80 segments on exactly the bright rows, repaired with 0 skipped. The `step` method is unchanged and remains the sharp/short-scar detector. Tests: test_adjacent_limit_is_column_aware, test_mad_detects_long_partial_segment_against_noise, test_mad_detects_whole_line_offset_not_just_internal_segments, and an updated (now matched-filter-aware) test_mad_outlier_method_repairs_only_segment. Co-Authored-By: Claude Opus 4.8 --- probeflow/processing/bad_lines.py | 166 ++++++++++++++++++++++++------ tests/test_processing.py | 94 +++++++++++++++-- 2 files changed, 221 insertions(+), 39 deletions(-) diff --git a/probeflow/processing/bad_lines.py b/probeflow/processing/bad_lines.py index b611d65..b969941 100644 --- a/probeflow/processing/bad_lines.py +++ b/probeflow/processing/bad_lines.py @@ -5,6 +5,7 @@ from dataclasses import dataclass import numpy as np +from scipy.ndimage import uniform_filter1d from ._image_utils import _nonnegative_finite @@ -38,6 +39,31 @@ class BadLineCorrectionInfo: max_adjacent_bad_lines: int = 1 +# Maximum sub-threshold gap (px) bridged inside a MAD segment so per-pixel +# noise holes don't shatter a coherent defect into fragments. Kept small and +# fixed (a noise-hole scale), independent of the minimum feature length. +_NOISE_GAP_TOL_PX = 4 + +# Along-row smoothing window (px) applied to the residual before MAD detection. +# Acts as a matched filter for *extended* scan-line defects: a coherent line is +# preserved while per-pixel noise averages down (SNR gain ~ sqrt(window)), so a +# faint-but-long bad line clears the threshold without texture flooding in. +_MAD_SMOOTH_WINDOW_PX = 11 + + +def _smooth_rows_nanaware(arr: np.ndarray, window: int) -> np.ndarray: + """Moving average along each row (axis 1), ignoring non-finite pixels.""" + if window <= 1: + return arr + finite = np.isfinite(arr) + vals = np.where(finite, arr, 0.0) + num = uniform_filter1d(vals, size=window, axis=1, mode="nearest") + den = uniform_filter1d(finite.astype(np.float64), size=window, axis=1, mode="nearest") + out = np.full_like(arr, np.nan, dtype=np.float64) + np.divide(num, den, out=out, where=den > 1e-6) + return out + + def _normalise_bad_segment_method(method: str) -> str: method = str(method or "step").lower().replace("-", "_") aliases = { @@ -110,30 +136,89 @@ def _contiguous_true_runs(mask: np.ndarray) -> list[tuple[int, int]]: return runs +def _close_small_gaps(mask: np.ndarray, gap_tol: int) -> np.ndarray: + """Bridge ``False`` gaps no longer than ``gap_tol`` between two ``True`` runs. + + Per-pixel noise punches sub-threshold holes through an otherwise coherent + bright/dark segment; closing those holes lets the segment be recovered as a + single run instead of many short fragments. + """ + if gap_tol <= 0: + return mask + out = np.asarray(mask, dtype=bool).copy() + runs = _contiguous_true_runs(out) + for (_s0, e0), (s1, _e1) in zip(runs, runs[1:]): + if s1 - e0 <= gap_tol: + out[e0:s1] = True + return out + + def _split_segments_by_adjacent_limit( segments: list[BadSegment] | tuple[BadSegment, ...], max_adjacent_bad_lines: int, ) -> tuple[list[BadSegment], list[BadSegment]]: - """Return (accepted, skipped) after applying adjacent-line safety limit.""" + """Return (accepted, skipped) after applying the adjacent-line safety limit. + + The limit guards against repairing a vertical *block* of bad pixels that has + no clean neighbour row to copy from. Crucially it is **column-aware**: two + segments on adjacent lines only belong to the same stack when their column + ranges overlap. Segments that merely share a line index with unrelated bad + lines elsewhere in the image are independently repairable and are kept. A + stack taller than ``max_adjacent_bad_lines`` consecutive overlapping lines is + skipped; everything else is repaired. + """ max_adjacent_bad_lines = max(1, int(max_adjacent_bad_lines)) - by_line: dict[int, list[BadSegment]] = {} - for seg in segments: - by_line.setdefault(int(seg.line_index), []).append(seg) + segs = list(segments) + n = len(segs) + if n == 0: + return [], [] + + parent = list(range(n)) + + def find(x: int) -> int: + root = x + while parent[root] != root: + root = parent[root] + while parent[x] != root: + parent[x], x = root, parent[x] + return root + + def union(a: int, b: int) -> None: + ra, rb = find(a), find(b) + if ra != rb: + parent[ra] = rb + + by_line: dict[int, list[int]] = {} + for i, seg in enumerate(segs): + by_line.setdefault(int(seg.line_index), []).append(i) + + def _overlaps(a: int, b: int) -> bool: + return (int(segs[a].start_col) < int(segs[b].end_col) + and int(segs[b].start_col) < int(segs[a].end_col)) + + # Link a segment to any column-overlapping segment on the next line; the + # connected components are the vertical stacks of contiguous bad pixels. + for i, seg in enumerate(segs): + for j in by_line.get(int(seg.line_index) + 1, ()): + if _overlaps(i, j): + union(i, j) + + components: dict[int, list[int]] = {} + for i in range(n): + components.setdefault(find(i), []).append(i) + accepted: list[BadSegment] = [] skipped: list[BadSegment] = [] - lines = sorted(by_line) - i = 0 - while i < len(lines): - group = [lines[i]] - i += 1 - while i < len(lines) and lines[i] == group[-1] + 1: - group.append(lines[i]) - i += 1 - group_segments = [seg for line in group for seg in by_line[line]] - if len(group) > max_adjacent_bad_lines: - skipped.extend(group_segments) - else: - accepted.extend(group_segments) + for members in components.values(): + lines = {int(segs[i].line_index) for i in members} + # Members are linked only through ±1-line overlaps, so the line set is a + # consecutive run; its span is the stack thickness. + stack_height = max(lines) - min(lines) + 1 + target = accepted if stack_height <= max_adjacent_bad_lines else skipped + target.extend(segs[i] for i in members) + + accepted.sort(key=lambda s: (int(s.line_index), int(s.start_col))) + skipped.sort(key=lambda s: (int(s.line_index), int(s.start_col))) return accepted, skipped @@ -203,27 +288,50 @@ def detect_bad_scanline_segments( k += 2 return segments - scale = _robust_scale(residuals) + # The residual is already referenced to the *neighbour rows*, so good pixels + # sit at ~0 and a bright/dark defect is a direct deviation from zero. We do + # NOT subtract each row's own median (the previous behaviour): once a defect + # covers more than ~half the row it becomes the row median, nets its own + # pixels to ~0, and vanishes — so lowering the threshold only surfaced noise. + # + # Scan-line defects are *extended*, so we smooth the residual along each row + # first (a matched filter): this preserves a coherent line's amplitude while + # averaging per-pixel noise down, letting a faint-but-long line clear the + # threshold without texture flooding in. ``step`` remains the detector for + # sharp, short scars. + window = max(1, min(_MAD_SMOOTH_WINDOW_PX, Nx)) + smoothed = _smooth_rows_nanaware(residuals, window) + scale = _robust_scale(smoothed) cutoff = threshold * scale if not np.isfinite(cutoff): return [] + + sign = 1.0 if polarity == "bright" else -1.0 segments = [] for row in range(Ny): - residual = residuals[row] - finite = np.isfinite(residual) + sm = smoothed[row] + finite = np.isfinite(sm) if not finite.any(): continue - centre = float(np.median(residual[finite])) - if polarity == "bright": - bad = finite & ((residual - centre) > cutoff) - else: - bad = finite & ((residual - centre) < -cutoff) + bad = finite & ((sign * sm) > cutoff) + # Bridge short noise-sized holes (independent of ``min_length`` so + # raising the minimum feature length to reject texture does not start + # merging unrelated noise spikes into spurious long runs). + bad = _close_small_gaps(bad, _NOISE_GAP_TOL_PX) & finite for start, end in _contiguous_true_runs(bad): - length = end - start - if not (min_length <= length <= max_len): + # No upper-length cap: a long (even full-width) coherent defect is + # exactly the target, and the neighbour-referenced baseline stays + # valid no matter how much of the row is affected. + if end - start < min_length: + continue + seg = sm[start:end] + seg = seg[np.isfinite(seg)] + if seg.size == 0: continue - local = np.abs(residual[start:end] - centre) - score = float(np.nanmedian(local) / max(scale, 1e-15)) + seg_median = float(np.median(seg)) + if (sign * seg_median) <= cutoff: + continue # not a defect on the median + score = float(abs(seg_median) / max(scale, 1e-15)) segments.append(BadSegment(row, start, end, score, method)) return segments diff --git a/tests/test_processing.py b/tests/test_processing.py index d4d02f4..eee16e1 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -94,21 +94,71 @@ def test_detects_and_repairs_only_injected_scanline_segment(self): np.testing.assert_allclose(out[7, 6:14], arr[7, 6:14], atol=1e-12) def test_mad_outlier_method_repairs_only_segment(self): - yy, xx = np.mgrid[:16, :24] - arr = (0.02 * yy + 0.01 * xx).astype(np.float64) + # The MAD method is a matched-filter detector for extended defects, so + # the segment must be longer than the smoothing window and the image + # needs realistic noise (a noise-free image collapses the robust scale). + rng = np.random.default_rng(0) + yy, xx = np.mgrid[:24, :48] + arr = (0.02 * yy + 0.01 * xx + rng.normal(0.0, 0.05, (24, 48))).astype(np.float64) damaged = arr.copy() - damaged[8, 5:11] -= 8.0 + damaged[12, 8:36] -= 8.0 # one long dark segment, rest of the row fine corrected, info = correct_bad_scanline_segments( damaged, threshold=5.0, method="mad", polarity="dark") - assert len(info.segments) == 1 - seg = info.segments[0] - assert (seg.line_index, seg.start_col, seg.end_col) == (8, 5, 11) - outside = np.ones(damaged.shape, dtype=bool) - outside[8, 5:11] = False - np.testing.assert_array_equal(corrected[outside], damaged[outside]) - np.testing.assert_allclose(corrected[8, 5:11], arr[8, 5:11], atol=1e-12) + assert {s.line_index for s in info.segments} == {12} # only the defect row + seg = max(info.segments, key=lambda s: s.end_col - s.start_col) + assert seg.start_col <= 10 and seg.end_col >= 34 # covers the defect + # the defect interior is repaired back toward the clean surface + np.testing.assert_allclose(corrected[12, 12:32], arr[12, 12:32], atol=0.3) + # rows away from the defect are untouched + np.testing.assert_array_equal(corrected[:11], damaged[:11]) + np.testing.assert_array_equal(corrected[14:], damaged[14:]) + + def test_mad_detects_whole_line_offset_not_just_internal_segments(self): + # A uniformly bright scan line has no internal step and nets to ~0 once a + # per-line median is subtracted, so the older segment-only logic missed it + # entirely. The augmented MAD method must flag the whole line. + rng = np.random.default_rng(1) + arr = (rng.normal(0.0, 0.3, (40, 60)) + + np.linspace(0.0, 2.0, 60)[None, :]).astype(np.float64) + damaged = arr.copy() + damaged[20] += 6.0 # whole-line bright offset, flat across the row + + segs = detect_bad_scanline_segments( + damaged, threshold=4.0, method="mad", polarity="bright") + by_line = {s.line_index: s for s in segs} + assert 20 in by_line, "whole bright line not detected by MAD method" + seg = by_line[20] + assert (seg.start_col, seg.end_col) == (0, 60) # full-width, not a scar + + # The step method targets internal edges and cannot see a flat offset. + step = detect_bad_scanline_segments( + damaged, threshold=4.0, method="step", polarity="bright") + assert 20 not in {s.line_index for s in step} + + def test_mad_detects_long_partial_segment_against_noise(self): + # The real artifact: a long bright segment (~200 of 256 px) while the + # rest of the row is fine, on a noisy image. Per-pixel thresholding + # shatters it into fragments and per-row-median centring hides it once it + # exceeds half the row; the detector must instead recover it as a single + # coherent run and ignore the surrounding noise. + rng = np.random.default_rng(7) + Ny, Nx = 60, 256 + arr = (rng.normal(0.0, 1.0, (Ny, Nx)) + + np.linspace(0.0, 3.0, Nx)[None, :]).astype(np.float64) + damaged = arr.copy() + damaged[30, 20:220] += 6.0 + + segs = detect_bad_scanline_segments( + damaged, threshold=4.0, method="mad", polarity="bright") + + assert {s.line_index for s in segs} == {30} # only the defect row + on_row = [s for s in segs if s.line_index == 30] + assert len(on_row) == 1 # one coherent run, not many fragments + seg = on_row[0] + assert seg.start_col <= 30 and seg.end_col >= 210 + assert (seg.end_col - seg.start_col) >= 180 def test_bright_polarity_detects_positive_not_negative_segment(self): arr = np.ones((14, 20), dtype=np.float64) @@ -180,6 +230,30 @@ def test_max_adjacent_bad_lines_skips_unsafe_group(self): assert info.corrected_segments == () np.testing.assert_array_equal(corrected, arr) + def test_adjacent_limit_is_column_aware(self): + # Two bad segments on adjacent lines in DIFFERENT columns are + # independently repairable; the adjacent-line limit must not discard + # them together just because their line indices are consecutive. + arr = np.ones((14, 40), dtype=np.float64) + arr[5, 2:14] += 10.0 + arr[6, 24:36] += 10.0 # adjacent line, non-overlapping columns + + corrected, info = correct_bad_scanline_segments( + arr, threshold=5.0, method="step", polarity="bright", + max_adjacent_bad_lines=1) + + assert {s.line_index for s in info.corrected_segments} == {5, 6} + assert info.skipped_segments == () + # A genuinely overlapping 2-high stack is still skipped at the limit. + arr2 = np.ones((14, 40), dtype=np.float64) + arr2[5, 4:20] += 10.0 + arr2[6, 4:20] += 10.0 + _c2, info2 = correct_bad_scanline_segments( + arr2, threshold=5.0, method="step", polarity="bright", + max_adjacent_bad_lines=1) + assert info2.corrected_segments == () + assert len(info2.skipped_segments) == 2 + def test_high_threshold_gives_no_correction(self): arr = np.ones((16, 16), dtype=np.float64) arr[7, 4:9] += 20.0