implement convex region properties

These properties are computed based on the image_convex subimages:

- area_convex
- feret_diameter_max
- solidity

python/cucim/src/cucim/skimage/_shared/utils.py

@@ -27,6 +27,28 @@
 ]
 
 
+# For n nonzero elements cupy.nonzero returns a tuple of length ndim where
+# each element is an array of size (n, ) corresponding to the coordinates on
+# a specific axis.
+#
+# Often for regionprops purposes we would rather have a single array of
+# size (n, ndim) instead of a the tuple of arrays.
+#
+# CuPy's `_ndarray_argwhere` (used internally by cupy.nonzero) already provides
+# this but is not part of the public API. To guard against potential future
+# change we provide a less efficient fallback implementation.
+try:
+    from cupy._core._routines_indexing import _ndarray_argwhere
+except ImportError:
+
+    def _ndarray_argwhere(a):
+        """Stack the result of cupy.nonzero into a single array
+
+        output shape will be (num_nonzero, ndim)
+        """
+        return cp.stack(cp.nonzero(a), axis=-1)
+
+
 def _count_wrappers(func):
     """Count the number of wrappers around `func`."""
     unwrapped = func

python/cucim/src/cucim/skimage/measure/_regionprops_gpu.py

@@ -23,6 +23,11 @@
     regionprops_num_perimeter_pixels,
     regionprops_num_pixels,
 )
+from ._regionprops_gpu_convex import (
+    convex_deps,
+    regionprops_area_convex,
+    regionprops_feret_diameter_max,
+)
 from ._regionprops_gpu_intensity_kernels import (
     intensity_deps,
     regionprops_intensity_mean,
@@ -35,10 +40,12 @@
     "equivalent_diameter_area",
     "regionprops_area",
     "regionprops_area_bbox",
+    "regionprops_area_convex",
     "regionprops_bbox_coords",
     "regionprops_coords",
     "regionprops_dict",
     "regionprops_extent",
+    "regionprops_feret_diameter_max",
     "regionprops_image",
     "regionprops_intensity_mean",
     "regionprops_intensity_min_max",
@@ -101,6 +108,7 @@
 # get_property_dependencies below).
 property_deps = dict()
 property_deps.update(basic_deps)
+property_deps.update(convex_deps)
 property_deps.update(intensity_deps)
 
 # set of properties that only supports 2D images
@@ -376,6 +384,21 @@ def regionprops_dict(
             offset_coordinates=True,
         )
 
+    if "area_convex" in required_props:
+        regionprops_area_convex(
+            out["image_convex"], max_label=max_label, props_dict=out
+        )
+
+    if "solidity" in required_props:
+        out["solidity"] = out["area"] / out["area_convex"]
+
+    if "feret_diameter_max" in required_props:
+        regionprops_feret_diameter_max(
+            out["image_convex"],
+            spacing=spacing,
+            props_dict=out,
+        )
+
     compute_coords = "coords" in required_props
     compute_coords_scaled = "coords_scaled" in required_props
     if compute_coords or compute_coords_scaled:

python/cucim/src/cucim/skimage/measure/_regionprops_gpu_basic_kernels.py

@@ -49,7 +49,6 @@
 basic_deps["equivalent_spherical_perimeter"] = ["equivalent_diameter_area"]
 basic_deps["extent"] = ["area", "area_bbox"]
 basic_deps["image"] = []
-basic_deps["image_convex"] = []
 basic_deps["image_filled"] = ["label_filled"]
 basic_deps["image_intensity"] = []
 basic_deps["label"] = []

python/cucim/src/cucim/skimage/measure/_regionprops_gpu_convex.py

@@ -0,0 +1,147 @@
+import math
+import warnings
+from collections.abc import Sequence
+
+import cupy as cp
+
+from cucim.skimage._shared.distance import pdist_max_blockwise
+from cucim.skimage._shared.utils import _ndarray_argwhere
+from cucim.skimage._vendored import ndimage as ndi
+
+# Store information on which other properties a given property depends on
+# This information will be used by `regionprops_dict` to make sure that when
+# a particular property is requested any dependent properties are computed
+# first.
+convex_deps = dict()
+convex_deps["image_convex"] = ["image"]  # computed by regionprops_image
+convex_deps["area_convex"] = ["image_convex"]
+convex_deps["feret_diameter_max"] = ["image_convex"]
+convex_deps["solidity"] = ["area", "area_convex"]
+
+
+def regionprops_area_convex(
+    images_convex,
+    max_label=None,
+    spacing=None,
+    area_dtype=cp.float64,
+    props_dict=None,
+):
+    """Compute the area of each convex image.
+
+    writes "area_convex" to props_dict
+
+    Parameters
+    ----------
+    images_convex : sequence of cupy.ndarray
+        Convex images for each region as produced by ``regionprops_image`` with
+        ``compute_convex=True``.
+    """
+    if max_label is None:
+        max_label = len(images_convex)
+    if not isinstance(images_convex, Sequence):
+        raise ValueError("Expected `images_convex` to be a sequence of images")
+    area_convex = cp.zeros((max_label,), dtype=area_dtype)
+    for i in range(max_label):
+        area_convex[i] = images_convex[i].sum()
+    if spacing is not None:
+        if isinstance(spacing, cp.ndarray):
+            pixel_area = cp.product(spacing)
+        else:
+            pixel_area = math.prod(spacing)
+        area_convex *= pixel_area
+    if props_dict is not None:
+        props_dict["area_convex"] = area_convex
+    return area_convex
+
+
+def _regionprops_coords_perimeter(
+    image,
+    connectivity=1,
+):
+    """
+    Takes an image of a single labeled region (e.g. one element of the tuple
+    resulting from regionprops_image) and returns the coordinates of the voxels
+    at the edge of that region.
+    """
+
+    # remove non-boundary pixels
+    binary_image = image > 0
+    footprint = ndi.generate_binary_structure(
+        binary_image.ndim, connectivity=connectivity
+    )
+    binary_image_eroded = ndi.binary_erosion(binary_image, footprint)
+    binary_edges = binary_image * ~binary_image_eroded
+    edge_coords = _ndarray_argwhere(binary_edges)
+    return edge_coords
+
+
+def _feret_diameter_max(image_convex, spacing=None, return_argmax=False):
+    """Compute the maximum Feret diameter of a single convex image region."""
+    if image_convex.size == 1:
+        warnings.warn(
+            "single element image, returning 0 for feret diameter", UserWarning
+        )
+        return 0
+    coords = _regionprops_coords_perimeter(image_convex, connectivity=1)
+    coords = coords.astype(cp.float32)
+
+    if spacing is not None:
+        if all(s == 1.0 for s in spacing):
+            spacing = None
+        else:
+            spacing = cp.asarray(spacing, dtype=cp.float32).reshape(1, -1)
+            coords *= spacing
+
+    out = pdist_max_blockwise(
+        coords,
+        metric="sqeuclidean",
+        compute_argmax=return_argmax,
+        coords_per_block=4000,
+    )
+    if return_argmax:
+        return math.sqrt(out[0]), out[1]
+    return math.sqrt(out[0])
+
+
+def regionprops_feret_diameter_max(
+    images_convex, spacing=None, props_dict=None
+):
+    """Compute the maximum Feret diameter of the convex hull of each image in
+    images_convex.
+
+    writes "feret_diameter_max" to props_dict
+
+    Parameters
+    ----------
+    images_convex : sequence of cupy.ndarray
+        Convex images for each region as produced by ``regionprops_image`` with
+        ``compute_convex=True``.
+    spacing : tuple of float, optional
+        The pixel spacing of the image.
+    props_dict : dict, optional
+        A dictionary to store the computed properties.
+
+    Notes
+    -----
+    The maximum Feret diameter is the maximum distance between any two
+    points on the convex hull of the region. The implementation here is based
+    on pairwise distances of all boundary coordinates rather than using
+    marching squares or marching cubes as in scikit-image. The implementation
+    here is n-dimensional.
+
+    The distance is between pixel centers and so may be approximately one pixel
+    width less than the one computed by scikit-image.
+    """
+    if not isinstance(images_convex, Sequence):
+        raise ValueError("Expected `images_convex` to be a sequence of images")
+    diameters = cp.asarray(
+        tuple(
+            _feret_diameter_max(
+                image_convex, spacing=spacing, return_argmax=False
+            )
+            for image_convex in images_convex
+        )
+    )
+    if props_dict is not None:
+        props_dict["feret_diameter_max"] = diameters
+    return diameters