From fe8f1e85385577ecd8d3370089e217c03431a6c5 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 04:29:59 +0000 Subject: [PATCH 01/11] add detection-level denoise to ShowDiffraction --- CHANGELOG.md | 5 + docs/api/showdiffraction.md | 4 +- src/quantem/widget/showdiffraction.py | 52 ++++++++++- tests/test_showdiffraction_denoise.py | 126 ++++++++++++++++++++++++++ 4 files changed, 181 insertions(+), 6 deletions(-) create mode 100644 tests/test_showdiffraction_denoise.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 79fbbf37..a87f953d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ new `rcN` heading when that rc is published to TestPyPI. ## Unreleased +- ShowDiffraction detection denoise: center refinement and spot/ring detection + now run on a denoised view of the frame (`detect_denoise`, default `"auto"`: + Anscombe for sparse counting data, light Gaussian for moderate-SNR data, + identity when clean). All fits and measurements keep using the raw frame, so + positions and radii are never biased by the smoothing. - Add `ChooseLattice`, an interactive 2D selector for choosing an ordered origin, a1, and a2 and exposing their `(row, col)` coordinates and derived lattice vectors for downstream analysis. diff --git a/docs/api/showdiffraction.md b/docs/api/showdiffraction.md index e985454c..32cff6e3 100644 --- a/docs/api/showdiffraction.md +++ b/docs/api/showdiffraction.md @@ -58,11 +58,11 @@ no console error, no NaN frame). | Contrast histogram (dual-thumb slider) | `dp_vmin_pct`, `dp_vmax_pct` | Drag either thumb (mouse or touch) for a live preview; traits update once on release | | Center mode dropdown | `center_mode` | `auto` re-detects the BF disk; `manual` enables click-to-set | | Click to set center (manual) | `center_row`, `center_col` | Crosshair moves; spot d-spacings recompute | -| Detect spots | `_detect_spots_request`, `spots` | Auto-finds every isolated peak with contrast at least 10% of the strongest (`min_relative`); no count cap | +| Detect spots | `_detect_spots_request`, `spots` | Auto-finds every isolated peak with contrast at least 10% of the strongest (`min_relative`); candidates come from the `detect_denoise` view, positions are refined on raw data; no count cap | | Add / remove spot (click) | `_spot_add_request`, `_spot_remove_request`, `spots` | Marker placed/removed; d-spacing updates | | Move spot (Move + drag) | `_spot_move_request`, `spots` | Re-picks the spot at the drop position; stale hkl clears | | Spot pick dropdown | `spot_refine`, `snap_enabled`, `snap_radius` | Clicked spots are Gaussian-fitted, snapped to the local maximum, or kept exactly as clicked | -| Detect rings | `_detect_rings_request`, `rings` | Auto-finds all Debye–Scherrer rings above the profile prominence threshold; no count cap | +| Detect rings | `_detect_rings_request`, `rings` | Auto-finds all Debye–Scherrer rings above the profile prominence threshold, on the `detect_denoise` view; ring fits stay on the raw profile; no count cap | | Add / remove ring | `_ring_add_request`, `_ring_remove_request`, `rings` | Ring overlay; ring d-spacing updates | | Calibrate from spot / ring | `_calibrate_from_spot_request`, `_calibrate_from_ring_request`, `k_pixel_size` | Sets k-space pixel size from a known d | | Auto button | `_auto_request`, `analysis_status` | Runs center, rings, calibration, fit, and indexing in one pass; status reports failed steps only | diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index 8f3cd422..aa79d354 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -23,6 +23,7 @@ from quantem.widget.export import ensure_mobile_viewport from quantem.widget.utils.array import to_numpy +from quantem.widget.utils.display_filter import apply_display_filter from quantem.widget.utils.state_io import resolve_widget_version, save_state_file, unwrap_state_payload from quantem.widget.utils.ui import UiMode, resolve_ui_mode @@ -1703,6 +1704,11 @@ class ShowDiffraction(anywidget.AnyWidget): Search radius in pixels for snapping / Gaussian refinement. spot_refine : bool, default True Sub-pixel refine spots with a 2D Gaussian fit on add. + detect_denoise : {"auto", "none", "gaussian", "anscombe"}, default "auto" + Denoise applied to the frame before center refinement and spot/ring + detection. "auto" estimates the noise level and picks a mode; fits + and measurements always run on the raw data, so positions are not + biased by the smoothing. dp_scale_mode : str, default "log" Diffraction display scaling ("linear", "log", "sqrt"). ui_mode : {"interactive", "presentation", "report", "minimal"}, default "interactive" @@ -1763,6 +1769,7 @@ class ShowDiffraction(anywidget.AnyWidget): "snap_enabled", "snap_radius", "spot_refine", + "detect_denoise", "center_mode", "calibration_source", "calibration_ref_d", @@ -1868,6 +1875,11 @@ class ShowDiffraction(anywidget.AnyWidget): spot_refine = traitlets.Bool(True).tag(sync=True) + # detection preprocessing; fits and measurements always use raw data + detect_denoise = traitlets.Enum( + ("auto", "none", "gaussian", "anscombe"), default_value="auto" + ).tag(sync=True) + # Indexing zone_axis = traitlets.Unicode("").tag(sync=True) phase_match = traitlets.Unicode("").tag(sync=True) @@ -1978,6 +1990,7 @@ def __init__( snap_enabled: bool = False, snap_radius: int = 5, spot_refine: bool = True, + detect_denoise: str = "auto", dp_scale_mode: str = "log", ui_mode: UiMode = "interactive", show_title: bool | None = None, @@ -2019,6 +2032,7 @@ def __init__( self.snap_enabled = snap_enabled self.snap_radius = snap_radius self.spot_refine = spot_refine + self.detect_denoise = detect_denoise ui = resolve_ui_mode( ui_mode, defaults={ @@ -2211,7 +2225,7 @@ def refine_center(self, *, method: str = "symmetry", search_radius: float = 8.0) raise ValueError(f"unknown refine method {method!r}") picked = pick_center( - self._displayed_frame().astype(np.float64), + self._detection_frame().astype(np.float64), method=method, mask=self._analysis_mask(), guess=(self.center_row, self.center_col), @@ -2256,6 +2270,34 @@ def _get_frame(self, idx: int) -> np.ndarray: def _displayed_frame(self) -> np.ndarray: return self._get_frame(self.frame_idx) + def _detection_frame(self) -> np.ndarray: + # candidate finding only; measurements keep the raw frame + frame = self._displayed_frame() + mode = self._resolve_detect_denoise(frame) + if mode == "none": + return frame + return apply_display_filter(frame, mode=mode, sigma=2.0) + + def _resolve_detect_denoise(self, frame: np.ndarray) -> str: + mode = self.detect_denoise + if mode != "auto": + return mode + + frame = np.asarray(frame, dtype=np.float64) + positive = frame[frame > 0] + if positive.size == 0: + return "none" + + if np.array_equal(positive, np.round(positive)): + # counting data: typical pixels in the shot noise regime, not the beam + return "anscombe" if float(np.median(positive)) <= 30.0 else "none" + + noise = 1.4826 * float(np.median(np.abs(frame - ndimage.median_filter(frame, size=3)))) + if noise <= 0.0: + return "none" + signal = float(np.percentile(frame, 99.5) - np.median(frame)) + return "gaussian" if signal < 50.0 * noise else "none" + def _update_frame(self, change=None): frame = self._displayed_frame() self.dp_stats = [ @@ -2356,7 +2398,7 @@ def detect_spots( replace: bool = True, ) -> Self: """Detect Bragg spots with contrast at least ``min_relative`` of the strongest peak.""" - frame = self._displayed_frame().astype(np.float64) + frame = self._detection_frame().astype(np.float64) n_rows, n_cols = frame.shape if exclude_radius is None: exclude_radius = max(self.bf_radius, 2.0 * float(min_distance)) @@ -2424,7 +2466,7 @@ def detect_rings( ) -> Self: """Detect Debye-Scherrer rings from radial profile peaks (max_rings=None keeps all).""" try: - radii_px, intensity = self._radial_profile() + radii_px, intensity = self._radial_profile(frame=self._detection_frame()) except Exception: return self y = np.asarray(intensity, dtype=np.float64) @@ -3040,8 +3082,10 @@ def _radial_profile( center: tuple[float, float] | None = None, angular_range: tuple[float, float] | None = None, frame_idx: int | None = None, + frame: np.ndarray | None = None, ) -> tuple[np.ndarray, np.ndarray]: - frame = self._displayed_frame() if frame_idx is None else self._get_frame(frame_idx) + if frame is None: + frame = self._displayed_frame() if frame_idx is None else self._get_frame(frame_idx) return radial_profile_px( frame, center=center or (self.center_row, self.center_col), diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py new file mode 100644 index 00000000..72166411 --- /dev/null +++ b/tests/test_showdiffraction_denoise.py @@ -0,0 +1,126 @@ +"""Detection-level denoise: candidates from a denoised view, measurements from raw data.""" + +import numpy as np +import pytest +import traitlets + +from quantem.widget.showdiffraction import ShowDiffraction + + +def _spot_dp(size=128, center=(64, 64), spacing=24.0, amp=40.0, sigma=2.0): + rows = np.arange(size, dtype=np.float64)[:, None] + cols = np.arange(size, dtype=np.float64)[None, :] + + def blob(r, c, a, s): + return a * np.exp(-((rows - r) ** 2 + (cols - c) ** 2) / (2 * s * s)) + + dp = blob(center[0], center[1], 300.0, 4.0) + truth = [ + (center[0], center[1] + spacing), + (center[0], center[1] - spacing), + (center[0] + spacing, center[1]), + (center[0] - spacing, center[1]), + ] + for r, c in truth: + dp = dp + blob(r, c, amp, sigma) + return dp, truth + + +def _ring_dp(radii, size=256, amp=30.0, sigma=2.5): + cen = (size // 2, size // 2) + rows = np.arange(size, dtype=np.float64)[:, None] + cols = np.arange(size, dtype=np.float64)[None, :] + r = np.hypot(rows - cen[0], cols - cen[1]) + dp = 200.0 * np.exp(-(r**2) / (2 * 5.0**2)) + for rr in radii: + dp = dp + amp * np.exp(-((r - rr) ** 2) / (2 * sigma**2)) + return dp, cen + + +def _true_spots_found(spots, truth, tol=2.0): + return sum( + any(abs(s["row"] - r) <= tol and abs(s["col"] - c) <= tol for r, c in truth) + for s in spots + ) + + +def _nearest(widget, r, c): + return min(np.hypot(s["row"] - r, s["col"] - c) for s in widget.spots) + + +def test_noisy_spot_detection_improves_with_denoise(): + dp, truth = _spot_dp() + rng = np.random.default_rng(7) + counts = rng.poisson(dp * 0.05).astype(np.float32) + + kwargs = dict(center=(64, 64), bf_radius=10, verbose=False) + raw = ShowDiffraction(counts, detect_denoise="none", **kwargs).detect_spots(max_spots=8) + auto = ShowDiffraction(counts, detect_denoise="auto", **kwargs).detect_spots(max_spots=8) + + # sparse counts: raw detection picks up shot noise, denoised finds just the lattice + assert len(auto.spots) == 4 + assert _true_spots_found(auto.spots, truth) == 4 + assert len(raw.spots) > len(auto.spots) + + +def test_denoise_does_not_shift_measured_positions(): + dp, truth = _spot_dp(amp=80.0) + rng = np.random.default_rng(11) + counts = rng.poisson(dp * 2.0).astype(np.float32) + + kwargs = dict(center=(64, 64), bf_radius=10, snap_radius=5, verbose=False) + raw = ShowDiffraction(counts, detect_denoise="none", **kwargs).detect_spots(max_spots=8) + auto = ShowDiffraction(counts, detect_denoise="anscombe", **kwargs).detect_spots(max_spots=8) + + # both find the lattice; refined positions agree because fits run on raw data + for widget in (raw, auto): + assert _true_spots_found(widget.spots, truth) == 4 + for r, c in truth: + p_raw = _nearest(raw, r, c) + p_auto = _nearest(auto, r, c) + assert p_auto <= 1.0 + assert abs(p_auto - p_raw) <= 0.3 + + +def test_noisy_ring_detection_and_raw_fit(): + radii = [50.0, 85.0] + dp, cen = _ring_dp(radii, amp=12.0) + rng = np.random.default_rng(3) + counts = rng.poisson(dp * 0.2).astype(np.float32) + + w = ShowDiffraction(counts, center=cen, bf_radius=15, detect_denoise="auto", verbose=False) + w.detect_rings(max_rings=4) + found = sorted(ring["radius_px"] for ring in w.rings) + for target in radii: + assert any(abs(f - target) <= 3.0 for f in found) + + # fit runs on the raw profile, so radii stay honest + w.fit_ring_profile() + for target in radii: + assert any( + r.get("fit_quality") is not None and abs(r["radius_px"] - target) <= 1.5 + for r in w.rings + ) + + +def test_auto_is_noop_on_clean_data(): + dp, truth = _spot_dp() + kwargs = dict(center=(64, 64), bf_radius=10, verbose=False) + raw = ShowDiffraction(dp.astype(np.float32), detect_denoise="none", **kwargs).detect_spots() + auto = ShowDiffraction(dp.astype(np.float32), detect_denoise="auto", **kwargs).detect_spots() + + assert len(raw.spots) == len(auto.spots) + for s_raw, s_auto in zip(raw.spots, auto.spots): + assert s_raw["row"] == s_auto["row"] and s_raw["col"] == s_auto["col"] + + +def test_detect_denoise_state_and_validation(): + dp, _ = _spot_dp() + w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) + assert w.state_dict()["detect_denoise"] == "gaussian" + + restored = ShowDiffraction(dp.astype(np.float32), verbose=False, state=w.state_dict()) + assert restored.detect_denoise == "gaussian" + + with pytest.raises(traitlets.TraitError): + ShowDiffraction(dp.astype(np.float32), detect_denoise="median", verbose=False) From d01727571e7f14088496d9ab0589b44b38699a33 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 05:00:34 +0000 Subject: [PATCH 02/11] add detection denoise tutorial section --- docs/tutorials/showdiffraction.ipynb | 42 ++++++++++++++++++---------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/docs/tutorials/showdiffraction.ipynb b/docs/tutorials/showdiffraction.ipynb index 0c5aa74b..853953e0 100644 --- a/docs/tutorials/showdiffraction.ipynb +++ b/docs/tutorials/showdiffraction.ipynb @@ -27,7 +27,7 @@ "quantem showdiffraction --demo # no data needed: real Fe3O4 SAED\n", "```\n", "\n", - "The same analysis is a few lines in a notebook. `showdiffraction_fe3o4` fetches the tutorial pattern (a packaged copy is used offline): Fe3O4 nanoparticles on an amorphous carbon support. The broad support halo near d = 3.6 Å is skipped with `exclude_radius`; Auto refines the center, detects rings, calibrates from the Fe3O4 phase, fits ring profiles, and indexes hkl." + "The same analysis is a few lines in a notebook. `showdiffraction_fe3o4` fetches the tutorial pattern (a packaged copy is used offline): Fe3O4 nanoparticles on an amorphous carbon support. The broad support halo near d = 3.6 \u00c5 is skipped with `exclude_radius`; Auto refines the center, detects rings, calibrates from the Fe3O4 phase, fits ring profiles, and indexes hkl." ] }, { @@ -51,7 +51,7 @@ "version_minor": 1 }, "text/plain": [ - "ShowDiffraction(shape=(1, 512, 512), sampling=(1.0 Å, 0.004992728220367598 1/Å), frame=0/1, title='Fe3O4 SAED (real)')" + "ShowDiffraction(shape=(1, 512, 512), sampling=(1.0 \u00c5, 0.004992728220367598 1/\u00c5), frame=0/1, title='Fe3O4 SAED (real)')" ] }, "execution_count": 1, @@ -96,7 +96,7 @@ "\n", "- One-click pipeline: Auto chains center refinement, ring detection, phase calibration, profile fitting, and hkl indexing; every step is also a plain Python call.\n", "- Exact crystallography: d-spacings and interplanar angles from the full metric tensor for any crystal system, systematic absences for the common structure types, and a phase library of 100+ standards with cited lattice constants.\n", - "- Honest identification: candidates are ranked by plain facts (matched lines, mean Δd, missing strong lines), and ranking stays a verification aid; nothing is applied silently.\n", + "- Honest identification: candidates are ranked by plain facts (matched lines, mean \u0394d, missing strong lines), and ranking stays a verification aid; nothing is applied silently.\n", "- Distortion-aware: a fitted elliptical distortion corrects every radius, profile, and calibration rather than being a separate pass.\n", "- Reproducible: JSON state save/load, CSV/JSON measurement tables, and the same analysis from the notebook, the exported HTML, or the `quantem` CLI." ] @@ -190,7 +190,7 @@ "version_minor": 1 }, "text/plain": [ - "ShowDiffraction(shape=(1, 256, 256), sampling=(1.0 Å, 0.018 1/Å), frame=0/1, spots=12, title='Single-crystal SAED')" + "ShowDiffraction(shape=(1, 256, 256), sampling=(1.0 \u00c5, 0.018 1/\u00c5), frame=0/1, spots=12, title='Single-crystal SAED')" ] }, "execution_count": 3, @@ -246,7 +246,7 @@ "version_minor": 1 }, "text/plain": [ - "ShowDiffraction(shape=(1, 512, 512), sampling=(1.0 Å, 0.003988000818768712 1/Å), frame=0/1, title='Magnetite-like rings')" + "ShowDiffraction(shape=(1, 512, 512), sampling=(1.0 \u00c5, 0.003988000818768712 1/\u00c5), frame=0/1, title='Magnetite-like rings')" ] }, "execution_count": 4, @@ -294,10 +294,10 @@ "name": "stdout", "output_type": "stream", "text": [ - "Fe3O4: 5/5 lines, mean Δd 0.82%, missing strong n/a\n", - "γ-Fe2O3: 5/5 lines, mean Δd 0.87%, missing strong n/a\n", - "α-Fe2O3 (hematite): 5/5 lines, mean Δd 1.26%, missing strong n/a\n", - "α-Fe: 1/5 lines, mean Δd 2.78%, missing strong n/a\n", + "Fe3O4: 5/5 lines, mean \u0394d 0.82%, missing strong n/a\n", + "\u03b3-Fe2O3: 5/5 lines, mean \u0394d 0.87%, missing strong n/a\n", + "\u03b1-Fe2O3 (hematite): 5/5 lines, mean \u0394d 1.26%, missing strong n/a\n", + "\u03b1-Fe: 1/5 lines, mean \u0394d 2.78%, missing strong n/a\n", "Fe3O4\n" ] } @@ -305,7 +305,7 @@ "source": [ "from quantem.widget import library_phase\n", "\n", - "expected = [library_phase(n) for n in (\"Fe3O4\", \"γ-Fe2O3\", \"α-Fe2O3 (hematite)\", \"α-Fe\")]\n", + "expected = [library_phase(n) for n in (\"Fe3O4\", \"\u03b3-Fe2O3\", \"\u03b1-Fe2O3 (hematite)\", \"\u03b1-Fe\")]\n", "verified = real.identify_phase(expected)\n", "\n", "for candidate in verified:\n", @@ -314,7 +314,7 @@ " missing = candidate[\"n_missing_strong\"]\n", " print(\n", " f\"{candidate['name']}: {candidate['matched']}/{candidate['n_obs']} lines, \"\n", - " f\"mean Δd {error_text}, \"\n", + " f\"mean \u0394d {error_text}, \"\n", " f\"missing strong {'n/a' if missing is None else missing}\"\n", " )\n", "\n", @@ -323,6 +323,20 @@ "print(candidates[0][\"name\"])" ] }, + { + "cell_type": "markdown", + "id": "29b98a84", + "metadata": {}, + "source": "## Detection denoise\n\nThe real Fe3O4 SAED thinned to a median of 0.5 counts per pixel, analyzed twice side by side. Detection and Auto run on a matched-filter denoised view (`detect_denoise`, default `auto`; Anscombe on counting data) while ring fits and measurements stay on the raw counts. The right panel also sets the display-only `denoise = \"nlm\"` so the pattern itself is readable; `show_detection_view` shows what detection saw instead." + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8edfbcc8", + "metadata": {}, + "outputs": [], + "source": "from ipywidgets import HBox, Layout\n\nfrom quantem.widget.data import showdiffraction_fe3o4\n\nfe3o4_saed = showdiffraction_fe3o4(verbose=False)\nscale = 0.5 / np.median(fe3o4_saed[fe3o4_saed > 0])\nsparse_saed = np.random.default_rng(2).poisson(np.clip(fe3o4_saed, 0, None) * scale).astype(np.float32)\n\n\ndef low_dose_auto(mode, title):\n w = ShowDiffraction(sparse_saed, title=title, offline=True, verbose=False, panel_width_px=430)\n w.detect_denoise = mode\n if mode != \"none\":\n w.denoise = \"nlm\"\n w.phase_name = \"Fe3O4\"\n w.run_auto()\n indexed = sum(1 for r in w.rings if r.get(\"hkl\"))\n print(f\"{title}: calibration rms {w.calibration_rms_px:.2f} px, {indexed}/{len(w.rings)} rings indexed\")\n return w\n\n\nraw = low_dose_auto(\"none\", \"raw\")\ndenoised = low_dose_auto(\"auto\", \"denoised\")\nHBox([raw, denoised], layout=Layout(overflow=\"auto\"))" + }, { "cell_type": "markdown", "id": "sddsave01", @@ -351,9 +365,9 @@ "output_type": "stream", "text": [ "Magnetite-like rings\n", - "════════════════════════════════\n", + "\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\u2550\n", "Frames: 1 (showing #0)\n", - "Detector: 512×512 (0.0040 1/Å/px)\n", + "Detector: 512\u00d7512 (0.0040 1/\u00c5/px)\n", "Calibration: phase (rms 0.62 px)\n", "Center: (255.5, 255.5) BF r=13.1 px\n", "Spots: 0\n", @@ -406,4 +420,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file From f4d6bcfbf71ae214bfd001791c37594dfcb60ca6 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 05:20:16 +0000 Subject: [PATCH 03/11] handle gain-scaled counting data in detect denoise auto --- src/quantem/widget/showdiffraction.py | 9 ++++++--- tests/test_showdiffraction_denoise.py | 9 +++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index aa79d354..a00c8e10 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -2288,9 +2288,12 @@ def _resolve_detect_denoise(self, frame: np.ndarray) -> str: if positive.size == 0: return "none" - if np.array_equal(positive, np.round(positive)): - # counting data: typical pixels in the shot noise regime, not the beam - return "anscombe" if float(np.median(positive)) <= 30.0 else "none" + # counting data may carry a detector gain: values then sit on integer + # multiples of the smallest positive value + counts = positive / float(positive.min()) + if np.allclose(counts, np.round(counts), atol=1e-3): + # typical pixels in the shot noise regime, not the beam + return "anscombe" if float(np.median(counts)) <= 30.0 else "none" noise = 1.4826 * float(np.median(np.abs(frame - ndimage.median_filter(frame, size=3)))) if noise <= 0.0: diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index 72166411..8d69751a 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -114,6 +114,15 @@ def test_auto_is_noop_on_clean_data(): assert s_raw["row"] == s_auto["row"] and s_raw["col"] == s_auto["col"] +def test_gain_scaled_counts_use_anscombe(): + dp, _ = _spot_dp() + counts = np.random.default_rng(5).poisson(dp * 0.05).astype(np.float32) + scaled = counts * 1.55e-5 + + w = ShowDiffraction(scaled, verbose=False) + assert w._resolve_detect_denoise(scaled) == "anscombe" + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From 83bfe32fc17940684322d4eb44d77b1e5732fabd Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 06:34:51 +0000 Subject: [PATCH 04/11] add show_detection_view display toggle --- src/quantem/widget/showdiffraction.py | 10 +++++++--- tests/test_showdiffraction_denoise.py | 15 +++++++++++++++ 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index a00c8e10..94ae4073 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -1708,7 +1708,8 @@ class ShowDiffraction(anywidget.AnyWidget): Denoise applied to the frame before center refinement and spot/ring detection. "auto" estimates the noise level and picks a mode; fits and measurements always run on the raw data, so positions are not - biased by the smoothing. + biased by the smoothing. Set the ``show_detection_view`` trait to + display the denoised view instead of the raw frame. dp_scale_mode : str, default "log" Diffraction display scaling ("linear", "log", "sqrt"). ui_mode : {"interactive", "presentation", "report", "minimal"}, default "interactive" @@ -1770,6 +1771,7 @@ class ShowDiffraction(anywidget.AnyWidget): "snap_radius", "spot_refine", "detect_denoise", + "show_detection_view", "center_mode", "calibration_source", "calibration_ref_d", @@ -1879,6 +1881,8 @@ class ShowDiffraction(anywidget.AnyWidget): detect_denoise = traitlets.Enum( ("auto", "none", "gaussian", "anscombe"), default_value="auto" ).tag(sync=True) + # display the denoised detection view instead of the raw frame + show_detection_view = traitlets.Bool(False).tag(sync=True) # Indexing zone_axis = traitlets.Unicode("").tag(sync=True) @@ -2116,7 +2120,7 @@ def _set_initial_geometry( self.auto_detect_center() def _observe_traits(self) -> None: - self.observe(self._update_frame, names=["frame_idx"]) + self.observe(self._update_frame, names=["frame_idx", "detect_denoise", "show_detection_view"]) self.observe(self._bake_offline_frames, names=["offline"]) self.observe(self._on_spot_add_request, names=["_spot_add_request"]) self.observe(self._on_spot_undo_request, names=["_spot_undo_request"]) @@ -2302,7 +2306,7 @@ def _resolve_detect_denoise(self, frame: np.ndarray) -> str: return "gaussian" if signal < 50.0 * noise else "none" def _update_frame(self, change=None): - frame = self._displayed_frame() + frame = self._detection_frame() if self.show_detection_view else self._displayed_frame() self.dp_stats = [ float(frame.mean()), float(frame.min()), diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index 8d69751a..f695d2e7 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -123,6 +123,21 @@ def test_gain_scaled_counts_use_anscombe(): assert w._resolve_detect_denoise(scaled) == "anscombe" +def test_show_detection_view_ships_denoised_frame(): + dp, _ = _spot_dp() + counts = np.random.default_rng(5).poisson(dp * 0.05).astype(np.float32) + + w = ShowDiffraction(counts, detect_denoise="anscombe", verbose=False) + raw_bytes = w.frame_bytes + w.show_detection_view = True + assert w.frame_bytes != raw_bytes + + # display only: stored data and measurements stay raw + assert np.array_equal(w._displayed_frame(), counts) + w.show_detection_view = False + assert w.frame_bytes == raw_bytes + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From f156225f6980baafa141a47a57edb2c4898d095d Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 07:07:12 +0000 Subject: [PATCH 05/11] add nlm display filter --- src/quantem/widget/utils/display_filter.py | 32 ++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/quantem/widget/utils/display_filter.py b/src/quantem/widget/utils/display_filter.py index d9e3d81f..6554862b 100644 --- a/src/quantem/widget/utils/display_filter.py +++ b/src/quantem/widget/utils/display_filter.py @@ -28,6 +28,7 @@ "anscombe", "bin2_anscombe", "bin4_anscombe", + "nlm", "tv", "denova", "denova_tv", @@ -120,6 +121,33 @@ def _bin2(image: np.ndarray, sigma: float | None = None) -> np.ndarray: ).astype(np.float32) +def _nlm(image: np.ndarray) -> np.ndarray: + """Poisson non-local means: variance-stabilize, patch-average, invert.""" + try: + from skimage.restoration import denoise_nl_means + except ImportError as exc: + raise ImportError("display filter 'nlm' requires scikit-image") from exc + + image = np.asarray(image, dtype=np.float32) + positive = image[image > 0] + if positive.size == 0: + return image.copy() + + # true counts when the data is gain-quantized, pseudo counts otherwise + gain = float(positive.min()) + ratio = positive / gain + if not np.allclose(ratio, np.round(ratio), atol=1e-3): + gain = float(np.percentile(positive, 99.5)) / 30.0 + + counts = np.clip(image, 0.0, None) / gain + stabilized = 2.0 * np.sqrt(counts + 3.0 / 8.0) # noise std near 1 + smoothed = denoise_nl_means( + stabilized, patch_size=5, patch_distance=6, h=0.8, sigma=1.0, fast_mode=True + ) + inverse = np.clip((smoothed * 0.5) ** 2 - 3.0 / 8.0, 0.0, None) + return (inverse * gain).astype(np.float32) + + def _denova(image: np.ndarray, method: str = "tv") -> np.ndarray: """Lab denova denoiser (auto-lambda). Requires the denova package.""" try: @@ -162,6 +190,8 @@ def apply_display_filter( - ``"anscombe"``: Anscombe transform, Gaussian, inverse; respects Poisson statistics of count data. With ``spatial_bin`` >= 2 the smoothing width becomes ``max(2, sigma*0.75)`` on the binned map. + - ``"nlm"``: Poisson non-local means, variance-stabilized patch + averaging that keeps spots sharp (requires scikit-image). - ``"tv"``: total-variation denoise (requires scikit-image). - ``"denova"`` / ``"denova_tv"`` / ``"denova_tv12"``: lab denova denoiser when the optional package is installed. @@ -235,6 +265,8 @@ def apply_display_filter( out = _bin2(out, None) if mode == "none": return out + if mode == "nlm": + return _nlm(out) if mode == "tv": try: from skimage.restoration import denoise_tv_chambolle From bc3d9843f79537030bc652181cb191d0cefe915b Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 07:07:12 +0000 Subject: [PATCH 06/11] add display-only denoise trait to ShowDiffraction --- src/quantem/widget/showdiffraction.py | 21 ++++++++++++++++++--- tests/test_showdiffraction_denoise.py | 14 ++++++++++++++ 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index 94ae4073..dc83b0a8 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -1709,7 +1709,9 @@ class ShowDiffraction(anywidget.AnyWidget): detection. "auto" estimates the noise level and picks a mode; fits and measurements always run on the raw data, so positions are not biased by the smoothing. Set the ``show_detection_view`` trait to - display the denoised view instead of the raw frame. + display the denoised view instead of the raw frame, or the + ``denoise`` trait for a display-only filter (``"nlm"`` keeps spots + sharp) that touches neither detection nor measurements. dp_scale_mode : str, default "log" Diffraction display scaling ("linear", "log", "sqrt"). ui_mode : {"interactive", "presentation", "report", "minimal"}, default "interactive" @@ -1772,6 +1774,7 @@ class ShowDiffraction(anywidget.AnyWidget): "spot_refine", "detect_denoise", "show_detection_view", + "denoise", "center_mode", "calibration_source", "calibration_ref_d", @@ -1883,6 +1886,10 @@ class ShowDiffraction(anywidget.AnyWidget): ).tag(sync=True) # display the denoised detection view instead of the raw frame show_detection_view = traitlets.Bool(False).tag(sync=True) + # display-only denoise for the shipped frame, stored data stays raw + denoise = traitlets.Enum( + ("none", "gaussian", "anscombe", "nlm", "tv"), default_value="none" + ).tag(sync=True) # Indexing zone_axis = traitlets.Unicode("").tag(sync=True) @@ -2120,7 +2127,10 @@ def _set_initial_geometry( self.auto_detect_center() def _observe_traits(self) -> None: - self.observe(self._update_frame, names=["frame_idx", "detect_denoise", "show_detection_view"]) + self.observe( + self._update_frame, + names=["frame_idx", "detect_denoise", "show_detection_view", "denoise"], + ) self.observe(self._bake_offline_frames, names=["offline"]) self.observe(self._on_spot_add_request, names=["_spot_add_request"]) self.observe(self._on_spot_undo_request, names=["_spot_undo_request"]) @@ -2306,7 +2316,12 @@ def _resolve_detect_denoise(self, frame: np.ndarray) -> str: return "gaussian" if signal < 50.0 * noise else "none" def _update_frame(self, change=None): - frame = self._detection_frame() if self.show_detection_view else self._displayed_frame() + if self.show_detection_view: + frame = self._detection_frame() + elif self.denoise != "none": + frame = apply_display_filter(self._displayed_frame(), mode=self.denoise, sigma=2.0) + else: + frame = self._displayed_frame() self.dp_stats = [ float(frame.mean()), float(frame.min()), diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index f695d2e7..520f3287 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -138,6 +138,20 @@ def test_show_detection_view_ships_denoised_frame(): assert w.frame_bytes == raw_bytes +def test_display_denoise_is_view_only(): + dp, _ = _spot_dp() + counts = np.random.default_rng(5).poisson(dp * 0.05).astype(np.float32) + + w = ShowDiffraction(counts, verbose=False) + raw_bytes = w.frame_bytes + w.denoise = "nlm" + assert w.frame_bytes != raw_bytes + assert np.array_equal(w._displayed_frame(), counts) + + w.denoise = "none" + assert w.frame_bytes == raw_bytes + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From 8675fd8af9e68e63e81c4b1dec1576fac4c0ab88 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 07:08:22 +0000 Subject: [PATCH 07/11] changelog for display denoise --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a87f953d..f88422e6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,10 @@ new `rcN` heading when that rc is published to TestPyPI. Anscombe for sparse counting data, light Gaussian for moderate-SNR data, identity when clean). All fits and measurements keep using the raw frame, so positions and radii are never biased by the smoothing. +- ShowDiffraction display denoise: a view-only `denoise` trait (including the + new Poisson non-local means `nlm` filter, which keeps spots sharp where the + detection blur softens them) and a `show_detection_view` toggle that + displays what detection saw; both leave stored data and measurements raw. - Add `ChooseLattice`, an interactive 2D selector for choosing an ordered origin, a1, and a2 and exposing their `(row, col)` coordinates and derived lattice vectors for downstream analysis. From d30552f42ceaa1f24ed7ffa7e34632ff7316bc8d Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 07:49:34 +0000 Subject: [PATCH 08/11] expose detect_spots noise floor as noise_sigma --- CHANGELOG.md | 3 +++ src/quantem/widget/showdiffraction.py | 10 ++++++++-- tests/test_showdiffraction_denoise.py | 21 +++++++++++++++++++++ 3 files changed, 32 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f88422e6..7c9e666a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,9 @@ new `rcN` heading when that rc is published to TestPyPI. Anscombe for sparse counting data, light Gaussian for moderate-SNR data, identity when clean). All fits and measurements keep using the raw frame, so positions and radii are never biased by the smoothing. +- ShowDiffraction `detect_spots` exposes its shot-noise contrast floor as + `noise_sigma`; lower it on frames whose diffuse scattering or detector + shadows inflate the robust noise estimate past real peak contrast. - ShowDiffraction display denoise: a view-only `denoise` trait (including the new Poisson non-local means `nlm` filter, which keeps spots sharp where the detection blur softens them) and a `show_detection_view` toggle that diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index dc83b0a8..773c215c 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -2417,9 +2417,15 @@ def detect_spots( min_distance: int = 6, min_relative: float = 0.1, exclude_radius: float | None = None, + noise_sigma: float = 5.0, replace: bool = True, ) -> Self: - """Detect Bragg spots with contrast at least ``min_relative`` of the strongest peak.""" + """Detect Bragg spots with contrast at least ``min_relative`` of the strongest peak. + + ``noise_sigma`` sets the shot-noise contrast floor in robust sigma + units; lower it on frames whose background structure inflates the + estimate (diffuse scattering, detector shadows). + """ frame = self._detection_frame().astype(np.float64) n_rows, n_cols = frame.shape if exclude_radius is None: @@ -2451,7 +2457,7 @@ def detect_spots( # contrast relative to the strongest peak, with a noise floor on noisy data contrast = np.expm1(prominence) sigma = 1.4826 * float(np.median(np.abs(work - np.median(work)))) - level = max(min_relative * float(contrast.max()), float(np.expm1(5.0 * sigma))) + level = max(min_relative * float(contrast.max()), float(np.expm1(noise_sigma * sigma))) keep = (prominence > 0) & (contrast >= level) coords, prominence = coords[keep], prominence[keep] if coords.size: diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index 520f3287..ff1d3905 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -152,6 +152,27 @@ def test_display_denoise_is_view_only(): assert w.frame_bytes == raw_bytes +def test_noise_sigma_floor_is_tunable(): + rows = np.arange(160, dtype=np.float64)[:, None] + cols = np.arange(160, dtype=np.float64)[None, :] + r = np.hypot(rows - 80, cols - 80) + + # wide dim spots on a diffuse cloud, the structured background case + dp = 80.0 * np.exp(-(r**2) / (2 * 20.0**2)) + truth = [(80, 125), (35, 80), (125, 80), (80, 35)] + for rr, cc in truth: + dp = dp + 6.0 * np.exp(-((rows - rr) ** 2 + (cols - cc) ** 2) / (2 * 4.0**2)) + counts = np.random.default_rng(9).poisson(dp).astype(np.float32) + + w = ShowDiffraction(counts, center=(80, 80), bf_radius=10, detect_denoise="anscombe", verbose=False) + w.detect_spots(min_distance=10, min_relative=0.3) + assert len(w.spots) == 0 + + w.detect_spots(min_distance=10, min_relative=0.3, noise_sigma=3.0) + near = sum(any(abs(s["row"] - a) <= 3 and abs(s["col"] - b) <= 3 for a, b in truth) for s in w.spots) + assert near >= 2 + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From ae13daa2d6d1dc5956dd344d4443036c7d76b862 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 07:56:26 +0000 Subject: [PATCH 09/11] test detected Au ring ratios index as fcc --- tests/test_showdiffraction_denoise.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index ff1d3905..6db94fce 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -173,6 +173,28 @@ def test_noise_sigma_floor_is_tunable(): assert near >= 2 +def test_detected_au_ratios_index_as_fcc(): + from quantem.widget import library_phase + + au = library_phase("Au") + k = 0.008 + hkls = ("111", "200", "220", "311") + radii = [1.0 / (au.d_spacing(tuple(int(c) for c in h)) * k) for h in hkls] + dp, cen = _ring_dp(radii, amp=25.0, sigma=2.0) + counts = np.random.default_rng(6).poisson(dp * 0.5).astype(np.float32) + + w = ShowDiffraction(counts, center=cen, bf_radius=15, verbose=False) + w.detect_rings(max_rings=4, exclude_radius=30) + w.calibrate_from_phase(au) + w.index_rings(au) + + # detected radius ratios must land on the allowed fcc reflections + ordered = sorted(w.rings, key=lambda r: r["radius_px"]) + assert [r["hkl"] for r in ordered] == list(hkls) + assert all(r["d_error"] < 0.01 for r in ordered) + assert abs(w.k_pixel_size - k) / k < 0.01 + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From ea7d49614260f076fd06c3e503accad3d98ae942 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 08:02:48 +0000 Subject: [PATCH 10/11] recover phase-predicted rings in run_auto --- src/quantem/widget/showdiffraction.py | 54 +++++++++++++++++++++++++++ tests/test_showdiffraction_denoise.py | 30 +++++++++++++++ 2 files changed, 84 insertions(+) diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index 773c215c..9233a370 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -2946,6 +2946,54 @@ def _all_phases(self, custom_only: bool = False) -> list[Phase]: continue return phases + def recover_predicted_rings( + self, phase: Phase, *, tol_px: float = 4.0, snr: float = 3.0 + ) -> list[float]: + """Add rings the detector missed at radii the calibrated phase predicts. + + For each allowed reflection with no ring within ``tol_px``, the raw + radial profile is tested near the predicted radius; a local maximum + at least ``snr`` robust sigmas above the detrended profile noise + becomes a ring. Returns the added radii in pixels. + """ + if not self.k_calibrated: + raise ValueError("recover_predicted_rings needs a k calibration") + from scipy.ndimage import gaussian_filter1d + + radii_px, intensity = self._radial_profile() + y_log = np.log1p(np.clip(intensity - intensity.min(), 0.0, None)) + detrended = y_log - gaussian_filter1d(y_log, sigma=max(3.0, y_log.size / 20.0)) + # successive differences ignore ring structure and detrending sidelobes + steps = np.diff(detrended) + noise = 1.4826 * float(np.median(np.abs(steps))) / np.sqrt(2.0) + if noise <= 0.0: + return [] + + existing = [r["radius_px"] for r in self.rings] + added = [] + for ref in phase.reflections(): + r_pred = 1.0 / (ref["d"] * self.k_pixel_size) + if r_pred <= self.bf_radius or r_pred >= float(radii_px[-1]): + continue + if any(abs(r_pred - r) <= tol_px for r in existing): + continue + window = np.abs(radii_px - r_pred) <= tol_px + if not window.any(): + continue + idx = np.flatnonzero(window) + peak = idx[np.argmax(detrended[idx])] + # interior local maximum with real prominence, not a window edge + if peak in (0, len(detrended) - 1) or peak in (idx[0], idx[-1]): + continue + if detrended[peak] < snr * noise: + continue + if detrended[peak] <= detrended[peak - 1] or detrended[peak] <= detrended[peak + 1]: + continue + self.add_ring(float(radii_px[peak])) + existing.append(float(radii_px[peak])) + added.append(float(radii_px[peak])) + return added + def run_auto( self, phase: Phase | None = None, @@ -2979,6 +3027,12 @@ def run_auto( problems.append(f"calibration failed ({exc})") phase = None if phase is not None and self.k_calibrated: + # rescue predicted reflections the detector under-detected + if self.recover_predicted_rings(phase): + try: + self.fit_ring_profile() + except (ValueError, ImportError): + pass try: self.index_rings(phase) if not any(r.get("hkl") for r in self.rings): diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index 6db94fce..82a1d03d 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -195,6 +195,36 @@ def test_detected_au_ratios_index_as_fcc(): assert abs(w.k_pixel_size - k) / k < 0.01 +def test_recover_predicted_rings_rescues_missed_reflection(): + from quantem.widget import library_phase + + au = library_phase("Au") + k = 0.008 + hkls = ("111", "200", "220", "311") + radii = [1.0 / (au.d_spacing(tuple(int(c) for c in h)) * k) for h in hkls] + + size = 256 + cen = (size // 2, size // 2) + rows = np.arange(size, dtype=np.float64)[:, None] + cols = np.arange(size, dtype=np.float64)[None, :] + r = np.hypot(rows - cen[0], cols - cen[1]) + dp = 200.0 * np.exp(-(r**2) / (2 * 5.0**2)) + for rr, amp in zip(radii, (25.0, 25.0, 4.0, 25.0)): + dp = dp + amp * np.exp(-((r - rr) ** 2) / (2 * 2.0**2)) + counts = np.random.default_rng(6).poisson(dp * 0.5).astype(np.float32) + + w = ShowDiffraction(counts, center=cen, bf_radius=15, verbose=False) + w.detect_rings(max_rings=3, exclude_radius=30) + assert len(w.rings) == 3 + + w.calibrate_from_phase(au) + added = w.recover_predicted_rings(au) + assert len(added) == 1 and abs(added[0] - radii[1]) <= 2.0 + + w.index_rings(au) + assert [x["hkl"] for x in sorted(w.rings, key=lambda y: y["radius_px"])] == list(hkls) + + def test_detect_denoise_state_and_validation(): dp, _ = _spot_dp() w = ShowDiffraction(dp.astype(np.float32), detect_denoise="gaussian", verbose=False) From b8ced224b79a86b2b51dc0fb63326d0ca0016778 Mon Sep 17 00:00:00 2001 From: henryhng Date: Mon, 10 Aug 2026 17:10:12 +0000 Subject: [PATCH 11/11] tidy denoise style --- src/quantem/widget/showdiffraction.py | 23 ++++++++++------------ src/quantem/widget/utils/display_filter.py | 4 ++-- tests/test_showdiffraction_denoise.py | 15 +++----------- 3 files changed, 15 insertions(+), 27 deletions(-) diff --git a/src/quantem/widget/showdiffraction.py b/src/quantem/widget/showdiffraction.py index 9233a370..7777a55c 100644 --- a/src/quantem/widget/showdiffraction.py +++ b/src/quantem/widget/showdiffraction.py @@ -1880,13 +1880,11 @@ class ShowDiffraction(anywidget.AnyWidget): spot_refine = traitlets.Bool(True).tag(sync=True) - # detection preprocessing; fits and measurements always use raw data + # Denoise detect_denoise = traitlets.Enum( ("auto", "none", "gaussian", "anscombe"), default_value="auto" ).tag(sync=True) - # display the denoised detection view instead of the raw frame show_detection_view = traitlets.Bool(False).tag(sync=True) - # display-only denoise for the shipped frame, stored data stays raw denoise = traitlets.Enum( ("none", "gaussian", "anscombe", "nlm", "tv"), default_value="none" ).tag(sync=True) @@ -2285,7 +2283,6 @@ def _displayed_frame(self) -> np.ndarray: return self._get_frame(self.frame_idx) def _detection_frame(self) -> np.ndarray: - # candidate finding only; measurements keep the raw frame frame = self._displayed_frame() mode = self._resolve_detect_denoise(frame) if mode == "none": @@ -2302,16 +2299,17 @@ def _resolve_detect_denoise(self, frame: np.ndarray) -> str: if positive.size == 0: return "none" - # counting data may carry a detector gain: values then sit on integer - # multiples of the smallest positive value + # Counting data counts = positive / float(positive.min()) if np.allclose(counts, np.round(counts), atol=1e-3): - # typical pixels in the shot noise regime, not the beam return "anscombe" if float(np.median(counts)) <= 30.0 else "none" - noise = 1.4826 * float(np.median(np.abs(frame - ndimage.median_filter(frame, size=3)))) + # Continuous data + residual = frame - ndimage.median_filter(frame, size=3) + noise = 1.4826 * float(np.median(np.abs(residual))) if noise <= 0.0: return "none" + signal = float(np.percentile(frame, 99.5) - np.median(frame)) return "gaussian" if signal < 50.0 * noise else "none" @@ -2958,12 +2956,11 @@ def recover_predicted_rings( """ if not self.k_calibrated: raise ValueError("recover_predicted_rings needs a k calibration") - from scipy.ndimage import gaussian_filter1d radii_px, intensity = self._radial_profile() y_log = np.log1p(np.clip(intensity - intensity.min(), 0.0, None)) - detrended = y_log - gaussian_filter1d(y_log, sigma=max(3.0, y_log.size / 20.0)) - # successive differences ignore ring structure and detrending sidelobes + detrended = y_log - ndimage.gaussian_filter1d(y_log, sigma=max(3.0, y_log.size / 20.0)) + # Profile noise estimate steps = np.diff(detrended) noise = 1.4826 * float(np.median(np.abs(steps))) / np.sqrt(2.0) if noise <= 0.0: @@ -2982,7 +2979,7 @@ def recover_predicted_rings( continue idx = np.flatnonzero(window) peak = idx[np.argmax(detrended[idx])] - # interior local maximum with real prominence, not a window edge + # Peak checks if peak in (0, len(detrended) - 1) or peak in (idx[0], idx[-1]): continue if detrended[peak] < snr * noise: @@ -3027,7 +3024,7 @@ def run_auto( problems.append(f"calibration failed ({exc})") phase = None if phase is not None and self.k_calibrated: - # rescue predicted reflections the detector under-detected + # Ring recovery if self.recover_predicted_rings(phase): try: self.fit_ring_profile() diff --git a/src/quantem/widget/utils/display_filter.py b/src/quantem/widget/utils/display_filter.py index 6554862b..a5bedf64 100644 --- a/src/quantem/widget/utils/display_filter.py +++ b/src/quantem/widget/utils/display_filter.py @@ -133,14 +133,14 @@ def _nlm(image: np.ndarray) -> np.ndarray: if positive.size == 0: return image.copy() - # true counts when the data is gain-quantized, pseudo counts otherwise + # Gain estimate gain = float(positive.min()) ratio = positive / gain if not np.allclose(ratio, np.round(ratio), atol=1e-3): gain = float(np.percentile(positive, 99.5)) / 30.0 counts = np.clip(image, 0.0, None) / gain - stabilized = 2.0 * np.sqrt(counts + 3.0 / 8.0) # noise std near 1 + stabilized = 2.0 * np.sqrt(counts + 3.0 / 8.0) smoothed = denoise_nl_means( stabilized, patch_size=5, patch_distance=6, h=0.8, sigma=1.0, fast_mode=True ) diff --git a/tests/test_showdiffraction_denoise.py b/tests/test_showdiffraction_denoise.py index 82a1d03d..127d6edc 100644 --- a/tests/test_showdiffraction_denoise.py +++ b/tests/test_showdiffraction_denoise.py @@ -4,6 +4,7 @@ import pytest import traitlets +from quantem.widget import library_phase from quantem.widget.showdiffraction import ShowDiffraction @@ -57,7 +58,6 @@ def test_noisy_spot_detection_improves_with_denoise(): raw = ShowDiffraction(counts, detect_denoise="none", **kwargs).detect_spots(max_spots=8) auto = ShowDiffraction(counts, detect_denoise="auto", **kwargs).detect_spots(max_spots=8) - # sparse counts: raw detection picks up shot noise, denoised finds just the lattice assert len(auto.spots) == 4 assert _true_spots_found(auto.spots, truth) == 4 assert len(raw.spots) > len(auto.spots) @@ -72,7 +72,6 @@ def test_denoise_does_not_shift_measured_positions(): raw = ShowDiffraction(counts, detect_denoise="none", **kwargs).detect_spots(max_spots=8) auto = ShowDiffraction(counts, detect_denoise="anscombe", **kwargs).detect_spots(max_spots=8) - # both find the lattice; refined positions agree because fits run on raw data for widget in (raw, auto): assert _true_spots_found(widget.spots, truth) == 4 for r, c in truth: @@ -94,7 +93,6 @@ def test_noisy_ring_detection_and_raw_fit(): for target in radii: assert any(abs(f - target) <= 3.0 for f in found) - # fit runs on the raw profile, so radii stay honest w.fit_ring_profile() for target in radii: assert any( @@ -132,7 +130,6 @@ def test_show_detection_view_ships_denoised_frame(): w.show_detection_view = True assert w.frame_bytes != raw_bytes - # display only: stored data and measurements stay raw assert np.array_equal(w._displayed_frame(), counts) w.show_detection_view = False assert w.frame_bytes == raw_bytes @@ -157,7 +154,7 @@ def test_noise_sigma_floor_is_tunable(): cols = np.arange(160, dtype=np.float64)[None, :] r = np.hypot(rows - 80, cols - 80) - # wide dim spots on a diffuse cloud, the structured background case + # Wide dim spots on a diffuse cloud dp = 80.0 * np.exp(-(r**2) / (2 * 20.0**2)) truth = [(80, 125), (35, 80), (125, 80), (80, 35)] for rr, cc in truth: @@ -169,13 +166,10 @@ def test_noise_sigma_floor_is_tunable(): assert len(w.spots) == 0 w.detect_spots(min_distance=10, min_relative=0.3, noise_sigma=3.0) - near = sum(any(abs(s["row"] - a) <= 3 and abs(s["col"] - b) <= 3 for a, b in truth) for s in w.spots) - assert near >= 2 + assert _true_spots_found(w.spots, truth, tol=3.0) >= 2 def test_detected_au_ratios_index_as_fcc(): - from quantem.widget import library_phase - au = library_phase("Au") k = 0.008 hkls = ("111", "200", "220", "311") @@ -188,7 +182,6 @@ def test_detected_au_ratios_index_as_fcc(): w.calibrate_from_phase(au) w.index_rings(au) - # detected radius ratios must land on the allowed fcc reflections ordered = sorted(w.rings, key=lambda r: r["radius_px"]) assert [r["hkl"] for r in ordered] == list(hkls) assert all(r["d_error"] < 0.01 for r in ordered) @@ -196,8 +189,6 @@ def test_detected_au_ratios_index_as_fcc(): def test_recover_predicted_rings_rescues_missed_reflection(): - from quantem.widget import library_phase - au = library_phase("Au") k = 0.008 hkls = ("111", "200", "220", "311")