bioimage_py.evaluation

32@dataclass(frozen=True)
33class ContingencyTable:
34    """The sparse contingency table between two segmentations.
35
36    All arrays use ``uint64`` (lossless for label ids and counts). Background (label ``0``) is kept;
37    ignoring it is a metric-level concern handled by the callers of this primitive.
38
39    Attributes:
40        pairs: The ``(N, 2)`` array of co-occurring label pairs ``[label_a, label_b]``, sorted
41            lexicographically by ``(label_a, label_b)``.
42        counts: The ``(N,)`` overlap count for each pair in `pairs`.
43        labels_a: The ``(Ka,)`` sorted unique labels present in the first segmentation.
44        sizes_a: The ``(Ka,)`` size (pixel count) of each label in `labels_a`.
45        labels_b: The ``(Kb,)`` sorted unique labels present in the second segmentation.
46        sizes_b: The ``(Kb,)`` size (pixel count) of each label in `labels_b`.
47        n_points: The total number of pixels counted (after any masking).
48    """
49
50    pairs: np.ndarray
51    counts: np.ndarray
52    labels_a: np.ndarray
53    sizes_a: np.ndarray
54    labels_b: np.ndarray
55    sizes_b: np.ndarray
56    n_points: int
57
58    def as_dicts(self) -> Tuple[Dict[int, int], Dict[int, int]]:
59        """Return the per-label sizes as ``({label_a: size}, {label_b: size})`` dictionaries.
60
61        Returns:
62            A dictionary mapping each label in the first segmentation to its size, and the analogous
63            dictionary for the second segmentation.
64        """
65        a_dict = {int(lab): int(cnt) for lab, cnt in zip(self.labels_a, self.sizes_a)}
66        b_dict = {int(lab): int(cnt) for lab, cnt in zip(self.labels_b, self.sizes_b)}
67        return a_dict, b_dict
68
69    def drop_ignore(self, ignore_a: Optional[Sequence[int]] = None,
70                    ignore_b: Optional[Sequence[int]] = None) -> "ContingencyTable":
71        """Return a copy with the voxels of the given ignore labels removed.
72
73        A pair is dropped if its A-label is in ``ignore_a`` **or** its B-label is in ``ignore_b`` (the
74        marginal sizes and ``n_points`` are recomputed over what remains). This is equivalent to
75        excluding those voxels before counting. Passing ``None`` (or an empty sequence) for a side is a
76        no-op for that side; dropping every present label yields the empty table.
77
78        Args:
79            ignore_a: Labels to ignore in the first segmentation.
80            ignore_b: Labels to ignore in the second segmentation.
81
82        Returns:
83            A new `ContingencyTable` without the ignored voxels.
84        """
85        if ignore_a is None and ignore_b is None:
86            return self
87        drop = np.zeros(self.pairs.shape[0], dtype=bool)
88        if ignore_a is not None:
89            drop |= np.isin(self.pairs[:, 0], np.asarray(list(ignore_a), dtype=self.pairs.dtype))
90        if ignore_b is not None:
91            drop |= np.isin(self.pairs[:, 1], np.asarray(list(ignore_b), dtype=self.pairs.dtype))
92        keep = ~drop
93        return _table_from_pairs(self.pairs[keep], self.counts[keep])

185def contingency_table(
186    seg_a: SourceLike,
187    seg_b: SourceLike,
188    num_workers: int = 1,
189    block_shape: Optional[Tuple[int, ...]] = None,
190    job_type: str = "local",
191    job_config: Optional[RunnerConfig] = None,
192    mask: Optional[SourceLike] = None,
193    block_ids: Optional[Sequence[int]] = None,
194    resume_from: Optional[str] = None,
195    pre_cleanup: Optional[Callable[[str], None]] = None,
196) -> ContingencyTable:
197    """Compute the contingency table (sparse overlap counts) between two segmentations.
198
199    The two segmentations are compared pixel-by-pixel and the overlap counts are accumulated
200    block-wise, so the result is exact regardless of how labels straddle block boundaries. The pairing
201    is symmetric: which input is the candidate and which is the ground truth is up to the caller.
202    Background (label ``0``) is included; ignoring labels is left to the metrics built on top of this.
203
204    Args:
205        seg_a: The first segmentation (a numpy/zarr/n5 array or a `Source`); must be integer-typed.
206        seg_b: The second segmentation; must be integer-typed and the same shape as `seg_a`.
207        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed
208            backends).
209        block_shape: Shape of the processing blocks. Defaults to the input chunk shape;
210            required for unchunked data.
211        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
212        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
213        mask: Optional binary mask; pixels outside the mask are excluded from the counts.
214        block_ids: Restrict processing to these block ids (e.g. to re-run previously failed blocks);
215            the table then reflects only those blocks.
216        resume_from: Distributed only; the preserved temp folder of a failed run to resume and
217            merge (see ``runner.run``). The returned table then covers the full volume (the
218            already-completed blocks merged with the re-run ones). Mutually exclusive with
219            ``block_ids``.
220        pre_cleanup: Optional ``pre_cleanup(tmp_folder)`` callback invoked on the orchestrating
221            process with the job temp folder right before it is deleted (distributed backends only).
222            Ignored for the ``local`` backend and for the direct (single-worker, unchunked) path.
223
224    Returns:
225        The merged `ContingencyTable`.
226    """
227    check_rerun_args(job_type, resume_from, block_ids)
228    src_a = as_source(seg_a)
229    src_b = as_source(seg_b)
230    for name, src in (("seg_a", src_a), ("seg_b", src_b)):
231        if not np.issubdtype(np.dtype(src.dtype), np.integer):
232            raise ValueError(f"contingency_table expects integer label images, got dtype {src.dtype} "
233                             f"for {name}.")
234
235    if check_direct(job_type, num_workers, block_shape, mask, block_ids):
236        table = _overlap_rows(src_a[full_roi(src_a.ndim)], src_b[full_roi(src_b.ndim)])
237        tables = [table] if table is not None else []
238    else:
239        runner = get_runner(job_type, job_config)
240        results = runner.run(_contingency_block, [seg_a, seg_b], num_workers=num_workers,
241                             block_shape=block_shape, mask=mask, block_ids=block_ids,
242                             resume_from=resume_from, has_return_val=True, name="contingency_table",
243                             pre_cleanup=pre_cleanup)
244        tables = [r for r in results if r is not None]
245
246    return _merge_tables(tables)

 89def variation_of_information(
 90    segmentation: SourceLike,
 91    groundtruth: SourceLike,
 92    *,
 93    ignore_seg: Optional[Sequence[int]] = None,
 94    ignore_gt: Optional[Sequence[int]] = None,
 95    use_log2: bool = True,
 96    num_workers: int = 1,
 97    block_shape: Optional[Tuple[int, ...]] = None,
 98    job_type: str = "local",
 99    job_config: Optional[RunnerConfig] = None,
100    mask: Optional[SourceLike] = None,
101) -> Tuple[float, float]:
102    """Compute the split and merge variation of information between two segmentations.
103
104    Args:
105        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
106        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
107        ignore_seg: Labels to ignore in the segmentation (their voxels are excluded).
108        ignore_gt: Labels to ignore in the groundtruth (their voxels are excluded).
109        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats).
110        num_workers: Number of parallel workers used to build the contingency table.
111        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
112        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
113        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
114        mask: Optional binary mask; voxels outside the mask are excluded.
115
116    Returns:
117        The split variation of information and the merge variation of information.
118    """
119    table = build_table(segmentation, groundtruth, ignore_seg=ignore_seg, ignore_gt=ignore_gt,
120                        num_workers=num_workers, block_shape=block_shape, job_type=job_type,
121                        job_config=job_config, mask=mask)
122    return vi_scores(table, use_log2=use_log2)

31def vi_scores(table: ContingencyTable, *, use_log2: bool = True) -> Tuple[float, float]:
32    """Compute the split and merge variation of information from a contingency table.
33
34    Args:
35        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
36        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats).
37
38    Returns:
39        The split variation of information (``H(seg | gt)``) and the merge variation of information
40        (``H(gt | seg)``).
41    """
42    n = table.n_points
43    if n == 0:
44        return 0.0, 0.0
45    log = np.log2 if use_log2 else np.log
46    counts = table.counts.astype("float64")
47    pa = table.sizes_a.astype("float64") / n
48    pb = table.sizes_b.astype("float64") / n
49    h_a = -np.sum(pa * log(pa))
50    h_b = -np.sum(pb * log(pb))
51    sa, sb = _pair_sizes(table)
52    mutual = np.sum(counts / n * log(n * counts / (sa * sb)))
53    return float(h_a - mutual), float(h_b - mutual)

125def object_vi(
126    segmentation: SourceLike,
127    groundtruth: SourceLike,
128    *,
129    ignore_seg: Optional[Sequence[int]] = None,
130    ignore_gt: Optional[Sequence[int]] = None,
131    use_log2: bool = True,
132    num_workers: int = 1,
133    block_shape: Optional[Tuple[int, ...]] = None,
134    job_type: str = "local",
135    job_config: Optional[RunnerConfig] = None,
136    mask: Optional[SourceLike] = None,
137) -> "pd.DataFrame":
138    """Compute the per-groundtruth-object variation of information between two segmentations.
139
140    Args:
141        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
142        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
143        ignore_seg: Labels to ignore in the segmentation (their voxels are excluded).
144        ignore_gt: Labels to ignore in the groundtruth (their voxels are excluded).
145        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats).
146        num_workers: Number of parallel workers used to build the contingency table.
147        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
148        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
149        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
150        mask: Optional binary mask; voxels outside the mask are excluded.
151
152    Returns:
153        A pandas DataFrame with one row per groundtruth object (columns ``label``, ``vi_split``,
154        ``vi_merge``).
155    """
156    table = build_table(segmentation, groundtruth, ignore_seg=ignore_seg, ignore_gt=ignore_gt,
157                        num_workers=num_workers, block_shape=block_shape, job_type=job_type,
158                        job_config=job_config, mask=mask)
159    return object_vi_scores(table, use_log2=use_log2)

56def object_vi_scores(table: ContingencyTable, *, use_log2: bool = True) -> "pd.DataFrame":
57    """Compute the per-groundtruth-object variation of information from a contingency table.
58
59    Based on https://arxiv.org/pdf/1708.02599.pdf (page 16).
60
61    Args:
62        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
63        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats).
64
65    Returns:
66        A pandas DataFrame with one row per groundtruth object, sorted by label, with columns
67        ``label`` (groundtruth id), ``vi_split`` and ``vi_merge``.
68    """
69    if table.pairs.shape[0] == 0:
70        return pd.DataFrame({"label": pd.Series(dtype="uint64"),
71                             "vi_split": pd.Series(dtype="float64"),
72                             "vi_merge": pd.Series(dtype="float64")})
73    log = np.log2 if use_log2 else np.log
74    counts = table.counts.astype("float64")
75    sa, sb = _pair_sizes(table)
76
77    # Group the pairs by their groundtruth (B) label.
78    order = np.argsort(table.pairs[:, 1], kind="stable")
79    b_sorted = table.pairs[:, 1][order]
80    c, sa_o, sb_o = counts[order], sa[order], sb[order]
81    starts = np.flatnonzero(np.concatenate(([True], b_sorted[1:] != b_sorted[:-1])))
82
83    vi_merge = np.add.reduceat(-(c / sb_o) * log(c / sb_o), starts)
84    vi_split = np.add.reduceat(-(c / sb_o) * log(c / sa_o), starts)
85    return pd.DataFrame({"label": b_sorted[starts].astype("uint64"),
86                         "vi_split": vi_split, "vi_merge": vi_merge}).reset_index(drop=True)

44def rand_index(
45    segmentation: SourceLike,
46    groundtruth: SourceLike,
47    *,
48    ignore_seg: Optional[Sequence[int]] = None,
49    ignore_gt: Optional[Sequence[int]] = None,
50    num_workers: int = 1,
51    block_shape: Optional[Tuple[int, ...]] = None,
52    job_type: str = "local",
53    job_config: Optional[RunnerConfig] = None,
54    mask: Optional[SourceLike] = None,
55) -> Tuple[float, float]:
56    """Compute the adapted rand error and the rand index between two segmentations.
57
58    Args:
59        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
60        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
61        ignore_seg: Labels to ignore in the segmentation (their voxels are excluded).
62        ignore_gt: Labels to ignore in the groundtruth (their voxels are excluded).
63        num_workers: Number of parallel workers used to build the contingency table.
64        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
65        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
66        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
67        mask: Optional binary mask; voxels outside the mask are excluded.
68
69    Returns:
70        The adapted rand error and the rand index.
71    """
72    table = build_table(segmentation, groundtruth, ignore_seg=ignore_seg, ignore_gt=ignore_gt,
73                        num_workers=num_workers, block_shape=block_shape, job_type=job_type,
74                        job_config=job_config, mask=mask)
75    return rand_scores(table)

21def rand_scores(table: ContingencyTable) -> Tuple[float, float]:
22    """Compute the adapted rand error and the rand index from a contingency table.
23
24    Args:
25        table: A contingency table for the two segmentations.
26
27    Returns:
28        The adapted rand error and the rand index.
29    """
30    n = table.n_points
31    if n == 0:
32        return 0.0, 1.0
33    sum_a2 = float(np.sum(table.sizes_a.astype("float64") ** 2))
34    sum_b2 = float(np.sum(table.sizes_b.astype("float64") ** 2))
35    sum_ab2 = float(np.sum(table.counts.astype("float64") ** 2))
36
37    precision = sum_ab2 / sum_a2
38    recall = sum_ab2 / sum_b2
39    adapted_rand_error = 1.0 - (2 * precision * recall) / (precision + recall)
40    rand_index = 1.0 - (sum_a2 + sum_b2 - 2 * sum_ab2) / (n * n)
41    return float(adapted_rand_error), float(rand_index)

40def cremi_score(
41    segmentation: SourceLike,
42    groundtruth: SourceLike,
43    *,
44    ignore_seg: Optional[Sequence[int]] = None,
45    ignore_gt: Optional[Sequence[int]] = None,
46    use_log2: bool = True,
47    num_workers: int = 1,
48    block_shape: Optional[Tuple[int, ...]] = None,
49    job_type: str = "local",
50    job_config: Optional[RunnerConfig] = None,
51    mask: Optional[SourceLike] = None,
52) -> Tuple[float, float, float, float]:
53    """Compute the CREMI score between two segmentations.
54
55    Args:
56        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
57        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
58        ignore_seg: Labels to ignore in the segmentation (their voxels are excluded).
59        ignore_gt: Labels to ignore in the groundtruth (their voxels are excluded).
60        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats) for the VI part.
61        num_workers: Number of parallel workers used to build the contingency table.
62        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
63        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
64        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
65        mask: Optional binary mask; voxels outside the mask are excluded.
66
67    Returns:
68        The split variation of information, the merge variation of information, the adapted rand error,
69        and the CREMI score.
70    """
71    table = build_table(segmentation, groundtruth, ignore_seg=ignore_seg, ignore_gt=ignore_gt,
72                        num_workers=num_workers, block_shape=block_shape, job_type=job_type,
73                        job_config=job_config, mask=mask)
74    return cremi_scores(table, use_log2=use_log2)

23def cremi_scores(table: ContingencyTable, *, use_log2: bool = True) -> Tuple[float, float, float, float]:
24    """Compute the CREMI score and its components from a single contingency table.
25
26    Args:
27        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
28        use_log2: Whether to use ``log2`` (bits) or natural ``log`` (nats) for the VI part.
29
30    Returns:
31        The split variation of information, the merge variation of information, the adapted rand error,
32        and the CREMI score (``sqrt(adapted_rand_error * (vi_split + vi_merge))``).
33    """
34    vi_split, vi_merge = vi_scores(table, use_log2=use_log2)
35    adapted_rand_error, _ = rand_scores(table)
36    cremi = float(np.sqrt(adapted_rand_error * (vi_split + vi_merge)))
37    return vi_split, vi_merge, adapted_rand_error, cremi

181def matching(
182    segmentation: SourceLike,
183    groundtruth: SourceLike,
184    *,
185    threshold: float = 0.5,
186    criterion: str = "iou",
187    ignore_label: Optional[int] = 0,
188    num_workers: int = 1,
189    block_shape: Optional[Tuple[int, ...]] = None,
190    job_type: str = "local",
191    job_config: Optional[RunnerConfig] = None,
192    mask: Optional[SourceLike] = None,
193) -> Dict[str, float]:
194    """Compute object-matching scores between two segmentations.
195
196    Args:
197        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
198        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
199        threshold: Overlap threshold for a match.
200        criterion: Matching criterion, one of ``"iou"``, ``"iot"`` or ``"iop"``.
201        ignore_label: Object label removed from both axes before matching (e.g. background). ``None``
202            keeps all objects.
203        num_workers: Number of parallel workers used to build the contingency table.
204        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
205        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
206        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
207        mask: Optional binary mask; voxels outside the mask are excluded.
208
209    Returns:
210        A mapping with keys ``precision``, ``recall``, ``segmentation_accuracy`` and ``f1``.
211    """
212    table = build_table(segmentation, groundtruth, num_workers=num_workers, block_shape=block_shape,
213                        job_type=job_type, job_config=job_config, mask=mask)
214    return matching_scores(table, threshold=threshold, criterion=criterion, ignore_label=ignore_label)

130def matching_scores(table: ContingencyTable, *, threshold: float = 0.5, criterion: str = "iou",
131                    ignore_label: Optional[int] = 0) -> Dict[str, float]:
132    """Compute object-matching scores from a contingency table.
133
134    Args:
135        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
136        threshold: Overlap threshold for a match.
137        criterion: Matching criterion, one of ``"iou"``, ``"iot"`` or ``"iop"``.
138        ignore_label: Object label removed from both axes before matching (e.g. background). ``None``
139            keeps all objects.
140
141    Returns:
142        A mapping with keys ``precision``, ``recall``, ``segmentation_accuracy`` and ``f1``.
143    """
144    n_true, n_matched, n_pred, scores = _compute_scores(table, criterion, ignore_label)
145    tp = _compute_tps(scores, n_matched, threshold)
146    fp, fn = n_pred - tp, n_true - tp
147    return {"precision": _precision(tp, fp, fn), "recall": _recall(tp, fp, fn),
148            "segmentation_accuracy": _segmentation_accuracy(tp, fp, fn), "f1": _f1(tp, fp, fn)}

217def mean_segmentation_accuracy(
218    segmentation: SourceLike,
219    groundtruth: SourceLike,
220    *,
221    thresholds: Optional[Sequence[float]] = None,
222    ignore_label: Optional[int] = 0,
223    return_accuracies: bool = False,
224    num_workers: int = 1,
225    block_shape: Optional[Tuple[int, ...]] = None,
226    job_type: str = "local",
227    job_config: Optional[RunnerConfig] = None,
228    mask: Optional[SourceLike] = None,
229) -> Union[float, Tuple[float, np.ndarray]]:
230    """Compute the mean segmentation accuracy between two segmentations.
231
232    Args:
233        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
234        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
235        thresholds: IoU thresholds to average over; defaults to ``np.arange(0.5, 1.0, 0.05)``.
236        ignore_label: Object label removed from both axes before matching. ``None`` keeps all objects.
237        return_accuracies: Whether to also return the per-threshold accuracies.
238        num_workers: Number of parallel workers used to build the contingency table.
239        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
240        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
241        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
242        mask: Optional binary mask; voxels outside the mask are excluded.
243
244    Returns:
245        The mean segmentation accuracy, and (only if ``return_accuracies``) the per-threshold accuracies.
246    """
247    table = build_table(segmentation, groundtruth, num_workers=num_workers, block_shape=block_shape,
248                        job_type=job_type, job_config=job_config, mask=mask)
249    return mean_segmentation_accuracy_scores(table, thresholds=thresholds, ignore_label=ignore_label,
250                                             return_accuracies=return_accuracies)

151def mean_segmentation_accuracy_scores(
152    table: ContingencyTable,
153    *,
154    thresholds: Optional[Sequence[float]] = None,
155    ignore_label: Optional[int] = 0,
156    return_accuracies: bool = False,
157) -> Union[float, Tuple[float, np.ndarray]]:
158    """Compute the mean segmentation accuracy (DSB-2018 style) from a contingency table.
159
160    Args:
161        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
162        thresholds: IoU thresholds to average over; defaults to ``np.arange(0.5, 1.0, 0.05)``.
163        ignore_label: Object label removed from both axes before matching. ``None`` keeps all objects.
164        return_accuracies: Whether to also return the per-threshold accuracies.
165
166    Returns:
167        The mean segmentation accuracy, and (only if ``return_accuracies``) the per-threshold accuracies.
168    """
169    n_true, n_matched, n_pred, scores = _compute_scores(table, "iou", ignore_label)
170    thresholds = np.arange(0.5, 1.0, 0.05) if thresholds is None else np.asarray(thresholds, "float64")
171    accuracies = np.array([
172        _segmentation_accuracy(tp, n_pred - tp, n_true - tp)
173        for tp in (_compute_tps(scores, n_matched, float(t)) for t in thresholds)
174    ])
175    mean_accuracy = float(np.mean(accuracies)) if accuracies.size else 0.0
176    if return_accuracies:
177        return mean_accuracy, accuracies
178    return mean_accuracy

126def centroid_matching(
127    segmentation: SourceLike,
128    groundtruth: SourceLike,
129    *,
130    distance_threshold: float,
131    resolution: Optional[Sequence[float]] = None,
132    ignore_label: Optional[int] = 0,
133    num_workers: int = 1,
134    block_shape: Optional[Tuple[int, ...]] = None,
135    job_type: str = "local",
136    job_config: Optional[RunnerConfig] = None,
137    mask: Optional[SourceLike] = None,
138) -> Dict[str, float]:
139    """Compute object-matching scores between two segmentations by centroid distance.
140
141    Derives each object's center of mass with :func:`bioimage_py.morphology.morphology` (the centroid
142    extraction runs block-wise / distributed via the runner arguments) and then matches the two centroid
143    sets with :func:`coordinate_matching`.
144
145    Args:
146        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`); must be
147            integer-typed.
148        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
149        distance_threshold: Maximum centroid distance (in voxels, or physical units if ``resolution``
150            is given) for two objects to count as a match.
151        resolution: Optional per-axis voxel spacing; when given, centroids are scaled by it before
152            distances are computed, so ``distance_threshold`` is interpreted in physical units.
153        ignore_label: Object label removed from both segmentations before matching (e.g. background).
154            ``None`` keeps all objects. Note ``morphology`` already excludes label ``0``.
155        num_workers: Number of parallel workers used to extract the centroids.
156        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
157        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
158        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
159        mask: Optional binary mask; voxels outside the mask are excluded.
160
161    Returns:
162        A mapping with keys ``precision``, ``recall``, ``segmentation_accuracy`` and ``f1``.
163    """
164    coords_a = _centroids(segmentation, ignore_label, num_workers, block_shape, job_type, job_config,
165                          mask)
166    coords_b = _centroids(groundtruth, ignore_label, num_workers, block_shape, job_type, job_config,
167                          mask)
168    return coordinate_matching(coords_a, coords_b, distance_threshold=distance_threshold,
169                               resolution=resolution)

 76def coordinate_matching(
 77    coordinates_a: object,
 78    coordinates_b: object,
 79    *,
 80    distance_threshold: float,
 81    resolution: Optional[Sequence[float]] = None,
 82) -> Dict[str, float]:
 83    """Match two point sets by centroid distance and compute the matching scores.
 84
 85    Each point in ``coordinates_a`` is matched to at most one point in ``coordinates_b`` (and vice
 86    versa) via the optimal one-to-one assignment, counting a pair as a match only if their Euclidean
 87    distance does not exceed ``distance_threshold``. ``coordinates_a`` plays the role of the prediction
 88    and ``coordinates_b`` the reference, so precision is computed over ``coordinates_a`` and recall over
 89    ``coordinates_b`` (matching the orientation of :func:`matching`).
 90
 91    Args:
 92        coordinates_a: The candidate points to evaluate, an array-like of shape ``(N, ndim)``.
 93        coordinates_b: The reference points, an array-like of shape ``(M, ndim)`` with the same
 94            ``ndim`` as ``coordinates_a``.
 95        distance_threshold: Maximum centroid distance (in the units of the coordinates, or of
 96            ``resolution`` if given) for two points to count as a match.
 97        resolution: Optional per-axis spacing; when given, coordinates are scaled by it before
 98            distances are computed, so ``distance_threshold`` is interpreted in physical units.
 99
100    Returns:
101        A mapping with keys ``precision``, ``recall``, ``segmentation_accuracy`` and ``f1``.
102    """
103    coords_a = _as_coordinates(coordinates_a, "coordinates_a")
104    coords_b = _as_coordinates(coordinates_b, "coordinates_b")
105    n_a, n_b = coords_a.shape[0], coords_b.shape[0]
106
107    if n_a and n_b:
108        if coords_a.shape[1] != coords_b.shape[1]:
109            raise ValueError(f"coordinates_a and coordinates_b must have the same ndim, got "
110                             f"{coords_a.shape[1]} and {coords_b.shape[1]}.")
111        if resolution is not None:
112            spacing = np.asarray(resolution, dtype="float64")
113            if spacing.shape != (coords_a.shape[1],):
114                raise ValueError(f"resolution must have one entry per axis ({coords_a.shape[1]}), got "
115                                 f"shape {spacing.shape}.")
116            coords_a, coords_b = coords_a * spacing, coords_b * spacing
117        tp = _count_matches(cdist(coords_a, coords_b), distance_threshold)
118    else:
119        tp = 0
120
121    fp, fn = n_a - tp, n_b - tp
122    return {"precision": _precision(tp, fp, fn), "recall": _recall(tp, fp, fn),
123            "segmentation_accuracy": _segmentation_accuracy(tp, fp, fn), "f1": _f1(tp, fp, fn)}

 60def dice_score(
 61    segmentation: SourceLike,
 62    groundtruth: SourceLike,
 63    *,
 64    threshold_seg: Optional[float] = 0,
 65    threshold_gt: Optional[float] = 0,
 66    num_workers: int = 1,
 67    block_shape: Optional[Tuple[int, ...]] = None,
 68    job_type: str = "local",
 69    job_config: Optional[RunnerConfig] = None,
 70    mask: Optional[SourceLike] = None,
 71) -> float:
 72    """Compute the dice score between a (binarized) segmentation and groundtruth.
 73
 74    To compare probability maps (values in ``[0, 1]``) pass ``threshold_seg=None`` /
 75    ``threshold_gt=None`` (a soft dice on the raw values); otherwise each input is binarized at its
 76    threshold.
 77
 78    Args:
 79        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
 80        groundtruth: The groundtruth; same shape as ``segmentation``.
 81        threshold_seg: Threshold applied to the segmentation, or ``None`` to use it as-is.
 82        threshold_gt: Threshold applied to the groundtruth, or ``None`` to use it as-is.
 83        num_workers: Number of parallel workers.
 84        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
 85        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
 86        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
 87        mask: Optional binary mask; voxels outside the mask are excluded.
 88
 89    Returns:
 90        The dice score.
 91    """
 92    if check_direct(job_type, num_workers, block_shape, mask, None):
 93        src_a, src_b = as_source(segmentation), as_source(groundtruth)
 94        intersection, sum_a, sum_b = _dice_sums(src_a[full_roi(src_a.ndim)], src_b[full_roi(src_b.ndim)],
 95                                                threshold_seg, threshold_gt)
 96    else:
 97        runner = get_runner(job_type, job_config)
 98        results = runner.run(_make_dice_compute(threshold_seg, threshold_gt),
 99                             [segmentation, groundtruth], num_workers=num_workers,
100                             block_shape=block_shape, mask=mask, has_return_val=True,
101                             name="dice_score")
102        results = [r for r in results if r is not None]
103        if not results:
104            return 0.0
105        intersection, sum_a, sum_b = np.array(results, dtype="float64").sum(axis=0)
106    return float(2.0 * intersection) / float(sum_a + sum_b + _EPS)

154def symmetric_best_dice_score(
155    segmentation: SourceLike,
156    groundtruth: SourceLike,
157    *,
158    ignore_label: Optional[int] = 0,
159    num_workers: int = 1,
160    block_shape: Optional[Tuple[int, ...]] = None,
161    job_type: str = "local",
162    job_config: Optional[RunnerConfig] = None,
163    mask: Optional[SourceLike] = None,
164) -> float:
165    """Compute the symmetric best dice score between two instance segmentations.
166
167    This metric is used in the CVPPP instance segmentation challenge.
168
169    Args:
170        segmentation: Candidate segmentation to evaluate (a numpy/zarr/n5 array or a `Source`).
171        groundtruth: The groundtruth segmentation; same shape as ``segmentation``.
172        ignore_label: Label excluded as an object on both sides (e.g. background). ``None`` keeps all.
173        num_workers: Number of parallel workers used to build the contingency table.
174        block_shape: Shape of the processing blocks. Defaults to the input chunk shape.
175        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
176        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
177        mask: Optional binary mask; voxels outside the mask are excluded.
178
179    Returns:
180        The symmetric best dice score.
181    """
182    table = build_table(segmentation, groundtruth, num_workers=num_workers, block_shape=block_shape,
183                        job_type=job_type, job_config=job_config, mask=mask)
184    return best_dice_scores(table, ignore_label=ignore_label)

123def best_dice_scores(table: ContingencyTable, *, ignore_label: Optional[int] = 0) -> float:
124    """Compute the symmetric best dice score from a contingency table.
125
126    For each object in one segmentation the best dice with any object in the other is taken; this is
127    averaged per segmentation and the smaller of the two averages is returned.
128
129    Args:
130        table: A contingency table built as ``contingency_table(segmentation, groundtruth)``.
131        ignore_label: Label excluded as an object on both sides (e.g. background). ``None`` keeps all.
132
133    Returns:
134        The symmetric best dice score.
135    """
136    if table.pairs.shape[0] == 0:
137        return 0.0
138    idx_a = np.searchsorted(table.labels_a, table.pairs[:, 0])
139    idx_b = np.searchsorted(table.labels_b, table.pairs[:, 1])
140    size_a = table.sizes_a.astype("float64")[idx_a]
141    size_b = table.sizes_b.astype("float64")[idx_b]
142    dice_pair = 2.0 * table.counts.astype("float64") / (size_a + size_b)
143
144    if ignore_label is None:
145        valid = np.ones(table.pairs.shape[0], dtype=bool)
146    else:
147        valid = (table.pairs[:, 0] != ignore_label) & (table.pairs[:, 1] != ignore_label)
148
149    dir_seg = _best_dice_direction(table.labels_a, idx_a, dice_pair, valid, ignore_label)
150    dir_gt = _best_dice_direction(table.labels_b, idx_b, dice_pair, valid, ignore_label)
151    return min(dir_seg, dir_gt)