diff --git a/docs/api.md b/docs/api.md index 690f7713a..a7da777a3 100644 --- a/docs/api.md +++ b/docs/api.md @@ -155,6 +155,7 @@ See the {doc}`extensibility guide ` for how to implement a custo experimental.im.calculate_image_features experimental.tl.calculate_tiling_qc experimental.tl.TilingQCParams + experimental.tl.SeamDetectionParams experimental.tl.assign_stitch_groups experimental.tl.StitchParams experimental.tl.make_stitched_labels diff --git a/src/squidpy/experimental/im/_tiling.py b/src/squidpy/experimental/im/_tiling.py index d23720990..f5371d9f1 100644 --- a/src/squidpy/experimental/im/_tiling.py +++ b/src/squidpy/experimental/im/_tiling.py @@ -351,12 +351,42 @@ def extract_labels_tile_lazy( ------- ``(crop_h, crop_w)`` numpy array with non-owned cells zeroed. """ - cy0, cx0, cy1, cx1 = spec.crop - tile_labels = _as_2d(_materialize(labels_da.isel(y=slice(cy0, cy1), x=slice(cx0, cx1))).copy()) + tile_labels = _crop_labels(labels_da, spec) _zero_non_owned(tile_labels, spec.owned_ids) return tile_labels +def _crop_labels(labels_da: xr.DataArray, spec: TileSpec) -> np.ndarray: + """Materialize a tile's crop region as a writable 2-D label array, before any masking.""" + cy0, cx0, cy1, cx1 = spec.crop + return _as_2d(_materialize(labels_da.isel(y=slice(cy0, cy1), x=slice(cx0, cx1))).copy()) + + +def extract_labels_tile_with_occupancy( + labels_da: xr.DataArray, + spec: TileSpec, +) -> tuple[np.ndarray, np.ndarray]: + """Extract a labels tile together with the *unmasked* occupancy of the same crop. + + ``tile_labels`` has non-owned cells zeroed, so every cell is scored by exactly one + tile. ``occupancy`` marks every labelled pixel in the crop, owned or not. + + Geometry that asks "is there tissue beyond this edge?" must use ``occupancy``: on the + masked array a cell whose neighbour is owned by the adjacent tile reads that neighbour + as background, which turns an ordinary inter-cell membrane into an apparent wide gap + along every tile border. Only one crop is materialized for both outputs. + + Returns + ------- + ``(tile_labels, occupancy)`` -- the owned-only label crop and the boolean occupancy of + the same region before masking. + """ + tile_labels = _crop_labels(labels_da, spec) + occupancy = tile_labels != 0 + _zero_non_owned(tile_labels, spec.owned_ids) + return tile_labels, occupancy + + def _zero_non_owned(tile_labels: np.ndarray, owned_ids: frozenset[int]) -> None: """Zero out labels not in *owned_ids* (in-place). diff --git a/src/squidpy/experimental/pl/_tiling_qc.py b/src/squidpy/experimental/pl/_tiling_qc.py index 9844af79e..c1f59b5be 100644 --- a/src/squidpy/experimental/pl/_tiling_qc.py +++ b/src/squidpy/experimental/pl/_tiling_qc.py @@ -20,6 +20,8 @@ def tiling_qc( "max_straight_edge_ratio", "cardinal_alignment_score", "is_outlier", + "is_seam_cut", + "seam_dist", ] = "nhood_outlier_fraction", cmap: str = "RdYlGn_r", figsize: tuple[float, float] | None = None, @@ -44,7 +46,9 @@ def tiling_qc( Which ``.obs`` column to colour by. One of ``"nhood_outlier_fraction"``, ``"smoothed_cut_score"``, ``"cut_score"``, ``"max_straight_edge_ratio"``, - ``"cardinal_alignment_score"``, ``"is_outlier"``. + ``"cardinal_alignment_score"``, ``"is_outlier"`` (MAD gate), + or -- from ``calculate_tiling_qc(detect_seams=True)`` -- the + emergent-seam columns ``"is_seam_cut"`` and ``"seam_dist"``. cmap Matplotlib colormap name. figsize @@ -70,9 +74,11 @@ def tiling_qc( "nhood_outlier_fraction": "Neighborhood outlier fraction", "smoothed_cut_score": "Smoothed cut score", "cut_score": "Cut score", - "is_outlier": "Outlier flag", + "is_outlier": "Outlier flag (MAD)", "max_straight_edge_ratio": "Max straight edge ratio", "cardinal_alignment_score": "Cardinal alignment score", + "is_seam_cut": "Seam-cut flag (emergent seam)", + "seam_dist": "Distance to seam", } show_kwargs: dict[str, object] = {"title": _TITLES.get(score_col, score_col)} diff --git a/src/squidpy/experimental/tl/__init__.py b/src/squidpy/experimental/tl/__init__.py index bd14f9e54..4ff3da5d4 100644 --- a/src/squidpy/experimental/tl/__init__.py +++ b/src/squidpy/experimental/tl/__init__.py @@ -1,7 +1,15 @@ from __future__ import annotations +from ._seam import SeamDetectionParams from ._stitched_labels import make_stitched_labels from ._tiling_qc import TilingQCParams, calculate_tiling_qc from ._tiling_stitch import StitchParams, assign_stitch_groups -__all__ = ["StitchParams", "TilingQCParams", "assign_stitch_groups", "calculate_tiling_qc", "make_stitched_labels"] +__all__ = [ + "SeamDetectionParams", + "StitchParams", + "TilingQCParams", + "assign_stitch_groups", + "calculate_tiling_qc", + "make_stitched_labels", +] diff --git a/src/squidpy/experimental/tl/_seam.py b/src/squidpy/experimental/tl/_seam.py new file mode 100644 index 000000000..8ea5ca941 --- /dev/null +++ b/src/squidpy/experimental/tl/_seam.py @@ -0,0 +1,457 @@ +"""Emergent-seam detection for tile-boundary cut cells. + +The MAD-based ``is_outlier`` gate in :func:`~squidpy.experimental.tl.calculate_tiling_qc` +flags cells whose boundary is unusually straight. In dense tissue with a wide inter-FOV +gap -- and when a cut leaves only one segmentable half -- that signal is swamped: most real +cut cells are missed and interior cells are flagged instead. + +This module adds a complementary, geometry-only detector that needs neither the FOV size +nor tile overlap, and works on single-sided cuts: + +1. For every cell, collect **all** long, cardinal (axis-aligned) flat boundary runs, together + with the depth of background lying just beyond each one. +2. Seam lines **emerge** as the coordinates where many such runs align -- a consensus over + many cells -- with a wide inter-FOV gap spreading a seam into a band. A peak counts as a + seam when its aligned-edge count is too large to arise from edges scattered uniformly + along the axis, so a weakly-populated seam is judged on its own evidence rather than + against the strongest peak present. +3. A cell is flagged ``is_seam_cut`` iff any of its cardinal edges lies in a detected seam + band and faces it. Membership does not require a wide gap: once the seam's location is + known, a close-gap or single-sided cut on that line counts too. + +**Which edges vote** depends on what the data supports. In packed tissue an ordinary facet +faces a neighbour one membrane away, so a wide background gap is rare and isolates the seam; +in sparse tissue nearly every edge faces open background and the gap says nothing, so all +edges vote and the alignment consensus carries detection alone. The split between the two +gap classes and the choice of channel are both read off the observed gap distribution +(:func:`gap_channel`), never set by hand. + +**No absolute-pixel thresholds.** Every length scale is expressed as a ratio and resolved +at runtime against the data's own length scale ``D`` (the median cell equivalent diameter) +and the observed gap distribution, so the same defaults transfer across resolutions, cell +sizes and FOV pitches. See :class:`SeamDetectionParams`. + +**Gaps must be probed on unmasked pixels.** The caller measures background beyond an edge on +an occupancy array that includes cells belonging to *other* processing tiles. Probing a +tile-masked array instead makes every neighbour dropped by masking look like open background, +which manufactures a seam along each processing tile border -- a detection result that would +depend on ``tile_size`` rather than on the data. +""" + +from __future__ import annotations + +from dataclasses import dataclass, fields +from statistics import median +from typing import Any + +import numpy as np +from scipy.signal import find_peaks +from scipy.stats import poisson +from skimage.filters import threshold_otsu + + +@dataclass(slots=True, frozen=True) +class SeamDetectionParams: + """Scale-invariant tuning knobs for emergent-seam cut-cell detection. + + No knob is an absolute pixel length. Each is either a **dimensionless ratio** resolved at + runtime against the data's length scale ``D`` (median cell equivalent diameter), a + **statistical criterion** (:attr:`alpha`, :attr:`gap_selectivity_max`) evaluated against + the data's own distributions, or a **pixel-grid constant** (:attr:`flat_tol`, + :attr:`bin_width`) that describes the raster rather than the tissue -- so the defaults + transfer across resolutions, cell sizes and FOV pitches. + + The thresholds that matter most are derived rather than set: the split between membrane + gaps and open background comes from the observed gap distribution, and the height a + histogram peak must reach comes from a null built on the edge count and axis extent. + """ + + edge_len_frac: float = 0.25 + """Minimum cardinal flat-run length, as a fraction of the median cell diameter ``D``. + + A cut rarely bisects a cell through its widest point, so requiring the straight edge to + span half the cell (the previous default) discards most real cuts. Measured on real + CosMx breast tissue against seam-adjacency ground truth, lowering this from 0.5 to 0.25 + raises recall 0.38 -> 0.54 at precision 0.98 -> 0.89 (F1 0.54 -> 0.67); on both synthetic + fixtures F1 moves by less than 0.01 either way, and the no-seam control stays empty at + every value. Raising it trades recall for precision.""" + + flat_tol: float = 1.5 + """Maximum deviation in pixels from a single coordinate for a run to count as flat. + A rasterisation / pixel-grid constant (~1px), independent of cell size or resolution -- + scaling it with the cell size would accept curved edges on large cells.""" + + probe_frac: float = 0.6 + """How far to probe for background beyond an edge, as a fraction of ``D``.""" + + gap_selectivity_max: float = 0.25 + """Use the wide-gap edge channel only while at most this fraction of edges passes the + (data-derived) gap split -- i.e. only while the split actually discriminates. + + In packed tissue an ordinary facet faces a neighbour one membrane away, so a wide gap is + rare and the split isolates the seam: measured 3-5% on real CosMx breast tissue and 5.5% + on the dense synthetic. In sparse tissue almost every edge faces open background (45-48% + measured), the gap says nothing about seams, and detection falls back to all edges. The + decision boundary sits an order of magnitude from either regime.""" + + bin_width: float = 2.0 + """Seam histogram bin width in pixels -- a raster resolution constant, independent of cell size.""" + + alpha: float = 0.01 + """Significance level for keeping a histogram peak, Bonferroni-corrected over the bins. + + A peak is a seam when its aligned-edge count is too large to come from edges scattered + uniformly along the axis. The null is computed from the data (edge count and axis + extent), so this transfers across datasets in a way a hand-set height ratio does not, and + a weakly-populated seam is judged on its own significance rather than against the + strongest peak on its axis.""" + + cluster_frac: float = 1.4 + """Histogram peaks within this multiple of ``D`` merge into one seam band (spans the gap).""" + + band_margin_frac: float = 0.15 + """Extra seam-band half-width, as a fraction of ``D``.""" + + flag_tol_frac: float = 0.25 + """Slack added to the band half-width when testing edge membership, as a fraction of ``D``.""" + + face_slack_frac: float = 0.15 + """Slack allowing an edge slightly past the band centre to still face the seam, fraction of ``D``.""" + + def __post_init__(self) -> None: + # Coerce every field (accepts numpy scalars cleanly); iterating the fields rather than a + # hand-written name list means a new knob cannot silently skip coercion. + for f in fields(self): + object.__setattr__(self, f.name, float(getattr(self, f.name))) + if self.edge_len_frac <= 0: + raise ValueError(f"edge_len_frac must be > 0, got {self.edge_len_frac}.") + if not 0.0 <= self.gap_selectivity_max <= 1.0: + raise ValueError(f"gap_selectivity_max must be in [0, 1], got {self.gap_selectivity_max}.") + if not 0.0 < self.alpha < 1.0: + raise ValueError(f"alpha must be in (0, 1), got {self.alpha}.") + if self.probe_frac <= 0: + raise ValueError(f"probe_frac must be > 0, got {self.probe_frac}.") + + def _resolve(self, diameter: float) -> SeamScale: + """Resolve every fraction against the data's length scale ``D`` (median cell diameter).""" + return SeamScale( + params=self, + diameter=float(diameter), + min_len=max(3, int(round(self.edge_len_frac * diameter))), + probe_depth=max(3, int(round(self.probe_frac * diameter))), + flag_tol=self.flag_tol_frac * diameter, + face_slack=self.face_slack_frac * diameter, + cluster_gap=self.cluster_frac * diameter, + band_margin=self.band_margin_frac * diameter, + ) + + +_SEAM_DEFAULTS = SeamDetectionParams() + + +@dataclass(frozen=True, slots=True) +class SeamScale: + """:class:`SeamDetectionParams` resolved against the data's length scale ``D`` (pixels). + + The fractions are dimensionless by design, but every consumer needs them in pixels. + Resolving once and passing this object keeps the arithmetic -- and the ``max(3, ...)`` + floors -- in a single place, so seam *detection* and the stitcher's edge *extraction* + cannot resolve the same knob to two different pixel values. + """ + + params: SeamDetectionParams + diameter: float + min_len: int + probe_depth: int + flag_tol: float + face_slack: float + cluster_gap: float + band_margin: float + + @property + def flat_tol(self) -> float: + """Flatness tolerance (px) -- a raster constant, so it is not scaled by ``D``.""" + return self.params.flat_tol + + +def _dominant_flat_line(extreme: np.ndarray, present: np.ndarray, flat_tol: float) -> tuple[float, int, int]: + """Longest flat run on one side of a cell, at *any* coordinate (not just the extreme). + + A cut leaves a straight boundary that is a plateau in the per-index extreme coordinate -- + but that plateau need not be at the cell's outermost point (a wide cell can have a partial + cut plus other geometry). Every distinct extreme value is a candidate line; the one whose + within-``flat_tol`` run of consecutive present indices is longest wins, returned as + ``(coord, start, length)``. + + All candidates are tested in one pass: the (candidate x index) membership matrix is + flattened with a False separator column, so a single run-length scan over the flat array + finds the longest run of any candidate. This sits in the innermost loop of the seam pass + -- four sides for every cell in the dataset -- so the per-candidate Python scan it replaces + was a measurable share of the runtime. + """ + cand = np.unique(np.round(extreme[present])) + if cand.size == 0: + return 0.0, 0, 0 + on = (np.abs(extreme[None, :] - cand[:, None]) <= flat_tol) & present[None, :] + padded = np.zeros((on.shape[0], on.shape[1] + 1), dtype=np.int8) + padded[:, :-1] = on + diffs = np.diff(np.concatenate(([np.int8(0)], padded.ravel()))) + starts = np.flatnonzero(diffs == 1) + if starts.size == 0: + return 0.0, 0, 0 + lengths = np.flatnonzero(diffs == -1) - starts + best = int(lengths.argmax()) # ties resolve to the lowest candidate, then leftmost run + row, col = divmod(int(starts[best]), padded.shape[1]) + return float(cand[row]), col, int(lengths[best]) + + +def _probe_gap( + occupancy: np.ndarray, axis: str, coord: int, side: int, run_lo: int, run_hi: int, probe_depth: int +) -> float: + """Median count of consecutive background pixels just outside a flat edge (tile-local coords). + + ``occupancy`` must mark **every** cell in the neighbourhood, not only the ones owned by the + current processing tile -- otherwise a neighbour dropped by tile masking reads as background + and an ordinary membrane is mistaken for a seam gap. + """ + height, width = occupancy.shape + step = 1 if side == -1 else -1 + # .tolist() keeps the sampled positions identical while making each `occupancy[y, x]` + # a fast int index rather than a numpy-scalar one. + idxs = np.linspace(run_lo, run_hi - 1, min(9, max(1, run_hi - run_lo))).astype(int).tolist() + depths = [] + for t in idxs: + d = 0 + for k in range(1, probe_depth + 1): + y, x = (t, coord + step * k) if axis == "v" else (coord + step * k, t) + if not (0 <= y < height and 0 <= x < width) or occupancy[y, x]: + break + d += 1 + depths.append(d) + return float(median(depths)) if depths else 0.0 + + +def cell_flat_edges( + mask: np.ndarray, + occupancy: np.ndarray, + bbox: tuple[int, int, int, int], + origin: tuple[int, int], + min_len: int, + flat_tol: float, + probe_depth: int, +) -> list[dict[str, Any]]: + """All cardinal flat boundary runs of one cell (multiple per cell; global coordinates). + + ``occupancy`` is the boolean occupancy of the surrounding crop **including cells not owned + by the current tile** -- see :func:`_probe_gap`. + + Returns a list of ``{"axis": "v"|"h", "coord": float, "span": int, "side": +1|-1, "gap": float}`` + -- one per side (right/left/bottom/top) whose longest flat run reaches ``min_len``. ``gap`` is + the background depth just beyond the edge (used later to separate seam cuts from touching facets). + Length filtering against the data scale and gap thresholding are applied by the caller. + """ + y0, x0, _, _ = bbox + oy, ox = origin + height, width = mask.shape + rows = mask.any(1) + cols = mask.any(0) + rightmost = np.where(rows, width - 1 - mask[:, ::-1].argmax(1), np.nan).astype(float) + leftmost = np.where(rows, mask.argmax(1), np.nan).astype(float) + bottommost = np.where(cols, height - 1 - mask[::-1, :].argmax(0), np.nan).astype(float) + topmost = np.where(cols, mask.argmax(0), np.nan).astype(float) + + out: list[dict[str, Any]] = [] + for axis, extreme, present, side in [ + ("v", rightmost, rows, -1), + ("v", leftmost, rows, +1), + ("h", bottommost, cols, -1), + ("h", topmost, cols, +1), + ]: + if int(present.sum()) < min_len: + continue + c, s, ln = _dominant_flat_line(extreme, present, flat_tol) + if ln < min_len: + continue + c_round = int(round(c)) + if axis == "v": # perp = x (column); parallel/along-seam = y (row) + perp_local, run_lo, run_hi = x0 + c_round, y0 + s, y0 + s + ln + coord_global, ext_lo, ext_hi = ox + perp_local, oy + run_lo, oy + run_hi + else: # perp = y (row); parallel/along-seam = x (column) + perp_local, run_lo, run_hi = y0 + c_round, x0 + s, x0 + s + ln + coord_global, ext_lo, ext_hi = oy + perp_local, ox + run_lo, ox + run_hi + gap = _probe_gap(occupancy, axis, perp_local, side, run_lo, run_hi, probe_depth) + out.append( + { + "axis": axis, + "coord": float(coord_global), + "span": int(ln), + "side": int(side), + "gap": gap, + "extent": (float(ext_lo), float(ext_hi)), + } + ) + return out + + +def _otsu_gap_split(gaps: np.ndarray) -> float: + """Otsu split of the observed edge-gap distribution (returns the threshold in pixels). + + The two classes are physical: an edge either faces a neighbour across a thin membrane, or + it faces open background. Taking the split from the data replaces a hand-set multiple of + an estimated membrane width, and adapts to imaging resolution and segmentation style. + + Gaps are small integer pixel depths, which skimage histograms one-bin-per-value, so there + is no bin count to tune. The ``+ 0.5`` puts the threshold between the two integer classes, + matching the ``>=`` test in :func:`gap_channel`. + """ + g = np.rint(gaps[np.isfinite(gaps)]).astype(int) + if g.size == 0: + return float("inf") + if np.unique(g).size < 2: + # A single distinct gap has no split; treat every edge as below threshold so the + # caller falls back to the all-edges channel rather than filtering on noise. + return float(g[0] + 1) + return float(threshold_otsu(g)) + 0.5 + + +def gap_channel(edges: list[dict[str, Any]], params: SeamDetectionParams = _SEAM_DEFAULTS) -> tuple[float, float, bool]: + """Choose the edge channel the data supports. + + Returns ``(threshold, selectivity, use_wide_gap)`` -- the data-derived gap split, the + fraction of edges above it, and whether that split is selective enough to filter on. + See :attr:`SeamDetectionParams.gap_selectivity_max`. + """ + gaps = np.array([e["gap"] for e in edges], dtype=float) + if gaps.size == 0: + return float("inf"), 0.0, False + thresh = _otsu_gap_split(gaps) + selectivity = float(np.mean(gaps >= thresh)) + return thresh, selectivity, selectivity <= params.gap_selectivity_max + + +def _min_significant_count(n_edges: int, n_bins: int, window: int, alpha: float) -> float: + """Smallest smoothed bin count not explainable by uniformly scattered edges. + + Under the null the edges are spread uniformly along the axis, so the count in a + ``2*window+1`` bin neighbourhood is Poisson with mean ``n_edges*(2*window+1)/n_bins``. + The threshold is that distribution's upper tail at ``alpha``, Bonferroni-corrected over + the bins tested. Both inputs come from the data, so the test transfers across datasets + where a fixed height ratio would not. + """ + if n_edges == 0 or n_bins == 0: + return float("inf") + lam = n_edges * (2 * window + 1) / n_bins + return float(max(3.0, poisson.isf(alpha / n_bins, lam) + 1)) + + +def detect_seams( + edges: list[dict[str, Any]], + extent_x: int, + extent_y: int, + scale: SeamScale, + channel: tuple[float, float, bool] | None = None, +) -> dict[str, list[tuple[float, float, int]]]: + """Locate seam bands per axis from the aligned edges the data supports. + + Two channels carry the seam signal in different tissue regimes, and the selectivity of the + gap split decides which one is used (:func:`gap_channel`): + + * **wide-gap edges** -- in packed tissue an ordinary facet faces a neighbour one membrane + away, so a wide gap is rare and isolates the seam cleanly. + * **all edges** -- in sparse tissue nearly every edge faces open background, so the gap + carries no information about seams and the alignment consensus carries detection alone. + + A peak becomes a seam when its aligned-edge count is too large to arise from edges + scattered uniformly along the axis (:func:`_min_significant_count`), so a weakly-populated + seam is judged on its own significance rather than against the strongest peak on its axis. + + Pass ``channel`` when the caller has already run :func:`gap_channel` (it is also useful as + a diagnostic), to avoid a second pass over every edge in the dataset. + + Returns ``{"v": [(centre, half_width, count), ...], "h": [...]}``. + """ + d = scale.diameter + params = scale.params + bin_w = params.bin_width + gap_thresh, _selectivity, use_wide_gap = channel if channel is not None else gap_channel(edges, params) + + out: dict[str, list[tuple[float, float, int]]] = {} + for axis, extent in (("v", extent_x), ("h", extent_y)): + coords = np.array( + [e["coord"] for e in edges if e["axis"] == axis and (not use_wide_gap or e["gap"] >= gap_thresh)] + ) + if coords.size == 0: + out[axis] = [] + continue + bins = np.arange(0, extent + bin_w, bin_w) + hist, _ = np.histogram(coords, bins=bins) + # Smooth only for sub-pixel wobble (small window); the CLUSTER step below is what + # widens a seam across the inter-FOV gap -- over-smoothing here collapses the gap's + # twin peaks into one narrow band and loses the far-side cells. + win = max(1, int(round((0.12 * d) / bin_w))) + smooth = np.convolve(hist, np.ones(2 * win + 1), mode="same") + height = _min_significant_count(coords.size, hist.size, win, params.alpha) + peaks, _ = find_peaks(smooth, height=height, distance=1) + if peaks.size == 0: + out[axis] = [] + continue + centres = bins[peaks] + bin_w / 2 + counts = smooth[peaks].astype(float) + order = np.argsort(centres) + centres, counts = centres[order], counts[order] + # Peaks closer together than `cluster_gap` belong to one seam: a wide inter-FOV gap + # puts one peak on each of its sides. + splits = np.flatnonzero(np.diff(centres) > scale.cluster_gap) + 1 + bands = [] + for cl in np.split(np.arange(centres.size), splits): + cc, ww = centres[cl], counts[cl] + bands.append( + ( + float(np.average(cc, weights=ww)), + float((cc.max() - cc.min()) / 2 + scale.band_margin), + int(ww.sum()), + ) + ) + out[axis] = bands + return out + + +def seam_offset( + edge: dict[str, Any], + bands: list[tuple[float, float, int]], + scale: SeamScale, +) -> float | None: + """Distance from ``edge`` to the seam band it lies on and faces, or ``None`` if it does neither. + + This is the definition of "this edge is a seam cut", and both stages must apply it + identically: :func:`flag_cells_on_seams` decides which cells `calculate_tiling_qc` marks, + and the stitcher decides which edges may be paired. Sharing one predicate is what makes + those two agree -- when it lived in both modules the copies drifted apart silently. + """ + best: float | None = None + for centre, half, _cnt in bands: + signed = centre - edge["coord"] + if abs(signed) > half + scale.flag_tol: + continue + faces = (signed >= -scale.face_slack) if edge["side"] == -1 else (signed <= scale.face_slack) + if faces and (best is None or abs(signed) < best): + best = abs(signed) + return best + + +def flag_cells_on_seams( + edges_by_cell: dict[int, list[dict[str, Any]]], + seams: dict[str, list[tuple[float, float, int]]], + scale: SeamScale, +) -> dict[int, float]: + """Return ``{cell_id: seam_dist}`` for every cell with a cardinal edge on and facing a seam. + + Once a seam is detected, membership does **not** require a wide gap -- a genuine two-sided + cut whose other half sits close still counts, because its edge lies on the seam consensus. + """ + flagged: dict[int, float] = {} + for cid, edges in edges_by_cell.items(): + dists = [d for e in edges if (d := seam_offset(e, seams[e["axis"]], scale)) is not None] + if dists: + flagged[cid] = min(dists) + return flagged diff --git a/src/squidpy/experimental/tl/_tiling_qc.py b/src/squidpy/experimental/tl/_tiling_qc.py index a8feb43c9..5feea8405 100644 --- a/src/squidpy/experimental/tl/_tiling_qc.py +++ b/src/squidpy/experimental/tl/_tiling_qc.py @@ -49,12 +49,21 @@ compute_cell_info_multiscale, compute_cell_info_tiled, extract_labels_tile_lazy, + extract_labels_tile_with_occupancy, ) +from squidpy.experimental.tl._seam import ( + SeamDetectionParams, + SeamScale, + cell_flat_edges, + flag_cells_on_seams, + gap_channel, +) +from squidpy.experimental.tl._seam import detect_seams as _detect_seam_bands from squidpy.experimental.tl._tiling_stitch import _STITCH_COLUMNS, _STITCH_PARAM_KEYS, StitchParams from squidpy.experimental.utils._labels import resolve_labels_array from squidpy.experimental.utils._params import resolve_params -__all__ = ["TilingQCParams", "calculate_tiling_qc"] +__all__ = ["SeamDetectionParams", "TilingQCParams", "calculate_tiling_qc"] @dataclass(slots=True, frozen=True) @@ -334,13 +343,22 @@ def _score_tile( min_area: int = _QC_DEFAULTS.min_area, downsample: int = 1, max_contour_points: int = _QC_DEFAULTS.max_contour_points, -) -> pd.DataFrame: + origin: tuple[int, int] = (0, 0), + seam_scale: SeamScale | None = None, + probe_occupancy: np.ndarray | None = None, +) -> tuple[pd.DataFrame, list[dict[str, Any]]]: """Compute tiling QC metrics for all cells in a numpy label tile. Parameters ---------- tile_labels ``(H, W)`` label array (background = 0, owned cells only). + seam_scale + Resolved seam thresholds. ``None`` skips seam-edge collection entirely. + probe_occupancy + ``(H, W)`` boolean occupancy of the same crop *before* non-owned cells were + masked out. Required alongside ``seam_scale``: seam-gap probing must see + neighbours owned by adjacent tiles, or every tile border looks like a seam. distance_tol Perpendicular distance tolerance for collinearity (pixels). min_area @@ -359,9 +377,12 @@ def _score_tile( """ regions = regionprops(tile_labels) if not regions: - return pd.DataFrame(columns=_TILE_SCORE_COLUMNS, dtype=float) + return pd.DataFrame(columns=_TILE_SCORE_COLUMNS, dtype=float), [] + if seam_scale is not None and probe_occupancy is None: + raise ValueError("seam_scale requires probe_occupancy (the unmasked crop occupancy).") rows: dict[int, dict[str, float]] = {} + edges: list[dict[str, Any]] = [] for region in regions: lid = region.label @@ -394,8 +415,22 @@ def _score_tile( "cardinal_alignment_score": cas, "cut_score": cs, } - - return pd.DataFrame.from_dict(rows, orient="index") + if seam_scale is not None: + # All cardinal flat boundary runs of this cell (full-resolution mask/tile, so + # independent of `downsample`). Seam detection + flagging happen globally later. + for e in cell_flat_edges( + region.image, + probe_occupancy, + region.bbox, + origin, + seam_scale.min_len, + seam_scale.flat_tol, + seam_scale.probe_depth, + ): + e["cell_id"] = lid + edges.append(e) + + return pd.DataFrame.from_dict(rows, orient="index"), edges # Centroid computation (shared logic with _feature.py) @@ -442,6 +477,8 @@ def calculate_tiling_qc( nmads_smoothed: float = 3, n_neighbors: int = 10, tiling_qc_params: TilingQCParams | Mapping[str, Any] | None = None, + detect_seams: bool = True, + seam_params: SeamDetectionParams | Mapping[str, Any] | None = None, n_jobs: int = -1, table_key_added: str | None = None, inplace: bool = True, @@ -500,6 +537,19 @@ def calculate_tiling_qc( a ``Mapping`` of its field names to values. See :class:`TilingQCParams` for each field's meaning and default. ``None`` (default) uses all defaults. + detect_seams + If ``True`` (default), also run emergent-seam detection: locate the + FOV seam grid from the collinear alignment of many cells' straight + cut-edges (no FOV size or tile overlap required) and flag cells whose + edge lies on a seam as ``is_seam_cut``. Unlike the MAD-based + ``is_outlier`` gate, this stays reliable in dense tissue, tolerates a + wide inter-FOV gap, and works when a cut leaves only one segmentable + half. The detected seams are recorded in ``.uns["tiling_qc"]["seams"]``. + seam_params + Advanced tuning knobs for seam detection as a + :class:`SeamDetectionParams` instance or a ``Mapping`` of its field + names. ``None`` (default) uses all defaults. Ignored when + ``detect_seams=False``. n_jobs Number of threads for tile processing. ``-1`` (default) uses all available CPUs; ``0`` and values below ``-1`` raise. @@ -533,6 +583,14 @@ def calculate_tiling_qc( neighbors that are smoothed-score outliers (MAD-based). Bounded [0, 1]; high values trace the tile grid. + When ``detect_seams=True`` two more columns are added: + + - ``is_seam_cut``: boolean, ``True`` for cells whose straight cardinal + edge lies on (and faces) a detected FOV seam. Recommended over + ``is_outlier`` for dense tissue / wide gaps / single-sided cuts. + - ``seam_dist``: distance (px) from the cell's edge to its seam-band + centre; ``NaN`` for cells that are not seam cuts. + Notes ----- Tile processing is parallelised via dask. When an active @@ -558,6 +616,7 @@ def calculate_tiling_qc( if n_neighbors < 1: raise ValueError(f"n_neighbors must be >= 1, got {n_neighbors}.") qc_params = _resolve_qc_params(tiling_qc_params) + resolved_seam = resolve_params(seam_params, SeamDetectionParams, label="seam_params") if detect_seams else None labels_da = resolve_labels_array(sdata, labels_key, scale) @@ -573,20 +632,36 @@ def calculate_tiling_qc( f"Tiling QC: {len(specs)} tiles ({tile_size}x{tile_size}, margin={overlap_margin}, downsample={downsample}x)." ) + # Data length scale D (median cell equivalent diameter) -> resolves every seam threshold to + # pixels once, so detection and the stitcher's edge extraction cannot disagree. + seam_diameter: float | None = None + seam_scale: SeamScale | None = None + if resolved_seam is not None: + _sizes = np.array([np.sqrt(ci.bbox_h * ci.bbox_w) for ci in cell_info.values()], dtype=float) + seam_diameter = float(np.median(_sizes)) if _sizes.size else 1.0 + seam_scale = resolved_seam._resolve(seam_diameter) + def _process_one(spec): - tile_lbl = extract_labels_tile_lazy(labels_da, spec) + if seam_scale is None: + tile_lbl, occupancy = extract_labels_tile_lazy(labels_da, spec), None + else: + tile_lbl, occupancy = extract_labels_tile_with_occupancy(labels_da, spec) return _score_tile( tile_lbl, distance_tol=qc_params.distance_tol, min_area=qc_params.min_area, downsample=downsample, max_contour_points=qc_params.max_contour_points, + origin=(spec.crop[0], spec.crop[1]), + seam_scale=seam_scale, + probe_occupancy=occupancy, ) # `_score_tile` is numba `nogil`, so threads scale (no process/pickle cost). results = _run_tiled(specs, _process_one, n_jobs=n_jobs, kind="threads", desc="tiles") - tile_dfs = [df for df in results if not df.empty] + tile_dfs = [df for df, _edges in results if not df.empty] + all_edges = [e for _df, edges in results for e in edges] if not tile_dfs: raise ValueError("No cells scored - labels may be empty or all below min_area.") @@ -649,6 +724,36 @@ def _process_one(spec): neighbor_outlier_frac = combined["is_outlier"].values[neighbor_idx].mean(axis=1) combined["nhood_outlier_fraction"] = neighbor_outlier_frac + # --- Emergent-seam cut-cell detection (FOV-size-agnostic, single-sided-friendly) --- + # Two stages: (1) locate the seam grid from wide-gap edges (a consensus over many cells; + # a lone straight membrane between two touching cells is off-grid and ignored); (2) flag + # any cell with a cardinal edge on a detected seam -- the gap is no longer required, so + # partial-edge and two-sided close-gap cuts are caught. Every threshold derives from the + # data: lengths from `seam_diameter`, the gap split and its channel choice from the observed + # gap distribution, and peak height from a Poisson null (note: not the `scale` argument). + seams_uns: dict[str, list[dict[str, float]]] = {"v": [], "h": []} + seam_gap_threshold: float | None = None + seam_gap_selectivity: float | None = None + if seam_scale is not None: + # One pass over the edge gaps: the channel decision is also the `.uns` diagnostic. + channel = gap_channel(all_edges, resolved_seam) + gap_thresh, selectivity, _use_wide = channel + seam_gap_threshold = float(gap_thresh) if np.isfinite(gap_thresh) else None + seam_gap_selectivity = float(selectivity) + seams = _detect_seam_bands(all_edges, W, H, seam_scale, channel) # stage 1 + edges_by_cell: dict[int, list[dict[str, Any]]] = {} + for e in all_edges: + edges_by_cell.setdefault(e["cell_id"], []).append(e) + flagged = flag_cells_on_seams(edges_by_cell, seams, seam_scale) # stage 2 + # `is_seam_cut` is by definition "has a seam distance", so derive it rather than + # building two parallel columns that have to stay in sync. + combined["seam_dist"] = combined.index.map(flagged).astype(float) + combined["is_seam_cut"] = combined["seam_dist"].notna() + seams_uns = { + axis: [{"coord": c, "half_width": h, "n_edges": int(k)} for c, h, k in bands] + for axis, bands in seams.items() + } + adata = ad.AnnData( X=np.empty((n_cells, 0), dtype=np.float32), ) @@ -683,6 +788,14 @@ def _process_one(spec): "nmads_smoothed": nmads_smoothed, "n_neighbors": n_neighbors, "tiling_qc_params": asdict(qc_params), + "detect_seams": resolved_seam is not None, + "seam_params": asdict(resolved_seam) if resolved_seam is not None else None, + "seams": seams_uns, + # Data length scale needed by assign_stitch_groups' seam-aware edge extraction. + "seam_diameter": seam_diameter, + # Diagnostics: where the gap split landed and whether it was selective enough to use. + "seam_gap_threshold": seam_gap_threshold, + "seam_gap_selectivity": seam_gap_selectivity, } if inplace: diff --git a/src/squidpy/experimental/tl/_tiling_stitch.py b/src/squidpy/experimental/tl/_tiling_stitch.py index 7df08c629..f587e0197 100644 --- a/src/squidpy/experimental/tl/_tiling_stitch.py +++ b/src/squidpy/experimental/tl/_tiling_stitch.py @@ -2,17 +2,28 @@ When segmentation is run tile-by-tile (Cellpose, Stardist, Mesmer, ...) cells that straddle tile boundaries get cut into 2-4 pieces with characteristic -straight, axis-aligned cut edges. :func:`~squidpy.experimental.tl.calculate_tiling_qc` flags these -as ``is_outlier=True``. This module pairs facing cut edges across boundaries -and assigns each candidate pair a heuristic geometric score in [0, 1]. - -The score is the flat (unweighted) mean of five dataset-independent geometric -features -- ``iou``, ``endpoint_match``, ``merge_compactness``, -``merge_solidity`` and ``gap_proximity`` -- computed from the cut-edge geometry -and the union mask after closing the seam gap. No model is fitted or shipped; -the features are recorded in ``.uns["tiling_stitch"]``. Users should tune -``min_confidence`` for their data; ``0.7`` is a reasonable starting point, not -a calibrated probability. +straight, axis-aligned cut edges. :func:`~squidpy.experimental.tl.calculate_tiling_qc` +flags these (``is_seam_cut`` when seam detection is enabled, else ``is_outlier``). +This module pairs facing cut edges across boundaries and assigns each candidate +pair a heuristic geometric score in [0, 1]. + +The score is the flat (unweighted) mean of four dataset-independent geometric +features -- ``iou``, ``endpoint_match``, ``merge_compactness`` and +``merge_solidity`` -- computed from the cut-edge geometry and the union mask +after closing the seam gap. No model is fitted or shipped; the features are +recorded in ``.uns["tiling_stitch"]``. Users should tune ``min_confidence`` +for their data; ``0.6`` is a reasonable starting point, not a calibrated +probability. + +Everything is scale-invariant: candidate enumeration is rank-based (k nearest +facing edges, no absolute-pixel search radius), cut-edge lengths are relative to +the data's median cell diameter ``D`` (read from +``.uns["tiling_qc"]["seam_diameter"]``), and how far apart two halves of one cut +may lie is bounded by the width of the seam band they sit on -- measured by +``calculate_tiling_qc``, not assumed. Cut edges are extracted only on +and facing the seam bands detected by ``calculate_tiling_qc`` -- so pairing no +longer merges touching interior cells, and both halves of a genuine cut are +recovered even when a 1-px segmentation fringe pushes the flat cut inward. The labels element is **never** modified here -- only ``.obs`` columns are written. Materialising a stitched labels element is opt-in via @@ -23,19 +34,19 @@ from collections.abc import Mapping from dataclasses import asdict, dataclass, field -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, Literal import numpy as np import spatialdata as sd import xarray as xr -from scipy.ndimage import binary_closing +from scipy.ndimage import distance_transform_edt from scipy.sparse import csr_matrix from scipy.sparse.csgraph import connected_components -from skimage.measure import find_contours, regionprops from skimage.measure import label as cc_label -from skimage.morphology import disk as morph_disk +from skimage.measure import regionprops from spatialdata._logging import logger as logg +from squidpy.experimental.tl._seam import SeamDetectionParams, SeamScale, cell_flat_edges, seam_offset from squidpy.experimental.utils._labels import iter_chunked_regionprops, resolve_labels_array from squidpy.experimental.utils._params import resolve_params @@ -46,8 +57,10 @@ __all__ = ["StitchParams", "assign_stitch_groups"] -# The geometric features whose flat mean is the stitch score. -_SCORE_FEATURES: tuple[str, ...] = ("iou", "endpoint_match", "merge_compactness", "merge_solidity", "gap_proximity") +# The geometric features whose flat mean is the stitch score. `gap_proximity` was intentionally +# dropped: with a per-pair `close_radius` that bridges each pair's own seam gap, it penalised +# wide-but-genuine seams and lowered recall without helping precision (validated on ground truth). +_SCORE_FEATURES: tuple[str, ...] = ("iou", "endpoint_match", "merge_compactness", "merge_solidity") # The subset computed by the expensive merge-union step; the rest are cheap # geometry features known before it, which drives the scoring early-prune. _SHAPE_FEATURES: tuple[str, ...] = ("merge_compactness", "merge_solidity") @@ -63,45 +76,30 @@ class StitchParams: advanced knobs -- the defaults rarely need changing. """ - distance_tol: float = 0.75 - """Sub-pixel tolerance for "lies on a bbox edge".""" - - min_edge_length: float = 5.0 - """Absolute floor on cut-edge length (pixels).""" - - min_edge_length_ratio: float = 0.4 - """Minimum cut-edge length relative to the cell's equivalent diameter.""" - - min_edge_coverage: float = 0.5 - """Minimum fraction of parallel-axis positions covered by near-edge contour points.""" - candidate_min_iou: float = 0.2 - """Loose 1-D IoU floor at candidate enumeration.""" + """Loose 1-D along-seam IoU floor for a facing edge to be a pair candidate.""" + + k_neighbors: int = 5 + """Rank-based candidate cap -- each cut edge is paired only with its ``k`` nearest + *facing* edges (by perpendicular gap). Replaces an absolute ``max_gap`` pixel + threshold, so the search adapts to the dataset's own seam-gap width.""" - close_radius: int = 3 - """Morphological closing disk radius for the union mask. Also the - length scale for ``gap_proximity`` (normalised by ``2 * close_radius``).""" + close_radius_min: int = 2 + """Floor for the per-pair morphological closing radius. The effective radius is + ``max(close_radius_min, ceil(gap / 2) + 1)`` so closing always bridges that pair's + own seam gap before the union's solidity/compactness are measured.""" def __post_init__(self) -> None: # Coerce numeric types (accept numpy scalars cleanly) and bounds-check. - self.distance_tol = float(self.distance_tol) - self.min_edge_length = float(self.min_edge_length) - self.min_edge_length_ratio = float(self.min_edge_length_ratio) - self.min_edge_coverage = float(self.min_edge_coverage) self.candidate_min_iou = float(self.candidate_min_iou) - self.close_radius = int(self.close_radius) - if self.distance_tol < 0: - raise ValueError(f"distance_tol must be >= 0, got {self.distance_tol}.") - if self.min_edge_length < 0: - raise ValueError(f"min_edge_length must be >= 0, got {self.min_edge_length}.") - if not 0.0 <= self.min_edge_length_ratio <= 1.0: - raise ValueError(f"min_edge_length_ratio must be in [0, 1], got {self.min_edge_length_ratio}.") - if not 0.0 <= self.min_edge_coverage <= 1.0: - raise ValueError(f"min_edge_coverage must be in [0, 1], got {self.min_edge_coverage}.") + self.k_neighbors = int(self.k_neighbors) + self.close_radius_min = int(self.close_radius_min) if not 0.0 <= self.candidate_min_iou <= 1.0: raise ValueError(f"candidate_min_iou must be in [0, 1], got {self.candidate_min_iou}.") - if self.close_radius < 0: - raise ValueError(f"close_radius must be >= 0, got {self.close_radius}.") + if self.k_neighbors < 1: + raise ValueError(f"k_neighbors must be >= 1, got {self.k_neighbors}.") + if self.close_radius_min < 0: + raise ValueError(f"close_radius_min must be >= 0, got {self.close_radius_min}.") def _resolve_stitch_params(stitch_params: StitchParams | Mapping[str, Any] | None) -> StitchParams: @@ -117,7 +115,7 @@ def _resolve_stitch_params(stitch_params: StitchParams | Mapping[str, Any] | Non # is the subset of top-level kwargs valid for re-running assign_stitch_groups # (the advanced tuning lives in a nested ``stitch_params`` dict). _STITCH_COLUMNS = ("stitch_group_id", "is_stitched", "n_pieces", "stitch_confidence") -_STITCH_PARAM_KEYS = frozenset({"min_confidence", "max_gap", "max_group_size"}) +_STITCH_PARAM_KEYS = frozenset({"min_confidence", "max_group_size"}) # Dataclasses @@ -168,7 +166,6 @@ class _StitchPair: confidence: float iou: float endpoint_match: float - gap_proximity: float merge_solidity: float merge_compactness: float edge_a: _CutEdge | None = field(default=None, repr=False) @@ -220,130 +217,81 @@ def _compute_outlier_bboxes( return bboxes -def _bbox_edge_run( - contour: np.ndarray, - perp_axis: int, - target: float, - distance_tol: float = _STITCH_DEFAULTS.distance_tol, - min_coverage: float = _STITCH_DEFAULTS.min_edge_coverage, -) -> tuple[float, float, float] | None: - """Find the extent of contour points lying near a single bbox edge. - - A genuine cut edge has many contour points clustered at the bbox boundary, - spanning a long parallel-axis range with high integer-position coverage. - A naturally curved cell only touches its bbox at a single point, which - fails either the count, length, or coverage check. - - Returns ``(ext_lo, ext_hi, length)`` if a substantial run is found. - """ - parallel_axis = 1 - perp_axis - near = np.abs(contour[:, perp_axis] - target) <= distance_tol - if near.sum() < 3: - return None - parallel_vals = contour[near, parallel_axis] - ext_lo = float(parallel_vals.min()) - ext_hi = float(parallel_vals.max()) - length = ext_hi - ext_lo - if length <= 0: - return None - width = max(int(np.ceil(length)), 1) - bins = np.zeros(width + 1, dtype=bool) - bins[np.clip((parallel_vals - ext_lo).astype(int), 0, width)] = True - coverage = float(bins.sum()) / (width + 1) - if coverage < min_coverage: - return None - return ext_lo, ext_hi, length - - def _extract_cut_edges( labels_da: xr.DataArray | np.ndarray, outlier_ids: Iterable[int], - bboxes: dict[int, tuple[int, int, int, int]] | None = None, - distance_tol: float = _STITCH_DEFAULTS.distance_tol, - min_edge_length: float = _STITCH_DEFAULTS.min_edge_length, - min_edge_length_ratio: float = _STITCH_DEFAULTS.min_edge_length_ratio, - min_edge_coverage: float = _STITCH_DEFAULTS.min_edge_coverage, + bboxes: dict[int, tuple[int, int, int, int]], + seams: dict[str, list[tuple[float, float, int]]], + scale: SeamScale, ) -> tuple[list[_CutEdge], dict[int, np.ndarray]]: - """Extract cardinal-aligned bbox-edge runs (cut-edge candidates) per outlier. - - For each outlier cell: - 1. Crop labels to its bbox + 1 px pad, build a binary mask. - 2. Trace its contour with :func:`skimage.measure.find_contours`. - 3. Check each of the 4 bbox-edge lines for a substantial straight run. + """Extract cut edges on and facing the detected seam bands, per outlier cell. - A piece cut at a tile boundary always has its cut on a bbox edge -- the - piece terminates exactly at the cut. Curved cells only touch the bbox - at a single contour point, which the density check rejects. + For each cut cell we take its dominant flat boundary line on each side + (:func:`~squidpy.experimental.tl._seam.cell_flat_edges` -- fringe-robust: the + flat line is found at *any* coordinate, not pinned to the bbox extreme, so a + 1-px segmentation fringe no longer drops the edge) and keep only edges that + lie on a detected seam band and face it. This mirrors the ``is_seam_cut`` + flagging in :func:`~squidpy.experimental.tl.calculate_tiling_qc`, so the two + stages agree on which edges are seam cuts. - Cells at a 4-tile corner produce 2 perpendicular edges; mid-stripe pieces - can produce 2 parallel edges. + ``scale`` is the detection scale rehydrated from ``.uns["tiling_qc"]``, so the + thresholds here are literally the ones the cells were flagged with. Returns ------- - The list of cut edges and, as a by-product of the per-cell crop already - read here, a ``{label_id -> boolean bbox mask}`` dict that lets the scoring - pass reconstruct merge unions in memory without re-reading the labels array. + The list of cut edges and, as a by-product of the per-cell crop already read + here, a ``{label_id -> boolean bbox mask}`` dict that lets the scoring pass + reconstruct merge unions in memory without re-reading the labels array. """ - outlier_list = [int(x) for x in outlier_ids] - if bboxes is None: - bboxes = _compute_outlier_bboxes(labels_da, outlier_list) + probe_depth = scale.probe_depth edges: list[_CutEdge] = [] outlier_crops: dict[int, np.ndarray] = {} - for lid in outlier_list: + for lid in [int(x) for x in outlier_ids]: bbox = bboxes.get(lid) if bbox is None: continue min_r, min_c, max_r, max_c = bbox - - crop_arr = _read_bbox_slice(labels_da, min_r, max_r, min_c, max_c) - cell_mask = crop_arr == lid # boolean bbox mask; reused by the scoring pass + # Read the bbox padded by the probe depth: cell_flat_edges probes background + # *beyond* the cut edge to measure the gap, so it needs the neighbourhood. + pad = probe_depth + 2 + r0 = max(0, min_r - pad) + c0 = max(0, min_c - pad) + crop = _read_bbox_slice(labels_da, r0, max_r + pad, c0, max_c + pad) + h = max_r - min_r + w = max_c - min_c + by0 = min_r - r0 + bx0 = min_c - c0 + # Slice to the cell's own bbox before comparing: the crop is padded by the probe depth, + # so comparing first would test several times the pixels that are kept. + cell_mask = crop[by0 : by0 + h, bx0 : bx0 + w] == lid # boolean bbox mask; reused by scoring if not cell_mask.any(): continue outlier_crops[lid] = cell_mask - # 1px zero-pad so cells filling their bbox still trace a closed contour. - mask = np.pad(cell_mask.astype(np.float32), 1, mode="constant", constant_values=0) - contours = find_contours(mask, 0.5) - if not contours: # degenerate mask traces nothing; skip it - continue - contour = max(contours, key=len) - contour_global = contour.copy() - contour_global[:, 0] += min_r - 1 - contour_global[:, 1] += min_c - 1 - - # Local centroid from the mask (avoids a second regionprops call). - ys, xs = np.where(mask) - cy = float(ys.mean()) + min_r - 1 - cx = float(xs.mean()) + min_c - 1 - area = float(mask.sum()) - eq_diameter = float(np.sqrt(4 * area / np.pi)) # diameter of the equal-area circle - min_len = max(min_edge_length, min_edge_length_ratio * eq_diameter) - - # find_contours places level set 0.5 outside the integer pixel boundary. - bbox_targets = [ - ("h", float(min_r) - 0.5), - ("h", float(max_r) - 0.5), - ("v", float(min_c) - 0.5), - ("v", float(max_c) - 0.5), - ] - for axis, target in bbox_targets: - perp_axis = 0 if axis == "h" else 1 - run = _bbox_edge_run(contour_global, perp_axis, target, distance_tol, min_edge_coverage) - if run is None: - continue - ext_lo, ext_hi, length = run - if length < min_len: + + flat = cell_flat_edges( + cell_mask, + crop != 0, # occupancy: this crop is read unmasked, so neighbours are visible + (by0, bx0, by0 + h, bx0 + w), + (r0, c0), + scale.min_len, + scale.flat_tol, + probe_depth, + ) + for e in flat: + # Keep only edges lying on a detected seam band and facing it -- the *same* + # predicate `calculate_tiling_qc` flags cells with, so the two stages agree by + # construction rather than by comment. + if seam_offset(e, seams.get(e["axis"], []), scale) is None: continue - cell_coord = cy if axis == "h" else cx - normal = 1 if cell_coord > target else -1 edges.append( _CutEdge( cell_id=lid, - axis=axis, - coord=target, - extent=(ext_lo, ext_hi), - normal_dir=normal, - length=float(length), + axis=e["axis"], + coord=e["coord"], + extent=e["extent"], + normal_dir=e["side"], + length=float(e["span"]), ) ) @@ -362,7 +310,7 @@ def _merge_shape_features( cell_b: int, bboxes: dict[int, tuple[int, int, int, int]], outlier_crops: dict[int, np.ndarray], - close_radius: int = _STITCH_DEFAULTS.close_radius, + close_radius: int, *, H: int, W: int, @@ -402,7 +350,11 @@ def _merge_shape_features( if not mask.any(): return zero - closed = binary_closing(mask, structure=morph_disk(close_radius)) + # Closing by a disk of radius r, via distance transforms: pixel-identical to + # `binary_closing(mask, disk(r))` but O(N) instead of O(N * r^2). The radius tracks each + # pair's own seam gap, so the explicit-footprint form got dramatically more expensive + # exactly on the wide-gap data this feature targets. + closed = distance_transform_edt(distance_transform_edt(~mask) <= close_radius) > close_radius cc = cc_label(closed, connectivity=2) if cc.max() == 0: return zero @@ -421,11 +373,15 @@ def _merge_shape_features( def _pair_geometry_features( e: _CutEdge, c: _CutEdge, - max_gap: float, candidate_min_iou: float = _STITCH_DEFAULTS.candidate_min_iou, ) -> dict[str, float] | None: """Compute geometry-only features for a candidate pair, returning ``None`` if the pair fails the basic facing/overlap/IoU filters. + + No absolute-pixel gap cutoff is applied here -- the raw perpendicular ``gap`` + is returned so the (rank-based) enumerator can pick each edge's nearest + facing partners, and the diameter-relative plausibility guard is applied in + :func:`_score_pairs`. """ if c.normal_dir == e.normal_dir: return None @@ -440,13 +396,9 @@ def _pair_geometry_features( if iou < candidate_min_iou: return None gap = abs(e.coord - c.coord) - if gap > max_gap: - return None endpoint_dist = abs(e.extent[0] - c.extent[0]) + abs(e.extent[1] - c.extent[1]) max_len = max(e.length, c.length) endpoint_match = max(0.0, 1.0 - endpoint_dist / max_len) if max_len > 0 else 0.0 - # Return the raw perpendicular gap; gap_proximity is derived later against - # the closing reach (2*close_radius), NOT against max_gap (a search radius). return { "iou": float(iou), "endpoint_match": float(endpoint_match), @@ -456,61 +408,47 @@ def _pair_geometry_features( def _enumerate_pair_candidates( edges: list[_CutEdge], - max_gap: float, + k_neighbors: int = _STITCH_DEFAULTS.k_neighbors, candidate_min_iou: float = _STITCH_DEFAULTS.candidate_min_iou, ) -> list[tuple[_CutEdge, _CutEdge, dict[str, float]]]: - """Find all (e, c) pairs of facing cut edges with their geometry features. + """Find candidate pairs of facing cut edges, rank-based (no absolute max_gap). - Returns one entry per surviving candidate. No selection / scoring yet. + For each edge, keep only its ``k_neighbors`` nearest *facing + overlapping* + edges by perpendicular gap. A rank-based cap adapts to the dataset's own + seam-gap width -- unlike an absolute pixel radius, which fails when the + inter-FOV gap is wider than a hand-tuned constant. Each unordered + ``(cell_a, cell_b, axis)`` pair is emitted once. No scoring yet. """ - out: list[tuple[_CutEdge, _CutEdge, dict[str, float]]] = [] + out: dict[tuple[int, int, str], tuple[_CutEdge, _CutEdge, dict[str, float]]] = {} by_axis: dict[str, list[_CutEdge]] = {"h": [], "v": []} for e in edges: by_axis[e.axis].append(e) for axis_edges in by_axis.values(): - axis_edges.sort(key=lambda e: e.coord) - coords = np.array([e.coord for e in axis_edges]) - for i, e in enumerate(axis_edges): - lo = int(np.searchsorted(coords, e.coord - max_gap, side="left")) - hi = int(np.searchsorted(coords, e.coord + max_gap, side="right")) - for j in range(lo, hi): - if j <= i: - continue # symmetry: emit each unordered pair once - c = axis_edges[j] - if c.cell_id == e.cell_id: + for e in axis_edges: + facing: list[tuple[float, _CutEdge, dict[str, float]]] = [] + for c in axis_edges: + if c is e or c.cell_id == e.cell_id: continue - feats = _pair_geometry_features(e, c, max_gap, candidate_min_iou=candidate_min_iou) + feats = _pair_geometry_features(e, c, candidate_min_iou=candidate_min_iou) if feats is None: continue - out.append((e, c, feats)) - return out + facing.append((feats["gap"], c, feats)) + facing.sort(key=lambda t: t[0]) + for _g, c, feats in facing[:k_neighbors]: + key = (min(e.cell_id, c.cell_id), max(e.cell_id, c.cell_id), e.axis) + if key not in out: + out[key] = (e, c, feats) + return list(out.values()) # Scoring -def _gap_proximity(gap: float, close_radius: int) -> float: - """Map the raw perpendicular gap to [0, 1] against the closing reach. - - Normalised by ``2 * close_radius`` -- the scale at which morphological - closing could actually bridge the seam -- so the feature is independent of - the ``max_gap`` search radius and only reaches 0 when the gap genuinely - exceeds what closing can join. When closing is disabled (``close_radius=0``) - the feature is inactive and returns ``1.0`` rather than collapsing the score. - """ - reach = 2 * close_radius - # gap<=0 (touching/overlapping) or reach<=0 (closing disabled, close_radius=0) - # -> the feature is inactive (neutral 1.0), never a silent score cliff. - if gap <= 0 or reach <= 0: - return 1.0 - return max(0.0, 1.0 - gap / reach) - - def _score_pair_features(features: dict[str, float]) -> float: """Return the heuristic stitch score in [0, 1]. - Flat (unweighted) mean of the five features in :data:`_SCORE_FEATURES`. + Flat (unweighted) mean of the four features in :data:`_SCORE_FEATURES`. The score is dataset-independent and not a calibrated probability -- users pick ``min_confidence`` based on their false-merge tolerance. """ @@ -527,29 +465,55 @@ def _max_achievable_score(known_features: dict[str, float]) -> float: return _score_pair_features({**known_features, **dict.fromkeys(_SHAPE_FEATURES, 1.0)}) +def _seam_span(seams: dict[str, list[tuple[float, float, int]]], axis: str, coord_a: float, coord_b: float) -> float: + """Widest plausible separation between two halves of one cut, from the measured seam band. + + Both halves sit on the inner edges of the same band, so their perpendicular separation + cannot exceed that band's width. Tying the bound to the detected seam rather than to a + multiple of the cell diameter also bounds the per-pair closing radius, which otherwise + grows with the gap until it bridges pieces that were never one cell. + """ + bands = seams.get(axis, []) + if not bands: + return 0.0 + mid = (coord_a + coord_b) / 2.0 + _centre, half, _cnt = min(bands, key=lambda b: abs(mid - b[0])) + return 2.0 * half + + def _score_pairs( candidates: list[tuple[_CutEdge, _CutEdge, dict[str, float]]], bboxes: dict[int, tuple[int, int, int, int]], outlier_crops: dict[int, np.ndarray], min_confidence: float, - close_radius: int = _STITCH_DEFAULTS.close_radius, + seams: dict[str, list[tuple[float, float, int]]], *, + close_radius_min: int = _STITCH_DEFAULTS.close_radius_min, H: int, W: int, ) -> list[_StitchPair]: """Compute shape features per candidate, score, and keep pairs >= min_confidence. - One entry per ``(cell_a, cell_b, axis)`` (keeping max confidence on duplicates). + The morphological closing radius is chosen *per pair* to bridge that pair's + own seam gap (``max(close_radius_min, ceil(gap / 2) + 1)``), and candidates + separated by more than the width of the seam band they lie on are discarded + (:func:`_seam_span`), which also bounds that radius. One entry per + ``(cell_a, cell_b, axis)`` (keeping max confidence on duplicates). """ scored: list[_StitchPair] = [] for e, c, geom in candidates: - known = {**geom, "gap_proximity": _gap_proximity(geom["gap"], close_radius)} + gap = geom["gap"] + if gap > _seam_span(seams, e.axis, e.coord, c.coord): # wider than its own seam: not one cut + continue + # Per-pair closing radius bridges this pair's own gap; the gap is already bounded by + # the seam band's width above, so the radius needs no separate cap. + close_radius = max(close_radius_min, int(np.ceil(gap / 2.0)) + 1) # Skip the costly union reconstruction when even the best case for the # deferred shape features can't reach min_confidence. - if _max_achievable_score(known) < min_confidence: + if _max_achievable_score(geom) < min_confidence: continue shape = _merge_shape_features(e.cell_id, c.cell_id, bboxes, outlier_crops, close_radius=close_radius, H=H, W=W) - feats = {**known, **shape} + feats = {**geom, **shape} confidence = _score_pair_features(feats) if confidence < min_confidence: continue @@ -566,7 +530,6 @@ def _score_pairs( confidence=confidence, iou=feats["iou"], endpoint_match=feats["endpoint_match"], - gap_proximity=feats["gap_proximity"], merge_solidity=feats["merge_solidity"], merge_compactness=feats["merge_compactness"], edge_a=ea, @@ -574,13 +537,7 @@ def _score_pairs( ) ) - # Deduplicate to one entry per (cell_a, cell_b, axis), keeping max confidence. - by_pair: dict[tuple[int, int, str], _StitchPair] = {} - for p in scored: - k = (p.cell_a, p.cell_b, p.axis) - if k not in by_pair or by_pair[k].confidence < p.confidence: - by_pair[k] = p - return sorted(by_pair.values(), key=lambda p: (-p.confidence, p.cell_a, p.cell_b)) + return sorted(scored, key=lambda p: (-p.confidence, p.cell_a, p.cell_b)) # Group assembly (union-find + validation) @@ -589,7 +546,7 @@ def _score_pairs( def _validate_group_geometry( pairs_in_group: list[_StitchPair], size: int, - max_gap: float, + gap_tol: float, ) -> bool: """Geometric sanity check for groups of size >= 3. @@ -598,7 +555,7 @@ def _validate_group_geometry( - **Corner group** (size 4, both axes present): the cut edges' endpoints must converge near a single junction point (one ``h`` cut crossing one ``v`` cut defines the junction). If the spread of edge extents from - the junction is greater than ``max_gap``, the group is implausible. + the junction is greater than ``gap_tol``, the group is implausible. - **Chain group** (size 3 or 4, all pairs share one axis): legitimate same-axis chains (e.g., a cell split by 3 horizontal seams into 4 @@ -606,6 +563,9 @@ def _validate_group_geometry( coordinates. Multiple pairs at the same seam coord would imply geometrically impossible "two cuts at the same seam" pairings -- a signature of a false-positive cluster -- so we reject. + + ``gap_tol`` is the width of the widest detected seam band, used consistently + with the per-pair bound applied during candidate scoring. """ h_pairs = [p for p in pairs_in_group if p.axis == "h"] v_pairs = [p for p in pairs_in_group if p.axis == "v"] @@ -616,10 +576,10 @@ def _validate_group_geometry( return True # 2-piece groups are trivially valid on one axis # Each pair's seam coord is roughly midway between its two edges. seam_coords = [round((p.edge_a.coord + p.edge_b.coord) / 2.0, 1) for p in pairs_in_group] - # Allow a max_gap-sized tolerance for "distinct" seams. + # Allow a gap_tol-sized tolerance for "distinct" seams. sorted_coords = sorted(seam_coords) for prev, cur in zip(sorted_coords, sorted_coords[1:], strict=False): - if cur - prev <= max_gap: + if cur - prev <= gap_tol: return False return True @@ -635,10 +595,10 @@ def _validate_group_geometry( junction_y = float(np.mean([e.coord for e in h_edges])) junction_x = float(np.mean([e.coord for e in v_edges])) for e in h_edges: - if min(abs(e.extent[0] - junction_x), abs(e.extent[1] - junction_x)) > max_gap: + if min(abs(e.extent[0] - junction_x), abs(e.extent[1] - junction_x)) > gap_tol: return False for e in v_edges: - if min(abs(e.extent[0] - junction_y), abs(e.extent[1] - junction_y)) > max_gap: + if min(abs(e.extent[0] - junction_y), abs(e.extent[1] - junction_y)) > gap_tol: return False return True @@ -647,7 +607,7 @@ def _assemble_groups( pairs: list[_StitchPair], candidate_ids: Iterable[int], max_group_size: int, - max_gap: float, + gap_tol: float, ) -> tuple[dict[int, int], dict[int, float]]: """Build stitch groups via union-find with size + corner validation. @@ -711,7 +671,7 @@ def _assemble_groups( # Geometric validation for 3+ piece groups: corner-junction for # mixed-axis 4-groups, chain (distinct seam coords) for same-axis 3+. - if size >= 3 and not _validate_group_geometry(group_pairs, size, max_gap): + if size >= 3 and not _validate_group_geometry(group_pairs, size, gap_tol): for m in mem: groups[m] = m confidences[m] = 1.0 @@ -738,29 +698,34 @@ def assign_stitch_groups( sdata: sd.SpatialData, labels_key: str, qc_table_key: str | None = None, - min_confidence: float = 0.7, - max_gap: float = 3.0, + min_confidence: float = 0.6, max_group_size: int = 4, + candidates: Literal["auto", "is_seam_cut", "is_outlier"] = "auto", stitch_params: StitchParams | Mapping[str, Any] | None = None, inplace: bool = True, ) -> ad.AnnData | None: """Assign tile-cut cell pieces to stitch groups. - Reads ``is_outlier=True`` cells flagged by - :func:`~squidpy.experimental.tl.calculate_tiling_qc`, pairs facing cut - edges across tile boundaries, scores each pair via a transparent geometric - composite, and assembles high-confidence pairs into stitch groups via - union-find. This only *annotates* which pieces belong together -- it does - **not** modify the labels element. Materialising a stitched labels element - is opt-in via :func:`!make_stitched_labels`. + Reads the cells flagged by :func:`~squidpy.experimental.tl.calculate_tiling_qc` + (``is_seam_cut`` by default), extracts each piece's cut edges *on and facing* + the detected FOV seam bands, pairs facing edges across boundaries (rank-based: + each edge with its ``k`` nearest facing partners -- no absolute-pixel search + radius), scores each pair via a transparent geometric composite, and assembles + high-confidence pairs into stitch groups via union-find. This only *annotates* + which pieces belong together -- it does **not** modify the labels element. + Materialising a stitched labels element is opt-in via :func:`!make_stitched_labels`. - The score per pair is the flat (unweighted) mean of five geometric features + The score per pair is the flat (unweighted) mean of four geometric features in [0, 1]: ``iou`` (1-D extent overlap), ``endpoint_match`` (chord endpoints - coincide), ``merge_compactness`` (``4*pi*A / P^2`` of the closed union mask), - ``merge_solidity`` (union area / convex hull area), and ``gap_proximity`` - (seam gap relative to the morphological closing reach). No coefficients are + coincide), ``merge_compactness`` (``4*pi*A / P^2`` of the closed union mask) + and ``merge_solidity`` (union area / convex hull area). No coefficients are fitted or shipped; the features are recorded in ``.uns["tiling_stitch"]``. + **Requires seam detection.** Run ``calculate_tiling_qc(..., detect_seams=True)`` + first: the seam bands and data length scale ``D`` it records in + ``.uns["tiling_qc"]`` localise cut-edge extraction to real seams, so pairing + no longer merges touching interior cells. + Parameters ---------- sdata @@ -771,16 +736,17 @@ def assign_stitch_groups( qc_table_key Key of the QC table. Defaults to ``"{labels_key}_qc"``. min_confidence - Threshold on ``stitch_confidence``. ``0.7`` (default) is a starting + Threshold on ``stitch_confidence``. ``0.6`` (default) is a starting point; raise it for stricter precision, lower for recall. Tune for your data -- the score is heuristic, not a calibrated probability. - max_gap - Maximum perpendicular distance (px) between facing cut edges for a pair - to be *considered* a candidate. This is a search radius only; it does - not scale the score. max_group_size Cap on group size; oversized groups (likely false merges) collapse to singletons. + candidates + Which QC column gates the cells considered for stitching. ``"auto"`` + (default) uses ``is_seam_cut`` when present -- these are localised to + detected FOV seams -- and otherwise falls back to ``is_outlier``. Set + explicitly to ``"is_seam_cut"`` or ``"is_outlier"`` to force one. stitch_params Advanced tuning knobs as a :class:`StitchParams` instance or a ``Mapping`` of its field names to values. See :class:`StitchParams` @@ -799,8 +765,6 @@ def assign_stitch_groups( raise ValueError(f"Labels key '{labels_key}' not found in sdata.labels.") if min_confidence < 0 or min_confidence > 1: raise ValueError(f"min_confidence must be in [0, 1], got {min_confidence}.") - if max_gap < 0: - raise ValueError(f"max_gap must be non-negative, got {max_gap}.") if max_group_size < 1: raise ValueError(f"max_group_size must be >= 1, got {max_group_size}.") params = _resolve_stitch_params(stitch_params) @@ -810,27 +774,53 @@ def assign_stitch_groups( raise ValueError(f"QC table '{table_key}' not found. Run calculate_tiling_qc first.") adata = sdata.tables[table_key].copy() - if "is_outlier" not in adata.obs.columns: - raise ValueError(f"QC table '{table_key}' is missing 'is_outlier'; re-run calculate_tiling_qc.") if "label_id" not in adata.obs.columns: raise ValueError(f"QC table '{table_key}' is missing 'label_id'.") + # Candidate gate: prefer the seam-aware `is_seam_cut` flag (localised to detected FOV seams, + # so pairing no longer merges touching interior cells); fall back to the MAD `is_outlier`. + # Seam data is mandatory below, so `is_seam_cut` is always present under "auto". + gate_col = "is_seam_cut" if candidates == "auto" else candidates + if gate_col not in adata.obs.columns: + raise ValueError( + f"QC table '{table_key}' is missing '{gate_col}'; re-run calculate_tiling_qc " + f"(with detect_seams=True for 'is_seam_cut')." + ) existing = [c for c in _STITCH_COLUMNS if c in adata.obs.columns] if existing: logg.warning(f"Overwriting existing stitch columns: {existing}.") adata.obs.drop(columns=existing, inplace=True) - # Resolve which labels DataArray was used at QC time (multi-scale aware). + # Seam contract from calculate_tiling_qc(detect_seams=True): the seam bands localise + # cut-edge extraction and `seam_diameter` sets every length scale. The extraction + # fractions mirror detection (read from seam_params) so both stages agree on seam cuts. qc_params = adata.uns.get("tiling_qc", {}) scale = qc_params.get("scale") + seams_uns = qc_params.get("seams") + diameter = qc_params.get("seam_diameter") + if not seams_uns or diameter is None: + raise ValueError( + f"QC table '{table_key}' has no seam detection results; assign_stitch_groups requires " + f"seam-localised cut edges. Re-run calculate_tiling_qc(..., detect_seams=True)." + ) + diameter = float(diameter) + seams = { + axis: [(float(b["coord"]), float(b["half_width"]), int(b["n_edges"])) for b in seams_uns.get(axis, [])] + for axis in ("v", "h") + } + # Rehydrate the exact parameters detection ran with, and resolve them against the same D. + seam_scale = resolve_params(qc_params.get("seam_params"), SeamDetectionParams, label="seam_params")._resolve( + diameter + ) + labels_da = resolve_labels_array(sdata, labels_key, scale) label_ids = adata.obs["label_id"].astype(int).to_numpy() - is_outlier = adata.obs["is_outlier"].to_numpy(dtype=bool) + is_outlier = adata.obs[gate_col].to_numpy(dtype=bool) outlier_ids = label_ids[is_outlier].tolist() n_outliers = len(outlier_ids) - logg.info(f"Stitching {n_outliers} outlier cells (out of {len(label_ids)} total).") + logg.info(f"Stitching {n_outliers} candidate cells ('{gate_col}', out of {len(label_ids)} total).") if n_outliers == 0: logg.warning("No outliers flagged; nothing to stitch.") @@ -846,21 +836,27 @@ def assign_stitch_groups( f"{len(missing)} outlier label_id(s) flagged in the QC table do not appear " f"in '{labels_key}' (e.g. {missing[:5]}); they will not be stitched." ) - edges, outlier_crops = _extract_cut_edges( - labels_da, - outlier_ids, - bboxes=bboxes, - distance_tol=params.distance_tol, - min_edge_length=params.min_edge_length, - min_edge_length_ratio=params.min_edge_length_ratio, - min_edge_coverage=params.min_edge_coverage, - ) + edges, outlier_crops = _extract_cut_edges(labels_da, outlier_ids, bboxes, seams, seam_scale) H, W = labels_da.shape[-2], labels_da.shape[-1] - candidates = _enumerate_pair_candidates(edges, max_gap=max_gap, candidate_min_iou=params.candidate_min_iou) + cand = _enumerate_pair_candidates( + edges, k_neighbors=params.k_neighbors, candidate_min_iou=params.candidate_min_iou + ) pairs = _score_pairs( - candidates, bboxes, outlier_crops, min_confidence, close_radius=params.close_radius, H=H, W=W + cand, + bboxes, + outlier_crops, + min_confidence, + seams, + close_radius_min=params.close_radius_min, + H=H, + W=W, + ) + # Tolerance for "distinct seams" / corner convergence: the widest detected band. + gap_tol = max( + (2.0 * half for bands in seams.values() for _c, half, _n in bands), + default=0.0, ) - groups, confidences = _assemble_groups(pairs, outlier_ids, max_group_size=max_group_size, max_gap=max_gap) + groups, confidences = _assemble_groups(pairs, outlier_ids, max_group_size=max_group_size, gap_tol=gap_tol) # Write .obs columns with three states distinguished by stitch_confidence: # - non-outlier cell -> own label_id, False, 1, NaN (not evaluated) @@ -902,8 +898,9 @@ def assign_stitch_groups( adata.uns[_METHOD_KEY] = { "min_confidence": float(min_confidence), - "max_gap": float(max_gap), "max_group_size": int(max_group_size), + "candidate_gate": gate_col, + "seam_diameter": float(diameter), "stitch_params": asdict(params), "n_outliers": int(n_outliers), "n_candidate_pairs": int(len(pairs)), diff --git a/tests/experimental/conftest.py b/tests/experimental/conftest.py index edc41a851..2027b387e 100644 --- a/tests/experimental/conftest.py +++ b/tests/experimental/conftest.py @@ -222,3 +222,89 @@ def make_clean_sdata() -> SpatialData: def sdata_clean() -> SpatialData: """Fixture wrapper around :func:`make_clean_sdata`.""" return make_clean_sdata() + + +# Dense-tissue + wide-gap + single-sided seam fixture (stresses emergent-seam detection) + +_DENSE_SIZE = 420 +_DENSE_BORDERS = (140, 280) # two internal seams per axis (3x3 tile grid) +_DENSE_GAP = 8 # wide inter-FOV gap (far wider than a touching-cell membrane) + + +@dataclass +class DenseSeamGroundTruth: + """Ground truth for the dense-tissue seam fixture.""" + + cut_cell_ids: frozenset[int] = field(default_factory=frozenset) + seam_coords: tuple[int, ...] = _DENSE_BORDERS + + +def make_dense_seam_sdata() -> tuple[SpatialData, DenseSeamGroundTruth]: + """Dense **touching** Voronoi cells cut by a wide seam, with single-sided drop-out. + + Mimics CosMx-style per-FOV segmentation on dense tissue: convex Voronoi cells with + thin membranes, a wide zeroed seam band, and small remnants dropped so many cuts leave + only one segmentable half. Convex cells never self-fragment, so the ground truth (a + cell is cut iff its pre-cut mask overlaps the seam band) is exact. + """ + from scipy.spatial import cKDTree + + rng = np.random.default_rng(0) + size, borders, gap = _DENSE_SIZE, _DENSE_BORDERS, _DENSE_GAP + pts = [(y + rng.integers(-4, 5), x + rng.integers(-4, 5)) for y in range(15, size, 15) for x in range(15, size, 15)] + tree = cKDTree(np.array(pts)) + yy, xx = np.mgrid[0:size, 0:size] + d, idx = tree.query(np.column_stack([yy.ravel(), xx.ravel()]), k=2) + cell = (idx[:, 0] + 1).reshape(size, size).astype(np.int32) + cell[((d[:, 1] - d[:, 0]) < 1.5).reshape(size, size)] = 0 # thin membranes + orig = cell.copy() + + half = gap // 2 + seam = np.zeros((size, size), bool) + for b in borders: + seam[b - half : b - half + gap, :] = True + seam[:, b - half : b - half + gap] = True + cut_orig = set(np.unique(orig[seam & (orig > 0)])) - {0} + orig_area = np.bincount(orig.ravel()) # one pass, indexed by label id + + seg = orig.copy() + seg[seam] = 0 + frag, nfrag = ndimage.label(seg > 0) + parent = np.zeros(nfrag + 1, np.int64) + frags_of: dict[int, list[int]] = {} + for f in range(1, nfrag + 1): + vals = orig[frag == f] + vals = vals[vals > 0] + parent[f] = np.bincount(vals).argmax() if vals.size else 0 + frags_of.setdefault(int(parent[f]), []).append(f) + + keep = np.ones(nfrag + 1, bool) + keep[0] = False + for o in cut_orig: + for f in frags_of.get(o, []): + if (frag == f).sum() / max(int(orig_area[o]), 1) < 0.45 and rng.random() < 0.75: + keep[f] = False + frag2, nf2 = ndimage.label(np.where(keep[frag], frag, 0) > 0) + + cut_ids: set[int] = set() + for f in range(1, nf2 + 1): + vals = orig[frag2 == f] + vals = vals[vals > 0] + if vals.size and int(np.bincount(vals).argmax()) in cut_orig: + cut_ids.add(f) + + labels_xr = xr.DataArray(da.from_array(frag2.astype(np.int32), chunks=(128, 128)), dims=["y", "x"]) + image_xr = xr.DataArray( + rng.integers(0, 255, (3, size, size), dtype=np.uint8), dims=["c", "y", "x"], coords={"c": ["R", "G", "B"]} + ) + sdata = SpatialData( + images={"image": Image2DModel.parse(image_xr)}, + labels={"labels": Labels2DModel.parse(labels_xr)}, + ) + return sdata, DenseSeamGroundTruth(cut_cell_ids=frozenset(cut_ids)) + + +@pytest.fixture() +def sdata_dense_seam() -> tuple[SpatialData, DenseSeamGroundTruth]: + """Fixture wrapper around :func:`make_dense_seam_sdata`.""" + return make_dense_seam_sdata() diff --git a/tests/experimental/test_seam.py b/tests/experimental/test_seam.py new file mode 100644 index 000000000..cd3269ff2 --- /dev/null +++ b/tests/experimental/test_seam.py @@ -0,0 +1,244 @@ +"""Tests for emergent-seam cut-cell detection in calculate_tiling_qc.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from squidpy.experimental.tl import SeamDetectionParams, calculate_tiling_qc +from squidpy.experimental.tl._seam import ( + _min_significant_count, + _otsu_gap_split, + cell_flat_edges, + detect_seams, + flag_cells_on_seams, + gap_channel, +) + + +def _assert_seams_match(seams, truth, context, tol=8): + """Assert each axis reports exactly the true seams, each within ``tol`` px.""" + for axis in ("v", "h"): + coords = sorted(b["coord"] for b in seams[axis]) + assert len(coords) == len(truth), ( + f"axis {axis} {context}: expected {list(truth)}, got {[round(c, 1) for c in coords]}" + ) + for want, got in zip(sorted(truth), coords, strict=True): + assert abs(got - want) <= tol, f"axis {axis} {context}: seam {got:.1f} too far from true {want}" + + +def _recall(pred_ids: set[int], truth_ids: frozenset[int]) -> float: + return len(pred_ids & truth_ids) / len(truth_ids) if truth_ids else 0.0 + + +def _precision(pred_ids: set[int], truth_ids: frozenset[int]) -> float: + return len(pred_ids & truth_ids) / len(pred_ids) if pred_ids else 0.0 + + +class TestSeamDetectionIntegration: + def test_seam_columns_added(self, sdata_dense_seam): + sdata, _ = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", detect_seams=True, inplace=False) + assert "is_seam_cut" in adata.obs + assert adata.obs["is_seam_cut"].dtype == bool + flagged = adata.obs["is_seam_cut"].to_numpy() + assert np.isfinite(adata.obs["seam_dist"].to_numpy()[flagged]).all() + assert np.isnan(adata.obs["seam_dist"].to_numpy()[~flagged]).all() + + def test_detect_seams_false_omits_columns(self, sdata_dense_seam): + sdata, _ = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", detect_seams=False, inplace=False) + assert "is_seam_cut" not in adata.obs + assert "seam_dist" not in adata.obs + assert adata.uns["tiling_qc"]["seams"] == {"v": [], "h": []} + + def test_seams_detected_near_true_border(self, sdata_dense_seam): + sdata, gt = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", detect_seams=True, inplace=False) + seams = adata.uns["tiling_qc"]["seams"] + found = [b["coord"] for b in seams["v"]] + [b["coord"] for b in seams["h"]] + assert found, "no seams detected" + for s in gt.seam_coords: + assert min(abs(c - s) for c in found) <= 8 + + def test_seam_recall_beats_outlier_on_dense_wide_gap(self, sdata_dense_seam): + sdata, gt = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", detect_seams=True, inplace=False) + lid = adata.obs["label_id"].to_numpy() + seam_ids = set(lid[adata.obs["is_seam_cut"].to_numpy()].tolist()) + outlier_ids = set(lid[adata.obs["is_outlier"].to_numpy()].tolist()) + # emergent-seam detection recovers most cuts where the MAD gate collapses + assert _recall(seam_ids, gt.cut_cell_ids) >= 0.7 + assert _recall(seam_ids, gt.cut_cell_ids) > _recall(outlier_ids, gt.cut_cell_ids) + 0.3 + assert _precision(seam_ids, gt.cut_cell_ids) >= 0.5 + + def test_seam_params_recorded_in_uns(self, sdata_dense_seam): + sdata, _ = sdata_dense_seam + adata = calculate_tiling_qc( + sdata, labels_key="labels", detect_seams=True, seam_params={"alpha": 0.005}, inplace=False + ) + assert adata.uns["tiling_qc"]["detect_seams"] is True + assert adata.uns["tiling_qc"]["seam_params"]["alpha"] == 0.005 + + def test_no_helper_columns_leak(self, sdata_dense_seam): + sdata, _ = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", detect_seams=True, inplace=False) + assert not any(c.startswith("_seam") for c in adata.obs.columns) + + +class TestSeamDetectionParams: + def test_defaults_are_dimensionless_or_pixel_constants(self): + p = SeamDetectionParams() + # length thresholds are fractions of the cell diameter (scale-invariant) + assert 0 < p.edge_len_frac <= 1 + # flat_tol / bin_width are pixel-grid constants (do not scale with cell size) + assert p.flat_tol >= 1 and p.bin_width >= 1 + + @pytest.mark.parametrize( + "kwargs,match", + [ + ({"edge_len_frac": 0.0}, "edge_len_frac"), + ({"gap_selectivity_max": 1.5}, "gap_selectivity_max"), + ({"alpha": 0.0}, "alpha"), + ({"probe_frac": 0.0}, "probe_frac"), + ], + ) + def test_invalid_raises(self, kwargs, match): + with pytest.raises(ValueError, match=match): + SeamDetectionParams(**kwargs) + + +class TestSeamUnits: + def test_flat_edge_detects_cut_with_gap(self): + # a cell filling the tile height, cut flat at col 10, wide background to the right + tile = np.zeros((20, 30), np.int32) + tile[0:20, 0:10] = 5 + from skimage.measure import regionprops + + rp = regionprops(tile)[0] + edges = cell_flat_edges(rp.image, tile != 0, rp.bbox, (0, 0), min_len=5, flat_tol=1.5, probe_depth=8) + vs = [e for e in edges if e["axis"] == "v" and e["side"] == -1] + assert vs and abs(vs[0]["coord"] - 9) <= 1 and vs[0]["gap"] >= 3 + + def test_flat_edge_gap_small_for_touching_neighbour(self): + # the flat side faces a neighbour ~1px away -> small gap (excluded from seam *detection*) + tile = np.zeros((20, 30), np.int32) + tile[0:20, 0:10] = 5 + tile[0:20, 11:21] = 6 + from skimage.measure import regionprops + + rp = next(r for r in regionprops(tile) if r.label == 5) + edges = cell_flat_edges(rp.image, tile != 0, rp.bbox, (0, 0), min_len=5, flat_tol=1.5, probe_depth=8) + vs = [e for e in edges if e["axis"] == "v" and e["side"] == -1] + assert vs and vs[0]["gap"] <= 2 # membrane, not a seam gap + + def test_otsu_splits_membranes_from_open_background(self): + # mostly 1px membranes plus a few wide gaps -> the split lands between the two modes + gaps = np.array([1, 1, 1, 1, 2, 8, 9, 10], dtype=float) + thr = _otsu_gap_split(gaps) + assert 2.0 < thr < 8.0 + + def test_gap_channel_filters_when_selective_and_falls_back_when_not(self): + # packed tissue: a few wide gaps among many membranes -> the split discriminates + packed = [{"gap": g} for g in [1] * 90 + [9] * 10] + _thr, selectivity, use_wide = gap_channel(packed) + assert use_wide and selectivity <= 0.25 + # sparse tissue: most edges face open background -> the gap says nothing, use all edges + sparse = [{"gap": g} for g in [0] * 50 + [8] * 50] + _thr, selectivity, use_wide = gap_channel(sparse) + assert not use_wide and selectivity > 0.25 + + def test_significance_threshold_scales_with_the_null_rate(self): + # denser scatter -> a larger count is needed before a peak is surprising + sparse_null = _min_significant_count(n_edges=100, n_bins=300, window=1, alpha=0.01) + dense_null = _min_significant_count(n_edges=3000, n_bins=300, window=1, alpha=0.01) + assert dense_null > sparse_null >= 3 + + def test_detect_and_flag_roundtrip(self): + rng = np.random.default_rng(0) + # a seam at x=100 (wide gaps) over scattered facet background (small gaps) + edges = [] + for c in 100.0 + rng.normal(0, 1, 30): + edges.append({"axis": "v", "coord": float(c), "span": 20, "side": -1, "gap": 9.0, "cell_id": -1}) + for c in rng.uniform(0, 200, 80): + edges.append({"axis": "v", "coord": float(c), "span": 20, "side": -1, "gap": 5.0, "cell_id": -1}) + scale = SeamDetectionParams()._resolve(20.0) + seams = detect_seams(edges, 200, 200, scale) + assert len(seams["v"]) == 1 and abs(seams["v"][0][0] - 100) <= 4 + # a left-body cell whose right edge lands on the seam is flagged; one far away is not + ebc = { + 1: [{"axis": "v", "coord": 100.0, "span": 20, "side": -1, "gap": 9.0}], + 2: [{"axis": "v", "coord": 20.0, "span": 20, "side": -1, "gap": 9.0}], + } + flagged = flag_cells_on_seams(ebc, seams, scale) + assert 1 in flagged and 2 not in flagged + + +class TestTileGridIndependence: + """Seam detection must depend on the data only, never on the QC tiling used to process it. + + ``tile_size`` is a throughput knob: it controls how the labels raster is chopped up for + parallel scoring. Two runs of the same data at different ``tile_size`` must therefore + agree on where the FOV seams are. + """ + + @pytest.mark.parametrize("tile_size", [128, 200, 420]) + def test_seams_match_truth_at_any_tile_size(self, sdata_dense_seam, tile_size): + sdata, gt = sdata_dense_seam + adata = calculate_tiling_qc(sdata, labels_key="labels", tile_size=tile_size, detect_seams=True, inplace=False) + _assert_seams_match(adata.uns["tiling_qc"]["seams"], gt.seam_coords, f"at tile_size={tile_size}") + + def test_no_seam_lands_on_the_processing_tile_grid(self, sdata_dense_seam): + """The QC tile borders are not seams; detecting one there is a processing artifact.""" + sdata, gt = sdata_dense_seam + tile_size = 200 + adata = calculate_tiling_qc(sdata, labels_key="labels", tile_size=tile_size, detect_seams=True, inplace=False) + d = adata.uns["tiling_qc"]["seam_diameter"] + extent = sdata.labels["labels"].shape[-1] + borders = [tile_size * k for k in range(1, 1 + extent // tile_size)] + seams = adata.uns["tiling_qc"]["seams"] + for axis in ("v", "h"): + for band in seams[axis]: + c = band["coord"] + near_true = min(abs(c - t) for t in gt.seam_coords) + near_border = min(abs(c - b) for b in borders) + assert near_true <= d or near_border > d, ( + f"axis {axis}: seam at {c:.1f} sits on the tile grid {borders} " + f"but not on a true seam {list(gt.seam_coords)}" + ) + + +class TestSparseTissueDetection: + """Detection must not depend on tissue density. + + In sparse tissue nearly every cell edge faces open background, so a wide-gap pre-filter + keeps almost every edge and discriminates nothing. The seam must then be found from the + alignment consensus alone -- many cells sharing one edge coordinate -- which is the + evidence that defines a seam in the first place. + """ + + @pytest.mark.parametrize("tile_size", [150, 200]) + def test_seams_found_in_sparse_tissue(self, sdata_tile_boundary, tile_size): + sdata, gt = sdata_tile_boundary + adata = calculate_tiling_qc(sdata, labels_key="labels", tile_size=tile_size, detect_seams=True, inplace=False) + seams = adata.uns["tiling_qc"]["seams"] + for axis, truth in (("v", gt.tile_borders_x), ("h", gt.tile_borders_y)): + coords = sorted(b["coord"] for b in seams[axis]) + for want in truth: + assert any(abs(c - want) <= 8 for c in coords), ( + f"axis {axis} at tile_size={tile_size}: true seam {want} not found " + f"in {[round(c, 1) for c in coords]}" + ) + + def test_seams_match_truth_at_any_overlap_margin(self, sdata_dense_seam): + """Gap probing reaches beyond a cell's own edge, so the crop must leave room for it. + + ``overlap_margin`` is sized to *contain* each owned cell; a cell sitting at the crop + edge would have its probe clipped and under-report its gap. Detection must not + depend on that either. + """ + sdata, gt = sdata_dense_seam + adata = calculate_tiling_qc( + sdata, labels_key="labels", tile_size=200, overlap_margin=2, detect_seams=True, inplace=False + ) + _assert_seams_match(adata.uns["tiling_qc"]["seams"], gt.seam_coords, "with overlap_margin=2") diff --git a/tests/experimental/test_tiling_stitch.py b/tests/experimental/test_tiling_stitch.py index c55b65548..6c15191b2 100644 --- a/tests/experimental/test_tiling_stitch.py +++ b/tests/experimental/test_tiling_stitch.py @@ -11,6 +11,7 @@ from spatialdata.models import Labels2DModel import squidpy as sq +from squidpy.experimental.tl import SeamDetectionParams from tests.conftest import DPI, PlotTester, PlotTesterMeta @@ -30,17 +31,19 @@ def test_columns_present(self, sdata_tile_boundary): assert col in adata.obs.columns def test_confidence_convention(self, sdata_tile_boundary): - # NaN = not evaluated (non-outlier), 1.0 = solo outlier, composite = stitched. + # NaN = not evaluated (non-candidate), 1.0 = solo candidate, composite = stitched. + # The candidate gate defaults to `is_seam_cut` when present, else `is_outlier`. sdata, _ = sdata_tile_boundary obs = _run_qc_and_stitch(sdata, min_confidence=0.5).obs + gate = "is_seam_cut" if "is_seam_cut" in obs.columns else "is_outlier" - non_outliers = ~obs["is_outlier"].astype(bool) - assert non_outliers.sum() > 0 - assert obs.loc[non_outliers, "stitch_confidence"].isna().all() - assert (obs.loc[non_outliers, "stitch_group_id"] == obs.loc[non_outliers, "label_id"]).all() - assert (obs.loc[non_outliers, "n_pieces"] == 1).all() + non_cands = ~obs[gate].astype(bool) + assert non_cands.sum() > 0 + assert obs.loc[non_cands, "stitch_confidence"].isna().all() + assert (obs.loc[non_cands, "stitch_group_id"] == obs.loc[non_cands, "label_id"]).all() + assert (obs.loc[non_cands, "n_pieces"] == 1).all() - solo = obs["is_outlier"].astype(bool) & ~obs["is_stitched"].astype(bool) + solo = obs[gate].astype(bool) & ~obs["is_stitched"].astype(bool) if solo.sum() > 0: assert (obs.loc[solo, "stitch_confidence"] == 1.0).all() @@ -93,9 +96,8 @@ def test_recovery_meets_quantitative_bounds(self, sdata_tile_boundary): def test_uns_records_params_and_features(self, sdata_tile_boundary): sdata, _ = sdata_tile_boundary - meta = _run_qc_and_stitch(sdata, min_confidence=0.7, max_gap=4.0).uns["tiling_stitch"] + meta = _run_qc_and_stitch(sdata, min_confidence=0.7).uns["tiling_stitch"] assert meta["min_confidence"] == 0.7 - assert meta["max_gap"] == 4.0 assert isinstance(meta["stitch_params"], dict) assert "model_coefficients" not in meta and "model_intercept" not in meta assert set(meta["score_features"]) == { @@ -103,7 +105,6 @@ def test_uns_records_params_and_features(self, sdata_tile_boundary): "endpoint_match", "merge_compactness", "merge_solidity", - "gap_proximity", } @pytest.mark.parametrize( @@ -170,6 +171,86 @@ def test_obs_and_uns_survive_zarr_roundtrip(self, sdata_tile_boundary, tmp_path) assert "tiling_stitch" in a2.uns +class TestPairingContract: + """Function-level tests of the seam-restricted, rank-based pairing pipeline. + + These exercise the stitcher directly on a hand-built two-sided cut (given the + seam bands + data scale that ``calculate_tiling_qc`` would supply), isolating + the pairing/scoring logic from detection. A dense synthetic fixture is a poor + end-to-end vehicle here -- its uniform geometry makes seam *detection* + over-flag -- so the two-sided merge contract is locked in deterministically at + this level instead. + """ + + @staticmethod + def _two_sided_labels(): + # One cell cut into a left half (id 1) and a right half (id 2) across a + # vertical seam at x~101, plus an off-seam distractor (id 3). bboxes use + # the skimage convention (max exclusive), as _compute_outlier_bboxes returns. + H, W = 60, 220 + arr = np.zeros((H, W), dtype=np.int32) + arr[20:40, 80:99] = 1 # left half: cols 80..98 + arr[20:40, 104:123] = 2 # right half: cols 104..122 (6 px seam gap) + arr[20:40, 160:180] = 3 # distractor, far from the seam + bboxes = {1: (20, 80, 40, 99), 2: (20, 104, 40, 123), 3: (20, 160, 40, 180)} + seams = {"v": [(101.0, 6.0, 50)], "h": []} + diameter = 20.0 + return arr, bboxes, seams, diameter, H, W + + def test_facing_halves_merge_off_seam_cell_ignored(self): + from squidpy.experimental.tl import _tiling_stitch as ts + + arr, bboxes, seams, diameter, H, W = self._two_sided_labels() + scale = SeamDetectionParams()._resolve(diameter) + edges, crops = ts._extract_cut_edges(arr, [1, 2, 3], bboxes, seams, scale) + # Both halves put a cut edge on the seam; the off-seam distractor does not. + assert {e.cell_id for e in edges} == {1, 2} + + cands = ts._enumerate_pair_candidates(edges, k_neighbors=5, candidate_min_iou=0.2) + pairs = ts._score_pairs(cands, bboxes, crops, 0.6, seams, close_radius_min=2, H=H, W=W) + merged = [p for p in pairs if {p.cell_a, p.cell_b} == {1, 2}] + assert len(merged) == 1 + assert merged[0].confidence >= 0.6 + assert 3 not in {p.cell_a for p in pairs} | {p.cell_b for p in pairs} + + def test_enumeration_is_rank_based_not_absolute_gap(self): + # The two halves sit 6 px apart -- far beyond the old 3 px max_gap default. + # Rank-based (k-NN) enumeration must still surface them as a candidate. + from squidpy.experimental.tl import _tiling_stitch as ts + + arr, bboxes, seams, diameter, _H, _W = self._two_sided_labels() + edges, _ = ts._extract_cut_edges(arr, [1, 2], bboxes, seams, SeamDetectionParams()._resolve(diameter)) + cands = ts._enumerate_pair_candidates(edges, k_neighbors=5, candidate_min_iou=0.2) + pair_ids = {(min(e.cell_id, c.cell_id), max(e.cell_id, c.cell_id)) for e, c, _ in cands} + assert (1, 2) in pair_ids + + def test_pieces_farther_apart_than_the_seam_are_not_merged(self): + """A cut's two halves cannot be separated by more than the seam band they lie on. + + The closing radius is scaled to each pair's own gap, so without a bound tied to the + measured seam any two aligned blobs get bridged by a disk large enough to join them + and then score as one compact, solid cell. + """ + from squidpy.experimental.tl import _tiling_stitch as ts + + H, W = 60, 260 + arr = np.zeros((H, W), dtype=np.int32) + arr[20:40, 80:99] = 1 # left piece, cols 80..98 + arr[20:40, 127:146] = 2 # right piece, cols 127..145 -> 29 px apart + bboxes = {1: (20, 80, 40, 99), 2: (20, 127, 40, 146)} + # a ~20 px wide seam band covering both edges: the pieces are still farther apart + # from each other than the seam itself is wide, so they are not one cut cell. + seams = {"v": [(112.5, 10.0, 50)], "h": []} + diameter = 20.0 + + edges, crops = ts._extract_cut_edges(arr, [1, 2], bboxes, seams, SeamDetectionParams()._resolve(diameter)) + cands = ts._enumerate_pair_candidates(edges, k_neighbors=5, candidate_min_iou=0.2) + pairs = ts._score_pairs(cands, bboxes, crops, 0.6, seams, close_radius_min=2, H=H, W=W) + assert not [p for p in pairs if {p.cell_a, p.cell_b} == {1, 2}], ( + "pieces 29 px apart across a 20 px seam were merged" + ) + + class TestStitchVisual(PlotTester, metaclass=PlotTesterMeta): _ZOOM = (150, 250, 250, 350) _SEAM_Y = 200