feat(examples): add segmentation & feature extraction examples

.github/workflows/notebooks.yml

@@ -0,0 +1,43 @@
+name: Convert notebooks to scripts
+
+on:
+  push:
+    branches: [main]
+    paths:
+      - "examples/notebooks/*.ipynb"
+
+jobs:
+  convert:
+    name: nbconvert to scripts
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: install uv
+        uses: astral-sh/setup-uv@v5
+
+      - name: set up python
+        run: uv python install 3.11
+
+      - name: install dependencies
+        run: uv sync --extra examples
+
+      - name: convert notebooks to scripts
+        run: |
+          for nb in examples/notebooks/*.ipynb; do
+            uv run jupyter nbconvert --to script "$nb" --output-dir examples/scripts
+          done
+
+      - name: commit generated scripts
+        run: |
+          git config user.name "github-actions[bot]"
+          git config user.email "github-actions[bot]@users.noreply.github.com"
+          git add examples/scripts/
+          if git diff --cached --quiet; then
+            echo "No changes to commit"
+          else
+            git commit -m "chore: auto-generate scripts from notebooks"
+            git push
+          fi

README.md

@@ -80,6 +80,8 @@ pull request on GitHub.
 | [Resolution Estimation (2D)](examples/notebooks/resolution_estimation_2d.ipynb) | FRC and DCR on STED microscopy data |
 | [Resolution Estimation (3D)](examples/notebooks/resolution_estimation_3d.ipynb) | FSC and DCR on 3D confocal pollen data |
 | [Deconvolution Iterations (3D)](examples/notebooks/deconvolution_iterations_3d.ipynb) | RL deconvolution stopping criteria via PSNR, SSIM, FSC, DCR |
+| [3D Monolayer Segmentation](examples/notebooks/segmentation_3d_monolayer.ipynb) | 3D nuclei and cell segmentation of hiPSC monolayer |
+| [3D Feature Extraction](examples/notebooks/feature_extraction_3d.ipynb) | GPU-accelerated regionprops on 3D fluorescence data |
 
 ## Citation
 If you use `cubic` in your research, please cite it:

cubic/feature/voxel.py

@@ -16,6 +16,9 @@ def regionprops(
     spacing: list[float] | None = None,
 ) -> list:
     """Extract region-based morphological features."""
+    # cucim requires spacing as tuple for kernel memoization
+    if spacing is not None:
+        spacing = tuple(spacing)
     return measure.regionprops(label_image, intensity_image, spacing=spacing)
 
 
@@ -30,6 +33,9 @@ def regionprops_table(
         properties = list(set(properties + ["label"]))
     else:
         properties = []
+    # cucim requires spacing as tuple for kernel memoization
+    if spacing is not None:
+        spacing = tuple(spacing)
     return measure.regionprops_table(
         label_image, intensity_image, properties=properties, spacing=spacing
     )

cubic/segmentation/segment_utils.py

@@ -20,51 +20,77 @@ def downscale_and_filter(
     downscale_anti_aliasing: bool = True,
     filter_size: int = 3,
     filter_shape: str = "square",
+    *,
+    downscale_xy_only: bool = True,
+    filter_mode: str = "nearest",
 ) -> npt.ArrayLike:
-    """Subsample and filter image prior to segmentiation.
+    """Subsample and filter image prior to segmentation.
 
     Parameters
     ----------
     image : npt.ArrayLike
         Image to be downsampled and filtered.
     downscale_factor : float, optional
         Factor by which to downscale the image, by default 0.5.
+    downscale_order : int, optional
+        Interpolation order for downscaling, by default 3.
+    downscale_anti_aliasing : bool, optional
+        Whether to apply anti-aliasing during downscaling, by default True.
+    downscale_xy_only : bool, optional
+        If True (default), only downscale XY dimensions, preserving Z for
+        3D images. If False, downscale all dimensions uniformly.
     filter_size : int, optional
         Size of median filter kernel, by default 3.
+    filter_shape : str, optional
+        Shape of the filter kernel: ``"square"`` (cube in 3D) uses
+        ``scipy.ndimage.median_filter(size=filter_size)`` which supports
+        boundary modes; ``"circular"`` (ball in 3D) uses
+        ``skimage.filters.median`` with a shaped footprint.
+    filter_mode : str, optional
+        Boundary mode for the median filter, by default ``"nearest"``.
+        Only used when ``filter_shape="square"``. Common values:
+        ``"constant"`` (zero-padding), ``"nearest"``, ``"reflect"``.
 
     Returns
     -------
     npt.ArrayLike
         Filtered and downsampled image.
 
     """
-    # cuCIM does not yet support rank-based median filter
-    # https://github.com/rapidsai/cucim/blob/main/python/cucim/src/cucim/skimage/filters/_median.py#L124
-    assert filter_shape in [
-        "square",
-        "circular",
-    ], "Filter shape must be 'square' or 'circular'."
+    from ..scipy import ndimage as _ndimage
+
+    if filter_shape not in ("square", "circular"):
+        raise ValueError("filter_shape must be 'square' or 'circular'.")
+
+    if downscale_factor < 1.0:
+        if downscale_xy_only:
+            from ..image_utils import rescale_xy
+
+            image = rescale_xy(
+                image,
+                scale=downscale_factor,
+                order=downscale_order,
+                anti_aliasing=downscale_anti_aliasing,
+            )
+        else:
+            image = transform.rescale(
+                image,
+                downscale_factor,
+                order=downscale_order,
+                anti_aliasing=downscale_anti_aliasing,
+            )
+
+    if filter_shape == "square":
+        return _ndimage.median_filter(image, size=filter_size, mode=filter_mode)
 
     if image.ndim == 2:
-        skimage_footprint = (
-            morphology.square if filter_shape == "square" else morphology.disk
-        )
+        footprint = morphology.disk(filter_size)
     elif image.ndim == 3:
-        skimage_footprint = (
-            morphology.cube if filter_shape == "square" else morphology.ball
-        )
+        footprint = morphology.ball(filter_size)
     else:
         raise ValueError("Image must be 2D or 3D.")
 
-    if downscale_factor < 1.0:
-        image = transform.rescale(
-            image,
-            downscale_factor,
-            order=downscale_order,
-            anti_aliasing=downscale_anti_aliasing,
-        )
-
-    return filters.median(image, footprint=skimage_footprint(filter_size))
+    return filters.median(image, footprint=footprint)
 
 
 def check_labeled_binary(image):
@@ -120,7 +146,7 @@ def cleanup_segmentation(
             )
             label_img[filled_mask] = label_id
 
-    return label(label_img).astype(np.uint8)
+    return label(label_img).astype(np.uint16)
 
 
 def find_objects(label_image, max_label=None):
@@ -253,26 +279,84 @@ def remove_thin_objects(label_image, min_z=2):
     return label_image
 
 
-def segment_watershed(image, markers=None, ball_size=15):
-    """Segment image using watershed algorithm."""
-    device = get_device(image)
+def segment_watershed(
+    image: npt.ArrayLike,
+    markers: npt.ArrayLike | None = None,
+    ball_size: int = 15,
+    *,
+    mask: npt.ArrayLike | None = None,
+    dilate_seeds: bool = False,
+) -> npt.ArrayLike:
+    """Segment image using watershed algorithm.
 
-    distance = distance_transform_edt(image)
-    coords = feature.peak_local_max(
-        distance, footprint=morphology.ball(ball_size), labels=image
-    )
+    When ``markers`` is None, computes a distance-based watershed:
+    EDT of the binary image is used to find peaks, which become markers,
+    and the watershed floods the negated distance.
+
+    When ``markers`` is provided, runs a marker-based watershed. By default
+    the image is used as both the landscape and mask. If ``mask`` is also
+    provided, the watershed uses the negated EDT of the mask as the
+    landscape (shape-based partitioning) and restricts flooding to the mask.
 
-    # https://github.com/rapidsai/cucim/issues/89
+    Parameters
+    ----------
+    image : npt.ArrayLike
+        Binary image to segment (distance-based) or intensity image
+        (marker-based when no mask is given).
+    markers : npt.ArrayLike or None, optional
+        Pre-computed markers for marker-based watershed. If None,
+        markers are generated from distance-transform peaks.
+    mask : npt.ArrayLike or None, optional
+        Binary mask restricting the watershed. When provided with markers,
+        the watershed landscape is the negated EDT of the mask (shape-based
+        partitioning). Only used when ``markers`` is not None.
+    ball_size : int, optional
+        Radius of the ball footprint for ``peak_local_max``, by default 15.
+        Only used when ``markers`` is None.
+    dilate_seeds : bool, optional
+        If True, dilate seed points with ``ball(1)`` before labeling.
+        This merges nearby peaks and reduces over-segmentation.
+        Only used when ``markers`` is None.
+
+    Returns
+    -------
+    npt.ArrayLike
+        Label image on the same device as the input.
+
+    """
+    from ..cuda import to_same_device
+
+    device = get_device(image)
+
+    # Distance-based watershed (no markers provided)
     if markers is None:
-        mask = np.zeros(distance.shape, dtype=bool)
-        mask[tuple(asnumpy(coords.T))] = True
-        markers = label(mask)
-        labels = watershed(-asnumpy(distance), markers, mask=asnumpy(image))
-    else:
-        labels = watershed(
-            asnumpy(image), markers=asnumpy(markers), mask=asnumpy(image)
-        )
-    # return on the same device as input
+        distance = distance_transform_edt(image)
+        footprint = morphology.ball(ball_size)
+        footprint = to_same_device(footprint, distance)
+        coords = feature.peak_local_max(distance, footprint=footprint, labels=image)
+
+        seed_mask = np.zeros(distance.shape, dtype=bool)
+        seed_mask[tuple(asnumpy(coords).T)] = True
+        seed_mask = to_device(seed_mask, device)
+        if dilate_seeds:
+            seed_mask = morphology.binary_dilation(
+                seed_mask, to_same_device(morphology.ball(1), seed_mask)
+            )
+        markers = label(seed_mask)
+        # watershed is not in cucim — run on CPU, return to original device
+        labels = watershed(-asnumpy(distance), asnumpy(markers), mask=asnumpy(image))
+        return to_device(labels, device)
+
+    # Marker-based watershed with explicit mask (shape-based partitioning)
+    if mask is not None:
+        distance = distance_transform_edt(asnumpy(mask))
+        ws_image = -distance
+        ws_image = ws_image - ws_image.min()
+        labels = watershed(ws_image, markers=asnumpy(markers), mask=asnumpy(mask))
+        return to_device(labels, device)
+
+    # Marker-based watershed without mask (image as landscape and mask)
+    labels = watershed(asnumpy(image), markers=asnumpy(markers), mask=asnumpy(image))
     return to_device(labels, device)
 
 

examples/data/README.md

@@ -32,3 +32,17 @@ Used in `examples/notebooks/resolution_estimation_3d.ipynb`.
 - **`40x_TAGoff_z_galvo.nd2`** — Pollen confocal (512x512x181, voxel 78x78x250 nm). Nikon A1 confocal, 40x/1.2 water, 488 nm excitation, GaAsP detector.
   [download from figshare](https://ndownloader.figshare.com/files/15203144)
   Source: Koho et al. (2019) *Nat. Commun.* 10:3103, [figshare dataset](https://doi.org/10.6084/m9.figshare.8159165.v1).
+
+#### 4. 3D cell monolayer segmentation
+
+Used in `examples/notebooks/segmentation_3d_monolayer.ipynb`. Auto-downloaded on first run using pooch.
+
+- **`3d_monolayer_xy1_ch0.tif`** — Membrane channel (60x256x256, uint16).
+- **`3d_monolayer_xy1_ch1.tif`** — Mitochondria channel (60x256x256, uint16).
+- **`3d_monolayer_xy1_ch2.tif`** — DNA channel (60x256x256, uint16).
+- **`3d_monolayer_xy1_ch2_NucleiLabels.tiff`** — CellProfiler nuclei reference labels.
+- **`3d_monolayer_xy1_ch0_CellsLabels.tiff`** — CellProfiler cell reference labels.
+
+Source: hiPSC data from the Allen Institute for Cell Science, provided with the
+[CellProfiler 3D monolayer tutorial](https://github.com/CellProfiler/tutorials/tree/master/3d_monolayer).
+[GitHub release assets](https://github.com/alxndrkalinin/cubic/releases/tag/v0.7.0a1).