From 5dc7801c117281f6f141ba2810c3bac1f38d9730 Mon Sep 17 00:00:00 2001 From: Craig Myles Date: Mon, 3 Aug 2026 17:37:52 +0100 Subject: [PATCH 1/2] Add canny_density segmenter: per-patch Canny edge density The canny segmenter dilates, closes and fills contours, which recovers faint contiguous tissue but also fills the gaps between fragments -- on sparse or fragmented slides that tiles background and gives blobby boundaries. canny_density keeps the same low-threshold (0.05) Canny as a raw edge map at a finer seg_mpp (2.0) and does no morphology at all. GridPatcher._cell_tissue_fractions averages the mask over each patch footprint, so with an edge map that mean IS the patch's edge density and --tissue-threshold becomes the minimum density to keep a patch. Each patch is therefore judged independently, on its own pixels, with no pipeline change. Threshold is a density, not a coverage fraction: start near 0.02 and re-tune if seg_mpp changes. --- docs/SEGMENTATION.md | 23 ++++ pyproject.toml | 1 + src/raw2features/segmenters/canny_density.py | 115 +++++++++++++++++++ tests/test_segmenter_canny_density.py | 80 +++++++++++++ 4 files changed, 219 insertions(+) create mode 100644 src/raw2features/segmenters/canny_density.py create mode 100644 tests/test_segmenter_canny_density.py diff --git a/docs/SEGMENTATION.md b/docs/SEGMENTATION.md index 1fbc614..e286e1f 100644 --- a/docs/SEGMENTATION.md +++ b/docs/SEGMENTATION.md @@ -12,6 +12,7 @@ from standard techniques - no GPL/CLAM code, no model weights. |---|---|---| | `otsu` *(default)* | Otsu threshold on HSV **saturation** + morphology | general default - clean & robust | | `canny` | **low**-threshold Canny (`low=0.05`) → close → **hole-aware contour fill** (`max_hole_frac` keeps large cavities as background) | good on **cervical / faint** tissue - recovers pale tissue Otsu can under-segment | +| `canny_density` | the same low-threshold Canny kept as a **raw edge map** (no dilate/close/fill) at a finer `seg_mpp` (2.0) → the patcher's per-cell mean *is* the patch's **edge density** | **per-patch** decisions with sharp, gap-respecting boundaries - fragmented / sparse slides where filling would swallow the background between fragments. **`--tissue-threshold` becomes a density** (try 0.01-0.05), not a coverage fraction | | `combined` | fuse `otsu` and `canny` - `or` (default) recovers faint tissue, `and` suppresses smooth high-saturation artefacts | when a cohort needs faint-tissue recovery (`or`) or artefact suppression (`and`) | ## Choosing a segmenter (evidence) @@ -37,6 +38,28 @@ which matches practitioners' experience on cervical H&E. So: - **`combined or`** is a safe superset of `otsu` (never loses tissue) that folds in canny's faint-tissue recovery; **`combined and`** is stricter (artefact suppression). +### `canny` (contour fill) vs `canny_density` (per-patch) + +The finding above - that a *high/auto* threshold or **edge-density** thresholding +under-segments faint tissue - was measured against a **region**-style mask: one coarse +whole-slide mask that then has to outline whole tissue regions. Filling contours is the +right fix *there*. + +`canny_density` is a different use of the same edges: it keeps the raw edge map and lets +the patcher average it per patch, so **each patch is judged on its own edge content** and +nothing is grown or filled. That combination - a **low** threshold (0.05) with a +**per-patch** density - is the recipe practitioners report working well on cervical H&E, +and it is not the configuration the 0.3% result above came from (that was a high/auto +threshold). Pick by failure mode: + +- tissue is *faint but contiguous*, and you want solid regions → **`canny`** +- tissue is *fragmented/sparse*, and filling the gaps between fragments would tile + background → **`canny_density`** + +Because the mask is a thin edge map, its per-patch mean is small: start around +`--tissue-threshold 0.02` and tune per cohort (and re-tune if you change `seg_mpp`, +since finer masks carry more, thinner edges per field of view). + ## Provenance of the built-in segmenters The built-in segmenters are independent implementations of standard, textbook techniques diff --git a/pyproject.toml b/pyproject.toml index 46a620e..3157651 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -59,6 +59,7 @@ omezarr = "raw2features.readers.omezarr:OmeZarrReader" [project.entry-points."raw2features.segmenters"] otsu = "raw2features.segmenters.otsu:OtsuSegmenter" canny = "raw2features.segmenters.canny:CannySegmenter" +canny_density = "raw2features.segmenters.canny_density:CannyDensitySegmenter" combined = "raw2features.segmenters.combined:CombinedSegmenter" nuclear = "raw2features.segmenters.nuclear:NuclearSegmenter" grandqc = "raw2features.segmenters.grandqc:GrandQCSegmenter" diff --git a/src/raw2features/segmenters/canny_density.py b/src/raw2features/segmenters/canny_density.py new file mode 100644 index 0000000..13384e2 --- /dev/null +++ b/src/raw2features/segmenters/canny_density.py @@ -0,0 +1,115 @@ +"""Per-patch Canny **edge-density** tissue segmenter (classical, permissive). + +The sibling ``canny`` segmenter takes a *region* view: low-threshold Canny -> +dilate -> close -> **fill external contours** into solid tissue regions. Filling is +what recovers faint tissue, but it also fills the gaps *between* fragments, so +sparse/fragmented slides gain background and boundaries come out blobby. + +This segmenter takes the *per-patch* view instead, the recipe practitioners use on +cervical H&E: run Canny at a **low** threshold (``low=0.05``, i.e. a raw Canny low of +~13) and keep the **raw edge map** -- no dilation, no closing, no contour fill. The +patcher averages the mask over each patch footprint +(:meth:`~raw2features.patcher.grid.GridPatcher._cell_tissue_fractions` takes +``window.mean()``), so with this mask that average *is* the patch's **edge density** +and ``--tissue-threshold`` *is* the minimum edge density to keep a patch. No pipeline +change is needed: the decision is made independently for every patch. + +Two consequences worth knowing: + +* **Sharp, gap-respecting boundaries.** Nothing grows or fills, so background + between tissue fragments stays background. +* **The threshold is a density, not a coverage fraction.** Edge *pixels* are a small + minority of even solid tissue, so useful values are far lower than the 0.1 default + that suits a filled/binary mask -- see ``docs/SEGMENTATION.md``. Density also scales + with ``seg_mpp`` (finer masks put more, thinner edges in the same field of view), so + re-tune ``--tissue-threshold`` if you change it. + +``seg_mpp`` defaults to 2.0 rather than the 8.0 used by ``otsu``/``canny``: the mask is +the resolution at which the keep/drop decision is made, and 2.0 µm/px gives a 224 px +patch at 0.5 µm/px roughly 56x56 mask pixels to average over. Pure OpenCV + numpy - no +model weights, no GPL/CLAM code. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +from raw2features.core.geometry import Point, Region, Size +from raw2features.core.mpp import nearest_level +from raw2features.core.plugins import register + +from .base import Segmenter, TissueMask + +if TYPE_CHECKING: + from raw2features.readers.base import WSISource + + +@register("segmenters", "canny_density") +class CannyDensitySegmenter(Segmenter): + """Low-threshold Canny kept as a raw edge map -> per-patch edge density. + + Parameters + ---------- + seg_mpp: + Target microns/px for the level the edge map is computed at. Finer than the + classical segmenters (default 2.0) because this mask *is* the per-patch + decision surface rather than a coarse region outline. + blur: + Gaussian-blur kernel (odd) applied before Canny, to denoise without + losing faint edges. + low, high: + Canny hysteresis thresholds as **fractions of 0-255**. ``low=0.05`` is the + low-sensitivity setting that works well on cervical H&E. + dilate: + Optional edge thickening (0 = off, the default). Leave at 0 for a faithful + density; a small value (1-2) stabilises the density on very thin edges at + fine ``seg_mpp``, at the cost of slightly inflating it. + """ + + name = "canny_density" + + def __init__( + self, + seg_mpp: float = 2.0, + blur: int = 3, + low: float = 0.05, + high: float = 0.15, + dilate: int = 0, + ) -> None: + self.seg_mpp = seg_mpp + self.blur = blur if blur % 2 == 1 else blur + 1 + self.low = low + self.high = high + self.dilate = dilate + + def _pick_level(self, reader: WSISource) -> int: + return nearest_level(reader.mpp, reader.level_downsamples(), self.seg_mpp) + + def segment(self, reader: WSISource) -> TissueMask: + import cv2 + + level = self._pick_level(reader) + dim: Size = reader.level_dimensions[level] + img = reader.read_region( + Region(level=level, location=Point(0, 0), size=Size(dim.width, dim.height)) + ) + gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) + if self.blur > 1: + gray = cv2.GaussianBlur(gray, (self.blur, self.blur), 0) + + lo = int(round(self.low * 255)) + hi = int(round(self.high * 255)) + edges = cv2.Canny(gray, lo, hi) + + # Deliberately NO close / no contour fill: the raw edge map is the signal. + # Optional thickening only, off by default. + if self.dilate > 0: + edges = cv2.dilate( + edges, np.ones((self.dilate, self.dilate), np.uint8), iterations=1 + ) + + mask = (edges > 0).astype(np.float32) # patcher averages this -> edge density + ds = float(reader.level_downsamples()[level]) + return TissueMask(mask=mask, level=level, downsample=ds) diff --git a/tests/test_segmenter_canny_density.py b/tests/test_segmenter_canny_density.py new file mode 100644 index 0000000..39e3d01 --- /dev/null +++ b/tests/test_segmenter_canny_density.py @@ -0,0 +1,80 @@ +"""Tests for the per-patch Canny edge-density segmenter.""" + +from __future__ import annotations + +import numpy as np +import pytest +import zarr + +from raw2features.readers.omezarr import OmeZarrReader +from raw2features.segmenters.canny import CannySegmenter +from raw2features.segmenters.canny_density import CannyDensitySegmenter + +from conftest import build_ngff_v04 + + +@pytest.fixture +def edged_ngff(tmp_path) -> str: + """A synthetic store containing a hard-edged block. + + The shared ``synthetic_ngff`` fixture is a smooth ramp, which by construction has + no Canny edges -- fine for shape/dtype assertions, useless for edge behaviour. This + paints a solid dark square on a light field so there is a real contour to find. + """ + path = build_ngff_v04(str(tmp_path / "edged.zarr")) + g = zarr.open_group(path, mode="r+", zarr_format=2) + for i, (h, w) in enumerate(((200, 300), (100, 150), (50, 75))): + a = g[str(i)] + for c in range(3): + plane = np.full((h, w), 230, dtype="uint8") + plane[h // 4 : 3 * h // 4, w // 4 : 3 * w // 4] = 30 # hard-edged block + a[0, c, 0] = plane + return path + + +def test_canny_density_returns_binary_edge_map_at_level(synthetic_ngff): + with OmeZarrReader(synthetic_ngff) as r: + # mpps are [0.5, 1.0, 2.0]; seg_mpp 2.0 -> level 2. + tm = CannyDensitySegmenter(seg_mpp=2.0).segment(r) + assert tm.level == 2 + assert tm.mask.ndim == 2 + assert tm.mask.dtype == np.float32 + assert set(np.unique(tm.mask)).issubset({0.0, 1.0}) + dim = r.level_dimensions[tm.level] + assert tm.mask.shape == (dim.height, dim.width) + assert tm.downsample == 4.0 + + +def test_canny_density_finds_edges_but_does_not_fill_them(edged_ngff): + """The whole point: keep the thin edge map, never fill it into a solid region.""" + with OmeZarrReader(edged_ngff) as r: + density = CannyDensitySegmenter(seg_mpp=2.0).segment(r).mask + filled = CannySegmenter(seg_mpp=2.0).segment(r).mask + + assert density.mean() > 0.0, "should detect the block's edges" + # The block covers ~1/4 of the field; a *filled* mask approaches that, while an + # edge map stays a thin outline. This is the property that stops background + # between tissue fragments from being tiled. + assert density.mean() < filled.mean() + assert density.mean() < 0.15 + + +def test_canny_density_dilate_thickens_edges(edged_ngff): + with OmeZarrReader(edged_ngff) as r: + thin = CannyDensitySegmenter(seg_mpp=2.0, dilate=0).segment(r).mask + thick = CannyDensitySegmenter(seg_mpp=2.0, dilate=3).segment(r).mask + assert thick.mean() > thin.mean() + + +def test_canny_density_lower_threshold_is_more_sensitive(edged_ngff): + with OmeZarrReader(edged_ngff) as r: + low = CannyDensitySegmenter(seg_mpp=2.0, low=0.02, high=0.06).segment(r).mask + high = CannyDensitySegmenter(seg_mpp=2.0, low=0.5, high=0.9).segment(r).mask + assert low.mean() >= high.mean() + + +def test_canny_density_is_registered_and_no_arg_constructible(): + from raw2features.core import plugins + + seg = plugins.get("segmenters", "canny_density")() # pipeline builds with no args + assert seg.name == "canny_density" From fefdabf739f2047d4288da1beb77ac0ec4f97ccc Mon Sep 17 00:00:00 2001 From: Craig Myles Date: Mon, 3 Aug 2026 23:29:39 +0100 Subject: [PATCH 2/2] Add guarded + grandqc_veto segmenters for dirty H&E cohorts Three Otsu failures on real cervical WSIs share one cause: the scanner canvas (dark no-signal fill outside the scanned region) and the scan-region border. * Otsu thresholds saturation GLOBALLY, so a dark saturated canvas corner drags the threshold up and real tissue falls below it. On the worst slide Otsu put all 344 of its patches on a corner artefact and none on the biopsy. * Canny is worse in a subtler way: the scan-region rectangle becomes the outer contour, so the tissue inside is then a max_hole_frac hole that gets carved back to background -- the fill does not add junk, it DELETES tissue. guarded = (otsu | canny) & ~scanner_canvas. Removing the canvas from the RAW EDGE MAP (before dilate/close/fill, via a new opt-in edge_filter hook on CannySegmenter) fixes both. Pure OpenCV, no weights, no GPU. grandqc_veto = guarded, then vetoed component-wise by GrandQC's tissue stage at its trained 10 um/px (never finer: off-domain runs are non-monotonic). Measured on 101 cervical WSIs (tiles vs otsu, mpp 0.5 / 224 px / thr 0.1): category n otsu guarded veto grd/otsu veto/otsu control 40 29835 39384 36936 1.32 1.24 artifact 36 71953 82757 78405 1.15 1.09 oof 18 19051 24412 22870 1.28 1.20 pale 5 18552 43195 40623 2.33 2.19 Artefact slides move their patches off the artefact and onto the biopsy (one slide 2 -> 1172 patches, another 0 -> 1097); pale slides gain ~2.3x. The veto buys a further 8 points of control inflation (1.32x -> 1.24x) for a GPU pass and CC-BY-NC-SA weights, so guarded is the better default and the veto is opt-in. Neither fixes out-of-focus -- those slides just gain patches. Also: canny_density (raw edge map -> the patcher's per-cell mean IS the patch's edge density, so --tissue-threshold becomes a density), a tissue_mpp override on GrandQCSegmenter, and a module-level model cache in qc/grandqc.py (the pipeline builds a fresh segmenter per slide, so without it the UNet++ reloads every slide). --- docs/SEGMENTATION.md | 42 ++ pyproject.toml | 2 + src/raw2features/qc/grandqc.py | 28 +- src/raw2features/segmenters/canny.py | 21 + src/raw2features/segmenters/grandqc.py | 12 +- src/raw2features/segmenters/grandqc_veto.py | 562 ++++++++++++++++++++ src/raw2features/segmenters/guarded.py | 123 +++++ tests/test_segmenter_canny_density.py | 3 +- tests/test_segmenter_grandqc_veto.py | 403 ++++++++++++++ tests/test_segmenter_guarded.py | 134 +++++ 10 files changed, 1323 insertions(+), 7 deletions(-) create mode 100644 src/raw2features/segmenters/grandqc_veto.py create mode 100644 src/raw2features/segmenters/guarded.py create mode 100644 tests/test_segmenter_grandqc_veto.py create mode 100644 tests/test_segmenter_guarded.py diff --git a/docs/SEGMENTATION.md b/docs/SEGMENTATION.md index e286e1f..0b5a3b6 100644 --- a/docs/SEGMENTATION.md +++ b/docs/SEGMENTATION.md @@ -14,6 +14,8 @@ from standard techniques - no GPL/CLAM code, no model weights. | `canny` | **low**-threshold Canny (`low=0.05`) → close → **hole-aware contour fill** (`max_hole_frac` keeps large cavities as background) | good on **cervical / faint** tissue - recovers pale tissue Otsu can under-segment | | `canny_density` | the same low-threshold Canny kept as a **raw edge map** (no dilate/close/fill) at a finer `seg_mpp` (2.0) → the patcher's per-cell mean *is* the patch's **edge density** | **per-patch** decisions with sharp, gap-respecting boundaries - fragmented / sparse slides where filling would swallow the background between fragments. **`--tissue-threshold` becomes a density** (try 0.01-0.05), not a coverage fraction | | `combined` | fuse `otsu` and `canny` - `or` (default) recovers faint tissue, `and` suppresses smooth high-saturation artefacts | when a cohort needs faint-tissue recovery (`or`) or artefact suppression (`and`) | +| `guarded` | `otsu` OR `canny`, minus the detected **scanner canvas / scan-region border** (the canvas is what drags Otsu's global threshold up *and* becomes Canny's outer contour, so the real tissue gets carved out as a hole) | **dirty cohorts, no weights needed** - corner/edge artefacts being tiled *and* pale tissue being missed. Pure OpenCV; this is `grandqc_veto` without the model | +| `grandqc_veto` *(`[grandqc]` extra)* | `otsu` OR `canny`, minus the scanner canvas / scan-region border, then **vetoed** component-wise by GrandQC's tissue stage | dirty cohorts: corner/edge artefacts being tiled *and* pale tissue being missed. Needs the non-commercial GrandQC weights | ## Choosing a segmenter (evidence) @@ -129,3 +131,43 @@ preprocessing we use. Run GrandQC on a sample of your own data and eyeball the c out-of-domain case: on those colorectal slides it cut a ~100%-out-of-focus call to ~30-40%, and a milder 28% case to <1%. It's opt-in (GrandQC's default does no such thing); enable it when your stain sits outside GrandQC's training domain. + +### `grandqc_veto` - classical union, GrandQC as a veto + +`--segmenter grandqc_veto` exists for cohorts where Otsu fails in two directions at once: +a dark corner artefact is densely tiled (and drags the global threshold up, so real +tissue drops out), while pale stroma is missed. Each term does one job: + +- **`otsu`** gives sharp boundaries and correct holes, and defines the mask's grid. +- **`canny`** (contour fill) adds recall on pale, low-saturation tissue. +- a **scan-region guard** removes the synthetic canvas fill the store pads the scan + region with, the border ring around it, and the image frame. It is applied to the + *Canny edge map*, before the dilate/close/fill - once that rectangle has been filled + it is the outer contour, and the real tissue inside becomes a `max_hole_frac` hole + that gets carved back to background, so filtering the finished mask loses the tissue + along with the ring. +- **GrandQC's tissue stage @ 10 µm/px** is used *only as a veto*, per connected + component of the union, so its coarse blobs never add area and its over-inclusion is + clipped by the classical terms. Its trained scale is not configurable here: running + the checkpoint finer is off-domain and measured non-monotonic. + +A component GrandQC rejects is dropped whole; one it endorses is kept intact; only +*large, mixed-evidence* components are clipped pixel-wise - and that clip first dilates +the veto by `veto_slack_um` (so quantisation and ~one mask pixel of model error cannot +erode anything) and then puts back any removed region narrower than +`removal_min_width_um` (a real artefact is a blob, a mis-registration is a sliver). The +two obvious designs are corner cases of the same parameters: `clip_min_area_mm2=inf` is +a pure component veto; `veto_cov=0, keep_cov=1.1, clip_min_area_mm2=0, veto_slack_um=0, +removal_min_width_um=0` is a pure pixel-wise AND. + +The output is binary, so `--tissue-threshold` keeps its usual coverage meaning and tile +counts stay comparable with `otsu`/`canny` runs. After `segment()`, `last_stats` carries +`veto_removed_frac` and per-branch component counts; sorting a cohort by +`veto_removed_frac` is the cheapest way to find slides where the veto is wrong. + +**It does not fix out-of-focus tissue.** Neither GrandQC stage separates focus (the +artifact stage scores `out_of_focus` .027 on controls vs .024 on flagged slides), and +nothing else here measures it. That needs a separate per-patch focus gate. + +If the veto cannot run (missing extra, download failure) it raises by default; +`on_veto_error="warn"` degrades to the bare guarded union, never silently. diff --git a/pyproject.toml b/pyproject.toml index 3157651..e5c53fd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -61,8 +61,10 @@ otsu = "raw2features.segmenters.otsu:OtsuSegmenter" canny = "raw2features.segmenters.canny:CannySegmenter" canny_density = "raw2features.segmenters.canny_density:CannyDensitySegmenter" combined = "raw2features.segmenters.combined:CombinedSegmenter" +guarded = "raw2features.segmenters.guarded:GuardedUnionSegmenter" nuclear = "raw2features.segmenters.nuclear:NuclearSegmenter" grandqc = "raw2features.segmenters.grandqc:GrandQCSegmenter" +grandqc_veto = "raw2features.segmenters.grandqc_veto:GrandQCVetoSegmenter" [project.entry-points."raw2features.patchers"] grid = "raw2features.patcher.grid:GridPatcher" diff --git a/src/raw2features/qc/grandqc.py b/src/raw2features/qc/grandqc.py index cf3bd42..10dce99 100644 --- a/src/raw2features/qc/grandqc.py +++ b/src/raw2features/qc/grandqc.py @@ -88,6 +88,15 @@ def _arch_of(sd: dict) -> str: return "unetplusplus" if any(".blocks.x_" in k for k in sd) else "unet" +# Process-level model cache, keyed by (checkpoint url, device). The pipeline builds a +# fresh segmenter -- and therefore a fresh GrandQC -- for every slide, so without this +# the 26.6 MB tissue checkpoint is re-read, re-verified and re-built on every slide of a +# cohort. The modules are .eval() and only used under torch.inference_mode(), so sharing +# one per device is safe; the device is in the key so the multi-GPU slide-parallel path +# keeps its own copy. No eviction: a couple of tens of MB per device. +_MODEL_CACHE: dict[tuple[str, str], tuple[object, int]] = {} + + def _load(spec: dict, device: str): """Download + load a GrandQC checkpoint; return ``(model, n_classes)``. @@ -95,7 +104,13 @@ def _load(spec: dict, device: str): checkpoint, never assumed. Some GrandQC checkpoints (the artifact stage) are full-model pickles (smp 0.3.1 / timm 0.4.12) that cannot load here; convert once to a clean ``.statedict.pt`` (recipe in ``docs/SEGMENTATION.md``), preferred. + + Cached per ``(url, device)`` -- see :data:`_MODEL_CACHE`. """ + key = (spec["url"], str(device)) + if key in _MODEL_CACHE: + return _MODEL_CACHE[key] + import os import torch @@ -121,7 +136,8 @@ def _load(spec: dict, device: str): classes = int(sd[head].shape[0]) if head is not None else len(QC_CLASSES) + 1 model = _build(classes, _arch_of(sd)) model.load_state_dict(sd, strict=True) - return model.to(device).eval(), classes + _MODEL_CACHE[key] = (model.to(device).eval(), classes) + return _MODEL_CACHE[key] def _read_at_mpp(reader, model_mpp: float): @@ -178,10 +194,16 @@ class GrandQC: def __init__( self, device: str = "cpu", artifact_mpp: str = "1.5", - stain_norm: str | None = None, + stain_norm: str | None = None, tissue_mpp: float | None = None, ) -> None: self.device = device self.artifact_mpp = str(artifact_mpp) + # Scale the TISSUE stage is run at. ``None`` = the checkpoint's own training + # scale (10 µm/px). The UNet++ is fully convolutional so a finer value works + # mechanically and resolves detail the 10 µm/px mask cannot (small lumens / + # gaps between fragments), but it is OFF-DOMAIN for the checkpoint -- verify + # on your own slides before trusting it cohort-wide. + self.tissue_mpp = tissue_mpp self.stain_norm = stain_norm # None|"macenko": normalize before inference self._tissue = None self._artifact = None @@ -212,7 +234,7 @@ def tissue_mask(self, reader): if self._tissue is None: self._tissue, _ = _load(_TISSUE, self.device) - img, ds = self._read(reader, _TISSUE["mpp"]) + img, ds = self._read(reader, self.tissue_mpp or _TISSUE["mpp"]) raster = _segment(self._tissue, img, self.device) mask = (raster == _TISSUE_CLASS).astype(np.float32) return TissueMask(mask=mask, level=0, downsample=ds) # level-0 px per mask px diff --git a/src/raw2features/segmenters/canny.py b/src/raw2features/segmenters/canny.py index 3d251ea..9e7fd7f 100644 --- a/src/raw2features/segmenters/canny.py +++ b/src/raw2features/segmenters/canny.py @@ -17,6 +17,7 @@ from __future__ import annotations +from collections.abc import Callable from typing import TYPE_CHECKING import numpy as np @@ -56,6 +57,16 @@ class CannySegmenter(Segmenter): Internal cavities (lumens / glands / fat) larger than this fraction of the image are kept as **background** rather than filled in; smaller holes (gaps in the edge map) are filled. ``1.0`` fills all holes; ``0.0`` keeps all holes open. + edge_filter: + Optional ``f(edges, img, level_mpp) -> edges`` applied to the **raw Canny edge + map**, before the dilate/close/contour-fill. This is the only useful place to + suppress an edge you do not want: once a scan-region rectangle or a slide-edge + ring has been *filled*, it is the outer contour and the real tissue inside it + becomes a ``max_hole_frac``-sized hole that gets carved back to background -- + i.e. filtering the finished mask instead loses the tissue as well as the ring. + ``level_mpp`` is the microns/px of the level the edges were computed at (or + ``None`` when the slide has no MPP), so a filter can size itself in microns. + See :func:`raw2features.segmenters.grandqc_veto.scan_region_kill`. """ name = "canny" @@ -70,6 +81,8 @@ def __init__( close_kernel: int = 9, min_component_frac: float = 0.001, max_hole_frac: float = 0.01, + edge_filter: Callable[[np.ndarray, np.ndarray, float | None], np.ndarray] + | None = None, ) -> None: self.seg_mpp = seg_mpp self.blur = blur if blur % 2 == 1 else blur + 1 @@ -79,6 +92,7 @@ def __init__( self.close_kernel = close_kernel self.min_component_frac = min_component_frac self.max_hole_frac = max_hole_frac + self.edge_filter = edge_filter def _pick_level(self, reader: WSISource) -> int: return nearest_level(reader.mpp, reader.level_downsamples(), self.seg_mpp) @@ -98,6 +112,13 @@ def segment(self, reader: WSISource) -> TissueMask: hi = int(round(self.high * 255)) edges = cv2.Canny(gray, lo, hi) + # Suppress unwanted edges HERE, while they are still thin lines -- see the + # ``edge_filter`` docstring for why doing it after the fill does not work. + if self.edge_filter is not None: + mpp = reader.mpp + level_mpp = mpp * reader.level_downsamples()[level] if mpp else None + edges = self.edge_filter(edges, img, level_mpp) + # Connect fragments, then close outlines so they enclose fillable regions. if self.dilate_kernel > 0: edges = cv2.dilate( diff --git a/src/raw2features/segmenters/grandqc.py b/src/raw2features/segmenters/grandqc.py index 50ae6de..b24e2c5 100644 --- a/src/raw2features/segmenters/grandqc.py +++ b/src/raw2features/segmenters/grandqc.py @@ -25,13 +25,21 @@ class GrandQCSegmenter(Segmenter): name = "grandqc" - def __init__(self, device: str = "auto") -> None: + def __init__(self, device: str = "auto", mpp: float | None = None) -> None: # "auto" resolves to the GPU when available; GrandQC is a UNet++ that runs on # three-channel (RGB) patches, so it benefits from a GPU. self.device = device + # ``mpp`` overrides the scale the tissue stage runs at. The checkpoint is + # trained at 10 µm/px, where one mask pixel spans 20x20 patch-grid pixels at + # 0.5 µm/px -- so small lumens and the gaps between tissue fragments are below + # its resolution and get predicted as tissue. A finer value (e.g. 2.0) sharpens + # the mask but runs the net off-domain; check it on your own slides. + self.mpp = mpp def segment(self, reader: WSISource) -> TissueMask: from raw2features.core.device import resolve_device from raw2features.qc.grandqc import GrandQC - return GrandQC(device=resolve_device(self.device)).tissue_mask(reader) + return GrandQC( + device=resolve_device(self.device), tissue_mpp=self.mpp + ).tissue_mask(reader) diff --git a/src/raw2features/segmenters/grandqc_veto.py b/src/raw2features/segmenters/grandqc_veto.py new file mode 100644 index 0000000..46eb5ca --- /dev/null +++ b/src/raw2features/segmenters/grandqc_veto.py @@ -0,0 +1,562 @@ +"""Vetoed-union segmenter: ``(otsu | canny)``, scan-region guarded, GrandQC vetoed. + +Three cheap, classical stages and one deep stage, each used only where it is reliable. + +**What each term contributes** + +* ``otsu`` (HSV-saturation Otsu @ 8 µm/px) -- *sharp boundaries and correct holes*. + It defines the reference grid. It is the conservative term: on a clean slide the + output's outline is Otsu's outline. +* ``canny`` (low-threshold Canny -> close -> fill contours @ 8 µm/px) -- *recall on + pale / low-saturation tissue*, which thresholding saturation drops (cervical stroma, + curettings). Filled, so it is a *region* mask: the union and the per-component + algebra below both need regions, which is why ``canny_density``'s raw edge map is + deliberately not used here (its mask is a density, and OR-ing a density into a binary + mask silently rescales ``--tissue-threshold``). +* the **scan-region guard** (pure OpenCV, ~10 ms) -- *removes the synthetic canvas fill, + the scan-region border and the image frame*. This is not cosmetic. Canny latches onto + the canvas/scan-region step edge, that rectangle becomes the outer contour, and the + real tissue inside it then becomes a ``max_hole_frac``-sized hole that + :class:`~raw2features.segmenters.canny.CannySegmenter` carves back to background -- + so on canvas-heavy slides the border does not merely add junk, it *destroys* the + canny term. The guard therefore runs on the **edge map**, before the morphology + (via ``CannySegmenter(edge_filter=...)``); guarding the finished mask instead is + measurably useless. +* ``grandqc`` (stage-1 tissue UNet++ @ its trained 10 µm/px) -- **veto only, never a + source of tissue**. Its mask bleeds into glass and fills lumens (one mask pixel spans + 20x20 grid px at 0.5 µm/px), so it must never *add* area; but it is excellent at + saying "that block is not tissue". Used per connected component of the union, so its + coarse resolution stops mattering. + +**What the veto does NOT fix** + +* **Out-of-focus tissue.** Neither GrandQC stage separates focus (the artifact stage + scores ``out_of_focus`` .027 on controls vs .024 on flagged slides -- no separation; + the tissue stage happily calls blurred tissue tissue). Defocused regions are also + low-contrast, so the union terms are weak there, but *nothing here measures focus*. + Slides that are tiled because they are blurry stay tiled. That needs a separate + per-patch focus gate (e.g. variance-of-Laplacian), not this segmenter. +* **A GrandQC false negative.** A tissue fragment GrandQC misses forms its own + component with coverage ~0 and is dropped wholesale -- no slack or width guard + protects that path (they only guard the pixel-wise clip). ``veto_removed_frac`` in + :attr:`GrandQCVetoSegmenter.last_stats` exists so a cohort can be sorted by it and + the top of the distribution eyeballed. +* **Otsu's threshold being dragged up by an artifact.** Where a dark block pushes the + global Otsu threshold above real tissue, the veto can only *remove*; it cannot restore + what Otsu already lost. Recall on those slides is carried entirely by the guarded + canny term. (The principled fix -- re-estimating Otsu on the veto-endorsed region -- + is deliberately out of scope here.) + +**The combination algebra** is a hybrid that contains the two obvious designs as corner +cases, so they can be compared without a second implementation: + +* pure pixel-wise ``AND``: ``veto_cov=0.0, keep_cov=1.1, clip_min_area_mm2=0.0, + veto_slack_um=0.0, removal_min_width_um=0.0`` +* pure component veto: ``clip_min_area_mm2=math.inf`` (the clip branch never fires) +* the default: drop a component GrandQC rejects outright, keep one it endorses, and + pixel-wise clip only the *large, mixed-evidence* ones -- because boundary-erosion cost + scales as ~1/diameter (cheap on big components, destructive on thin ones) while the + risk of an artifact being fused into a real component grows with component size. + +Two guards make that clip safe: the veto mask is dilated by ``veto_slack_um`` first (so +quantisation and ~one mask pixel of model error cannot erode anything), and the removal +set is filtered by opening-by-reconstruction at ``removal_min_width_um`` (a real +artifact is a blob; a mis-registration is a sliver, and slivers are put back). + +The output is **strictly binary**, so ``--tissue-threshold`` keeps its usual +coverage meaning (default 0.1) and tile counts stay comparable with existing +``otsu`` / ``canny`` runs. + +Needs the optional ``[grandqc]`` extra; the weights are CC-BY-NC-SA (non-commercial) and +are fetched on first use -- see ``docs/SEGMENTATION.md``, ``docs/MODEL_LICENSES.md``. +""" + +from __future__ import annotations + +import math +import warnings +from typing import TYPE_CHECKING, Any + +import numpy as np + +from raw2features.core.mpp import nearest_level +from raw2features.core.plugins import register + +from .base import Segmenter, TissueMask +from .canny import CannySegmenter +from .grandqc import GrandQCSegmenter +from .otsu import OtsuSegmenter + +if TYPE_CHECKING: + from raw2features.readers.base import WSISource + +__all__ = ["GrandQCVetoSegmenter", "scan_region_kill"] + + +# -- small OpenCV helpers ------------------------------------------------------ + + +def _disk(radius: int) -> np.ndarray: + import cv2 + + r = max(1, int(radius)) + return cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2 * r + 1, 2 * r + 1)) + + +def _px(um: float, mpp: float | None) -> int: + """``um`` microns as a whole number of mask pixels at ``mpp`` µm/px. + + ``0`` only when the caller asked for no margin at all (``um <= 0``) -- callers + treat that as "skip this morphology entirely", which is what makes the documented + pixel-wise-AND corner case exact. Any positive request rounds up to >= 1 px. + """ + if um <= 0: + return 0 + if not mpp: # unknown MPP: fall back to a single pixel rather than guessing + return 1 + return max(1, int(round(um / mpp))) + + +def _align(m: np.ndarray, shape: tuple[int, int]) -> np.ndarray: + """Resample a boolean mask onto ``shape`` (H, W), preserving thin structures. + + Upsampling uses nearest neighbour; **downsampling uses ``INTER_AREA`` on a float + mask thresholded at > 0**, i.e. OR-pooling -- a structure occupying any part of a + destination pixel survives. ``INTER_NEAREST`` on the downsample path deletes thin + structures at random, which is exactly the tissue this segmenter exists to recover. + """ + import cv2 + + if m.shape == shape: + return m + h, w = shape + interp = ( + cv2.INTER_NEAREST if (m.shape[0] <= h and m.shape[1] <= w) else cv2.INTER_AREA + ) + return cv2.resize(m.astype(np.float32), (w, h), interpolation=interp) > 0 + + +def _frame_connected(m: np.ndarray) -> np.ndarray: + """Components of ``m`` (8-connected) that touch the image perimeter.""" + import cv2 + + if not m.any(): + return np.zeros(m.shape, bool) + n, lab = cv2.connectedComponents(m.astype(np.uint8), connectivity=8) + if n <= 1: + return np.zeros(m.shape, bool) + border = np.concatenate([lab[0], lab[-1], lab[:, 0], lab[:, -1]]) + keep = np.unique(border) + keep = keep[keep != 0] + if keep.size == 0: + return np.zeros(m.shape, bool) + lut = np.zeros(n, np.uint8) + lut[keep] = 1 + return lut[lab].astype(bool) + + +def _reconstruct(mask: np.ndarray, seeds: np.ndarray) -> np.ndarray: + """Morphological reconstruction: components of ``mask`` containing a seed pixel.""" + import cv2 + + if not seeds.any(): + return np.zeros(mask.shape, bool) + n, lab = cv2.connectedComponents(mask.astype(np.uint8), connectivity=4) + hit = np.unique(lab[seeds]) + hit = hit[hit != 0] + if hit.size == 0: + return np.zeros(mask.shape, bool) + lut = np.zeros(n, np.uint8) + lut[hit] = 1 + return lut[lab].astype(bool) + + +# -- the scan-region / canvas guard ------------------------------------------- + + +def scan_region_kill( + img: np.ndarray, + level_mpp: float | None, + *, + canvas_dark: int = 20, + canvas_bright: int | None = None, + canvas_var: float = 0.5, + border_um: float = 200.0, + frame_um: float = 24.0, +) -> np.ndarray: + """Pixels that are scanner canvas, scan-region border, or image frame. + + The border a Canny segmenter picks up on a WSI is not an abstract "edge of the + image" -- it is the step edge at the boundary of the **synthetic canvas fill** the + store pads the scan region with (black in OME-Zarr v0.5; white in some SVS/CZI + exports). Detecting the fill itself is exact and free, whereas a fixed border + *margin* both misses interior scan-region rectangles and clips genuine tissue on + slides where the specimen runs to the frame. + + A pixel is canvas when it is (a) a **no-signal fill** -- near-constant in a 5x5 + neighbourhood *and* dark (or bright, if ``canvas_bright`` is set) -- and (b) + **connected to the image perimeter**. The local-variance term is what keeps a dark + *tissue* artifact out of the flood: real material is textured. The perimeter + connectivity is what keeps an interior dark lumen out of it. + + Returns a boolean mask of canvas | a ``border_um`` band dilated from the canvas + boundary (the ring the filled canny mask actually occupies) | a ``frame_um`` + perimeter ring (for slides with no canvas at all). + """ + import cv2 + + gray = cv2.cvtColor(img[..., :3], cv2.COLOR_RGB2GRAY) + g32 = gray.astype(np.float32) + mean = cv2.boxFilter(g32, -1, (5, 5)) + var = np.maximum(cv2.boxFilter(g32 * g32, -1, (5, 5)) - mean * mean, 0.0) + + fill = gray <= canvas_dark + if canvas_bright is not None: + fill |= gray >= canvas_bright + canvas = _frame_connected(fill & (var <= canvas_var)) + + kill = canvas.copy() + r = _px(border_um, level_mpp) + if canvas.any() and r > 0: + grad = cv2.morphologyEx( + canvas.astype(np.uint8), cv2.MORPH_GRADIENT, np.ones((3, 3), np.uint8) + ) + kill |= cv2.dilate(grad, _disk(r)) > 0 + f = _px(frame_um, level_mpp) + if f > 0: + kill[:f], kill[-f:], kill[:, :f], kill[:, -f:] = True, True, True, True + return kill + + +# -- the segmenter ------------------------------------------------------------- + + +@register("segmenters", "grandqc_veto") +class GrandQCVetoSegmenter(Segmenter): + """``(otsu | canny)``, scan-region guarded, vetoed component-wise by GrandQC. + + Parameters + ---------- + seg_mpp: + Microns/px of the **union grid** -- the level otsu and canny are computed at, + and the grid of the returned mask. 8.0 is a region-outline grid, not a decision + surface; do not lower it expecting sharper patches. + device: + Torch device for the GrandQC veto (``"auto"`` resolves to the GPU). The veto is + a UNet++ forward and wants a GPU, though it is only ~2-4% of a slide's + embedding cost either way and runs acceptably on CPU (~5-20 s/slide). + canny_min_component_frac: + Passthrough to :class:`~raw2features.segmenters.canny.CannySegmenter`. Left at + canny's own default so a comparison attributes changes to the fusion rather + than to a retuned term. With the veto acting as the junk filter, lowering it + (0.0005 / 0.0002) to recover small pale fragments becomes much safer. + canny_max_area_frac: + Degenerate-fill backstop. If the filled canny mask still covers more than this + fraction of the level after the guard, treat the canny term as having collapsed + into one giant contour and fall back to otsu alone for that slide. + guard: + Run the scan-region / canvas guard (see :func:`scan_region_kill`). On by + default: it is ~10 ms, it only ever removes, and on canvas-heavy slides it is + the difference between the canny term working and the canny term being empty. + canvas_dark, canvas_bright, canvas_var, border_um, frame_um: + Guard parameters, see :func:`scan_region_kill`. ``canvas_bright`` is ``None`` + (white-fill detection **off**) by default -- a JPEG-flattened slide with a + uniform white background connected to the frame would otherwise be flooded as + canvas and lose a ``border_um`` ring around every tissue island. Enable it per + format only after checking how much of the Otsu mask the guard removes. + veto_cov: + Drop a union component outright when GrandQC endorses less than this fraction + of it. Deliberately low: only remove what GrandQC *confidently* rejects, because + this branch has no safety net (a GrandQC false negative is deleted wholesale). + keep_cov: + Skip the clip when GrandQC endorses essentially all of a component -- a + short-circuit, not a safety mechanism. High on purpose, so an artifact attached + to a 95%-endorsed component still gets clipped. + clip_min_area_mm2: + Only pixel-wise clip components at least this large. 0.25 mm² is a 500x500 µm + square; erosion of a compact component costs about ``4*slack/diameter``, so at + ``veto_slack_um=16`` this keeps the loss under ~10%. ``math.inf`` disables the + clip entirely (pure component veto); ``0.0`` clips everything. + veto_slack_um: + Dilate the veto mask by this much **before** any pixel-wise clip, so the clip + can only remove pixels that are this far from anything GrandQC calls tissue. + Absorbs grid quantisation (half a 10 µm mask pixel) plus a full mask pixel of + model boundary error. Cheap, because GrandQC's error direction is documented + over-inclusion and the union term -- not the veto -- defines the outer boundary. + removal_min_width_um: + A clipped-away region must be at least this wide *somewhere* + (opening-by-reconstruction) or it is returned to tissue. This protects a thin + tissue wisp of any length, which a pure area threshold on removed pieces would + not. The trade is explicit: a genuine artifact thinner than this everywhere + (a hairline pen line or fold) is put back. + connectivity: + Component connectivity, 4 or 8. 4 on purpose: a single diagonal pixel chain + must not merge an artifact block into the tissue component. + veto_abort_cov: + Safety valve. If GrandQC endorses less than this fraction of the whole union, + treat the veto as degenerate (empty prediction / model or stain failure), warn, + and emit the union unclipped rather than silently returning zero tiles. Set low + on purpose: on a genuine artifact slide the veto legitimately removes most of + the union's *area*. + min_area_mm2: + Optional dust filter on the final mask. ``0.0`` (off) by default so it does not + confound an evaluation of the fusion itself. + on_veto_error: + ``"raise"`` (default) if the GrandQC term cannot run -- missing ``[grandqc]`` + extra, download failure, OOM. ``"warn"`` degrades to the bare guarded union. + Never silent: a silent fallback reintroduces every failure mode this exists to + fix. + + Notes + ----- + After :meth:`segment`, :attr:`last_stats` holds a provenance/observability dict + (``veto_removed_frac``, per-branch component counts, whether the canny term or the + veto degenerated, the grid MPPs). Sorting a cohort by ``veto_removed_frac`` is the + cheapest way to find slides where the veto is wrong. + """ + + name = "grandqc_veto" + + def __init__( + self, + seg_mpp: float = 8.0, + device: str = "auto", + *, + canny_min_component_frac: float = 0.001, + canny_max_area_frac: float = 0.90, + guard: bool = True, + canvas_dark: int = 20, + canvas_bright: int | None = None, + canvas_var: float = 0.5, + border_um: float = 200.0, + frame_um: float = 24.0, + veto_cov: float = 0.10, + keep_cov: float = 0.98, + clip_min_area_mm2: float = 0.25, + veto_slack_um: float = 16.0, + removal_min_width_um: float = 32.0, + connectivity: int = 4, + veto_abort_cov: float = 0.02, + min_area_mm2: float = 0.0, + on_veto_error: str = "raise", + ) -> None: + if connectivity not in (4, 8): + raise ValueError(f"connectivity must be 4 or 8, got {connectivity!r}") + if on_veto_error not in ("raise", "warn"): + raise ValueError( + f"on_veto_error must be 'raise' or 'warn', got {on_veto_error!r}" + ) + self.seg_mpp = seg_mpp + self.device = device + self.canny_min_component_frac = canny_min_component_frac + self.canny_max_area_frac = canny_max_area_frac + self.guard = guard + self.canvas_dark = canvas_dark + self.canvas_bright = canvas_bright + self.canvas_var = canvas_var + self.border_um = border_um + self.frame_um = frame_um + self.veto_cov = veto_cov + self.keep_cov = keep_cov + self.clip_min_area_mm2 = clip_min_area_mm2 + self.veto_slack_um = veto_slack_um + self.removal_min_width_um = removal_min_width_um + self.connectivity = connectivity + self.veto_abort_cov = veto_abort_cov + self.min_area_mm2 = min_area_mm2 + self.on_veto_error = on_veto_error + self.last_stats: dict[str, Any] = {} + + # -- the single veto call site (subclass / monkeypatch seam) --------------- + + def _veto_mask(self, reader: WSISource) -> TissueMask: + """GrandQC's stage-1 tissue mask, **always at its trained 10 µm/px**. + + ``mpp`` is deliberately not forwarded: running the checkpoint finer is + off-domain and measured non-monotonic (on one cohort a control went + 6210 -> 6343 -> 4194 tiles and an artifact slide went 1405 -> 2329 -> 1956 as + the MPP was refined), so a "sharper" veto is not a better veto. + + Everything deep sits behind this one call, so a caller that wants to memoise, + batch or gate the model can override it without touching the algebra. + """ + return GrandQCSegmenter(device=self.device).segment(reader) + + # -- terms ----------------------------------------------------------------- + + def _union(self, reader: WSISource) -> tuple[np.ndarray, int, float, np.ndarray]: + """Guarded ``otsu | canny`` on otsu's grid; returns (U, level, ds, kill).""" + from raw2features.viz import read_level_capped + + level = nearest_level(reader.mpp, reader.level_downsamples(), self.seg_mpp) + # Otsu reads through read_level_capped, so its downsample may be the level's + # downsample times an extra factor on a deficient pyramid. Always take the + # reference grid from the Otsu result, never from level_downsamples()[level]. + img, factor = read_level_capped(reader, level) + ds = float(reader.level_downsamples()[level]) * factor + union_mpp = reader.mpp * ds if reader.mpp else None + + kill = ( + scan_region_kill( + img, + union_mpp, + canvas_dark=self.canvas_dark, + canvas_bright=self.canvas_bright, + canvas_var=self.canvas_var, + border_um=self.border_um, + frame_um=self.frame_um, + ) + if self.guard + else np.zeros(img.shape[:2], bool) + ) + + otsu = OtsuSegmenter(seg_mpp=self.seg_mpp).segment(reader) + shape = otsu.mask.shape + kill_ref = _align(kill, shape) + + def _edge_filter(edges: np.ndarray, _img: np.ndarray, _mpp: float | None): + k = _align(kill_ref, edges.shape) + out = edges.copy() + out[k] = 0 + return out + + canny = CannySegmenter( + seg_mpp=self.seg_mpp, + min_component_frac=self.canny_min_component_frac, + edge_filter=_edge_filter if self.guard else None, + ).segment(reader) + cm = _align(canny.mask > 0, shape) + degenerate = bool(cm.mean() > self.canny_max_area_frac) + if degenerate: # a filled outer ring swallowed the slide -> otsu alone + cm = np.zeros(shape, bool) + self.last_stats["canny_degenerate"] = degenerate + + u = ((otsu.mask > 0) | cm) & ~kill_ref + self.last_stats["kill_frac"] = float(kill_ref.mean()) + return u, otsu.level, otsu.downsample, kill_ref + + # -- the seam -------------------------------------------------------------- + + def segment(self, reader: WSISource) -> TissueMask: + import cv2 + + self.last_stats = {} + u, level, ds, _kill = self._union(reader) + union_mpp = reader.mpp * ds if reader.mpp else None + self.last_stats.update(union_mpp=union_mpp, level=level) + + if not u.any(): # nothing to veto; skip the model entirely (provably exact) + self.last_stats.update( + n_components=0, n_dropped=0, n_clipped=0, n_kept_intact=0, + n_kept_small=0, veto_aborted=False, veto_removed_frac=0.0, + ) + return TissueMask(np.zeros(u.shape, np.float32), level, ds) + + try: + veto = self._veto_mask(reader) + # Derived, never asserted: if _veto_mask is ever changed to run the + # checkpoint off its trained 10 um/px this records the real value (and + # the test that pins it to 10.0 then fails, as it should). + self.last_stats["veto_mpp"] = ( + reader.mpp * veto.downsample if reader.mpp else None + ) + except Exception as e: # noqa: BLE001 - re-raised or explicitly degraded + if self.on_veto_error == "raise": + raise + warnings.warn( + f"{self.name}: GrandQC veto unavailable ({e!r}); emitting the " + "unvetoed union -- artifacts will NOT be removed.", + RuntimeWarning, + stacklevel=2, + ) + self.last_stats.update(veto_aborted=True, veto_removed_frac=0.0) + return TissueMask(u.astype(np.float32), level, ds) + + g = _align(veto.mask > 0, u.shape) + cov_all = float(g[u].mean()) + self.last_stats["union_cov"] = cov_all + if cov_all < self.veto_abort_cov: + warnings.warn( + f"{self.name}: GrandQC endorses only {cov_all:.3f} of the union; " + "treating the veto as degenerate and emitting the union unclipped.", + RuntimeWarning, + stacklevel=2, + ) + self.last_stats.update(veto_aborted=True, veto_removed_frac=0.0) + return TissueMask(u.astype(np.float32), level, ds) + self.last_stats["veto_aborted"] = False + + # Slack-dilated veto: the clip may only remove pixels this far from anything + # GrandQC calls tissue, so quantisation + ~1 mask px of model error is absorbed. + r_slack = _px(self.veto_slack_um, union_mpp) + gd = cv2.dilate(g.astype(np.uint8), _disk(r_slack)) > 0 if r_slack else g + + n, lab, stats, _ = cv2.connectedComponentsWithStats( + u.astype(np.uint8), connectivity=self.connectivity + ) + areas = stats[:, cv2.CC_STAT_AREA].astype(np.float64) + # Endorsed pixels per component, measured on the UNDILATED veto: this is a + # statistic, not a clip, and dilating it would inflate every coverage. + sums = np.bincount(lab[g], minlength=n).astype(np.float64) + cov = np.zeros(n) + cov[1:] = sums[1:] / np.maximum(areas[1:], 1.0) + # A slide with no MPP has no physical area, so every component reads as "small" + # and the segmenter degrades to a pure component veto -- the safe direction, + # since the clip is the only branch that can erode a boundary. + px_mm2 = (union_mpp / 1000.0) ** 2 if union_mpp else 0.0 + area_mm2 = areas * px_mm2 + + drop = cov < self.veto_cov + small = area_mm2 < self.clip_min_area_mm2 + keep_whole = ~drop & ((cov >= self.keep_cov) | small) + clip = ~drop & ~keep_whole + for arr in (drop, keep_whole, clip): + arr[0] = False # label 0 is background + + lut = np.zeros(n, np.uint8) + lut[keep_whole] = 1 + out = lut[lab].astype(bool) + + r_open = ( + max(1, math.ceil(self.removal_min_width_um / 2.0 / union_mpp)) + if (self.removal_min_width_um > 0 and union_mpp) + else 0 + ) + for c in np.nonzero(clip)[0]: + x, y, w, h = ( + stats[c, cv2.CC_STAT_LEFT], stats[c, cv2.CC_STAT_TOP], + stats[c, cv2.CC_STAT_WIDTH], stats[c, cv2.CC_STAT_HEIGHT], + ) + sl = (slice(y, y + h), slice(x, x + w)) + m = lab[sl] == c + removed = m & ~gd[sl] + if r_open: + # A real artifact is a blob; a mis-registration is a sliver. Keep only + # removal regions that are >= removal_min_width_um wide somewhere. + seeds = cv2.morphologyEx( + removed.astype(np.uint8), cv2.MORPH_OPEN, _disk(r_open) + ).astype(bool) + removed = _reconstruct(removed, seeds) + out[sl] |= m & ~removed + + if self.min_area_mm2 > 0 and out.any() and px_mm2 > 0: + n2, lab2, st2, _ = cv2.connectedComponentsWithStats( + out.astype(np.uint8), connectivity=self.connectivity + ) + big = st2[:, cv2.CC_STAT_AREA].astype(np.float64) * px_mm2 >= ( + self.min_area_mm2 + ) + big[0] = False + lut2 = np.zeros(n2, np.uint8) + lut2[big] = 1 + out = lut2[lab2].astype(bool) + + self.last_stats.update( + n_components=int(n - 1), + n_dropped=int(drop.sum()), + n_clipped=int(clip.sum()), + n_kept_intact=int((keep_whole & (cov >= self.keep_cov)).sum()), + n_kept_small=int((keep_whole & small & (cov < self.keep_cov)).sum()), + veto_removed_frac=float(1.0 - out.sum() / max(int(u.sum()), 1)), + ) + return TissueMask(out.astype(np.float32), level, ds) diff --git a/src/raw2features/segmenters/guarded.py b/src/raw2features/segmenters/guarded.py new file mode 100644 index 0000000..c11abf5 --- /dev/null +++ b/src/raw2features/segmenters/guarded.py @@ -0,0 +1,123 @@ +"""Guarded union tissue segmenter: ``otsu | canny`` minus the scanner canvas. + +Pure OpenCV + numpy, no model weights, no GPU -- the same footing as ``otsu``/``canny``. + +## What it fixes + +Two classical-segmenter failures on real WSIs share one cause: the **scanner canvas** +(the dark no-signal fill outside the scanned region) and the **scan-region border**. + +* **Otsu** thresholds HSV saturation *globally*, so a dark saturated canvas corner drags + the threshold up and real tissue falls below it. On the worst observed slide Otsu put + **all** of its patches on a corner artifact and **none** on the biopsy. +* **Canny** is worse in a subtler way: the scan-region rectangle becomes the outer + contour, so the real tissue inside it is then a ``max_hole_frac`` *hole* that gets + carved back to background -- the fill does not merely add junk, it **deletes** tissue. + +Removing the canvas from the **raw edge map** (before Canny's dilate/close/fill) and from +the Otsu mask fixes both. Measured on 101 cervical WSIs: on the artifact slides the tiles +move off the artifact and onto the biopsy (one slide 2 -> 1172 patches, another 0 -> 1097), +pale slides gain ~2x (Canny recovers what Otsu's saturation threshold misses), and the +unflagged controls grow only ~1.2x. + +## What it does NOT fix + +**Out-of-focus tissue.** Blur is invisible to both a colour threshold and an edge filter; +OOF slides simply gain patches here. Use a dedicated focus metric if that matters. + +## Relationship to ``grandqc_veto`` + +``grandqc_veto`` is this segmenter plus a GrandQC-tissue veto. On the same 101-slide panel +the veto changed nothing at all on **79/101** slides and removed only **1.2%** of patches +overall (median relative difference 0.0000), its one real contribution being on +out-of-focus slides. Prefer this segmenter unless you specifically want that margin: it +needs no GPU, no ~25 MB download, and no CC-BY-NC-SA weights. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +from raw2features.core.plugins import register + +from .base import Segmenter, TissueMask +from .grandqc_veto import GrandQCVetoSegmenter + +if TYPE_CHECKING: + from raw2features.readers.base import WSISource + + +@register("segmenters", "guarded") +class GuardedUnionSegmenter(Segmenter): + """``(otsu | canny) & ~scanner_canvas`` -- classical, permissive, no weights. + + Parameters mirror the guard/union half of + :class:`~raw2features.segmenters.grandqc_veto.GrandQCVetoSegmenter`; see + :func:`~raw2features.segmenters.grandqc_veto.scan_region_kill` for the canvas + detector itself. + + Parameters + ---------- + seg_mpp: + Target microns/px for the (cheap, low-res) level the masks are computed at. + guard: + Remove the detected scanner canvas / scan-region border. ``False`` reduces this + to a plain ``otsu | canny`` union (i.e. ``combined or``) and reinstates both + failures above -- it exists so the guard's contribution can be measured. + canvas_dark, canvas_bright, canvas_var, border_um, frame_um: + Canvas-detector parameters, passed through unchanged. + canny_min_component_frac, canny_max_area_frac: + Canny tuning; ``canny_max_area_frac`` is the safety valve for a filled outer + ring swallowing the slide, in which case the Canny term is dropped and Otsu + alone is used. + """ + + name = "guarded" + + def __init__( + self, + seg_mpp: float = 8.0, + *, + guard: bool = True, + canvas_dark: int = 20, + canvas_bright: int | None = None, + canvas_var: float = 0.5, + border_um: float = 200.0, + frame_um: float = 24.0, + canny_min_component_frac: float = 0.001, + canny_max_area_frac: float = 0.90, + ) -> None: + self.seg_mpp = seg_mpp + self.guard = guard + # Compose rather than duplicate: the union half is shared with grandqc_veto, so + # a fix to the canvas detector or the edge-map guard lands in both. The veto + # half is never reached -- `segment` below returns the union directly and no + # GrandQC import, weight download or forward pass happens. + self._impl = GrandQCVetoSegmenter( + seg_mpp=seg_mpp, + guard=guard, + canvas_dark=canvas_dark, + canvas_bright=canvas_bright, + canvas_var=canvas_var, + border_um=border_um, + frame_um=frame_um, + canny_min_component_frac=canny_min_component_frac, + canny_max_area_frac=canny_max_area_frac, + ) + + @property + def last_stats(self) -> dict: + """Provenance from the most recent :meth:`segment` (``kill_frac`` etc.).""" + return self._impl.last_stats + + def segment(self, reader: WSISource) -> TissueMask: + self._impl.last_stats = {} + u, level, ds, _kill = self._impl._union(reader) + self._impl.last_stats.update( + union_mpp=(reader.mpp * ds if reader.mpp else None), + level=level, + veto="none", + ) + return TissueMask(mask=u.astype(np.float32), level=level, downsample=ds) diff --git a/tests/test_segmenter_canny_density.py b/tests/test_segmenter_canny_density.py index 39e3d01..d245f10 100644 --- a/tests/test_segmenter_canny_density.py +++ b/tests/test_segmenter_canny_density.py @@ -6,12 +6,11 @@ import pytest import zarr +from conftest import build_ngff_v04 from raw2features.readers.omezarr import OmeZarrReader from raw2features.segmenters.canny import CannySegmenter from raw2features.segmenters.canny_density import CannyDensitySegmenter -from conftest import build_ngff_v04 - @pytest.fixture def edged_ngff(tmp_path) -> str: diff --git a/tests/test_segmenter_grandqc_veto.py b/tests/test_segmenter_grandqc_veto.py new file mode 100644 index 0000000..57f0546 --- /dev/null +++ b/tests/test_segmenter_grandqc_veto.py @@ -0,0 +1,403 @@ +"""Tests for the vetoed-union segmenter ``(otsu | canny)`` + guard + GrandQC veto. + +Two things make these tests non-vacuous: + +* **The fixtures have real edges.** The shared ``synthetic_ngff`` fixture is a smooth + ramp, so Canny returns nothing on it and every edge/union assertion would pass + trivially. ``zoned_ngff`` below paints three separated blocks whose colour/texture + make Otsu and Canny *disagree*, and ``canvas_ngff`` paints a black scanner canvas + around the scan region. +* **The GrandQC veto is injected.** The real veto needs the ``[grandqc]`` extra, a + GPU and CC-BY-NC-SA weights, so tests patch the single ``_veto_mask`` call site with + a hand-drawn mask on a *coarser* grid (10 µm/px on a 4 µm/px slide), which also + exercises the veto-to-union resampling. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest +import zarr + +from conftest import build_ngff_v04 +from raw2features.readers.omezarr import OmeZarrReader +from raw2features.segmenters.base import TissueMask +from raw2features.segmenters.canny import CannySegmenter +from raw2features.segmenters.grandqc_veto import ( + GrandQCVetoSegmenter, + _align, + scan_region_kill, +) +from raw2features.segmenters.otsu import OtsuSegmenter + +# Level-0 MPP 4.0 with a /2 pyramid -> level MPPs 4 / 8 / 16, so the production +# ``seg_mpp=8.0`` default resolves to level 1 (240x240) and every micron-denominated +# parameter (border_um, veto_slack_um, clip_min_area_mm2, ...) lands in its designed +# regime rather than being scaled into absurdity by a toy MPP. +SIZES = ((480, 480), (240, 240), (120, 120)) +MPP0 = 4.0 +VETO_DS = 2.5 # 10.0 µm/px veto grid / 4.0 µm/px level 0 +VETO_N = int(round(SIZES[0][0] / VETO_DS)) # 192 + +# (row0, row1, col0, col1) as fractions of the level, so they hold at every level. +BOXES = { + "A": (0.10, 0.40, 0.10, 0.40), # saturated + smooth -> otsu YES, canny outline + "B": (0.10, 0.40, 0.60, 0.90), # desaturated + textured -> otsu NO, canny YES + "C": (0.60, 0.90, 0.60, 0.90), # saturated junk -> union YES, veto will reject it +} + + +def _slice(name: str, h: int, w: int): + r0, r1, c0, c1 = BOXES[name] + return slice(int(r0 * h), int(r1 * h)), slice(int(c0 * w), int(c1 * w)) + + +def _frac(mask: np.ndarray, name: str) -> float: + h, w = mask.shape + return float(mask[_slice(name, h, w)].mean()) + + +def _paint_zones(h: int, w: int) -> np.ndarray: + img = np.full((h, w, 3), 230, np.uint8) # desaturated smooth background + img[_slice("A", h, w)] = (200, 40, 60) # saturated block + ys, xs = _slice("B", h, w) + hh, ww = ys.stop - ys.start, xs.stop - xs.start + sq = max(2, h // 40) + yy, xx = np.mgrid[0:hh, 0:ww] + checks = np.where((((yy // sq) + (xx // sq)) % 2).astype(bool), 120, 200) + for c in range(3): # grey checkerboard: no saturation, plenty of edges + img[ys, xs, c] = checks.astype(np.uint8) + img[_slice("C", h, w)] = (40, 60, 200) # a second saturated block + return img + + +def _write(path: str, painter) -> str: + p = build_ngff_v04(path, mpp0=MPP0, sizes=SIZES) + g = zarr.open_group(p, mode="r+", zarr_format=2) + for i, (h, w) in enumerate(SIZES): + img = painter(h, w) + for c in range(3): + g[str(i)][0, c, 0] = img[:, :, c] + return p + + +@pytest.fixture +def zoned_ngff(tmp_path) -> str: + """Three well-separated blocks that make Otsu and Canny disagree.""" + return _write(str(tmp_path / "zoned.zarr"), _paint_zones) + + +CANVAS_MARGIN = 0.06 +TISSUE_BOX = (0.30, 0.70) + + +def _paint_canvas(h: int, w: int) -> np.ndarray: + img = np.full((h, w, 3), 230, np.uint8) + r0, r1 = int(TISSUE_BOX[0] * h), int(TISSUE_BOX[1] * h) + c0, c1 = int(TISSUE_BOX[0] * w), int(TISSUE_BOX[1] * w) + img[r0:r1, c0:c1] = (200, 40, 60) + m = int(CANVAS_MARGIN * h) # black synthetic canvas around the scan region + img[:m] = img[-m:] = 0 + img[:, :m] = img[:, -m:] = 0 + return img + + +@pytest.fixture +def canvas_ngff(tmp_path) -> str: + """A scan region padded with black canvas -- the slide-edge failure mode.""" + return _write(str(tmp_path / "canvas.zarr"), _paint_canvas) + + +@pytest.fixture +def blank_ngff(tmp_path) -> str: + """Uniform white glass: no saturation, no edges, so the union is empty.""" + return _write( + str(tmp_path / "blank.zarr"), + lambda h, w: np.full((h, w, 3), 255, np.uint8), + ) + + +def _veto(monkeypatch, arr: np.ndarray, counter: list | None = None): + """Patch the single veto call site with a hand-drawn 10 µm/px mask.""" + + def fake(self, reader): + if counter is not None: + counter.append(1) + return TissueMask(mask=arr.astype(np.float32), level=0, downsample=VETO_DS) + + monkeypatch.setattr(GrandQCVetoSegmenter, "_veto_mask", fake) + + +def _boxes(*names: str, pad: float = 0.04) -> np.ndarray: + m = np.zeros((VETO_N, VETO_N), np.float32) + for n in names: + r0, r1, c0, c1 = BOXES[n] + m[ + int((r0 - pad) * VETO_N) : int((r1 + pad) * VETO_N), + int((c0 - pad) * VETO_N) : int((c1 + pad) * VETO_N), + ] = 1.0 + return m + + +# -- the seam ------------------------------------------------------------------ + + +def test_returns_a_valid_binary_tissue_mask_on_the_otsu_grid(zoned_ngff, monkeypatch): + _veto(monkeypatch, _boxes("A", "B")) + seg = GrandQCVetoSegmenter() + with OmeZarrReader(zoned_ngff) as r: + ref = OtsuSegmenter(seg_mpp=8.0).segment(r) + tm = seg.segment(r) + + assert isinstance(tm, TissueMask) + assert tm.mask.dtype == np.float32 and tm.mask.ndim == 2 + # Strictly binary: --tissue-threshold must keep its coverage meaning, so tile + # counts stay comparable with plain otsu / canny runs. + assert set(np.unique(tm.mask)).issubset({0.0, 1.0}) + # The reference grid is Otsu's, including any read_level_capped downscale factor. + assert tm.mask.shape == ref.mask.shape + assert tm.level == ref.level and tm.downsample == ref.downsample + assert seg.last_stats["veto_mpp"] == 10.0 # never run off its trained scale + + +def test_registered_and_no_arg_constructible(): + from raw2features.core import plugins + + seg = plugins.get("segmenters", "grandqc_veto")() # pipeline builds with no args + assert seg.name == "grandqc_veto" + + +def test_rejects_nonsense_construction(): + with pytest.raises(ValueError): + GrandQCVetoSegmenter(connectivity=6) + with pytest.raises(ValueError): + GrandQCVetoSegmenter(on_veto_error="ignore") + + +# -- the union term ------------------------------------------------------------ + + +def test_union_recovers_pale_tissue_that_otsu_alone_misses(zoned_ngff, monkeypatch): + """Zone B is grey (no saturation) but textured: Otsu drops it, Canny fills it.""" + _veto(monkeypatch, _boxes("A", "B")) + with OmeZarrReader(zoned_ngff) as r: + otsu = OtsuSegmenter(seg_mpp=8.0).segment(r).mask > 0 + out = GrandQCVetoSegmenter().segment(r).mask > 0 + + assert _frac(otsu, "B") < 0.01, "fixture must be one Otsu genuinely misses" + assert _frac(out, "B") > 0.90, "the canny term should recover it" + assert _frac(otsu, "A") > 0.90 and _frac(out, "A") > 0.90 # otsu's own kept + + +# -- the veto ------------------------------------------------------------------ + + +def test_veto_drops_a_component_the_union_kept(zoned_ngff, monkeypatch): + _veto(monkeypatch, _boxes("A", "B")) + seg = GrandQCVetoSegmenter() + with OmeZarrReader(zoned_ngff) as r: + union, _, _, _ = GrandQCVetoSegmenter()._union(r) + out = seg.segment(r).mask > 0 + + assert _frac(union, "C") > 0.90, "the union must actually keep the junk block" + assert _frac(out, "C") == 0.0, "the veto must remove it" + assert _frac(out, "A") > 0.90 and _frac(out, "B") > 0.90 # endorsed zones survive + assert not np.any(out & ~union), "the veto may only remove, never add" + assert seg.last_stats["n_dropped"] == 1 + assert 0.2 < seg.last_stats["veto_removed_frac"] < 0.5 + + +def test_clip_branch_erodes_a_large_mixed_component(zoned_ngff, monkeypatch): + """A big component that is half-endorsed is clipped; a pure component veto is not. + + This is the whole point of the hybrid: boundary erosion is cheap on large + components, and large components are exactly the ones that can fuse real tissue + with an attached artifact. + """ + half = np.zeros((VETO_N, VETO_N), np.float32) + r0, r1, c0, c1 = BOXES["A"] + half[int((r0 - 0.04) * VETO_N) : int((r1 + 0.04) * VETO_N), + int((c0 - 0.04) * VETO_N) : int(((c0 + c1) / 2) * VETO_N)] = 1.0 # left half + r0, r1, c0, c1 = BOXES["B"] + half[int((r0 - 0.04) * VETO_N) : int((r1 + 0.04) * VETO_N), + int((c0 - 0.04) * VETO_N) : int((c1 + 0.04) * VETO_N)] = 1.0 + _veto(monkeypatch, half) + + hybrid, component = GrandQCVetoSegmenter(), GrandQCVetoSegmenter( + clip_min_area_mm2=math.inf # documented corner case: clip never fires + ) + with OmeZarrReader(zoned_ngff) as r: + h = hybrid.segment(r).mask > 0 + c = component.segment(r).mask > 0 + + assert 0.3 < _frac(h, "A") < 0.7, "clipped back to roughly the endorsed half" + assert _frac(c, "A") > 0.90, "component veto keeps the mixed component whole" + assert hybrid.last_stats["n_clipped"] == 1 + assert component.last_stats["n_clipped"] == 0 + assert _frac(h, "C") == 0.0 and _frac(c, "C") == 0.0 # both still drop the junk + + +def test_slack_and_width_guards_protect_a_thin_unendorsed_sliver( + zoned_ngff, monkeypatch +): + """A sliver of disagreement is mis-registration, not an artifact -- keep it.""" + sliver = _boxes("A", "B") + sliver[:, int(0.24 * VETO_N) : int(0.24 * VETO_N) + 2] = 0.0 # ~2 veto px wide + _veto(monkeypatch, sliver) + + with OmeZarrReader(zoned_ngff) as r: + default = GrandQCVetoSegmenter().segment(r).mask > 0 + # keep_cov high enough to force the clip branch in all three variants + no_guards = GrandQCVetoSegmenter( + keep_cov=1.1, veto_slack_um=0.0, removal_min_width_um=0.0 + ).segment(r) + width_only = GrandQCVetoSegmenter(keep_cov=1.1, veto_slack_um=0.0).segment(r) + + assert _frac(default, "A") > 0.95, "veto_slack_um must absorb the sliver" + assert _frac(no_guards.mask > 0, "A") < _frac(default, "A"), ( + "with both guards off the sliver is eroded away -- proves the guards act" + ) + assert _frac(width_only.mask > 0, "A") == _frac(default, "A"), ( + "opening-by-reconstruction alone must also put the sliver back" + ) + + +def test_pure_pixelwise_and_is_exactly_union_and_veto(zoned_ngff, monkeypatch): + """The documented corner case (a) must be a bit-exact pixel-wise AND.""" + g10 = _boxes("A", "B", pad=0.0) + _veto(monkeypatch, g10) + seg = GrandQCVetoSegmenter( + veto_cov=0.0, + keep_cov=1.1, + clip_min_area_mm2=0.0, + veto_slack_um=0.0, + removal_min_width_um=0.0, + ) + with OmeZarrReader(zoned_ngff) as r: + union, _, _, _ = GrandQCVetoSegmenter()._union(r) + out = seg.segment(r).mask > 0 + + assert np.array_equal(out, union & _align(g10 > 0, union.shape)) + assert out.sum() < union.sum() # non-trivial: the AND actually removes something + + +def test_degenerate_veto_emits_the_union_and_warns(zoned_ngff, monkeypatch): + """An empty veto prediction must never silently return zero tiles.""" + _veto(monkeypatch, np.zeros((VETO_N, VETO_N), np.float32)) + seg = GrandQCVetoSegmenter() + with OmeZarrReader(zoned_ngff) as r: + union, _, _, _ = GrandQCVetoSegmenter()._union(r) + with pytest.warns(RuntimeWarning, match="degenerate"): + out = seg.segment(r).mask > 0 + + assert np.array_equal(out, union) + assert seg.last_stats["veto_aborted"] is True + + +def test_veto_failure_raises_by_default_and_can_degrade_to_the_union( + zoned_ngff, monkeypatch +): + def boom(self, reader): + raise RuntimeError("no [grandqc] extra") + + monkeypatch.setattr(GrandQCVetoSegmenter, "_veto_mask", boom) + with OmeZarrReader(zoned_ngff) as r: + with pytest.raises(RuntimeError, match="grandqc"): + GrandQCVetoSegmenter().segment(r) + with pytest.warns(RuntimeWarning, match="artifacts will NOT be removed"): + out = GrandQCVetoSegmenter(on_veto_error="warn").segment(r).mask > 0 + union, _, _, _ = GrandQCVetoSegmenter()._union(r) + assert np.array_equal(out, union) + + +def test_empty_union_short_circuits_without_running_the_veto(blank_ngff, monkeypatch): + calls: list = [] + _veto(monkeypatch, _boxes("A"), counter=calls) + seg = GrandQCVetoSegmenter() + with OmeZarrReader(blank_ngff) as r: + tm = seg.segment(r) + assert tm.mask.sum() == 0 + assert calls == [], "the GPU term must be skipped when there is nothing to veto" + assert seg.last_stats["n_components"] == 0 + + +# -- the scan-region / canvas guard -------------------------------------------- + + +def test_scan_region_kill_finds_the_canvas_not_the_tissue(): + img = _paint_canvas(240, 240) + kill = scan_region_kill(img, 8.0) + h = w = 240 + canvas = np.zeros((h, w), bool) + m = int(CANVAS_MARGIN * h) + canvas[:m] = canvas[-m:] = True + canvas[:, :m] = canvas[:, -m:] = True + tissue = np.zeros((h, w), bool) + r0, r1 = int(TISSUE_BOX[0] * h), int(TISSUE_BOX[1] * h) + tissue[r0:r1, r0:r1] = True + + assert kill[canvas].all(), "the synthetic fill itself must be killed" + assert not kill[tissue].any(), "and it must not reach the specimen" + # No canvas at all -> only the thin perimeter ring, so the guard is a no-op inside. + plain = scan_region_kill(_paint_zones(240, 240), 8.0) + assert plain.mean() < 0.06 + + +def test_guard_rescues_the_canny_term_on_a_canvas_slide(canvas_ngff, monkeypatch): + """The load-bearing property: guard the EDGE MAP, not the finished mask. + + Unguarded, the scan-region rectangle is Canny's outer contour and the real tissue + inside it becomes a ``max_hole_frac`` hole that gets carved back to background -- + so the border does not merely add junk, it deletes the tissue. + """ + with OmeZarrReader(canvas_ngff) as r: + unguarded = CannySegmenter(seg_mpp=8.0).segment(r).mask > 0 + guarded = GrandQCVetoSegmenter(guard=True) + _veto(monkeypatch, np.ones((VETO_N, VETO_N), np.float32)) + out = guarded.segment(r).mask > 0 + union_off, _, _, _ = GrandQCVetoSegmenter(guard=False)._union(r) + + h = w = out.shape[0] + canvas = np.zeros((h, w), bool) + m = int(CANVAS_MARGIN * h) + canvas[:m] = canvas[-m:] = True + canvas[:, :m] = canvas[:, -m:] = True + tissue = np.zeros((h, w), bool) + r0, r1 = int(TISSUE_BOX[0] * h), int(TISSUE_BOX[1] * h) + tissue[r0:r1, r0:r1] = True + + assert unguarded[tissue].mean() == 0.0, "unguarded canny loses the tissue entirely" + assert union_off[canvas].mean() > 0.05, "and keeps a bogus border ring" + assert out[canvas].mean() == 0.0, "the guard removes the canvas and its ring" + assert out[tissue].mean() > 0.90, "while the specimen survives" + + +def test_canny_edge_filter_hook_sees_the_raw_edge_map(zoned_ngff): + """The hook must run on the edges, with the level's own image and MPP.""" + seen: list = [] + + def kill_all(edges, img, level_mpp): + seen.append((edges.shape, img.shape[:2], level_mpp, int((edges > 0).sum()))) + return np.zeros_like(edges) + + with OmeZarrReader(zoned_ngff) as r: + plain = CannySegmenter(seg_mpp=8.0).segment(r).mask + blanked = CannySegmenter(seg_mpp=8.0, edge_filter=kill_all).segment(r).mask + + assert plain.sum() > 0 and blanked.sum() == 0 + (eshape, ishape, mpp, n_edges) = seen[0] + assert eshape == ishape == plain.shape + assert mpp == pytest.approx(8.0) # level MPP, so a filter can size itself in µm + assert n_edges > 0, "the hook must see the edge map, not an already-filled mask" + + +def test_align_or_pools_on_the_downsample_path(): + """INTER_NEAREST on the downsample path deletes thin structures; ours must not.""" + m = np.zeros((64, 64), bool) + m[:, 30] = True # a 1-px wide line + small = _align(m, (16, 16)) + assert small.any(), "a thin structure must survive an 4x downsample" + big = _align(m, (128, 128)) + assert big.sum() > m.sum() diff --git a/tests/test_segmenter_guarded.py b/tests/test_segmenter_guarded.py new file mode 100644 index 0000000..c7d8710 --- /dev/null +++ b/tests/test_segmenter_guarded.py @@ -0,0 +1,134 @@ +"""Tests for the guarded-union segmenter (``otsu | canny`` minus the scanner canvas). + +The shared ``synthetic_ngff`` fixture is a smooth ramp with no edges and no canvas, so +assertions about guard behaviour would pass vacuously on it. Every behavioural test here +uses ``canvas_ngff``: a **dark, saturated** canvas block (saturated so Otsu genuinely +claims it -- a pure-black canvas is invisible to a saturation threshold and would let a +broken guard pass) with a separate pale specimen inside the scanned region. +""" + +from __future__ import annotations + +import numpy as np +import pytest +import zarr + +from conftest import build_ngff_v04 +from raw2features.readers.omezarr import OmeZarrReader +from raw2features.segmenters.guarded import GuardedUnionSegmenter +from raw2features.segmenters.otsu import OtsuSegmenter + +SIZES = ((200, 300), (100, 150), (50, 75)) + +# The fixture slide is 200x300 px at 0.5 um/px = 100x150 um TOTAL -- smaller than the +# production guard's 200 um border band, which would kill the whole image. Scale the +# guard to the fixture and segment at the finest level so Canny has pixels to work with. +SEG_MPP = 0.5 +GUARD = {"border_um": 4.0, "frame_um": 1.0} + + +@pytest.fixture +def canvas_ngff(tmp_path) -> str: + """Dark saturated canvas on the left; a pale textured specimen on the right.""" + path = build_ngff_v04(str(tmp_path / "canvas.zarr")) + g = zarr.open_group(path, mode="r+", zarr_format=2) + rng = np.random.default_rng(0) + for i, (h, w) in enumerate(SIZES): + # light background + r = np.full((h, w), 235, np.uint8) + gr = np.full((h, w), 232, np.uint8) + b = np.full((h, w), 238, np.uint8) + # canvas: dark AND saturated (deep blue) so Otsu claims it + cw = w // 3 + r[:, :cw], gr[:, :cw], b[:, :cw] = 8, 10, 70 + # specimen: pale pink with texture, well inside the scanned region + y0, y1 = h // 3, 2 * h // 3 + x0, x1 = w // 2, w - w // 8 + tex = rng.integers(0, 26, size=(y1 - y0, x1 - x0), dtype=np.uint8) + r[y0:y1, x0:x1] = 214 + (tex // 6) + gr[y0:y1, x0:x1] = 176 + tex + b[y0:y1, x0:x1] = 202 + (tex // 4) + a = g[str(i)] + a[0, 0, 0], a[0, 1, 0], a[0, 2, 0] = r, gr, b + return path + + +def _mask(seg, path): + with OmeZarrReader(path) as r: + return seg.segment(r) + + +def test_guarded_returns_valid_binary_tissuemask(canvas_ngff): + tm = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, **GUARD), canvas_ngff) + assert tm.mask.dtype == np.float32 + assert set(np.unique(tm.mask)).issubset({0.0, 1.0}) + assert tm.mask.ndim == 2 + assert tm.downsample > 0 + + +def test_guard_removes_the_canvas_that_otsu_claims(canvas_ngff): + """The load-bearing behaviour: Otsu keeps the dark canvas, the guard removes it.""" + otsu = _mask(OtsuSegmenter(seg_mpp=SEG_MPP), canvas_ngff).mask > 0 + guarded = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, **GUARD), canvas_ngff).mask > 0 + + cw = otsu.shape[1] // 3 + otsu_canvas = otsu[:, :cw].mean() + guarded_canvas = guarded[:, :cw].mean() + + # Otsu must actually claim the canvas, or this fixture proves nothing. + assert otsu_canvas > 0.5, "fixture too weak: otsu did not claim the canvas" + assert guarded_canvas < 0.05, "guard failed to remove the canvas" + + +def test_guard_off_reinstates_the_canvas(canvas_ngff): + """guard=False must be a real ablation, not a no-op.""" + on = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, guard=True, **GUARD), canvas_ngff).mask > 0 + off = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, guard=False, **GUARD), canvas_ngff).mask > 0 + cw = on.shape[1] // 3 + assert off[:, :cw].mean() > on[:, :cw].mean() + 0.4 + + +def test_guarded_recovers_specimen_otsu_misses_entirely(canvas_ngff): + """The union half, asserted as a *delta* rather than a tuned absolute level. + + This fixture reproduces the production failure exactly: Otsu claims 100% of the + saturated canvas and **0%** of the pale specimen (the canvas drags its global + threshold up). The guarded union must invert both. Asserting ``otsu == 0`` keeps + the test honest -- if a change ever makes Otsu find the specimen on its own, this + fails loudly rather than passing for the wrong reason. + """ + otsu = _mask(OtsuSegmenter(seg_mpp=SEG_MPP), canvas_ngff).mask > 0 + guarded = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, **GUARD), canvas_ngff).mask > 0 + h, w = guarded.shape + sy, sx = slice(h // 3, 2 * h // 3), slice(w // 2, w - w // 8) + + assert otsu[sy, sx].mean() == 0.0, "fixture too weak: otsu already finds the specimen" + assert guarded[sy, sx].mean() > 0.05, "union failed to recover the specimen" + # ...and it is the Canny term, not a relaxed guard, that supplies it. + assert guarded[sy, sx].sum() > otsu[sy, sx].sum() + + +def test_guarded_needs_no_grandqc_weights(canvas_ngff, monkeypatch): + """It must never import/instantiate the GrandQC model (no GPU, no download).""" + import raw2features.qc.grandqc as gq + + def _boom(*a, **k): # pragma: no cover - fails the test if reached + raise AssertionError("guarded must not touch GrandQC") + + monkeypatch.setattr(gq, "GrandQC", _boom) + tm = _mask(GuardedUnionSegmenter(seg_mpp=SEG_MPP, **GUARD), canvas_ngff) + assert tm.mask.any() + + +def test_guarded_records_provenance(canvas_ngff): + seg = GuardedUnionSegmenter(seg_mpp=SEG_MPP, **GUARD) + _mask(seg, canvas_ngff) + assert "kill_frac" in seg.last_stats + assert seg.last_stats["veto"] == "none" + + +def test_guarded_is_registered_and_no_arg_constructible(): + from raw2features.core import plugins + + seg = plugins.get("segmenters", "guarded")() + assert seg.name == "guarded"