feat(dynacell/eval): whole-cell + nucleus instance average-precision

Add instance mAP (cubic average_precision, pred vs GT) to the eval
pipeline for two targets, gated by compute_instance_ap:

- nucleus (backend=cellpose): each side segments its own nucleus channel
  into instance labels independently.
- whole-cell (backend=cellpose_watershed, membrane): shared GT-nuclei
  Cellpose seeds -> membrane EDT watershed (segment_whole_cell, a cubic
  port of the .tmp/wholecell cellgrid prototype) -> cytoplasm-only labels.

Both share one AP wrapper (instance_metrics), a uint16 instance-label
cache (cache.py write/read_instance_mask + pipeline_cache _fov_instances
with manifest auto-invalidation + a whole-cell seed preflight), the
merged-column metric output (AP_<th>/mAP/n_gt/n_pred/instance_*@0.50 in
mask_metrics), and the _process_one_fov instance branch. 2D (default,
single slice at frac/sharpest) and 3D (full-volume do_3D) are supported.

A bidirectional guard validates the backend/target/toggle combination,
the cached-final gate requires AP columns when compute_instance_ap is
set, and the instance toggles are grouped-run invariants.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: alxndrkalinin

applications/dynacell/src/dynacell/evaluation/_configs/eval.yaml

@@ -18,6 +18,44 @@ defaults:
 benchmark: null
 
 target_name: ???
+
+# Nucleus/membrane organelle-mask segmenter. `supermodel` = segmenter-model-zoo
+# SuperModel (default, unchanged behavior). `cellpose` = GPU Cellpose-SAM + fnet
+# pre/post (nucleus only). `cellpose_watershed` = whole-cell (membrane only):
+# GPU Cellpose-SAM nuclei seeds + membrane EDT watershed. The cellpose /
+# cellpose_watershed backends with compute_instance_ap=true add instance mAP
+# (cubic average_precision) to mask_metrics; labels cache under
+# instance_masks/<target>__<backend>.zarr.
+segmentation:
+  backend: supermodel  # supermodel | cellpose (nucleus) | cellpose_watershed (whole-cell)
+  dimension: 2d        # 2d | 3d — instance backends only (3d = full-volume do_3D)
+  slice_selection: frac    # frac | sharpest — picks the 2d plane (dimension=2d)
+  slice_fraction: 0.30     # fractional Z depth when slice_selection=frac
+  nuclei_channel_name: null  # GT-plate channel for watershed seeds (cellpose_watershed only)
+  # Nucleus Cellpose-SAM params (the cellpose path + the cellpose_watershed seeds).
+  cellpose:
+    target_voxel_um: 0.58  # isotropic grid the nuclei are segmented on (always)
+    cellprob_threshold: 0.0
+    flow_threshold: 0.4
+    min_obj_size: 30       # drop nuclei < N px/vox at target_voxel_um; 2D default — raise (~500) for dimension=3d
+  # Whole-cell watershed stage (cellpose_watershed only). Physical-um params,
+  # converted to working-grid pixels per cell_voxel_um (so they are grid-independent).
+  watershed:
+    cell_voxel_um: 0.3     # isotropic cell-stage grid (null = native resolution)
+    close_um: 2.5          # grayscale-closing disk radius
+    wall_sigma_um: 0.35    # gaussian sigma before wall thresholding
+    wall_min_um: 1.0       # drop membrane-wall specks (um^ndim)
+    hole_um: 3.0           # fill cell-mask holes (um^ndim)
+    min_cell_um: 15.0      # drop cells below this (2D um^2 default; set ~50 um^3 for dimension=3d)
+    memb_clahe: true       # CLAHE the membrane channel
+    subtract_nuclei: true  # carve nuclei out -> cytoplasmic-shape-only metrics
+
+# Whole-cell / nucleus instance average-precision. Requires an instance backend:
+# (cellpose_watershed, membrane) or (cellpose, nucleus).
+compute_instance_ap: false
+instance_metrics:
+  iou_thresholds: [0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95]
+
 io:
   pred_path: ???
   gt_path: ???
@@ -110,6 +148,8 @@ force_recompute:
   pred_dinov3: false
   pred_dynaclr: false
   pred_celldino: false
+  gt_instances: false
+  pred_instances: false
   final_metrics: false
 
 save:

applications/dynacell/src/dynacell/evaluation/cache.py

@@ -45,10 +45,28 @@ class CachePaths:
     manifest: Path
     masks_dir: Path
     features_dir: Path
+    instance_masks_dir: Path
 
-    def mask_plate(self, target_name: str) -> Path:
-        """Return the HCS OME-Zarr plate for masks of *target_name*."""
-        return self.masks_dir / f"{target_name}.zarr"
+    def mask_plate(self, target_name: str, backend: str = "supermodel") -> Path:
+        """Return the HCS OME-Zarr plate for masks of *target_name*.
+
+        Non-default segmentation backends get a ``__{backend}`` filename infix so
+        masks from different segmenters never collide in one cache dir; the
+        default ``supermodel`` keeps the bare ``{target_name}.zarr`` path so
+        pre-existing caches are unaffected.
+        """
+        stem = target_name if backend == "supermodel" else f"{target_name}__{backend}"
+        return self.masks_dir / f"{stem}.zarr"
+
+    def instance_mask_plate(self, target_name: str, backend: str) -> Path:
+        """Return the HCS OME-Zarr plate for instance (uint16) labels.
+
+        Instance-label caches always carry a ``__{backend}`` infix
+        (``nucleus__cellpose.zarr`` / ``membrane__cellpose_watershed.zarr``) and
+        live under a separate directory from the binary ``organelle_masks`` so
+        the two never collide.
+        """
+        return self.instance_masks_dir / f"{target_name}__{backend}.zarr"
 
     def cp_features(self) -> Path:
         """Return the zarr group path for CP regionprops features."""
@@ -75,6 +93,7 @@ def cache_paths(cache_dir: Path | str) -> CachePaths:
         manifest=root / "manifest.yaml",
         masks_dir=root / "organelle_masks",
         features_dir=root / "features",
+        instance_masks_dir=root / "instance_masks",
     )
 
 
@@ -242,7 +261,7 @@ def built_at_now() -> str:
     return datetime.now(timezone.utc).isoformat(timespec="seconds")
 
 
-def read_mask(paths: CachePaths, target_name: str, pos_name: str) -> np.ndarray | None:
+def read_mask(paths: CachePaths, target_name: str, pos_name: str, backend: str = "supermodel") -> np.ndarray | None:
     """Read cached organelle masks for a single position.
 
     Returns
@@ -251,7 +270,7 @@ def read_mask(paths: CachePaths, target_name: str, pos_name: str) -> np.ndarray
         Bool array of shape ``(T, D, H, W)``, or ``None`` if the plate or
         position is absent.
     """
-    plate_path = paths.mask_plate(target_name)
+    plate_path = paths.mask_plate(target_name, backend)
     if not plate_path.exists():
         return None
     with open_ome_zarr(plate_path, mode="r") as plate:
@@ -297,6 +316,7 @@ def write_mask(
     masks: np.ndarray,
     *,
     channel_name: str = _MASK_CHANNEL,
+    backend: str = "supermodel",
 ) -> None:
     """Append masks for a single position to the ``{target_name}.zarr`` plate.
 
@@ -315,7 +335,7 @@ def write_mask(
     """
     if masks.ndim != 4:
         raise ValueError(f"masks must be 4-D (T, D, H, W); got shape {masks.shape}")
-    plate_path = paths.mask_plate(target_name)
+    plate_path = paths.mask_plate(target_name, backend)
     plate_path.parent.mkdir(parents=True, exist_ok=True)
     data = masks.astype(bool)[:, None]  # (T, 1, D, H, W)
     if plate_path.exists() and _is_position_malformed(plate_path, pos_name):
@@ -341,6 +361,86 @@ def write_mask(
         position.create_image("0", data)
 
 
+_INSTANCE_MASK_CHANNEL = "instance_seg"
+
+
+def read_instance_mask(paths: CachePaths, target_name: str, pos_name: str, backend: str) -> np.ndarray | None:
+    """Read cached instance labels for a single position.
+
+    Returns
+    -------
+    numpy.ndarray | None
+        uint16 array of shape ``(T, D, H, W)`` (2-D runs are stored with
+        ``D=1``), or ``None`` if the plate or position is absent.
+    """
+    plate_path = paths.instance_mask_plate(target_name, backend)
+    if not plate_path.exists():
+        return None
+    with open_ome_zarr(plate_path, mode="r") as plate:
+        try:
+            position = plate[pos_name]
+        except KeyError:
+            return None
+        data = np.asarray(position.data[:, 0]).astype(np.uint16)
+    return data
+
+
+def write_instance_mask(
+    paths: CachePaths,
+    target_name: str,
+    pos_name: str,
+    labels: np.ndarray,
+    *,
+    channel_name: str = _INSTANCE_MASK_CHANNEL,
+    backend: str,
+) -> None:
+    """Append uint16 instance labels for a single position to the instance plate.
+
+    Mirrors :func:`write_mask` but preserves integer labels (no bool coercion).
+
+    Parameters
+    ----------
+    paths
+        Cache paths.
+    target_name
+        Organelle name (mask plate filename stem, with the ``__{backend}`` infix).
+    pos_name
+        HCS position name in ``row/col/fov`` form.
+    labels
+        uint16 array of shape ``(T, D, H, W)`` (2-D runs use ``D=1``).
+    channel_name
+        OME-Zarr channel label to write for this instance plate.
+    backend
+        Segmentation backend (selects the plate filename infix).
+    """
+    if labels.ndim != 4:
+        raise ValueError(f"labels must be 4-D (T, D, H, W); got shape {labels.shape}")
+    plate_path = paths.instance_mask_plate(target_name, backend)
+    plate_path.parent.mkdir(parents=True, exist_ok=True)
+    data = labels.astype(np.uint16)[:, None]  # (T, 1, D, H, W)
+    if plate_path.exists() and _is_position_malformed(plate_path, pos_name):
+        _rewrite_inner_array(plate_path / pos_name, data)
+        return
+    mode = "r+" if plate_path.exists() else "w"
+    with open_ome_zarr(
+        plate_path,
+        mode=mode,
+        layout="hcs",
+        channel_names=[channel_name],
+        version="0.5",
+    ) as plate:
+        row, col, fov = pos_name.split("/")
+        try:
+            position = plate[pos_name]
+        except KeyError:
+            position = plate.create_position(row, col, fov)
+        try:
+            del position["0"]
+        except KeyError:
+            pass
+        position.create_image("0", data)
+
+
 def _features_group_path(
     paths: CachePaths,
     kind: FeatureKind,

applications/dynacell/src/dynacell/evaluation/instance_metrics.py

@@ -0,0 +1,87 @@
+"""Instance-segmentation average-precision metrics (Cellpose-style).
+
+Thin numpy wrapper over :func:`cubic.metrics.average_precision` (TP / (TP + FP +
+FN) per IoU threshold, the Cellpose definition). Used for both the nucleus and
+whole-cell instance-AP eval paths. The returned dict merges directly into the
+per-(FOV, t) mask-metric rows, so every per-threshold AP becomes its own
+``mask_metrics.csv`` column / ``.npy`` key.
+"""
+
+import numpy as np
+from cubic.metrics import average_precision
+
+DEFAULT_IOU_THRESHOLDS = (0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95)
+"""IoU thresholds for the AP sweep (Cellpose / StarDist standard 0.50..0.95)."""
+
+_PRIMARY_THRESHOLD = 0.50
+"""IoU threshold whose TP/FP/FN counts are persisted (the standard 0.50 operating
+point); stored under ``instance_{TP,FP,FN}@0.50`` so the threshold is explicit."""
+
+
+def _relabel_sequential(labels: np.ndarray) -> np.ndarray:
+    """Relabel an integer label image to a dense ``0, 1..K`` (background stays 0).
+
+    ``cubic.average_precision`` derives object counts from ``labels.max()``, so
+    labels must be gap-free for FP/FN to be correct. Uses ``np.unique`` with
+    ``return_inverse`` (not connected-component relabeling — disjoint pieces that
+    share an id stay one object).
+    """
+    labels = np.asarray(labels)
+    uniq, inv = np.unique(labels, return_inverse=True)
+    inv = inv.reshape(labels.shape)
+    # uniq is sorted ascending; if a 0 background exists it maps to 0, real ids to
+    # 1..K. If no 0 is present (no background), shift so ids become 1..K.
+    return inv if uniq[0] == 0 else inv + 1
+
+
+def instance_average_precision(
+    labels_pred: np.ndarray,
+    labels_gt: np.ndarray,
+    iou_thresholds=DEFAULT_IOU_THRESHOLDS,
+) -> dict:
+    """Average precision of predicted vs ground-truth instance labels.
+
+    Parameters
+    ----------
+    labels_pred, labels_gt : numpy.ndarray
+        Predicted and ground-truth integer instance-label images (same shape).
+    iou_thresholds : sequence of float
+        IoU thresholds for the AP sweep.
+
+    Returns
+    -------
+    dict
+        ``AP_<th>`` (one per threshold), ``mAP`` (their mean), ``n_gt``,
+        ``n_pred``, and ``instance_{TP,FP,FN}@0.50``. Both sides empty → all-NaN
+        AP/mAP; exactly one side empty → AP/mAP 0.0.
+    """
+    thresholds = [float(t) for t in iou_thresholds]
+    pred = _relabel_sequential(labels_pred)
+    gt = _relabel_sequential(labels_gt)
+    n_pred = int(pred.max())
+    n_gt = int(gt.max())
+
+    # cubic.average_precision returns 0 (not NaN) on empty inputs, so pre-check the
+    # degenerate cases here and call it only when both sides have objects.
+    if n_gt == 0 and n_pred == 0:
+        ap_vals = [float("nan")] * len(thresholds)
+        tp = fp = fn = float("nan")
+    elif n_gt == 0 or n_pred == 0:
+        ap_vals = [0.0] * len(thresholds)
+        tp, fp, fn = 0.0, float(n_pred), float(n_gt)
+    else:
+        ap, tp_arr, fp_arr, fn_arr = average_precision(gt, pred, thresholds)
+        ap_vals = [float(a) for a in np.atleast_1d(ap)]
+        idx = thresholds.index(_PRIMARY_THRESHOLD) if _PRIMARY_THRESHOLD in thresholds else 0
+        tp = float(np.atleast_1d(tp_arr)[idx])
+        fp = float(np.atleast_1d(fp_arr)[idx])
+        fn = float(np.atleast_1d(fn_arr)[idx])
+
+    result = {"n_gt": n_gt, "n_pred": n_pred}
+    for th, a in zip(thresholds, ap_vals):
+        result[f"AP_{th:.2f}"] = a
+    result["mAP"] = float(np.mean(ap_vals))
+    result[f"instance_TP@{_PRIMARY_THRESHOLD:.2f}"] = tp
+    result[f"instance_FP@{_PRIMARY_THRESHOLD:.2f}"] = fp
+    result[f"instance_FN@{_PRIMARY_THRESHOLD:.2f}"] = fn
+    return result

applications/dynacell/src/dynacell/evaluation/instance_metrics_test.py

@@ -0,0 +1,95 @@
+"""Tests for the instance average-precision wrapper (real cubic, CPU)."""
+
+import numpy as np
+
+from dynacell.evaluation.instance_metrics import (
+    DEFAULT_IOU_THRESHOLDS,
+    _relabel_sequential,
+    instance_average_precision,
+)
+
+
+def _two_squares(shape=(16, 16)) -> np.ndarray:
+    """A label image with two well-separated square instances (ids 1, 2)."""
+    lab = np.zeros(shape, dtype=np.uint16)
+    lab[2:6, 2:6] = 1
+    lab[10:14, 10:14] = 2
+    return lab
+
+
+def test_identical_labels_map_one():
+    """Identical label images score mAP 1.0 and AP 1.0 at every threshold."""
+    gt = _two_squares()
+    result = instance_average_precision(gt.copy(), gt)
+    assert result["mAP"] == 1.0
+    assert result["n_gt"] == 2 and result["n_pred"] == 2
+    for th in DEFAULT_IOU_THRESHOLDS:
+        assert result[f"AP_{th:.2f}"] == 1.0
+    assert result["instance_TP@0.50"] == 2.0
+    assert result["instance_FP@0.50"] == 0.0
+    assert result["instance_FN@0.50"] == 0.0
+
+
+def test_disjoint_labels_map_zero():
+    """Predicting objects nowhere near the GT scores mAP 0.0."""
+    gt = _two_squares()
+    pred = np.zeros_like(gt)
+    pred[2:6, 10:14] = 1  # overlaps neither GT square
+    result = instance_average_precision(pred, gt)
+    assert result["mAP"] == 0.0
+
+
+def test_one_side_empty_is_zero():
+    """One empty side → AP 0.0 (not NaN); FP/FN reflect the non-empty side."""
+    gt = _two_squares()
+    empty = np.zeros_like(gt)
+    miss = instance_average_precision(empty, gt)
+    assert miss["mAP"] == 0.0
+    assert miss["instance_FN@0.50"] == 2.0
+    assert miss["instance_FP@0.50"] == 0.0
+    extra = instance_average_precision(gt.copy(), empty)
+    assert extra["mAP"] == 0.0
+    assert extra["instance_FP@0.50"] == 2.0
+    assert extra["instance_FN@0.50"] == 0.0
+
+
+def test_both_empty_is_nan():
+    """Both sides empty → AP undefined (NaN), not a spurious 0.0."""
+    empty = np.zeros((16, 16), dtype=np.uint16)
+    result = instance_average_precision(empty, empty)
+    assert np.isnan(result["mAP"])
+    assert result["n_gt"] == 0 and result["n_pred"] == 0
+    for th in DEFAULT_IOU_THRESHOLDS:
+        assert np.isnan(result[f"AP_{th:.2f}"])
+
+
+def test_arg_swap_moves_fp_fn_not_map():
+    """Swapping pred/gt leaves mAP unchanged but swaps FP and FN."""
+    gt = _two_squares()
+    pred = gt.copy()
+    pred[10:14, 10:14] = 0  # pred has only 1 of the 2 objects
+    a = instance_average_precision(pred, gt)  # 1 pred, 2 gt -> 1 FN
+    b = instance_average_precision(gt, pred)  # 2 pred, 1 gt -> 1 FP
+    assert a["mAP"] == b["mAP"]
+    assert a["instance_FN@0.50"] == b["instance_FP@0.50"] == 1.0
+    assert a["instance_FP@0.50"] == b["instance_FN@0.50"] == 0.0
+
+
+def test_gapped_relabel_no_merge():
+    """``_relabel_sequential`` densifies ids without merging disjoint objects."""
+    lab = np.zeros((8, 8), dtype=np.uint16)
+    lab[0:2, 0:2] = 5
+    lab[5:7, 5:7] = 9  # gap in ids (5, 9) -> should become (1, 2)
+    out = _relabel_sequential(lab)
+    assert sorted(np.unique(out).tolist()) == [0, 1, 2]
+    # the two original objects stay distinct (no merge)
+    assert out[0, 0] != out[5, 5]
+
+
+def test_relabel_preserves_object_count_for_ap():
+    """AP object counts come from labels.max(); relabel keeps n correct."""
+    lab = np.zeros((8, 8), dtype=np.uint16)
+    lab[0:2, 0:2] = 7  # single object with a non-1 id
+    result = instance_average_precision(lab.copy(), lab)
+    assert result["n_gt"] == 1 and result["n_pred"] == 1
+    assert result["mAP"] == 1.0