diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6564ba372..54bc258ac 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -201,12 +201,14 @@ repos: ^python/cucim/src/cucim/skimage/morphology/_medial_axis_lookup[.]py$| ^python/cucim/src/cucim/skimage/morphology/binary[.]py$| ^python/cucim/src/cucim/skimage/morphology/convex_hull[.]py$| + ^python/cucim/src/cucim/skimage/morphology/extrema[.]py$| ^python/cucim/src/cucim/skimage/morphology/footprints[.]py$| ^python/cucim/src/cucim/skimage/morphology/gray[.]py$| ^python/cucim/src/cucim/skimage/morphology/isotropic[.]py$| ^python/cucim/src/cucim/skimage/morphology/misc[.]py$| ^python/cucim/src/cucim/skimage/morphology/tests/test_binary[.]py$| ^python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull[.]py$| + ^python/cucim/src/cucim/skimage/morphology/tests/test_extrema[.]py$| ^python/cucim/src/cucim/skimage/morphology/tests/test_footprints[.]py$| ^python/cucim/src/cucim/skimage/morphology/tests/test_gray[.]py$| ^python/cucim/src/cucim/skimage/morphology/tests/test_isotropic[.]py$| @@ -483,12 +485,14 @@ repos: python/cucim/src/cucim/skimage/morphology/_medial_axis_lookup[.]py$| python/cucim/src/cucim/skimage/morphology/binary[.]py$| python/cucim/src/cucim/skimage/morphology/convex_hull[.]py$| + python/cucim/src/cucim/skimage/morphology/extrema[.]py$| python/cucim/src/cucim/skimage/morphology/footprints[.]py$| python/cucim/src/cucim/skimage/morphology/gray[.]py$| python/cucim/src/cucim/skimage/morphology/isotropic[.]py$| python/cucim/src/cucim/skimage/morphology/misc[.]py$| python/cucim/src/cucim/skimage/morphology/tests/test_binary[.]py$| python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull[.]py$| + python/cucim/src/cucim/skimage/morphology/tests/test_extrema[.]py$| python/cucim/src/cucim/skimage/morphology/tests/test_footprints[.]py$| python/cucim/src/cucim/skimage/morphology/tests/test_gray[.]py$| python/cucim/src/cucim/skimage/morphology/tests/test_isotropic[.]py$| diff --git a/benchmarks/skimage/cucim_morphology_bench.py b/benchmarks/skimage/cucim_morphology_bench.py index 58664beaf..92400c116 100644 --- a/benchmarks/skimage/cucim_morphology_bench.py +++ b/benchmarks/skimage/cucim_morphology_bench.py @@ -16,10 +16,15 @@ import skimage.data import skimage.morphology from _image_bench import ImageBench +from packaging.version import Version import cucim.skimage import cucim.skimage.morphology +# scikit-image >= 0.26 renamed min_size -> max_size in remove_small_objects +# and area_threshold -> max_size in remove_small_holes +_skimage_has_max_size = Version(skimage.__version__) >= Version("0.26") + class BinaryMorphologyBench(ImageBench): def __init__( @@ -107,15 +112,26 @@ def _init_test_data(self, dtype): def set_args(self, dtype): a = self._init_test_data(dtype) - self.args_cpu = (cp.asnumpy(a), 6) - self.args_gpu = (a, 6) + self.args_cpu = (cp.asnumpy(a),) + self.args_gpu = (a,) + # scikit-image < 0.26 uses min_size, >= 0.26 uses max_size + if _skimage_has_max_size: + self.fixed_kwargs_cpu["max_size"] = 6 + else: + self.fixed_kwargs_cpu["min_size"] = 6 + self.fixed_kwargs_gpu["max_size"] = 6 class RemoveSmallHolesBench(RemoveSmallObjectsBench): def set_args(self, dtype): a = ~self._init_test_data(dtype) - self.args_cpu = (cp.asnumpy(a), 5) - self.args_gpu = (a, 5) + self.args_cpu = (cp.asnumpy(a),) + self.args_gpu = (a,) + if _skimage_has_max_size: + self.fixed_kwargs_cpu["max_size"] = 5 + else: + self.fixed_kwargs_cpu["area_threshold"] = 5 + self.fixed_kwargs_gpu["max_size"] = 5 def main(args): @@ -161,6 +177,37 @@ def main(args): ("thin", dict(), dict(), False, True), # grayreconstruct.py ("reconstruction", dict(), dict(), False, True), + # extrema.py + ( + "local_maxima", + dict(), + dict(connectivity=list(range(1, len(shape) + 1))), + False, + True, + ), + ( + "local_minima", + dict(), + dict(connectivity=list(range(1, len(shape) + 1))), + False, + True, + ), + # h values scaled to image range: integers use camera (0-255), + # floats use camera/255 (0-1) + ( + "h_maxima", + dict(), + dict(h=[10, 40] if np.issubdtype(dtypes[0], np.integer) else [0.04, 0.16]), + False, + True, + ), + ( + "h_minima", + dict(), + dict(h=[10, 40] if np.issubdtype(dtypes[0], np.integer) else [0.04, 0.16]), + False, + True, + ), # footprints.py # OMIT the functions from this file (each creates a structuring element) ]: @@ -298,7 +345,8 @@ def main(args): 'binary_closing', 'isotropic_erosion', 'isotropic_dilation', 'isotropic_opening', 'isotropic_closing','remove_small_objects', 'remove_small_holes', 'erosion', 'dilation', 'opening', 'closing', - 'white_tophat', 'black_tophat', 'thin', 'medial_axis', 'reconstruction' + 'white_tophat', 'black_tophat', 'thin', 'medial_axis', 'reconstruction', + 'local_maxima', 'local_minima', 'h_maxima', 'h_minima', ] # fmt: on dtype_choices = [ diff --git a/benchmarks/skimage/run-nv-bench-morphology.sh b/benchmarks/skimage/run-nv-bench-morphology.sh index 324c59749..f88339094 100755 --- a/benchmarks/skimage/run-nv-bench-morphology.sh +++ b/benchmarks/skimage/run-nv-bench-morphology.sh @@ -6,7 +6,7 @@ MAX_DURATION="${CUCIM_BENCHMARK_MAX_DURATION:-10}" param_shape=("512,512" "3840,2160" "3840,2160,3" "192,192,192") -param_filt=(binary_erosion binary_dilation binary_opening binary_closing isotropic_erosion isotropic_dilation isotropic_opening isotropic_closing remove_small_objects remove_small_holes erosion dilation opening closing white_tophat black_tophat medial_axis thin reconstruction) +param_filt=(binary_erosion binary_dilation binary_opening binary_closing isotropic_erosion isotropic_dilation isotropic_opening isotropic_closing remove_small_objects remove_small_holes erosion dilation opening closing white_tophat black_tophat medial_axis thin reconstruction h_maxima h_minima local_maxima local_minima) param_dt=(uint8) for shape in "${param_shape[@]}"; do for filt in "${param_filt[@]}"; do @@ -18,7 +18,7 @@ done # Note: Omit binary_*, medial_axis and thin from floating point benchmarks. # (these functions only take binary input). -param_filt_float=(remove_small_objects remove_small_holes erosion dilation opening closing white_tophat black_tophat reconstruction) +param_filt_float=(remove_small_objects remove_small_holes erosion dilation opening closing white_tophat black_tophat reconstruction h_maxima h_minima local_maxima local_minima) param_dt=(float32) for shape in "${param_shape[@]}"; do for filt in "${param_filt_float[@]}"; do diff --git a/python/cucim/src/cucim/skimage/morphology/__init__.py b/python/cucim/src/cucim/skimage/morphology/__init__.py index 42936cfe5..24de21c34 100644 --- a/python/cucim/src/cucim/skimage/morphology/__init__.py +++ b/python/cucim/src/cucim/skimage/morphology/__init__.py @@ -1,10 +1,11 @@ # SPDX-FileCopyrightText: 2009-2022 the scikit-image team -# SPDX-FileCopyrightText: Copyright (c) 2021-2025, NVIDIA CORPORATION. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2021-2026, NVIDIA CORPORATION. All rights reserved. # SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause """Morphological algorithms, e.g., closing, opening, skeletonization.""" from ._skeletonize import medial_axis, thin +from .extrema import h_maxima, h_minima, local_maxima, local_minima from .binary import ( binary_closing, binary_dilation, @@ -50,6 +51,10 @@ "binary_dilation", "binary_opening", "binary_closing", + "h_maxima", + "h_minima", + "local_maxima", + "local_minima", "isotropic_dilation", "isotropic_erosion", "isotropic_opening", diff --git a/python/cucim/src/cucim/skimage/morphology/extrema.py b/python/cucim/src/cucim/skimage/morphology/extrema.py new file mode 100644 index 000000000..db9673017 --- /dev/null +++ b/python/cucim/src/cucim/skimage/morphology/extrema.py @@ -0,0 +1,710 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +""" +Local extrema detection functions. + +These functions identify local maxima and minima in n-dimensional arrays, +with proper handling of plateaus (connected regions of equal value). +""" + +import warnings + +import cupy as cp +import numpy as np + +import cucim.skimage._vendored.ndimage as ndi +from cucim.skimage._shared.utils import warn +from cucim.skimage.morphology import grayreconstruct +from cucim.skimage.util import dtype_limits, invert + +__all__ = ["h_maxima", "h_minima", "local_maxima", "local_minima"] + + +def _add_constant_clip(image, const_value): + """Add constant to the image while handling overflow issues gracefully.""" + min_dtype, max_dtype = dtype_limits(image, clip_negative=False) + + if const_value > (max_dtype - min_dtype): + raise ValueError( + "The added constant is not compatible with the image data type." + ) + + result = image + image.dtype.type(const_value) + clip_mask = image > image.dtype.type(max_dtype - const_value) + result[clip_mask] = image.dtype.type(max_dtype) + return result + + +def _subtract_constant_clip(image, const_value): + """Subtract constant from image while handling underflow issues.""" + min_dtype, max_dtype = dtype_limits(image, clip_negative=False) + + if const_value > (max_dtype - min_dtype): + raise ValueError( + "The subtracted constant is not compatible with the image data type." + ) + + result = image - image.dtype.type(const_value) + clip_mask = image < image.dtype.type(const_value + min_dtype) + result[clip_mask] = image.dtype.type(min_dtype) + return result + + +def h_maxima(image, h, footprint=None, *, reconstruct_on_cpu="auto"): + """Determine all maxima of the image with height >= h. + + The local maxima are defined as connected sets of pixels with equal + gray level strictly greater than the gray level of all pixels in direct + neighborhood of the set. + + A local maximum M of height h is a local maximum for which + there is at least one path joining M with an equal or higher local maximum + on which the minimal value is f(M) - h (i.e. the values along the path + are not decreasing by more than h with respect to the maximum's value) + and no path to an equal or higher local maximum for which the minimal + value is greater. + + The global maxima of the image are also found by this function. + + Parameters + ---------- + image : ndarray + The input image for which the maxima are to be calculated. + h : unsigned integer + The minimal height of all extracted maxima. + footprint : ndarray, optional + The neighborhood expressed as an n-D array of 1's and 0's. + Default is the ball of radius 1 according to the maximum norm + (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.) + reconstruct_on_cpu : bool or {'auto'}, optional + If ``'auto'`` (default), the internal morphological reconstruction uses + the GPU path when supported by the footprint and offset. If False, it + always runs entirely on the GPU. If True, it falls back to + scikit-image's CPU-based reconstruction loop. See + `cucim.skimage.morphology.reconstruction` for details. + + Returns + ------- + h_max : ndarray + The local maxima of height >= h and the global maxima. + The resulting image is a binary image, where pixels belonging to + the determined maxima take value 1, the others take value 0. + + See Also + -------- + cucim.skimage.morphology.h_minima + cucim.skimage.morphology.local_maxima + cucim.skimage.morphology.local_minima + + References + ---------- + .. [1] Soille, P., "Morphological Image Analysis: Principles and + Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883. + + Examples + -------- + >>> import cupy as cp + >>> from cucim.skimage.morphology import extrema + + We create an image (quadratic function with a maximum in the center and + 4 additional constant maxima. + The heights of the maxima are: 1, 21, 41, 61, 81 + + >>> w = 10 + >>> x, y = cp.mgrid[0:w,0:w] + >>> f = 20 - 0.2*((x - w/2)**2 + (y-w/2)**2) + >>> f[2:4,2:4] = 40; f[2:4,7:9] = 60; f[7:9,2:4] = 80; f[7:9,7:9] = 100 + >>> f = f.astype(int) + + We can calculate all maxima with a height of at least 40: + + >>> maxima = extrema.h_maxima(f, 40) + + The resulting image will contain 3 local maxima. + """ + + # Check for h value that is larger than the range of the image. If this + # is True then there are no h-maxima in the image. + if h > cp.ptp(image): + return cp.zeros(image.shape, dtype=np.uint8) + + # Check for floating point h value. For this to work properly + # we need to explicitly convert image to float64. + # + # FIXME: This could give incorrect results if image is int64 and + # has a very high dynamic range. The dtype of image is + # changed to float64, and different integer values could + # become the same float due to rounding. + # + # see: https://github.com/rapidsai/cucim/issues/1090 + # + # >>> ii64 = np.iinfo(np.int64) + # >>> a = np.array([ii64.max, ii64.max - 2]) + # >>> a[0] == a[1] + # False + # >>> b = a.astype(np.float64) + # >>> b[0] == b[1] + # True + # + if np.issubdtype(type(h), np.floating) and np.issubdtype( + image.dtype, np.integer + ): + if (h % 1) != 0: + warn( + "possible precision loss converting image to " + "floating point. To silence this warning, " + "ensure image and h have same data type.", + stacklevel=2, + ) + image = image.astype(float) + else: + h = image.dtype.type(h) + + if h == 0: + raise ValueError("h = 0 is ambiguous, use local_maxima() instead?") + + if np.issubdtype(image.dtype, np.floating): + # The purpose of the resolution variable is to allow for the + # small rounding errors that inevitably occur when doing + # floating point arithmetic. We want shifted_img to be + # guaranteed to be h less than image. If we only subtract h + # there may be pixels where shifted_img ends up being + # slightly greater than image - h. + # + # The resolution is scaled based on the pixel values in the + # image because floating point precision is relative. A + # very large value of 1.0e10 will have a large precision, + # say +-1.0e4, and a very small value of 1.0e-10 will have + # a very small precision, say +-1.0e-16. + # + resolution = (2 * cp.finfo(image.dtype).resolution) * cp.abs(image) + shifted_img = image - h - resolution + else: + shifted_img = _subtract_constant_clip(image, h) + + rec_img = grayreconstruct.reconstruction( + shifted_img, + image, + method="dilation", + footprint=footprint, + reconstruct_on_cpu=reconstruct_on_cpu, + ) + residue_img = image - rec_img + return (residue_img >= h).astype(cp.uint8) + + +def h_minima(image, h, footprint=None, *, reconstruct_on_cpu="auto"): + """Determine all minima of the image with depth >= h. + + The local minima are defined as connected sets of pixels with equal + gray level strictly smaller than the gray levels of all pixels in direct + neighborhood of the set. + + A local minimum M of depth h is a local minimum for which + there is at least one path joining M with an equal or lower local minimum + on which the maximal value is f(M) + h (i.e. the values along the path + are not increasing by more than h with respect to the minimum's value) + and no path to an equal or lower local minimum for which the maximal + value is smaller. + + The global minima of the image are also found by this function. + + Parameters + ---------- + image : ndarray + The input image for which the minima are to be calculated. + h : unsigned integer + The minimal depth of all extracted minima. + footprint : ndarray, optional + The neighborhood expressed as an n-D array of 1's and 0's. + Default is the ball of radius 1 according to the maximum norm + (i.e. a 3x3 square for 2D images, a 3x3x3 cube for 3D images, etc.) + reconstruct_on_cpu : bool or {'auto'}, optional + If ``'auto'`` (default), the internal morphological reconstruction uses + the GPU path when supported by the footprint and offset. If False, it + always runs entirely on the GPU. If True, it falls back to + scikit-image's CPU-based reconstruction loop. See + `cucim.skimage.morphology.reconstruction` for details. + + See Also + -------- + skimage.morphology.h_maxima + skimage.morphology.local_maxima + skimage.morphology.local_minima + + References + ---------- + .. [1] Soille, P., "Morphological Image Analysis: Principles and + Applications" (Chapter 6), 2nd edition (2003), ISBN 3540429883. + + Examples + -------- + >>> import cupy as cp + >>> from cucim.skimage.morphology import extrema + + We create an image (quadratic function with a minimum in the center and + 4 additional constant maxima. + The depth of the minima are: 1, 21, 41, 61, 81 + + >>> w = 10 + >>> x, y = cp.mgrid[0:w,0:w] + >>> f = 180 + 0.2*((x - w/2)**2 + (y-w/2)**2) + >>> f[2:4,2:4] = 160; f[2:4,7:9] = 140; f[7:9,2:4] = 120; f[7:9,7:9] = 100 + >>> f = f.astype(int) + + We can calculate all minima with a depth of at least 40: + + >>> minima = extrema.h_minima(f, 40) + + The resulting image will contain 3 local minima. + """ + if h > cp.ptp(image): + return cp.zeros(image.shape, dtype=cp.uint8) + + if np.issubdtype(type(h), np.floating) and np.issubdtype( + image.dtype, np.integer + ): + if (h % 1) != 0: + warn( + "possible precision loss converting image to " + "floating point. To silence this warning, " + "ensure image and h have same data type.", + stacklevel=2, + ) + image = image.astype(float) + else: + h = image.dtype.type(h) + + if h == 0: + raise ValueError("h = 0 is ambiguous, use local_minima() instead?") + + if np.issubdtype(image.dtype, np.floating): + resolution = 2 * cp.finfo(image.dtype).resolution * cp.abs(image) + shifted_img = image + h + resolution + else: + shifted_img = _add_constant_clip(image, h) + + rec_img = grayreconstruct.reconstruction( + shifted_img, + image, + method="erosion", + footprint=footprint, + reconstruct_on_cpu=reconstruct_on_cpu, + ) + residue_img = rec_img - image + return (residue_img >= h).astype(cp.uint8) + + +def _resolve_footprint(footprint, connectivity, ndim): + """ + Resolve footprint from footprint or connectivity parameter. + + Parameters + ---------- + footprint : ndarray or None + Explicit footprint array. + connectivity : int or None + Connectivity value for generating footprint. + ndim : int + Number of dimensions. + + Returns + ------- + footprint : cupy.ndarray + Boolean footprint array with shape (3,) * ndim. + + Raises + ------ + ValueError + If footprint has wrong number of dimensions or wrong size. + """ + if footprint is None: + if connectivity is None: + # Default: full connectivity + connectivity = ndim + # Clamp connectivity to valid range + connectivity = min(connectivity, ndim) + # Generate binary structure: returns (3,) * ndim shaped array + footprint = ndi.generate_binary_structure(ndim, connectivity) + else: + # Convert footprint to cupy boolean array + footprint = cp.asarray(footprint, dtype=bool) + + # Validate footprint dimensions + if footprint.ndim != ndim: + raise ValueError( + f"footprint must have the same number of dimensions as image, " + f"got footprint.ndim={footprint.ndim} but image.ndim={ndim}" + ) + + # Validate footprint size (must be 3 in each dimension) + for i, s in enumerate(footprint.shape): + if s != 3: + raise ValueError( + f"Each dimension size of footprint must be 3, " + f"got footprint.shape[{i}]={s}" + ) + + return footprint + + +def local_maxima( + image, footprint=None, connectivity=None, indices=False, allow_borders=True +): + """ + Find local maxima of n-dimensional array. + + This function finds local maxima in an image, where a local maximum is + defined as a connected set of pixels (a plateau) with equal gray level + that is strictly greater than ALL pixels in its external neighborhood. + + Parameters + ---------- + image : ndarray + Input array to find local maxima in. + footprint : ndarray, optional + Boolean array defining the neighborhood. If None, a full connectivity + footprint is generated based on the ``connectivity`` parameter. + ``footprint`` and ``connectivity`` are mutually exclusive. + connectivity : int, optional + The connectivity defining the neighborhood of a pixel. Used to generate + a footprint if ``footprint`` is None. A connectivity of 1 means pixels + sharing an edge (4-connectivity in 2D, 6-connectivity in 3D). A + connectivity equal to ``image.ndim`` means pixels sharing a corner + (8-connectivity in 2D, 26-connectivity in 3D). + Default is ``image.ndim`` (full connectivity). + indices : bool, optional + If True, return a tuple of arrays containing the coordinates of local + maxima. If False (default), return a boolean mask with the same shape + as ``image``. + allow_borders : bool, optional + If True (default), local maxima at the image border are detected. + If False, pixels at the image border are excluded from consideration. + + Returns + ------- + maxima : ndarray or tuple of ndarrays + If ``indices=False``, returns a boolean array of the same shape as + ``image`` where True indicates pixels that are part of a local maximum. + If ``indices=True``, returns a tuple of 1-D arrays containing the + coordinates of local maxima (as returned by ``cp.nonzero``). + + See Also + -------- + local_minima : Find local minima. + skimage.feature.peak_local_max : Find peaks with additional filtering. + + Notes + ----- + This function identifies strict local maxima: a connected plateau is + only considered a local maximum if ALL pixels in its external neighborhood + have strictly lower values (not equal). When a plateau is identified as + a local maximum, ALL pixels belonging to that plateau are marked as True. + + Integer inputs use ``float64`` intermediates while validating plateaus. + For high-range ``int64`` or ``uint64`` inputs, distinct adjacent values may + round to the same ``float64`` value. In that case, extrema differing by + less than the local ``float64`` spacing may be missed. + + The algorithm works as follows: + 1. Find candidate pixels using a hollow (center-excluded) maximum filter + 2. Label connected regions of candidates + 3. Validate each labeled region has uniform intensity + 4. Check that each valid plateau's external boundary has strictly lower + values + + Examples + -------- + >>> import cupy as cp + >>> from cucim.skimage.morphology import local_maxima + >>> image = cp.array([[0, 0, 0, 0, 0], + ... [0, 1, 1, 1, 0], + ... [0, 1, 2, 1, 0], + ... [0, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0]], dtype=cp.float32) + >>> local_maxima(image) + array([[False, False, False, False, False], + [False, False, False, False, False], + [False, False, True, False, False], + [False, False, False, False, False], + [False, False, False, False, False]]) + + For a plateau maximum, all plateau pixels are marked True: + + >>> image = cp.array([[0, 0, 0], + ... [0, 5, 5], + ... [0, 5, 5]], dtype=cp.float32) + >>> local_maxima(image) + array([[False, False, False], + [False, True, True], + [False, True, True]]) + """ + image = cp.asarray(image) + + if image.size == 0: + if indices: + return tuple(cp.array([], dtype=cp.intp) for _ in range(image.ndim)) + return cp.zeros(image.shape, dtype=bool) + + # Validate mutually exclusive parameters + if footprint is not None and connectivity is not None: + raise ValueError( + "footprint and connectivity are mutually exclusive, " + "only provide one." + ) + + # Call the internal implementation + return _local_maxima( + image, + footprint=footprint, + connectivity=connectivity, + indices=indices, + allow_borders=allow_borders, + ) + + +def local_minima( + image, footprint=None, connectivity=None, indices=False, allow_borders=True +): + """ + Find local minima of n-dimensional array. + + This function finds local minima in an image, where a local minimum is + defined as a connected set of pixels (a plateau) with equal gray level + that is strictly less than ALL pixels in its external neighborhood. + + Parameters + ---------- + image : ndarray + Input array to find local minima in. + footprint : ndarray, optional + Boolean array defining the neighborhood. If None, a full connectivity + footprint is generated based on the ``connectivity`` parameter. + ``footprint`` and ``connectivity`` are mutually exclusive. + connectivity : int, optional + The connectivity defining the neighborhood of a pixel. Used to generate + a footprint if ``footprint`` is None. A connectivity of 1 means pixels + sharing an edge (4-connectivity in 2D, 6-connectivity in 3D). A + connectivity equal to ``image.ndim`` means pixels sharing a corner + (8-connectivity in 2D, 26-connectivity in 3D). + Default is ``image.ndim`` (full connectivity). + indices : bool, optional + If True, return a tuple of arrays containing the coordinates of local + minima. If False (default), return a boolean mask with the same shape + as ``image``. + allow_borders : bool, optional + If True (default), local minima at the image border are detected. + If False, pixels at the image border are excluded from consideration. + + Returns + ------- + minima : ndarray or tuple of ndarrays + If ``indices=False``, returns a boolean array of the same shape as + ``image`` where True indicates pixels that are part of a local minimum. + If ``indices=True``, returns a tuple of 1-D arrays containing the + coordinates of local minima (as returned by ``cp.nonzero``). + + See Also + -------- + local_maxima : Find local maxima. + + Notes + ----- + This function identifies strict local minima: a connected plateau is + only considered a local minimum if ALL pixels in its external neighborhood + have strictly higher values (not equal). When a plateau is identified as + a local minimum, ALL pixels belonging to that plateau are marked as True. + + Integer inputs use ``float64`` intermediates while validating plateaus. + For high-range ``int64`` or ``uint64`` inputs, distinct adjacent values may + round to the same ``float64`` value. In that case, extrema differing by + less than the local ``float64`` spacing may be missed. + + Examples + -------- + >>> import cupy as cp + >>> from cucim.skimage.morphology import local_minima + >>> image = cp.array([[5, 5, 5, 5, 5], + ... [5, 4, 4, 4, 5], + ... [5, 4, 3, 4, 5], + ... [5, 4, 4, 4, 5], + ... [5, 5, 5, 5, 5]], dtype=cp.float32) + >>> local_minima(image) + array([[False, False, False, False, False], + [False, False, False, False, False], + [False, False, True, False, False], + [False, False, False, False, False], + [False, False, False, False, False]]) + """ + return local_maxima( + invert(image, signed_float=True), + footprint=footprint, + connectivity=connectivity, + indices=indices, + allow_borders=allow_borders, + ) + + +def _local_maxima(image, footprint, connectivity, indices, allow_borders): + """ + Internal implementation of local_maxima. + + Parameters + ---------- + image : cupy.ndarray + Input image (already validated as cupy array). + footprint : cupy.ndarray or None + Footprint defining the neighborhood. + connectivity : int or None + Connectivity for generating footprint. + indices : bool + Whether to return indices instead of mask. + allow_borders : bool + Whether to allow maxima at borders. + + Returns + ------- + result : cupy.ndarray or tuple of cupy.ndarray + Boolean mask or tuple of coordinate arrays. + """ + # Validate dtype - float16 is not supported + if image.dtype == np.float16: + raise TypeError( + f"dtype {image.dtype} which is not supported, try using " + "float32 or float64 instead" + ) + + ndim = image.ndim + original_shape = image.shape + + # Resolve footprint from footprint or connectivity parameter + footprint = _resolve_footprint(footprint, connectivity, ndim) + + # Handle small arrays (any dimension < 3) + # When allow_borders=False, we cannot apply a 3x3 footprint properly + if not allow_borders and any(s < 3 for s in image.shape): + warnings.warn( + "maxima can't exist with any dimension smaller 3 " + "if borders are excluded", + UserWarning, + stacklevel=4, + ) + if indices: + return tuple(cp.array([], dtype=cp.intp) for _ in range(ndim)) + return cp.zeros(original_shape, dtype=bool) + + # Pad image if allow_borders=True + # This allows border pixels to be compared against "virtual" neighbors + # with the minimum value, making them potential maxima + if allow_borders: + # Get minimum value for padding + if np.issubdtype(image.dtype, np.floating): + pad_value = float(image.min()) + else: + pad_value = int(image.min()) + + # Pad with 1 pixel on each side + image = cp.pad( + image, pad_width=1, mode="constant", constant_values=pad_value + ) + + # Step 3: Find candidates using maximum filter + # Apply standard maximum filter (more efficient than hollow footprint + # since rectangular footprints use separable 1D filters internally) + # image_max[i,j] = max value in neighborhood of (i,j) including center + image_max = ndi.maximum_filter(image, footprint=footprint) + + # Candidates are pixels where value equals the local maximum + # i.e., pixels that are >= all neighbors in their neighborhood + # These could be single-pixel maxima or part of a plateau + candidates = image == image_max + + # Step 4: Label candidate regions and validate + # Connected candidates form potential plateaus + labels, num_labels = ndi.label(candidates, structure=footprint) + + if num_labels == 0: + # No candidates found + result = cp.zeros(image.shape, dtype=bool) + else: + # Step 5: Validate plateau boundaries using two checks: + # + # Check 1: No equal-value non-candidate neighbors + # If a candidate pixel has a neighbor with the same value that is NOT + # a candidate, that neighbor must see something higher. This means the + # plateau's extended neighborhood includes higher values → invalid. + # + # We compute max of non-candidate neighbor values. For candidates, + # we use a very low value so they don't contribute. + if np.issubdtype(image.dtype, np.floating): + low_val = float(image.min()) - 1.0 + else: + # For integers, use float to allow going below min + low_val = float(image.min()) - 1.0 + + non_candidate_image = cp.where( + candidates, low_val, image.astype(cp.float64) + ) + max_non_candidate_neighbor = ndi.maximum_filter( + non_candidate_image, footprint=footprint + ) + # For each candidate, check if any non-candidate neighbor has value + # >= candidate value. If so, this candidate (and its region) is invalid. + no_bad_neighbor = image.astype(cp.float64) > max_non_candidate_neighbor + + # Check 2: Region must have at least one strictly lower neighbor + # This ensures the region is a STRICT maximum (not a constant plateau). + # For constant images, no pixel has a strictly lower neighbor. + min_neighbor = ndi.minimum_filter(image, footprint=footprint) + has_lower_neighbor = image > min_neighbor + + # Find labels where ANY pixel fails check 1 (has bad neighbor) + # These labels are invalid + bad_labels_mask = candidates & ~no_bad_neighbor + if bad_labels_mask.any(): + bad_labels = cp.unique(labels[bad_labels_mask]) + bad_labels = bad_labels[bad_labels > 0] + else: + bad_labels = cp.array([], dtype=labels.dtype) + + # Find labels where AT LEAST ONE pixel passes check 2 (has lower) + # These labels could be valid (if they also pass check 1) + good_labels_mask = candidates & has_lower_neighbor + if good_labels_mask.any(): + labels_with_lower = cp.unique(labels[good_labels_mask]) + else: + labels_with_lower = cp.array([], dtype=labels.dtype) + + # Build validity lookup: + # Valid = has at least one lower neighbor AND no bad neighbors + valid_labels_lookup = cp.zeros(num_labels + 1, dtype=bool) + valid_labels_lookup[labels_with_lower] = True # start with these + valid_labels_lookup[bad_labels] = False # remove invalid + valid_labels_lookup[0] = False # background is never valid + + # When allow_borders=False, invalidate labels that touch the border + # (not just border pixels, but entire connected regions touching border) + if not allow_borders: + # Lazy import to avoid circular import with segmentation module + from cucim.skimage.segmentation._clear_border import ( + _get_border_labels, + ) + + # Find labels that touch the border and mark them invalid + border_labels = _get_border_labels(labels) + border_labels = border_labels[border_labels > 0] + valid_labels_lookup[border_labels] = False + + result = valid_labels_lookup[labels] + + # Remove padding from result if we added it + if allow_borders: + # Extract the center region (remove padding) + slices = tuple(slice(1, -1) for _ in range(ndim)) + result = result[slices] + + if indices: + return cp.nonzero(result) + return result diff --git a/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py b/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py index 5442caa38..cee9e8b8f 100644 --- a/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py +++ b/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py @@ -1,7 +1,7 @@ # SPDX-FileCopyrightText: Copyright (c) 2003-2009 Massachusetts Institute of Technology # SPDX-FileCopyrightText: Copyright (c) 2009-2011 Broad Institute # SPDX-FileCopyrightText: 2003 Lee Kamentsky -# SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION. All rights reserved. # SPDX-License-Identifier: Apache-2.0 AND (GPL-2.0-only OR BSD-3-Clause) """ @@ -21,10 +21,20 @@ import skimage from packaging.version import Version +import cucim.skimage._vendored.ndimage as ndi + old_reconstruction_pyx = Version(skimage.__version__) < Version("0.20.0") -def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): +def reconstruction( + seed, + mask, + method="dilation", + footprint=None, + offset=None, + *, + reconstruct_on_cpu="auto", +): """Perform a morphological reconstruction of an image. Morphological reconstruction by dilation is similar to basic morphological @@ -66,6 +76,12 @@ def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): The coordinates of the center of the footprint. Default is located on the geometrical center of the footprint, in that case footprint dimensions must be odd. + reconstruct_on_cpu : bool or {'auto'}, optional + If ``'auto'`` (default), use the iterative GPU-native algorithm when + the footprint and offset can be represented by the GPU filter path, + otherwise fall back to scikit-image's Cython reconstruction_loop on + the CPU. If False, always use the GPU-native algorithm. If True, + always use the CPU algorithm (requires host-device transfers). Returns ------- @@ -124,8 +140,34 @@ def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): Notes ----- - The algorithm is taken from [1]_. Applications for grayscale reconstruction - are discussed in [2]_ and [3]_. + Two algorithm variants are available, selected by the + ``reconstruct_on_cpu`` parameter: + + When ``reconstruct_on_cpu=True``, the CPU-based algorithm from [1]_ is used. This + processes pixels in rank order using a linked-list traversal via + scikit-image's Cython ``reconstruction_loop``. Data is transferred to + the host for the inner loop and back to the device afterward. + + When ``reconstruct_on_cpu=False``, an iterative parallel algorithm runs + entirely on the GPU. Each iteration applies a conditional dilation (or + erosion) via ``maximum_filter`` (or ``minimum_filter``) and clamps by the + mask, repeating until convergence. With the default + ``reconstruct_on_cpu='auto'``, this GPU path is used except in rare cases + not supported by the GPU algorithm, such as a non-centered offset for an + even-sized footprint that cannot be represented by ``ndimage``'s filter + origin. This is highly efficient when the seed is close to the mask (e.g., + in ``h_maxima`` / ``h_minima``), typically converging in few iterations + with each iteration fully parallelized. However, in pathological cases + where a single seed value must propagate across a long path (e.g., a + single-pixel seed with a uniform mask), convergence can require O(N) + iterations and the CPU path may be faster. + + An alternative GPU algorithm based on irregular wavefront propagation [4]_ + would handle such pathological cases efficiently by tiling the image into + shared-memory-sized blocks and only re-processing tiles whose boundaries + change. This approach is not currently implemented in cuCIM. + + Applications for grayscale reconstruction are discussed in [2]_ and [3]_. References ---------- @@ -136,49 +178,198 @@ def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): on Image Processing (1993) .. [3] Soille, P., "Morphological Image Analysis: Principles and Applications", Chapter 6, 2nd edition (2003), ISBN 3540429883. + .. [4] Teodoro, G., Banerjee, T., Kurc, T.M., Sussman, A., Pan, T., and + Saltz, J.H., "Efficient irregular wavefront propagation algorithms + on hybrid CPU-GPU machines", Parallel Computing, 39(4-5), + pp. 189-211, 2013. """ - from ..filters._rank_order import rank_order + seed = cp.asarray(seed) + mask = cp.asarray(mask) + + if reconstruct_on_cpu not in (False, True, "auto"): + raise ValueError( + "reconstruct_on_cpu must be one of False, True, or 'auto'" + ) + + if seed.shape != mask.shape: + raise ValueError( + f"seed and mask must have the same shape, " + f"got {seed.shape} and {mask.shape}" + ) + + if method not in ("dilation", "erosion"): + raise ValueError( + "Reconstruction method can be one of 'erosion' " + f"or 'dilation'. Got '{method}'." + ) - assert tuple(seed.shape) == tuple(mask.shape) if method == "dilation" and cp.any(seed > mask): # synchronize! raise ValueError( "Intensity of seed image must be less than that " "of the mask image for reconstruction by dilation." ) - elif method == "erosion" and cp.any(seed < mask): # synchronize! raise ValueError( "Intensity of seed image must be greater than that " "of the mask image for reconstruction by erosion." ) - try: - from skimage.morphology._grayreconstruct import reconstruction_loop - except ImportError: - try: - from skimage.morphology._greyreconstruct import reconstruction_loop - except ImportError: - raise ImportError("reconstruction requires scikit-image") - if footprint is None: - footprint = np.ones([3] * seed.ndim, dtype=bool) + footprint_shape = (3,) * seed.ndim + elif isinstance(footprint, cp.ndarray): + footprint_shape = footprint.shape else: - if isinstance(footprint, cp.ndarray): - footprint = cp.asnumpy(footprint) - footprint = footprint.astype(bool, copy=True) + footprint_shape = np.asarray(footprint).shape if offset is None: - if not all([d % 2 == 1 for d in footprint.shape]): + if not all([d % 2 == 1 for d in footprint_shape]): raise ValueError("Footprint dimensions must all be odd") - offset = np.array([d // 2 for d in footprint.shape]) + offset = np.array([d // 2 for d in footprint_shape]) else: if isinstance(offset, cp.ndarray): offset = cp.asnumpy(offset) - if offset.ndim != footprint.ndim: + else: + offset = np.asarray(offset) + if offset.ndim != len(footprint_shape): raise ValueError("Offset and footprint ndims must be equal.") - if not all([(0 <= o < d) for o, d in zip(offset, footprint.shape)]): + if not all([(0 <= o < d) for o, d in zip(offset, footprint_shape)]): raise ValueError("Offset must be included inside footprint") + if reconstruct_on_cpu == "auto": + reconstruct_on_cpu = not _gpu_reconstruction_supports_offset( + footprint_shape, offset + ) + elif not reconstruct_on_cpu and not _gpu_reconstruction_supports_offset( + footprint_shape, offset + ): + raise ValueError( + "The GPU reconstruction path cannot represent this footprint " + "offset. Use reconstruct_on_cpu=True or reconstruct_on_cpu='auto'." + ) + + if reconstruct_on_cpu: + if footprint is None: + footprint = np.ones(footprint_shape, dtype=bool) + elif isinstance(footprint, cp.ndarray): + footprint = cp.asnumpy(footprint).astype(bool, copy=True) + else: + footprint = np.asarray(footprint, dtype=bool).copy() + return _reconstruction_cpu(seed, mask, method, footprint, offset) + else: + if footprint is None: + footprint = cp.ones(footprint_shape, dtype=bool) + else: + footprint = cp.asarray(footprint, dtype=bool).copy() + return _reconstruction_gpu(seed, mask, method, footprint, offset) + + +def _gpu_reconstruction_supports_offset(footprint_shape, offset): + """Return whether ndimage's origin can express this footprint offset.""" + for size, off in zip(footprint_shape, offset): + origin = int(size // 2 - off) + if origin < -(size // 2) or origin > (size - 1) // 2: + return False + return True + + +def _extreme_cval(dtype, *, low): + """Return a boundary value that cannot win a min/max filter.""" + dtype = np.dtype(dtype) + if np.issubdtype(dtype, np.integer): + info = np.iinfo(dtype) + return info.min if low else info.max + if np.issubdtype(dtype, np.floating): + return -np.inf if low else np.inf + if np.issubdtype(dtype, np.bool_): + return False if low else True + raise TypeError(f"dtype {dtype} is not supported") + + +def _reconstruction_gpu(seed, mask, method, footprint, offset): + """Iterative GPU-native morphological reconstruction. + + Repeatedly applies conditional dilation (or erosion) until convergence. + Each iteration is a bulk GPU operation via maximum_filter/minimum_filter. + """ + # Work in a dtype that can represent the reconstruction + work_dtype = np.promote_types(seed.dtype, mask.dtype) + if method == "erosion" and np.issubdtype(work_dtype, np.unsignedinteger): + work_dtype = np.promote_types(work_dtype, np.int8) + + current = seed.astype(work_dtype, copy=True) + mask_work = mask.astype(work_dtype, copy=False) + + # Convert reconstruction's `offset` to ndimage's `origin` parameter. + # The CPU reconstruction defines offset as the footprint anchor: pixel i + # propagates TO neighbors at positions defined by (footprint_pos - offset). + # In the iterative GPU approach, maximum_filter has pixel i RECEIVE from + # its window. To match the CPU's propagation direction, mirror the + # footprint and reverse the relationship: + # origin = shape // 2 - offset (per dimension). + # For the default centered case (offset = shape // 2), origin = 0. + origin = tuple(int(d // 2 - o) for d, o in zip(footprint.shape, offset)) + footprint = cp.ascontiguousarray( + footprint[(slice(None, None, -1),) * footprint.ndim] + ) + + # Adaptive batch size: run several iterations between convergence checks + # to amortize the cost of the device synchronization in .any(). + check_every = 4 + max_check_every = 64 + max_iter = max(seed.size, 1000) # safety limit + + is_dilation = method == "dilation" + cval = _extreme_cval(current.dtype, low=is_dilation) + filter_kwargs = dict( + footprint=footprint, origin=origin, mode="constant", cval=cval + ) + n_iter = 0 + while n_iter < max_iter: + for _ in range(check_every): + if is_dilation: + dilated = ndi.maximum_filter(current, **filter_kwargs) + cp.minimum(dilated, mask_work, out=current) + else: + eroded = ndi.minimum_filter(current, **filter_kwargs) + cp.maximum(eroded, mask_work, out=current) + n_iter += 1 + + # Check convergence: would one more iteration change anything? + if is_dilation: + candidate = ndi.maximum_filter(current, **filter_kwargs) + cp.minimum(candidate, mask_work, out=candidate) + else: + candidate = ndi.minimum_filter(current, **filter_kwargs) + cp.maximum(candidate, mask_work, out=candidate) + + if not (candidate != current).any(): + break + + current = candidate + n_iter += 1 + + # Ramp up batch size if convergence is slow + check_every = min(check_every * 2, max_check_every) + + return current + + +def _reconstruction_cpu(seed, mask, method, footprint, offset): + """CPU-based morphological reconstruction using scikit-image's Cython loop. + + Data is transferred to host, processed single-threaded via + scikit-image's reconstruction_loop, then transferred back. + """ + from ..filters._rank_order import rank_order + + try: + from skimage.morphology._grayreconstruct import reconstruction_loop + except ImportError: + try: + from skimage.morphology._greyreconstruct import reconstruction_loop + except ImportError: + raise ImportError("reconstruction requires scikit-image") + # Cross out the center of the footprint footprint[tuple(slice(d, d + 1) for d in offset)] = False @@ -191,13 +382,9 @@ def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): # we can interleave image and mask pixels when sorting. if method == "dilation": pad_value = cp.min(seed).item() - elif method == "erosion": - pad_value = cp.max(seed).item() else: - raise ValueError( - "Reconstruction method can be one of 'erosion' " - f"or 'dilation'. Got '{method}'." - ) + pad_value = cp.max(seed).item() + # CuPy Backend: modified to allow images_dtype based on input dtype # instead of float64 images_dtype = np.promote_types(seed.dtype, mask.dtype) @@ -260,7 +447,6 @@ def reconstruction(seed, mask, method="dilation", footprint=None, offset=None): value_rank, value_map = rank_order(-images) value_map = -value_map - # TODO: implement reconstruction_loop on the GPU? For now, run it on host. start = index_sorted[0] value_rank = cp.asnumpy(value_rank.astype(unsigned_int_dtype, copy=False)) reconstruction_loop(value_rank, prev, next, nb_strides, start, image_stride) diff --git a/python/cucim/src/cucim/skimage/morphology/tests/test_extrema.py b/python/cucim/src/cucim/skimage/morphology/tests/test_extrema.py new file mode 100644 index 000000000..c8c686af2 --- /dev/null +++ b/python/cucim/src/cucim/skimage/morphology/tests/test_extrema.py @@ -0,0 +1,754 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +import math + +import cupy as cp +import numpy as np +import pytest + +from cucim.skimage.morphology import extrema, local_maxima, local_minima + +eps = 1e-12 + + +def diff(a, b): + a = cp.asarray(a, dtype=cp.float64) + b = cp.asarray(b, dtype=cp.float64) + t = float(((a - b) ** 2).sum()) + return math.sqrt(t) + + +class TestExtrema: + def test_saturated_arithmetic(self): + """Adding/subtracting a constant and clipping""" + # Test for unsigned integer + data = cp.array( + [[250, 251, 5, 5], [100, 200, 253, 252], [4, 10, 1, 3]], + dtype=cp.uint8, + ) + # adding the constant + img_constant_added = extrema._add_constant_clip(data, 4) + expected = cp.array( + [[254, 255, 9, 9], [104, 204, 255, 255], [8, 14, 5, 7]], + dtype=cp.uint8, + ) + error = diff(img_constant_added, expected) + assert error < eps + img_constant_subtracted = extrema._subtract_constant_clip(data, 4) + expected = cp.array( + [[246, 247, 1, 1], [96, 196, 249, 248], [0, 6, 0, 0]], + dtype=cp.uint8, + ) + error = diff(img_constant_subtracted, expected) + assert error < eps + + # Test for signed integer + data = cp.array([[32767, 32766], [-32768, -32767]], dtype=cp.int16) + img_constant_added = extrema._add_constant_clip(data, 1) + expected = cp.array([[32767, 32767], [-32767, -32766]], dtype=cp.int16) + error = diff(img_constant_added, expected) + assert error < eps + img_constant_subtracted = extrema._subtract_constant_clip(data, 1) + expected = cp.array([[32766, 32765], [-32768, -32768]], dtype=cp.int16) + error = diff(img_constant_subtracted, expected) + assert error < eps + + def test_h_maxima(self): + """h-maxima for various data types""" + + data = cp.array( + [ + [10, 11, 13, 14, 14, 15, 14, 14, 13, 11], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + [13, 15, 40, 40, 18, 18, 18, 60, 60, 15], + [14, 16, 40, 40, 19, 19, 19, 60, 60, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [15, 16, 18, 19, 19, 20, 19, 19, 18, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [14, 16, 80, 80, 19, 19, 19, 100, 100, 16], + [13, 15, 80, 80, 18, 18, 18, 100, 100, 15], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + ], + dtype=cp.uint8, + ) + + expected_result = cp.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=cp.uint8, + ) + for dtype in [cp.uint8, cp.uint64, cp.int8, cp.int64]: + data = data.astype(dtype) + out = extrema.h_maxima(data, 40) + + error = diff(expected_result, out) + assert error < eps + + def test_h_minima(self): + """h-minima for various data types""" + + data = cp.array( + [ + [10, 11, 13, 14, 14, 15, 14, 14, 13, 11], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + [13, 15, 40, 40, 18, 18, 18, 60, 60, 15], + [14, 16, 40, 40, 19, 19, 19, 60, 60, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [15, 16, 18, 19, 19, 20, 19, 19, 18, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [14, 16, 80, 80, 19, 19, 19, 100, 100, 16], + [13, 15, 80, 80, 18, 18, 18, 100, 100, 15], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + ], + dtype=cp.uint8, + ) + data = 100 - data + expected_result = cp.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=cp.uint8, + ) + for dtype in [cp.uint8, cp.uint64, cp.int8, cp.int64]: + data = data.astype(dtype) + out = extrema.h_minima(data, 40) + + error = diff(expected_result, out) + assert error < eps + assert out.dtype == expected_result.dtype + + def test_extrema_float(self): + """specific tests for float type""" + data = cp.array( + [ + [0.10, 0.11, 0.13, 0.14, 0.14, 0.15, 0.14, 0.14, 0.13, 0.11], + [0.11, 0.13, 0.15, 0.16, 0.16, 0.16, 0.16, 0.16, 0.15, 0.13], + [0.13, 0.15, 0.40, 0.40, 0.18, 0.18, 0.18, 0.60, 0.60, 0.15], + [0.14, 0.16, 0.40, 0.40, 0.19, 0.19, 0.19, 0.60, 0.60, 0.16], + [0.14, 0.16, 0.18, 0.19, 0.19, 0.19, 0.19, 0.19, 0.18, 0.16], + [0.15, 0.182, 0.18, 0.19, 0.204, 0.20, 0.19, 0.19, 0.18, 0.16], + [0.14, 0.16, 0.18, 0.19, 0.19, 0.19, 0.19, 0.19, 0.18, 0.16], + [0.14, 0.16, 0.80, 0.80, 0.19, 0.19, 0.19, 1.0, 1.0, 0.16], + [0.13, 0.15, 0.80, 0.80, 0.18, 0.18, 0.18, 1.0, 1.0, 0.15], + [0.11, 0.13, 0.15, 0.16, 0.16, 0.16, 0.16, 0.16, 0.15, 0.13], + ], + dtype=cp.float32, + ) + inverted_data = 1.0 - data + + out = extrema.h_maxima(data, 0.003) + expected_result = cp.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=cp.uint8, + ) + + error = diff(expected_result, out) + assert error < eps + + out = extrema.h_minima(inverted_data, 0.003) + error = diff(expected_result, out) + assert error < eps + + def test_h_maxima_float_image(self): + """specific tests for h-maxima float image type""" + w = 10 + x, y = cp.mgrid[0:w, 0:w] + data = 20 - 0.2 * ((x - w / 2) ** 2 + (y - w / 2) ** 2) + data[2:4, 2:4] = 40 + data[2:4, 7:9] = 60 + data[7:9, 2:4] = 80 + data[7:9, 7:9] = 100 + data = data.astype(cp.float32) + + expected_result = cp.zeros_like(data) + expected_result[(data > 19.9)] = 1.0 + + for h in [1.0e-12, 1.0e-6, 1.0e-3, 1.0e-2, 1.0e-1, 0.1]: + out = extrema.h_maxima(data, h) + error = diff(expected_result, out) + assert error < eps + + def test_h_maxima_float_h(self): + """specific tests for h-maxima float h parameter""" + data = cp.array( + [ + [0, 0, 0, 0, 0], + [0, 3, 3, 3, 0], + [0, 3, 4, 3, 0], + [0, 3, 3, 3, 0], + [0, 0, 0, 0, 0], + ], + dtype=cp.uint8, + ) + + h_vals = np.linspace(1.0, 2.0, 100) + failures = 0 + + for h in h_vals: + if h % 1 != 0: + with pytest.warns( + UserWarning, + match="possible precision loss converting image", + ): + maxima = extrema.h_maxima(data, h) + else: + maxima = extrema.h_maxima(data, h) + + if maxima[2, 2] == 0: + failures += 1 + + assert failures == 0 + + def test_h_maxima_large_h(self): + """test that h-maxima works correctly for large h""" + data = cp.array( + [ + [10, 10, 10, 10, 10], + [10, 13, 13, 13, 10], + [10, 13, 14, 13, 10], + [10, 13, 13, 13, 10], + [10, 10, 10, 10, 10], + ], + dtype=cp.uint8, + ) + + maxima = extrema.h_maxima(data, 5) + assert cp.sum(maxima) == 0 + + data = cp.array( + [ + [10, 10, 10, 10, 10], + [10, 13, 13, 13, 10], + [10, 13, 14, 13, 10], + [10, 13, 13, 13, 10], + [10, 10, 10, 10, 10], + ], + dtype=cp.float32, + ) + + maxima = extrema.h_maxima(data, 5.0) + assert cp.sum(maxima) == 0 + + def test_h_minima_float_image(self): + """specific tests for h-minima float image type""" + w = 10 + x, y = cp.mgrid[0:w, 0:w] + data = 180 + 0.2 * ((x - w / 2) ** 2 + (y - w / 2) ** 2) + data[2:4, 2:4] = 160 + data[2:4, 7:9] = 140 + data[7:9, 2:4] = 120 + data[7:9, 7:9] = 100 + data = data.astype(cp.float32) + + expected_result = cp.zeros_like(data) + expected_result[(data < 180.1)] = 1.0 + + for h in [1.0e-12, 1.0e-6, 1.0e-3, 1.0e-2, 1.0e-1, 0.1]: + out = extrema.h_minima(data, h) + error = diff(expected_result, out) + assert error < eps + + def test_h_minima_float_h(self): + """specific tests for h-minima float h parameter""" + data = cp.array( + [ + [4, 4, 4, 4, 4], + [4, 1, 1, 1, 4], + [4, 1, 0, 1, 4], + [4, 1, 1, 1, 4], + [4, 4, 4, 4, 4], + ], + dtype=cp.uint8, + ) + + h_vals = np.linspace(1.0, 2.0, 100) + failures = 0 + for h in h_vals: + if h % 1 != 0: + with pytest.warns( + UserWarning, + match="possible precision loss converting image", + ): + minima = extrema.h_minima(data, h) + else: + minima = extrema.h_minima(data, h) + + if minima[2, 2] == 0: + failures += 1 + + assert failures == 0 + + def test_h_minima_large_h(self): + """test that h-minima works correctly for large h""" + data = cp.array( + [ + [14, 14, 14, 14, 14], + [14, 11, 11, 11, 14], + [14, 11, 10, 11, 14], + [14, 11, 11, 11, 14], + [14, 14, 14, 14, 14], + ], + dtype=cp.uint8, + ) + + maxima = extrema.h_minima(data, 5) + assert cp.sum(maxima) == 0 + + data = cp.array( + [ + [14, 14, 14, 14, 14], + [14, 11, 11, 11, 14], + [14, 11, 10, 11, 14], + [14, 11, 11, 11, 14], + [14, 14, 14, 14, 14], + ], + dtype=cp.float32, + ) + + maxima = extrema.h_minima(data, 5.0) + assert cp.sum(maxima) == 0 + + +class TestLocalMaxima: + """Some tests for local_minima are included as well.""" + + supported_dtypes = [ + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.int8, + np.int16, + np.int32, + np.int64, + np.float32, + np.float64, + ] + image = cp.asarray( + [ + [1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 2, 0, 0, 3, 3, 0, 0, 4, 0, 2, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 4, 4, 0, 3, 0, 0, 0], + [0, 2, 0, 1, 0, 2, 1, 0, 0, 0, 0, 3, 0, 0, 0], + [0, 0, 2, 0, 2, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0], + ], + dtype=np.uint8, + ) + # Connectivity 2, maxima can touch border, returned with default values + expected_default = cp.asarray( + [ + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + # Connectivity 1 (cross), maxima can touch border + expected_cross = cp.asarray( + [ + [1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0], + [0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + + def test_empty(self): + """Test result with empty image.""" + result = local_maxima(cp.asarray([[]]), indices=False) + assert result.size == 0 + assert result.dtype == bool + assert result.shape == (1, 0) + + result = local_maxima(cp.asarray([]), indices=True) + assert isinstance(result, tuple) + assert len(result) == 1 + assert result[0].size == 0 + assert result[0].dtype == np.intp + + result = local_maxima(cp.asarray([[]]), indices=True) + assert isinstance(result, tuple) + assert len(result) == 2 + assert result[0].size == 0 + assert result[0].dtype == np.intp + assert result[1].size == 0 + assert result[1].dtype == np.intp + + def test_dtypes(self): + """Test results with default configuration for all supported dtypes.""" + for dtype in self.supported_dtypes: + result = local_maxima(self.image.astype(dtype)) + assert result.dtype == bool + cp.testing.assert_array_equal(result, self.expected_default) + + def test_dtypes_old(self): + """ + Test results with default configuration and data copied from old unit + tests for all supported dtypes. + """ + data = cp.asarray( + [ + [10, 11, 13, 14, 14, 15, 14, 14, 13, 11], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + [13, 15, 40, 40, 18, 18, 18, 60, 60, 15], + [14, 16, 40, 40, 19, 19, 19, 60, 60, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [15, 16, 18, 19, 19, 20, 19, 19, 18, 16], + [14, 16, 18, 19, 19, 19, 19, 19, 18, 16], + [14, 16, 80, 80, 19, 19, 19, 100, 100, 16], + [13, 15, 80, 80, 18, 18, 18, 100, 100, 15], + [11, 13, 15, 16, 16, 16, 16, 16, 15, 13], + ], + dtype=np.uint8, + ) + expected = cp.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + for dtype in self.supported_dtypes: + image = data.astype(dtype) + result = local_maxima(image) + assert result.dtype == bool + cp.testing.assert_array_equal(result, expected) + + def test_connectivity(self): + """Test results if footprint is a scalar.""" + # Connectivity 1: generates cross shaped footprint + result_conn1 = local_maxima(self.image, connectivity=1) + assert result_conn1.dtype == bool + cp.testing.assert_array_equal(result_conn1, self.expected_cross) + + # Connectivity 2: generates square shaped footprint + result_conn2 = local_maxima(self.image, connectivity=2) + assert result_conn2.dtype == bool + cp.testing.assert_array_equal(result_conn2, self.expected_default) + + # Connectivity 3: generates square shaped footprint + result_conn3 = local_maxima(self.image, connectivity=3) + assert result_conn3.dtype == bool + cp.testing.assert_array_equal(result_conn3, self.expected_default) + + def test_footprint(self): + """Test results if footprint is given.""" + footprint_cross = cp.asarray( + [[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=bool + ) + result_footprint_cross = local_maxima( + self.image, footprint=footprint_cross + ) + assert result_footprint_cross.dtype == bool + cp.testing.assert_array_equal( + result_footprint_cross, self.expected_cross + ) + + for footprint in [ + ((True,) * 3,) * 3, + cp.ones((3, 3), dtype=np.float64), + cp.ones((3, 3), dtype=np.uint8), + cp.ones((3, 3), dtype=bool), + ]: + # Test different dtypes for footprint which expects a boolean array + # but will accept and convert other types if possible + result_footprint_square = local_maxima( + self.image, footprint=footprint + ) + assert result_footprint_square.dtype == bool + cp.testing.assert_array_equal( + result_footprint_square, self.expected_default + ) + + footprint_x = cp.asarray([[1, 0, 1], [0, 1, 0], [1, 0, 1]], dtype=bool) + expected_footprint_x = cp.asarray( + [ + [1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + result_footprint_x = local_maxima(self.image, footprint=footprint_x) + assert result_footprint_x.dtype == bool + cp.testing.assert_array_equal(result_footprint_x, expected_footprint_x) + + def test_indices(self): + """Test output if indices of peaks are desired.""" + # Connectivity 1 + expected_conn1 = cp.nonzero(self.expected_cross) + result_conn1 = local_maxima(self.image, connectivity=1, indices=True) + assert len(result_conn1) == len(expected_conn1) + for r, e in zip(result_conn1, expected_conn1): + cp.testing.assert_array_equal(r, e) + + # Connectivity 2 + expected_conn2 = cp.nonzero(self.expected_default) + result_conn2 = local_maxima(self.image, connectivity=2, indices=True) + assert len(result_conn2) == len(expected_conn2) + for r, e in zip(result_conn2, expected_conn2): + cp.testing.assert_array_equal(r, e) + + def test_allow_borders(self): + """Test maxima detection at the image border.""" + # Use connectivity 1 to allow many maxima, only filtering at border is + # of interest + result_with_border = local_maxima( + self.image, connectivity=1, allow_borders=True + ) + assert result_with_border.dtype == bool + cp.testing.assert_array_equal(result_with_border, self.expected_cross) + + expected_without_border = cp.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0], + [0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + result_without_border = local_maxima( + self.image, connectivity=1, allow_borders=False + ) + assert result_with_border.dtype == bool + cp.testing.assert_array_equal( + result_without_border, expected_without_border + ) + + def test_nd(self): + """Test one- and three-dimensional case.""" + # One-dimension + x_1d = cp.asarray([1, 1, 0, 1, 2, 3, 0, 2, 1, 2, 0]) + expected_1d = np.asarray([1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0], dtype=bool) + result_1d = local_maxima(x_1d) + assert result_1d.dtype == bool + cp.testing.assert_array_equal(result_1d, expected_1d) + + # 3-dimensions (adapted from old unit test) + x_3d = cp.zeros((8, 8, 8), dtype=np.uint8) + expected_3d = cp.zeros((8, 8, 8), dtype=bool) + # first maximum: only one pixel + x_3d[1, 1:3, 1:3] = 100 + x_3d[2, 2, 2] = 200 + x_3d[3, 1:3, 1:3] = 100 + expected_3d[2, 2, 2] = 1 + # second maximum: three pixels in z-direction + x_3d[5:8, 1, 1] = 200 + expected_3d[5:8, 1, 1] = 1 + # third: two maxima in 0 and 3. + x_3d[0, 5:8, 5:8] = 200 + x_3d[1, 6, 6] = 100 + x_3d[2, 5:7, 5:7] = 200 + x_3d[0:3, 5:8, 5:8] += 50 + expected_3d[0, 5:8, 5:8] = 1 + expected_3d[2, 5:7, 5:7] = 1 + # four : one maximum in the corner of the square + x_3d[6:8, 6:8, 6:8] = 200 + x_3d[7, 7, 7] = 255 + expected_3d[7, 7, 7] = 1 + result_3d = local_maxima(x_3d) + assert result_3d.dtype == bool + cp.testing.assert_array_equal(result_3d, expected_3d) + + def test_constant(self): + """Test behaviour for 'flat' images.""" + const_image = cp.full((7, 6), 42, dtype=np.uint8) + expected = cp.zeros((7, 6), dtype=np.uint8) + for dtype in self.supported_dtypes: + const_image = const_image.astype(dtype) + # test for local maxima + result = local_maxima(const_image) + assert result.dtype == bool + cp.testing.assert_array_equal(result, expected) + # test for local minima + result = local_minima(const_image) + assert result.dtype == bool + cp.testing.assert_array_equal(result, expected) + + def test_extrema_float(self): + """Specific tests for float type.""" + # Copied from old unit test for local_maxima + image = cp.asarray( + [ + [0.10, 0.11, 0.13, 0.14, 0.14, 0.15, 0.14, 0.14, 0.13, 0.11], + [0.11, 0.13, 0.15, 0.16, 0.16, 0.16, 0.16, 0.16, 0.15, 0.13], + [0.13, 0.15, 0.40, 0.40, 0.18, 0.18, 0.18, 0.60, 0.60, 0.15], + [0.14, 0.16, 0.40, 0.40, 0.19, 0.19, 0.19, 0.60, 0.60, 0.16], + [0.14, 0.16, 0.18, 0.19, 0.19, 0.19, 0.19, 0.19, 0.18, 0.16], + [0.15, 0.182, 0.18, 0.19, 0.204, 0.20, 0.19, 0.19, 0.18, 0.16], + [0.14, 0.16, 0.18, 0.19, 0.19, 0.19, 0.19, 0.19, 0.18, 0.16], + [0.14, 0.16, 0.80, 0.80, 0.19, 0.19, 0.19, 1.0, 1.0, 0.16], + [0.13, 0.15, 0.80, 0.80, 0.18, 0.18, 0.18, 1.0, 1.0, 0.15], + [0.11, 0.13, 0.15, 0.16, 0.16, 0.16, 0.16, 0.16, 0.15, 0.13], + ], + dtype=np.float32, + ) + inverted_image = 1.0 - image + expected_result = cp.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=bool, + ) + + # Test for local maxima with automatic step calculation + result = local_maxima(image) + assert result.dtype == bool + cp.testing.assert_array_equal(result, expected_result) + + # Test for local minima with automatic step calculation + result = local_minima(inverted_image) + assert result.dtype == bool + cp.testing.assert_array_equal(result, expected_result) + + def test_extrema_small_float(self): + image = cp.asarray( + [ + [ + 9.89232736e-20, + 8.78543302e-20, + 8.78543302e-20, + 9.89232736e-20, + ], + [ + 8.78543302e-20, + 6.38842355e-20, + 6.38842355e-20, + 8.78543302e-20, + ], + [ + 8.78543302e-20, + 6.38842355e-20, + 6.38842355e-20, + 8.78543302e-20, + ], + [ + 9.89232736e-20, + 8.78543302e-20, + 8.78543302e-20, + 9.89232736e-20, + ], + ] + ) + + result = local_minima(image) + + expected_result = cp.asarray( + [ + [False, False, False, False], + [False, True, True, False], + [False, True, True, False], + [False, False, False, False], + ] + ) + + cp.testing.assert_array_equal(result, expected_result) + + def test_exceptions(self): + """Test if input validation triggers correct exceptions.""" + # Mismatching number of dimensions + with pytest.raises(ValueError, match="number of dimensions"): + local_maxima(self.image, footprint=cp.ones((3, 3, 3), dtype=bool)) + with pytest.raises(ValueError, match="number of dimensions"): + local_maxima(self.image, footprint=cp.ones((3,), dtype=bool)) + + # All dimensions in footprint must be of size 3 + with pytest.raises(ValueError, match="dimension size"): + local_maxima(self.image, footprint=cp.ones((2, 3), dtype=bool)) + with pytest.raises(ValueError, match="dimension size"): + local_maxima(self.image, footprint=cp.ones((5, 5), dtype=bool)) + + with pytest.raises(TypeError, match="float16 which is not supported"): + local_maxima(cp.empty(1, dtype=np.float16)) + + def test_small_array(self): + """Test output for arrays with dimension smaller 3. + + If any dimension of an array is smaller than 3 and `allow_borders` is + false a footprint, which has at least 3 elements in each + dimension, can't be applied. This is an implementation detail so + `local_maxima` should still return valid output (see gh-3261). + + If `allow_borders` is true the array is padded internally and there is + no problem. + """ + warning_msg = "maxima can't exist .* any dimension smaller 3 .*" + x = cp.asarray([0, 1]) + local_maxima(x, allow_borders=True) # no warning + with pytest.warns(UserWarning, match=warning_msg): + result = local_maxima(x, allow_borders=False) + cp.testing.assert_array_equal(result, [0, 0]) + assert result.dtype == bool + + x = cp.asarray([[1, 2], [2, 2]]) + local_maxima(x, allow_borders=True, indices=True) # no warning + with pytest.warns(UserWarning, match=warning_msg): + result = local_maxima(x, allow_borders=False, indices=True) + assert isinstance(result, tuple) + assert len(result) == 2 + assert result[0].size == 0 + assert result[1].size == 0 + assert result[0].dtype == np.intp + assert result[1].dtype == np.intp diff --git a/python/cucim/src/cucim/skimage/morphology/tests/test_reconstruction.py b/python/cucim/src/cucim/skimage/morphology/tests/test_reconstruction.py index 072399524..a12fb3464 100644 --- a/python/cucim/src/cucim/skimage/morphology/tests/test_reconstruction.py +++ b/python/cucim/src/cucim/skimage/morphology/tests/test_reconstruction.py @@ -1,6 +1,6 @@ # SPDX-FileCopyrightText: Copyright (c) 2003-2009 Massachusetts Institute of Technology # SPDX-FileCopyrightText: Copyright (c) 2009-2011 Broad Institute -# SPDX-FileCopyrightText: Copyright (c) 2021-2025, NVIDIA CORPORATION. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2021-2026, NVIDIA CORPORATION. All rights reserved. # SPDX-License-Identifier: Apache-2.0 AND (GPL-2.0-only OR BSD-3-Clause) """ @@ -25,36 +25,55 @@ from cucim.skimage.morphology import reconstruction -def test_zeros(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_zeros(reconstruct_on_cpu): """Test reconstruction with image and mask of zeros""" assert_array_almost_equal( - reconstruction(cp.zeros((5, 7)), cp.zeros((5, 7))), 0 + reconstruction( + cp.zeros((5, 7)), + cp.zeros((5, 7)), + reconstruct_on_cpu=reconstruct_on_cpu, + ), + 0, ) -def test_image_equals_mask(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_image_equals_mask(reconstruct_on_cpu): """Test reconstruction where the image and mask are the same""" assert_array_almost_equal( - reconstruction(cp.ones((7, 5)), cp.ones((7, 5))), 1 + reconstruction( + cp.ones((7, 5)), + cp.ones((7, 5)), + reconstruct_on_cpu=reconstruct_on_cpu, + ), + 1, ) -def test_image_less_than_mask(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_image_less_than_mask(reconstruct_on_cpu): """Test reconstruction where the image is uniform and less than mask""" image = cp.ones((5, 5)) mask = cp.ones((5, 5)) * 2 - assert_array_almost_equal(reconstruction(image, mask), 1) + assert_array_almost_equal( + reconstruction(image, mask, reconstruct_on_cpu=reconstruct_on_cpu), 1 + ) -def test_one_image_peak(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_one_image_peak(reconstruct_on_cpu): """Test reconstruction with one peak pixel""" image = cp.ones((5, 5)) image[2, 2] = 2 mask = cp.ones((5, 5)) * 3 - assert_array_almost_equal(reconstruction(image, mask), 2) + assert_array_almost_equal( + reconstruction(image, mask, reconstruct_on_cpu=reconstruct_on_cpu), 2 + ) -def test_two_image_peaks(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_two_image_peaks(reconstruct_on_cpu): """Test reconstruction with two peak pixels isolated by the mask""" # fmt: off image = cp.array([[1, 1, 1, 1, 1, 1, 1, 1], @@ -78,20 +97,31 @@ def test_two_image_peaks(): [1, 1, 1, 1, 1, 3, 3, 3], [1, 1, 1, 1, 1, 3, 3, 3]]) # fmt: on - assert_array_almost_equal(reconstruction(image, mask), expected) + assert_array_almost_equal( + reconstruction(image, mask, reconstruct_on_cpu=reconstruct_on_cpu), + expected, + ) -def test_zero_image_one_mask(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_zero_image_one_mask(reconstruct_on_cpu): """Test reconstruction with an image of all zeros and a mask that's not""" - result = reconstruction(cp.zeros((10, 10)), cp.ones((10, 10))) + result = reconstruction( + cp.zeros((10, 10)), + cp.ones((10, 10)), + reconstruct_on_cpu=reconstruct_on_cpu, + ) assert_array_almost_equal(result, 0) -def test_fill_hole(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_fill_hole(reconstruct_on_cpu): """Test reconstruction by erosion, which should fill holes in mask.""" seed = cp.array([0, 8, 8, 8, 8, 8, 8, 8, 8, 0]) mask = cp.array([0, 3, 6, 2, 1, 1, 1, 4, 2, 0]) - result = reconstruction(seed, mask, method="erosion") + result = reconstruction( + seed, mask, method="erosion", reconstruct_on_cpu=reconstruct_on_cpu + ) assert_array_almost_equal(result, cp.array([0, 3, 6, 4, 4, 4, 4, 4, 2, 0])) @@ -104,10 +134,11 @@ def test_invalid_seed(): reconstruction(seed * 0.5, mask, method="erosion") +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) @pytest.mark.parametrize( "unsigned_dtype", [cp.uint8, cp.uint16, cp.uint32, cp.uint64] ) -def test_erosion_uint8_input(unsigned_dtype): +def test_erosion_uint8_input(unsigned_dtype, reconstruct_on_cpu): image = cp.asarray(data.moon()) # Rescale image intensity so that we can see dim features. image = rescale_intensity(image, in_range=(50, 200)) @@ -117,7 +148,9 @@ def test_erosion_uint8_input(unsigned_dtype): mask = image image2 = image2.astype(unsigned_dtype) - image2_erosion = reconstruction(image2, mask, method="erosion") + image2_erosion = reconstruction( + image2, mask, method="erosion", reconstruct_on_cpu=reconstruct_on_cpu + ) expected_out_type = cp.promote_types(image2.dtype, cp.int8) assert image2_erosion.dtype == expected_out_type @@ -137,24 +170,44 @@ def test_erosion_uint8_input(unsigned_dtype): ) -def test_invalid_footprint(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_invalid_footprint(reconstruct_on_cpu): seed = cp.ones((5, 5)) mask = cp.ones((5, 5)) with pytest.raises(ValueError): - reconstruction(seed, mask, footprint=np.ones((4, 4))) + reconstruction( + seed, + mask, + footprint=cp.ones((4, 4)), + reconstruct_on_cpu=reconstruct_on_cpu, + ) with pytest.raises(ValueError): - reconstruction(seed, mask, footprint=np.ones((3, 4))) - reconstruction(seed, mask, footprint=np.ones((3, 3))) + reconstruction( + seed, + mask, + footprint=cp.ones((3, 4)), + reconstruct_on_cpu=reconstruct_on_cpu, + ) + reconstruction( + seed, + mask, + footprint=cp.ones((3, 3)), + reconstruct_on_cpu=reconstruct_on_cpu, + ) -def test_invalid_method(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_invalid_method(reconstruct_on_cpu): seed = cp.array([0, 8, 8, 8, 8, 8, 8, 8, 8, 0]) mask = cp.array([0, 3, 6, 2, 1, 1, 1, 4, 2, 0]) with pytest.raises(ValueError): - reconstruction(seed, mask, method="foo") + reconstruction( + seed, mask, method="foo", reconstruct_on_cpu=reconstruct_on_cpu + ) -def test_invalid_offset_not_none(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_invalid_offset_not_none(reconstruct_on_cpu): """Test reconstruction with invalid not None offset parameter""" # fmt: off image = cp.array([[1, 1, 1, 1, 1, 1, 1, 1], @@ -178,10 +231,12 @@ def test_invalid_offset_not_none(): method="dilation", footprint=cp.ones((3, 3)), offset=cp.array([3, 0]), + reconstruct_on_cpu=reconstruct_on_cpu, ) -def test_offset_not_none(): +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_offset_not_none(reconstruct_on_cpu): """Test reconstruction with valid offset parameter""" seed = cp.array([0, 3, 6, 2, 1, 1, 1, 4, 2, 0]) mask = cp.array([0, 8, 6, 8, 8, 8, 8, 4, 4, 0]) @@ -194,12 +249,130 @@ def test_offset_not_none(): method="dilation", footprint=cp.ones(3), offset=cp.array([0]), + reconstruct_on_cpu=reconstruct_on_cpu, + ), + expected, + ) + + +@pytest.mark.parametrize("reconstruct_on_cpu", ["auto", False, True]) +def test_noncentered_offset_matches_skimage(reconstruct_on_cpu): + seed = cp.array([0, 1, 0, 2, 0, 2]) + mask = cp.array([2, 3, 0, 3, 1, 3]) + footprint = cp.ones(3) + offset = cp.array([0]) + expected = cp.array([0, 1, 0, 2, 1, 2]) + + assert_array_almost_equal( + reconstruction( + seed, + mask, + footprint=footprint, + offset=offset, + reconstruct_on_cpu=reconstruct_on_cpu, ), expected, ) -def test_reconstruction_float_inputs(): +def test_asymmetric_footprint_matches_skimage_gpu_path(): + seed = cp.array([0, 2, 0, 0, 0]) + mask = cp.array([0, 3, 3, 3, 0]) + footprint = cp.array([0, 1, 1]) + offset = cp.array([1]) + + result = reconstruction( + seed, + mask, + footprint=footprint, + offset=offset, + reconstruct_on_cpu=False, + ) + expected = reconstruction_cpu( + cp.asnumpy(seed), + cp.asnumpy(mask), + footprint=cp.asnumpy(footprint), + offset=cp.asnumpy(offset), + ) + + cp.testing.assert_allclose( + result, cp.asarray(expected).astype(result.dtype) + ) + + +def test_auto_falls_back_for_unsupported_gpu_offset(): + seed = cp.array([0, 1, 0, 2, 0, 2]) + mask = cp.array([2, 3, 0, 3, 1, 3]) + footprint = cp.ones(2) + offset = cp.array([0]) + + result = reconstruction( + seed, + mask, + footprint=footprint, + offset=offset, + reconstruct_on_cpu="auto", + ) + expected = reconstruction_cpu( + cp.asnumpy(seed), + cp.asnumpy(mask), + footprint=cp.asnumpy(footprint), + offset=cp.asnumpy(offset), + ) + cp.testing.assert_allclose( + result, cp.asarray(expected).astype(result.dtype) + ) + + with pytest.raises(ValueError, match="GPU reconstruction path"): + reconstruction( + seed, + mask, + footprint=footprint, + offset=offset, + reconstruct_on_cpu=False, + ) + + +@pytest.mark.parametrize( + "dtype", + [ + cp.int8, + cp.uint8, + cp.int16, + cp.uint16, + cp.int32, + cp.uint32, + cp.int64, + cp.uint64, + cp.float32, + cp.float64, + ], +) +@pytest.mark.parametrize("method", ["dilation", "erosion"]) +def test_gpu_cval_dtype_extremes(dtype, method): + dtype = cp.dtype(dtype) + if dtype.kind in "iu": + info = cp.iinfo(dtype) + low = info.min + high = info.max + else: + low = -1.0 + high = 1.0 + + if method == "dilation": + seed = cp.array([low, low, low], dtype=dtype) + mask = cp.array([low, high, low], dtype=dtype) + else: + seed = cp.array([high, high, high], dtype=dtype) + mask = cp.array([high, low, high], dtype=dtype) + + result = reconstruction(seed, mask, method=method, reconstruct_on_cpu=False) + expected_value = low if method == "dilation" else high + cp.testing.assert_array_equal(result, cp.full_like(result, expected_value)) + + +@pytest.mark.parametrize("reconstruct_on_cpu", [False, True]) +def test_reconstruction_float_inputs(reconstruct_on_cpu): """Verifies fix for: https://github.com/rapidsai/cuci/issues/36 Run the 2D example from the reconstruction docstring and compare the output @@ -211,5 +384,9 @@ def test_reconstruction_float_inputs(): h = 0.3 seed = bumps - h background_cpu = reconstruction_cpu(seed, bumps) - background = reconstruction(cp.asarray(seed), cp.asarray(bumps)) + background = reconstruction( + cp.asarray(seed), + cp.asarray(bumps), + reconstruct_on_cpu=reconstruct_on_cpu, + ) cp.testing.assert_allclose(background, background_cpu) diff --git a/python/cucim/src/cucim/skimage/segmentation/_clear_border.py b/python/cucim/src/cucim/skimage/segmentation/_clear_border.py index c83cb3d91..a7daacf17 100644 --- a/python/cucim/src/cucim/skimage/segmentation/_clear_border.py +++ b/python/cucim/src/cucim/skimage/segmentation/_clear_border.py @@ -1,11 +1,44 @@ # SPDX-FileCopyrightText: 2009-2022 the scikit-image team -# SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION. All rights reserved. # SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause import cupy as cp from ..measure import label + +def _get_border_labels(labels, buffer_size=0): + """Return unique labels that touch the image border. + + Parameters + ---------- + labels : cupy.ndarray + Label image. + buffer_size : int, optional + The width of the border examined. By default (0), only labels + that touch the outermost pixels are returned. + + Returns + ------- + border_labels : cupy.ndarray + 1-D array of unique label values that touch the border region. + """ + # Create border mask using slice-based approach + border_mask = cp.zeros(labels.shape, dtype=bool) + ext = buffer_size + 1 + slstart = slice(ext) + slend = slice(-ext, None) + slices = [slice(None)] * labels.ndim + for d in range(labels.ndim): + slices[d] = slstart + border_mask[tuple(slices)] = True + slices[d] = slend + border_mask[tuple(slices)] = True + slices[d] = slice(None) + + return cp.unique(labels[border_mask]) + + _clear_border_labels = cp.ElementwiseKernel( in_params="raw X labels, raw X borders_indices, int32 nvals, Y bgval", out_params="Y out", @@ -89,6 +122,10 @@ def clear_border(labels, buffer_size=0, bgval=0, mask=None, *, out=None): else: out = labels.copy() + # Re-label, in case we are dealing with a binary out + # and to get consistent labeling + relabeled, number = label(out, background=0, return_num=True) + if mask is not None: err_msg = ( f"labels and mask should have the same shape but " @@ -98,29 +135,14 @@ def clear_border(labels, buffer_size=0, bgval=0, mask=None, *, out=None): raise ValueError(err_msg) if mask.dtype != bool: raise TypeError("mask should be of type bool.") + # For mask case, border is where mask is False borders = ~mask + borders_indices = cp.unique(relabeled[borders]) else: - # create borders with buffer_size - borders = cp.zeros_like(out, dtype=bool) - ext = buffer_size + 1 - slstart = slice(ext) - slend = slice(-ext, None) - slices = [slice(None) for _ in out.shape] - for d in range(out.ndim): - slices[d] = slstart - borders[tuple(slices)] = True - slices[d] = slend - borders[tuple(slices)] = True - slices[d] = slice(None) - - # Re-label, in case we are dealing with a binary out - # and to get consistent labeling - labels, number = label(out, background=0, return_num=True) - - # determine all objects that are connected to borders - borders_indices = cp.unique(labels[borders]) + # Use helper function to get labels touching border + borders_indices = _get_border_labels(relabeled, buffer_size=buffer_size) _clear_border_labels( - labels, borders_indices, borders_indices.size, bgval, out + relabeled, borders_indices, borders_indices.size, bgval, out ) return out