bioimage_py.morphology.regionprops

  1"""Advanced per-object morphology: numpy moment features + surface area and centroid correction.
  2
  3A *second pass* on top of :func:`bioimage_py.morphology.morphology`: for each labeled object it crops
  4the sub-volume by the precomputed bounding box, masks ``== label``, and computes shape features
  5directly with numpy (in physical units via ``resolution``) — the physical volume (``area``), the
  6``extent``, the ``equivalent_diameter_area`` and the major/minor ``axis_*_length`` from the object's
  7second moments — plus an optional marching-cubes ``surface_area`` (3D) and a corrected centroid (the
  8center-of-mass if it lies inside the object, else the deepest-interior voxel via the Euclidean
  9distance transform).
 10
 11The moment features reproduce the corresponding scikit-image ``regionprops`` definitions exactly, but
 12without constructing a ``RegionProperties`` object per label (and without the expensive convex-hull
 13``solidity`` / topological ``euler_number``), so the pass scales linearly in the number of objects.
 14Work is mapped one task per object with the generic
 15:meth:`bioimage_py.runner.base.Runner.map`, so it runs identically across ``local`` / ``subprocess``
 16/ ``slurm``; for distributed backends the base table is serialized to a temp file the workers read.
 17"""
 18from __future__ import annotations
 19
 20import functools
 21import os
 22from typing import Any, Callable, Dict, List, Optional, Sequence, Union
 23
 24import bioimage_cpp as bic
 25import numpy as np
 26import pandas as pd
 27
 28from ..runner import get_runner
 29from ..runner.config import RunnerConfig
 30from ..sources import Source, SourceLike, as_source
 31from ..util import check_rerun_args
 32from .morphology import _axis_names
 33
 34__all__ = ["regionprops"]
 35
 36# Per-worker cache so each task reopens the segmentation `SourceSpec` once, not per object (keyed by
 37# the spec's stable fields). The consumed table columns travel in the closure as numpy arrays, so
 38# they need no worker-side caching.
 39_SEG_CACHE: Dict[Any, Source] = {}
 40
 41
 42def _read_table(path: str) -> "pd.DataFrame":
 43    """Read a serialized table from a ``.csv`` / ``.xlsx`` path."""
 44    lower = str(path).lower()
 45    if lower.endswith((".xlsx", ".xls")):
 46        return pd.read_excel(path)
 47    if lower.endswith(".csv"):
 48        return pd.read_csv(path)
 49    raise ValueError(f"Unsupported table format: {path!r} (use .csv or .xlsx/.xls).")
 50
 51
 52def _load_table(table: Union[str, "pd.DataFrame"]) -> "pd.DataFrame":
 53    """Return ``table`` as a DataFrame (pass-through, or read from a path)."""
 54    if isinstance(table, pd.DataFrame):
 55        return table
 56    return _read_table(str(table))
 57
 58
 59def _required_columns(axes: Sequence[str]) -> List[str]:
 60    """The base-morphology columns this op consumes."""
 61    return (["label"] + [f"com_{a}" for a in axes]
 62            + [f"bb_min_{a}" for a in axes] + [f"bb_max_{a}" for a in axes])
 63
 64
 65def _resolve_seg(seg: Union[Source, Any]) -> Source:
 66    """Return an opened segmentation source, reopening (and caching) a `SourceSpec` if needed."""
 67    if isinstance(seg, Source):
 68        return seg
 69    key = (seg.kind, seg.path, seg.internal_path)
 70    src = _SEG_CACHE.get(key)
 71    if src is None:
 72        from ..sources.dispatch import from_spec
 73        src = from_spec(seg)
 74        _SEG_CACHE[key] = src
 75    return src
 76
 77
 78def _column_arrays(df: "pd.DataFrame", axes: Sequence[str]) -> Dict[str, np.ndarray]:
 79    """Extract the consumed columns as numpy arrays (built once, indexed by row position).
 80
 81    Indexing these arrays by position avoids per-object pandas scalar access, which dominates the
 82    cost when there are many objects. Returns ``label`` (``(N,)``) and ``com`` / ``bb_min`` /
 83    ``bb_max`` (each ``(N, ndim)``).
 84    """
 85    return {
 86        "label": df["label"].to_numpy(dtype="int64"),
 87        "com": np.stack([df[f"com_{a}"].to_numpy(dtype="float64") for a in axes], axis=1),
 88        "bb_min": np.stack([df[f"bb_min_{a}"].to_numpy(dtype="int64") for a in axes], axis=1),
 89        "bb_max": np.stack([df[f"bb_max_{a}"].to_numpy(dtype="int64") for a in axes], axis=1),
 90    }
 91
 92
 93def _surface_area(mask: np.ndarray, resolution: Sequence[float]) -> float:
 94    """Physical surface area (nm² etc.) of a 3D mask via marching cubes; 0.0 if degenerate."""
 95    from skimage.measure import marching_cubes, mesh_surface_area
 96    if not mask.any():
 97        return 0.0
 98    # Pad by 1 so border voxels produce closed faces; spacing makes the area physical.
 99    padded = np.pad(mask, 1).astype("float32")
100    try:
101        verts, faces, _, _ = marching_cubes(padded, level=0.5, spacing=tuple(resolution))
102        return float(mesh_surface_area(verts, faces))
103    except (RuntimeError, ValueError):
104        # Objects too thin/small for a closed surface.
105        return 0.0
106
107
108def _axis_lengths(coords: np.ndarray, spacing: np.ndarray, ndim: int) -> "tuple[float, float]":
109    """Major/minor axis lengths (physical) of the inertia-equivalent ellipse/ellipsoid.
110
111    Mirrors scikit-image: the inertia tensor is ``trace(C) * I - C`` for the population covariance
112    ``C`` of the (physical) voxel coordinates; the axis lengths derive from its eigenvalues (sorted
113    descending). Returns ``(0.0, 0.0)`` for an empty object.
114    """
115    n = coords.shape[0]
116    if n == 0:
117        return 0.0, 0.0
118    x = coords.astype("float64") * spacing
119    xc = x - x.mean(axis=0)
120    cov = (xc.T @ xc) / n  # population covariance == central second moments / area.
121    inertia = np.trace(cov) * np.eye(ndim) - cov
122    ev = np.sort(np.clip(np.linalg.eigvalsh(inertia), 0.0, None))[::-1]  # descending.
123    if ndim == 3:
124        major = float(np.sqrt(max(0.0, 10.0 * (ev[0] + ev[1] - ev[2]))))
125        minor = float(np.sqrt(max(0.0, 10.0 * (-ev[0] + ev[1] + ev[2]))))
126        return major, minor
127    # 2D (and the nD fallback): the ellipse axes are 4 * sqrt(extreme eigenvalues).
128    return float(4.0 * np.sqrt(ev[0])), float(4.0 * np.sqrt(ev[-1]))
129
130
131def _corrected_centroid(mask: np.ndarray, com_global_vox: np.ndarray, origin: np.ndarray,
132                        resolution: Sequence[float]) -> np.ndarray:
133    """Return the corrected centroid in physical units.
134
135    Keeps the center-of-mass when it lands inside the object; otherwise uses the deepest-interior
136    voxel (argmax of the Euclidean distance transform), then scales to physical units.
137    """
138    spacing = np.asarray(resolution, dtype="float64")
139    if not mask.any():
140        return com_global_vox * spacing
141    com_local = com_global_vox - origin
142    rounded = np.round(com_local).astype(int)
143    inside = (np.all(rounded >= 0) and np.all(rounded < np.asarray(mask.shape))
144              and bool(mask[tuple(rounded)]))
145    if inside:
146        centroid_vox = com_global_vox
147    else:
148        # Deepest-interior point: argmax of the (anisotropic) exact EDT from bioimage_cpp.
149        dt = bic.distance.distance_transform(mask.astype("uint8"), sampling=tuple(resolution))
150        idx = np.unravel_index(int(np.argmax(dt)), mask.shape)
151        centroid_vox = np.asarray(idx, dtype="float64") + origin
152    return centroid_vox * spacing
153
154
155def _object_features(index: int, ctx: Dict[str, Any]) -> Dict[str, Any]:
156    """Compute the feature row for a single object (the per-item function handed to ``runner.map``)."""
157    seg = _resolve_seg(ctx["seg"])
158    axes, ndim = ctx["axes"], ctx["ndim"]
159    res = np.asarray(ctx["resolution"], dtype="float64")
160
161    label = int(ctx["label"][index])
162    com = ctx["com"][index]
163    bb_min = ctx["bb_min"][index]
164    bb_max = ctx["bb_max"][index]
165
166    crop = np.asarray(seg[tuple(slice(int(lo), int(hi)) for lo, hi in zip(bb_min, bb_max))])
167    mask = crop == label
168    coords = np.argwhere(mask)  # crop-local voxel coordinates of the object.
169    n_voxels = int(coords.shape[0])
170    area = n_voxels * float(np.prod(res))  # physical volume under spacing (== skimage's 'area').
171
172    bbox_voxels = int(np.prod(bb_max - bb_min))
173    major, minor = _axis_lengths(coords, res, ndim)
174    result: Dict[str, Any] = {
175        "label": label,
176        "n_voxels": n_voxels,
177        "area": area,
178        "extent": (n_voxels / bbox_voxels) if bbox_voxels > 0 else float("nan"),
179        "equivalent_diameter_area": ((2 * ndim * area / np.pi) ** (1.0 / ndim)) if area > 0 else 0.0,
180        "axis_major_length": major,
181        "axis_minor_length": minor,
182    }
183
184    if ctx["compute_surface"] and ndim == 3:
185        result["surface_area"] = _surface_area(mask, ctx["resolution"])
186
187    centroid = _corrected_centroid(mask, com, bb_min.astype("float64"), ctx["resolution"])
188    for a, ax in enumerate(axes):
189        result[f"centroid_{ax}"] = float(centroid[a])
190    for a, ax in enumerate(axes):
191        result[f"bb_min_{ax}"] = int(bb_min[a])
192        result[f"bb_max_{ax}"] = int(bb_max[a])
193    return result
194
195
196def _order_columns(axes: Sequence[str], compute_surface: bool, ndim: int) -> List[str]:
197    """The fixed output column order: identity, shape features, geometry, then surface."""
198    cols = ["label", "n_voxels", "area", "extent", "equivalent_diameter_area",
199            "axis_major_length", "axis_minor_length"]
200    cols += [f"centroid_{a}" for a in axes]
201    cols += [f"bb_min_{a}" for a in axes] + [f"bb_max_{a}" for a in axes]
202    if compute_surface and ndim == 3:
203        cols.append("surface_area")
204    return cols
205
206
207def _write_table(out: "pd.DataFrame", output_path: str) -> None:
208    """Write the result table to a ``.csv`` / ``.xlsx`` path (creating parent dirs)."""
209    parent = os.path.dirname(os.path.abspath(str(output_path)))
210    if parent:
211        os.makedirs(parent, exist_ok=True)
212    lower = str(output_path).lower()
213    if lower.endswith((".xlsx", ".xls")):
214        out.to_excel(output_path, index=False)
215    elif lower.endswith(".csv"):
216        out.to_csv(output_path, index=False)
217    else:
218        raise ValueError(f"Unsupported output format: {output_path!r} (use .csv or .xlsx/.xls).")
219
220
221def regionprops(
222    input: SourceLike,
223    table: Union[str, "pd.DataFrame"],
224    *,
225    resolution: Optional[Sequence[float]] = None,
226    compute_surface: bool = False,
227    output_path: Optional[str] = None,
228    num_workers: int = 1,
229    job_type: str = "local",
230    job_config: Optional[RunnerConfig] = None,
231    item_ids: Optional[Sequence[int]] = None,
232    resume_from: Optional[str] = None,
233    pre_cleanup: Optional[Callable[[str], None]] = None,
234    verbose: bool = True,
235) -> "pd.DataFrame":
236    """Compute per-object morphology features for a labeled volume, one task per object.
237
238    For each object listed in ``table`` (the output of :func:`morphology`), the sub-volume is cropped by
239    its bounding box, masked to the label, and described with numpy in physical units (via
240    ``resolution``): the physical volume ``area``, the ``extent`` (filled fraction of the bounding box),
241    the ``equivalent_diameter_area`` (diameter of the ball with the same volume) and the major/minor
242    ``axis_*_length`` (from the object's second moments). These reproduce the corresponding
243    scikit-image ``regionprops`` definitions exactly. Optionally a marching-cubes ``surface_area`` (3D
244    only) and a corrected centroid (the center-of-mass when it lies inside the object, otherwise the
245    deepest-interior voxel — the argmax of the Euclidean distance transform) are added.
246
247    Args:
248        input: The labeled segmentation (a numpy/zarr/n5 array or a `Source`); integer-typed. For the
249            ``subprocess``/``slurm`` backends it must be file-backed (zarr/n5).
250        table: The base morphology table — a pandas DataFrame or a path to a ``.csv`` / ``.xlsx`` file.
251            Must contain ``label``, ``com_<axis>``, ``bb_min_<axis>`` and ``bb_max_<axis>`` (``bb_max``
252            is the exclusive slice stop, as produced by :func:`morphology`).
253        resolution: Per-axis physical voxel size in array (e.g. z, y, x) order. Defaults to ones (voxel
254            units).
255        compute_surface: Whether to add a marching-cubes ``surface_area`` (3D inputs only). This is the
256            most expensive per-object step, so it is off by default.
257        output_path: Optional ``.csv`` / ``.xlsx`` path to also write the result to.
258        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed backends).
259        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
260        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
261        item_ids: Restrict processing to these item indices (rows of ``table``), e.g. to re-run
262            previously failed objects; the result then covers only those rows. An item id indexes a
263            *row* of ``table``, so the same table (same row order) must be passed as in the original
264            run. Mutually exclusive with ``resume_from``.
265        resume_from: Distributed only; the preserved temp folder of a failed run to resume and merge
266            (see ``runner.run``). The result then covers all objects (the already-completed merged
267            with the re-run ones). The recommended rerun path for regionprops. Mutually exclusive
268            with ``item_ids``.
269        pre_cleanup: Optional ``pre_cleanup(tmp_folder)`` callback invoked on the orchestrating
270            process with the job temp folder right before it is deleted (distributed backends only).
271            Use it to read out the per-task timing files under ``tmp_folder/timings/`` before cleanup.
272            Ignored for the ``local`` backend (no temp folder).
273        verbose: Whether to show a progress bar for the ``local`` backend. The distributed backends
274            always report their progress, since the name also labels their jobs.
275
276    Returns:
277        A pandas DataFrame with one row per object, sorted by ``label``: ``label``, ``n_voxels`` (raw
278        voxel count), ``area`` (physical volume), ``extent``, ``equivalent_diameter_area``,
279        ``axis_major_length``, ``axis_minor_length``, ``centroid_<axis>`` (corrected, physical),
280        ``bb_min_<axis>``/``bb_max_<axis>`` (global voxels), and ``surface_area`` (only when
281        ``compute_surface`` and the input is 3D).
282    """
283    check_rerun_args(job_type, resume_from, item_ids, subset_name="item_ids")
284    src = as_source(input)
285    if not np.issubdtype(np.dtype(src.dtype), np.integer):
286        raise ValueError(f"regionprops expects an integer label image, got dtype {src.dtype}.")
287    ndim = src.ndim
288    axes = _axis_names(ndim)
289
290    if resolution is None:
291        resolution = tuple(1.0 for _ in range(ndim))
292    else:
293        resolution = tuple(float(r) for r in resolution)
294        if len(resolution) != ndim:
295            raise ValueError(f"resolution {resolution} does not match the input ndim {ndim}.")
296
297    df = _load_table(table)
298    missing = [c for c in _required_columns(axes) if c not in df.columns]
299    if missing:
300        raise ValueError(
301            f"table is missing required columns {missing}; pass the output of "
302            "bioimage_py.morphology.morphology (label / com_* / bb_min_* / bb_max_*)."
303        )
304
305    n = len(df)
306    if n == 0:
307        cols = _order_columns(axes, compute_surface, ndim)
308        return pd.DataFrame({c: pd.Series(dtype="float64") for c in cols})
309
310    seg_arg: Any = src
311    if job_type != "local":
312        try:
313            seg_arg = src.to_spec()
314        except ValueError as err:
315            raise ValueError(
316                f"Distributed regionprops requires a file-backed (zarr/n5) segmentation. {err}"
317            ) from err
318
319    # The consumed columns travel with the closure as numpy arrays (built once): for distributed
320    # backends they are cloudpickled into the single shared payload the workers read.
321    ctx = {
322        "seg": seg_arg, "resolution": resolution, "axes": tuple(axes), "ndim": ndim,
323        "compute_surface": bool(compute_surface), **_column_arrays(df, axes),
324    }
325    runner = get_runner(job_type, job_config)
326    # The local runner only shows a progress bar for a non-empty name.
327    name = "regionprops" if verbose or job_type != "local" else ""
328    results = runner.map(functools.partial(_object_features, ctx=ctx), n,
329                         item_ids=item_ids, resume_from=resume_from,
330                         num_workers=num_workers, has_return_val=True, name=name,
331                         pre_cleanup=pre_cleanup)
332
333    out = pd.DataFrame(results)
334    out = out[_order_columns(axes, compute_surface, ndim)].sort_values("label")
335    out = out.reset_index(drop=True)
336    if output_path is not None:
337        _write_table(out, output_path)
338    return out