wip

docs/SEGMENTATION.md

@@ -12,7 +12,10 @@ 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`) |
+| `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)
 
@@ -37,6 +40,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
@@ -106,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.

pyproject.toml

@@ -59,9 +59,12 @@ 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"
+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"

src/raw2features/qc/grandqc.py

@@ -88,14 +88,29 @@ 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)``.
 
     The architecture (UNet++ tissue / UNet artifact) and class count are read from the
     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 ``<name>.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

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(

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)

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)