diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 54bc258ac..ffeaae699 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -149,6 +149,11 @@ repos: ^python/cucim/src/cucim/skimage/filters/_window[.]py$| ^python/cucim/src/cucim/skimage/filters/cuda/histogram_median[.]cu$| ^python/cucim/src/cucim/skimage/filters/lpi_filter[.]py$| + ^python/cucim/src/cucim/skimage/filters/rank/__init__[.]py$| + ^python/cucim/src/cucim/skimage/filters/rank/_generic[.]py$| + ^python/cucim/src/cucim/skimage/filters/rank/_percentile[.]py$| + ^python/cucim/src/cucim/skimage/filters/rank/_rank_filter[.]py$| + ^python/cucim/src/cucim/skimage/filters/rank/tests/test_rank[.]py$| ^python/cucim/src/cucim/skimage/filters/ridges[.]py$| ^python/cucim/src/cucim/skimage/filters/tests/test_correlate[.]py$| ^python/cucim/src/cucim/skimage/filters/tests/test_edges[.]py$| @@ -434,6 +439,11 @@ repos: python/cucim/src/cucim/skimage/filters/_window[.]py$| python/cucim/src/cucim/skimage/filters/lpi_filter[.]py$| python/cucim/src/cucim/skimage/filters/ridges[.]py$| + python/cucim/src/cucim/skimage/filters/rank/__init__[.]py$| + python/cucim/src/cucim/skimage/filters/rank/_generic[.]py$| + python/cucim/src/cucim/skimage/filters/rank/_percentile[.]py$| + python/cucim/src/cucim/skimage/filters/rank/_rank_filter[.]py$| + python/cucim/src/cucim/skimage/filters/rank/tests/test_rank[.]py$| python/cucim/src/cucim/skimage/filters/tests/test_correlate[.]py$| python/cucim/src/cucim/skimage/filters/tests/test_edges[.]py$| python/cucim/src/cucim/skimage/filters/tests/test_fft_based[.]py$| diff --git a/benchmarks/skimage/_image_bench.py b/benchmarks/skimage/_image_bench.py index aab6b84ba..0a1672f6d 100644 --- a/benchmarks/skimage/_image_bench.py +++ b/benchmarks/skimage/_image_bench.py @@ -40,11 +40,17 @@ def __init__( module_gpu=cupyx.scipy.ndimage, function_is_generator=False, run_cpu=True, + fixed_kwargs_cpu=None, + fixed_kwargs_gpu=None, ): self.shape = shape self.function_name = function_name - self.fixed_kwargs_cpu = self._update_kwargs_arrays(fixed_kwargs, "cpu") - self.fixed_kwargs_gpu = self._update_kwargs_arrays(fixed_kwargs, "gpu") + if fixed_kwargs_cpu is None: + fixed_kwargs_cpu = fixed_kwargs + if fixed_kwargs_gpu is None: + fixed_kwargs_gpu = fixed_kwargs + self.fixed_kwargs_cpu = self._update_kwargs_arrays(fixed_kwargs_cpu, "cpu") + self.fixed_kwargs_gpu = self._update_kwargs_arrays(fixed_kwargs_gpu, "gpu") self.var_kwargs = var_kwargs self.index_str = index_str # self.set_args_kwargs = set_args_kwargs diff --git a/benchmarks/skimage/cucim_filters_rank_bench.py b/benchmarks/skimage/cucim_filters_rank_bench.py new file mode 100644 index 000000000..9d1380669 --- /dev/null +++ b/benchmarks/skimage/cucim_filters_rank_bench.py @@ -0,0 +1,213 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +import argparse +import os + +import numpy as np +import pandas as pd +import skimage.filters.rank +from _image_bench import ImageBench +from skimage.morphology import disk + +import cucim.skimage.filters.rank + +RANK_FILTERS = [ + # generic.py + ("autolevel", dict(), dict()), + ("enhance_contrast", dict(), dict()), + ("entropy", dict(), dict()), + ("equalize", dict(), dict()), + ("geometric_mean", dict(), dict()), + ("gradient", dict(), dict()), + ("majority", dict(), dict()), + ("maximum", dict(), dict()), + ("mean", dict(), dict()), + ("median", dict(), dict()), + ("minimum", dict(), dict()), + ("modal", dict(), dict()), + ("noise_filter", dict(), dict()), + ("pop", dict(), dict()), + ("subtract_mean", dict(), dict()), + ("sum", dict(), dict()), + ("threshold", dict(), dict()), + # percentile.py + ("autolevel_percentile", dict(p0=0.1, p1=0.9), dict()), + ("enhance_contrast_percentile", dict(p0=0.1, p1=0.9), dict()), + ("gradient_percentile", dict(p0=0.1, p1=0.9), dict()), + ("mean_percentile", dict(p0=0.1, p1=0.9), dict()), + ("percentile", dict(p0=0.5), dict()), + ("pop_percentile", dict(p0=0.1, p1=0.9), dict()), + ("subtract_mean_percentile", dict(p0=0.1, p1=0.9), dict()), + ("sum_percentile", dict(p0=0.1, p1=0.9), dict()), + ("threshold_percentile", dict(p0=0.5), dict()), + # bilateral.py + ("mean_bilateral", dict(s0=10, s1=10), dict()), + ("pop_bilateral", dict(s0=10, s1=10), dict()), + ("sum_bilateral", dict(s0=10, s1=10), dict()), +] + + +def _parse_shape(img_size): + return tuple(list(map(int, img_size.split(",")))) + + +def _parse_radii(radii): + return [int(r) for r in radii.split(",")] + + +def _parse_footprint_sizes(footprint_sizes): + return [int(size) for size in footprint_sizes.split(",")] + + +def _make_footprints(args): + if args.footprint_shape == "disk": + return [disk(radius).astype(bool) for radius in _parse_radii(args.radii)] + + footprint_sizes = _parse_footprint_sizes(args.footprint_sizes) + if any(size <= 1 or size % 2 == 0 for size in footprint_sizes): + raise ValueError("rank filter benchmark footprint sizes must be odd and > 1") + + return [np.ones((size, size), dtype=bool) for size in footprint_sizes] + + +def main(args): + cfile = "cucim_filters_rank_results.csv" + if getattr(args, "no_resume", False) or not os.path.exists(cfile): + all_results = pd.DataFrame() + else: + all_results = pd.read_csv(cfile, index_col=0) + dtypes = [np.dtype(args.dtype)] + + shape = _parse_shape(args.img_size) + if len(shape) != 2: + raise ValueError("rank filter benchmarks use 2D images") + + footprints = _make_footprints(args) + + for function_name, fixed_kwargs, var_kwargs in RANK_FILTERS: + if function_name != args.func_name: + continue + + var_kwargs = dict(var_kwargs) + var_kwargs["footprint"] = footprints + fixed_kwargs_gpu = dict(fixed_kwargs) + fixed_kwargs_gpu["backend"] = args.backend + + B = ImageBench( + function_name=function_name, + shape=shape, + dtypes=dtypes, + fixed_kwargs=fixed_kwargs, + fixed_kwargs_gpu=fixed_kwargs_gpu, + var_kwargs=var_kwargs, + module_cpu=skimage.filters.rank, + module_gpu=cucim.skimage.filters.rank, + run_cpu=not args.no_cpu, + ) + results = B.run_benchmark(duration=args.duration) + all_results = pd.concat([all_results, results["full"]]) + + fbase = os.path.splitext(cfile)[0] + all_results.to_csv(cfile, index=True) + with open(fbase + ".md", "w") as f: + f.write(all_results.to_markdown()) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Benchmarking cuCIM rank filters") + func_name_choices = [filter_spec[0] for filter_spec in RANK_FILTERS] + dtype_choices = [ + "float16", + "float32", + "float64", + "int8", + "int16", + "int32", + "int64", + "uint8", + "uint16", + "uint32", + "uint64", + ] + parser.add_argument( + "-i", "--img_size", type=str, help="Size of input image", required=True + ) + parser.add_argument( + "-d", + "--dtype", + type=str, + help="Dtype of input image", + choices=dtype_choices, + required=True, + ) + parser.add_argument( + "-f", + "--func_name", + type=str, + help="function to benchmark", + choices=func_name_choices, + required=True, + ) + parser.add_argument( + "-t", + "--duration", + type=int, + help="time to run benchmark", + required=True, + ) + parser.add_argument( + "--footprint_sizes", + type=str, + help=( + "Comma-separated odd square footprint side lengths to benchmark. " + "The all-ones rectangular footprints exercise the histogram fast " + "path for supported uint8 2D rank filters." + ), + default="3,7,15,31", + ) + parser.add_argument( + "--radii", + type=str, + help=( + "Comma-separated disk footprint radii to benchmark when " + "--footprint_shape=disk." + ), + default="1,3,7,15", + ) + parser.add_argument( + "--footprint_shape", + type=str, + choices=["rectangle", "disk"], + help=( + "Footprint family to benchmark. rectangle uses all-ones square " + "footprints and is the default." + ), + default="rectangle", + ) + parser.add_argument( + "--backend", + type=str, + choices=["auto", "histogram", "elementwise"], + help=( + "cuCIM rank backend to benchmark. auto uses automatic dispatch, " + "histogram requires the uint8 2D rectangular histogram backend, " + "and elementwise forces the generic per-output-pixel backend." + ), + default="auto", + ) + parser.add_argument( + "--no_cpu", + action="store_true", + help="disable cpu measurements", + default=False, + ) + parser.add_argument( + "--no_resume", + action="store_true", + help="do not load existing results CSV; save only this run's results (overwrite)", + default=False, + ) + + args = parser.parse_args() + main(args) diff --git a/benchmarks/skimage/run-nv-bench-filters-rank.sh b/benchmarks/skimage/run-nv-bench-filters-rank.sh new file mode 100755 index 000000000..d9c18c496 --- /dev/null +++ b/benchmarks/skimage/run-nv-bench-filters-rank.sh @@ -0,0 +1,23 @@ +#!/bin/bash +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +# Use env var if set/non-empty, otherwise default to 3 +MAX_DURATION="${CUCIM_BENCHMARK_MAX_DURATION:-3}" + +# Use env var if set/non-empty, otherwise default to "3,7,15,31" +FOOTPRINT_SIZES="${CUCIM_BENCHMARK_RANK_FOOTPRINT_SIZES:-3,7,15,31}" + +# Use env var if set/non-empty, otherwise default to automatic dispatch +BACKEND="${CUCIM_BENCHMARK_RANK_BACKEND:-auto}" + +param_shape=("512,512" "1920,1080") +param_filt=(autolevel enhance_contrast entropy equalize geometric_mean gradient majority maximum mean median minimum modal noise_filter pop subtract_mean sum threshold autolevel_percentile enhance_contrast_percentile gradient_percentile mean_percentile percentile pop_percentile subtract_mean_percentile sum_percentile threshold_percentile mean_bilateral pop_bilateral sum_bilateral) +param_dt=(uint8) +for shape in "${param_shape[@]}"; do + for filt in "${param_filt[@]}"; do + for dt in "${param_dt[@]}"; do + python cucim_filters_rank_bench.py -f "$filt" -i "$shape" -d "$dt" -t "$MAX_DURATION" --footprint_sizes "$FOOTPRINT_SIZES" --backend "$BACKEND" + done + done +done diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py index 06d16d4c2..64792b8ea 100644 --- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py +++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py @@ -169,6 +169,7 @@ def _call_kernel( weights, output, structure=None, + mask=None, weights_dtype=numpy.float64, structure_dtype=numpy.float64, ): @@ -200,6 +201,9 @@ def _call_kernel( if structure is not None: structure = cupy.ascontiguousarray(structure, structure_dtype) args.append(structure) + if mask is not None: + mask = cupy.ascontiguousarray(mask, bool) + args.append(mask) output = _util._get_output(output, input, None, complex_output) # noqa needs_temp = cupy.shares_memory(output, input, "MAY_SHARE_BOUNDS") if needs_temp: diff --git a/python/cucim/src/cucim/skimage/filters/__init__.pyi b/python/cucim/src/cucim/skimage/filters/__init__.pyi index 3bdf5051a..c1e110cf3 100644 --- a/python/cucim/src/cucim/skimage/filters/__init__.pyi +++ b/python/cucim/src/cucim/skimage/filters/__init__.pyi @@ -1,5 +1,5 @@ # SPDX-FileCopyrightText: 2009-2022 the scikit-image team -# SPDX-FileCopyrightText: Copyright (c) 2023-2025, NVIDIA CORPORATION. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION. All rights reserved. # SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause __all__ = [ @@ -24,6 +24,7 @@ __all__ = [ "prewitt", "prewitt_h", "prewitt_v", + "rank", "rank_order", "roberts", "roberts_neg_diag", @@ -52,6 +53,7 @@ __all__ = [ "window", ] +from . import rank from ._fft_based import butterworth from ._gabor import gabor, gabor_kernel from ._gaussian import difference_of_gaussians, gaussian diff --git a/python/cucim/src/cucim/skimage/filters/_median.py b/python/cucim/src/cucim/skimage/filters/_median.py index 31e4176de..9f3a541b7 100644 --- a/python/cucim/src/cucim/skimage/filters/_median.py +++ b/python/cucim/src/cucim/skimage/filters/_median.py @@ -1,5 +1,5 @@ # 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 import math @@ -34,9 +34,8 @@ def median( If ``None``, ``footprint`` will be a N-D array with 3 elements for each dimension (e.g., vector, square, cube, etc.). If `footprint` is a tuple of integers, it will be an array of ones with the given shape. - Otherwise, if ``behavior=='rank'``, ``footprint`` is a 2-D array of 1's - and 0's. If ``behavior=='ndimage'``, ``footprint`` is a N-D array of - 1's and 0's with the same number of dimension as ``image``. + Otherwise, ``footprint`` is an N-D array of 1's and 0's with the same + number of dimensions as ``image``. Note that upstream scikit-image currently does not support supplying a tuple for `footprint`. It is added here to avoid overhead of generating a small weights array in cases where it is not needed. @@ -44,26 +43,19 @@ def median( If None, a new array is allocated. mode : {'reflect', 'constant', 'nearest', 'mirror','‘wrap'}, optional The mode parameter determines how the array borders are handled, where - ``cval`` is the value when mode is equal to 'constant'. + ``cval`` is the value when mode is equal to 'constant'. ``mode`` is + only used when ``behavior='ndimage'``. Default is 'nearest'. - - .. versionadded:: 0.15 - ``mode`` is used when ``behavior='ndimage'``. cval : scalar, optional - Value to fill past edges of input if mode is 'constant'. Default is 0.0 - - .. versionadded:: 0.15 - ``cval`` was added in 0.15 is used when ``behavior='ndimage'``. + Value to fill past edges of input if mode is 'constant'. ``cval`` is + only used when ``behavior='ndimage'``. Default is 0.0. behavior : {'ndimage', 'rank'}, optional - Either to use the old behavior (i.e., < 0.15) or the new behavior. - The old behavior will call the :func:`skimage.filters.rank.median`. - The new behavior will call the :func:`scipy.ndimage.median_filter`. - Default is 'ndimage'. - - .. versionadded:: 0.15 - ``behavior`` is introduced in 0.15 - .. versionchanged:: 0.16 - Default ``behavior`` has been changed from 'rank' to 'ndimage' + Behavior 'ndimage' behaves like `cupyx.scipy.ndimage.median_filter`, + while 'rank' uses :func:`cucim.skimage.filters.rank.median`. cuCIM rank + filters use reflected boundary extension with a constant footprint + size and support N-D images. This differs from scikit-image rank + filters, which crop neighborhoods near image boundaries and support + only 2-D and 3-D images. Default is 'ndimage'. Other Parameters ---------------- @@ -81,8 +73,8 @@ def median( Returns ------- - out : 2-D array (same dtype as input image) - Output image. + out : N-D array + Output image with the same shape as the input image. See also -------- @@ -124,9 +116,15 @@ def median( "otherwise.", stacklevel=2, ) - raise NotImplementedError("rank behavior not currently implemented") - # TODO: implement median rank filter - # return generic.median(image, footprint=footprint, out=out) + from .rank import median as rank_median + + if isinstance(footprint, tuple): + if len(footprint) != image.ndim: + raise ValueError( + "tuple footprint must have ndim matching image" + ) + footprint = cp.ones(footprint, dtype=bool) + return rank_median(image, footprint=footprint, out=out) if footprint is None: footprint_shape = (3,) * image.ndim diff --git a/python/cucim/src/cucim/skimage/filters/rank/__init__.py b/python/cucim/src/cucim/skimage/filters/rank/__init__.py new file mode 100644 index 000000000..cf398a4a2 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/__init__.py @@ -0,0 +1,222 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +"""GPU-accelerated rank filters. + +This module provides GPU (CuPy/CUDA) implementations of most of the local +rank filters from ``skimage.filters.rank``, including generic, percentile, and +bilateral variants. The only unimplemented functions are: + +1. ``otsu`` (local Otsu thresholding via between-class variance maximization) +2. ``windowed_histogram`` (returns the full local histogram per pixel) + +These two operations did not map cleanly to the common patterns shared across +the other filters implemented here. + +Implementation approach +----------------------- + +The general elementwise backend computes output pixels independently: each +output pixel is assigned to a single GPU thread that gathers its local +neighborhood and applies the requested operation. Within this backend, +operations that do not require sorted values use streaming reductions, while +operations requiring rank ordering use a sorted-neighborhood kernel. + +For compatible 2-D uint8 inputs with fully populated, odd-sized rectangular +footprints (including square footprints), selected filters can instead use a +cooperative sliding-window histogram fast path. This backend partitions the +output rows among CUDA blocks, whose threads cooperatively maintain and +evaluate local histograms while traversing multiple output pixels. + +scikit-image, by contrast, always uses a sliding-window histogram approach that +incrementally updates a histogram as it moves across the image. Its current CPU +implementation relies on sequential neighborhood updates and is restricted to +2-D inputs (with 3-D support for a subset of filters). The elementwise backend +in cuCIM does not have this dimensionality restriction and all filters are +available in N-D, although the histogram-based GPU fast path is restricted to +2-D uint8 inputs. + +For small window sizes, the general elementwise approach is faster on the GPU, +but a histogram-based approach becomes much faster at large window sizes. The +default setting of ``backend='auto'`` uses thresholds derived from performance +measurements on an RTX A6000 to choose the best approach for a given window +size. This choice can be overridden by explicitly choosing +``backend='elementwise'`` or ``backend='histogram'``. + +Note that the elementwise approach is more general and supports features such +as arbitrary footprint shape and an image ``mask``, which are not available +for the histogram-based approach. + +Note that behavior at image boundaries differs from scikit-image. The GPU +implementations in cuCIM extend the boundary by reflection, so the footprint +does not shrink at image edges. A supplied mask can still reduce the number of +samples contributing to a neighborhood. The scikit-image implementation does +not extend the image and instead crops the footprint to remain within the image +edges. + +cuCIM vs scikit-image +--------------------- + +The table below summarizes known behavioral differences. Results are otherwise +expected to match. + +| Feature | scikit-image (CPU) | cuCIM (GPU) | +|--------------------------|-----------------------------------------------------------------|-------------| +| Dimensions | 2-D (3-D for generic filters) | N-dimensional | +| Supported dtypes | uint8 and uint16 natively; other real inputs are converted to uint8 | Integer and real floating-point dtypes; non-uint8 inputs are converted to uint8 by default | +| Output dtype | Usually the processed input dtype; entropy returns float64 | Usually the processed input dtype; entropy defaults to float32 | +| Algorithm | Sliding-window histogram | Streaming reductions, sorted neighborhoods, or uint8 2-D histogram fast path | +| Boundary handling | Excludes out-of-bounds pixels (population decreases at borders) | SciPy ``ndimage``-style reflected boundary extension; no boundary-induced population decrease | + +By default, non-uint8 inputs are converted to uint8 with +``img_as_ubyte`` before rank filtering. Set ``cast_to_uint8=False`` to opt +out of this conversion and use the elementwise kernels on the native dtype. +These kernels have the following behavioral differences. The cuCIM behavior is +generally preferable in these cases. + +| Feature | scikit-image (CPU) | cuCIM (GPU) | +|--------------------------|-----------------------------------------------------------------|-------------| +| ``mean_percentile`` | Can return zero when no histogram bin lies in the percentile interval | Selects samples by rank, avoiding this empty-bin artifact | +| ``subtract_mean_percentile`` | Can return zero when no histogram bin lies in the percentile interval | Selects samples by rank, avoiding this empty-bin artifact | +| ``sum`` | uint8/uint16 output can overflow | Preserves native input dtype; use a wider dtype to avoid overflow | +| ``sum_bilateral`` | uint8/uint16 output can overflow | Preserves native input dtype; use a wider dtype to avoid overflow | +| ``sum_percentile`` | uint8/uint16 output can overflow | Preserves native input dtype; use a wider dtype to avoid overflow | + +See the ``_percentile``, ``_generic``, and ``_bilateral`` modules for +additional per-function notes on dtype handling and behavioral differences. + +Histogram fast path +------------------- + +At larger window sizes, when the input is uint8 (or converted to uint8 via +``cast_to_uint8``, which is enabled by default), a histogram-based approach is +often beneficial. + +A uint8 2-D sliding-histogram backend is selected automatically for these rank +filters when all compatibility conditions below are met and the fully +populated rectangular footprint (including square footprints) is at least the +operation-specific benchmark-derived cutoff size. Sparse or arbitrarily shaped +footprints require the elementwise backend. + +* ``percentile`` +* ``median`` (implemented as ``percentile(p0=0.5)``) +* ``threshold_percentile`` +* ``mean_percentile`` with a non-full percentile range +* ``sum_percentile`` with a non-full percentile range +* ``pop_percentile`` with a non-full percentile range +* ``gradient_percentile`` with a non-full percentile range +* ``autolevel_percentile`` with a non-full percentile range +* ``enhance_contrast_percentile`` with a non-full percentile range +* ``subtract_mean_percentile`` with a non-full percentile range +* ``equalize`` +* ``geometric_mean`` +* ``mean_bilateral`` +* ``modal`` +* ``majority`` (alias for ``modal``) +* ``pop_bilateral`` +* ``sum_bilateral`` +* ``entropy`` + +The compatibility conditions are: + +* input image is 2-D and either has dtype ``uint8`` or is converted to + ``uint8`` before backend selection with ``cast_to_uint8=True``, the default +* footprint is a fully populated rectangular footprint with odd side lengths + greater than 1, for example ``cupy.ones((15, 15), dtype=bool)`` +* output dtype is one of uint8, uint16, float32, or float64; ``entropy`` + specifically requires float32 or float64 output +* no ``mask`` is provided +* no footprint shift is requested (``shift_x == shift_y == 0`` and no nonzero + ``shifts``) +* the internal boundary mode is ``reflect`` (the public rank wrappers use this + mode). This is the SciPy ``ndimage`` meaning of ``reflect``, which repeats + the edge value and is equivalent to ``numpy.pad(..., mode='symmetric')``: + ``d c b a | a b c d | d c b a``. The naming differs across libraries, so + this should not be confused with NumPy's ``mode='reflect'``. It also differs + from scikit-image rank filters, which do not extend the image and therefore + use smaller cropped neighborhoods near edges and corners. +* footprint half-width does not exceed the corresponding image extent + +With ``backend='auto'``, any unsupported case falls back to the generic +elementwise GPU implementation. For compatible calls, automatic dispatch also +requires a benchmark-derived minimum footprint area. Smaller footprints stay +on the generic per-output-pixel backend unless ``backend='histogram'`` is +requested explicitly. + +The automatic selection can be overridden with the keyword-only ``backend`` +parameter accepted by rank filters: + +* ``backend='auto'`` keeps the automatic selection behavior. +* ``backend='histogram'`` requires the histogram backend and raises + ``ValueError`` if the call is not compatible. +* ``backend='elementwise'`` forces the generic per-output-pixel backend. + +""" # noqa: E501 + +from ._generic import ( + autolevel, + enhance_contrast, + entropy, + equalize, + geometric_mean, + gradient, + majority, + maximum, + mean, + median, + minimum, + modal, + noise_filter, + pop, + subtract_mean, + sum, + threshold, +) +from ._percentile import ( + autolevel_percentile, + enhance_contrast_percentile, + gradient_percentile, + mean_percentile, + percentile, + pop_percentile, + subtract_mean_percentile, + sum_percentile, + threshold_percentile, +) +from ._bilateral import mean_bilateral, pop_bilateral, sum_bilateral + +__all__ = [ + 'autolevel', + 'autolevel_percentile', + 'enhance_contrast', + 'enhance_contrast_percentile', + 'entropy', + 'equalize', + 'geometric_mean', + 'gradient', + 'gradient_percentile', + 'majority', + 'maximum', + 'mean', + 'mean_bilateral', + 'mean_percentile', + 'median', + 'minimum', + 'modal', + 'noise_filter', + 'percentile', + 'pop', + 'pop_bilateral', + 'pop_percentile', + 'subtract_mean', + 'subtract_mean_percentile', + 'sum', + 'sum_bilateral', + 'sum_percentile', + 'threshold', + 'threshold_percentile', + # --- Not yet implemented --- + # 'otsu', + # 'windowed_histogram', +] diff --git a/python/cucim/src/cucim/skimage/filters/rank/_bilateral.py b/python/cucim/src/cucim/skimage/filters/rank/_bilateral.py new file mode 100644 index 000000000..be446b6c1 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/_bilateral.py @@ -0,0 +1,215 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +"""Bilateral rank filters for GPU (CuPy/CUDA). + +Bilateral filters include only neighborhood pixels whose grayvalue is within +a specified range around the center pixel: ``g - s1 < value < g + s0``, where +``g`` is the center pixel grayvalue and ``s0``, ``s1`` define the range. + +See ``cucim.skimage.filters.rank`` for a summary of differences between the +cuCIM and scikit-image implementations. +""" + +from ._percentile import ( + _apply, + _doc_boundary_note, + _doc_cast_to_uint8_param, + _doc_common_params, +) + +__all__ = [ + "mean_bilateral", + "pop_bilateral", + "sum_bilateral", +] + +# --- Docstring fragments --- + +_doc_s0_s1_params = """ + s0, s1 : int, optional + Define the bilateral range ``(g - s1, g + s0)`` around the center + pixel grayvalue ``g``. Only neighborhood pixels with values strictly + inside this interval are included. Default is 10.""" + +_doc_shifts_param = """ + shifts : sequence of int, optional (keyword-only) + N-dimensional offsets. If provided, shift_x and shift_y must be 0. + Length must match image.ndim.""" + +_doc_backend_param = """ + backend : {'auto', 'histogram', 'elementwise'}, optional (keyword-only) + Algorithm backend. ``'auto'`` selects the best compatible backend, + ``'histogram'`` requires a uint8 2-D image with a fully populated, + odd-sized rectangular footprint, and ``'elementwise'`` forces the + generic per-output-pixel backend. ``'histogram'`` raises ``ValueError`` + for an incompatible call.""" + +_doc_returns = """ + Returns + ------- + out : cupy.ndarray + Output image with the same shape as the input. The default output dtype + is uint8 for non-uint8 inputs converted with + ``cast_to_uint8=True``; otherwise it follows the input dtype unless + ``out`` controls it. +""" + + +def _build_bilateral_docstring(summary): + """Build a docstring for a bilateral rank filter.""" + return ( + summary + + "\n\n Parameters\n ----------" + + _doc_common_params + + _doc_s0_s1_params + + _doc_shifts_param + + _doc_backend_param + + _doc_cast_to_uint8_param + + "\n" + + _doc_returns + + _doc_boundary_note + ) + + +def mean_bilateral( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + s0=10, + s1=10, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "bilateral_mean", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=0, + p1=0, + shifts=shifts, + s0=s0, + s1=s1, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +mean_bilateral.__doc__ = _build_bilateral_docstring( + """Apply a flat kernel bilateral filter. + + This is an edge-preserving and noise reducing denoising filter. It averages + pixels based on their spatial closeness and radiometric similarity. + + Spatial closeness is measured by considering only the local pixel + neighborhood given by a footprint (structuring element). + + Radiometric similarity is defined by the graylevel interval ``(g-s1, g+s0)`` + where ``g`` is the current pixel graylevel. + + Only pixels belonging to the footprint and having a graylevel inside this + interval are averaged.""", +) + + +def pop_bilateral( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + s0=10, + s1=10, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "bilateral_pop", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=0, + p1=0, + shifts=shifts, + s0=s0, + s1=s1, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +pop_bilateral.__doc__ = _build_bilateral_docstring( + """Return the local number (population) of pixels in the bilateral range. + + The number of pixels is defined as the number of pixels which are included + in the footprint and the mask, and additionally have a graylevel inside + the interval ``(g-s1, g+s0)`` where ``g`` is the grayvalue of the center + pixel.""", +) + + +def sum_bilateral( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + s0=10, + s1=10, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "bilateral_sum", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=0, + p1=0, + shifts=shifts, + s0=s0, + s1=s1, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +sum_bilateral.__doc__ = _build_bilateral_docstring( + """Return the local sum of pixels in the bilateral range. + + Only pixels belonging to the footprint AND having a graylevel inside the + interval ``(g-s1, g+s0)`` are summed, where ``g`` is the current pixel + graylevel. + + The sum may overflow in a narrow output dtype. To accumulate into a wider + dtype, either provide a wider ``out`` array or promote the input and set + ``cast_to_uint8=False``. + + .. note:: + + scikit-image processes uint8 and uint16 inputs natively and returns + the sum in that dtype, so sufficiently large sums can overflow. cuCIM + can use a wider input or output dtype to avoid overflow.""", +) diff --git a/python/cucim/src/cucim/skimage/filters/rank/_generic.py b/python/cucim/src/cucim/skimage/filters/rank/_generic.py new file mode 100644 index 000000000..0559130cf --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/_generic.py @@ -0,0 +1,860 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +"""Generic rank filters for GPU. + +These are equivalent to the corresponding ``*_percentile`` functions called +with ``p0=0, p1=1`` (full range). This GPU implementation supports integer and +floating-point dtypes and N-D images; Boolean images are not supported. +scikit-image processes uint8 and uint16 natively and is restricted to 2-D/3-D +inputs. + +""" + +import cupy as cp +import numpy as np + +from ._percentile import ( + _apply, + _doc_boundary_note, + _doc_cast_to_uint8_param, + _doc_common_params, +) + +__all__ = [ + "autolevel", + "enhance_contrast", + "entropy", + "equalize", + "geometric_mean", + "gradient", + "majority", + "maximum", + "mean", + "median", + "minimum", + "modal", + "noise_filter", + "pop", + "subtract_mean", + "sum", + "threshold", +] + +# --- Docstring fragments for generic (no p0/p1) functions --- + +_doc_shifts_param_generic = """ + shift_z : int, optional + Additional footprint-center offset for 3-D images. For general N-D + offsets, use ``shifts`` instead. Default is 0. + shifts : sequence of int, optional (keyword-only) + N-dimensional offsets. If provided, shift_x, shift_y and shift_z + must be 0. Length must match image.ndim.""" + +_doc_backend_param = """ + backend : {'auto', 'histogram', 'elementwise'}, optional (keyword-only) + Algorithm backend. ``'auto'`` selects the best compatible backend, + ``'histogram'`` requires a supported operation on a uint8 2-D image + with a fully populated, odd-sized rectangular footprint, and + ``'elementwise'`` forces the generic per-output-pixel backend. + ``'histogram'`` raises ``ValueError`` for an incompatible call.""" + +_doc_returns = """ + Returns + ------- + out : cupy.ndarray + Output image with the same shape as the input. The default output dtype + is uint8 for non-uint8 inputs converted with + ``cast_to_uint8=True``; otherwise it follows the input dtype unless + ``out`` controls it. ``entropy`` defaults to float32 for non-floating + inputs when ``out`` is not provided. +""" + +_doc_common_params_median = _doc_common_params.replace( + """ footprint : cupy.ndarray + The neighborhood expressed as an array of 1's and 0's.""", + """ footprint : cupy.ndarray or None, optional + The neighborhood expressed as an array of 1's and 0's. If None, a + full footprint with shape ``(3,) * image.ndim`` is used.""", +) + + +def _build_generic_docstring(summary): + """Build a docstring for a generic (no p0/p1) rank filter.""" + return ( + summary + + "\n\n Parameters\n ----------" + + _doc_common_params + + _doc_shifts_param_generic + + _doc_backend_param + + _doc_cast_to_uint8_param + + "\n" + + _doc_returns + + _doc_boundary_note + ) + + +def _build_median_docstring(summary): + """Build a docstring for median, which has a default footprint.""" + return ( + summary + + "\n\n Parameters\n ----------" + + _doc_common_params_median + + _doc_shifts_param_generic + + _doc_backend_param + + _doc_cast_to_uint8_param + + "\n" + + _doc_returns + + _doc_boundary_note + ) + + +def _apply_generic( + operation, + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + p0=0, + p1=1, + backend="auto", + cast_to_uint8=True, + out_dtype=None, +): + """Apply a generic rank filter (defaults to full range p0=0, p1=1).""" + if not isinstance(image, cp.ndarray): + raise ValueError("image must be a CuPy array") + + # Convert shift_z into the N-D shifts parameter + if shifts is not None: + if shift_x != 0 or shift_y != 0 or shift_z != 0: + raise ValueError( + "shift_x, shift_y and shift_z must be 0 when shifts " + "is specified" + ) + elif shift_z != 0: + if image.ndim < 3: + raise ValueError( + "shift_z is only valid for 3D or higher dimensional images" + ) + shifts_list = [0] * image.ndim + shifts_list[0] = shift_z + shifts_list[1] = shift_y + shifts_list[2] = shift_x + shifts = tuple(shifts_list) + + return _apply( + operation, + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x if shifts is None else 0, + shift_y=shift_y if shifts is None else 0, + p0=p0, + p1=p1, + shifts=shifts, + out_dtype=out_dtype, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +def autolevel( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "autolevel", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +autolevel.__doc__ = _build_generic_docstring( + """Auto-level image using local histogram. + + This filter locally stretches the histogram of gray values to cover the + entire range of values from "white" to "black".""", +) + + +def gradient( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "gradient", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +gradient.__doc__ = _build_generic_docstring( + """Return local gradient of an image (i.e. local maximum - local + minimum).""", +) + + +def mean( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "mean", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +mean.__doc__ = _build_generic_docstring( + """Return local mean of an image.""", +) + + +def subtract_mean( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + result = _apply_generic( + "subtract_mean", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + # The generic version uses an offset of mid_bin - 1 (127 for uint8), + # while the percentile version uses mid_bin (128 for uint8). Adjust by + # subtracting 1 for integer dtypes to match scikit-image's generic + # subtract_mean. + if np.issubdtype(result.dtype, np.integer): + result -= 1 + return result + + +subtract_mean.__doc__ = _build_generic_docstring( + """Return image subtracted from its local mean. + + The output is:: + + out = (g - mean) * 0.5 + mid_bin - 1 + + where ``mean`` is the local neighborhood mean, ``g`` is the center pixel + value, and ``mid_bin`` is ``(dtype_max + 1) / 2`` (128 for uint8), so the + effective offset is 127 for uint8. + + .. note:: + + This function uses an output offset of ``(dtype_max + 1) / 2 - 1`` + (127 for uint8), matching scikit-image's ``subtract_mean``. The + percentile variant ``subtract_mean_percentile`` uses an offset of + ``(dtype_max + 1) / 2`` (128 for uint8), matching scikit-image's + ``subtract_mean_percentile``.""", +) + + +def enhance_contrast( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "enhance_contrast", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +enhance_contrast.__doc__ = _build_generic_docstring( + """Enhance contrast of an image. + + This replaces each pixel by the local maximum if the pixel gray value is + closer to the local maximum than the local minimum. Otherwise it is + replaced by the local minimum.""", +) + + +def pop( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "pop", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +pop.__doc__ = _build_generic_docstring( + """Return the local number (population) of pixels. + + The number of pixels is defined as the number of pixels which are included + in the footprint and the mask. + + .. note:: + + The output is constant across the entire image (equal to the number of + active footprint elements), except when a mask is provided. Unlike + scikit-image, the GPU implementation does not reduce the count at image + borders because the underlying kernel uses reflected boundary + extension. In scikit-image, the population decreases at borders + because the sliding window excludes out-of-bounds pixels.""", +) + + +def sum( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "sum", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +sum.__doc__ = _build_generic_docstring( + """Return the local sum of pixels. + + The sum may overflow in a narrow output dtype. To accumulate into a wider + dtype, either provide a wider ``out`` array or promote the input and set + ``cast_to_uint8=False``. + + .. note:: + + scikit-image processes uint8 and uint16 inputs natively and returns + the sum in that dtype, so sufficiently large sums can overflow. cuCIM + can use a wider input or output dtype to avoid overflow.""", +) + + +def minimum( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "minimum", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +minimum.__doc__ = _build_generic_docstring( + """Return the local minimum of an image. + + .. note:: + + This uses a streaming reduction over the neighborhood. If mask support + is not needed, ``cupyx.scipy.ndimage.minimum_filter`` may be faster.""", +) + + +def maximum( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "maximum", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +maximum.__doc__ = _build_generic_docstring( + """Return the local maximum of an image. + + .. note:: + + This uses a streaming reduction over the neighborhood. If mask support + is not needed, ``cupyx.scipy.ndimage.maximum_filter`` may be faster.""", +) + + +def median( + image, + footprint=None, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + if footprint is None and isinstance(image, cp.ndarray): + footprint = cp.ones((3,) * image.ndim, dtype=bool) + return _apply_generic( + "percentile", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + p0=0.5, + cast_to_uint8=cast_to_uint8, + ) + + +median.__doc__ = _build_median_docstring( + """Return the local median of an image. + + .. note:: + + This is implemented via ``percentile(p0=0.5)`` to ensure consistent + neighborhood-level mask handling with other rank filters. If mask + support is not needed, ``cupyx.scipy.ndimage.median_filter`` may + be faster.""", +) + + +def threshold( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "threshold_mean", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +threshold.__doc__ = _build_generic_docstring( + """Local threshold of an image. + + The output is 1 if the grayvalue of the center pixel is greater than the + local mean and 0 otherwise:: + + out = 1 if g > mean else 0 + + where ``g`` is the center pixel value and ``mean`` is the mean of all + neighborhood values. + + .. note:: + + This differs from ``threshold_percentile``, which compares to the + p0-th percentile value and outputs 0 or ``dtype_max`` (e.g. 0/255 + for uint8). The generic ``threshold`` compares to the local **mean** + and outputs 0 or 1.""", +) + + +def equalize( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "equalize", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +equalize.__doc__ = _build_generic_docstring( + """Equalize image using local histogram. + + The output is the rank of the center pixel within its local neighborhood, + scaled to the full output range:: + + out = dtype_max * rank(g) / N + + where ``rank(g)`` is the number of neighborhood values <= ``g``, ``N`` + is the neighborhood population, and ``dtype_max`` is the maximum value + for the output dtype (255 for uint8, 1.0 for float).""", +) + + +def geometric_mean( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "geometric_mean", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +geometric_mean.__doc__ = _build_generic_docstring( + """Return the local geometric mean of an image. + + The output is:: + + out = round(exp(mean(log(values + 1))) - 1) + + The ``+1`` / ``-1`` offset ensures that zero-valued pixels are handled + correctly (``log(0)`` is undefined).""", +) + + +def noise_filter( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "noise_filter", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +noise_filter.__doc__ = _build_generic_docstring( + """Noise feature filter. + + Returns 0 if the center pixel value appears among its neighbors (i.e. it + is not isolated noise). Otherwise returns the minimum absolute distance + to the nearest neighbor value. Higher values indicate more isolated + pixels.""", +) + + +def modal( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "modal", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +modal.__doc__ = _build_generic_docstring( + """Return the local modal (most frequent) value of an image.""", +) + + +def majority( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply_generic( + "modal", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +majority.__doc__ = _build_generic_docstring( + """Return the local majority value of an image. + + This is an alias for ``modal`` — it returns the most frequent value + in the local neighborhood.""", +) + + +def entropy( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + shift_z=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + out_dtype = None + if ( + out is None + and isinstance(image, cp.ndarray) + and (image.dtype.kind != "f" or cast_to_uint8) + ): + out_dtype = cp.float32 + return _apply_generic( + "entropy", + image, + footprint, + out, + mask, + shift_x, + shift_y, + shift_z, + shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + out_dtype=out_dtype, + ) + + +entropy.__doc__ = _build_generic_docstring( + """Return the local Shannon entropy of an image. + + The output is the entropy in bits of the local grayvalue distribution:: + + out = -sum(p * log2(p)) + + where the sum is over unique values in the neighborhood and + ``p = count / N`` is the probability of each value. + + .. note:: + + The output is a floating-point quantity (entropy in bits). When ``out`` + is not provided, integer inputs produce a floating-point output. + Explicit integer ``out`` arrays are respected and will truncate + fractional entropy values.""", +) diff --git a/python/cucim/src/cucim/skimage/filters/rank/_histogram.py b/python/cucim/src/cucim/skimage/filters/rank/_histogram.py new file mode 100644 index 000000000..dc6700d89 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/_histogram.py @@ -0,0 +1,275 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 + +import os + +import cupy as cp + +from ..._shared.utils import _to_np_mode +from ..._vendored import pad + +_HISTOGRAM_OPS = { + "percentile": 0, + "threshold": 1, + "mean": 2, + "sum": 3, + "pop": 4, + "gradient": 5, + "autolevel": 6, + "entropy": 7, + "enhance_contrast": 8, + "subtract_mean": 9, + "equalize": 10, + "bilateral_mean": 11, + "bilateral_pop": 12, + "bilateral_sum": 13, + "modal": 14, + "geometric_mean": 15, +} + +# Performance-tuning thresholds derived from benchmarks on an RTX A6000. +_HISTOGRAM_MIN_FOOTPRINT_AREA = { + "sum": 15 * 15, + "enhance_contrast": 17 * 17, + "gradient": 17 * 17, + "percentile": 17 * 17, + "pop": 17 * 17, + "threshold": 17 * 17, + "mean": 19 * 19, + "subtract_mean": 19 * 19, + "autolevel": 21 * 21, + "entropy": 25 * 25, + "bilateral_sum": 27 * 27, + "bilateral_pop": 29 * 29, + "bilateral_mean": 33 * 33, + "modal": 15 * 15, + "geometric_mean": 15 * 15, + "equalize": 91 * 91, +} + +_DEFAULT_SCRATCH_MB = 256 +_DEFAULT_MAX_PARTITIONS = 256 +_INT16_MAX = 32767 +_HISTOGRAM_COUNTER_TYPES = { + "int16": (cp.int16, "short"), + "int32": (cp.int32, "int"), +} +_HISTOGRAM_OUTPUT_TYPES = { + "float32": (cp.float32, "float"), + "float64": (cp.float64, "double"), + "uint8": (cp.uint8, "unsigned char"), + "uint16": (cp.uint16, "unsigned short"), +} + + +def _can_use_rank_histogram( + image, + footprint_shape, + output, + mask, + modes, + origins, + *, + has_weights, + operation, + p0, + p1, +): + """Return True for the restricted uint8 2-D histogram backend. + + This backend is intentionally narrow. It is selected only for supported + rank operations on 2-D uint8 images with an all-ones odd rectangular + footprint, no mask, reflect mode and zero origin. Unsupported cases fall + back to the generic rank implementation. + """ + if operation not in _HISTOGRAM_OPS: + return False + if ( + operation + in { + "autolevel", + "enhance_contrast", + "mean", + "subtract_mean", + "sum", + "pop", + "gradient", + } + and p0 <= 0 + and p1 >= 100 + ): + return False + if image.ndim != 2 or image.dtype != cp.uint8: + return False + if output is not None and output.dtype.name not in _HISTOGRAM_OUTPUT_TYPES: + return False + if ( + operation == "entropy" + and output is not None + and ( + output.dtype.kind != "f" + or output.dtype.name not in _HISTOGRAM_OUTPUT_TYPES + ) + ): + return False + if operation == "entropy" and output is None: + return False + if mask is not None or has_weights: + return False + if tuple(modes) != ("reflect", "reflect"): + return False + if any(origin != 0 for origin in origins): + return False + if len(footprint_shape) != 2: + return False + if any(size <= 1 or size % 2 == 0 for size in footprint_shape): + return False + radii = tuple(size // 2 for size in footprint_shape) + if any(radius > size for radius, size in zip(radii, image.shape)): + return False + return True + + +def _should_use_rank_histogram(operation, footprint_shape): + """Return True when benchmarks favor histogram over elementwise.""" + min_area = _HISTOGRAM_MIN_FOOTPRINT_AREA.get(operation) + if min_area is None: + return False + return footprint_shape[0] * footprint_shape[1] >= min_area + + +def _get_env_int(name, default): + value = os.environ.get(name) + if value is None or value == "": + return default + value = int(value) + if value <= 0: + raise ValueError(f"{name} must be a positive integer") + return value + + +def _get_histogram_counter_dtype(footprint_shape): + """Return the narrowest safe column-histogram counter dtype.""" + footprint_area = footprint_shape[0] * footprint_shape[1] + if footprint_area <= _INT16_MAX: + return cp.int16 + return cp.int32 + + +def _get_rank_histogram_partitions( + out_rows, cols, partitions=None, *, counter_dtype=cp.int32 +): + """Choose row partitions for the sliding-histogram backend. + + More partitions expose more row-band parallelism, but scratch memory grows + linearly as ``partitions * cols * 256 * sizeof(counter_dtype)``. + """ + if partitions is not None: + return min(max(1, int(partitions)), out_rows) + + partitions = os.environ.get("CUCIM_RANK_HISTOGRAM_PARTITIONS") + if partitions not in (None, ""): + return min(max(1, int(partitions)), out_rows) + + scratch_mb = _get_env_int( + "CUCIM_RANK_HISTOGRAM_SCRATCH_MB", _DEFAULT_SCRATCH_MB + ) + max_partitions = _get_env_int( + "CUCIM_RANK_HISTOGRAM_MAX_PARTITIONS", _DEFAULT_MAX_PARTITIONS + ) + bytes_per_partition = cols * 256 * cp.dtype(counter_dtype).itemsize + partitions_by_memory = max(1, (scratch_mb << 20) // bytes_per_partition) + + return min(max(1, out_rows // 2), max_partitions, partitions_by_memory) + + +@cp.memoize(for_each_device=True) +def _get_histogram_rank_kernel( + operation, counter_dtype_name, output_dtype_name +): + kernel_directory = os.path.join(os.path.dirname(__file__), "cuda") + with open(os.path.join(kernel_directory, "histogram_rank.cu")) as f: + code = "\n".join(f.readlines()) + + _, counter_type = _HISTOGRAM_COUNTER_TYPES[counter_dtype_name] + _, output_type = _HISTOGRAM_OUTPUT_TYPES[output_dtype_name] + code = ( + f"#define RANK_HIST_OP {_HISTOGRAM_OPS[operation]}\n" + f"#define HIST_COUNTER_T {counter_type}\n" + f"#define RANK_HIST_OUTPUT_T {output_type}\n" + code + ) + return cp.RawKernel(code=code, name="cuRankHistogram2DUint8") + + +def _rank_histogram( + image, + footprint_shape, + operation, + *, + output=None, + mode="reflect", + cval=0, + p0=0, + p1=100, + s0=0, + s1=0, + dtype_max=255, + partitions=None, +): + """Apply a uint8 2-D rectangular rank filter using a sliding histogram.""" + image = cp.ascontiguousarray(image) + radii = tuple(size // 2 for size in footprint_shape) + npad = tuple((radius, radius) for radius in radii) + np_mode = _to_np_mode(mode) + if np_mode == "constant": + pad_kwargs = dict(mode=np_mode, constant_values=cval) + else: + pad_kwargs = dict(mode=np_mode) + padded = pad(image, npad, **pad_kwargs) + + out_dtype = output.dtype if output is not None else padded.dtype + if cp.dtype(out_dtype).name not in _HISTOGRAM_OUTPUT_TYPES: + raise ValueError(f"unsupported histogram output dtype: {out_dtype}") + out = cp.empty(padded.shape, dtype=out_dtype) + rows, cols = padded.shape + out_rows = image.shape[0] + counter_dtype = _get_histogram_counter_dtype(footprint_shape) + counter_dtype_name = cp.dtype(counter_dtype).name + output_dtype_name = cp.dtype(out.dtype).name + partitions = _get_rank_histogram_partitions( + out_rows, cols, partitions=partitions, counter_dtype=counter_dtype + ) + + hist = cp.zeros((partitions * cols * 256,), dtype=counter_dtype) + op_code = _HISTOGRAM_OPS[operation] + kernel = _get_histogram_rank_kernel( + operation, counter_dtype_name, output_dtype_name + ) + window_size = footprint_shape[0] * footprint_shape[1] + kernel( + (partitions,), + (256,), + ( + padded, + out, + hist, + radii[0], + radii[1], + float(p0), + float(p1), + float(s0), + float(s1), + float(dtype_max), + op_code, + window_size, + rows, + cols, + ), + ) + + out_sl = tuple(slice(radius, -radius) for radius in radii) + result = out[out_sl] + if output is not None: + output[...] = result + return output + return result diff --git a/python/cucim/src/cucim/skimage/filters/rank/_percentile.py b/python/cucim/src/cucim/skimage/filters/rank/_percentile.py new file mode 100644 index 000000000..99bda5761 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/_percentile.py @@ -0,0 +1,762 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +"""Percentile rank filters for GPU (CuPy/CUDA). + +Inferior and superior ranks, provided by the user, are passed to the kernel +function to provide a softer version of the rank filters. E.g. +``autolevel_percentile`` will stretch image levels between percentile [p0, p1] +instead of using [min, max]. It means that isolated bright or dark pixels will +not produce halos. + +See ``cucim.skimage.filters.rank`` for a summary of differences between the +cuCIM and scikit-image implementations. + +Dtype notes +----------- + +Some operations use a ``dtype_max`` value based on the output dtype that +affects output scaling (``autolevel_percentile``, ``threshold_percentile``, +``subtract_mean_percentile``). With the default ``cast_to_uint8=True``, +non-uint8 inputs are converted before filtering and the default output dtype +is uint8. When processing other dtypes natively with ``cast_to_uint8=False``: + +- **Unsigned integers** (uint8, uint16, etc.): ``dtype_max`` is the type + maximum (255, 65535, ...). This matches scikit-image's behavior. +- **Float** (float32, float64): ``dtype_max = 1.0``, assuming normalized + images in [0, 1]. Operations like ``autolevel_percentile`` will scale + output to [0.0, 1.0] and ``threshold_percentile`` will output 0.0 or 1.0. +- **Signed integers** (int8, int16, etc.): ``dtype_max`` uses the positive + maximum (127 for int8, 32767 for int16). This means operations like + ``autolevel_percentile`` scale to [0, 127] (positive half only) and + ``subtract_mean_percentile`` centers at 64 (not 0). Signed integer inputs + are accepted but may not give intuitive results for these operations. + +""" + +import warnings + +import cupy as cp + +from ...util import img_as_ubyte +from ._rank_filter import _skimage_rank_filter + +__all__ = [ + "autolevel_percentile", + "enhance_contrast_percentile", + "gradient_percentile", + "mean_percentile", + "percentile", + "pop_percentile", + "subtract_mean_percentile", + "sum_percentile", + "threshold_percentile", +] + +_ZERO_FOR_EMPTY_FOOTPRINT_OPS = { + "geometric_mean", + "maximum", + "mean", + "minimum", +} + +# --- Common docstring fragments --- + +_doc_common_params = """ + image : cupy.ndarray + N-D input image with an integer or floating-point dtype. Boolean + images are not supported. + footprint : cupy.ndarray + The neighborhood expressed as an array of 1's and 0's. + out : cupy.ndarray, optional + If None, a new array is allocated. + mask : cupy.ndarray, optional + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int, optional + Footprint-center offsets along axes 1 and 0, respectively. For general + N-D offsets, use ``shifts`` instead. Default is 0.""" + +_doc_p0_p1_params = """ + p0, p1 : float, optional, in interval [0, 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. Defaults are 0 and 1, respectively.""" + +_doc_p0_only_param = """ + p0 : float, optional, in interval [0, 1] + Set the percentile value. Default is 0.""" + +_doc_shifts_param = """ + shifts : sequence of int, optional (keyword-only) + N-dimensional offsets. If provided, shift_x and shift_y must be 0. + Length must match image.ndim.""" + +_doc_backend_param = """ + backend : {'auto', 'histogram', 'elementwise'}, optional (keyword-only) + Algorithm backend. ``'auto'`` selects the best compatible backend, + ``'histogram'`` requires a supported operation on a uint8 2-D image + with a fully populated, odd-sized rectangular footprint, and + ``'elementwise'`` forces the generic per-output-pixel backend. + ``'histogram'`` raises ``ValueError`` for an incompatible call.""" + +_doc_boundary_note = """ + + Notes + ----- + Rank filters use reflected boundary extension. The name ``reflect`` can be + confusing because padding libraries use different conventions. Here it + follows the SciPy ``ndimage`` convention, which repeats the edge value + (equivalent to ``numpy.pad(..., mode='symmetric')``):: + + d c b a | a b c d | d c b a + + This also differs from scikit-image's rank filters, which do not extend + the image at the boundary. scikit-image uses cropped neighborhoods at + edges and corners, so the effective footprint population can be smaller + near the image border. A supplied mask can still reduce the number of + samples contributing to a cuCIM neighborhood.""" + +_doc_cast_to_uint8_param = """ + cast_to_uint8 : bool, optional (keyword-only) + If True, non-uint8 image inputs are converted to uint8 with + ``img_as_ubyte`` before backend selection. This matches scikit-image's + conversion behavior for many other input dtypes and can enable the + uint8 histogram backend for compatible inputs. Set this to False to + process a wider input dtype natively. Default is True.""" + +_doc_returns = """ + Returns + ------- + out : cupy.ndarray + Output image with same shape as input. The default output dtype is + uint8 for non-uint8 inputs converted with ``cast_to_uint8=True``; + otherwise it follows the input dtype unless ``out`` controls it. +""" + + +def _build_docstring(summary, *, p0_only=False): + """Build a docstring from common fragments.""" + pct_params = _doc_p0_only_param if p0_only else _doc_p0_p1_params + return ( + summary + + "\n\n Parameters\n ----------" + + _doc_common_params + + pct_params + + _doc_shifts_param + + _doc_backend_param + + _doc_cast_to_uint8_param + + "\n" + + _doc_returns + + _doc_boundary_note + ) + + +def _preprocess_input( + image, + footprint=None, + out=None, + mask=None, + out_dtype=None, + shifts=None, + cast_to_uint8=True, +): + """Preprocess and verify input for filters.rank methods (GPU version).""" + if not isinstance(image, cp.ndarray): + raise ValueError("image must be a CuPy array") + + if isinstance(out, cp.ndarray) and cp.shares_memory( + image, out, "MAY_SHARE_BOUNDS" + ): + raise NotImplementedError("Cannot perform rank operation in place.") + + input_dtype = image.dtype + if input_dtype == bool or out_dtype == bool: + raise ValueError("dtype cannot be bool.") + + if cast_to_uint8 and image.dtype != cp.dtype(cp.uint8): + if image.dtype.kind == "f": + warnings.warn( + f"Possible precision loss converting image of type " + f"{image.dtype} to uint8 as required by rank filters. " + f"Convert manually using cucim.skimage.util.img_as_ubyte to " + f"silence this warning.", + stacklevel=3, + ) + image = img_as_ubyte(image) + else: + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message=( + r"Downcasting .* to uint8 without scaling because " + r"max value .* fits in uint8" + ), + category=UserWarning, + ) + image = img_as_ubyte(image) + input_dtype = image.dtype + + # Convert footprint to boolean CuPy array + if footprint is not None: + if not isinstance(footprint, cp.ndarray): + raise ValueError("footprint must be a CuPy array") + footprint = cp.ascontiguousarray(footprint > 0, dtype=bool) + if footprint.ndim != image.ndim: + raise ValueError( + "Image dimensions and footprint dimensions do not match" + ) + + # Ensure image is contiguous + if not image.flags.c_contiguous: + image = cp.ascontiguousarray(image) + + # Handle mask + if mask is not None: + if not isinstance(mask, cp.ndarray): + raise ValueError("mask must be a CuPy array") + mask = cp.ascontiguousarray(mask > 0, dtype=bool) + if mask.shape != image.shape: + raise ValueError("Mask shape must match image shape") + + # Handle output array + if out is None: + if out_dtype is None: + out_dtype = image.dtype + out = cp.empty(image.shape, dtype=out_dtype) + else: + if not isinstance(out, cp.ndarray): + raise ValueError("out must be a CuPy array") + if out.shape != image.shape: + raise ValueError("out shape must match image shape") + + # Handle shifts parameter + origin = 0 # Default origin + if shifts is not None: + if not hasattr(shifts, "__len__"): + raise ValueError("shifts must be a sequence") + if len(shifts) != image.ndim: + raise ValueError( + f"shifts length ({len(shifts)}) must match image.ndim " + f"({image.ndim})" + ) + # Convert shifts to origin (shifts are offsets from center) + # Note: In ndimage, origin shifts the filter in the opposite + # direction. For now, map shifts directly to origin + if any(s != 0 for s in shifts): + # origin = tuple(-s for s in shifts) # Negate for opposite + origin = tuple(shifts) # Or use directly + # TODO: Verify the sign convention matches scikit-image + + return image, footprint, out, mask, origin + + +def _apply( + operation, + image, + footprint, + out, + mask, + shift_x, + shift_y, + p0, + p1, + out_dtype=None, + shifts=None, + s0=0, + s1=0, + backend="auto", + cast_to_uint8=True, +): + """Apply percentile range filter with specified operation.""" + if not isinstance(image, cp.ndarray): + raise ValueError("image must be a CuPy array") + + # Handle shift_x, shift_y vs shifts + if shifts is not None: + if shift_x != 0 or shift_y != 0: + raise ValueError( + "shift_x and shift_y must be 0 when shifts is specified" + ) + else: + # Convert shift_x, shift_y to shifts for 2D compatibility + if image.ndim >= 2: + # For 2D+: shift_y applies to axis 0, shift_x to axis 1 + shifts = [0] * image.ndim + shifts[0] = shift_y + shifts[1] = shift_x + shifts = tuple(shifts) + elif shift_x != 0 or shift_y != 0: + raise ValueError( + "shift_x and shift_y are only valid for 2D or higher " + "dimensional images" + ) + + image, footprint, out, mask, origin = _preprocess_input( + image, + footprint, + out, + mask, + out_dtype, + shifts=shifts, + cast_to_uint8=cast_to_uint8, + ) + + if ( + operation in _ZERO_FOR_EMPTY_FOOTPRINT_OPS + and footprint is not None + and not bool(cp.any(footprint)) + ): + out.fill(0) + return out + + # Convert percentiles from [0, 1] to [0, 100] for our implementation + p0_pct = p0 * 100.0 + p1_pct = p1 * 100.0 + + # Call the GPU implementation + result = _skimage_rank_filter( + image, + p0=p0_pct, + p1=p1_pct, + operation=operation, + footprint=footprint, + output=out, + mode="reflect", + cval=0.0, + origin=origin, + axes=None, + mask=mask, + s0=s0, + s1=s1, + backend=backend, + ) + + return result + + +def autolevel_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "autolevel", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +autolevel_percentile.__doc__ = _build_docstring( + """Return grayscale local autolevel of an image. + + This filter locally stretches the histogram of grayvalues to cover the + entire range of values from "white" to "black". + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is:: + + out = dtype_max * (clamp(g, v_p0, v_p1) - v_p0) / (v_p1 - v_p0) + + where ``v_p0`` and ``v_p1`` are the local values at percentiles p0 and + p1, ``g`` is the center pixel value, and ``dtype_max`` is the maximum + value for the output dtype (255 for uint8, 65535 for uint16, 1.0 for + float). See the module-level dtype notes for signed integer behavior.""", +) + + +def gradient_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "gradient", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +gradient_percentile.__doc__ = _build_docstring( + """Return local gradient of an image (i.e. local maximum - local minimum). + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is:: + + out = v_p1 - v_p0 + + where ``v_p0`` and ``v_p1`` are the local values at percentiles p0 and + p1.""", +) + + +def mean_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "mean", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +mean_percentile.__doc__ = _build_docstring( + """Return local mean of an image. + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is the arithmetic mean of all neighborhood values + whose sorted position falls within the [p0, p1] percentile range. + + .. note:: + + scikit-image's histogram-based implementation can produce spurious + zero outputs in low-variance neighborhoods where no histogram bin + falls entirely within the percentile window. Both cuCIM backends + ensure that a nonempty neighborhood contributes at least one value to + the selected percentile range, avoiding such artifacts.""", +) + + +def subtract_mean_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "subtract_mean", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +subtract_mean_percentile.__doc__ = _build_docstring( + """Return image subtracted from its local mean. + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is:: + + out = (g - mean_p) * 0.5 + mid_bin + + where ``mean_p`` is the mean of neighborhood values in the [p0, p1] + percentile range, ``g`` is the center pixel value, and ``mid_bin`` is + ``(dtype_max + 1) / 2`` (128 for uint8). + + .. note:: + + scikit-image's histogram-based implementation can produce spurious + zero outputs in low-variance neighborhoods where no histogram bin + falls entirely within the percentile window. Both cuCIM backends + ensure that a nonempty neighborhood contributes at least one value to + the selected percentile range, avoiding such artifacts. + + .. note:: + + This function uses an output offset of ``(dtype_max + 1) / 2`` (128 + for uint8), matching scikit-image's ``subtract_mean_percentile``. The + non-percentile ``subtract_mean`` in ``filters.rank`` uses an offset of + ``(dtype_max + 1) / 2 - 1`` (127 for uint8), matching scikit-image's + ``subtract_mean``.""", +) + + +def enhance_contrast_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "enhance_contrast", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +enhance_contrast_percentile.__doc__ = _build_docstring( + """Enhance contrast of an image. + + This replaces each pixel by the local maximum if the pixel grayvalue is + closer to the local maximum than the local minimum. Otherwise it is + replaced by the local minimum. + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is:: + + out = v_p1 if (v_p1 - g) < (g - v_p0) else v_p0 + + where ``v_p0`` and ``v_p1`` are the local values at percentiles p0 and + p1, and ``g`` is the center pixel value.""", +) + + +def percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "percentile", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p0, # p1 not used for single percentile + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +percentile.__doc__ = _build_docstring( + """Return local percentile of an image. + + Returns the value of the p0 lower percentile of the local grayvalue + distribution. The output is the value at position + ``floor(p0 * N)`` in the sorted neighborhood, where N is the + neighborhood population.""", + p0_only=True, +) + + +def pop_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "pop", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +pop_percentile.__doc__ = _build_docstring( + """Return the local number (population) of pixels. + + The number of pixels is defined as the number of pixels which are included + in the footprint and the mask. + + Only grayvalues between percentiles [p0, p1] are considered in the + filter. The output is the count of neighborhood pixels whose values fall + in histogram bins where the cumulative count is in [p0 * N, p1 * N], + where N is the neighborhood population.""", +) + + +def sum_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + p1=1, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "sum", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p1, + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +sum_percentile.__doc__ = _build_docstring( + """Return the local sum of pixels. + + Only grayvalues between percentiles [p0, p1] are considered in the filter. + + The sum may overflow in a narrow output dtype. To accumulate into a wider + dtype, either provide a wider ``out`` array or promote the input and set + ``cast_to_uint8=False``. + + .. note:: + + scikit-image processes uint8 and uint16 inputs natively and returns + the sum in that dtype, so sufficiently large sums can overflow. cuCIM + can use a wider input or output dtype to avoid overflow.""", +) + + +def threshold_percentile( + image, + footprint, + out=None, + mask=None, + shift_x=0, + shift_y=0, + p0=0, + *, + shifts=None, + backend="auto", + cast_to_uint8=True, +): + return _apply( + "threshold", + image, + footprint, + out=out, + mask=mask, + shift_x=shift_x, + shift_y=shift_y, + p0=p0, + p1=p0, # p1 not used for threshold + shifts=shifts, + backend=backend, + cast_to_uint8=cast_to_uint8, + ) + + +threshold_percentile.__doc__ = _build_docstring( + """Local threshold of an image. + + The output is ``dtype_max`` if the grayvalue of the center pixel is greater + than or equal to the value at the p0 percentile, and 0 otherwise:: + + out = dtype_max if g >= v_p0 else 0 + + where ``v_p0`` is the local value at percentile p0, ``g`` is the center + pixel value, and ``dtype_max`` is the maximum value for the output dtype + (255 for uint8, 65535 for uint16, 1.0 for float). See the module-level + dtype notes for signed integer behavior. + + .. note:: + + This differs from generic ``threshold``, which compares to the local + **mean** and outputs 0 or 1 (not 0 or ``dtype_max``). + ``threshold_percentile`` compares to the **p0-th percentile** value + and outputs 0 or ``dtype_max``.""", + p0_only=True, +) diff --git a/python/cucim/src/cucim/skimage/filters/rank/_rank_filter.py b/python/cucim/src/cucim/skimage/filters/rank/_rank_filter.py new file mode 100644 index 000000000..97dea7438 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/_rank_filter.py @@ -0,0 +1,1133 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +"""Shared implementation machinery for GPU rank filters. + +This internal module normalizes neighborhood arguments, generates elementwise +CUDA kernels, and dispatches compatible operations to the sliding-histogram +backend. It is used by the public generic, percentile, and bilateral rank-filter +APIs defined in ``_generic``, ``_percentile``, and ``_bilateral``. +""" + +import math + +import cupy as cp +import numpy as np + +import cucim.skimage._vendored._ndimage_filters_core as _filters_core +import cucim.skimage._vendored._ndimage_util as _util +from cucim.skimage._vendored._ndimage_filters import ( + __SHELL_SORT, + _get_shell_gap, +) + +from ._histogram import ( + _can_use_rank_histogram, + _rank_histogram, + _should_use_rank_histogram, +) + + +def _get_streaming_rank_kernel( + p0, + p1, + operation, + modes, + w_shape, + offsets, + cval, + int_type, + has_weights, + *, + has_mask=False, + dtype_max=255, + s0=0.0, + s1=0.0, +): + """Generate a rank kernel that reduces during neighborhood traversal. + + This path is for operations that do not need sorted neighborhood values. + It avoids allocating ``values`` and calling ``sort`` in the generated + kernel, which is much cheaper for large footprints. + """ + if operation in ("minimum", "maximum"): + best_var = "min_val" if operation == "minimum" else "max_val" + comparator = "<" if operation == "minimum" else ">" + pre = f"int n_vals = 0;\nX {best_var};" + update = f""" + X v = {{value}}; + if (n_vals == 0 || v {comparator} {best_var}) {{ + {best_var} = v; + }} + n_vals++; + """ + post = f""" + if (n_vals == 0) {{ + y = cast(x[i]); + return; + }} + y = cast({best_var}); + """ + elif operation in ("gradient", "enhance_contrast", "autolevel"): + pre = "int n_vals = 0;\nX min_val;\nX max_val;" + update = """ + X v = {value}; + if (n_vals == 0) { + min_val = v; + max_val = v; + } else { + if (v < min_val) min_val = v; + if (v > max_val) max_val = v; + } + n_vals++; + """ + if operation == "gradient": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(max_val - min_val); + """ + elif operation == "enhance_contrast": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + X g = x[i]; + if (max_val - g < g - min_val) { + y = cast(max_val); + } else { + y = cast(min_val); + } + """ + else: + post = f""" + if (n_vals == 0) {{ + y = cast(x[i]); + return; + }} + X g = x[i]; + X clamped = (g < min_val) ? min_val : + ((g > max_val) ? max_val : g); + double delta = static_cast(max_val - min_val); + if (delta > 0) {{ + double scaled = (static_cast(clamped - min_val) + / delta) * static_cast({dtype_max}); + y = cast(scaled); + }} else {{ + y = cast(0); + }} + """ + elif operation in ("mean", "sum", "subtract_mean", "threshold_mean"): + pre = "int n_vals = 0;\ndouble sum = 0.0;" + update = """ + X v = {value}; + sum += static_cast(v); + n_vals++; + """ + if operation == "mean": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(sum / n_vals); + """ + elif operation == "sum": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(sum); + """ + elif operation == "subtract_mean": + _mid_bin = (dtype_max + 1) // 2 + post = f""" + if (n_vals == 0) {{ + y = cast(x[i]); + return; + }} + double mean = sum / n_vals; + X g = x[i]; + y = cast((static_cast(g) - mean) * 0.5 + {_mid_bin}); + """ + else: + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + double mean = sum / n_vals; + X g = x[i]; + y = (static_cast(g) > mean) ? cast(1) : cast(0); + """ + elif operation == "pop": + pre = "int n_vals = 0;" + update = "n_vals++;" + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(n_vals); + """ + elif operation == "equalize": + pre = "int n_vals = 0;\nint eq_rank = 0;\nX g = x[i];" + update = """ + X v = {value}; + if (v <= g) eq_rank++; + n_vals++; + """ + post = f""" + if (n_vals == 0) {{ + y = cast(x[i]); + return; + }} + y = cast(static_cast({dtype_max}) * eq_rank / n_vals); + """ + elif operation == "geometric_mean": + pre = "int n_vals = 0;\ndouble log_sum = 0.0;" + update = """ + X v = {value}; + log_sum += log(static_cast(v) + 1.0); + n_vals++; + """ + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(round(exp(log_sum / n_vals) - 1.0)); + """ + elif operation == "noise_filter": + pre = """ + int n_vals = 0; + bool nf_found = false; + bool nf_has_dist = false; + typename RankNoiseDistance::type nf_min_dist = 0; + X g = x[i]; + """ + update = """ + X v = {value}; + if (v == g) { + nf_found = true; + } else { + typedef typename RankNoiseDistance::type DistanceT; + DistanceT vd = static_cast(v); + DistanceT gd = static_cast(g); + DistanceT d = (v > g) ? vd - gd : gd - vd; + if (!nf_has_dist || d < nf_min_dist) nf_min_dist = d; + nf_has_dist = true; + } + n_vals++; + """ + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = nf_found ? cast(0) : cast(nf_min_dist); + """ + elif operation in ("bilateral_mean", "bilateral_pop", "bilateral_sum"): + pre = """ + int n_vals = 0; + int bilat_pop = 0; + double bilat_sum = 0.0; + X g = x[i]; + double gd = static_cast(g); + """ + update = f""" + X v = {{value}}; + double vd = static_cast(v); + if (gd > (vd - {s0}) && gd < (vd + {s1})) {{ + bilat_pop++; + bilat_sum += vd; + }} + n_vals++; + """ + if operation == "bilateral_mean": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = (bilat_pop > 0) ? cast(bilat_sum / bilat_pop) : cast(0); + """ + elif operation == "bilateral_pop": + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = cast(bilat_pop); + """ + else: + post = """ + if (n_vals == 0) { + y = cast(x[i]); + return; + } + y = (bilat_pop > 0) ? cast(bilat_sum) : cast(0); + """ + else: + raise ValueError(f"Unsupported streaming operation: {operation}") + + update = update.replace("{value}", "__VALUE__") + update = update.replace("{", "{{").replace("}", "}}") + update = update.replace("__VALUE__", "{value}") + + if has_mask: + ndim = len(w_shape) + index_expr = " + ".join([f"ix_{j}" for j in range(ndim)]) + found = ( + "{{ ptrdiff_t _neighbor_idx = (" + index_expr + ") / sizeof(X); " + "if ((bool)mask[_neighbor_idx]) {{ " + update + " }} }}" + ) + else: + found = update + + op_name = operation.replace("_", "") + mask_str = "_masked" if has_mask else "" + preamble = "" + if operation == "noise_filter": + preamble = r""" +template ::value> +struct RankNoiseDistance { + typedef typename std::make_unsigned::type type; +}; + +template +struct RankNoiseDistance { + typedef double type; +}; +""" + return _filters_core._generate_nd_kernel( + f"rank_stream_{op_name}_{int(p0)}_{int(p1)}{mask_str}", + pre, + found, + post, + modes, + w_shape, + int_type, + offsets, + cval, + has_weights=has_weights, + has_mask=has_mask, + preamble=preamble, + ) + + +@cp.memoize(for_each_device=True) +def _get_percentile_range_kernel( + filter_size, + p0, + p1, + operation, + modes, + w_shape, + offsets, + cval, + int_type, + has_weights, + *, + has_mask=False, + dtype_max=255, + s0=0.0, + s1=0.0, +): + """Generate a kernel for computing statistics on a percentile range. + + Parameters + ---------- + filter_size : int + Total number of values in the neighborhood (when mask is not used). + p0 : float + Lower percentile (0-100). + p1 : float + Upper percentile (0-100). + operation : str + Operation to perform. Supported values are ``'autolevel'``, + ``'bilateral_mean'``, ``'bilateral_pop'``, ``'bilateral_sum'``, + ``'enhance_contrast'``, ``'entropy'``, ``'equalize'``, + ``'geometric_mean'``, ``'gradient'``, ``'maximum'``, ``'mean'``, + ``'minimum'``, ``'modal'``, ``'noise_filter'``, ``'percentile'``, + ``'pop'``, ``'subtract_mean'``, ``'sum'``, ``'threshold'``, and + ``'threshold_mean'``. + modes : tuple of str + Boundary handling modes. + w_shape : tuple of int + Shape of the footprint/kernel. + offsets : tuple of int + Offsets for the footprint origin. + cval : float + Constant value for 'constant' mode. + int_type : str + Integer type to use for indexing. + has_weights : bool + Whether an explicit footprint array is used. + has_mask : bool + Whether an image mask is used to filter neighborhood pixels. + dtype_max : scalar + Maximum representable output value used by scaled operations. + s0, s1 : float + Bilateral graylevel-range parameters. + + Returns + ------- + kernel : cupy.ElementwiseKernel + The compiled CUDA kernel. + """ + # Convert percentiles to array indices + # Note: Matches scikit-image's histogram-based percentile approach + # where values are included if cumsum is in [p0*pop, p1*pop] + _single_percentile_op = operation in ("percentile", "threshold") + _bilateral_op = operation in ( + "bilateral_mean", + "bilateral_pop", + "bilateral_sum", + ) + _full_range = p0 <= 0 and p1 >= 100 + _streaming_ops = { + "bilateral_mean", + "bilateral_pop", + "bilateral_sum", + "equalize", + "geometric_mean", + "minimum", + "maximum", + "noise_filter", + "threshold_mean", + } + _full_range_streaming_ops = { + "autolevel", + "enhance_contrast", + "gradient", + "mean", + "pop", + "subtract_mean", + "sum", + } + _skip_idx = _single_percentile_op or _bilateral_op + if not _bilateral_op: + if p0 < 0 or p0 > 100: + raise ValueError("Percentiles must be in range [0, 100]") + if not _single_percentile_op: + if p1 < 0 or p1 > 100: + raise ValueError("Percentiles must be in range [0, 100]") + if p0 >= p1: + raise ValueError("p0 must be less than p1") + + if operation in _streaming_ops or ( + _full_range and operation in _full_range_streaming_ops + ): + return _get_streaming_rank_kernel( + p0, + p1, + operation, + modes, + w_shape, + offsets, + cval, + int_type, + has_weights, + has_mask=has_mask, + dtype_max=dtype_max, + s0=s0, + s1=s1, + ) + + # When has_mask is True, we need to dynamically calculate indices based on + # the actual number of values collected (which depends on the mask). + # We'll use runtime calculation in the CUDA code. + # "percentile", "threshold", and bilateral ops compute their own indices, + # so idx_start/idx_end are not needed. + if not has_mask and not _skip_idx: + # Calculate indices for the percentile range + # (pre-computed at compile time) + # Matches scikit-image's histogram-based approach where value at + # index i is included if: (i + 1) >= p0 * N AND (i + 1) <= p1 * N + # Rearranging: p0 * N - 1 <= i <= p1 * N - 1 + # + # For integer indices in a sorted array (using Python's range with + # exclusive upper bound): + # idx_start = ceil(p0 * N - 1) = ceil(p0 * N) - 1 + # idx_end = floor(p1 * N - 1) + 1 = floor(p1 * N) + idx_start = max(0, int(math.ceil(p0 * filter_size / 100.0)) - 1) + # int() gives floor for positive values + idx_end = int(p1 * filter_size / 100.0) + + # Ensure at least one value is included + if idx_end <= idx_start: + idx_end = idx_start + 1 + idx_end = min(idx_end, filter_size) + n_values = idx_end - idx_start + + # Always use full sorting for percentile ranges + array_size = filter_size + sorter = __SHELL_SORT.format(gap=_get_shell_gap(filter_size)) + + if has_mask: + # Runtime calculation of indices based on actual count + if operation in ("percentile", "threshold", "pop") or _bilateral_op: + post = """ + if (iv == 0) {{ + y = cast(x[i]); // No valid values, keep original + return; + }} + sort(values, iv);""" + else: + post = f""" + if (iv == 0) {{ + y = cast(x[i]); // No valid values, keep original + return; + }} + sort(values, iv); + int actual_start = max(0, (int)ceil({p0 / 100.0} * iv) - 1); + int actual_end = (int)({p1 / 100.0} * iv); + if (actual_end <= actual_start) actual_end = actual_start + 1; + if (actual_end > iv) actual_end = iv;""" + else: + post = f""" + sort(values, {filter_size});""" + + # Generate the post-processing code based on the operation + if operation == "mean": + # Standard mean of values in percentile range + if has_mask: + # Runtime calculation of indices based on actual count + post += """ + int n_vals = actual_end - actual_start; + double sum = 0.0; + for (int j = actual_start; j < actual_end; j++) {{ + sum += static_cast(values[j]); + }} + y = cast(sum / n_vals); + """ + else: + post += f""" + double sum = 0.0; + for (int j = {idx_start}; j < {idx_end}; j++) {{ + sum += static_cast(values[j]); + }} + y = cast(sum / {n_values}); + """ + elif operation == "sum": + # Sum of values in percentile range + if has_mask: + post += """ + double sum = 0.0; + for (int j = actual_start; j < actual_end; j++) {{ + sum += static_cast(values[j]); + }} + y = cast(sum); + """ + else: + post += f""" + double sum = 0.0; + for (int j = {idx_start}; j < {idx_end}; j++) {{ + sum += static_cast(values[j]); + }} + y = cast(sum); + """ + elif operation == "gradient": + # Gradient: max - min in percentile range + if has_mask: + post += """ + X min_val = values[actual_start]; + X max_val = values[actual_end - 1]; + y = cast(max_val - min_val); + """ + else: + post += f""" + X min_val = values[{idx_start}]; + X max_val = values[{idx_end - 1}]; + y = cast(max_val - min_val); + """ + elif operation == "subtract_mean": + # Subtract mean: scikit-image formula: + # (g - mean) * 0.5 + mid_bin + # where mid_bin = n_bins / 2 (128 for uint8, 32768 for uint16). + # This centers the result so that g == mean maps to mid_bin. + _mid_bin = (dtype_max + 1) // 2 + if has_mask: + post += f""" + int n_vals = actual_end - actual_start; + double sum = 0.0; + for (int j = actual_start; j < actual_end; j++) {{ + sum += static_cast(values[j]); + }} + double mean = sum / n_vals; + X g = x[i]; + y = cast((static_cast(g) - mean) * 0.5 + {_mid_bin}); + """ + else: + post += f""" + double sum = 0.0; + for (int j = {idx_start}; j < {idx_end}; j++) {{ + sum += static_cast(values[j]); + }} + double mean = sum / {n_values}; + X g = x[i]; + y = cast((static_cast(g) - mean) * 0.5 + {_mid_bin}); + """ + elif operation == "enhance_contrast": + # Enhance contrast: replace with closer extreme (min or max) + if has_mask: + post += """ + X min_val = values[actual_start]; + X max_val = values[actual_end - 1]; + X g = x[i]; + // Replace with whichever extreme is closer + if (max_val - g < g - min_val) {{ + y = cast(max_val); + }} else {{ + y = cast(min_val); + }} + """ + else: + post += f""" + X min_val = values[{idx_start}]; + X max_val = values[{idx_end - 1}]; + X g = x[i]; + if (max_val - g < g - min_val) {{ + y = cast(max_val); + }} else {{ + y = cast(min_val); + }} + """ + elif operation == "percentile": + # Single percentile value (p0 determines which percentile) + # Note: This returns the value AT the p0 percentile + if has_mask: + post += f""" + int percentile_idx; + if ({p0 / 100.0} == 1.0) {{ + // p0 = 100%: return maximum + percentile_idx = iv - 1; + }} else {{ + // Find index where cumsum > p0 * pop + percentile_idx = (int)({p0 / 100.0} * iv); + if (percentile_idx >= iv) percentile_idx = iv - 1; + }} + y = cast(values[percentile_idx]); + """ + else: + # For no mask, we can use precomputed idx_start + post += f""" + int percentile_idx; + if ({p0 / 100.0} == 1.0) {{ + percentile_idx = {filter_size - 1}; + }} else {{ + percentile_idx = (int)({p0 / 100.0} * {filter_size}); + if (percentile_idx >= {filter_size}) {{ + percentile_idx = {filter_size - 1}; + }} + }} + y = cast(values[percentile_idx]); + """ + elif operation == "pop": + # Population: count of pixels in percentile range. + # Must match scikit-image's histogram-bin grouping: groups of equal + # values are included/excluded as a whole based on whether the + # cumulative count after adding the group falls in [p0*pop, p1*pop]. + if has_mask: + post += f""" + int pop_n = 0; + int pop_cumsum = 0; + int pop_j = 0; + while (pop_j < iv) {{ + X pop_val = values[pop_j]; + int pop_gs = 0; + while (pop_j < iv && values[pop_j] == pop_val) {{ + pop_j++; + pop_gs++; + }} + pop_cumsum += pop_gs; + if ((double)pop_cumsum >= {p0 / 100.0} * iv && + (double)pop_cumsum <= {p1 / 100.0} * iv) {{ + pop_n += pop_gs; + }} + }} + y = cast(pop_n); + """ + else: + post += f""" + int pop_n = 0; + int pop_cumsum = 0; + int pop_j = 0; + while (pop_j < {filter_size}) {{ + X pop_val = values[pop_j]; + int pop_gs = 0; + while (pop_j < {filter_size} && values[pop_j] == pop_val) {{ + pop_j++; + pop_gs++; + }} + pop_cumsum += pop_gs; + if ((double)pop_cumsum >= {p0 / 100.0} * {filter_size} && + (double)pop_cumsum <= {p1 / 100.0} * {filter_size}) {{ + pop_n += pop_gs; + }} + }} + y = cast(pop_n); + """ + elif operation == "threshold": + # Threshold: binary output comparing center pixel to p0 percentile. + # scikit-image uses (n_bins - 1) * (g >= threshold), which gives the + # dtype max value (e.g. 255 for uint8) or 0. + if has_mask: + post += f""" + int threshold_idx = (int)({p0 / 100.0} * iv); + if (threshold_idx >= iv) threshold_idx = iv - 1; + X threshold_val = values[threshold_idx]; + X g = x[i]; + y = (g >= threshold_val) ? cast({dtype_max}) : cast(0); + """ + else: + post += f""" + int threshold_idx = (int)({p0 / 100.0} * {filter_size}); + if (threshold_idx >= {filter_size}) {{ + threshold_idx = {filter_size - 1}; + }} + X threshold_val = values[threshold_idx]; + X g = x[i]; + y = (g >= threshold_val) ? cast({dtype_max}) : cast(0); + """ + elif operation == "autolevel": + # Autolevel: stretch pixel values to full dtype range based on local + # percentile min/max. scikit-image formula: + # (n_bins - 1) * (clamp(g, imin, imax) - imin) / (imax - imin) + # Scales output to [0, dtype_max], NOT [0, local_max]. + if has_mask: + post += f""" + X min_val = values[actual_start]; + X max_val = values[actual_end - 1]; + X g = x[i]; + X clamped = (g < min_val) ? min_val : \ +((g > max_val) ? max_val : g); + double delta = static_cast(max_val - min_val); + if (delta > 0) {{ + double scaled = (static_cast(clamped - min_val) \ +/ delta) * static_cast({dtype_max}); + y = cast(scaled); + }} else {{ + y = cast(0); + }} + """ + else: + post += f""" + X min_val = values[{idx_start}]; + X max_val = values[{idx_end - 1}]; + X g = x[i]; + X clamped = (g < min_val) ? min_val : \ +((g > max_val) ? max_val : g); + double delta = static_cast(max_val - min_val); + if (delta > 0) {{ + double scaled = (static_cast(clamped - min_val) \ +/ delta) * static_cast({dtype_max}); + y = cast(scaled); + }} else {{ + y = cast(0); + }} + """ + elif operation == "modal": + # Modal: most frequent value (mode) in the neighborhood. + # Scan sorted array for the longest run of equal values. + if has_mask: + post += """ + X mode_val = values[0]; + int mode_max = 1; + int mode_cur = 1; + for (int j = 1; j < iv; j++) { + if (values[j] == values[j - 1]) { + mode_cur++; + } else { + if (mode_cur > mode_max) { + mode_max = mode_cur; + mode_val = values[j - 1]; + } + mode_cur = 1; + } + } + if (mode_cur > mode_max) mode_val = values[iv - 1]; + y = cast(mode_val); + """ + else: + post += f""" + X mode_val = values[0]; + int mode_max = 1; + int mode_cur = 1; + for (int j = 1; j < {filter_size}; j++) {{ + if (values[j] == values[j - 1]) {{ + mode_cur++; + }} else {{ + if (mode_cur > mode_max) {{ + mode_max = mode_cur; + mode_val = values[j - 1]; + }} + mode_cur = 1; + }} + }} + if (mode_cur > mode_max) mode_val = values[{filter_size} - 1]; + y = cast(mode_val); + """ + elif operation == "entropy": + # Shannon entropy in bits: -sum(p * log2(p)) where p = count / N. + # Uses run-length counting on the sorted array to get value counts. + # log2(p) = log(p) / log(2); 0.6931471805599453 = log(2). + if has_mask: + post += """ + double ent = 0.0; + int ent_j = 0; + while (ent_j < iv) { + int ent_count = 1; + while (ent_j + ent_count < iv && + values[ent_j + ent_count] == values[ent_j]) + ent_count++; + double p = static_cast(ent_count) / iv; + ent -= p * log(p) / 0.6931471805599453; + ent_j += ent_count; + } + y = cast(ent); + """ + else: + post += f""" + double ent = 0.0; + int ent_j = 0; + while (ent_j < {filter_size}) {{ + int ent_count = 1; + while (ent_j + ent_count < {filter_size} && + values[ent_j + ent_count] == values[ent_j]) + ent_count++; + double p = static_cast(ent_count) / {filter_size}; + ent -= p * log(p) / 0.6931471805599453; + ent_j += ent_count; + }} + y = cast(ent); + """ + else: + raise ValueError( + f"Unsupported operation: {operation}. " + "Supported sorted operations: 'mean', 'sum', 'gradient', " + "'subtract_mean', 'enhance_contrast', 'percentile', 'pop', " + "'threshold', 'autolevel', 'modal', 'entropy'" + ) + + # Sanitize operation name for kernel name (replace special chars) + op_name = operation.replace("_", "") + + # Build the pre string and found string with neighborhood-level masking + pre = f"int iv = 0;\nX values[{array_size}];" + + if has_mask: + # Neighborhood-level masking: check mask at each neighbor's position + # Calculate the neighbor's element index from the byte offset. + # Note: ix_{j} are calculated in _generate_nd_kernel (in + # _ndimage_filters_core.py) within the neighborhood iteration loops. + # Each ix_{j} = coordinate_{j} * xstride_{j}, where xstride_{j} comes + # from x.strides()[{j}] (byte stride). The sum (ix_0 + ix_1 + ...) + # gives the total byte offset to the neighbor pixel from the start + # of the array. + ndim = len(w_shape) + index_expr = " + ".join([f"ix_{j}" for j in range(ndim)]) + # Use string concatenation (not f-string) so that {{ / }} are + # correctly interpreted as literal braces by .format() later. + found = ( + "{{ ptrdiff_t _neighbor_idx = (" + index_expr + ") / sizeof(X); " + "if ((bool)mask[_neighbor_idx]) {{ " + "values[iv++] = {value}; " + "}} }}" + ) + else: + found = "values[iv++] = {value};" + + mask_str = "_masked" if has_mask else "" + return _filters_core._generate_nd_kernel( + f"percentile_range_{filter_size}_{int(p0)}_{int(p1)}_{op_name}{mask_str}", + pre, + found, + post, + modes, + w_shape, + int_type, + offsets, + cval, + has_weights=has_weights, + has_mask=has_mask, + preamble=sorter, + ) + + +def _is_decomposed_footprint(footprint): + """Return True for morphology footprint decomposition sequences.""" + if not isinstance(footprint, (tuple, list)) or len(footprint) == 0: + return False + + for item in footprint: + if not isinstance(item, (tuple, list)) or len(item) != 2: + return False + footprint_part, num_iter = item + if not hasattr(footprint_part, "ndim"): + return False + if not isinstance(num_iter, (int, np.integer)): + return False + + return True + + +def _skimage_rank_filter( + input, + p0, + p1, + operation="mean", + size=None, + footprint=None, + output=None, + mode="reflect", + cval=0.0, + origin=0, + axes=None, + *, + mask=None, + s0=0, + s1=0, + backend="auto", +): + """Internal helper for percentile range filters. + + This function computes statistics (mean, sum, etc.) on values within + a specified percentile range [p0, p1] of the neighborhood. + + Parameters + ---------- + input : cupy.ndarray + The input array. + p0 : float + Lower percentile (0-100). + p1 : float + Upper percentile (0-100). + operation : str, optional + Operation to perform. Supported values are ``'autolevel'``, + ``'bilateral_mean'``, ``'bilateral_pop'``, ``'bilateral_sum'``, + ``'enhance_contrast'``, ``'entropy'``, ``'equalize'``, + ``'geometric_mean'``, ``'gradient'``, ``'maximum'``, ``'mean'``, + ``'minimum'``, ``'modal'``, ``'noise_filter'``, ``'percentile'``, + ``'pop'``, ``'subtract_mean'``, ``'sum'``, ``'threshold'``, and + ``'threshold_mean'``. Default is ``'mean'``. + size : int or sequence of int, optional + Size of the neighborhood. One of `size` or `footprint` must be provided. + footprint : cupy.ndarray, optional + Boolean array specifying the neighborhood shape. + output : cupy.ndarray, dtype or None, optional + The array in which to place the output. + mode : str or sequence of str, optional + Boundary handling mode. Default is ``'reflect'``. + cval : scalar, optional + Value to fill past edges if mode is 'constant'. Default is 0.0. + origin : int or sequence of int, optional + Origin of the footprint. Default is 0. + axes : tuple of int or None, optional + Axes along which to apply the filter. Default is None (all axes). + mask : cupy.ndarray or None, optional + If provided, only neighbor pixels where mask is True are included + when computing statistics. This matches scikit-image's filters.rank + behavior where the mask filters which pixels in the local neighborhood + contribute to the computation. Output is computed for all pixels, but + each uses a different set of neighbors based on the mask. + s0, s1 : float, optional + Bilateral graylevel-range parameters. Default is 0. + backend : {'auto', 'histogram', 'elementwise'}, optional + ``'auto'`` selects the histogram backend for compatible calls above + the tuned footprint-size threshold and otherwise uses the elementwise + backend. ``'histogram'`` requires a compatible uint8 2-D call and + raises ``ValueError`` otherwise. Default is ``'auto'``. + + Returns + ------- + output : cupy.ndarray + The filtered array. + + Notes + ----- + This function is for internal use as a common implementation for filters + under cucim.skimage.filters.rank. + """ + if backend not in ("auto", "histogram", "elementwise"): + raise ValueError( + "backend must be one of 'auto', 'histogram' or 'elementwise'" + ) + + ndim = input.ndim + axes = _util._check_axes(axes, ndim) + num_axes = len(axes) + default_footprint = footprint is None + if _is_decomposed_footprint(footprint): + raise ValueError( + "decomposed footprint sequences are not supported by rank filters" + ) + sizes, footprint, _ = _filters_core._check_size_footprint_structure( + num_axes, + size, + footprint, + None, + force_footprint=operation == "noise_filter", + ) + if cval is cp.nan: + raise NotImplementedError("NaN cval is unsupported") + + # Validate percentiles + p0 = float(p0) + p1 = float(p1) + _bilateral_op = operation in ( + "bilateral_mean", + "bilateral_pop", + "bilateral_sum", + ) + _single_percentile_op = operation in ("percentile", "threshold") + if not _bilateral_op: + if p0 < 0 or p0 > 100: + raise ValueError("Percentiles must be in range [0, 100]") + if not _single_percentile_op: + if p1 < 0 or p1 > 100: + raise ValueError("Percentiles must be in range [0, 100]") + if p0 >= p1: + raise ValueError("p0 must be less than p1") + + has_weights = True + if sizes is not None: + has_weights = False + filter_size = math.prod(sizes) + if filter_size == 0: + return cp.zeros_like(input) + footprint_shape = tuple(sizes) + ( + axes, + footprint, + origins, + modes, + int_type, + ) = _filters_core._check_nd_args( + input, + None, + mode, + origin, + "footprint", + axes=axes, + sizes=footprint_shape, + ) + else: + if footprint.size == 0: + return cp.zeros_like(input) + + ( + axes, + footprint, + origins, + modes, + int_type, + ) = _filters_core._check_nd_args( + input, footprint, mode, origin, "footprint", axes=axes + ) + + if operation == "noise_filter": + # The footprint anchor addresses the center input pixel. Exclude + # it so that it cannot make every pixel appear non-isolated. + footprint = footprint.copy() + anchor = _filters_core._origins_to_offsets(origins, footprint.shape) + footprint[anchor] = False + + if default_footprint: + filter_size = footprint.size + else: + footprint_shape = footprint.shape + filter_size = int(footprint.sum()) + if filter_size == footprint.size: + # can omit passing the footprint if it is all ones + sizes = footprint.shape + has_weights = False + + if not has_weights: + footprint = None + + offsets = _filters_core._origins_to_offsets(origins, footprint_shape) + if num_axes < ndim and not has_weights: + offsets = tuple(_util._expand_origin(ndim, axes, offsets)) + modes = tuple(_util._expand_mode(ndim, axes, modes)) + footprint_shape_temp = [1] * ndim + for s, ax in zip(footprint_shape, axes): + footprint_shape_temp[ax] = s + footprint_shape = tuple(footprint_shape_temp) + + has_mask = mask is not None + + # Compute dtype max for threshold operation (binary output needs + # the type's max value, not the local neighborhood max). + _out_dtype = output.dtype if output is not None else input.dtype + if np.issubdtype(_out_dtype, np.integer): + _dtype_max = int(np.iinfo(_out_dtype).max) + else: + _dtype_max = 1.0 + + can_use_histogram = _can_use_rank_histogram( + input, + footprint_shape, + output, + mask, + modes, + origins, + has_weights=has_weights, + operation=operation, + p0=p0, + p1=p1, + ) + if backend == "histogram" and not can_use_histogram: + raise ValueError( + "backend='histogram' requires a supported uint8 2D rank " + "operation, compatible output, no mask, zero shifts, reflect " + "mode, and an all-ones odd rectangular footprint" + ) + + if backend == "histogram" or ( + backend == "auto" + and can_use_histogram + and _should_use_rank_histogram(operation, footprint_shape) + ): + return _rank_histogram( + input, + footprint_shape, + operation, + output=output, + mode=modes[0], + cval=cval, + p0=p0, + p1=p1, + s0=s0, + s1=s1, + dtype_max=_dtype_max, + ) + + kernel = _get_percentile_range_kernel( + filter_size, + p0, + p1, + operation, + modes, + footprint_shape, + offsets, + float(cval), + int_type, + has_weights=has_weights, + has_mask=has_mask, + dtype_max=_dtype_max, + s0=float(s0), + s1=float(s1), + ) + kwargs = dict(weights_dtype=bool) + if has_mask: + kwargs["mask"] = mask + return _filters_core._call_kernel( + kernel, input, footprint, output, **kwargs + ) diff --git a/python/cucim/src/cucim/skimage/filters/rank/cuda/histogram_rank.cu b/python/cucim/src/cucim/skimage/filters/rank/cuda/histogram_rank.cu new file mode 100644 index 000000000..5d1f2a71d --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/cuda/histogram_rank.cu @@ -0,0 +1,755 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ + +#define OP_PERCENTILE 0 +#define OP_THRESHOLD 1 +#define OP_MEAN 2 +#define OP_SUM 3 +#define OP_POP 4 +#define OP_GRADIENT 5 +#define OP_AUTOLEVEL 6 +#define OP_ENTROPY 7 +#define OP_ENHANCE_CONTRAST 8 +#define OP_SUBTRACT_MEAN 9 +#define OP_EQUALIZE 10 +#define OP_BILATERAL_MEAN 11 +#define OP_BILATERAL_POP 12 +#define OP_BILATERAL_SUM 13 +#define OP_MODAL 14 +#define OP_GEOMETRIC_MEAN 15 + +#ifndef HIST_COUNTER_T +# define HIST_COUNTER_T int +#endif + +#ifndef RANK_HIST_OUTPUT_T +# define RANK_HIST_OUTPUT_T unsigned char +#endif + +#if RANK_HIST_OP == OP_GEOMETRIC_MEAN +__device__ __constant__ double geometricMeanLogLut[256] = { 0, + 0.69314718055994529, + 1.0986122886681098, + 1.3862943611198906, + 1.6094379124341003, + 1.791759469228055, + 1.9459101490553132, + 2.0794415416798357, + 2.1972245773362196, + 2.3025850929940459, + 2.3978952727983707, + 2.4849066497880004, + 2.5649493574615367, + 2.6390573296152584, + 2.7080502011022101, + 2.7725887222397811, + 2.8332133440562162, + 2.8903717578961645, + 2.9444389791664403, + 2.9957322735539909, + 3.044522437723423, + 3.0910424533583161, + 3.1354942159291497, + 3.1780538303479458, + 3.2188758248682006, + 3.2580965380214821, + 3.2958368660043291, + 3.3322045101752038, + 3.3672958299864741, + 3.4011973816621555, + 3.4339872044851463, + 3.4657359027997265, + 3.4965075614664802, + 3.5263605246161616, + 3.5553480614894135, + 3.5835189384561099, + 3.6109179126442243, + 3.6375861597263857, + 3.6635616461296463, + 3.6888794541139363, + 3.713572066704308, + 3.7376696182833684, + 3.7612001156935624, + 3.784189633918261, + 3.8066624897703196, + 3.8286413964890951, + 3.8501476017100584, + 3.8712010109078911, + 3.8918202981106265, + 3.912023005428146, + 3.9318256327243257, + 3.9512437185814275, + 3.970291913552122, + 3.9889840465642745, + 4.0073331852324712, + 4.0253516907351496, + 4.0430512678345503, + 4.0604430105464191, + 4.0775374439057197, + 4.0943445622221004, + 4.1108738641733114, + 4.1271343850450917, + 4.1431347263915326, + 4.1588830833596715, + 4.1743872698956368, + 4.1896547420264252, + 4.2046926193909657, + 4.219507705176107, + 4.2341065045972597, + 4.2484952420493594, + 4.2626798770413155, + 4.2766661190160553, + 4.290459441148391, + 4.3040650932041702, + 4.3174881135363101, + 4.3307333402863311, + 4.3438054218536841, + 4.3567088266895917, + 4.3694478524670215, + 4.3820266346738812, + 4.3944491546724391, + 4.4067192472642533, + 4.4188406077965983, + 4.4308167988433134, + 4.4426512564903167, + 4.4543472962535073, + 4.4659081186545837, + 4.4773368144782069, + 4.4886363697321396, + 4.499809670330265, + 4.5108595065168497, + 4.5217885770490405, + 4.5325994931532563, + 4.5432947822700038, + 4.5538768916005408, + 4.5643481914678361, + 4.5747109785033828, + 4.5849674786705723, + 4.5951198501345898, + 4.6051701859880918, + 4.6151205168412597, + 4.6249728132842707, + 4.6347289882296359, + 4.6443908991413725, + 4.6539603501575231, + 4.6634390941120669, + 4.6728288344619058, + 4.6821312271242199, + 4.6913478822291435, + 4.7004803657924166, + 4.7095302013123339, + 4.7184988712950942, + 4.7273878187123408, + 4.7361984483944957, + 4.7449321283632502, + 4.7535901911063645, + 4.7621739347977563, + 4.7706846244656651, + 4.7791234931115296, + 4.7874917427820458, + 4.7957905455967413, + 4.8040210447332568, + 4.8121843553724171, + 4.8202815656050371, + 4.8283137373023015, + 4.836281906951478, + 4.8441870864585912, + 4.8520302639196169, + 4.8598124043616719, + 4.8675344504555822, + 4.8751973232011512, + 4.8828019225863706, + 4.8903491282217537, + 4.8978397999509111, + 4.9052747784384296, + 4.9126548857360524, + 4.9199809258281251, + 4.9272536851572051, + 4.9344739331306915, + 4.9416424226093039, + 4.9487598903781684, + 4.9558270576012609, + 4.962844630259907, + 4.9698132995760007, + 4.9767337424205742, + 4.9836066217083363, + 4.990432586778736, + 4.9972122737641147, + 5.0039463059454592, + 5.0106352940962555, + 5.0172798368149243, + 5.0238805208462765, + 5.0304379213924353, + 5.0369526024136295, + 5.0434251169192468, + 5.0498560072495371, + 5.0562458053483077, + 5.0625950330269669, + 5.0689042022202315, + 5.0751738152338266, + 5.0814043649844631, + 5.0875963352323836, + 5.0937502008067623, + 5.0998664278241987, + 5.1059454739005803, + 5.1119877883565437, + 5.1179938124167554, + 5.1239639794032588, + 5.1298987149230735, + 5.1357984370502621, + 5.1416635565026603, + 5.1474944768134527, + 5.1532915944977793, + 5.1590552992145291, + 5.1647859739235145, + 5.1704839950381514, + 5.1761497325738288, + 5.181783550292085, + 5.1873858058407549, + 5.1929568508902104, + 5.1984970312658261, + 5.2040066870767951, + 5.2094861528414214, + 5.2149357576089859, + 5.2203558250783244, + 5.2257466737132017, + 5.2311086168545868, + 5.2364419628299492, + 5.2417470150596426, + 5.2470240721604862, + 5.2522734280466299, + 5.2574953720277815, + 5.2626901889048856, + 5.2678581590633282, + 5.2729995585637468, + 5.2781146592305168, + 5.2832037287379885, + 5.2882670306945352, + 5.2933048247244923, + 5.2983173665480363, + 5.3033049080590757, + 5.3082676974012051, + 5.3132059790417872, + 5.3181199938442161, + 5.3230099791384085, + 5.3278761687895813, + 5.3327187932653688, + 5.3375380797013179, + 5.3423342519648109, + 5.3471075307174685, + 5.3518581334760666, + 5.3565862746720123, + 5.3612921657094255, + 5.3659760150218512, + 5.3706380281276624, + 5.3752784076841653, + 5.3798973535404597, + 5.3844950627890888, + 5.389071729816501, + 5.393627546352362, + 5.3981627015177525, + 5.4026773818722793, + 5.4071717714601188, + 5.4116460518550396, + 5.4161004022044201, + 5.4205349992722862, + 5.4249500174814029, + 5.4293456289544411, + 5.43372200355424, + 5.4380793089231956, + 5.4424177105217932, + 5.4467373716663099, + 5.4510384535657002, + 5.4553211153577017, + 5.4595855141441589, + 5.4638318050256105, + 5.4680601411351315, + 5.472270673671475, + 5.476463551931511, + 5.4806389233419912, + 5.4847969334906548, + 5.4889377261566867, + 5.4930614433405482, + 5.4971682252932021, + 5.5012582105447274, + 5.5053315359323625, + 5.5093883366279774, + 5.5134287461649825, + 5.5174528964647074, + 5.521460917862246, + 5.5254529391317835, + 5.5294290875114234, + 5.5333894887275203, + 5.5373342670185366, + 5.5412635451584258, + 5.5451774444795623 }; +#endif + +__device__ void histogramPrefixScan256(int* hist, int* scan) +{ + int tx = threadIdx.x; + if (tx < 256) + { + scan[tx] = hist[tx]; + } + __syncthreads(); + + for (int offset = 1; offset < 256; offset <<= 1) + { + int v = 0; + if (tx >= offset && tx < 256) + { + v = scan[tx - offset]; + } + __syncthreads(); + if (tx >= offset && tx < 256) + { + scan[tx] += v; + } + __syncthreads(); + } +} + +__device__ void reduceSum256(int* values) +{ + int tx = threadIdx.x; + for (int stride = 128; stride > 0; stride >>= 1) + { + if (tx < stride) + { + values[tx] += values[tx + stride]; + } + __syncthreads(); + } +} + +__device__ void histogramWeightedPrefixScan256(int* hist, int* scan) +{ + int tx = threadIdx.x; + if (tx < 256) + { + scan[tx] = hist[tx] * tx; + } + __syncthreads(); + + for (int offset = 1; offset < 256; offset <<= 1) + { + int v = 0; + if (tx >= offset && tx < 256) + { + v = scan[tx - offset]; + } + __syncthreads(); + if (tx >= offset && tx < 256) + { + scan[tx] += v; + } + __syncthreads(); + } +} + +__device__ RANK_HIST_OUTPUT_T histogramRankValue(int* hist, + int* scan, + int* tmp0, + int* tmp1, + double* dtmp, + int op, + int window_size, + double p0, + double p1, + double s0, + double s1, + double dtype_max, + unsigned char center) +{ + int tx = threadIdx.x; + __shared__ int result; + __shared__ int range_start; + __shared__ int range_end; + // clang-format off +#if RANK_HIST_OP == OP_MEAN || RANK_HIST_OP == OP_SUM || RANK_HIST_OP == OP_SUBTRACT_MEAN || RANK_HIST_OP == OP_BILATERAL_MEAN || RANK_HIST_OP == OP_BILATERAL_POP || RANK_HIST_OP == OP_BILATERAL_SUM + // clang-format on + __shared__ int range_start_sum; + __shared__ int range_end_sum; +#endif + +#if RANK_HIST_OP == OP_ENTROPY + double ent = 0.0; + if (tx < 256 && hist[tx] > 0) + { + double p = ((double)hist[tx]) / window_size; + ent = -p * log(p) / 0.6931471805599453; + } + dtmp[tx] = ent; + __syncthreads(); + for (int stride = 128; stride > 0; stride >>= 1) + { + if (tx < stride) + { + dtmp[tx] += dtmp[tx + stride]; + } + __syncthreads(); + } + return static_cast(dtmp[0]); +#elif RANK_HIST_OP == OP_MODAL + tmp0[tx] = hist[tx]; + tmp1[tx] = tx; + __syncthreads(); + for (int stride = 128; stride > 0; stride >>= 1) + { + if (tx < stride) + { + int other_count = tmp0[tx + stride]; + int other_value = tmp1[tx + stride]; + if (other_count > tmp0[tx] || (other_count == tmp0[tx] && other_value < tmp1[tx])) + { + tmp0[tx] = other_count; + tmp1[tx] = other_value; + } + } + __syncthreads(); + } + return static_cast(tmp1[0]); +#elif RANK_HIST_OP == OP_GEOMETRIC_MEAN + double log_sum = 0.0; + if (tx < 256 && hist[tx] > 0) + { + log_sum = ((double)hist[tx]) * geometricMeanLogLut[tx]; + } + dtmp[tx] = log_sum; + __syncthreads(); + for (int stride = 128; stride > 0; stride >>= 1) + { + if (tx < stride) + { + dtmp[tx] += dtmp[tx + stride]; + } + __syncthreads(); + } + return static_cast(round(exp(dtmp[0] / window_size) - 1.0)); +#else + op = RANK_HIST_OP; + histogramPrefixScan256(hist, scan); + int pop = scan[255]; + + if (tx == 0) + { + result = 0; + int start = max(0, (int)ceil(p0 * pop / 100.0) - 1); + int end = (int)(p1 * pop / 100.0); + if (end <= start) + { + end = start + 1; + } + if (end > pop) + { + end = pop; + } + range_start = start; + range_end = end; + } + __syncthreads(); + + if (op == OP_PERCENTILE || op == OP_THRESHOLD) + { + int target; + if (p0 == 100.0) + { + target = pop - 1; + } + else + { + target = (int)(p0 * pop / 100.0); + if (target >= pop) + { + target = pop - 1; + } + } + + if (tx < 256 && hist[tx] > 0) + { + int bin_start = scan[tx] - hist[tx]; + if (bin_start <= target && scan[tx] > target) + { + result = tx; + } + } + __syncthreads(); + + if (op == OP_THRESHOLD) + { + return (center >= result) ? static_cast(dtype_max) : static_cast(0); + } + return static_cast(result); + } + +# if RANK_HIST_OP == OP_EQUALIZE + return static_cast(dtype_max * ((double)scan[center]) / pop); +# endif + + // clang-format off +#if RANK_HIST_OP == OP_MEAN || RANK_HIST_OP == OP_SUM || RANK_HIST_OP == OP_SUBTRACT_MEAN || RANK_HIST_OP == OP_BILATERAL_MEAN || RANK_HIST_OP == OP_BILATERAL_POP || RANK_HIST_OP == OP_BILATERAL_SUM + // clang-format on + histogramWeightedPrefixScan256(hist, tmp1); + if (tx == 0) + { + range_start_sum = 0; + range_end_sum = 0; + } + __syncthreads(); + +# if RANK_HIST_OP == OP_BILATERAL_MEAN || RANK_HIST_OP == OP_BILATERAL_POP || RANK_HIST_OP == OP_BILATERAL_SUM + if (tx == 0) + { + int start_bin = max(0, (int)floor((double)center - s1) + 1); + int stop_bin = min(256, (int)ceil((double)center + s0)); + if (stop_bin <= start_bin) + { + range_start = 0; + range_end = 0; + range_start_sum = 0; + range_end_sum = 0; + } + else + { + range_start = (start_bin > 0) ? scan[start_bin - 1] : 0; + range_end = scan[stop_bin - 1]; + range_start_sum = (start_bin > 0) ? tmp1[start_bin - 1] : 0; + range_end_sum = tmp1[stop_bin - 1]; + } + } + __syncthreads(); +# else + if (tx < 256 && hist[tx] > 0) + { + int bin_end = scan[tx]; + int bin_start = bin_end - hist[tx]; + int weighted_end = tmp1[tx]; + int weighted_start = weighted_end - hist[tx] * tx; + + if (range_start > 0 && bin_start < range_start && bin_end >= range_start) + { + range_start_sum = weighted_start + (range_start - bin_start) * tx; + } + if (range_end > 0 && bin_start < range_end && bin_end >= range_end) + { + range_end_sum = weighted_start + (range_end - bin_start) * tx; + } + } + __syncthreads(); +# endif + + int selected_count_total = range_end - range_start; + int selected_sum_total = range_end_sum - range_start_sum; + if (op == OP_BILATERAL_POP) + { + return static_cast(selected_count_total); + } + if (selected_count_total <= 0) + { + return static_cast(0); + } + if (op == OP_BILATERAL_MEAN) + { + return static_cast(((double)selected_sum_total) / selected_count_total); + } + if (op == OP_BILATERAL_SUM) + { + return static_cast(selected_sum_total); + } + if (op == OP_MEAN) + { + return static_cast(((double)selected_sum_total) / selected_count_total); + } + if (op == OP_SUBTRACT_MEAN) + { + double mean = ((double)selected_sum_total) / selected_count_total; + return static_cast(((double)center - mean) * 0.5 + floor((dtype_max + 1.0) / 2.0)); + } + return static_cast(selected_sum_total); +# endif + + int selected_count = 0; + int selected_sum = 0; + if (tx < 256) + { + int bin_start = scan[tx] - hist[tx]; + int bin_end = scan[tx]; + selected_count = min(bin_end, range_end) - max(bin_start, range_start); + if (selected_count < 0) + { + selected_count = 0; + } + selected_sum = selected_count * tx; + } + + if (op == OP_POP) + { + double low = p0 * pop / 100.0; + double high = p1 * pop / 100.0; + int count = 0; + if (tx < 256 && hist[tx] > 0 && (double)scan[tx] >= low && (double)scan[tx] <= high) + { + count = hist[tx]; + } + tmp0[tx] = count; + __syncthreads(); + reduceSum256(tmp0); + return static_cast(tmp0[0]); + } + + if (op == OP_GRADIENT || op == OP_AUTOLEVEL || op == OP_ENHANCE_CONTRAST) + { + tmp0[tx] = selected_count > 0 ? tx : 255; + tmp1[tx] = selected_count > 0 ? tx : 0; + __syncthreads(); + for (int stride = 128; stride > 0; stride >>= 1) + { + if (tx < stride) + { + tmp0[tx] = min(tmp0[tx], tmp0[tx + stride]); + tmp1[tx] = max(tmp1[tx], tmp1[tx + stride]); + } + __syncthreads(); + } + if (op == OP_GRADIENT) + { + return static_cast(tmp1[0] - tmp0[0]); + } + + int min_val = tmp0[0]; + int max_val = tmp1[0]; + if (op == OP_ENHANCE_CONTRAST) + { + return (max_val - center < center - min_val) ? static_cast(max_val) : + static_cast(min_val); + } + + int clamped = min(max((int)center, min_val), max_val); + int delta = max_val - min_val; + if (delta > 0) + { + return static_cast(((double)(clamped - min_val) / delta) * dtype_max); + } + return static_cast(0); + } + + // clang-format off +#if RANK_HIST_OP != OP_MEAN && RANK_HIST_OP != OP_SUM && RANK_HIST_OP != OP_SUBTRACT_MEAN && RANK_HIST_OP != OP_BILATERAL_MEAN && RANK_HIST_OP != OP_BILATERAL_POP && RANK_HIST_OP != OP_BILATERAL_SUM + // clang-format on + tmp0[tx] = selected_count; + tmp1[tx] = selected_sum; + __syncthreads(); + reduceSum256(tmp0); + reduceSum256(tmp1); + + if (op == OP_MEAN) + { + return static_cast(((double)tmp1[0]) / tmp0[0]); + } + if (op == OP_SUBTRACT_MEAN) + { + double mean = ((double)tmp1[0]) / tmp0[0]; + return static_cast(((double)center - mean) * 0.5 + floor((dtype_max + 1.0) / 2.0)); + } + return static_cast(tmp1[0]); +# endif +#endif +} + +extern "C" __global__ void cuRankHistogram2DUint8(const unsigned char* src, + RANK_HIST_OUTPUT_T* dest, + HIST_COUNTER_T* histPar, + int r0, + int r1, + double p0, + double p1, + double s0, + double s1, + double dtype_max, + int op, + int window_size, + int rows, + int cols) +{ + __shared__ int H[256]; + __shared__ int Hscan[256]; + __shared__ int tmp0[256]; + __shared__ int tmp1[256]; + __shared__ double dtmp[256]; + + int tx = threadIdx.x; + int out_rows = rows - 2 * r0; + int rows_per_block = (out_rows + gridDim.x - 1) / gridDim.x; + int start_out = blockIdx.x * rows_per_block; + int stop_out = min(out_rows, start_out + rows_per_block); + + if (start_out >= stop_out) + { + return; + } + + int start_row = r0 + start_out; + int stop_row = r0 + stop_out; + HIST_COUNTER_T* hist = histPar + blockIdx.x * cols * 256; + + for (int col = tx; col < cols; col += blockDim.x) + { + HIST_COUNTER_T* col_hist = hist + col * 256; + for (int row = start_row - r0; row <= start_row + r0; row++) + { + col_hist[src[row * cols + col]]++; + } + } + __syncthreads(); + + for (int row = start_row; row < stop_row; row++) + { + if (tx < 256) + { + int total = 0; + for (int col = 0; col <= 2 * r1; col++) + { + total += (int)hist[col * 256 + tx]; + } + H[tx] = total; + } + __syncthreads(); + + for (int col = r1; col < cols - r1; col++) + { + unsigned char center = src[row * cols + col]; + RANK_HIST_OUTPUT_T value = + histogramRankValue(H, Hscan, tmp0, tmp1, dtmp, op, window_size, p0, p1, s0, s1, dtype_max, center); + + if (tx == 0) + { + dest[row * cols + col] = value; + } + __syncthreads(); + + if (col < cols - r1 - 1 && tx < 256) + { + int sub_col = col - r1; + int add_col = col + r1 + 1; + H[tx] += (int)hist[add_col * 256 + tx] - (int)hist[sub_col * 256 + tx]; + } + __syncthreads(); + } + + if (row < stop_row - 1) + { + int sub_row = row - r0; + int add_row = row + r0 + 1; + for (int col = tx; col < cols; col += blockDim.x) + { + HIST_COUNTER_T* col_hist = hist + col * 256; + col_hist[src[sub_row * cols + col]]--; + col_hist[src[add_row * cols + col]]++; + } + __syncthreads(); + } + } +} diff --git a/python/cucim/src/cucim/skimage/filters/rank/tests/test_rank.py b/python/cucim/src/cucim/skimage/filters/rank/tests/test_rank.py new file mode 100644 index 000000000..52ee82e97 --- /dev/null +++ b/python/cucim/src/cucim/skimage/filters/rank/tests/test_rank.py @@ -0,0 +1,1913 @@ +# SPDX-FileCopyrightText: 2009-2022 the scikit-image team +# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 AND BSD-3-Clause + +import warnings + +import cupy as cp +import numpy as np +import pytest +from skimage import data +from skimage._shared.testing import fetch + +from cucim.skimage import morphology, util +from cucim.skimage._shared.testing import expected_warnings +from cucim.skimage.filters import rank +from cucim.skimage.filters.rank import ( + __all__ as all_rank_filters, + subtract_mean, +) +from cucim.skimage.filters.rank._histogram import ( + _get_histogram_counter_dtype, + _get_rank_histogram_partitions, + _should_use_rank_histogram, +) +from cucim.skimage.morphology import ball, disk, gray +from cucim.skimage.util import img_as_float, img_as_ubyte + + +def _reflect_index(index, size): + if index < 0: + index = -1 - index + index %= 2 * size + return min(index, 2 * size - 1 - index) + + +def _cast_uint8(value): + return np.uint8(int(value) % 256) + + +def _rank_filter_brute_force_uint8( + image, + footprint, + operation, + *, + mask=None, + s0=10, + s1=10, +): + radius = tuple(s // 2 for s in footprint.shape) + out = np.empty_like(image) + dtype_max = 255 + for row in range(image.shape[0]): + for col in range(image.shape[1]): + values = [] + for frow in range(footprint.shape[0]): + for fcol in range(footprint.shape[1]): + if not footprint[frow, fcol]: + continue + if ( + operation == "noise_filter" + and ( + frow, + fcol, + ) + == radius + ): + continue + irow = _reflect_index( + row + frow - radius[0], image.shape[0] + ) + icol = _reflect_index( + col + fcol - radius[1], image.shape[1] + ) + if mask is not None and not mask[irow, icol]: + continue + values.append(int(image[irow, icol])) + + center = int(image[row, col]) + if not values: + out[row, col] = image[row, col] + elif operation == "minimum": + out[row, col] = min(values) + elif operation == "maximum": + out[row, col] = max(values) + elif operation == "mean": + out[row, col] = _cast_uint8(sum(values) / len(values)) + elif operation == "sum": + out[row, col] = _cast_uint8(sum(values)) + elif operation == "subtract_mean": + mean = sum(values) / len(values) + out[row, col] = _cast_uint8( + _cast_uint8((center - mean) * 0.5 + 128) - 1 + ) + elif operation == "pop": + out[row, col] = _cast_uint8(len(values)) + elif operation == "threshold": + out[row, col] = _cast_uint8( + center > (sum(values) / len(values)) + ) + elif operation == "gradient": + out[row, col] = _cast_uint8(max(values) - min(values)) + elif operation == "entropy": + counts = np.bincount(values, minlength=256) + probabilities = counts[counts > 0] / len(values) + out[row, col] = _cast_uint8( + -(probabilities * np.log2(probabilities)).sum() + ) + elif operation == "autolevel": + min_val = min(values) + max_val = max(values) + delta = max_val - min_val + if delta > 0: + clamped = min(max(center, min_val), max_val) + out[row, col] = _cast_uint8( + (clamped - min_val) / delta * dtype_max + ) + else: + out[row, col] = 0 + elif operation == "enhance_contrast": + min_val = min(values) + max_val = max(values) + if max_val - center < center - min_val: + out[row, col] = max_val + else: + out[row, col] = min_val + elif operation == "equalize": + rank = sum(value <= center for value in values) + out[row, col] = _cast_uint8(dtype_max * rank / len(values)) + elif operation == "geometric_mean": + log_sum = sum(np.log(value + 1.0) for value in values) + out[row, col] = _cast_uint8( + round(np.exp(log_sum / len(values)) - 1.0) + ) + elif operation in {"modal", "majority"}: + counts = np.bincount(values, minlength=256) + out[row, col] = int(np.argmax(counts)) + elif operation == "noise_filter": + if center in values: + out[row, col] = 0 + else: + out[row, col] = min(abs(value - center) for value in values) + elif operation in { + "mean_bilateral", + "pop_bilateral", + "sum_bilateral", + }: + bilateral_values = [ + value + for value in values + if center > value - s0 and center < value + s1 + ] + if operation == "pop_bilateral": + out[row, col] = _cast_uint8(len(bilateral_values)) + elif not bilateral_values: + out[row, col] = 0 + elif operation == "mean_bilateral": + out[row, col] = _cast_uint8( + sum(bilateral_values) / len(bilateral_values) + ) + else: + out[row, col] = _cast_uint8(sum(bilateral_values)) + else: + raise ValueError(f"unsupported operation: {operation}") + return out + + +def _rank_percentile_brute_force_uint8( + image, + footprint, + operation, + *, + p0=0, + p1=1, +): + radius = tuple(s // 2 for s in footprint.shape) + out = np.empty_like(image) + for row in range(image.shape[0]): + for col in range(image.shape[1]): + values = [] + for frow in range(footprint.shape[0]): + for fcol in range(footprint.shape[1]): + if not footprint[frow, fcol]: + continue + irow = _reflect_index( + row + frow - radius[0], image.shape[0] + ) + icol = _reflect_index( + col + fcol - radius[1], image.shape[1] + ) + values.append(int(image[irow, icol])) + values.sort() + center = int(image[row, col]) + pop = len(values) + + if operation == "percentile": + if p0 == 1: + percentile_idx = pop - 1 + else: + percentile_idx = int(p0 * pop) + if percentile_idx >= pop: + percentile_idx = pop - 1 + out[row, col] = values[percentile_idx] + continue + + if operation == "threshold_percentile": + threshold_idx = int(p0 * pop) + if threshold_idx >= pop: + threshold_idx = pop - 1 + out[row, col] = 255 if center >= values[threshold_idx] else 0 + continue + + idx_start = max(0, int(np.ceil(p0 * pop)) - 1) + idx_end = int(p1 * pop) + if idx_end <= idx_start: + idx_end = idx_start + 1 + if idx_end > pop: + idx_end = pop + selected = values[idx_start:idx_end] + + if operation == "mean_percentile": + out[row, col] = _cast_uint8(sum(selected) / len(selected)) + elif operation == "sum_percentile": + out[row, col] = _cast_uint8(sum(selected)) + elif operation == "gradient_percentile": + out[row, col] = _cast_uint8(selected[-1] - selected[0]) + elif operation == "subtract_mean_percentile": + mean = sum(selected) / len(selected) + out[row, col] = _cast_uint8((center - mean) * 0.5 + 128) + elif operation == "enhance_contrast_percentile": + min_val = selected[0] + max_val = selected[-1] + if max_val - center < center - min_val: + out[row, col] = max_val + else: + out[row, col] = min_val + elif operation == "autolevel_percentile": + min_val = selected[0] + max_val = selected[-1] + delta = max_val - min_val + if delta > 0: + clamped = min(max(center, min_val), max_val) + out[row, col] = _cast_uint8( + (clamped - min_val) / delta * 255 + ) + else: + out[row, col] = 0 + elif operation == "pop_percentile": + count = 0 + cumsum = 0 + i = 0 + while i < pop: + value = values[i] + group_size = 0 + while i < pop and values[i] == value: + group_size += 1 + i += 1 + cumsum += group_size + if cumsum >= p0 * pop and cumsum <= p1 * pop: + count += group_size + out[row, col] = _cast_uint8(count) + else: + raise ValueError(f"unsupported operation: {operation}") + return out + + +@pytest.mark.parametrize("dtype", [np.uint8, np.uint16]) +def test_subtract_mean_underflow_correction(dtype): + # Input: [10, 10, 10] + footprint = cp.ones((1, 3)) + arr = cp.array([[10, 10, 10]], dtype=dtype) + result = subtract_mean(arr, footprint, cast_to_uint8=(dtype == np.uint8)) + + if dtype == np.uint8: + expected_val = 127 + else: + expected_val = 32767 + # note: scikit-image expected_val for uint16 was + # expected_val = (arr.max() + 1) // 2 - 1 + + assert cp.all(result == expected_val) + + +@pytest.mark.parametrize( + "filter_name", + [ + "minimum", + "maximum", + "mean", + "sum", + "subtract_mean", + "pop", + "threshold", + "gradient", + "autolevel", + "enhance_contrast", + "equalize", + "geometric_mean", + "noise_filter", + "mean_bilateral", + "pop_bilateral", + "sum_bilateral", + ], +) +@pytest.mark.parametrize("use_mask", [False, True]) +def test_streaming_rank_filter_ops_uint8(filter_name, use_mask): + image = np.array( + [ + [0, 5, 10, 50, 90], + [3, 8, 20, 60, 120], + [7, 11, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + footprint = np.array( + [ + [1, 0, 1], + [1, 1, 0], + [0, 1, 1], + ], + dtype=bool, + ) + mask = np.array( + [ + [1, 1, 0, 1, 1], + [1, 0, 1, 1, 0], + [0, 1, 1, 0, 1], + [1, 1, 1, 1, 1], + ], + dtype=bool, + ) + kwargs = {} + if "bilateral" in filter_name: + kwargs.update(s0=6, s1=9) + + result = getattr(rank, filter_name)( + cp.asarray(image), + cp.asarray(footprint), + mask=cp.asarray(mask) if use_mask else None, + **kwargs, + ) + expected = _rank_filter_brute_force_uint8( + image, + footprint, + filter_name, + mask=mask if use_mask else None, + **kwargs, + ) + cp.testing.assert_array_equal(result, cp.asarray(expected)) + + +@pytest.mark.parametrize( + "shift_kwargs", + [{}, {"shift_x": 1}, {"shifts": (-1, 0)}], +) +def test_noise_filter_excludes_shifted_anchor(shift_kwargs): + image = cp.zeros((7, 7), dtype=cp.uint8) + image[3, 3] = 10 + footprint = cp.ones((3, 3), dtype=bool) + footprint_before = footprint.copy() + + result = rank.noise_filter( + image, + footprint, + backend="elementwise", + **shift_kwargs, + ) + + assert result[3, 3] == 10 + cp.testing.assert_array_equal(footprint, footprint_before) + + +@pytest.mark.parametrize( + "dtype, values, expected", + [ + (cp.float32, [0.5, 0.6, 0.8], 0.1), + (cp.uint64, [0, 2**32, 2**33], 2**32), + ], +) +def test_noise_filter_preserves_native_distance_precision( + dtype, values, expected +): + image = cp.asarray([values], dtype=dtype) + footprint = cp.ones((1, 3), dtype=bool) + + result = rank.noise_filter( + image, + footprint, + backend="elementwise", + cast_to_uint8=False, + ) + + if dtype == cp.float32: + cp.testing.assert_allclose(result[0, 1], expected, rtol=1e-6) + else: + cp.testing.assert_array_equal(result[0, 1], expected) + + +@pytest.mark.parametrize( + "filter_name, kwargs", + [ + ("percentile", dict(p0=0.5)), + ("percentile", dict(p0=1.0)), + ("threshold_percentile", dict(p0=0.5)), + ("mean_percentile", dict(p0=0.25, p1=0.75)), + ("sum_percentile", dict(p0=0.25, p1=0.75)), + ("pop_percentile", dict(p0=0.25, p1=0.75)), + ("gradient_percentile", dict(p0=0.25, p1=0.75)), + ("autolevel_percentile", dict(p0=0.25, p1=0.75)), + ("enhance_contrast_percentile", dict(p0=0.25, p1=0.75)), + ("subtract_mean_percentile", dict(p0=0.25, p1=0.75)), + ], +) +def test_histogram_rank_percentile_ops_uint8_rectangular(filter_name, kwargs): + image = np.array( + [ + [0, 5, 10, 50, 90], + [3, 8, 20, 60, 120], + [7, 11, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + footprint = np.ones((3, 3), dtype=bool) + result = getattr(rank, filter_name)( + cp.asarray(image), + cp.asarray(footprint), + backend="histogram", + **kwargs, + ) + expected = _rank_percentile_brute_force_uint8( + image, + footprint, + filter_name, + **kwargs, + ) + cp.testing.assert_array_equal(result, cp.asarray(expected)) + + +def test_histogram_rank_entropy_uint8_rectangular(): + image = np.array( + [ + [0, 0, 10, 50, 90], + [0, 8, 20, 60, 120], + [7, 8, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + footprint = np.ones((3, 3), dtype=bool) + result = rank.entropy( + cp.asarray(image), cp.asarray(footprint), backend="histogram" + ) + expected = rank.entropy( + cp.asarray(image), cp.asarray(footprint), backend="elementwise" + ) + assert result.dtype.kind == "f" + cp.testing.assert_allclose(result, expected) + + +@pytest.mark.parametrize( + "filter_name, kwargs", + [ + ("equalize", {}), + ("geometric_mean", {}), + ("mean_bilateral", dict(s0=6, s1=9)), + ("modal", {}), + ("majority", {}), + ("pop_bilateral", dict(s0=6, s1=9)), + ("sum_bilateral", dict(s0=6, s1=9)), + ], +) +def test_histogram_rank_prefix_ops_uint8_rectangular(filter_name, kwargs): + image = np.array( + [ + [0, 5, 10, 50, 90], + [3, 8, 20, 60, 120], + [7, 11, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + footprint = np.ones((3, 3), dtype=bool) + result = getattr(rank, filter_name)( + cp.asarray(image), + cp.asarray(footprint), + backend="histogram", + **kwargs, + ) + expected = _rank_filter_brute_force_uint8( + image, + footprint, + filter_name, + **kwargs, + ) + cp.testing.assert_array_equal(result, cp.asarray(expected)) + + +def test_rank_backend_override_histogram_and_elementwise(): + image = cp.asarray( + np.array( + [ + [0, 5, 10, 50, 90], + [3, 8, 20, 60, 120], + [7, 11, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + ) + footprint = cp.ones((3, 3), dtype=bool) + + automatic = rank.percentile(image, footprint, p0=0.5, backend="auto") + histogram = rank.percentile(image, footprint, p0=0.5, backend="histogram") + elementwise = rank.percentile( + image, footprint, p0=0.5, backend="elementwise" + ) + + cp.testing.assert_array_equal(histogram, automatic) + cp.testing.assert_array_equal(elementwise, automatic) + + +def test_rank_backend_auto_uses_elementwise_below_histogram_cutoff(): + image = cp.asarray( + np.array( + [ + [0, 5, 10, 50, 90], + [3, 8, 20, 60, 120], + [7, 11, 30, 70, 150], + [13, 17, 40, 80, 180], + ], + dtype=np.uint8, + ) + ) + footprint = cp.ones((3, 3), dtype=bool) + + automatic = rank.percentile(image, footprint, p0=0.5, backend="auto") + elementwise = rank.percentile( + image, footprint, p0=0.5, backend="elementwise" + ) + histogram = rank.percentile(image, footprint, p0=0.5, backend="histogram") + + cp.testing.assert_array_equal(automatic, elementwise) + cp.testing.assert_array_equal(histogram, elementwise) + + +def test_rank_backend_histogram_rejects_incompatible_input(): + image = cp.asarray(np.arange(25, dtype=np.uint16).reshape(5, 5)) + footprint = cp.ones((3, 3), dtype=bool) + + with pytest.raises(ValueError, match="backend='histogram' requires"): + rank.percentile( + image, + footprint, + p0=0.5, + backend="histogram", + cast_to_uint8=False, + ) + + +@pytest.mark.parametrize("out_dtype", [cp.float32, cp.uint16]) +def test_rank_backend_histogram_supports_non_uint8_output(out_dtype): + image = cp.asarray(np.arange(32 * 32, dtype=np.uint8).reshape(32, 32)) + footprint = cp.ones((17, 17), dtype=bool) + out = cp.empty(image.shape, dtype=out_dtype) + + result = rank.percentile( + image, footprint, p0=0.5, out=out, backend="histogram" + ) + expected = rank.percentile( + image, footprint, p0=0.5, backend="histogram" + ).astype(out_dtype) + + assert result is out + assert result.dtype == out_dtype + cp.testing.assert_array_equal(result, expected) + + auto = rank.percentile( + image, footprint, p0=0.5, out=cp.empty_like(out), backend="auto" + ) + cp.testing.assert_array_equal(auto, expected) + + +@pytest.mark.parametrize( + "filter_name, kwargs", + [ + ("sum_percentile", dict(p0=0.01, p1=1.0)), + ("pop_percentile", dict(p0=0.01, p1=1.0)), + ("sum_bilateral", dict(s0=1, s1=1)), + ], +) +def test_histogram_rank_preserves_uint16_results(filter_name, kwargs): + image = cp.full((32, 32), 255, dtype=cp.uint8) + footprint = cp.ones((29, 29), dtype=bool) + + histogram = getattr(rank, filter_name)( + image, + footprint, + out=cp.empty(image.shape, dtype=cp.uint16), + backend="histogram", + **kwargs, + ) + elementwise = getattr(rank, filter_name)( + image, + footprint, + out=cp.empty(image.shape, dtype=cp.uint16), + backend="elementwise", + **kwargs, + ) + automatic = getattr(rank, filter_name)( + image, + footprint, + out=cp.empty(image.shape, dtype=cp.uint16), + backend="auto", + **kwargs, + ) + + cp.testing.assert_array_equal(histogram, elementwise) + cp.testing.assert_array_equal(automatic, elementwise) + + +@pytest.mark.parametrize("out_dtype", [cp.float32, cp.uint16]) +@pytest.mark.parametrize( + "filter_name, kwargs", + [ + ("threshold_percentile", dict(p0=0.5)), + ("equalize", {}), + ("autolevel_percentile", dict(p0=0.1, p1=0.9)), + ("subtract_mean_percentile", dict(p0=0.1, p1=0.9)), + ], +) +def test_histogram_rank_uses_output_dtype_scale(filter_name, kwargs, out_dtype): + image = cp.arange(24 * 24, dtype=cp.uint8).reshape(24, 24) + footprint = cp.ones((17, 17), dtype=bool) + + histogram = getattr(rank, filter_name)( + image, + footprint, + out=cp.empty(image.shape, dtype=out_dtype), + backend="histogram", + **kwargs, + ) + elementwise = getattr(rank, filter_name)( + image, + footprint, + out=cp.empty(image.shape, dtype=out_dtype), + backend="elementwise", + **kwargs, + ) + + cp.testing.assert_allclose(histogram, elementwise) + + +def test_histogram_rank_rejects_unsupported_output_dtype(): + image = cp.arange(24 * 24, dtype=cp.uint8).reshape(24, 24) + footprint = cp.ones((17, 17), dtype=bool) + out = cp.empty(image.shape, dtype=cp.uint32) + + with pytest.raises(ValueError, match="backend='histogram' requires"): + rank.sum_percentile( + image, + footprint, + p0=0.01, + p1=1.0, + out=out, + backend="histogram", + ) + + automatic = rank.sum_percentile( + image, + footprint, + p0=0.01, + p1=1.0, + out=out, + backend="auto", + ) + elementwise = rank.sum_percentile( + image, + footprint, + p0=0.01, + p1=1.0, + out=cp.empty_like(out), + backend="elementwise", + ) + cp.testing.assert_array_equal(automatic, elementwise) + + +def test_rank_backend_invalid_value_raises(): + image = cp.asarray(np.arange(25, dtype=np.uint8).reshape(5, 5)) + footprint = cp.ones((3, 3), dtype=bool) + + with pytest.raises(ValueError, match="backend must be one of"): + rank.percentile(image, footprint, p0=0.5, backend="bad") + + +def test_rank_requires_cupy_inputs(): + image = cp.asarray(np.arange(25, dtype=np.uint8).reshape(5, 5)) + footprint = cp.ones((3, 3), dtype=bool) + mask = cp.ones_like(image, dtype=bool) + + with pytest.raises(ValueError, match="image must be a CuPy array"): + rank.percentile(cp.asnumpy(image), footprint, p0=0.5) + + with pytest.raises(ValueError, match="footprint must be a CuPy array"): + rank.percentile(image, cp.asnumpy(footprint), p0=0.5) + + with pytest.raises(ValueError, match="mask must be a CuPy array"): + rank.percentile(image, footprint, mask=cp.asnumpy(mask), p0=0.5) + + +@pytest.mark.parametrize("alias", ["view", "transpose"]) +def test_rank_rejects_overlapping_output_alias(alias): + image = cp.arange(25, dtype=cp.uint8).reshape(5, 5) + footprint = cp.ones((3, 3), dtype=bool) + out = image.view() if alias == "view" else image.T + + with pytest.raises( + NotImplementedError, match="Cannot perform rank operation in place" + ): + rank.percentile( + image, + footprint, + p0=0.5, + out=out, + backend="elementwise", + ) + + +def test_rank_median_default_footprint(): + image = cp.asarray(np.arange(25, dtype=np.uint8).reshape(5, 5)) + expected = rank.median(image, cp.ones((3, 3), dtype=bool)) + result = rank.median(image) + + cp.testing.assert_array_equal(result, expected) + + +def test_rank_default_cast_to_uint8_matches_explicit_float_conversion(): + image = cp.linspace(0, 1, 25, dtype=cp.float32).reshape(5, 5) + footprint = cp.ones((3, 3), dtype=bool) + image_u8 = img_as_ubyte(image) + + expected = rank.percentile( + image_u8, footprint, p0=0.5, backend="elementwise" + ) + with expected_warnings(["Possible precision loss"]): + result = rank.percentile( + image, + footprint, + p0=0.5, + backend="elementwise", + ) + + assert result.dtype == cp.uint8 + cp.testing.assert_array_equal(result, expected) + + +def test_rank_default_cast_to_uint8_matches_explicit_uint16_conversion(): + image = cp.linspace(0, 65535, 25, dtype=cp.uint16).reshape(5, 5) + footprint = cp.ones((3, 3), dtype=bool) + image_u8 = img_as_ubyte(image) + + expected = rank.percentile( + image_u8, footprint, p0=0.5, backend="elementwise" + ) + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always") + result = rank.percentile( + image, + footprint, + p0=0.5, + backend="elementwise", + ) + assert not record + + assert result.dtype == cp.uint8 + cp.testing.assert_array_equal(result, expected) + + +def test_rank_cast_to_uint8_before_histogram_backend_selection(): + image = cp.linspace(0, 1, 25 * 25, dtype=cp.float32).reshape(25, 25) + footprint = cp.ones((17, 17), dtype=bool) + image_u8 = img_as_ubyte(image) + + with pytest.raises(ValueError, match="backend='histogram' requires"): + rank.percentile( + image, + footprint, + p0=0.5, + backend="histogram", + cast_to_uint8=False, + ) + + expected = rank.percentile(image_u8, footprint, p0=0.5, backend="histogram") + with expected_warnings(["Possible precision loss"]): + result = rank.percentile( + image, + footprint, + p0=0.5, + backend="histogram", + ) + + assert result.dtype == cp.uint8 + cp.testing.assert_array_equal(result, expected) + + +def test_rank_uint16_elementwise_does_not_warn_about_bins(): + image = cp.asarray(np.arange(25, dtype=np.uint16).reshape(5, 5)) + image[-1, -1] = 2048 + footprint = cp.ones((3, 3), dtype=bool) + + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always") + rank.percentile( + image, + footprint, + p0=0.5, + backend="elementwise", + cast_to_uint8=False, + ) + assert not record + + with pytest.raises(ValueError, match="backend='histogram' requires"): + rank.percentile( + image, + footprint, + p0=0.5, + backend="histogram", + cast_to_uint8=False, + ) + + +def test_rank_histogram_partitions_default_and_env(monkeypatch): + monkeypatch.delenv("CUCIM_RANK_HISTOGRAM_PARTITIONS", raising=False) + monkeypatch.delenv("CUCIM_RANK_HISTOGRAM_SCRATCH_MB", raising=False) + monkeypatch.delenv("CUCIM_RANK_HISTOGRAM_MAX_PARTITIONS", raising=False) + + assert ( + _get_rank_histogram_partitions(1080, 1080, counter_dtype=cp.int32) + == 242 + ) + assert ( + _get_rank_histogram_partitions(1080, 1080, counter_dtype=cp.int16) + == 256 + ) + + monkeypatch.setenv("CUCIM_RANK_HISTOGRAM_MAX_PARTITIONS", "64") + assert ( + _get_rank_histogram_partitions(1080, 1080, counter_dtype=cp.int32) == 64 + ) + + monkeypatch.setenv("CUCIM_RANK_HISTOGRAM_PARTITIONS", "32") + assert ( + _get_rank_histogram_partitions(1080, 1080, counter_dtype=cp.int32) == 32 + ) + + +def test_rank_histogram_counter_dtype(): + assert _get_histogram_counter_dtype((181, 181)) == cp.int16 + assert _get_histogram_counter_dtype((181, 183)) == cp.int32 + + +def test_rank_histogram_auto_cutoffs(): + assert not _should_use_rank_histogram("percentile", (15, 15)) + assert _should_use_rank_histogram("percentile", (17, 17)) + assert not _should_use_rank_histogram("mean", (17, 17)) + assert _should_use_rank_histogram("mean", (19, 19)) + assert not _should_use_rank_histogram("entropy", (23, 23)) + assert _should_use_rank_histogram("entropy", (25, 25)) + assert not _should_use_rank_histogram("bilateral_mean", (31, 31)) + assert _should_use_rank_histogram("bilateral_mean", (33, 33)) + assert not _should_use_rank_histogram("geometric_mean", (13, 13)) + assert _should_use_rank_histogram("geometric_mean", (15, 15)) + assert not _should_use_rank_histogram("modal", (13, 13)) + assert _should_use_rank_histogram("modal", (15, 15)) + assert not _should_use_rank_histogram("equalize", (71, 71)) + assert _should_use_rank_histogram("equalize", (91, 91)) + + +# # Note: Explicitly read all values into a dict. Otherwise, stochastic test +# # failures related to I/O can occur during parallel test cases. +ref_data = dict(np.load(fetch("data/rank_filter_tests.npz"))) +ref_data_3d = dict(np.load(fetch("data/rank_filters_tests_3d.npz"))) + + +class TestRank: + def setup_method(self): + np.random.seed(0) + # This image is used along with @run_in_parallel + # to ensure that the same seed is used for each thread. + self.image = cp.asarray(np.random.rand(25, 25)) + np.random.seed(0) + self.volume = cp.asarray(np.random.rand(10, 10, 10)) + # Set again the seed for the other tests. + np.random.seed(0) + self.footprint = morphology.disk(1) + self.footprint_3d = morphology.ball(1) + self.refs = ref_data + self.refs_3d = ref_data_3d + + # Filters where the only differences vs scikit-image are at image + # borders (due to reflected boundary extension vs excluded pixels). + # For these, we compare only interior pixels. + _border_differences_allowed = { + "entropy", + "equalize", + "geometric_mean", + "majority", + "mean", + "mean_bilateral", + "mean_percentile", + "median", + "modal", + "pop", + "pop_bilateral", + "pop_percentile", + "subtract_mean", + "subtract_mean_percentile", + "sum", + "sum_bilateral", + "sum_percentile", + } + + # Filters with known algorithmic differences that are documented and + # expected. These are tested separately or skipped here. + _xfail_filters = { + # gradient_percentile: scikit-image's histogram p1-inversion quirk + # makes imax=255 always; our sorted-array computes correct max-min. + "gradient_percentile", + # noise_filter: scikit-image treats the dtype endpoints as implicit + # neighbors when all actual neighbors lie on one side of the center. + # Our implementation returns the nearest actual neighbor distance. + "noise_filter", + } + + @pytest.mark.parametrize("filter", all_rank_filters) + def test_rank_filter(self, filter): + """Test rank filters with uint8 input against scikit-image reference. + + The reference data in rank_filter_tests.npz was generated by + scikit-image which internally converts float images to uint8. We + keep the same default conversion behavior for closer compatibility. + """ + if filter in self._xfail_filters: + pytest.skip( + f"{filter}: known algorithmic difference vs scikit-image (not a bug)" + ) + expected = cp.asarray(self.refs[filter]) + with expected_warnings(["Possible precision loss"]): + result = getattr(rank, filter)( + self.image, self.footprint, cast_to_uint8=True + ) + if filter in self._border_differences_allowed: + # Only compare interior pixels — borders differ due to + # reflected boundary extension (GPU) vs excluded pixels + # (scikit-image). + expected = expected[1:-1, 1:-1] + result = result[1:-1, 1:-1] + if filter == "subtract_mean_percentile": + # Allow off-by-1 due to documented mid_bin offset difference: + # percentile variant uses (dtype_max + 1) / 2 = 128 for uint8, + # scikit-image's percentile variant uses the same, but the + # histogram-based integer arithmetic can round differently, + # giving a systematic -1 offset on some pixels. + cp.testing.assert_allclose(expected, result, atol=1) + else: + cp.testing.assert_allclose(expected, result) + + @pytest.mark.parametrize("filter", all_rank_filters) + def test_rank_filter_footprint_sequence_unsupported(self, filter): + footprint_sequence = morphology.diamond(3, decomposition="sequence") + with pytest.raises(ValueError): + getattr(rank, filter)( + self.image.astype(np.uint8), footprint_sequence + ) + + @pytest.mark.parametrize("outdt", [None]) # , cp.float32, cp.float64]) + @pytest.mark.parametrize( + "filter", + [ + "autolevel", + "equalize", + "gradient", + "majority", + "maximum", + "mean", + "geometric_mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "sum", + "threshold", + "noise_filter", + "entropy", + ], + ) + def test_rank_filters_3D(self, filter, outdt): + if filter in self._xfail_filters: + pytest.skip( + f"{filter}: known algorithmic difference vs scikit-image (not a bug)" + ) + expected = cp.asarray(self.refs_3d[filter]) + if outdt is not None: + out = cp.zeros_like(expected, dtype=outdt) + else: + out = None + with expected_warnings(["Possible precision loss"]): + result = getattr(rank, filter)( + self.volume, self.footprint_3d, out=out, cast_to_uint8=True + ) + if outdt is not None: + # Avoid rounding issues comparing to expected result + if filter == "sum": + # sum test data seems to be 8-bit disguised as 16-bit + datadt = cp.uint8 + else: + datadt = expected.dtype + # Take modulus first to avoid undefined behavior for + # float->uint8 conversions. + result = cp.mod(result, 256.0).astype(datadt) + if filter in self._border_differences_allowed: + # Only compare interior pixels — borders differ due to + # reflected boundary extension (GPU) vs excluded pixels + # (scikit-image). + expected = expected[1:-1, 1:-1, 1:-1] + result = result[1:-1, 1:-1, 1:-1] + if filter == "subtract_mean_percentile": + # Allow off-by-1 due to documented mid_bin offset difference: + # percentile variant uses (dtype_max + 1) / 2 = 128 for uint8, + # scikit-image's percentile variant uses the same, but the + # histogram-based integer arithmetic can round differently, + # giving a systematic -1 offset on some pixels. + cp.testing.assert_allclose(expected, result, atol=1) + cp.testing.assert_array_almost_equal(expected, result) + + def test_random_sizes(self): + # make sure the size is not a problem + elem = cp.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=cp.uint8) + for m, n in np.random.randint(1, 101, size=(10, 2)): + mask = cp.ones((m, n), dtype=cp.uint8) + + image8 = cp.ones((m, n), dtype=cp.uint8) + out8 = cp.empty_like(image8) + rank.mean( + image=image8, + footprint=elem, + mask=mask, + out=out8, + shift_x=0, + shift_y=0, + ) + assert image8.shape == out8.shape + rank.mean( + image=image8, + footprint=elem, + mask=mask, + out=out8, + shift_x=+1, + shift_y=+1, + ) + assert image8.shape == out8.shape + + rank.geometric_mean( + image=image8, + footprint=elem, + mask=mask, + out=out8, + shift_x=0, + shift_y=0, + ) + assert image8.shape == out8.shape + + rank.geometric_mean( + image=image8, + footprint=elem, + mask=mask, + out=out8, + shift_x=+1, + shift_y=+1, + ) + assert image8.shape == out8.shape + + image16 = cp.ones((m, n), dtype=cp.uint16) + out16 = cp.empty_like(image8, dtype=cp.uint16) + rank.mean( + image=image16, + footprint=elem, + mask=mask, + out=out16, + shift_x=0, + shift_y=0, + ) + assert image16.shape == out16.shape + rank.mean( + image=image16, + footprint=elem, + mask=mask, + out=out16, + shift_x=+1, + shift_y=+1, + ) + assert image16.shape == out16.shape + + rank.geometric_mean( + image=image16, + footprint=elem, + mask=mask, + out=out16, + shift_x=0, + shift_y=0, + ) + assert image16.shape == out16.shape + rank.geometric_mean( + image=image16, + footprint=elem, + mask=mask, + out=out16, + shift_x=+1, + shift_y=+1, + ) + assert image16.shape == out16.shape + + rank.mean_percentile( + image=image16, + mask=mask, + out=out16, + footprint=elem, + shift_x=0, + shift_y=0, + p0=0.1, + p1=0.9, + ) + assert image16.shape == out16.shape + rank.mean_percentile( + image=image16, + mask=mask, + out=out16, + footprint=elem, + shift_x=+1, + shift_y=+1, + p0=0.1, + p1=0.9, + ) + assert image16.shape == out16.shape + + @pytest.mark.parametrize("r", list(range(3, 20, 2))) + def test_compare_with_gray_dilation(self, r): + # compare the result of maximum filter with dilate + + image = (cp.random.rand(100, 100) * 256).astype(cp.uint8) + out = cp.empty_like(image) + mask = cp.ones(image.shape, dtype=cp.uint8) + + elem = cp.ones((r, r), dtype=np.uint8) + rank.maximum(image=image, footprint=elem, out=out, mask=mask) + cm = gray.dilation(image, elem) + cp.testing.assert_array_equal(out, cm) + + @pytest.mark.parametrize("r", list(range(3, 20, 2))) + def test_compare_with_gray_erosion(self, r): + # compare the result of maximum filter with erode + + image = (cp.random.rand(100, 100) * 256).astype(cp.uint8) + out = cp.empty_like(image) + mask = cp.ones(image.shape, dtype=cp.uint8) + + elem = cp.ones((r, r), dtype=np.uint8) + rank.minimum(image=image, footprint=elem, out=out, mask=mask) + cm = gray.erosion(image, elem) + cp.testing.assert_array_equal(out, cm) + + def test_population(self): + # check the number of valid pixels in the neighborhood + image = cp.zeros((5, 5), dtype=np.uint8) + elem = cp.ones((3, 3), dtype=np.uint8) + out = cp.empty_like(image) + mask = cp.ones(image.shape, dtype=np.uint8) + + rank.pop(image=image, footprint=elem, out=out, mask=mask) + r = cp.array( + [ + [4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4], + ] + ) + # Note: omit boundaries due to known difference in boundary handling + cp.testing.assert_array_equal(r[1:-1, 1:-1], out[1:-1, 1:-1]) + + def test_structuring_element8(self): + # check the output for a custom footprint + + r = cp.array( + [ + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 255, 0, 0, 0], + [0, 0, 255, 255, 255, 0], + [0, 0, 0, 255, 255, 0], + [0, 0, 0, 0, 0, 0], + ] + ) + + # 8-bit + image = cp.zeros((6, 6), dtype=np.uint8) + image[2, 2] = 255 + elem = cp.asarray([[1, 1, 0], [1, 1, 1], [0, 0, 1]], dtype=np.uint8) + out = cp.empty_like(image) + mask = cp.ones(image.shape, dtype=np.uint8) + + rank.maximum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=1, + shift_y=1, + ) + cp.testing.assert_array_equal(r, out) + + # 16-bit + image = cp.zeros((6, 6), dtype=np.uint16) + image[2, 2] = 255 + out = np.empty_like(image) + + rank.maximum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=1, + shift_y=1, + ) + cp.testing.assert_array_equal(r, out) + + def test_inplace_output(self): + # rank filters are not supposed to filter inplace + + footprint = disk(20, decomposition=None) + image = (cp.random.rand(500, 500) * 256).astype(np.uint8) + out = image + with pytest.raises(NotImplementedError): + rank.mean(image, footprint, out=out) + + def test_compare_autolevels(self): + # compare autolevel and percentile autolevel with p0=0.0 and p1=1.0 + # should returns the same arrays + + image = util.img_as_ubyte(cp.asarray(data.camera())) + + footprint = disk(20, decomposition=None) + loc_autolevel = rank.autolevel(image, footprint=footprint) + loc_perc_autolevel = rank.autolevel_percentile( + image, footprint=footprint, p0=0.0, p1=1.0 + ) + + cp.testing.assert_array_equal(loc_autolevel, loc_perc_autolevel) + + def test_compare_autolevels_16bit(self): + # compare autolevel(16-bit) and percentile autolevel(16-bit) with + # p0=0.0 and p1=1.0 should returns the same arrays + + image = cp.asarray(data.camera()).astype(np.uint16) * 4 + + footprint = disk(20, decomposition=None) + loc_autolevel = rank.autolevel(image, footprint=footprint) + loc_perc_autolevel = rank.autolevel_percentile( + image, footprint=footprint, p0=0.0, p1=1.0 + ) + + cp.testing.assert_array_equal(loc_autolevel, loc_perc_autolevel) + + @pytest.mark.parametrize( + "method", + [ + "autolevel", + "equalize", + "gradient", + "threshold", + "subtract_mean", + "enhance_contrast", + "pop", + ], + ) + def test_compare_ubyte_vs_float(self, method): + # Create signed int8 image that and convert it to uint8 + image_uint = img_as_ubyte(cp.asarray(data.camera()[:50, :50])) + image_float = img_as_float(image_uint) + + disk3 = disk(3, decomposition=None) + func = getattr(rank, method) + out_u = func(image_uint, disk3) + with expected_warnings(["Possible precision loss"]): + out_f = func(image_float, disk3, cast_to_uint8=True) + cp.testing.assert_array_equal(out_u, out_f) + + @pytest.mark.parametrize( + "method", + [ + "equalize", + "autolevel", + "gradient", + "majority", + "maximum", + "mean", + "geometric_mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "sum", + "threshold", + "noise_filter", + "entropy", + ], + ) + def test_compare_ubyte_vs_float_3d(self, method): + # Create signed int8 volume that and convert it to uint8 + np.random.seed(0) + volume_uint = np.random.randint( + 0, high=256, size=(10, 20, 30), dtype=np.uint8 + ) + volume_uint = cp.asarray(volume_uint) + volume_float = img_as_float(volume_uint) + + ball3 = ball(3, decomposition=None) + func = getattr(rank, method) + out_u = func(volume_uint, ball3) + with expected_warnings(["Possible precision loss"]): + out_f = func(volume_float, ball3, cast_to_uint8=True) + cp.testing.assert_array_equal(out_u, out_f) + + @pytest.mark.parametrize( + "method", + [ + "autolevel", + "equalize", + "gradient", + "maximum", + "mean", + "geometric_mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "threshold", + ], + ) + def test_compare_8bit_unsigned_vs_signed(self, method): + # filters applied on 8-bit image or 16-bit image (having only real 8-bit + # of dynamic) should be identical + + # Create signed int8 image that and convert it to uint8 + image = img_as_ubyte(cp.asarray(data.camera()))[::2, ::2] + image[image > 127] = 0 + image_s = image.astype(np.int8) + image_u = img_as_ubyte(image_s) + cp.testing.assert_array_equal(image_u, img_as_ubyte(image_s)) + func = getattr(rank, method) + disk3 = disk(3, decomposition=None) + out_u = func(image_u, disk3) + # with expected_warnings(["Possible precision loss"]): + out_s = func(image_s, disk3, cast_to_uint8=True) + cp.testing.assert_array_equal(out_u, out_s) + + @pytest.mark.parametrize( + "method", + [ + "equalize", + "autolevel", + "gradient", + "majority", + "maximum", + "mean", + "geometric_mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "sum", + "threshold", + "noise_filter", + "entropy", + ], + ) + def test_compare_8bit_unsigned_vs_signed_3d(self, method): + # filters applied on 8-bit volume or 16-bit volume (having only real 8-bit + # of dynamic) should be identical + + # Create signed int8 volume that and convert it to uint8 + np.random.seed(0) + volume_s = np.random.randint( + 0, high=127, size=(10, 20, 30), dtype=np.int8 + ) + volume_s = cp.asarray(volume_s) + volume_u = img_as_ubyte(volume_s) + cp.testing.assert_array_equal(volume_u, img_as_ubyte(volume_s)) + + ball3 = ball(3, decomposition=None) + func = getattr(rank, method) + out_u = func(volume_u, ball3) + # with expected_warnings(["Possible precision loss"]): + out_s = func(volume_s, ball3, cast_to_uint8=True) + cp.testing.assert_array_equal(out_u, out_s) + + @pytest.mark.parametrize( + "method", + [ + "autolevel", + "equalize", + "gradient", + "maximum", + "mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "threshold", + ], + ) + def test_compare_8bit_vs_16bit(self, method): + # filters applied on 8-bit image or 16-bit image (having only real 8-bit + # of dynamic) should be identical + image8 = util.img_as_ubyte(cp.asarray(data.camera())[::2, ::2]) + image16 = image8.astype(cp.uint16) + cp.testing.assert_array_equal(image8, image16) + + func = getattr(rank, method) + + disk3 = disk(3, decomposition=None) + f8 = func(image8, disk3) + f16 = func(image16, disk3, cast_to_uint8=True) + cp.testing.assert_array_equal(f8, f16) + + @pytest.mark.parametrize( + "method", + [ + "equalize", + "autolevel", + "gradient", + "majority", + "maximum", + "mean", + "geometric_mean", + "subtract_mean", + "median", + "minimum", + "modal", + "enhance_contrast", + "pop", + "sum", + "threshold", + "noise_filter", + "entropy", + ], + ) + def test_compare_8bit_vs_16bit_3d(self, method): + np.random.seed(0) + volume8 = np.random.randint( + 128, high=256, size=(10, 10, 10), dtype=np.uint8 + ) + volume8 = cp.asarray(volume8) + volume16 = volume8.astype(cp.uint16) + + func = getattr(rank, method) + + ball3 = ball(3, decomposition=None) + f8 = func(volume8, ball3) + f16 = func(volume16, ball3, cast_to_uint8=True) + cp.testing.assert_array_equal(f8, f16) + + @pytest.mark.parametrize("dtype", [cp.uint8, cp.uint16]) + def test_trivial_footprint8(self, dtype): + # check that min, max and mean returns identity if footprint + # contains only central pixel + + image = cp.zeros((5, 5), dtype=dtype) + out = cp.zeros_like(image) + mask = cp.ones_like(image, dtype=cp.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = cp.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=cp.uint8) + rank.mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.geometric_mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.minimum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.maximum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + + @pytest.mark.parametrize("dtype", [cp.uint8, cp.uint16]) + def test_smallest_footprint8(self, dtype): + # check that min, max and mean returns identity if footprint + # contains only central pixel + + image = cp.zeros((5, 5), dtype=dtype) + out = cp.zeros_like(image) + mask = cp.ones_like(image, dtype=cp.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = cp.array([[1]], dtype=cp.uint8) + rank.mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.minimum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.maximum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + + def test_empty_footprint(self): + image = cp.zeros((5, 5), dtype=np.uint16) + out = cp.zeros_like(image) + mask = cp.ones_like(image, dtype=np.uint8) + res = cp.zeros_like(image) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + elem = cp.array([[0, 0, 0], [0, 0, 0]], dtype=np.uint8) + + rank.mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(res, out) + rank.geometric_mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(res, out) + rank.minimum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(res, out) + rank.maximum( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(res, out) + + def test_entropy(self): + # verify that entropy is coherent with bitdepth of the input data + + footprint = cp.ones((16, 16), dtype=cp.uint8) + # 1 bit per pixel + data = cp.tile(cp.asarray([0, 1]), (100, 100)).astype(cp.uint8) + assert cp.max(rank.entropy(data, footprint)) == 1 + + # 2 bit per pixel + data = cp.tile(cp.asarray([[0, 1], [2, 3]]), (10, 10)).astype(cp.uint8) + assert cp.max(rank.entropy(data, footprint)) == 2 + + # 3 bit per pixel + data = cp.tile( + cp.asarray([[0, 1, 2, 3], [4, 5, 6, 7]]), (10, 10) + ).astype(cp.uint8) + assert cp.max(rank.entropy(data, footprint)) == 3 + + # 4 bit per pixel + data = cp.tile(cp.reshape(cp.arange(16), (4, 4)), (10, 10)).astype( + cp.uint8 + ) + assert cp.max(rank.entropy(data, footprint)) == 4 + + # 6 bit per pixel + data = cp.tile(cp.reshape(cp.arange(64), (8, 8)), (10, 10)).astype( + cp.uint8 + ) + assert cp.max(rank.entropy(data, footprint)) == 6 + + # 8-bit per pixel + data = cp.tile(cp.reshape(cp.arange(256), (16, 16)), (10, 10)).astype( + cp.uint8 + ) + assert cp.max(rank.entropy(data, footprint)) == 8 + + # 12 bit per pixel + footprint = cp.ones((64, 64), dtype=cp.uint8) + data = cp.zeros((65, 65), dtype=cp.uint16) + data[:64, :64] = cp.reshape(cp.arange(4096), (64, 64)) + assert cp.max(rank.entropy(data, footprint, cast_to_uint8=False)) == 12 + + # make sure output is floating point + # with expected_warnings(['Bad rank filter performance']): + out = rank.entropy(data, cp.ones((16, 16), dtype=cp.uint8)) + assert out.dtype.kind == "f" + + def test_footprint_dtypes(self): + image = cp.zeros((5, 5), dtype=cp.uint8) + out = cp.zeros_like(image) + mask = cp.ones_like(image, dtype=cp.uint8) + image[2, 2] = 255 + image[2, 3] = 128 + image[1, 2] = 16 + + for dtype in ( + bool, + cp.uint8, + cp.uint16, + cp.int32, + cp.int64, + cp.float32, + cp.float64, + ): + elem = cp.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], dtype=dtype) + rank.mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.geometric_mean( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + rank.mean_percentile( + image=image, + footprint=elem, + out=out, + mask=mask, + shift_x=0, + shift_y=0, + ) + cp.testing.assert_array_equal(image, out) + + def test_16bit(self): + image = cp.zeros((21, 21), dtype=np.uint16) + footprint = cp.ones((3, 3), dtype=np.uint8) + + for bitdepth in range(17): + value = 2**bitdepth - 1 + image[10, 10] = value + expected = [] + # if bitdepth >= 11: + # expected = ['Bad rank filter performance'] + with expected_warnings(expected): + assert ( + rank.minimum(image, footprint, cast_to_uint8=False)[10, 10] + == 0 + ) + assert ( + rank.maximum(image, footprint, cast_to_uint8=False)[10, 10] + == value + ) + mean_val = rank.mean(image, footprint, cast_to_uint8=False)[ + 10, 10 + ] + assert mean_val == int(value / footprint.size) + + def test_bilateral(self): + image = cp.zeros((21, 21), dtype=cp.uint16) + footprint = cp.ones((3, 3), dtype=cp.uint8) + + image[10, 10] = 1000 + image[10, 11] = 1010 + image[10, 9] = 900 + + kwargs = dict(s0=1, s1=1, cast_to_uint8=False) + assert ( + rank.mean_bilateral(image, footprint, **kwargs)[10, 10].get() + == 1000 + ) + assert rank.pop_bilateral(image, footprint, **kwargs)[10, 10].get() == 1 + kwargs = dict(s0=11, s1=11, cast_to_uint8=False) + assert ( + rank.mean_bilateral(image, footprint, **kwargs)[10, 10].get() + == 1005 + ) + assert rank.pop_bilateral(image, footprint, **kwargs)[10, 10].get() == 2 + + def test_percentile_min(self): + # check that percentile p0 = 0 is identical to local min + img = cp.asarray(data.camera()) + img16 = img.astype(cp.uint16) + footprint = disk(15, decomposition=None) + # check for 8bit + img_p0 = rank.percentile(img, footprint=footprint, p0=0) + img_min = rank.minimum(img, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_min) + # check for 16bit + img_p0 = rank.percentile(img16, footprint=footprint, p0=0) + img_min = rank.minimum(img16, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_min) + + def test_percentile_max(self): + # check that percentile p0 = 1 is identical to local max + img = cp.asarray(data.camera()) + img16 = img.astype(cp.uint16) + footprint = disk(15, decomposition=None) + # check for 8bit + img_p0 = rank.percentile(img, footprint=footprint, p0=1.0) + img_max = rank.maximum(img, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_max) + # check for 16bit + img_p0 = rank.percentile(img16, footprint=footprint, p0=1.0) + img_max = rank.maximum(img16, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_max) + + def test_percentile_median(self): + # check that percentile p0 = 0.5 is identical to local median + img = cp.asarray(data.camera()) + img16 = img.astype(cp.uint16) + footprint = disk(15, decomposition=None) + # check for 8bit + img_p0 = rank.percentile(img, footprint=footprint, p0=0.5) + img_max = rank.median(img, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_max) + # check for 16bit + img_p0 = rank.percentile(img16, footprint=footprint, p0=0.5) + img_max = rank.median(img16, footprint=footprint) + cp.testing.assert_array_equal(img_p0, img_max) + + def test_sum(self): + # check the number of valid pixels in the neighborhood + + image8 = cp.array( + [ + [0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0], + ], + dtype=cp.uint8, + ) + image16 = 400 * cp.array( + [ + [0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0], + ], + dtype=cp.uint16, + ) + elem = cp.ones((3, 3), dtype=cp.uint8) + out8 = cp.empty_like(image8) + out16 = cp.empty_like(image16) + mask = cp.ones(image8.shape, dtype=cp.uint8) + + r = cp.array( + [ + [1, 2, 3, 2, 1], + [2, 4, 6, 4, 2], + [3, 6, 9, 6, 3], + [2, 4, 6, 4, 2], + [1, 2, 3, 2, 1], + ], + dtype=cp.uint8, + ) + rank.sum(image=image8, footprint=elem, out=out8, mask=mask) + cp.testing.assert_array_equal(r, out8) + rank.sum_percentile( + image=image8, footprint=elem, out=out8, mask=mask, p0=0.0, p1=1.0 + ) + cp.testing.assert_array_equal(r, out8) + rank.sum_bilateral( + image=image8, footprint=elem, out=out8, mask=mask, s0=255, s1=255 + ) + cp.testing.assert_array_equal(r, out8) + + r = 400 * cp.array( + [ + [1, 2, 3, 2, 1], + [2, 4, 6, 4, 2], + [3, 6, 9, 6, 3], + [2, 4, 6, 4, 2], + [1, 2, 3, 2, 1], + ], + dtype=cp.uint16, + ) + rank.sum( + image=image16, + footprint=elem, + out=out16, + mask=mask, + cast_to_uint8=False, + ) + cp.testing.assert_array_equal(r, out16) + rank.sum_percentile( + image=image16, + footprint=elem, + out=out16, + mask=mask, + p0=0.0, + p1=1.0, + cast_to_uint8=False, + ) + cp.testing.assert_array_equal(r, out16) + rank.sum_bilateral( + image=image16, + footprint=elem, + out=out16, + mask=mask, + s0=1000, + s1=1000, + cast_to_uint8=False, + ) + cp.testing.assert_array_equal(r, out16) + + def test_median_default_value(self): + a = cp.zeros((3, 3), dtype=cp.uint8) + a[1] = 1 + full_footprint = cp.ones((3, 3), dtype=cp.uint8) + cp.testing.assert_array_equal( + rank.median(a), rank.median(a, full_footprint) + ) + assert rank.median(a)[1, 1].get() == 0 + assert rank.median(a, disk(1, decomposition=None))[1, 1].get() == 1 + + def test_output_same_dtype(self): + image = (cp.random.rand(100, 100) * 256).astype(cp.uint8) + out = cp.empty_like(image) + mask = cp.ones(image.shape, dtype=cp.uint8) + elem = cp.ones((3, 3), dtype=cp.uint8) + rank.maximum(image=image, footprint=elem, out=out, mask=mask) + cp.testing.assert_array_equal(image.dtype, out.dtype) + + def test_input_boolean_dtype(self): + image = (cp.random.rand(100, 100) * 256).astype(bool) + elem = cp.ones((3, 3), dtype=bool) + with pytest.raises(ValueError): + rank.maximum(image=image, footprint=elem) diff --git a/python/cucim/src/cucim/skimage/filters/tests/test_median.py b/python/cucim/src/cucim/skimage/filters/tests/test_median.py index b6c0abebc..66f35ddbe 100644 --- a/python/cucim/src/cucim/skimage/filters/tests/test_median.py +++ b/python/cucim/src/cucim/skimage/filters/tests/test_median.py @@ -1,5 +1,5 @@ # 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 import math @@ -12,7 +12,7 @@ from cucim.skimage import morphology from cucim.skimage._shared.testing import expected_warnings -from cucim.skimage.filters import median +from cucim.skimage.filters import median, rank @pytest.fixture @@ -34,14 +34,12 @@ def camera(): return cp.array(data.camera()) -# TODO: mode='rank' disabled until it has been implemented @pytest.mark.parametrize( "mode, cval, behavior, warning_type", [ ("nearest", 0.0, "ndimage", []), - # ('constant', 0.0, 'rank', (UserWarning,)), - # ('nearest', 0.0, 'rank', []), - ("nearest", 0.0, "ndimage", []), + ("constant", 0.0, "rank", (UserWarning,)), + ("nearest", 0.0, "rank", []), ], ) def test_median_warning(image, mode, cval, behavior, warning_type): @@ -52,7 +50,6 @@ def test_median_warning(image, mode, cval, behavior, warning_type): median(image, mode=mode, behavior=behavior) -# TODO: update if rank.median implemented @pytest.mark.parametrize( "behavior, func", [("ndimage", ndimage.median_filter)], @@ -95,6 +92,34 @@ def test_median_behavior( ) +@pytest.mark.parametrize("footprint_tuple", (False, True)) +@pytest.mark.parametrize("out", [None, "array"]) +def test_median_rank_behavior_matches_rank_median(camera, footprint_tuple, out): + footprint_shape = (3, 5) + if footprint_tuple: + footprint = footprint_shape + rank_footprint = cp.ones(footprint_shape, dtype=bool) + else: + footprint = cp.ones(footprint_shape, dtype=bool) + rank_footprint = footprint + cam2 = camera[:64, :75] + out_arg = cp.zeros_like(cam2) if out == "array" else None + expected_out = cp.zeros_like(cam2) if out == "array" else None + + result = median( + cam2, + footprint, + behavior="rank", + out=out_arg, + ) + expected = rank.median(cam2, footprint=rank_footprint, out=expected_out) + + if out_arg is not None: + assert result is out_arg + assert expected is expected_out + assert_allclose(result, expected) + + @pytest.mark.parametrize( "mode", ["reflect", "mirror", "nearest", "constant", "wrap"] ) @@ -262,18 +287,12 @@ def test_median_preserve_dtype(image, dtype): assert median_image.dtype == dtype -# TODO: update if rank.median implemented -# def test_median_error_ndim(): -# img = cp.random.randint(0, 10, size=(5, 5, 5), dtype=cp.uint8) -# with pytest.raises(ValueError): -# median(img, behavior='rank') - - -# TODO: update if rank.median implemented @pytest.mark.parametrize( "img, behavior", - # (cp.random.randint(0, 10, size=(3, 3), dtype=cp.uint8), 'rank'), [ + (cp.random.randint(0, 10, size=(3, 3), dtype=cp.uint8), "rank"), + # note: upstream scikit-image is 2D-only in rank mode, but cuCIM is nD + (cp.random.randint(0, 10, size=(3, 3, 3), dtype=cp.uint8), "rank"), (cp.random.randint(0, 10, size=(3, 3), dtype=cp.uint8), "ndimage"), (cp.random.randint(0, 10, size=(3, 3, 3), dtype=cp.uint8), "ndimage"), ],