Implement upstream changes from scikit-image 0.23 (part 2 of 2: morph…

…ology) (#728)

There were a number of non-trivial changes to the morphology module, so I broke those out from the other changes in #727. Please review and merge that MR first before reviewing this one.

Highlights from upstream are:
- binary morphology functions have a new `mode` argument that controls how values outside the image boundaries are interpreted
- grayscale morphology functions have new `mode` and `cval` arguments that control how boundaries are extended (these were already available in `scipy.ndimage`/`cupyx.scipy.ndimage`, they just weren't exposed via the `skimage`/`cuCIM` APIs)
- binary and grayscale morphology functions have bug fixes in the case of even-sized/non-symmetric footprints
  - additional corresponding test cases were added

Aside from the upstream changes, novel changes in this MR are:
- refactored utility functions to mirror and pad the footprints to allow use with the cuCIM-specific optimization of passing a tuple for rectangular footprints instead of explicitly allocating a GPU footprint array
- refactored some test cases to better use `pytest.mark.parametrize`
- some grayscale tests now compare directly to `skimage` CPU outputs instead fetching previously saved values
- bumped our version pinning for scikit-image to allow 0.23.x to be installed

I marked as "non-breaking" as the existing behavior has not changed except in the case of the bug fixes for even-sized footprints.

Authors:
  - Gregory Lee (https://github.com/grlee77)

Approvers:
  - Gigon Bae (https://github.com/gigony)
  - Jake Awe (https://github.com/AyodeAwe)

URL: #728

conda/environments/all_cuda-118_arch-x86_64.yaml

@@ -43,7 +43,7 @@ dependencies:
 - python>=3.8,<3.12
 - pywavelets>=1.0
 - recommonmark
-- scikit-image>=0.19.0,<0.23.0a0
+- scikit-image>=0.19.0,<0.24.0a0
 - scipy>=1.6.0
 - sphinx<6
 - sysroot_linux-64==2.17

conda/environments/all_cuda-122_arch-x86_64.yaml

@@ -42,7 +42,7 @@ dependencies:
 - python>=3.8,<3.12
 - pywavelets>=1.0
 - recommonmark
-- scikit-image>=0.19.0,<0.23.0a0
+- scikit-image>=0.19.0,<0.24.0a0
 - scipy>=1.6.0
 - sphinx<6
 - sysroot_linux-64==2.17

dependencies.yaml

@@ -247,7 +247,7 @@ dependencies:
           - click
           - lazy_loader>=0.1
           - numpy>=1.23.4,<2.0a0
-          - scikit-image>=0.19.0,<0.23.0a0
+          - scikit-image>=0.19.0,<0.24.0a0
           - scipy>=1.6.0
       - output_types: conda
         packages:

python/cucim/pyproject.toml

@@ -28,7 +28,7 @@ dependencies = [
     "cupy-cuda11x>=12.0.0",
     "lazy_loader>=0.1",
     "numpy>=1.23.4,<2.0a0",
-    "scikit-image>=0.19.0,<0.23.0a0",
+    "scikit-image>=0.19.0,<0.24.0a0",
     "scipy>=1.6.0",
 ] # This list was generated by `rapids-dependency-file-generator`. To make changes, edit ../../dependencies.yaml and run `rapids-dependency-file-generator`.
 classifiers = [

python/cucim/src/cucim/skimage/morphology/_skeletonize.py

@@ -56,7 +56,9 @@ def thin(image, max_num_iter=None):
     Parameters
     ----------
     image : binary (M, N) ndarray
-        The image to be thinned.
+        The image to thin. If this input isn't already a binary image,
+        it gets converted into one: In this case, zero values are considered
+        background (False), nonzero values are considered foreground (True).
     max_num_iter : int, number of iterations, optional
         Regardless of the value of this parameter, the thinned image
         is returned immediately if an iteration produces no change.
@@ -94,10 +96,10 @@ def thin(image, max_num_iter=None):
 
     Examples
     --------
-    >>> square = np.zeros((7, 7), dtype=np.uint8)
+    >>> square = np.zeros((7, 7), dtype=bool)
     >>> square[1:-1, 2:-2] = 1
     >>> square[0, 1] =  1
-    >>> square
+    >>> square.view(cp.uint8)
     array([[0, 1, 0, 0, 0, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
@@ -106,7 +108,7 @@ def thin(image, max_num_iter=None):
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
     >>> skel = thin(square)
-    >>> skel.astype(np.uint8)
+    >>> skel.view(np.uint8)
     array([[0, 1, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, 1, 0, 0, 0],
@@ -119,11 +121,11 @@ def thin(image, max_num_iter=None):
     check_nD(image, 2)
 
     # convert image to uint8 with values in {0, 1}
-    skel = cp.asarray(image, dtype=bool).astype(cp.uint8)
+    skel = cp.asarray(image, dtype=bool).view(cp.uint8)
 
     # neighborhood mask
     mask = cp.asarray(
-        [[8, 4, 2], [16, 0, 1], [32, 64, 128]], dtype=cp.uint8  # noqa  # noqa
+        [[8, 4, 2], [16, 0, 1], [32, 64, 128]], dtype=cp.uint8  # noqa
     )
 
     G123_LUT = cp.asarray(_G123_LUT)
@@ -177,7 +179,10 @@ def medial_axis(
     Parameters
     ----------
     image : binary ndarray, shape (M, N)
-        The image of the shape to be skeletonized.
+        The image of the shape to skeletonize. If this input isn't already a
+        binary image, it gets converted into one: In this case, zero values are
+        considered background (False), nonzero values are considered
+        foreground (True).
     mask : binary ndarray, shape (M, N), optional
         If a mask is given, only those elements in `image` with a true
         value in `mask` are used for computing the medial axis.
@@ -205,7 +210,7 @@ def medial_axis(
 
     See Also
     --------
-    skeletonize
+    skeletonize, thin
 
     Notes
     -----
@@ -232,17 +237,17 @@ def medial_axis(
 
     Examples
     --------
-    >>> square = np.zeros((7, 7), dtype=np.uint8)
+    >>> square = np.zeros((7, 7), dtype=bool)
     >>> square[1:-1, 2:-2] = 1
-    >>> square
+    >>> square.view(cp.uint8)
     array([[0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 1, 1, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
-    >>> medial_axis(square).astype(np.uint8)
+    >>> medial_axis(square).view(cp.uint8)
     array([[0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 1, 0, 0],
            [0, 0, 0, 1, 0, 0, 0],

python/cucim/src/cucim/skimage/morphology/binary.py

@@ -1,17 +1,16 @@
 """
 Binary morphological operations
 """
-import functools
 
 import cupy as cp
 
 import cucim.skimage._vendored.ndimage as ndi
 
-from .footprints import _footprint_is_sequence
+from .footprints import _footprint_is_sequence, pad_footprint
 from .misc import default_footprint
 
 
-def _iterate_binary_func(binary_func, image, footprint, out):
+def _iterate_binary_func(binary_func, image, footprint, out, border_value):
     """Helper to call `binary_func` for each footprint in a sequence.
 
     binary_func is a binary morphology function that accepts "structure",
@@ -24,7 +23,12 @@ def _iterate_binary_func(binary_func, image, footprint, out):
     #     `iterations > 1` is added.
     fp, num_iter = footprint[0]
     binary_func(
-        image, structure=fp, output=out, iterations=num_iter, brute_force=True
+        image,
+        structure=fp,
+        output=out,
+        iterations=num_iter,
+        border_value=border_value,
+        brute_force=True,
     )
     for fp, num_iter in footprint[1:]:
         # Note: out.copy() because the computation cannot be in-place!
@@ -34,6 +38,7 @@ def _iterate_binary_func(binary_func, image, footprint, out):
             structure=fp,
             output=out,
             iterations=num_iter,
+            border_value=border_value,
             brute_force=True,
         )
     return out
@@ -43,7 +48,7 @@ def _iterate_binary_func(binary_func, image, footprint, out):
 # default with the same dimension as the input image and size 3 along each
 # axis.
 @default_footprint
-def binary_erosion(image, footprint=None, out=None):
+def binary_erosion(image, footprint=None, out=None, *, mode="ignore"):
     """Return fast binary morphological erosion of an image.
 
     This function returns the same result as grayscale erosion but performs
@@ -65,6 +70,15 @@ def binary_erosion(image, footprint=None, out=None):
     out : ndarray of bool, optional
         The array to store the result of the morphology. If None is
         passed, a new array will be allocated.
+    mode : str, optional
+        The `mode` parameter determines how the array borders are handled.
+        Valid modes are: 'max', 'min', 'ignore'.
+        If 'max' or 'ignore', pixels outside the image domain are assumed
+        to be `True`, which causes them to not influence the result.
+        Default is 'ignore'.
+
+        .. versionadded:: 24.06
+            `mode` was added in 24.06.
 
     Returns
     -------
@@ -81,25 +95,37 @@ def binary_erosion(image, footprint=None, out=None):
     would apply a 9x1 footprint followed by a 1x9 footprint resulting in a net
     effect that is the same as ``footprint=cp.ones((9, 9))``, but with lower
     computational cost. Most of the builtin footprints such as
-    ``skimage.morphology.disk`` provide an option to automatically generate a
-    footprint sequence of this type.
+    :func:`skimage.morphology.disk` provide an option to automatically generate
+    a footprint sequence of this type.
 
+    For even-sized footprints, :func:`skimage.morphology.erosion` and
+    this function produce an output that differs: one is shifted by one pixel
+    compared to the other.
     """
     if out is None:
         out = cp.empty(image.shape, dtype=bool)
 
-    if _footprint_is_sequence(footprint):
-        binary_func = functools.partial(ndi.binary_erosion, border_value=True)
-        return _iterate_binary_func(binary_func, image, footprint, out)
+    if mode not in {"max", "min", "ignore"}:
+        raise ValueError(f"unsupported mode, got {mode!r}")
+    border_value = False if mode == "min" else True
+
+    footprint = pad_footprint(footprint, pad_end=True)
+
+    if not _footprint_is_sequence(footprint):
+        footprint = [(footprint, 1)]
 
-    ndi.binary_erosion(
-        image, structure=footprint, output=out, border_value=True
+    out = _iterate_binary_func(
+        binary_func=ndi.binary_erosion,
+        image=image,
+        footprint=footprint,
+        out=out,
+        border_value=border_value,
     )
     return out
 
 
 @default_footprint
-def binary_dilation(image, footprint=None, out=None):
+def binary_dilation(image, footprint=None, out=None, *, mode="ignore"):
     """Return fast binary morphological dilation of an image.
 
     This function returns the same result as grayscale dilation but performs
@@ -121,6 +147,15 @@ def binary_dilation(image, footprint=None, out=None):
     out : ndarray of bool, optional
         The array to store the result of the morphology. If None is
         passed, a new array will be allocated.
+    mode : str, optional
+        The `mode` parameter determines how the array borders are handled.
+        Valid modes are: 'max', 'min', 'ignore'.
+        If 'min' or 'ignore', pixels outside the image domain are assumed
+        to be `False`, which causes them to not influence the result.
+        Default is 'ignore'.
+
+        .. versionadded:: 24.06
+            `mode` was added in 24.06.
 
     Returns
     -------
@@ -137,22 +172,37 @@ def binary_dilation(image, footprint=None, out=None):
     would apply a 9x1 footprint followed by a 1x9 footprint resulting in a net
     effect that is the same as ``footprint=cp.ones((9, 9))``, but with lower
     computational cost. Most of the builtin footprints such as
-    ``skimage.morphology.disk`` provide an option to automatically generate a
-    footprint sequence of this type.
+    :func:`skimage.morphology.disk` provide an option to automatically generate
+    a footprint sequence of this type.
 
+    For non-symmetric footprints, :func:`skimage.morphology.binary_dilation`
+    and :func:`skimage.morphology.dilation` produce an output that differs:
+    `binary_dilation` mirrors the footprint, whereas `dilation` does not.
     """
     if out is None:
         out = cp.empty(image.shape, dtype=bool)
 
-    if _footprint_is_sequence(footprint):
-        return _iterate_binary_func(ndi.binary_dilation, image, footprint, out)
+    if mode not in {"max", "min", "ignore"}:
+        raise ValueError(f"unsupported mode, got {mode!r}")
+    border_value = True if mode == "max" else False
+
+    footprint = pad_footprint(footprint, pad_end=True)
+
+    if not _footprint_is_sequence(footprint):
+        footprint = [(footprint, 1)]
 
-    ndi.binary_dilation(image, structure=footprint, output=out)
+    out = _iterate_binary_func(
+        binary_func=ndi.binary_dilation,
+        image=image,
+        footprint=footprint,
+        out=out,
+        border_value=border_value,
+    )
     return out
 
 
 @default_footprint
-def binary_opening(image, footprint=None, out=None):
+def binary_opening(image, footprint=None, out=None, *, mode="ignore"):
     """Return fast binary morphological opening of an image.
 
     This function returns the same result as grayscale opening but performs
@@ -175,6 +225,15 @@ def binary_opening(image, footprint=None, out=None):
     out : ndarray of bool, optional
         The array to store the result of the morphology. If None
         is passed, a new array will be allocated.
+    mode : str, optional
+        The `mode` parameter determines how the array borders are handled.
+        Valid modes are: 'max', 'min', 'ignore'.
+        If 'ignore', pixels outside the image domain are assumed to be `True`
+        for the erosion and `False` for the dilation, which causes them to not
+        influence the result. Default is 'ignore'.
+
+        .. versionadded:: 24.06
+            `mode` was added in 24.06
 
     Returns
     -------
@@ -190,17 +249,16 @@ def binary_opening(image, footprint=None, out=None):
     would apply a 9x1 footprint followed by a 1x9 footprint resulting in a net
     effect that is the same as ``footprint=cp.ones((9, 9))``, but with lower
     computational cost. Most of the builtin footprints such as
-    ``skimage.morphology.disk`` provide an option to automatically generate a
-    footprint sequence of this type.
-
+    :func:`skimage.morphology.disk` provide an option to automatically generate
+    a footprint sequence of this type.
     """
-    eroded = binary_erosion(image, footprint)
-    out = binary_dilation(eroded, footprint, out=out)
+    tmp = binary_erosion(image, footprint, mode=mode)
+    out = binary_dilation(tmp, footprint, out=out, mode=mode)
     return out
 
 
 @default_footprint
-def binary_closing(image, footprint=None, out=None):
+def binary_closing(image, footprint=None, out=None, *, mode="ignore"):
     """Return fast binary morphological closing of an image.
 
     This function returns the same result as grayscale closing but performs
@@ -223,6 +281,15 @@ def binary_closing(image, footprint=None, out=None):
     out : ndarray of bool, optional
         The array to store the result of the morphology. If None,
         is passed, a new array will be allocated.
+    mode : str, optional
+        The `mode` parameter determines how the array borders are handled.
+        Valid modes are: 'max', 'min', 'ignore'.
+        If 'ignore', pixels outside the image domain are assumed to be `True`
+        for the erosion and `False` for the dilation, which causes them to not
+        influence the result. Default is 'ignore'.
+
+        .. versionadded:: 24.06
+            `mode` was added in 24.06.
 
     Returns
     -------
@@ -238,10 +305,9 @@ def binary_closing(image, footprint=None, out=None):
     would apply a 9x1 footprint followed by a 1x9 footprint resulting in a net
     effect that is the same as ``footprint=cp.ones((9, 9))``, but with lower
     computational cost. Most of the builtin footprints such as
-    ``skimage.morphology.disk`` provide an option to automatically generate a
-    footprint sequence of this type.
-
+    :func:`skimage.morphology.disk` provide an option to automatically generate
+    a footprint sequence of this type.
     """
-    dilated = binary_dilation(image, footprint)
-    out = binary_erosion(dilated, footprint, out=out)
+    tmp = binary_dilation(image, footprint, mode=mode)
+    out = binary_erosion(tmp, footprint, out=out, mode=mode)
     return out