bioimage_py.morphology

141def morphology(
142    input: SourceLike,
143    num_workers: int = 1,
144    block_shape: Optional[Tuple[int, ...]] = None,
145    job_type: str = "local",
146    job_config: Optional[RunnerConfig] = None,
147    mask: Optional[SourceLike] = None,
148    block_ids: Optional[Sequence[int]] = None,
149    resume_from: Optional[str] = None,
150    pre_cleanup: Optional[Callable[[str], None]] = None,
151) -> "pd.DataFrame":
152    """Compute per-label morphology (size, center of mass, bounding box) of a labeled volume.
153
154    Statistics are computed block-wise and merged so the result is exact regardless of how labels
155    straddle block boundaries. The background label ``0`` is excluded.
156
157    Args:
158        input: The input label image (a numpy/zarr/n5 array or a `Source`); must be integer-typed.
159        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed
160            backends).
161        block_shape: Shape of the processing blocks. Defaults to the input chunk shape;
162            required for unchunked data.
163        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
164        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
165        mask: Optional binary mask; voxels outside the mask are excluded from the computation.
166        block_ids: Restrict processing to these block ids (e.g. to re-run previously failed blocks);
167            the table then reflects only those blocks.
168        resume_from: Distributed only; the preserved temp folder of a failed run to resume and
169            merge (see ``runner.run``). The returned table then covers the full volume (the
170            already-completed blocks merged with the re-run ones). Mutually exclusive with
171            ``block_ids``.
172        pre_cleanup: Optional ``pre_cleanup(tmp_folder)`` callback invoked on the orchestrating
173            process with the job temp folder right before it is deleted (distributed backends only).
174            Use it to read out the per-task timing files under ``tmp_folder/timings/`` before cleanup.
175            Ignored for the ``local`` backend and for the direct (single-worker, unchunked) path, which
176            have no temp folder.
177
178    Returns:
179        A pandas DataFrame with one row per label, sorted by label, with columns ``label``, ``size``,
180        ``com_<axis>`` (center of mass), ``bb_min_<axis>`` and ``bb_max_<axis>``. The bounding box is
181        slice-ready: ``bb_min`` is inclusive and ``bb_max`` is the exclusive stop, so the object's box
182        is ``tuple(slice(bb_min_<axis>, bb_max_<axis>) for axis in axes)``. Axis names are ``z``/``y``/
183        ``x`` for 2D/3D data and ``axis0`` ... otherwise.
184    """
185    check_rerun_args(job_type, resume_from, block_ids)
186    src = as_source(input)
187    if not np.issubdtype(np.dtype(src.dtype), np.integer):
188        raise ValueError(f"morphology expects an integer label image, got dtype {src.dtype}.")
189    ndim = src.ndim
190
191    if check_direct(job_type, num_workers, block_shape, mask, block_ids):
192        table = _block_table(src[full_roi(src.ndim)], [0] * ndim)
193        tables = [table] if table is not None else []
194    else:
195        runner = get_runner(job_type, job_config)
196        results = runner.run(_morphology_block, [input], num_workers=num_workers,
197                             block_shape=block_shape, mask=mask, block_ids=block_ids,
198                             resume_from=resume_from, has_return_val=True, name="morphology",
199                             pre_cleanup=pre_cleanup)
200        tables = [r for r in results if r is not None]
201
202    return _to_dataframe(*_merge_tables(tables, ndim), ndim)

222def regionprops(
223    input: SourceLike,
224    table: Union[str, "pd.DataFrame"],
225    *,
226    resolution: Optional[Sequence[float]] = None,
227    compute_surface: bool = False,
228    output_path: Optional[str] = None,
229    num_workers: int = 1,
230    job_type: str = "local",
231    job_config: Optional[RunnerConfig] = None,
232    item_ids: Optional[Sequence[int]] = None,
233    resume_from: Optional[str] = None,
234    pre_cleanup: Optional[Callable[[str], None]] = None,
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
274    Returns:
275        A pandas DataFrame with one row per object, sorted by ``label``: ``label``, ``n_voxels`` (raw
276        voxel count), ``area`` (physical volume), ``extent``, ``equivalent_diameter_area``,
277        ``axis_major_length``, ``axis_minor_length``, ``centroid_<axis>`` (corrected, physical),
278        ``bb_min_<axis>``/``bb_max_<axis>`` (global voxels), and ``surface_area`` (only when
279        ``compute_surface`` and the input is 3D).
280    """
281    check_rerun_args(job_type, resume_from, item_ids, subset_name="item_ids")
282    src = as_source(input)
283    if not np.issubdtype(np.dtype(src.dtype), np.integer):
284        raise ValueError(f"regionprops expects an integer label image, got dtype {src.dtype}.")
285    ndim = src.ndim
286    axes = _axis_names(ndim)
287
288    if resolution is None:
289        resolution = tuple(1.0 for _ in range(ndim))
290    else:
291        resolution = tuple(float(r) for r in resolution)
292        if len(resolution) != ndim:
293            raise ValueError(f"resolution {resolution} does not match the input ndim {ndim}.")
294
295    df = _load_table(table)
296    missing = [c for c in _required_columns(axes) if c not in df.columns]
297    if missing:
298        raise ValueError(
299            f"table is missing required columns {missing}; pass the output of "
300            "bioimage_py.morphology.morphology (label / com_* / bb_min_* / bb_max_*)."
301        )
302
303    n = len(df)
304    if n == 0:
305        cols = _order_columns(axes, compute_surface, ndim)
306        return pd.DataFrame({c: pd.Series(dtype="float64") for c in cols})
307
308    seg_arg: Any = src
309    if job_type != "local":
310        try:
311            seg_arg = src.to_spec()
312        except ValueError as err:
313            raise ValueError(
314                f"Distributed regionprops requires a file-backed (zarr/n5) segmentation. {err}"
315            ) from err
316
317    # The consumed columns travel with the closure as numpy arrays (built once): for distributed
318    # backends they are cloudpickled into the single shared payload the workers read.
319    ctx = {
320        "seg": seg_arg, "resolution": resolution, "axes": tuple(axes), "ndim": ndim,
321        "compute_surface": bool(compute_surface), **_column_arrays(df, axes),
322    }
323    runner = get_runner(job_type, job_config)
324    results = runner.map(functools.partial(_object_features, ctx=ctx), n,
325                         item_ids=item_ids, resume_from=resume_from,
326                         num_workers=num_workers, has_return_val=True, name="regionprops",
327                         pre_cleanup=pre_cleanup)
328
329    out = pd.DataFrame(results)
330    out = out[_order_columns(axes, compute_surface, ndim)].sort_values("label")
331    out = out.reset_index(drop=True)
332    if output_path is not None:
333        _write_table(out, output_path)
334    return out

 80def distance_transform(
 81    input: SourceLike,
 82    halo: Sequence[int],
 83    *,
 84    sampling: Sampling = None,
 85    return_distances: bool = True,
 86    return_indices: bool = False,
 87    distances: Optional[SourceLike] = None,
 88    indices: Optional[SourceLike] = None,
 89    block_shape: Optional[Tuple[int, ...]] = None,
 90    job_type: str = "local",
 91    job_config: Optional[RunnerConfig] = None,
 92    num_workers: int = 1,
 93    block_ids: Optional[Sequence[int]] = None,
 94    resume_from: Optional[str] = None,
 95) -> Union[SourceLike, Tuple[SourceLike, SourceLike]]:
 96    """Compute the (halo-based) distance transform of 2D/3D data, block-wise.
 97
 98    Each block is computed over a halo-padded region via ``bioimage_cpp.distance.distance_transform``
 99    and the halo-free inner block is written back, so the result is exact only up to the block
100    boundary plus ``halo``. Single-stage, so it supports ``block_ids`` / ``resume_from``.
101
102    Args:
103        input: The input data (a numpy/zarr/n5 array or a `Source`); 2D or 3D.
104        halo: Per-axis halo enlarging each block; **required** for the block-wise path (choose it to
105            cover the maximum distance of interest). Ignored by the direct (single-block) path.
106        sampling: The per-axis voxel spacing passed to the distance transform (isotropic ``1`` by
107            default).
108        return_distances: Whether to compute the distance map.
109        return_indices: Whether to compute the index (feature) transform, an ``(ndim, *shape)`` array
110            holding, per voxel, the global coordinates of the nearest background voxel.
111        distances: Output array for the distances (``float32``, shape ``input.shape``). Optional for
112            local execution -- allocated and returned if omitted; **required** for distributed
113            execution.
114        indices: Output array for the indices (``int32``, shape ``(ndim, *input.shape)``). Optional
115            for local execution; **required** for distributed execution if ``return_indices``.
116        block_shape: Shape of the processing blocks. Defaults to the input chunk shape; required
117            for unchunked data.
118        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
119        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
120        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed
121            backends).
122        block_ids: Restrict processing to these block ids (e.g. to re-run previously failed blocks).
123            Mutually exclusive with ``resume_from``.
124        resume_from: Distributed only; the preserved temp folder of a failed run to resume (see
125            ``runner.run``). Mutually exclusive with ``block_ids``.
126
127    Returns:
128        The distances array if only ``return_distances``; the indices array if only
129        ``return_indices``; a ``(distances, indices)`` tuple if both.
130    """
131    check_rerun_args(job_type, resume_from, block_ids)
132    src = as_source(input)
133    ndim = src.ndim
134    if ndim not in (2, 3):
135        raise ValueError(f"distance_transform supports 2D or 3D data, got {ndim}D.")
136    if not (return_distances or return_indices):
137        raise ValueError("At least one of 'return_distances' or 'return_indices' must be True.")
138
139    direct = (is_direct(job_type, num_workers, block_shape)
140              and block_ids is None and resume_from is None)
141
142    dist_out = (_allocate_output(distances, tuple(src.shape), "float32", job_type, "distances")
143                if return_distances else None)
144    idx_out = (_allocate_output(indices, (ndim,) + tuple(src.shape), "int32", job_type, "indices")
145               if return_indices else None)
146    outputs: List[SourceLike] = []
147    if return_distances:
148        outputs.append(dist_out)
149    if return_indices:
150        outputs.append(idx_out)
151
152    if direct:
153        ret = bic.distance.distance_transform(
154            src[full_roi(ndim)], sampling=sampling, return_distances=return_distances,
155            return_indices=return_indices, number_of_threads=1,
156        )
157        if return_distances and return_indices:
158            dist, ind = ret
159        elif return_distances:
160            dist, ind = ret, None
161        else:
162            dist, ind = None, ret
163        if return_distances:
164            as_source(dist_out)[full_roi(ndim)] = dist
165        if return_indices:
166            as_source(idx_out)[full_roi(ndim + 1)] = ind
167    else:
168        runner = get_runner(job_type, job_config)
169        runner.run(_make_distance_block(sampling, return_distances, return_indices, ndim),
170                   [input], outputs=outputs, halo=normalize_halo(halo, ndim), block_shape=block_shape,
171                   num_workers=num_workers, block_ids=block_ids, resume_from=resume_from,
172                   name="distance_transform")
173
174    if return_distances and return_indices:
175        return dist_out, idx_out
176    return dist_out if return_distances else idx_out

218def map_points_to_objects(
219    segmentation: SourceLike,
220    points: np.ndarray,
221    block_shape: Tuple[int, ...],
222    *,
223    halo: Optional[Sequence[int]] = None,
224    sampling: Sampling = None,
225    num_workers: int = 1,
226    job_type: str = "local",
227    job_config: Optional[RunnerConfig] = None,
228    block_ids: Optional[Sequence[int]] = None,
229    resume_from: Optional[str] = None,
230) -> Tuple[np.ndarray, np.ndarray]:
231    """Map point coordinates to the nearest object in a segmentation and measure the distance.
232
233    Each (halo-padded) block computes the index transform of its background and maps the points it
234    contains to the nearest object. Choose ``halo`` to cover the maximum distance of interest; points
235    near a block boundary are resolved to the assignment with the smallest distance across overlapping
236    blocks.
237
238    Args:
239        segmentation: The label image (a numpy/zarr/n5 array or a `Source`).
240        points: The integer point coordinates, an ``(n_points, ndim)`` array.
241        block_shape: Shape of the processing blocks.
242        halo: Per-axis halo enlarging each block; choose it large enough to cover the maximum
243            distance of interest. ``None`` uses non-overlapping blocks (distances are then only
244            correct within each block).
245        sampling: The per-axis voxel spacing passed to the distance transform.
246        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed
247            backends).
248        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
249        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
250        block_ids: Restrict processing to these block ids (e.g. to re-run previously failed blocks).
251            Mutually exclusive with ``resume_from``.
252        resume_from: Distributed only; the preserved temp folder of a failed run to resume (see
253            ``runner.run``). Mutually exclusive with ``block_ids``.
254
255    Returns:
256        A ``(object_ids, object_distances)`` tuple, each of length ``n_points``: the id of the nearest
257        object per point (``0`` if none was found) and the corresponding distance (``inf`` if none).
258    """
259    check_rerun_args(job_type, resume_from, block_ids)
260    seg_src = as_source(segmentation)
261    points = np.asarray(points)
262    if points.ndim != 2 or points.shape[1] != seg_src.ndim:
263        raise ValueError(
264            f"points must be an (n_points, {seg_src.ndim}) array, got shape {points.shape}."
265        )
266    n_points = len(points)
267
268    has_halo = halo is not None
269    runner = get_runner(job_type, job_config)
270    results = runner.run(_make_map_points(points, sampling, has_halo, np.dtype(seg_src.dtype)),
271                         [segmentation], block_shape=block_shape,
272                         halo=normalize_halo(halo, seg_src.ndim) if has_halo else None,
273                         num_workers=num_workers, block_ids=block_ids, resume_from=resume_from,
274                         has_return_val=True, name="map_points_to_objects")
275
276    object_ids = np.zeros(n_points, dtype=seg_src.dtype)
277    object_distances = np.full(n_points, np.inf, dtype="float32")
278    for res in results:
279        if res is None:
280            continue
281        pids, oids, dists = res
282        take = dists < object_distances[pids]  # overlapping blocks: keep the closest assignment.
283        object_ids[pids[take]] = oids[take]
284        object_distances[pids[take]] = dists[take]
285    return object_ids, object_distances

48def find_local_maxima(
49    input: SourceLike,
50    *,
51    min_distance: int = 1,
52    threshold_abs: Optional[float] = None,
53    threshold_rel: Optional[float] = None,
54    block_shape: Optional[Tuple[int, ...]] = None,
55    job_type: str = "local",
56    job_config: Optional[RunnerConfig] = None,
57    num_workers: int = 1,
58    block_ids: Optional[Sequence[int]] = None,
59    resume_from: Optional[str] = None,
60) -> np.ndarray:
61    """Find local maxima of the data, block-wise (based on ``skimage.feature.peak_local_max``).
62
63    Args:
64        input: The input data (a numpy/zarr/n5 array or a `Source`).
65        min_distance: The minimum allowed distance between maxima (the non-maximum suppression
66            radius); also drives the block halo.
67        threshold_abs: Minimum intensity of a maximum. Defaults to the data minimum.
68        threshold_rel: Minimum intensity of a maximum, as a fraction of the data maximum.
69        block_shape: Shape of the processing blocks. Defaults to the input chunk shape; required
70            for unchunked data.
71        job_type: Execution backend: one of ``"local"``, ``"subprocess"`` or ``"slurm"``.
72        job_config: Backend configuration (a `RunnerConfig` / `SlurmConfig`).
73        num_workers: Number of parallel workers (threads for ``local``, tasks for distributed
74            backends).
75        block_ids: Restrict processing to these block ids; the maxima of just those blocks are
76            returned. Mutually exclusive with ``resume_from``.
77        resume_from: Distributed only; the preserved temp folder of a failed run to resume (see
78            ``runner.run``). Mutually exclusive with ``block_ids``.
79
80    Returns:
81        An ``(n_maxima, ndim)`` array of the detected maxima coordinates.
82    """
83    check_rerun_args(job_type, resume_from, block_ids)
84    src = as_source(input)
85    ndim = src.ndim
86    if check_direct(job_type, num_workers, block_shape, None, block_ids):
87        return peak_local_max(src[full_roi(ndim)], min_distance=min_distance,
88                              threshold_abs=threshold_abs, threshold_rel=threshold_rel)
89
90    halo = [min_distance + 8] * ndim
91    runner = get_runner(job_type, job_config)
92    results = runner.run(_make_local_maxima(min_distance, threshold_abs, threshold_rel),
93                         [input], block_shape=block_shape, halo=halo, num_workers=num_workers,
94                         block_ids=block_ids, resume_from=resume_from, has_return_val=True,
95                         name="find_local_maxima")
96    results = [r for r in results if r is not None]
97    if not results:
98        return np.zeros((0, ndim), dtype="int64")
99    return np.concatenate(results, axis=0)