Skip to content

Seam report

patchworks.seam_report(labels: Union[str, Path, 'zarr.Array'], tile_shape: Sequence[int], *, component: str = '0', min_voxels: int = 4, max_faces: int = 64, warn_ratio: float = 2.0) -> dict[str, Any]

Measure whether tile seams leave artifacts in a merged label image.

For every axis, compares the orphan rate at tile seams -- the fraction of labels touching one side of a seam with nothing continuing on the other -- against the same rate on control planes halfway through each tile. A seam rate well above the interior rate means the tiling shows: raise overlap (it should cover about one object), or try stitch="iou".

Parameters:

Name Type Description Default
labels (str, Path or Array)

The merged labels: a zarr array, or a group path plus component (default "0", a label group's full resolution).

 required
tile_shape sequence of int

The tile shape the segmentation ran with.

 required
component str

Array inside labels when it is a path.

 '0'
min_voxels int

Ignore labels with fewer voxels than this on the slice (edge grazes).

 4
max_faces int

Cap on seam faces (and as many control faces) read per axis, spread evenly over the image (default 64). Each face is a two-voxel slab of one tile's cross-section, so this bounds the I/O on a huge store: 64 z-faces of a 1024 x 1024 int32 tile read ~1 GB.

 64
warn_ratio float

Log a warning when an axis' seam rate exceeds its interior rate by this factor.

 2.0

Returns:

Type Description
dict

{"axes": {axis: {...}}, "worst_seams": [...]}. Per axis: seam_labels, seam_orphans, seam_rate, interior_rate (None when tiles are one voxel thick there, leaving no interior plane), and ratio. worst_seams lists the faces with the most orphans, with their position, for a look in the viewer.

Examples:

>>> seam_report("scan.zarr/labels/cells", (16, 1024, 1024))
{'axes': {0: {'seam_rate': 0.11, 'interior_rate': 0.10, ...}, ...}, ...}
Source code in src/patchworks/_seams.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
def seam_report(
    labels: Union[str, Path, "zarr.Array"],
    tile_shape: Sequence[int],
    *,
    component: str = "0",
    min_voxels: int = 4,
    max_faces: int = 64,
    warn_ratio: float = 2.0,
) -> dict[str, Any]:
    """Measure whether tile seams leave artifacts in a merged label image.

    For every axis, compares the **orphan rate** at tile seams -- the fraction
    of labels touching one side of a seam with nothing continuing on the
    other -- against the same rate on control planes halfway through each
    tile. A seam rate well above the interior rate means the tiling shows:
    raise ``overlap`` (it should cover about one object), or try
    ``stitch="iou"``.

    Parameters
    ----------
    labels : str, Path or zarr.Array
        The merged labels: a zarr array, or a group path plus *component*
        (default ``"0"``, a label group's full resolution).
    tile_shape : sequence of int
        The tile shape the segmentation ran with.
    component : str, optional
        Array inside *labels* when it is a path.
    min_voxels : int, optional
        Ignore labels with fewer voxels than this on the slice (edge grazes).
    max_faces : int, optional
        Cap on seam faces (and as many control faces) read per axis, spread
        evenly over the image (default 64). Each face is a two-voxel slab of
        one tile's cross-section, so this bounds the I/O on a huge store:
        64 z-faces of a 1024 x 1024 int32 tile read ~1 GB.
    warn_ratio : float, optional
        Log a warning when an axis' seam rate exceeds its interior rate by
        this factor.

    Returns
    -------
    dict
        ``{"axes": {axis: {...}}, "worst_seams": [...]}``. Per axis:
        ``seam_labels``, ``seam_orphans``, ``seam_rate``, ``interior_rate``
        (None when tiles are one voxel thick there, leaving no interior
        plane), and ``ratio``. ``worst_seams`` lists the faces with the most
        orphans, with their position, for a look in the viewer.

    Examples
    --------
    >>> seam_report("scan.zarr/labels/cells", (16, 1024, 1024))  # doctest: +SKIP
    {'axes': {0: {'seam_rate': 0.11, 'interior_rate': 0.10, ...}, ...}, ...}
    """
    arr = (
        zarr.open_group(str(labels), mode="r")[component]
        if isinstance(labels, (str, Path))
        else labels
    )
    shape = arr.shape
    tile = tuple(int(t) for t in tile_shape)
    if len(tile) != len(shape):
        raise ValueError(
            f"tile_shape {tile} has {len(tile)} axes; the labels have "
            f"{len(shape)}"
        )

    axes: dict[int, dict[str, Any]] = {}
    worst: list[dict[str, Any]] = []
    for ax, (n, t) in enumerate(zip(shape, tile)):
        positions = list(range(t, n, t))
        if not positions:
            continue
        other = [a for a in range(len(shape)) if a != ax]
        columns = list(_iproduct(*[range(0, shape[a], tile[a]) for a in other]))
        faces = [(p, c) for p in positions for c in columns]
        if len(faces) > max_faces:
            pick = np.linspace(0, len(faces) - 1, max_faces).astype(int)
            faces = [faces[i] for i in pick]

        def slab(pos: int, col: tuple[int, ...]) -> np.ndarray:
            sl: list[slice] = [slice(None)] * len(shape)
            sl[ax] = slice(pos - 1, pos + 1)
            for a, off in zip(other, col):
                sl[a] = slice(off, min(off + tile[a], shape[a]))
            return np.moveaxis(np.asarray(arr[tuple(sl)]), ax, 0)

        seam_n = seam_o = inner_n = inner_o = 0
        for pos, col in faces:
            s = slab(pos, col)
            # Both directions: an orphan on either side is a seam mark.
            n1, o1 = _orphans(s[0], s[1], min_voxels)
            n2, o2 = _orphans(s[1], s[0], min_voxels)
            seam_n += n1 + n2
            seam_o += o1 + o2
            if o1 + o2:
                worst.append(
                    {
                        "axis": ax,
                        "position": pos,
                        "offset": list(col),
                        "orphans": o1 + o2,
                        "labels": n1 + n2,
                    }
                )
            if t >= 2:
                c = slab(pos - t // 2, col)  # halfway into the tile below
                m1, p1 = _orphans(c[0], c[1], min_voxels)
                m2, p2 = _orphans(c[1], c[0], min_voxels)
                inner_n += m1 + m2
                inner_o += p1 + p2

        seam_rate = seam_o / seam_n if seam_n else 0.0
        interior = (inner_o / inner_n if inner_n else 0.0) if t >= 2 else None
        ratio = (
            seam_rate / interior
            if interior
            else (float("inf") if seam_rate and interior == 0.0 else None)
        )
        axes[ax] = {
            "seams": len(positions),
            "faces_read": len(faces),
            "seam_labels": seam_n,
            "seam_orphans": seam_o,
            "seam_rate": seam_rate,
            "interior_rate": interior,
            "ratio": ratio,
        }
        if ratio is not None and ratio > warn_ratio and seam_o:
            logger.warning(
                "seam report: axis %d seams orphan %.1f%% of labels vs "
                "%.1f%% inside tiles (%.1fx) -- the tiling shows. Raise "
                "overlap to about one object, or try stitch='iou'.",
                ax,
                100 * seam_rate,
                100 * (interior or 0.0),
                ratio,
            )
        else:
            logger.info(
                "seam report: axis %d seam orphan rate %.1f%%, interior %s",
                ax,
                100 * seam_rate,
                "n/a" if interior is None else f"{100 * interior:.1f}%",
            )

    worst.sort(key=lambda f: (-f["orphans"], f["axis"], f["position"]))
    return {"axes": axes, "worst_seams": worst[:20]}

Choosing the overlap

patchworks.suggest_overlap(image: Union[da.Array, np.ndarray], fn: Callable[[np.ndarray], np.ndarray], tile_shape: Sequence[int], *, candidates: Sequence[int] = (0, 4, 8, 16, 32, 64), region: 'tuple[slice, ...] | None' = None, crop_tiles: int = 2, target: float = 0.99, stitch: str = 'touch') -> dict[str, Any]

Smallest overlap whose tiled result matches an untiled one.

Segments a crop spanning crop_tiles tiles per axis once without tiling, then with tile_process at each candidate overlap, scoring each by :func:object_f1 against the untiled reference. Candidates are tried in increasing order and the search stops at the first that reaches target.

Parameters:

Name Type Description Default
image dask or NumPy array

The image as tile_process would see it (channel already chosen).

 required
fn callable

The segmentation function. It is run on the crop in one piece, so the crop must fit it (crop_tiles=2 means a 2x2(x2) block of tiles).

 required
tile_shape sequence of int

The tile shape the real run will use.

 required
candidates sequence of int

Overlaps to try, in voxels (applied on every axis the tile allows).

 (0, 4, 8, 16, 32, 64)
region tuple of slice

The crop to use; default a centred block of tiles. Pick one with typical objects -- an empty crop agrees at any overlap.

 None
crop_tiles int

Tiles per axis in the default crop (default 2).

 2
target float

F1 counted as "tiling makes no difference" (default 0.99).

 0.99
stitch str

Stitching mode of the real run ("touch" or "iou").

 'touch'

Returns:

Type Description
dict

{"overlap": chosen or None, "scores": {overlap: f1}, "region": crop, "reference_objects": n}. overlap is None when no candidate reached target (raise the candidates, or the crop has objects larger than any of them).

Examples:

>>> from patchworks import suggest_overlap
>>> suggest_overlap(img, fn, (16, 512, 512))["overlap"]
16
Source code in src/patchworks/_autotune.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def suggest_overlap(
    image: Union[da.Array, np.ndarray],
    fn: Callable[[np.ndarray], np.ndarray],
    tile_shape: Sequence[int],
    *,
    candidates: Sequence[int] = (0, 4, 8, 16, 32, 64),
    region: "tuple[slice, ...] | None" = None,
    crop_tiles: int = 2,
    target: float = 0.99,
    stitch: str = "touch",
) -> dict[str, Any]:
    """Smallest overlap whose tiled result matches an untiled one.

    Segments a crop spanning ``crop_tiles`` tiles per axis once without
    tiling, then with ``tile_process`` at each candidate overlap, scoring
    each by :func:`object_f1` against the untiled reference. Candidates are
    tried in increasing order and the search stops at the first that reaches
    *target*.

    Parameters
    ----------
    image : dask or NumPy array
        The image as ``tile_process`` would see it (channel already chosen).
    fn : callable
        The segmentation function. It is run on the crop in one piece, so the
        crop must fit it (``crop_tiles=2`` means a 2x2(x2) block of tiles).
    tile_shape : sequence of int
        The tile shape the real run will use.
    candidates : sequence of int, optional
        Overlaps to try, in voxels (applied on every axis the tile allows).
    region : tuple of slice, optional
        The crop to use; default a centred block of tiles. Pick one with
        typical objects -- an empty crop agrees at any overlap.
    crop_tiles : int, optional
        Tiles per axis in the default crop (default 2).
    target : float, optional
        F1 counted as "tiling makes no difference" (default 0.99).
    stitch : str, optional
        Stitching mode of the real run (``"touch"`` or ``"iou"``).

    Returns
    -------
    dict
        ``{"overlap": chosen or None, "scores": {overlap: f1}, "region":
        crop, "reference_objects": n}``. ``overlap`` is None when no
        candidate reached *target* (raise the candidates, or the crop has
        objects larger than any of them).

    Examples
    --------
    >>> from patchworks import suggest_overlap  # doctest: +SKIP
    >>> suggest_overlap(img, fn, (16, 512, 512))["overlap"]  # doctest: +SKIP
    16
    """
    from ._core import tile_process

    tile = tuple(int(t) for t in tile_shape)
    arr = image if isinstance(image, da.Array) else da.from_array(image)
    region = region or _centre_crop(arr.shape, tile, crop_tiles)
    crop = np.asarray(arr[region])
    if crop.size == 0:
        raise ValueError(f"empty crop {region}")
    reference = np.asarray(fn(crop))
    n_ref = len(np.unique(reference)) - (1 if (reference == 0).any() else 0)
    if n_ref == 0:
        logger.warning(
            "suggest_overlap: the crop %s holds no objects, so every overlap "
            "agrees; pass region= with typical objects",
            region,
        )

    scores: dict[int, float] = {}
    chosen = None
    scratch = tempfile.mkdtemp(prefix="pws_overlap_")
    try:
        lazy = da.from_array(crop, chunks=tile)
        for ov in sorted(int(c) for c in candidates):
            tiled = tile_process(
                lazy,
                fn,
                overlap=ov,
                stitch=stitch,
                write_to=f"{scratch}/ov{ov}.zarr",
                progress=False,
                log_file=False,
            )
            scores[ov] = object_f1(reference, np.asarray(tiled))
            logger.info(
                "suggest_overlap: overlap %d -> F1 %.3f vs untiled",
                ov,
                scores[ov],
            )
            if scores[ov] >= target:
                chosen = ov
                break
    finally:
        shutil.rmtree(scratch, ignore_errors=True)
    if chosen is None:
        logger.warning(
            "suggest_overlap: no candidate reached F1 %.2f (best %s); try "
            "larger overlaps",
            target,
            max(scores.items(), key=lambda kv: kv[1]) if scores else None,
        )
    return {
        "overlap": chosen,
        "scores": scores,
        "region": [(s.start, s.stop) for s in region],
        "reference_objects": int(n_ref),
    }

patchworks.object_f1(a: np.ndarray, b: np.ndarray, iou: float = 0.5) -> float

Object-level F1 between two label images at an IoU threshold.

Two objects match when their IoU exceeds iou (above 0.5 a match is necessarily one-to-one). 1.0 means every object in one has its counterpart in the other; background-only images agree perfectly.

Source code in src/patchworks/_autotune.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def object_f1(a: np.ndarray, b: np.ndarray, iou: float = 0.5) -> float:
    """Object-level F1 between two label images at an IoU threshold.

    Two objects match when their IoU exceeds *iou* (above 0.5 a match is
    necessarily one-to-one). 1.0 means every object in one has its
    counterpart in the other; background-only images agree perfectly.
    """
    from ._merge import _pair_stats

    pairs, inter, a_area, b_area = _pair_stats(a, b)
    na, nb = len(a_area), len(b_area)
    if na + nb == 0:
        return 1.0
    tp = sum(
        1
        for (ai, bi), c in zip(pairs.tolist(), inter.tolist())
        if c / (a_area[ai] + b_area[bi] - c) > iou
    )
    return 2 * tp / (na + nb)