Repository: rapidsai/cucim (language: python)

Issue:
[FEA] refactor separable convolution for improved performance
Benchmarking vs. [VPI](https://docs.nvidia.com/vpi/) indicates that a several-fold improvement may be possible. I will try to work on an implementation making better use of shared memory. Ideally this would be an n-d version of the row and column filters implemented here: https://github.com/opencv/opencv_contrib/tree/4.x/modules/cudafilters/src/cuda




Discussion / hints:

Fix patch:
diff --git a/benchmarks/skimage/bench_convolve.py b/benchmarks/skimage/bench_convolve.py
new file mode 100644
index 000000000..60007b8a3
--- /dev/null
+++ b/benchmarks/skimage/bench_convolve.py
@@ -0,0 +1,139 @@
+"""
+Benchmark locally modified ndimage functions vs. their CuPy counterparts
+"""
+import cupy as cp
+import cupyx.scipy.ndimage as ndi
+import pytest
+from cupyx.profiler import benchmark
+
+from cucim.skimage._vendored.ndimage import (
+    convolve1d, correlate1d, gaussian_filter, gaussian_filter1d,
+    gaussian_gradient_magnitude, gaussian_laplace, laplace, prewitt, sobel,
+    uniform_filter, uniform_filter1d,
+)
+
+d = cp.cuda.Device()
+
+
+def _get_image(shape, dtype, seed=123):
+
+    rng = cp.random.default_rng(seed)
+    dtype = cp.dtype(dtype)
+    if dtype.kind == 'b':
+        image = rng.integers(0, 1, shape, dtype=cp.uint8).astype(bool)
+    elif dtype.kind in 'iu':
+        image = rng.integers(0, 128, shape, dtype=dtype)
+    elif dtype.kind in 'c':
+        real_dtype = cp.asarray([], dtype=dtype).real.dtype
+        image = rng.standard_normal(shape, dtype=real_dtype)
+        image = image + 1j * rng.standard_normal(shape, dtype=real_dtype)
+    else:
+        if dtype == cp.float16:
+            image = rng.standard_normal(shape).astype(dtype)
+        else:
+            image = rng.standard_normal(shape, dtype=dtype)
+    return image
+
+
+def _compare_implementations(
+    shape, kernel_size, axis, dtype, mode, cval=0.0, origin=0,
+    output_dtype=None, kernel_dtype=None, output_preallocated=False,
+    function=convolve1d, max_duration=1
+):
+    dtype = cp.dtype(dtype)
+    if kernel_dtype is None:
+        kernel_dtype = dtype
+    image = _get_image(shape, dtype)
+    kernel = _get_image((kernel_size,), kernel_dtype)
+    kwargs = dict(axis=axis, mode=mode, cval=cval, origin=origin)
+    if output_dtype is not None:
+        output_dtype = cp.dtype(output_dtype)
+    function_ref = getattr(ndi, function.__name__)
+    if output_preallocated:
+        if output_dtype is None:
+            output_dtype = image.dtype
+        output1 = cp.empty(image.shape, dtype=output_dtype)
+        output2 = cp.empty(image.shape, dtype=output_dtype)
+        kwargs.update(dict(output=output1))
+        perf1 = benchmark(function_ref, (image, kernel), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+        kwargs.update(dict(output=output2, algorithm='shared_memory'))
+        perf2 = benchmark(function, (image, kernel), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+        return perf1, perf2
+    kwargs.update(dict(output=output_dtype))
+    perf1 = benchmark(function_ref, (image, kernel), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+    kwargs.update(dict(output=output_dtype, algorithm='shared_memory'))
+    perf2 = benchmark(function, (image, kernel), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+    return perf1, perf2
+
+
+def _compare_implementations_other(
+    shape, dtype, mode, cval=0.0,
+    output_dtype=None, kernel_dtype=None, output_preallocated=False,
+    function=convolve1d, func_kwargs={}, max_duration=1,
+):
+    dtype = cp.dtype(dtype)
+    image = _get_image(shape, dtype)
+    kwargs = dict(mode=mode, cval=cval)
+    if func_kwargs:
+        kwargs.update(func_kwargs)
+    if output_dtype is not None:
+        output_dtype = cp.dtype(output_dtype)
+    function_ref = getattr(ndi, function.__name__)
+    if output_preallocated:
+        if output_dtype is None:
+            output_dtype = image.dtype
+        output1 = cp.empty(image.shape, dtype=output_dtype)
+        output2 = cp.empty(image.shape, dtype=output_dtype)
+        kwargs.update(dict(output=output1))
+        perf1 = benchmark(function_ref, (image,), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+        kwargs.update(dict(output=output1, algorithm='shared_memory'))
+        perf2 = benchmark(function, (image,), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+        return perf1, perf2
+    kwargs.update(dict(output=output_dtype))
+    perf1 = benchmark(function_ref, (image,), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+    kwargs.update(dict(output=output_dtype, algorithm='shared_memory'))
+    perf2 = benchmark(function, (image,), kwargs=kwargs, n_warmup=10, n_repeat=10000, max_duration=max_duration)
+    return perf1, perf2
+
+
+print("\n\n")
+print("function | shape | dtype | mode | kernel size | preallocated | axis | dur (ms), CuPy | dur (ms), cuCIM | acceleration ")
+print("---------|-------|-------|------|-------------|--------------|------|----------------|-----------------|--------------")
+for function in [convolve1d]:
+    for shape in [(512, 512), (3840, 2160), (64, 64, 64), (256, 256, 256)]:
+        for dtype in [cp.float32, cp.uint8]:
+            for mode in ['nearest']:
+                for kernel_size in [3, 7, 11, 41]:
+                    for output_preallocated in [False]:  # , True]:
+                        for axis in range(len(shape)):
+                            output_dtype = dtype
+                            perf1, perf2 = _compare_implementations(shape=shape, kernel_size=kernel_size, mode=mode, axis=axis, dtype=dtype, output_dtype=output_dtype, output_preallocated=output_preallocated, function=function)
+                            t_elem = perf1.gpu_times * 1000.
+                            t_shared = perf2.gpu_times * 1000.
+                            print(f"{function.__name__} | {shape} | {cp.dtype(dtype).name} | {mode} | {kernel_size=} | prealloc={output_preallocated} | {axis=} | {t_elem.mean():0.3f} +/- {t_elem.std():0.3f}  | {t_shared.mean():0.3f} +/- {t_shared.std():0.3f} | {t_elem.mean() / t_shared.mean():0.3f}")
+
+
+print("function | kwargs | shape | dtype | mode | preallocated | dur (ms), CuPy | dur (ms), cuCIM | acceleration ")
+print("---------|--------|-------|-------|------|--------------|----------------|-----------------|--------------")
+for function, func_kwargs in [
+    # (gaussian_filter1d, dict(sigma=1.0, axis=0)),
+    # (gaussian_filter1d, dict(sigma=1.0, axis=-1)),
+    # (gaussian_filter1d, dict(sigma=4.0, axis=0)),
+    # (gaussian_filter1d, dict(sigma=4.0, axis=-1)),
+    (gaussian_filter, dict(sigma=1.0)),
+    (gaussian_filter, dict(sigma=4.0)),
+    (uniform_filter, dict(size=11)),
+    (prewitt, dict(axis=0)),
+    (sobel, dict(axis=0)),
+    (prewitt, dict(axis=-1)),
+    (sobel, dict(axis=-1)),
+]:
+    for shape in [(512, 512), (3840, 2160), (64, 64, 64), (256, 256, 256)]:
+        for (dtype, output_dtype) in [(cp.float32, cp.float32), (cp.uint8, cp.float32)]:
+            for mode in ['nearest']:
+                for output_preallocated in [False, True]:
+                    perf1, perf2 = _compare_implementations_other(shape=shape, mode=mode, dtype=dtype, output_dtype=output_dtype, output_preallocated=output_preallocated, function=function, func_kwargs=func_kwargs)
+                    t_elem = perf1.gpu_times * 1000.
+                    t_shared = perf2.gpu_times * 1000.
+                    print(f"{function.__name__} | {func_kwargs} | {shape} | {cp.dtype(dtype).name} | {mode} | {output_preallocated} | {t_elem.mean():0.3f} +/- {t_elem.std():0.3f}  | {t_shared.mean():0.3f} +/- {t_shared.std():0.3f} | {t_elem.mean() / t_shared.mean():0.3f}")
+
diff --git a/python/cucim/setup.cfg b/python/cucim/setup.cfg
index 517f33fff..7743540c2 100644
--- a/python/cucim/setup.cfg
+++ b/python/cucim/setup.cfg
@@ -65,6 +65,6 @@ line_length = 80
 known_first_party = cucim
 default_section = THIRDPARTY
 forced_separate = test_cucim
-skip = .tox,.eggs,ci/templates,build,dist,versioneer.py
+skip = .tox,.eggs,ci/templates,build,dist,versioneer.py,ndimage.py
 sections = FUTURE,STDLIB,THIRDPARTY,FIRSTPARTY,LOCALFOLDER
 multi_line_output = GRID
diff --git a/python/cucim/src/cucim/skimage/_shared/filters.py b/python/cucim/src/cucim/skimage/_shared/filters.py
index 2f68221b4..56cee8000 100644
--- a/python/cucim/src/cucim/skimage/_shared/filters.py
+++ b/python/cucim/src/cucim/skimage/_shared/filters.py
@@ -7,7 +7,8 @@
 from collections.abc import Iterable
 
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
+
 
 from .._shared import utils
 from .._shared.utils import _supported_float_type, convert_to_float
diff --git a/python/cucim/src/cucim/skimage/_vendored/_internal.py b/python/cucim/src/cucim/skimage/_vendored/_internal.py
index b6d163a36..196486d94 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_internal.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_internal.py
@@ -1,3 +1,7 @@
+import math
+from functools import reduce
+from operator import mul
+
 import cupy
 import numpy
 
@@ -61,3 +65,11 @@ def _normalize_axis_indices(axes, ndim):  # NOQA
             res.append(axis)
 
         return tuple(sorted(res))
+
+
+if hasattr(math, 'prod'):
+    prod = math.prod
+else:
+
+    def prod(iterable, *, start=1):
+        return reduce(mul, iterable, start)
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
index 1ce475371..e4dc1311d 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
@@ -1,8 +1,224 @@
 """A vendored subset of cupyx.scipy.ndimage._filters"""
+import warnings
 
 import cupy
+import numpy
 
+from cucim.skimage._vendored import _internal as internal
 from cucim.skimage._vendored import _ndimage_filters_core as _filters_core
+from cucim.skimage._vendored import _ndimage_util as _util
+
+
+def correlate(input, weights, output=None, mode='reflect', cval=0.0, origin=0):
+    """Multi-dimensional correlate.
+
+    The array is correlated with the given kernel.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        weights (cupy.ndarray): Array of weights, same number of dimensions as
+            input
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of correlate.
+
+    .. seealso:: :func:`scipy.ndimage.correlate`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    return _correlate_or_convolve(input, weights, output, mode, cval, origin)
+
+
+def convolve(input, weights, output=None, mode='reflect', cval=0.0, origin=0):
+    """Multi-dimensional convolution.
+
+    The array is convolved with the given kernel.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        weights (cupy.ndarray): Array of weights, same number of dimensions as
+            input
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of convolution.
+
+    .. seealso:: :func:`scipy.ndimage.convolve`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    return _correlate_or_convolve(input, weights, output, mode, cval, origin,
+                                  True)
+
+
+def correlate1d(input, weights, axis=-1, output=None, mode="reflect", cval=0.0,
+                origin=0, *, algorithm=None):
+    """One-dimensional correlate.
+
+    The array is correlated with the given kernel.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        weights (cupy.ndarray): One-dimensional array of weights
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int): The origin parameter controls the placement of the
+            filter, relative to the center of the current element of the
+            input. Default is ``0``.
+
+    Returns:
+        cupy.ndarray: The result of the 1D correlation.
+
+    .. seealso:: :func:`scipy.ndimage.correlate1d`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    return _correlate_or_convolve1d(
+        input, weights, axis, output, mode, cval, origin, False, algorithm
+    )
+
+
+def convolve1d(input, weights, axis=-1, output=None, mode="reflect", cval=0.0,
+               origin=0, *, algorithm=None):
+    """One-dimensional convolution.
+
+    The array is convolved with the given kernel.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        weights (cupy.ndarray): One-dimensional array of weights
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int): The origin parameter controls the placement of the
+            filter, relative to the center of the current element of the
+            input. Default is ``0``.
+    Returns:
+        cupy.ndarray: The result of the 1D convolution.
+
+    .. seealso:: :func:`scipy.ndimage.convolve1d`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    return _correlate_or_convolve1d(
+        input, weights, axis, output, mode, cval, origin, True, algorithm
+    )
+
+
+def _correlate_or_convolve(input, weights, output, mode, cval, origin,
+                           convolution=False):
+    origins, int_type = _filters_core._check_nd_args(input, weights,
+                                                     mode, origin)
+    if weights.size == 0:
+        return cupy.zeros_like(input)
+
+    _util._check_cval(mode, cval, _util._is_integer_output(output, input))
+
+    if convolution:
+        weights = weights[tuple([slice(None, None, -1)] * weights.ndim)]
+        origins = list(origins)
+        for i, wsize in enumerate(weights.shape):
+            origins[i] = -origins[i]
+            if wsize % 2 == 0:
+                origins[i] -= 1
+        origins = tuple(origins)
+    elif weights.dtype.kind == "c":
+        # numpy.correlate conjugates weights rather than input.
+        weights = weights.conj()
+    weights_dtype = _util._get_weights_dtype(input, weights, use_cucim_casting=True)  # noqa
+    offsets = _filters_core._origins_to_offsets(origins, weights.shape)
+    kernel = _get_correlate_kernel(mode, weights.shape, int_type,
+                                   offsets, cval)
+    output = _filters_core._call_kernel(kernel, input, weights, output,
+                                        weights_dtype=weights_dtype)
+    return output
+
+
+def _correlate_or_convolve1d(input, weights, axis, output, mode, cval, origin,
+                             convolution=False, algorithm=None):
+    # Calls fast shared-memory convolution when possible, otherwise falls back
+    # to the vendored elementwise _correlate_or_convolve
+    if algorithm is None:
+        if input.ndim == 2 and weights.size <= 256:
+            algorithm = 'shared_memory'
+        else:
+            algorithm = 'elementwise'
+    elif algorithm not in ['shared_memory', 'elementwise']:
+        raise ValueError(
+            "algorithm must be 'shared_memory', 'elementwise' or None"
+        )
+    if mode == 'wrap':
+        mode = 'grid-wrap'
+    if algorithm == 'shared_memory':
+        from cucim.skimage.filters._separable_filtering import (
+            ResourceLimitError, _shmem_convolve1d)
+        if input.ndim not in [2, 3]:
+            raise NotImplementedError(
+                f"shared_memory not implemented for ndim={input.ndim}"
+            )
+        try:
+            out = _shmem_convolve1d(input, weights, axis=axis, output=output,
+                                    mode=mode, cval=cval, origin=origin,
+                                    convolution=convolution)
+            return out
+        except ResourceLimitError:
+            # fallback to elementwise if inadequate shared memory available
+            warnings.warn(
+                "Inadequate resources for algorithm='shared_memory: "
+                "falling back to the elementwise implementation"
+            )
+            algorithm = 'elementwise'
+    if algorithm == 'elementwise':
+        weights, origins = _filters_core._convert_1d_args(
+            input.ndim, weights, origin, axis
+        )
+        return _correlate_or_convolve(
+            input, weights, output, mode, cval, origins, convolution
+        )
 
 
 @cupy.memoize(for_each_device=True)
@@ -13,3 +229,949 @@ def _get_correlate_kernel(mode, w_shape, int_type, offsets, cval):
         'sum += cast<W>({value}) * wval;',
         'y = cast<Y>(sum);',
         mode, w_shape, int_type, offsets, cval, ctype='W')
+
+
+def _run_1d_correlates(input, params, get_weights, output, mode, cval,
+                       origin=0, **filter_kwargs):
+    """
+    Enhanced version of _run_1d_filters that uses correlate1d as the filter
+    function. The params are a list of values to pass to the get_weights
+    callable given. If duplicate param values are found, the weights are
+    reused from the first invocation of get_weights. The get_weights callable
+    must return a 1D array of weights to give to correlate1d.
+    """
+    wghts = {}
+    for param in params:
+        if param not in wghts:
+            wghts[param] = get_weights(param)
+    wghts = [wghts[param] for param in params]
+    return _filters_core._run_1d_filters(
+        [None if w is None else correlate1d for w in wghts],
+        input, wghts, output, mode, cval, origin, **filter_kwargs)
+
+
+def uniform_filter1d(input, size, axis=-1, output=None, mode="reflect",
+                     cval=0.0, origin=0, *, algorithm=None):
+    """One-dimensional uniform filter along the given axis.
+
+    The lines of the array along the given axis are filtered with a uniform
+    filter of the given size.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int): Length of the uniform filter.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int): The origin parameter controls the placement of the
+            filter, relative to the center of the current element of the
+            input. Default is ``0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.uniform_filter1d`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+    weights = cupy.full(size, 1 / size, dtype=weights_dtype)
+    return correlate1d(input, weights, axis, output, mode, cval,
+                       origin, algorithm=algorithm)
+
+
+def uniform_filter(input, size=3, output=None, mode="reflect", cval=0.0,
+                   origin=0, *, algorithm=None):
+    """Multi-dimensional uniform filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int or sequence of int): Lengths of the uniform filter for each
+            dimension. A single value applies to all axes.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of ``0`` is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.uniform_filter`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    sizes = _util._fix_sequence_arg(size, input.ndim, 'size', int)
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+
+    def get(size):
+        return None if size <= 1 else cupy.full(size, 1 / size, dtype=weights_dtype)  # noqa
+
+    return _run_1d_correlates(
+        input, sizes, get, output, mode, cval, origin, algorithm=algorithm
+    )
+
+
+def gaussian_filter1d(input, sigma, axis=-1, order=0, output=None,
+                      mode="reflect", cval=0.0, truncate=4.0, *,
+                      algorithm=None):
+    """One-dimensional Gaussian filter along the given axis.
+
+    The lines of the array along the given axis are filtered with a Gaussian
+    filter of the given standard deviation.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        sigma (scalar): Standard deviation for Gaussian kernel.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        order (int): An order of ``0``, the default, corresponds to convolution
+            with a Gaussian kernel. A positive order corresponds to convolution
+            with that derivative of a Gaussian.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        truncate (float): Truncate the filter at this many standard deviations.
+            Default is ``4.0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.gaussian_filter1d`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    radius = int(float(truncate) * float(sigma) + 0.5)
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+    weights = _gaussian_kernel1d(sigma, int(order), radius, weights_dtype)
+    return correlate1d(
+        input, weights, axis, output, mode, cval, algorithm=algorithm
+    )
+
+
+def gaussian_filter(input, sigma, order=0, output=None, mode="reflect",
+                    cval=0.0, truncate=4.0, *, algorithm=None):
+    """Multi-dimensional Gaussian filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        sigma (scalar or sequence of scalar): Standard deviations for each axis
+            of Gaussian kernel. A single value applies to all axes.
+        order (int or sequence of scalar): An order of ``0``, the default,
+            corresponds to convolution with a Gaussian kernel. A positive order
+            corresponds to convolution with that derivative of a Gaussian. A
+            single value applies to all axes.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        truncate (float): Truncate the filter at this many standard deviations.
+            Default is ``4.0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.gaussian_filter`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    sigmas = _util._fix_sequence_arg(sigma, input.ndim, 'sigma', float)
+    orders = _util._fix_sequence_arg(order, input.ndim, 'order', int)
+    truncate = float(truncate)
+    weights_dtype = cupy.promote_types(input, cupy.float32)
+
+    def get(param, dtype=weights_dtype):
+        sigma, order = param
+        radius = int(truncate * float(sigma) + 0.5)
+        if radius <= 0:
+            return None
+        return _gaussian_kernel1d(sigma, order, radius, dtype)
+
+    return _run_1d_correlates(input, list(zip(sigmas, orders)), get, output,
+                              mode, cval, 0, algorithm=algorithm)
+
+
+def _gaussian_kernel1d(sigma, order, radius, dtype=cupy.float64):
+    """
+    Computes a 1-D Gaussian correlation kernel.
+    """
+    if order < 0:
+        raise ValueError('order must be non-negative')
+    sigma2 = sigma * sigma
+    x = numpy.arange(-radius, radius + 1)
+    phi_x = numpy.exp(-0.5 / sigma2 * x ** 2)
+    phi_x /= phi_x.sum()
+
+    if order == 0:
+        return cupy.asarray(phi_x)
+
+    # f(x) = q(x) * phi(x) = q(x) * exp(p(x))
+    # f'(x) = (q'(x) + q(x) * p'(x)) * phi(x)
+    # p'(x) = -1 / sigma ** 2
+    # Implement q'(x) + q(x) * p'(x) as a matrix operator and apply to the
+    # coefficients of q(x)
+    exponent_range = numpy.arange(order + 1)
+    q = numpy.zeros(order + 1)
+    q[0] = 1
+    D = numpy.diag(exponent_range[1:], 1)  # D @ q(x) = q'(x)
+    P = numpy.diag(numpy.ones(order) / -sigma2, -1)  # P @ q(x) = q(x) * p'(x)
+    Q_deriv = D + P
+    for _ in range(order):
+        q = Q_deriv.dot(q)
+    q = (x[:, None] ** exponent_range).dot(q)
+    return cupy.asarray((q * phi_x)[::-1], order='C', dtype=dtype)
+
+
+def prewitt(input, axis=-1, output=None, mode="reflect", cval=0.0, *,
+            algorithm=None):
+    """Compute a Prewitt filter along the given axis.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.prewitt`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+    smooth = cupy.ones(3, dtype=weights_dtype)
+    return _prewitt_or_sobel(
+        input, axis, output, mode, cval, smooth, algorithm
+    )
+
+
+def sobel(input, axis=-1, output=None, mode="reflect", cval=0.0, *,
+          algorithm=None):
+    """Compute a Sobel filter along the given axis.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.sobel`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+    smooth = cupy.array([1, 2, 1], dtype=weights_dtype)
+    return _prewitt_or_sobel(
+        input, axis, output, mode, cval, smooth, algorithm
+    )
+
+
+def _prewitt_or_sobel(input, axis, output, mode, cval, weights, algorithm):
+    axis = internal._normalize_axis_index(axis, input.ndim)
+
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+
+    def get(is_diff, dtype=weights_dtype):
+        return cupy.array([-1, 0, 1], dtype=dtype) if is_diff else weights  # noqa
+
+    return _run_1d_correlates(input, [a == axis for a in range(input.ndim)],
+                              get, output, mode, cval, algorithm=algorithm)
+
+
+def generic_laplace(input, derivative2, output=None, mode="reflect",
+                    cval=0.0, extra_arguments=(), extra_keywords=None):
+    """Multi-dimensional Laplace filter using a provided second derivative
+    function.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        derivative2 (callable): Function or other callable with the following
+            signature that is called once per axis::
+
+                derivative2(input, axis, output, mode, cval,
+                            *extra_arguments, **extra_keywords)
+
+            where ``input`` and ``output`` are ``cupy.ndarray``, ``axis`` is an
+            ``int`` from ``0`` to the number of dimensions, and ``mode``,
+            ``cval``, ``extra_arguments``, ``extra_keywords`` are the values
+            given to this function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        extra_arguments (sequence, optional):
+            Sequence of extra positional arguments to pass to ``derivative2``.
+        extra_keywords (dict, optional):
+            dict of extra keyword arguments to pass ``derivative2``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.generic_laplace`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    if extra_keywords is None:
+        extra_keywords = {}
+    ndim = input.ndim
+    modes = _util._fix_sequence_arg(mode, ndim, 'mode',
+                                    _util._check_mode)
+    output = _util._get_output(output, input)
+    if ndim == 0:
+        output[:] = input
+        return output
+    derivative2(input, 0, output, modes[0], cval,
+                *extra_arguments, **extra_keywords)
+    if ndim > 1:
+        tmp = _util._get_output(output.dtype, input)
+        for i in range(1, ndim):
+            derivative2(input, i, tmp, modes[i], cval,
+                        *extra_arguments, **extra_keywords)
+            output += tmp
+    return output
+
+
+def laplace(input, output=None, mode="reflect", cval=0.0, *, algorithm=None):
+    """Multi-dimensional Laplace filter based on approximate second
+    derivatives.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.laplace`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    weights_dtype = cupy.promote_types(input.dtype, cupy.float32)
+    weights = cupy.array([1, -2, 1], dtype=weights_dtype)
+
+    def derivative2(input, axis, output, mode, cval):
+        return correlate1d(
+            input, weights, axis, output, mode, cval, algorithm=algorithm
+        )
+
+    return generic_laplace(input, derivative2, output, mode, cval)
+
+
+def gaussian_laplace(input, sigma, output=None, mode="reflect",
+                     cval=0.0, *, algorithm=None, **kwargs):
+    """Multi-dimensional Laplace filter using Gaussian second derivatives.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        sigma (scalar or sequence of scalar): Standard deviations for each axis
+            of Gaussian kernel. A single value applies to all axes.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        kwargs (dict, optional):
+            dict of extra keyword arguments to pass ``gaussian_filter()``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.gaussian_laplace`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    def derivative2(input, axis, output, mode, cval):
+        order = [0] * input.ndim
+        order[axis] = 2
+        return gaussian_filter(input, sigma, order, output, mode, cval,
+                               algorithm=algorithm, **kwargs)
+    return generic_laplace(input, derivative2, output, mode, cval)
+
+
+def generic_gradient_magnitude(input, derivative, output=None,
+                               mode="reflect", cval=0.0,
+                               extra_arguments=(), extra_keywords=None):
+    """Multi-dimensional gradient magnitude filter using a provided derivative
+    function.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        derivative (callable): Function or other callable with the following
+            signature that is called once per axis::
+
+                derivative(input, axis, output, mode, cval,
+                           *extra_arguments, **extra_keywords)
+
+            where ``input`` and ``output`` are ``cupy.ndarray``, ``axis`` is an
+            ``int`` from ``0`` to the number of dimensions, and ``mode``,
+            ``cval``, ``extra_arguments``, ``extra_keywords`` are the values
+            given to this function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        extra_arguments (sequence, optional):
+            Sequence of extra positional arguments to pass to ``derivative2``.
+        extra_keywords (dict, optional):
+            dict of extra keyword arguments to pass ``derivative2``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.generic_gradient_magnitude`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    if extra_keywords is None:
+        extra_keywords = {}
+    ndim = input.ndim
+    modes = _util._fix_sequence_arg(mode, ndim, 'mode',
+                                    _util._check_mode)
+    output = _util._get_output(output, input)
+    if ndim == 0:
+        output[:] = input
+        return output
+    derivative(input, 0, output, modes[0], cval,
+               *extra_arguments, **extra_keywords)
+    output *= output
+    if ndim > 1:
+        tmp = _util._get_output(output.dtype, input)
+        for i in range(1, ndim):
+            derivative(input, i, tmp, modes[i], cval,
+                       *extra_arguments, **extra_keywords)
+            tmp *= tmp
+            output += tmp
+    return cupy.sqrt(output, output, casting='unsafe')
+
+
+def gaussian_gradient_magnitude(input, sigma, output=None, mode="reflect",
+                                cval=0.0, *, algorithm=None, **kwargs):
+    """Multi-dimensional gradient magnitude using Gaussian derivatives.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        sigma (scalar or sequence of scalar): Standard deviations for each axis
+            of Gaussian kernel. A single value applies to all axes.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        kwargs (dict, optional):
+            dict of extra keyword arguments to pass ``gaussian_filter()``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.gaussian_gradient_magnitude`
+
+    .. note::
+        When the output data type is integral (or when no output is provided
+        and input is integral) the results may not perfectly match the results
+        from SciPy due to floating-point rounding of intermediate results.
+    """
+    def derivative(input, axis, output, mode, cval):
+        order = [0] * input.ndim
+        order[axis] = 1
+        return gaussian_filter(input, sigma, order, output, mode, cval,
+                               algorithm=algorithm, **kwargs)
+    return generic_gradient_magnitude(input, derivative, output, mode, cval)
+
+
+def minimum_filter(input, size=None, footprint=None, output=None,
+                   mode="reflect", cval=0.0, origin=0):
+    """Multi-dimensional minimum filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int or sequence of int): One of ``size`` or ``footprint`` must be
+            provided. If ``footprint`` is given, ``size`` is ignored. Otherwise
+            ``footprint = cupy.ones(size)`` with ``size`` automatically made to
+            match the number of dimensions in ``input``.
+        footprint (cupy.ndarray): a boolean array which specifies which of the
+            elements within this shape will get passed to the filter function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.minimum_filter`
+    """
+    return _min_or_max_filter(input, size, footprint, None, output, mode,
+                              cval, origin, 'min')
+
+
+def maximum_filter(input, size=None, footprint=None, output=None,
+                   mode="reflect", cval=0.0, origin=0):
+    """Multi-dimensional maximum filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int or sequence of int): One of ``size`` or ``footprint`` must be
+            provided. If ``footprint`` is given, ``size`` is ignored. Otherwise
+            ``footprint = cupy.ones(size)`` with ``size`` automatically made to
+            match the number of dimensions in ``input``.
+        footprint (cupy.ndarray): a boolean array which specifies which of the
+            elements within this shape will get passed to the filter function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.maximum_filter`
+    """
+    return _min_or_max_filter(input, size, footprint, None, output, mode,
+                              cval, origin, 'max')
+
+
+def _min_or_max_filter(input, size, ftprnt, structure, output, mode, cval,
+                       origin, func):
+    # structure is used by morphology.grey_erosion() and grey_dilation()
+    # and not by the regular min/max filters
+
+    sizes, ftprnt, structure = _filters_core._check_size_footprint_structure(
+        input.ndim, size, ftprnt, structure)
+    if cval is cupy.nan:
+        raise NotImplementedError("NaN cval is unsupported")
+
+    if sizes is not None:
+        # Seperable filter, run as a series of 1D filters
+        fltr = minimum_filter1d if func == 'min' else maximum_filter1d
+        return _filters_core._run_1d_filters(
+            [fltr if size > 1 else None for size in sizes],
+            input, sizes, output, mode, cval, origin)
+
+    origins, int_type = _filters_core._check_nd_args(input, ftprnt,
+                                                     mode, origin, 'footprint')
+    if structure is not None and structure.ndim != input.ndim:
+        raise RuntimeError('structure array has incorrect shape')
+
+    if ftprnt.size == 0:
+        return cupy.zeros_like(input)
+    offsets = _filters_core._origins_to_offsets(origins, ftprnt.shape)
+    kernel = _get_min_or_max_kernel(mode, ftprnt.shape, func,
+                                    offsets, float(cval), int_type,
+                                    has_structure=structure is not None,
+                                    has_central_value=bool(ftprnt[offsets]))
+    return _filters_core._call_kernel(kernel, input, ftprnt, output,
+                                      structure, weights_dtype=bool)
+
+
+def minimum_filter1d(input, size, axis=-1, output=None, mode="reflect",
+                     cval=0.0, origin=0):
+    """Compute the minimum filter along a single axis.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int): Length of the minimum filter.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int): The origin parameter controls the placement of the
+            filter, relative to the center of the current element of the
+            input. Default is ``0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.minimum_filter1d`
+    """
+    return _min_or_max_1d(input, size, axis, output, mode, cval, origin, 'min')
+
+
+def maximum_filter1d(input, size, axis=-1, output=None, mode="reflect",
+                     cval=0.0, origin=0):
+    """Compute the maximum filter along a single axis.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int): Length of the maximum filter.
+        axis (int): The axis of input along which to calculate. Default is -1.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int): The origin parameter controls the placement of the
+            filter, relative to the center of the current element of the
+            input. Default is ``0``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.maximum_filter1d`
+    """
+    return _min_or_max_1d(input, size, axis, output, mode, cval, origin, 'max')
+
+
+def _min_or_max_1d(input, size, axis=-1, output=None, mode="reflect", cval=0.0,
+                   origin=0, func='min'):
+    ftprnt = cupy.ones(size, dtype=bool)
+    ftprnt, origin = _filters_core._convert_1d_args(input.ndim, ftprnt,
+                                                    origin, axis)
+    origins, int_type = _filters_core._check_nd_args(input, ftprnt,
+                                                     mode, origin, 'footprint')
+    offsets = _filters_core._origins_to_offsets(origins, ftprnt.shape)
+    kernel = _get_min_or_max_kernel(mode, ftprnt.shape, func, offsets,
+                                    float(cval), int_type, has_weights=False)
+    return _filters_core._call_kernel(kernel, input, None, output,
+                                      weights_dtype=bool)
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_min_or_max_kernel(mode, w_shape, func, offsets, cval, int_type,
+                           has_weights=True, has_structure=False,
+                           has_central_value=True):
+    # When there are no 'weights' (the footprint, for the 1D variants) then
+    # we need to make sure intermediate results are stored as doubles for
+    # consistent results with scipy.
+    ctype = 'X' if has_weights else 'double'
+    value = '{value}'
+    if not has_weights:
+        value = 'cast<double>({})'.format(value)
+
+    # Having a non-flat structure biases the values
+    if has_structure:
+        value += ('-' if func == 'min' else '+') + 'cast<X>(sval)'
+
+    if has_central_value:
+        pre = '{} value = x[i];'
+        found = 'value = {func}({value}, value);'
+    else:
+        # If the central pixel is not included in the footprint we cannot
+        # assume `x[i]` is not below the min or above the max and thus cannot
+        # seed with that value. Instead we keep track of having set `value`.
+        pre = '{} value; bool set = false;'
+        found = 'value = set ? {func}({value}, value) : {value}; set=true;'
+
+    return _filters_core._generate_nd_kernel(
+        func, pre.format(ctype),
+        found.format(func=func, value=value), 'y = cast<Y>(value);',
+        mode, w_shape, int_type, offsets, cval, ctype=ctype,
+        has_weights=has_weights, has_structure=has_structure)
+
+
+def rank_filter(input, rank, size=None, footprint=None, output=None,
+                mode="reflect", cval=0.0, origin=0):
+    """Multi-dimensional rank filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        rank (int): The rank of the element to get. Can be negative to count
+            from the largest value, e.g. ``-1`` indicates the largest value.
+        size (int or sequence of int): One of ``size`` or ``footprint`` must be
+            provided. If ``footprint`` is given, ``size`` is ignored. Otherwise
+            ``footprint = cupy.ones(size)`` with ``size`` automatically made to
+            match the number of dimensions in ``input``.
+        footprint (cupy.ndarray): a boolean array which specifies which of the
+            elements within this shape will get passed to the filter function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.rank_filter`
+    """
+    rank = int(rank)
+    return _rank_filter(input, lambda fs: rank + fs if rank < 0 else rank,
+                        size, footprint, output, mode, cval, origin)
+
+
+def median_filter(input, size=None, footprint=None, output=None,
+                  mode="reflect", cval=0.0, origin=0):
+    """Multi-dimensional median filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (int or sequence of int): One of ``size`` or ``footprint`` must be
+            provided. If ``footprint`` is given, ``size`` is ignored. Otherwise
+            ``footprint = cupy.ones(size)`` with ``size`` automatically made to
+            match the number of dimensions in ``input``.
+        footprint (cupy.ndarray): a boolean array which specifies which of the
+            elements within this shape will get passed to the filter function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.median_filter`
+    """
+    return _rank_filter(input, lambda fs: fs // 2,
+                        size, footprint, output, mode, cval, origin)
+
+
+def percentile_filter(input, percentile, size=None, footprint=None,
+                      output=None, mode="reflect", cval=0.0, origin=0):
+    """Multi-dimensional percentile filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        percentile (scalar): The percentile of the element to get (from ``0``
+            to ``100``). Can be negative, thus ``-20`` equals ``80``.
+        size (int or sequence of int): One of ``size`` or ``footprint`` must be
+            provided. If ``footprint`` is given, ``size`` is ignored. Otherwise
+            ``footprint = cupy.ones(size)`` with ``size`` automatically made to
+            match the number of dimensions in ``input``.
+        footprint (cupy.ndarray): a boolean array which specifies which of the
+            elements within this shape will get passed to the filter function.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output. Default is is same dtype as the input.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``'constant'``. Default is ``0.0``.
+        origin (int or sequence of int): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of the filtering.
+
+    .. seealso:: :func:`scipy.ndimage.percentile_filter`
+    """
+    percentile = float(percentile)
+    if percentile < 0.0:
+        percentile += 100.0
+    if percentile < 0.0 or percentile > 100.0:
+        raise RuntimeError('invalid percentile')
+    if percentile == 100.0:
+        def get_rank(fs):
+            return fs - 1
+    else:
+        def get_rank(fs):
+            return int(float(fs) * percentile / 100.0)
+    return _rank_filter(input, get_rank,
+                        size, footprint, output, mode, cval, origin)
+
+
+def _rank_filter(input, get_rank, size=None, footprint=None, output=None,
+                 mode="reflect", cval=0.0, origin=0):
+    _, footprint, _ = _filters_core._check_size_footprint_structure(
+        input.ndim, size, footprint, None, force_footprint=True)
+    if cval is cupy.nan:
+        raise NotImplementedError("NaN cval is unsupported")
+    origins, int_type = _filters_core._check_nd_args(input, footprint,
+                                                     mode, origin, 'footprint')
+    if footprint.size == 0:
+        return cupy.zeros_like(input)
+    filter_size = int(footprint.sum())
+    rank = get_rank(filter_size)
+    if rank < 0 or rank >= filter_size:
+        raise RuntimeError('rank not within filter footprint size')
+    if rank == 0:
+        return _min_or_max_filter(input, None, footprint, None, output, mode,
+                                  cval, origins, 'min')
+    if rank == filter_size - 1:
+        return _min_or_max_filter(input, None, footprint, None, output, mode,
+                                  cval, origins, 'max')
+    offsets = _filters_core._origins_to_offsets(origins, footprint.shape)
+    kernel = _get_rank_kernel(filter_size, rank, mode, footprint.shape,
+                              offsets, float(cval), int_type)
+    return _filters_core._call_kernel(kernel, input, footprint, output,
+                                      weights_dtype=bool)
+
+
+__SHELL_SORT = '''
+__device__ void sort(X *array, int size) {{
+    int gap = {gap};
+    while (gap > 1) {{
+        gap /= 3;
+        for (int i = gap; i < size; ++i) {{
+            X value = array[i];
+            int j = i - gap;
+            while (j >= 0 && value < array[j]) {{
+                array[j + gap] = array[j];
+                j -= gap;
+            }}
+            array[j + gap] = value;
+        }}
+    }}
+}}'''
+
+
+@cupy._util.memoize()
+def _get_shell_gap(filter_size):
+    gap = 1
+    while gap < filter_size:
+        gap = 3 * gap + 1
+    return gap
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_rank_kernel(filter_size, rank, mode, w_shape, offsets, cval,
+                     int_type):
+    s_rank = min(rank, filter_size - rank - 1)
+    # The threshold was set based on the measurements on a V100
+    # TODO(leofang, anaruse): Use Optuna to automatically tune the threshold,
+    # as it may vary depending on the GPU in use, compiler version, dtype,
+    # filter size, etc.
+    if s_rank <= 80:
+        # When s_rank is small and register usage is low, this partial
+        # selection sort approach is faster than general sorting approach
+        # using shell sort.
+        if s_rank == rank:
+            comp_op = '<'
+        else:
+            comp_op = '>'
+        array_size = s_rank + 2
+        found_post = '''
+            if (iv > {rank} + 1) {{{{
+                int target_iv = 0;
+                X target_val = values[0];
+                for (int jv = 1; jv <= {rank} + 1; jv++) {{{{
+                    if (target_val {comp_op} values[jv]) {{{{
+                        target_val = values[jv];
+                        target_iv = jv;
+                    }}}}
+                }}}}
+                if (target_iv <= {rank}) {{{{
+                    values[target_iv] = values[{rank} + 1];
+                }}}}
+                iv = {rank} + 1;
+            }}}}'''.format(rank=s_rank, comp_op=comp_op)
+        post = '''
+            X target_val = values[0];
+            for (int jv = 1; jv <= {rank}; jv++) {{
+                if (target_val {comp_op} values[jv]) {{
+                    target_val = values[jv];
+                }}
+            }}
+            y=cast<Y>(target_val);'''.format(rank=s_rank, comp_op=comp_op)
+        sorter = ''
+    else:
+        array_size = filter_size
+        found_post = ''
+        post = 'sort(values,{});\ny=cast<Y>(values[{}]);'.format(
+            filter_size, rank)
+        sorter = __SHELL_SORT.format(gap=_get_shell_gap(filter_size))
+
+    return _filters_core._generate_nd_kernel(
+        'rank_{}_{}'.format(filter_size, rank),
+        'int iv = 0;\nX values[{}];'.format(array_size),
+        'values[iv++] = {value};' + found_post, post,
+        mode, w_shape, int_type, offsets, cval, preamble=sorter)
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 60aec68cf..469dff9d9 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py
@@ -1,24 +1,172 @@
 """A vendored subset of cupyx.scipy.ndimage._filters_core"""
+import warnings
 
 import cupy
 import numpy
 
+from cucim.skimage._vendored import _internal as internal
 from cucim.skimage._vendored import _ndimage_util as _util
 
-includes = r'''
-// workaround for HIP: line begins with #include
+
+def _origins_to_offsets(origins, w_shape):
+    return tuple(x // 2 + o for x, o in zip(w_shape, origins))
+
+
+def _check_size_footprint_structure(ndim, size, footprint, structure,
+                                    stacklevel=3, force_footprint=False):
+    if structure is None and footprint is None:
+        if size is None:
+            raise RuntimeError("no footprint or filter size provided")
+        sizes = _util._fix_sequence_arg(size, ndim, 'size', int)
+        if force_footprint:
+            return None, cupy.ones(sizes, bool), None
+        return sizes, None, None
+    if size is not None:
+        warnings.warn("ignoring size because {} is set".format(
+            'structure' if footprint is None else 'footprint'),
+            UserWarning, stacklevel=stacklevel + 1)
+
+    if footprint is not None:
+        footprint = cupy.array(footprint, bool, True, 'C')
+        if not footprint.any():
+            raise ValueError("all-zero footprint is not supported")
+
+    if structure is None:
+        if not force_footprint and footprint.all():
+            if footprint.ndim != ndim:
+                raise RuntimeError("size must have length equal to input rank")
+            return footprint.shape, None, None
+        return None, footprint, None
+
+    structure = cupy.ascontiguousarray(structure)
+    if footprint is None:
+        footprint = cupy.ones(structure.shape, bool)
+    return None, footprint, structure
+
+
+def _convert_1d_args(ndim, weights, origin, axis):
+    if weights.ndim != 1 or weights.size < 1:
+        raise RuntimeError('incorrect filter size')
+    axis = internal._normalize_axis_index(axis, ndim)
+    w_shape = [1] * ndim
+    w_shape[axis] = weights.size
+    weights = weights.reshape(w_shape)
+    origins = [0] * ndim
+    origins[axis] = _util._check_origin(origin, weights.size)
+    return weights, tuple(origins)
+
+
+def _check_nd_args(input, weights, mode, origin, wghts_name='filter weights'):
+    _util._check_mode(mode)
+    # Weights must always be less than 2 GiB
+    if weights.nbytes >= (1 << 31):
+        raise RuntimeError('weights must be 2 GiB or less, use FFTs instead')
+    weight_dims = [x for x in weights.shape if x != 0]
+    if len(weight_dims) != input.ndim:
+        raise RuntimeError('{} array has incorrect shape'.format(wghts_name))
+    origins = _util._fix_sequence_arg(origin, len(weight_dims), 'origin', int)
+    for origin, width in zip(origins, weight_dims):
+        _util._check_origin(origin, width)
+    return tuple(origins), _util._get_inttype(input)
+
+
+def _run_1d_filters(filters, input, args, output, mode, cval, origin=0,
+                    **filter_kwargs):
+    """
+    Runs a series of 1D filters forming an nd filter. The filters must be a
+    list of callables that take input, arg, axis, output, mode, cval, origin.
+    The args is a list of values that are passed for the arg value to the
+    filter. Individual filters can be None causing that axis to be skipped.
+    """
+    output = _util._get_output(output, input)
+    modes = _util._fix_sequence_arg(mode, input.ndim, 'mode',
+                                    _util._check_mode)
+    # for filters, "wrap" is a synonym for "grid-wrap".
+    modes = ['grid-wrap' if m == 'wrap' else m for m in modes]
+    origins = _util._fix_sequence_arg(origin, input.ndim, 'origin', int)
+    n_filters = sum(filter is not None for filter in filters)
+    if n_filters == 0:
+        output[:] = input
+        return output
+    # We can't operate in-place efficiently, so use a 2-buffer system
+    temp = _util._get_output(output.dtype, input) if n_filters > 1 else None   # noqa
+    iterator = zip(filters, args, modes, origins)
+    for axis, (fltr, arg, mode, origin) in enumerate(iterator):
+        if fltr is None:
+            continue
+        else:
+            break
+    if n_filters % 2 == 0:
+        fltr(input, arg, axis, temp, mode, cval, origin, **filter_kwargs)
+        input = temp
+    else:
+        fltr(input, arg, axis, output, mode, cval, origin, **filter_kwargs)
+        if n_filters == 1:
+            return output
+        input, output = output, temp
+    for axis, (fltr, arg, mode, origin) in enumerate(iterator, start=axis + 1):
+        if fltr is None:
+            continue
+        fltr(input, arg, axis, output, mode, cval, origin, **filter_kwargs)
+        input, output = output, input
+    return input
+
+
+def _call_kernel(kernel, input, weights, output, structure=None,
+                 weights_dtype=numpy.float64, structure_dtype=numpy.float64):
+    """
+    Calls a constructed ElementwiseKernel. The kernel must take an input image,
+    an optional array of weights, an optional array for the structure, and an
+    output array.
+
+    weights and structure can be given as None (structure defaults to None) in
+    which case they are not passed to the kernel at all. If the output is given
+    as None then it will be allocated in this function.
+
+    This function deals with making sure that the weights and structure are
+    contiguous and float64 (or bool for weights that are footprints)*, that the
+    output is allocated and appriopately shaped. This also deals with the
+    situation that the input and output arrays overlap in memory.
+
+    * weights is always cast to float64 or bool in order to get an output
+    compatible with SciPy, though float32 might be sufficient when input dtype
+    is low precision. If weights_dtype is passed as weights.dtype then no
+    dtype conversion will occur. The input and output are never converted.
+    """
+    args = [input]
+    complex_output = input.dtype.kind == 'c'
+    if weights is not None:
+        weights = cupy.ascontiguousarray(weights, weights_dtype)
+        complex_output = complex_output or weights.dtype.kind == 'c'
+        args.append(weights)
+    if structure is not None:
+        structure = cupy.ascontiguousarray(structure, structure_dtype)
+        args.append(structure)
+    output = _util._get_output(output, input, None, complex_output)  # noqa
+    needs_temp = cupy.shares_memory(output, input, 'MAY_SHARE_BOUNDS')
+    if needs_temp:
+        output, temp = _util._get_output(output.dtype, input, None, complex_output), output  # noqa
+    args.append(output)
+    kernel(*args)
+    if needs_temp:
+        output[:] = temp
+        output = temp
+    return output
+
+
+_ndimage_includes = r'''
 #include <type_traits>  // let Jitify handle this
 #include <cupy/math_constants.h>
-'''
 
-
-_CAST_FUNCTION = """
-// Implements a casting function to make it compatible with scipy
-// Use like cast<to_type>(value)
 template<> struct std::is_floating_point<float16> : std::true_type {};
 template<> struct std::is_signed<float16> : std::true_type {};
 template<class T> struct std::is_signed<complex<T>> : std::is_signed<T> {};
+'''
 
+
+_ndimage_CAST_FUNCTION = """
+// Implements a casting function to make it compatible with scipy
+// Use like cast<to_type>(value)
 template <class B, class A>
 __device__ __forceinline__
 typename std::enable_if<(!std::is_floating_point<A>::value
@@ -148,7 +296,7 @@ def _generate_nd_kernel(name, pre, found, post, mode, w_shape, int_type,
         name += '_with_structure'
     if has_mask:
         name += '_with_mask'
-    preamble = includes + _CAST_FUNCTION + preamble
+    preamble = _ndimage_includes + _ndimage_CAST_FUNCTION + preamble
     options += ('--std=c++11', '-DCUPY_USE_JITIFY')
     return cupy.ElementwiseKernel(in_params, out_params, operation, name,
                                   reduce_dims=False, preamble=preamble,
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_interp_kernels.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interp_kernels.py
new file mode 100644
index 000000000..fa65eed36
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interp_kernels.py
@@ -0,0 +1,598 @@
+import cupy
+import numpy
+
+from cucim.skimage._vendored import \
+    _ndimage_spline_kernel_weights as _spline_kernel_weights
+from cucim.skimage._vendored import \
+    _ndimage_spline_prefilter_core as _spline_prefilter_core
+from cucim.skimage._vendored import _ndimage_util as _util
+
+math_constants_preamble = r'''
+// workaround for HIP: line begins with #include
+#include <cupy/math_constants.h>
+'''
+
+spline_weights_inline = _spline_kernel_weights.spline_weights_inline
+
+
+def _get_coord_map(ndim, nprepad=0):
+    """Extract target coordinate from coords array (for map_coordinates).
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        coords (ndarray): array of shape (ncoords, ndim) containing the target
+            coordinates.
+        c_j: variables to hold the target coordinates
+
+    computes::
+
+        c_j = coords[i + j * ncoords];
+
+    ncoords is determined by the size of the output array, y.
+    y will be indexed by the CIndexer, _ind.
+    Thus ncoords = _ind.size();
+
+    """
+    ops = []
+    ops.append('ptrdiff_t ncoords = _ind.size();')
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = coords[i + {j} * ncoords]{pre};''')
+    return ops
+
+
+def _get_coord_zoom_and_shift(ndim, nprepad=0):
+    """Compute target coordinate based on a shift followed by a zoom.
+
+    This version zooms from the center of the edge pixels.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        in_coord[ndim]: array containing the source coordinate
+        zoom[ndim]: array containing the zoom for each axis
+        shift[ndim]: array containing the zoom for each axis
+
+    computes::
+
+        c_j = zoom[j] * (in_coord[j] - shift[j])
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = zoom[{j}] * ((W)in_coord[{j}] - shift[{j}]){pre};''')
+    return ops
+
+
+def _get_coord_zoom_and_shift_grid(ndim, nprepad=0):
+    """Compute target coordinate based on a shift followed by a zoom.
+
+    This version zooms from the outer edges of the grid pixels.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        in_coord[ndim]: array containing the source coordinate
+        zoom[ndim]: array containing the zoom for each axis
+        shift[ndim]: array containing the zoom for each axis
+
+    computes::
+
+        c_j = zoom[j] * (in_coord[j] - shift[j] + 0.5) - 0.5
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = zoom[{j}] * ((W)in_coord[{j}] - shift[j] + 0.5) - 0.5{pre};''')
+    return ops
+
+
+def _get_coord_zoom(ndim, nprepad=0):
+    """Compute target coordinate based on a zoom.
+
+    This version zooms from the center of the edge pixels.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        in_coord[ndim]: array containing the source coordinate
+        zoom[ndim]: array containing the zoom for each axis
+
+    computes::
+
+        c_j = zoom[j] * in_coord[j]
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = zoom[{j}] * (W)in_coord[{j}]{pre};''')
+    return ops
+
+
+def _get_coord_zoom_grid(ndim, nprepad=0):
+    """Compute target coordinate based on a zoom (grid_mode=True version).
+
+    This version zooms from the outer edges of the grid pixels.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        in_coord[ndim]: array containing the source coordinate
+        zoom[ndim]: array containing the zoom for each axis
+
+    computes::
+
+        c_j = zoom[j] * (in_coord[j] + 0.5) - 0.5
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = zoom[{j}] * ((W)in_coord[{j}] + 0.5) - 0.5{pre};''')
+    return ops
+
+
+def _get_coord_shift(ndim, nprepad=0):
+    """Compute target coordinate based on a shift.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        in_coord[ndim]: array containing the source coordinate
+        shift[ndim]: array containing the zoom for each axis
+
+    computes::
+
+        c_j = in_coord[j] - shift[j]
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    for j in range(ndim):
+        ops.append(f'''
+    W c_{j} = (W)in_coord[{j}] - shift[{j}]{pre};''')
+    return ops
+
+
+def _get_coord_affine(ndim, nprepad=0):
+    """Compute target coordinate based on a homogeneous transformation matrix.
+
+    The homogeneous matrix has shape (ndim, ndim + 1). It corresponds to
+    affine matrix where the last row of the affine is assumed to be:
+    ``[0] * ndim + [1]``.
+
+    Notes
+    -----
+    Assumes the following variables have been initialized on the device::
+
+        mat(array): array containing the (ndim, ndim + 1) transform matrix.
+        in_coords(array): coordinates of the input
+
+    For example, in 2D:
+
+        c_0 = mat[0] * in_coords[0] + mat[1] * in_coords[1] + aff[2];
+        c_1 = mat[3] * in_coords[0] + mat[4] * in_coords[1] + aff[5];
+
+    """
+    ops = []
+    pre = f" + (W){nprepad}" if nprepad > 0 else ''
+    ncol = ndim + 1
+    for j in range(ndim):
+        ops.append(f'''
+            W c_{j} = (W)0.0;''')
+        for k in range(ndim):
+            ops.append(f'''
+            c_{j} += mat[{ncol * j + k}] * (W)in_coord[{k}];''')
+        ops.append(f'''
+            c_{j} += mat[{ncol * j + ndim}]{pre};''')
+    return ops
+
+
+def _unravel_loop_index(shape, uint_t='unsigned int'):
+    """
+    declare a multi-index array in_coord and unravel the 1D index, i into it.
+    This code assumes that the array is a C-ordered array.
+    """
+    ndim = len(shape)
+    code = [f'''
+        {uint_t} in_coord[{ndim}];
+        {uint_t} s, t, idx = i;''']
+    for j in range(ndim - 1, 0, -1):
+        code.append(f'''
+        s = {shape[j]};
+        t = idx / s;
+        in_coord[{j}] = idx - t * s;
+        idx = t;''')
+    code.append('''
+        in_coord[0] = idx;''')
+    return '\n'.join(code)
+
+
+def _generate_interp_custom(coord_func, ndim, large_int, yshape, mode, cval,
+                            order, name='', integer_output=False, nprepad=0,
+                            omit_in_coord=False):
+    """
+    Args:
+        coord_func (function): generates code to do the coordinate
+            transformation. See for example, `_get_coord_shift`.
+        ndim (int): The number of dimensions.
+        large_int (bool): If true use Py_ssize_t instead of int for indexing.
+        yshape (tuple): Shape of the output array.
+        mode (str): Signal extension mode to use at the array boundaries
+        cval (float): constant value used when `mode == 'constant'`.
+        name (str): base name for the interpolation kernel
+        integer_output (bool): boolean indicating whether the output has an
+            integer type.
+        nprepad (int): integer indicating the amount of prepadding at the
+            boundaries.
+
+    Returns:
+        operation (str): code body for the ElementwiseKernel
+        name (str): name for the ElementwiseKernel
+    """
+
+    ops = []
+    internal_dtype = 'double' if integer_output else 'Y'
+    ops.append(f'{internal_dtype} out = 0.0;')
+
+    if large_int:
+        uint_t = 'size_t'
+        int_t = 'ptrdiff_t'
+    else:
+        uint_t = 'unsigned int'
+        int_t = 'int'
+
+    # determine strides for x along each axis
+    for j in range(ndim):
+        ops.append(f'const {int_t} xsize_{j} = x.shape()[{j}];')
+    ops.append(f'const {uint_t} sx_{ndim - 1} = 1;')
+    for j in range(ndim - 1, 0, -1):
+        ops.append(f'const {uint_t} sx_{j - 1} = sx_{j} * xsize_{j};')
+
+    if not omit_in_coord:
+        # create in_coords array to store the unraveled indices
+        ops.append(_unravel_loop_index(yshape, uint_t))
+
+    # compute the transformed (target) coordinates, c_j
+    ops = ops + coord_func(ndim, nprepad)
+
+    if cval is numpy.nan:
+        cval = '(Y)CUDART_NAN'
+    elif cval == numpy.inf:
+        cval = '(Y)CUDART_INF'
+    elif cval == -numpy.inf:
+        cval = '(Y)(-CUDART_INF)'
+    else:
+        cval = f'({internal_dtype}){cval}'
+
+    if mode == 'constant':
+        # use cval if coordinate is outside the bounds of x
+        _cond = ' || '.join(
+            [f'(c_{j} < 0) || (c_{j} > xsize_{j} - 1)' for j in range(ndim)])
+        ops.append(f'''
+        if ({_cond})
+        {{
+            out = {cval};
+        }}
+        else
+        {{''')
+
+    if order == 0:
+        if mode == 'wrap':
+            ops.append('double dcoord;')  # mode 'wrap' requires this to work
+        for j in range(ndim):
+            # determine nearest neighbor
+            if mode == 'wrap':
+                ops.append(f'''
+                dcoord = c_{j};''')
+            else:
+                ops.append(f'''
+                {int_t} cf_{j} = ({int_t})floor((double)c_{j} + 0.5);''')
+
+            # handle boundary
+            if mode != 'constant':
+                if mode == 'wrap':
+                    ixvar = 'dcoord'
+                    float_ix = True
+                else:
+                    ixvar = f'cf_{j}'
+                    float_ix = False
+                ops.append(
+                    _util._generate_boundary_condition_ops(
+                        mode, ixvar, f'xsize_{j}', int_t, float_ix))
+                if mode == 'wrap':
+                    ops.append(f'''
+                {int_t} cf_{j} = ({int_t})floor(dcoord + 0.5);''')
+
+            # sum over ic_j will give the raveled coordinate in the input
+            ops.append(f'''
+            {int_t} ic_{j} = cf_{j} * sx_{j};''')
+        _coord_idx = ' + '.join([f'ic_{j}' for j in range(ndim)])
+        if mode == 'grid-constant':
+            _cond = ' || '.join([f'(ic_{j} < 0)' for j in range(ndim)])
+            ops.append(f'''
+            if ({_cond}) {{
+                out = {cval};
+            }} else {{
+                out = ({internal_dtype})x[{_coord_idx}];
+            }}''')
+        else:
+            ops.append(f'''
+            out = ({internal_dtype})x[{_coord_idx}];''')
+
+    elif order == 1:
+        for j in range(ndim):
+            # get coordinates for linear interpolation along axis j
+            ops.append(f'''
+            {int_t} cf_{j} = ({int_t})floor((double)c_{j});
+            {int_t} cc_{j} = cf_{j} + 1;
+            {int_t} n_{j} = (c_{j} == cf_{j}) ? 1 : 2;  // points needed
+            ''')
+
+            if mode == 'wrap':
+                ops.append(f'''
+                double dcoordf = c_{j};
+                double dcoordc = c_{j} + 1;''')
+            else:
+                # handle boundaries for extension modes.
+                ops.append(f'''
+                {int_t} cf_bounded_{j} = cf_{j};
+                {int_t} cc_bounded_{j} = cc_{j};''')
+
+            if mode != 'constant':
+                if mode == 'wrap':
+                    ixvar = 'dcoordf'
+                    float_ix = True
+                else:
+                    ixvar = f'cf_bounded_{j}'
+                    float_ix = False
+                ops.append(
+                    _util._generate_boundary_condition_ops(
+                        mode, ixvar, f'xsize_{j}', int_t, float_ix))
+
+                ixvar = 'dcoordc' if mode == 'wrap' else f'cc_bounded_{j}'
+                ops.append(
+                    _util._generate_boundary_condition_ops(
+                        mode, ixvar, f'xsize_{j}', int_t, float_ix))
+                if mode == 'wrap':
+                    ops.append(
+                        f'''
+                    {int_t} cf_bounded_{j} = ({int_t})floor(dcoordf);;
+                    {int_t} cc_bounded_{j} = ({int_t})floor(dcoordf + 1);;
+                    '''
+                    )
+
+            ops.append(f'''
+            for (int s_{j} = 0; s_{j} < n_{j}; s_{j}++)
+                {{
+                    W w_{j};
+                    {int_t} ic_{j};
+                    if (s_{j} == 0)
+                    {{
+                        w_{j} = (W)cc_{j} - c_{j};
+                        ic_{j} = cf_bounded_{j} * sx_{j};
+                    }} else
+                    {{
+                        w_{j} = c_{j} - (W)cf_{j};
+                        ic_{j} = cc_bounded_{j} * sx_{j};
+                    }}''')
+    elif order > 1:
+        if mode == 'grid-constant':
+            spline_mode = 'constant'
+        elif mode == 'nearest':
+            spline_mode = 'nearest'
+        else:
+            spline_mode = _spline_prefilter_core._get_spline_mode(mode)
+
+        # wx, wy are temporary variables used during spline weight computation
+        ops.append(f'''
+            W wx, wy;
+            {int_t} start;''')
+        for j in range(ndim):
+            # determine weights along the current axis
+            ops.append(f'''
+            W weights_{j}[{order + 1}];''')
+            ops.append(spline_weights_inline[order].format(j=j, order=order))
+
+            # get starting coordinate for spline interpolation along axis j
+            if mode in ['wrap']:
+                ops.append(f'double dcoord = c_{j};')
+                coord_var = 'dcoord'
+                ops.append(
+                    _util._generate_boundary_condition_ops(
+                        mode, coord_var, f'xsize_{j}', int_t, True))
+            else:
+                coord_var = f'(double)c_{j}'
+
+            if order & 1:
+                op_str = '''
+                start = ({int_t})floor({coord_var}) - {order_2};'''
+            else:
+                op_str = '''
+                start = ({int_t})floor({coord_var} + 0.5) - {order_2};'''
+            ops.append(
+                op_str.format(
+                    int_t=int_t, coord_var=coord_var, order_2=order // 2
+                ))
+
+            # set of coordinate values within spline footprint along axis j
+            ops.append(f'''{int_t} ci_{j}[{order + 1}];''')
+            for k in range(order + 1):
+                ixvar = f'ci_{j}[{k}]'
+                ops.append(f'''
+                {ixvar} = start + {k};''')
+                ops.append(
+                    _util._generate_boundary_condition_ops(
+                        spline_mode, ixvar, f'xsize_{j}', int_t))
+
+            # loop over the order + 1 values in the spline filter
+            ops.append(f'''
+            W w_{j};
+            {int_t} ic_{j};
+            for (int k_{j} = 0; k_{j} <= {order}; k_{j}++)
+                {{
+                    w_{j} = weights_{j}[k_{j}];
+                    ic_{j} = ci_{j}[k_{j}] * sx_{j};
+            ''')
+
+    if order > 0:
+
+        _weight = ' * '.join([f'w_{j}' for j in range(ndim)])
+        _coord_idx = ' + '.join([f'ic_{j}' for j in range(ndim)])
+        if mode == 'grid-constant' or (order > 1 and mode == 'constant'):
+            _cond = ' || '.join([f'(ic_{j} < 0)' for j in range(ndim)])
+            ops.append(f'''
+            if ({_cond}) {{
+                out += {cval} * ({internal_dtype})({_weight});
+            }} else {{
+                {internal_dtype} val = ({internal_dtype})x[{_coord_idx}];
+                out += val * ({internal_dtype})({_weight});
+            }}''')
+        else:
+            ops.append(f'''
+            {internal_dtype} val = ({internal_dtype})x[{_coord_idx}];
+            out += val * ({internal_dtype})({_weight});''')
+
+        ops.append('}' * ndim)
+
+    if mode == 'constant':
+        ops.append('}')
+
+    if integer_output:
+        ops.append('y = (Y)rint((double)out);')
+    else:
+        ops.append('y = (Y)out;')
+    operation = '\n'.join(ops)
+
+    mode_str = mode.replace('-', '_')  # avoid hyphen in kernel name
+    name = 'cupyx_scipy_ndimage_interpolate_{}_order{}_{}_{}d_y{}'.format(
+        name, order, mode_str, ndim, '_'.join([f'{j}' for j in yshape]),
+    )
+    if uint_t == 'size_t':
+        name += '_i64'
+    return operation, name
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_map_kernel(ndim, large_int, yshape, mode, cval=0.0, order=1,
+                    integer_output=False, nprepad=0):
+    in_params = 'raw X x, raw W coords'
+    out_params = 'Y y'
+    operation, name = _generate_interp_custom(
+        coord_func=_get_coord_map,
+        ndim=ndim,
+        large_int=large_int,
+        yshape=yshape,
+        mode=mode,
+        cval=cval,
+        order=order,
+        name='map',
+        integer_output=integer_output,
+        nprepad=nprepad,
+        omit_in_coord=True,  # input image coordinates are not needed
+    )
+    return cupy.ElementwiseKernel(in_params, out_params, operation, name,
+                                  preamble=math_constants_preamble)
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_shift_kernel(ndim, large_int, yshape, mode, cval=0.0, order=1,
+                      integer_output=False, nprepad=0):
+    in_params = 'raw X x, raw W shift'
+    out_params = 'Y y'
+    operation, name = _generate_interp_custom(
+        coord_func=_get_coord_shift,
+        ndim=ndim,
+        large_int=large_int,
+        yshape=yshape,
+        mode=mode,
+        cval=cval,
+        order=order,
+        name='shift',
+        integer_output=integer_output,
+        nprepad=nprepad,
+    )
+    return cupy.ElementwiseKernel(in_params, out_params, operation, name,
+                                  preamble=math_constants_preamble)
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_zoom_shift_kernel(ndim, large_int, yshape, mode, cval=0.0, order=1,
+                           integer_output=False, grid_mode=False, nprepad=0):
+    in_params = 'raw X x, raw W shift, raw W zoom'
+    out_params = 'Y y'
+    if grid_mode:
+        zoom_shift_func = _get_coord_zoom_and_shift_grid
+    else:
+        zoom_shift_func = _get_coord_zoom_and_shift
+    operation, name = _generate_interp_custom(
+        coord_func=zoom_shift_func,
+        ndim=ndim,
+        large_int=large_int,
+        yshape=yshape,
+        mode=mode,
+        cval=cval,
+        order=order,
+        name="zoom_shift_grid" if grid_mode else "zoom_shift",
+        integer_output=integer_output,
+        nprepad=nprepad,
+    )
+    return cupy.ElementwiseKernel(in_params, out_params, operation, name,
+                                  preamble=math_constants_preamble)
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_zoom_kernel(ndim, large_int, yshape, mode, cval=0.0, order=1,
+                     integer_output=False, grid_mode=False, nprepad=0):
+    in_params = 'raw X x, raw W zoom'
+    out_params = 'Y y'
+    operation, name = _generate_interp_custom(
+        coord_func=_get_coord_zoom_grid if grid_mode else _get_coord_zoom,
+        ndim=ndim,
+        large_int=large_int,
+        yshape=yshape,
+        mode=mode,
+        cval=cval,
+        order=order,
+        name="zoom_grid" if grid_mode else "zoom",
+        integer_output=integer_output,
+        nprepad=nprepad,
+    )
+    return cupy.ElementwiseKernel(in_params, out_params, operation, name,
+                                  preamble=math_constants_preamble)
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_affine_kernel(ndim, large_int, yshape, mode, cval=0.0, order=1,
+                       integer_output=False, nprepad=0):
+    in_params = 'raw X x, raw W mat'
+    out_params = 'Y y'
+    operation, name = _generate_interp_custom(
+        coord_func=_get_coord_affine,
+        ndim=ndim,
+        large_int=large_int,
+        yshape=yshape,
+        mode=mode,
+        cval=cval,
+        order=order,
+        name='affine',
+        integer_output=integer_output,
+        nprepad=nprepad,
+    )
+    return cupy.ElementwiseKernel(in_params, out_params, operation, name,
+                                  preamble=math_constants_preamble)
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py
new file mode 100644
index 000000000..18275213b
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py
@@ -0,0 +1,757 @@
+import cmath
+import math
+import warnings
+
+import cupy
+import numpy
+from cupy import _core
+from cupy.cuda import runtime
+
+from cucim.skimage._vendored import _ndimage_interp_kernels as _interp_kernels
+from cucim.skimage._vendored import \
+    _ndimage_spline_prefilter_core as _spline_prefilter_core
+from cucim.skimage._vendored import _ndimage_util as _util
+from cucim.skimage._vendored._internal import _normalize_axis_index, prod
+
+
+def _check_parameter(func_name, order, mode):
+    if order is None:
+        warnings.warn(f'Currently the default order of {func_name} is 1. In a '
+                      'future release this may change to 3 to match '
+                      'scipy.ndimage ')
+    elif order < 0 or 5 < order:
+        raise ValueError('spline order is not supported')
+
+    if mode not in ('constant', 'grid-constant', 'nearest', 'mirror',
+                    'reflect', 'grid-mirror', 'wrap', 'grid-wrap', 'opencv',
+                    '_opencv_edge'):
+        raise ValueError('boundary mode ({}) is not supported'.format(mode))
+
+
+def _get_spline_output(input, output):
+    """Create workspace array, temp, and the final dtype for the output.
+
+    Differs from SciPy by not always forcing the internal floating point dtype
+    to be double precision.
+    """
+    complex_data = input.dtype.kind == 'c'
+    if complex_data:
+        min_float_dtype = cupy.complex64
+    else:
+        min_float_dtype = cupy.float32
+    if isinstance(output, cupy.ndarray):
+        if complex_data and output.dtype.kind != 'c':
+            raise ValueError(
+                'output must have complex dtype for complex inputs'
+            )
+        float_dtype = cupy.promote_types(output.dtype, min_float_dtype)
+        output_dtype = output.dtype
+    else:
+        if output is None:
+            output = output_dtype = input.dtype
+        else:
+            output_dtype = cupy.dtype(output)
+        float_dtype = cupy.promote_types(output, min_float_dtype)
+
+    if (isinstance(output, cupy.ndarray)
+            and output.dtype == float_dtype == output_dtype
+            and output.flags.c_contiguous):
+        if output is not input:
+            _core.elementwise_copy(input, output)
+        temp = output
+    else:
+        temp = input.astype(float_dtype, copy=False)
+        temp = cupy.ascontiguousarray(temp)
+        if cupy.shares_memory(temp, input, 'MAY_SHARE_BOUNDS'):
+            temp = temp.copy()
+    return temp, float_dtype, output_dtype
+
+
+def spline_filter1d(input, order=3, axis=-1, output=cupy.float64,
+                    mode='mirror'):
+    """
+    Calculate a 1-D spline filter along the given axis.
+
+    The lines of the array along the given axis are filtered by a
+    spline filter. The order of the spline must be >= 2 and <= 5.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        axis (int): The axis along which the spline filter is applied. Default
+            is the last axis.
+        output (cupy.ndarray or dtype, optional): The array in which to place
+            the output, or the dtype of the returned array. Default is
+            ``numpy.float64``.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+
+    Returns:
+        cupy.ndarray: The result of prefiltering the input.
+
+    .. seealso:: :func:`scipy.spline_filter1d`
+    """
+    if order < 0 or order > 5:
+        raise RuntimeError('spline order not supported')
+    x = input
+    ndim = x.ndim
+    axis = _normalize_axis_index(axis, ndim)
+
+    # order 0, 1 don't require reshaping as no CUDA kernel will be called
+    # scalar or size 1 arrays also don't need to be filtered
+    run_kernel = not (order < 2 or x.ndim == 0 or x.shape[axis] == 1)
+    if not run_kernel:
+        output = _util._get_output(output, input)
+        _core.elementwise_copy(x, output)
+        return output
+
+    temp, data_dtype, output_dtype = _get_spline_output(x, output)
+    data_type = cupy._core._scalar.get_typename(temp.dtype)
+    pole_type = cupy._core._scalar.get_typename(temp.real.dtype)
+
+    index_type = _util._get_inttype(input)
+    index_dtype = cupy.int32 if index_type == 'int' else cupy.int64
+
+    n_samples = x.shape[axis]
+    n_signals = x.size // n_samples
+    info = cupy.array((n_signals, n_samples) + x.shape, dtype=index_dtype)
+
+    # empirical choice of block size that seemed to work well
+    block_size = max(2 ** math.ceil(numpy.log2(n_samples / 32)), 8)
+    kern = _spline_prefilter_core.get_raw_spline1d_kernel(
+        axis,
+        ndim,
+        mode,
+        order=order,
+        index_type=index_type,
+        data_type=data_type,
+        pole_type=pole_type,
+        block_size=block_size,
+    )
+
+    # Due to recursive nature, a given line of data must be processed by a
+    # single thread. n_signals lines will be processed in total.
+    block = (block_size,)
+    grid = ((n_signals + block[0] - 1) // block[0],)
+
+    # apply prefilter gain
+    poles = _spline_prefilter_core.get_poles(order=order)
+    temp *= _spline_prefilter_core.get_gain(poles)
+
+    # apply caual + anti-causal IIR spline filters
+    kern(grid, block, (temp, info))
+
+    if isinstance(output, cupy.ndarray) and temp is not output:
+        # copy kernel output into the user-provided output array
+        _core.elementwise_copy(temp, output)
+        return output
+    return temp.astype(output_dtype, copy=False)
+
+
+def spline_filter(input, order=3, output=cupy.float64, mode='mirror'):
+    """Multidimensional spline filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        output (cupy.ndarray or dtype, optional): The array in which to place
+            the output, or the dtype of the returned array. Default is
+            ``numpy.float64``.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+
+    Returns:
+        cupy.ndarray: The result of prefiltering the input.
+
+    .. seealso:: :func:`scipy.spline_filter1d`
+    """
+    if order < 2 or order > 5:
+        raise RuntimeError('spline order not supported')
+
+    x = input
+    temp, data_dtype, output_dtype = _get_spline_output(x, output)
+    if order not in [0, 1] and input.ndim > 0:
+        for axis in range(x.ndim):
+            spline_filter1d(x, order, axis, output=temp, mode=mode)
+            x = temp
+    if isinstance(output, cupy.ndarray):
+        _core.elementwise_copy(temp, output)
+    else:
+        output = temp
+    if output.dtype != output_dtype:
+        output = output.astype(output_dtype)
+    return output
+
+
+def _check_coordinates(coordinates, order, allow_float32=True):
+    if coordinates.dtype.kind == 'f':
+        if allow_float32:
+            coord_dtype = cupy.promote_types(coordinates.dtype, cupy.float32)
+        else:
+            coord_dtype = cupy.promote_types(coordinates.dtype, cupy.float64)
+        coordinates = coordinates.astype(coord_dtype, copy=False)
+    elif coordinates.dtype.kind in 'iu':
+        if order > 1:
+            # order > 1 (spline) kernels require floating-point coordinates
+            if allow_float32:
+                coord_dtype = cupy.promote_types(
+                    coordinates.dtype, cupy.float32
+                )
+            else:
+                coord_dtype = cupy.promote_types(
+                    coordinates.dtype, cupy.float64
+                )
+            coordinates = coordinates.astype(coord_dtype)
+    else:
+        raise ValueError('coordinates should have floating point dtype')
+    if not coordinates.flags.c_contiguous:
+        coordinates = cupy.ascontiguousarray(coordinates)
+    return coordinates
+
+
+def _prepad_for_spline_filter(input, mode, cval):
+    if mode in ['nearest', 'grid-constant']:
+        # these modes need padding to get accurate boundary values
+        npad = 12  # empirical factor chosen by SciPy
+        if mode == 'grid-constant':
+            kwargs = dict(mode='constant', constant_values=cval)
+        else:
+            kwargs = dict(mode='edge')
+        padded = cupy.pad(input, npad, **kwargs)
+    else:
+        npad = 0
+        padded = input
+    return padded, npad
+
+
+def _filter_input(image, prefilter, mode, cval, order):
+    """Perform spline prefiltering when needed.
+
+    Spline orders > 1 need a prefiltering stage to preserve resolution.
+
+    For boundary modes without analytical spline boundary conditions, some
+    prepadding of the input with cupy.pad is used to maintain accuracy.
+    ``npad`` is an integer corresponding to the amount of padding at each edge
+    of the array.
+    """
+    if not prefilter or order < 2:
+        return (cupy.ascontiguousarray(image), 0)
+    padded, npad = _prepad_for_spline_filter(image, mode, cval)
+    float_dtype = cupy.promote_types(image.dtype, cupy.float32)
+    filtered = spline_filter(padded, order, output=float_dtype, mode=mode)
+    return cupy.ascontiguousarray(filtered), npad
+
+
+def map_coordinates(input, coordinates, output=None, order=3,
+                    mode='constant', cval=0.0, prefilter=True):
+    """Map the input array to new coordinates by interpolation.
+
+    The array of coordinates is used to find, for each point in the output, the
+    corresponding coordinates in the input. The value of the input at those
+    coordinates is determined by spline interpolation of the requested order.
+
+    The shape of the output is derived from that of the coordinate array by
+    dropping the first axis. The values of the array along the first axis are
+    the coordinates in the input array at which the output value is found.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        coordinates (array_like): The coordinates at which ``input`` is
+            evaluated.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+        cval (scalar): Value used for points outside the boundaries of
+            the input if ``mode='constant'`` or ``mode='opencv'``. Default is
+            0.0
+        prefilter (bool): It is not used yet. It just exists for compatibility
+            with :mod:`scipy.ndimage`.
+
+    Returns:
+        cupy.ndarray:
+            The result of transforming the input. The shape of the output is
+            derived from that of ``coordinates`` by dropping the first axis.
+
+    .. seealso:: :func:`scipy.ndimage.map_coordinates`
+    """
+
+    _check_parameter('map_coordinates', order, mode)
+
+    if mode == 'opencv' or mode == '_opencv_edge':
+        input = cupy.pad(input, [(1, 1)] * input.ndim, 'constant',
+                         constant_values=cval)
+        coordinates = cupy.add(coordinates, 1)
+        mode = 'constant'
+
+    ret = _util._get_output(output, input, coordinates.shape[1:])
+    integer_output = ret.dtype.kind in 'iu'
+    _util._check_cval(mode, cval, integer_output)
+
+    if input.dtype.kind in 'iu':
+        input = input.astype(cupy.float32)
+    coordinates = _check_coordinates(coordinates, order)
+    filtered, nprepad = _filter_input(input, prefilter, mode, cval, order)
+    large_int = max(prod(input.shape), coordinates.shape[0]) > 1 << 31
+    kern = _interp_kernels._get_map_kernel(
+        input.ndim, large_int, yshape=coordinates.shape, mode=mode, cval=cval,
+        order=order, integer_output=integer_output, nprepad=nprepad)
+    kern(filtered, coordinates, ret)
+    return ret
+
+
+def affine_transform(input, matrix, offset=0.0, output_shape=None, output=None,
+                     order=3, mode='constant', cval=0.0, prefilter=True, *,
+                     texture_memory=False):
+    """Apply an affine transformation.
+
+    Given an output image pixel index vector ``o``, the pixel value is
+    determined from the input image at position
+    ``cupy.dot(matrix, o) + offset``.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        matrix (cupy.ndarray): The inverse coordinate transformation matrix,
+            mapping output coordinates to input coordinates. If ``ndim`` is the
+            number of dimensions of ``input``, the given matrix must have one
+            of the following shapes:
+
+                - ``(ndim, ndim)``: the linear transformation matrix for each
+                  output coordinate.
+                - ``(ndim,)``: assume that the 2D transformation matrix is
+                  diagonal, with the diagonal specified by the given value.
+                - ``(ndim + 1, ndim + 1)``: assume that the transformation is
+                  specified using homogeneous coordinates. In this case, any
+                  value passed to ``offset`` is ignored.
+                - ``(ndim, ndim + 1)``: as above, but the bottom row of a
+                  homogeneous transformation matrix is always
+                  ``[0, 0, ..., 1]``, and may be omitted.
+
+        offset (float or sequence): The offset into the array where the
+            transform is applied. If a float, ``offset`` is the same for each
+            axis. If a sequence, ``offset`` should contain one value for each
+            axis.
+        output_shape (tuple of ints): Shape tuple.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+        cval (scalar): Value used for points outside the boundaries of
+            the input if ``mode='constant'`` or ``mode='opencv'``. Default is
+            0.0
+        prefilter (bool): It is not used yet. It just exists for compatibility
+            with :mod:`scipy.ndimage`.
+        texture_memory (bool): If True, uses GPU texture memory. Supports only:
+
+            - 2D and 3D float32 arrays as input
+            - ``(ndim + 1, ndim + 1)`` homogeneous float32 transformation
+                matrix
+            - ``mode='constant'`` and ``mode='nearest'``
+            - ``order=0`` (nearest neighbor) and ``order=1`` (linear
+                interpolation)
+            - NVIDIA CUDA GPUs
+
+    Returns:
+        cupy.ndarray or None:
+            The transformed input. If ``output`` is given as a parameter,
+            ``None`` is returned.
+
+    .. seealso:: :func:`scipy.ndimage.affine_transform`
+    """
+
+    if texture_memory:
+        # _texture only available in CuPy 10.x so delay the import
+        # We do not use this texture-based implementation in cuCIM.
+        from cucim.skimage._vendored import _texture
+
+        if runtime.is_hip:
+            raise RuntimeError(
+                'HIP currently does not support texture acceleration')
+        tm_interp = 'linear' if order > 0 else 'nearest'
+        return _texture.affine_transformation(data=input,
+                                              transformation_matrix=matrix,
+                                              output_shape=output_shape,
+                                              output=output,
+                                              interpolation=tm_interp,
+                                              mode=mode,
+                                              border_value=cval)
+
+    _check_parameter('affine_transform', order, mode)
+
+    offset = _util._fix_sequence_arg(offset, input.ndim, 'offset', float)
+
+    if matrix.ndim not in [1, 2] or matrix.shape[0] < 1:
+        raise RuntimeError('no proper affine matrix provided')
+    if matrix.ndim == 2:
+        if matrix.shape[0] == matrix.shape[1] - 1:
+            offset = matrix[:, -1]
+            matrix = matrix[:, :-1]
+        elif matrix.shape[0] == input.ndim + 1:
+            offset = matrix[:-1, -1]
+            matrix = matrix[:-1, :-1]
+        if matrix.shape != (input.ndim, input.ndim):
+            raise RuntimeError('improper affine shape')
+
+    if mode == 'opencv':
+        m = cupy.zeros((input.ndim + 1, input.ndim + 1))
+        m[:-1, :-1] = matrix
+        m[:-1, -1] = offset
+        m[-1, -1] = 1
+        m = cupy.linalg.inv(m)
+        m[:2] = cupy.roll(m[:2], 1, axis=0)
+        m[:2, :2] = cupy.roll(m[:2, :2], 1, axis=1)
+        matrix = m[:-1, :-1]
+        offset = m[:-1, -1]
+
+    if output_shape is None:
+        output_shape = input.shape
+
+    if mode == 'opencv' or mode == '_opencv_edge':
+        if matrix.ndim == 1:
+            matrix = cupy.diag(matrix)
+        coordinates = cupy.indices(output_shape, dtype=cupy.float64)
+        coordinates = cupy.dot(matrix, coordinates.reshape((input.ndim, -1)))
+        coordinates += cupy.expand_dims(cupy.asarray(offset), -1)
+        ret = _util._get_output(output, input, shape=output_shape)
+        ret[:] = map_coordinates(input, coordinates, ret.dtype, order, mode,
+                                 cval, prefilter).reshape(output_shape)
+        return ret
+
+    matrix = matrix.astype(cupy.float64, copy=False)
+    ndim = input.ndim
+    output = _util._get_output(output, input, shape=output_shape)
+    if input.dtype.kind in 'iu':
+        input = input.astype(cupy.float32)
+    filtered, nprepad = _filter_input(input, prefilter, mode, cval, order)
+
+    integer_output = output.dtype.kind in 'iu'
+    _util._check_cval(mode, cval, integer_output)
+    large_int = max(prod(input.shape), prod(output_shape)) > 1 << 31
+    if matrix.ndim == 1:
+        offset = cupy.asarray(offset, dtype=cupy.float64)
+        offset = -offset / matrix
+        kern = _interp_kernels._get_zoom_shift_kernel(
+            ndim, large_int, output_shape, mode, cval=cval, order=order,
+            integer_output=integer_output, nprepad=nprepad)
+        kern(filtered, offset, matrix, output)
+    else:
+        kern = _interp_kernels._get_affine_kernel(
+            ndim, large_int, output_shape, mode, cval=cval, order=order,
+            integer_output=integer_output, nprepad=nprepad)
+        m = cupy.zeros((ndim, ndim + 1), dtype=cupy.float64)
+        m[:, :-1] = matrix
+        m[:, -1] = cupy.asarray(offset, dtype=cupy.float64)
+        kern(filtered, m, output)
+    return output
+
+
+def _minmax(coor, minc, maxc):
+    if coor[0] < minc[0]:
+        minc[0] = coor[0]
+    elif coor[0] > maxc[0]:
+        maxc[0] = coor[0]
+    if coor[1] < minc[1]:
+        minc[1] = coor[1]
+    elif coor[1] > maxc[1]:
+        maxc[1] = coor[1]
+    return minc, maxc
+
+
+def rotate(input, angle, axes=(1, 0), reshape=True, output=None, order=3,
+           mode='constant', cval=0.0, prefilter=True):
+    """Rotate an array.
+
+    The array is rotated in the plane defined by the two axes given by the
+    ``axes`` parameter using spline interpolation of the requested order.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        angle (float): The rotation angle in degrees.
+        axes (tuple of 2 ints): The two axes that define the plane of rotation.
+            Default is the first two axes.
+        reshape (bool): If ``reshape`` is True, the output shape is adapted so
+            that the input array is contained completely in the output. Default
+            is True.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+        cval (scalar): Value used for points outside the boundaries of
+            the input if ``mode='constant'`` or ``mode='opencv'``. Default is
+            0.0
+        prefilter (bool): It is not used yet. It just exists for compatibility
+            with :mod:`scipy.ndimage`.
+
+    Returns:
+        cupy.ndarray or None:
+            The rotated input.
+
+    .. seealso:: :func:`scipy.ndimage.rotate`
+    """
+
+    _check_parameter('rotate', order, mode)
+
+    if mode == 'opencv':
+        mode = '_opencv_edge'
+
+    input_arr = input
+    axes = list(axes)
+    if axes[0] < 0:
+        axes[0] += input_arr.ndim
+    if axes[1] < 0:
+        axes[1] += input_arr.ndim
+    if axes[0] > axes[1]:
+        axes = [axes[1], axes[0]]
+    if axes[0] < 0 or input_arr.ndim <= axes[1]:
+        raise ValueError('invalid rotation plane specified')
+
+    ndim = input_arr.ndim
+    rad = math.radians(angle)
+    sincos = cmath.rect(1, rad)
+    cos, sin = sincos.real, sincos.imag
+
+    # determine offsets and output shape as in scipy.ndimage.rotate
+    rot_matrix = numpy.array([[cos, sin],
+                              [-sin, cos]])
+
+    img_shape = numpy.asarray(input_arr.shape)
+    in_plane_shape = img_shape[axes]
+    if reshape:
+        # Compute transformed input bounds
+        iy, ix = in_plane_shape
+        out_bounds = rot_matrix @ [[0, 0, iy, iy],
+                                   [0, ix, 0, ix]]
+        # Compute the shape of the transformed input plane
+        out_plane_shape = (out_bounds.ptp(axis=1) + 0.5).astype(cupy.int64)
+    else:
+        out_plane_shape = img_shape[axes]
+
+    out_center = rot_matrix @ ((out_plane_shape - 1) / 2)
+    in_center = (in_plane_shape - 1) / 2
+
+    output_shape = img_shape
+    output_shape[axes] = out_plane_shape
+    output_shape = tuple(output_shape)
+
+    matrix = numpy.identity(ndim)
+    matrix[axes[0], axes[0]] = cos
+    matrix[axes[0], axes[1]] = sin
+    matrix[axes[1], axes[0]] = -sin
+    matrix[axes[1], axes[1]] = cos
+
+    offset = numpy.zeros(ndim, dtype=cupy.float64)
+    offset[axes] = in_center - out_center
+
+    matrix = cupy.asarray(matrix)
+    offset = cupy.asarray(offset)
+
+    return affine_transform(input, matrix, offset, output_shape, output, order,
+                            mode, cval, prefilter)
+
+
+def shift(input, shift, output=None, order=3, mode='constant', cval=0.0,
+          prefilter=True):
+    """Shift an array.
+
+    The array is shifted using spline interpolation of the requested order.
+    Points outside the boundaries of the input are filled according to the
+    given mode.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        shift (float or sequence): The shift along the axes. If a float,
+            ``shift`` is the same for each axis. If a sequence, ``shift``
+            should contain one value for each axis.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+        cval (scalar): Value used for points outside the boundaries of
+            the input if ``mode='constant'`` or ``mode='opencv'``. Default is
+            0.0
+        prefilter (bool): It is not used yet. It just exists for compatibility
+            with :mod:`scipy.ndimage`.
+
+    Returns:
+        cupy.ndarray or None:
+            The shifted input.
+
+    .. seealso:: :func:`scipy.ndimage.shift`
+    """
+
+    _check_parameter('shift', order, mode)
+
+    shift = _util._fix_sequence_arg(shift, input.ndim, 'shift', float)
+
+    if mode == 'opencv':
+        mode = '_opencv_edge'
+
+        output = affine_transform(
+            input,
+            cupy.ones(input.ndim, input.dtype),
+            cupy.negative(cupy.asarray(shift)),
+            None,
+            output,
+            order,
+            mode,
+            cval,
+            prefilter,
+        )
+    else:
+        output = _util._get_output(output, input)
+        if input.dtype.kind in 'iu':
+            input = input.astype(cupy.float32)
+        filtered, nprepad = _filter_input(input, prefilter, mode, cval, order)
+        integer_output = output.dtype.kind in 'iu'
+        _util._check_cval(mode, cval, integer_output)
+        large_int = prod(input.shape) > 1 << 31
+        kern = _interp_kernels._get_shift_kernel(
+            input.ndim, large_int, input.shape, mode, cval=cval, order=order,
+            integer_output=integer_output, nprepad=nprepad)
+        shift = cupy.asarray(shift, dtype=cupy.float64, order='C')
+        if shift.ndim != 1:
+            raise ValueError('shift must be 1d')
+        if shift.size != filtered.ndim:
+            raise ValueError('len(shift) must equal input.ndim')
+        kern(filtered, shift, output)
+    return output
+
+
+def zoom(input, zoom, output=None, order=3, mode='constant', cval=0.0,
+         prefilter=True, *, grid_mode=False):
+    """Zoom an array.
+
+    The array is zoomed using spline interpolation of the requested order.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        zoom (float or sequence): The zoom factor along the axes. If a float,
+            ``zoom`` is the same for each axis. If a sequence, ``zoom`` should
+            contain one value for each axis.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array.
+        order (int): The order of the spline interpolation, default is 3. Must
+            be in the range 0-5.
+        mode (str): Points outside the boundaries of the input are filled
+            according to the given mode (``'constant'``, ``'nearest'``,
+            ``'mirror'``, ``'reflect'``, ``'wrap'``, ``'grid-mirror'``,
+            ``'grid-wrap'``, ``'grid-constant'`` or ``'opencv'``).
+        cval (scalar): Value used for points outside the boundaries of
+            the input if ``mode='constant'`` or ``mode='opencv'``. Default is
+            0.0
+        prefilter (bool): It is not used yet. It just exists for compatibility
+            with :mod:`scipy.ndimage`.
+        grid_mode (bool, optional): If False, the distance from the pixel
+            centers is zoomed. Otherwise, the distance including the full pixel
+            extent is used. For example, a 1d signal of length 5 is considered
+            to have length 4 when ``grid_mode`` is False, but length 5 when
+            ``grid_mode`` is True. See the following visual illustration:
+
+            .. code-block:: text
+
+                    | pixel 1 | pixel 2 | pixel 3 | pixel 4 | pixel 5 |
+                         |<-------------------------------------->|
+                                            vs.
+                    |<----------------------------------------------->|
+
+            The starting point of the arrow in the diagram above corresponds to
+            coordinate location 0 in each mode.
+
+    Returns:
+        cupy.ndarray or None:
+            The zoomed input.
+
+    .. seealso:: :func:`scipy.ndimage.zoom`
+    """
+
+    _check_parameter('zoom', order, mode)
+
+    zoom = _util._fix_sequence_arg(zoom, input.ndim, 'zoom', float)
+
+    output_shape = []
+    for s, z in zip(input.shape, zoom):
+        output_shape.append(int(round(s * z)))
+    output_shape = tuple(output_shape)
+
+    if mode == 'opencv':
+        zoom = []
+        offset = []
+        for in_size, out_size in zip(input.shape, output_shape):
+            if out_size > 1:
+                zoom.append(float(in_size) / out_size)
+                offset.append((zoom[-1] - 1) / 2.0)
+            else:
+                zoom.append(0)
+                offset.append(0)
+        mode = 'nearest'
+
+        output = affine_transform(
+            input,
+            cupy.asarray(zoom),
+            offset,
+            output_shape,
+            output,
+            order,
+            mode,
+            cval,
+            prefilter,
+        )
+    else:
+        if grid_mode:
+
+            # warn about modes that may have surprising behavior
+            suggest_mode = None
+            if mode == 'constant':
+                suggest_mode = 'grid-constant'
+            elif mode == 'wrap':
+                suggest_mode = 'grid-wrap'
+            if suggest_mode is not None:
+                warnings.warn(
+                    f'It is recommended to use mode = {suggest_mode} instead '
+                    f'of {mode} when grid_mode is True.')
+
+        zoom = []
+        for in_size, out_size in zip(input.shape, output_shape):
+            if grid_mode and out_size > 0:
+                zoom.append(in_size / out_size)
+            elif out_size > 1:
+                zoom.append((in_size - 1) / (out_size - 1))
+            else:
+                zoom.append(0)
+
+        output = _util._get_output(output, input, shape=output_shape)
+        if input.dtype.kind in 'iu':
+            input = input.astype(cupy.float32)
+        filtered, nprepad = _filter_input(input, prefilter, mode, cval, order)
+        integer_output = output.dtype.kind in 'iu'
+        _util._check_cval(mode, cval, integer_output)
+        large_int = max(prod(input.shape), prod(output_shape)) > 1 << 31
+        kern = _interp_kernels._get_zoom_kernel(
+            input.ndim, large_int, output_shape, mode, order=order,
+            integer_output=integer_output, grid_mode=grid_mode,
+            nprepad=nprepad)
+        zoom = cupy.asarray(zoom, dtype=cupy.float64)
+        kern(filtered, zoom, output)
+    return output
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_morphology.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_morphology.py
new file mode 100644
index 000000000..5c6462cff
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_morphology.py
@@ -0,0 +1,1016 @@
+import operator
+import warnings
+
+import cupy
+import numpy
+from cupy import _core
+
+from cucim.skimage._vendored import _ndimage_filters as _filters
+from cucim.skimage._vendored import _ndimage_filters_core as _filters_core
+from cucim.skimage._vendored import _ndimage_util as _util
+
+
+@cupy.memoize(for_each_device=True)
+def _get_binary_erosion_kernel(
+    w_shape, int_type, offsets, center_is_true, border_value, invert, masked,
+    all_weights_nonzero
+):
+    if invert:
+        border_value = int(not border_value)
+        true_val = 0
+        false_val = 1
+    else:
+        true_val = 1
+        false_val = 0
+
+    if masked:
+        pre = """
+            bool mv = (bool)mask[i];
+            bool _in = (bool)x[i];
+            if (!mv) {{
+                y = cast<Y>(_in);
+                return;
+            }} else if ({center_is_true} && _in == {false_val}) {{
+                y = cast<Y>(_in);
+                return;
+            }}""".format(center_is_true=int(center_is_true),
+                         false_val=false_val)
+    else:
+        pre = """
+            bool _in = (bool)x[i];
+            if ({center_is_true} && _in == {false_val}) {{
+                y = cast<Y>(_in);
+                return;
+            }}""".format(center_is_true=int(center_is_true),
+                         false_val=false_val)
+    pre = pre + """
+            y = cast<Y>({true_val});""".format(true_val=true_val)
+
+    # {{{{ required because format is called again within _generate_nd_kernel
+    found = """
+        if ({{cond}}) {{{{
+            if (!{border_value}) {{{{
+                y = cast<Y>({false_val});
+                return;
+            }}}}
+        }}}} else {{{{
+            bool nn = {{value}} ? {true_val} : {false_val};
+            if (!nn) {{{{
+                y = cast<Y>({false_val});
+                return;
+            }}}}
+        }}}}""".format(true_val=int(true_val),
+                       false_val=int(false_val),
+                       border_value=int(border_value),)
+
+    name = 'binary_erosion'
+    if false_val:
+        name += '_invert'
+    return _filters_core._generate_nd_kernel(
+        name,
+        pre,
+        found,
+        '',
+        'constant', w_shape, int_type, offsets, 0, ctype='Y', has_weights=True,
+        has_structure=False, has_mask=masked, binary_morphology=True,
+        all_weights_nonzero=all_weights_nonzero)
+
+
+def _center_is_true(structure, origin):
+    coor = tuple([oo + ss // 2 for ss, oo in zip(structure.shape, origin)])
+    return bool(structure[coor])  # device synchronization
+
+
+def iterate_structure(structure, iterations, origin=None):
+    """Iterate a structure by dilating it with itself.
+
+    Args:
+        structure(array_like): Structuring element (an array of bools,
+            for example), to be dilated with itself.
+        iterations(int): The number of dilations performed on the structure
+            with itself.
+        origin(int or tuple of int, optional): If origin is None, only the
+            iterated structure is returned. If not, a tuple of the iterated
+            structure and the modified origin is returned.
+
+    Returns:
+        cupy.ndarray: A new structuring element obtained by dilating
+        ``structure`` (``iterations`` - 1) times with itself.
+
+    .. seealso:: :func:`scipy.ndimage.iterate_structure`
+    """
+    if iterations < 2:
+        return structure.copy()
+    ni = iterations - 1
+    shape = [ii + ni * (ii - 1) for ii in structure.shape]
+    pos = [ni * (structure.shape[ii] // 2) for ii in range(len(shape))]
+    slc = tuple(
+        slice(pos[ii], pos[ii] + structure.shape[ii], None)
+        for ii in range(len(shape))
+    )
+    out = cupy.zeros(shape, bool)
+    out[slc] = structure != 0
+    out = binary_dilation(out, structure, iterations=ni)
+    if origin is None:
+        return out
+    else:
+        origin = _util._fix_sequence_arg(origin, structure.ndim, 'origin', int)
+        origin = [iterations * o for o in origin]
+        return out, origin
+
+
+def generate_binary_structure(rank, connectivity):
+    """Generate a binary structure for binary morphological operations.
+
+    Args:
+        rank(int): Number of dimensions of the array to which the structuring
+            element will be applied, as returned by ``np.ndim``.
+        connectivity(int): ``connectivity`` determines which elements of the
+            output array belong to the structure, i.e., are considered as
+            neighbors of the central element. Elements up to a squared distance
+            of ``connectivity`` from the center are considered neighbors.
+            ``connectivity`` may range from 1 (no diagonal elements are
+            neighbors) to ``rank`` (all elements are neighbors).
+
+    Returns:
+        cupy.ndarray: Structuring element which may be used for binary
+        morphological operations, with ``rank`` dimensions and all
+        dimensions equal to 3.
+
+    .. seealso:: :func:`scipy.ndimage.generate_binary_structure`
+    """
+    if connectivity < 1:
+        connectivity = 1
+    if rank < 1:
+        return cupy.asarray(True, dtype=bool)
+    output = numpy.fabs(numpy.indices([3] * rank) - 1)
+    output = numpy.add.reduce(output, 0)
+    output = output <= connectivity
+    return cupy.asarray(output)
+
+
+def _binary_erosion(input, structure, iterations, mask, output, border_value,
+                    origin, invert, brute_force=True):
+    try:
+        iterations = operator.index(iterations)
+    except TypeError:
+        raise TypeError('iterations parameter should be an integer')
+
+    if input.dtype.kind == 'c':
+        raise TypeError('Complex type not supported')
+    if structure is None:
+        structure = generate_binary_structure(input.ndim, 1)
+        all_weights_nonzero = input.ndim == 1
+        center_is_true = True
+        default_structure = True
+    else:
+        structure = structure.astype(dtype=bool, copy=False)
+        # transfer to CPU for use in determining if it is fully dense
+        # structure_cpu = cupy.asnumpy(structure)
+        default_structure = False
+    if structure.ndim != input.ndim:
+        raise RuntimeError('structure and input must have same dimensionality')
+    if not structure.flags.c_contiguous:
+        structure = cupy.ascontiguousarray(structure)
+    if structure.size < 1:
+        raise RuntimeError('structure must not be empty')
+
+    if mask is not None:
+        if mask.shape != input.shape:
+            raise RuntimeError('mask and input must have equal sizes')
+        if not mask.flags.c_contiguous:
+            mask = cupy.ascontiguousarray(mask)
+        masked = True
+    else:
+        masked = False
+    origin = _util._fix_sequence_arg(origin, input.ndim, 'origin', int)
+
+    if isinstance(output, cupy.ndarray):
+        if output.dtype.kind == 'c':
+            raise TypeError('Complex output type not supported')
+    else:
+        output = bool
+    output = _util._get_output(output, input)
+    temp_needed = cupy.shares_memory(output, input, 'MAY_SHARE_BOUNDS')
+    if temp_needed:
+        # input and output arrays cannot share memory
+        temp = output
+        output = _util._get_output(output.dtype, input)
+    if structure.ndim == 0:
+        # kernel doesn't handle ndim=0, so special case it here
+        if float(structure):
+            output[...] = cupy.asarray(input, dtype=bool)
+        else:
+            output[...] = ~cupy.asarray(input, dtype=bool)
+        return output
+    origin = tuple(origin)
+    int_type = _util._get_inttype(input)
+    offsets = _filters_core._origins_to_offsets(origin, structure.shape)
+    if not default_structure:
+        # synchronize required to determine if all weights are non-zero
+        nnz = int(cupy.count_nonzero(structure))
+        all_weights_nonzero = nnz == structure.size
+        if all_weights_nonzero:
+            center_is_true = True
+        else:
+            center_is_true = _center_is_true(structure, origin)
+
+    erode_kernel = _get_binary_erosion_kernel(
+        structure.shape, int_type, offsets, center_is_true, border_value,
+        invert, masked, all_weights_nonzero,
+    )
+
+    if iterations == 1:
+        if masked:
+            output = erode_kernel(input, structure, mask, output)
+        else:
+            output = erode_kernel(input, structure, output)
+    elif center_is_true and not brute_force:
+        raise NotImplementedError(
+            'only brute_force iteration has been implemented'
+        )
+    else:
+        if cupy.shares_memory(output, input, 'MAY_SHARE_BOUNDS'):
+            raise ValueError('output and input may not overlap in memory')
+        tmp_in = cupy.empty_like(input, dtype=output.dtype)
+        tmp_out = output
+        if iterations >= 1 and not iterations & 1:
+            tmp_in, tmp_out = tmp_out, tmp_in
+        if masked:
+            tmp_out = erode_kernel(input, structure, mask, tmp_out)
+        else:
+            tmp_out = erode_kernel(input, structure, tmp_out)
+        # TODO: kernel doesn't return the changed status, so determine it here
+        changed = not (input == tmp_out).all()  # synchronize!
+        ii = 1
+        while ii < iterations or ((iterations < 1) and changed):
+            tmp_in, tmp_out = tmp_out, tmp_in
+            if masked:
+                tmp_out = erode_kernel(tmp_in, structure, mask, tmp_out)
+            else:
+                tmp_out = erode_kernel(tmp_in, structure, tmp_out)
+            changed = not (tmp_in == tmp_out).all()
+            ii += 1
+            if not changed and (not ii & 1):  # synchronize!
+                # can exit early if nothing changed
+                # (only do this after even number of tmp_in/out swaps)
+                break
+        output = tmp_out
+    if temp_needed:
+        _core.elementwise_copy(output, temp)
+        output = temp
+    return output
+
+
+def binary_erosion(input, structure=None, iterations=1, mask=None, output=None,
+                   border_value=0, origin=0, brute_force=False):
+    """Multidimensional binary erosion with a given structuring element.
+
+    Binary erosion is a mathematical morphology operation used for image
+    processing.
+
+    Args:
+        input(cupy.ndarray): The input binary array_like to be eroded.
+            Non-zero (True) elements form the subset to be eroded.
+        structure(cupy.ndarray, optional): The structuring element used for the
+            erosion. Non-zero elements are considered True. If no structuring
+            element is provided an element is generated with a square
+            connectivity equal to one. (Default value = None).
+        iterations(int, optional): The erosion is repeated ``iterations`` times
+            (one, by default). If iterations is less than 1, the erosion is
+            repeated until the result does not change anymore. Only an integer
+            of iterations is accepted.
+        mask(cupy.ndarray or None, optional): If a mask is given, only those
+            elements with a True value at the corresponding mask element are
+            modified at each iteration. (Default value = None)
+        output(cupy.ndarray, optional): Array of the same shape as input, into
+            which the output is placed. By default, a new array is created.
+        border_value(int (cast to 0 or 1), optional): Value at the
+            border in the output array. (Default value = 0)
+        origin(int or tuple of ints, optional): Placement of the filter, by
+            default 0.
+        brute_force(boolean, optional): Memory condition: if False, only the
+            pixels whose value was changed in the last iteration are tracked as
+            candidates to be updated (eroded) in the current iteration; if
+            True all pixels are considered as candidates for erosion,
+            regardless of what happened in the previous iteration.
+
+    Returns:
+        cupy.ndarray: The result of binary erosion.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_erosion`
+    """
+    return _binary_erosion(input, structure, iterations, mask, output,
+                           border_value, origin, 0, brute_force)
+
+
+def binary_dilation(input, structure=None, iterations=1, mask=None,
+                    output=None, border_value=0, origin=0, brute_force=False):
+    """Multidimensional binary dilation with the given structuring element.
+
+    Args:
+        input(cupy.ndarray): The input binary array_like to be dilated.
+            Non-zero (True) elements form the subset to be dilated.
+        structure(cupy.ndarray, optional): The structuring element used for the
+            dilation. Non-zero elements are considered True. If no structuring
+            element is provided an element is generated with a square
+            connectivity equal to one. (Default value = None).
+        iterations(int, optional): The dilation is repeated ``iterations``
+            times (one, by default). If iterations is less than 1, the dilation
+            is repeated until the result does not change anymore. Only an
+            integer of iterations is accepted.
+        mask(cupy.ndarray or None, optional): If a mask is given, only those
+            elements with a True value at the corresponding mask element are
+            modified at each iteration. (Default value = None)
+        output(cupy.ndarray, optional): Array of the same shape as input, into
+            which the output is placed. By default, a new array is created.
+        border_value(int (cast to 0 or 1), optional): Value at the
+            border in the output array. (Default value = 0)
+        origin(int or tuple of ints, optional): Placement of the filter, by
+            default 0.
+        brute_force(boolean, optional): Memory condition: if False, only the
+            pixels whose value was changed in the last iteration are tracked as
+            candidates to be updated (dilated) in the current iteration; if
+            True all pixels are considered as candidates for dilation,
+            regardless of what happened in the previous iteration.
+
+    Returns:
+        cupy.ndarray: The result of binary dilation.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_dilation`
+    """
+    if structure is None:
+        structure = generate_binary_structure(input.ndim, 1)
+    origin = _util._fix_sequence_arg(origin, input.ndim, 'origin', int)
+    structure = structure[tuple([slice(None, None, -1)] * structure.ndim)]
+    for ii in range(len(origin)):
+        origin[ii] = -origin[ii]
+        if not structure.shape[ii] & 1:
+            origin[ii] -= 1
+    return _binary_erosion(input, structure, iterations, mask, output,
+                           border_value, origin, 1, brute_force)
+
+
+def binary_opening(input, structure=None, iterations=1, output=None, origin=0,
+                   mask=None, border_value=0, brute_force=False):
+    """
+    Multidimensional binary opening with the given structuring element.
+
+    The *opening* of an input image by a structuring element is the
+    *dilation* of the *erosion* of the image by the structuring element.
+
+    Args:
+        input(cupy.ndarray): The input binary array to be opened.
+            Non-zero (True) elements form the subset to be opened.
+        structure(cupy.ndarray, optional): The structuring element used for the
+            opening. Non-zero elements are considered True. If no structuring
+            element is provided an element is generated with a square
+            connectivity equal to one. (Default value = None).
+        iterations(int, optional): The opening is repeated ``iterations`` times
+            (one, by default). If iterations is less than 1, the opening is
+            repeated until the result does not change anymore. Only an integer
+            of iterations is accepted.
+        output(cupy.ndarray, optional): Array of the same shape as input, into
+            which the output is placed. By default, a new array is created.
+        origin(int or tuple of ints, optional): Placement of the filter, by
+            default 0.
+        mask(cupy.ndarray or None, optional): If a mask is given, only those
+            elements with a True value at the corresponding mask element are
+            modified at each iteration. (Default value = None)
+        border_value(int (cast to 0 or 1), optional): Value at the
+            border in the output array. (Default value = 0)
+        brute_force(boolean, optional): Memory condition: if False, only the
+            pixels whose value was changed in the last iteration are tracked as
+            candidates to be updated (dilated) in the current iteration; if
+            True all pixels are considered as candidates for opening,
+            regardless of what happened in the previous iteration.
+
+    Returns:
+        cupy.ndarray: The result of binary opening.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_opening`
+    """
+    if structure is None:
+        rank = input.ndim
+        structure = generate_binary_structure(rank, 1)
+    tmp = binary_erosion(input, structure, iterations, mask, None,
+                         border_value, origin, brute_force)
+    return binary_dilation(tmp, structure, iterations, mask, output,
+                           border_value, origin, brute_force)
+
+
+def binary_closing(input, structure=None, iterations=1, output=None, origin=0,
+                   mask=None, border_value=0, brute_force=False):
+    """
+    Multidimensional binary closing with the given structuring element.
+
+    The *closing* of an input image by a structuring element is the
+    *erosion* of the *dilation* of the image by the structuring element.
+
+    Args:
+        input(cupy.ndarray): The input binary array to be closed.
+            Non-zero (True) elements form the subset to be closed.
+        structure(cupy.ndarray, optional): The structuring element used for the
+            closing. Non-zero elements are considered True. If no structuring
+            element is provided an element is generated with a square
+            connectivity equal to one. (Default value = None).
+        iterations(int, optional): The closing is repeated ``iterations`` times
+            (one, by default). If iterations is less than 1, the closing is
+            repeated until the result does not change anymore. Only an integer
+            of iterations is accepted.
+        output(cupy.ndarray, optional): Array of the same shape as input, into
+            which the output is placed. By default, a new array is created.
+        origin(int or tuple of ints, optional): Placement of the filter, by
+            default 0.
+        mask(cupy.ndarray or None, optional): If a mask is given, only those
+            elements with a True value at the corresponding mask element are
+            modified at each iteration. (Default value = None)
+        border_value(int (cast to 0 or 1), optional): Value at the
+            border in the output array. (Default value = 0)
+        brute_force(boolean, optional): Memory condition: if False, only the
+            pixels whose value was changed in the last iteration are tracked as
+            candidates to be updated (dilated) in the current iteration; if
+            True all pixels are considered as candidates for closing,
+            regardless of what happened in the previous iteration.
+
+    Returns:
+        cupy.ndarray: The result of binary closing.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_closing`
+    """
+    if structure is None:
+        rank = input.ndim
+        structure = generate_binary_structure(rank, 1)
+    tmp = binary_dilation(input, structure, iterations, mask, None,
+                          border_value, origin, brute_force)
+    return binary_erosion(tmp, structure, iterations, mask, output,
+                          border_value, origin, brute_force)
+
+
+def binary_hit_or_miss(input, structure1=None, structure2=None, output=None,
+                       origin1=0, origin2=None):
+    """
+    Multidimensional binary hit-or-miss transform.
+
+    The hit-or-miss transform finds the locations of a given pattern
+    inside the input image.
+
+    Args:
+        input (cupy.ndarray): Binary image where a pattern is to be detected.
+        structure1 (cupy.ndarray, optional): Part of the structuring element to
+            be fitted to the foreground (non-zero elements) of ``input``. If no
+            value is provided, a structure of square connectivity 1 is chosen.
+        structure2 (cupy.ndarray, optional): Second part of the structuring
+            element that has to miss completely the foreground. If no value is
+            provided, the complementary of ``structure1`` is taken.
+        output (cupy.ndarray, dtype or None, optional): Array of the same shape
+            as input, into which the output is placed. By default, a new array
+            is created.
+        origin1 (int or tuple of ints, optional): Placement of the first part
+            of the structuring element ``structure1``, by default 0 for a
+            centered structure.
+        origin2 (int or tuple of ints or None, optional): Placement of the
+            second part of the structuring element ``structure2``, by default 0
+            for a centered structure. If a value is provided for ``origin1``
+            and not for ``origin2``, then ``origin2`` is set to ``origin1``.
+
+    Returns:
+        cupy.ndarray: Hit-or-miss transform of ``input`` with the given
+        structuring element (``structure1``, ``structure2``).
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_hit_or_miss`
+    """
+    if structure1 is None:
+        structure1 = generate_binary_structure(input.ndim, 1)
+    if structure2 is None:
+        structure2 = cupy.logical_not(structure1)
+    origin1 = _util._fix_sequence_arg(origin1, input.ndim, 'origin1', int)
+    if origin2 is None:
+        origin2 = origin1
+    else:
+        origin2 = _util._fix_sequence_arg(origin2, input.ndim, 'origin2', int)
+
+    tmp1 = _binary_erosion(input, structure1, 1, None, None, 0, origin1, 0,
+                           False)
+    inplace = isinstance(output, cupy.ndarray)
+    result = _binary_erosion(input, structure2, 1, None, output, 0, origin2, 1,
+                             False)
+    if inplace:
+        cupy.logical_not(output, output)
+        cupy.logical_and(tmp1, output, output)
+    else:
+        cupy.logical_not(result, result)
+        return cupy.logical_and(tmp1, result)
+
+
+def binary_propagation(input, structure=None, mask=None, output=None,
+                       border_value=0, origin=0):
+    """
+    Multidimensional binary propagation with the given structuring element.
+
+    Args:
+        input (cupy.ndarray): Binary image to be propagated inside ``mask``.
+        structure (cupy.ndarray, optional): Structuring element used in the
+            successive dilations. The output may depend on the structuring
+            element, especially if ``mask`` has several connex components. If
+            no structuring element is provided, an element is generated with a
+            squared connectivity equal to one.
+        mask (cupy.ndarray, optional): Binary mask defining the region into
+            which ``input`` is allowed to propagate.
+        output (cupy.ndarray, optional): Array of the same shape as input, into
+            which the output is placed. By default, a new array is created.
+        border_value (int, optional): Value at the border in the output array.
+            The value is cast to 0 or 1.
+        origin (int or tuple of ints, optional): Placement of the filter.
+
+    Returns:
+        cupy.ndarray : Binary propagation of ``input`` inside ``mask``.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_propagation`
+    """
+    return binary_dilation(input, structure, -1, mask, output, border_value,
+                           origin, brute_force=True)
+
+
+def binary_fill_holes(input, structure=None, output=None, origin=0):
+    """Fill the holes in binary objects.
+
+    Args:
+        input (cupy.ndarray): N-D binary array with holes to be filled.
+        structure (cupy.ndarray, optional):  Structuring element used in the
+            computation; large-size elements make computations faster but may
+            miss holes separated from the background by thin regions. The
+            default element (with a square connectivity equal to one) yields
+            the intuitive result where all holes in the input have been filled.
+        output (cupy.ndarray, dtype or None, optional): Array of the same shape
+            as input, into which the output is placed. By default, a new array
+            is created.
+        origin (int, tuple of ints, optional): Position of the structuring
+            element.
+
+    Returns:
+        cupy.ndarray: Transformation of the initial image ``input`` where holes
+        have been filled.
+
+    .. warning::
+
+        This function may synchronize the device.
+
+    .. seealso:: :func:`scipy.ndimage.binary_fill_holes`
+    """
+    mask = cupy.logical_not(input)
+    tmp = cupy.zeros(mask.shape, bool)
+    inplace = isinstance(output, cupy.ndarray)
+    # TODO (grlee77): set brute_force=False below once implemented
+    if inplace:
+        binary_dilation(tmp, structure, -1, mask, output, 1, origin,
+                        brute_force=True)
+        cupy.logical_not(output, output)
+    else:
+        output = binary_dilation(tmp, structure, -1, mask, None, 1, origin,
+                                 brute_force=True)
+        cupy.logical_not(output, output)
+        return output
+
+
+def grey_erosion(input, size=None, footprint=None, structure=None, output=None,
+                 mode='reflect', cval=0.0, origin=0):
+    """Calculates a greyscale erosion.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the greyscale erosion. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for greyscale erosion. Non-zero values
+            give the set of neighbors of the center over which minimum is
+            chosen.
+        structure (array of ints): Structuring element used for the greyscale
+            erosion. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of greyscale erosion.
+
+    .. seealso:: :func:`scipy.ndimage.grey_erosion`
+    """
+
+    if size is None and footprint is None and structure is None:
+        raise ValueError('size, footprint or structure must be specified')
+
+    return _filters._min_or_max_filter(input, size, footprint, structure,
+                                       output, mode, cval, origin, 'min')
+
+
+def grey_dilation(input, size=None, footprint=None, structure=None,
+                  output=None, mode='reflect', cval=0.0, origin=0):
+    """Calculates a greyscale dilation.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the greyscale dilation. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for greyscale dilation. Non-zero values
+            give the set of neighbors of the center over which maximum is
+            chosen.
+        structure (array of ints): Structuring element used for the greyscale
+            dilation. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of greyscale dilation.
+
+    .. seealso:: :func:`scipy.ndimage.grey_dilation`
+    """
+
+    if size is None and footprint is None and structure is None:
+        raise ValueError('size, footprint or structure must be specified')
+    if structure is not None:
+        structure = cupy.array(structure)
+        structure = structure[tuple([slice(None, None, -1)] * structure.ndim)]
+    if footprint is not None:
+        footprint = cupy.array(footprint)
+        footprint = footprint[tuple([slice(None, None, -1)] * footprint.ndim)]
+
+    origin = _util._fix_sequence_arg(origin, input.ndim, 'origin', int)
+    for i in range(len(origin)):
+        origin[i] = -origin[i]
+        if footprint is not None:
+            sz = footprint.shape[i]
+        elif structure is not None:
+            sz = structure.shape[i]
+        elif numpy.isscalar(size):
+            sz = size
+        else:
+            sz = size[i]
+        if sz % 2 == 0:
+            origin[i] -= 1
+
+    return _filters._min_or_max_filter(input, size, footprint, structure,
+                                       output, mode, cval, origin, 'max')
+
+
+def grey_closing(input, size=None, footprint=None, structure=None,
+                 output=None, mode='reflect', cval=0.0, origin=0):
+    """Calculates a multi-dimensional greyscale closing.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the greyscale closing. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for greyscale closing. Non-zero values
+            give the set of neighbors of the center over which closing is
+            chosen.
+        structure (array of ints): Structuring element used for the greyscale
+            closing. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of greyscale closing.
+
+    .. seealso:: :func:`scipy.ndimage.grey_closing`
+    """
+    if (size is not None) and (footprint is not None):
+        warnings.warn('ignoring size because footprint is set', UserWarning,
+                      stacklevel=2)
+    tmp = grey_dilation(input, size, footprint, structure, None, mode, cval,
+                        origin)
+    return grey_erosion(tmp, size, footprint, structure, output, mode, cval,
+                        origin)
+
+
+def grey_opening(input, size=None, footprint=None, structure=None,
+                 output=None, mode='reflect', cval=0.0, origin=0):
+    """Calculates a multi-dimensional greyscale opening.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the greyscale opening. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for greyscale opening. Non-zero values
+            give the set of neighbors of the center over which opening is
+            chosen.
+        structure (array of ints): Structuring element used for the greyscale
+            opening. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The result of greyscale opening.
+
+    .. seealso:: :func:`scipy.ndimage.grey_opening`
+    """
+    if (size is not None) and (footprint is not None):
+        warnings.warn('ignoring size because footprint is set', UserWarning,
+                      stacklevel=2)
+    tmp = grey_erosion(input, size, footprint, structure, None, mode, cval,
+                       origin)
+    return grey_dilation(tmp, size, footprint, structure, output, mode, cval,
+                         origin)
+
+
+def morphological_gradient(
+    input,
+    size=None,
+    footprint=None,
+    structure=None,
+    output=None,
+    mode='reflect',
+    cval=0.0,
+    origin=0,
+):
+    """
+    Multidimensional morphological gradient.
+
+    The morphological gradient is calculated as the difference between a
+    dilation and an erosion of the input with a given structuring element.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the morphological gradient. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for morphological gradient. Non-zero
+            values give the set of neighbors of the center over which opening
+            is chosen.
+        structure (array of ints): Structuring element used for the
+            morphological gradient. ``structure`` may be a non-flat
+            structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The morphological gradient of the input.
+
+    .. seealso:: :func:`scipy.ndimage.morphological_gradient`
+    """
+    tmp = grey_dilation(
+        input, size, footprint, structure, None, mode, cval, origin
+    )
+    if isinstance(output, cupy.ndarray):
+        grey_erosion(
+            input, size, footprint, structure, output, mode, cval, origin
+        )
+        return cupy.subtract(tmp, output, output)
+    else:
+        return tmp - grey_erosion(
+            input, size, footprint, structure, None, mode, cval, origin
+        )
+
+
+def morphological_laplace(
+    input,
+    size=None,
+    footprint=None,
+    structure=None,
+    output=None,
+    mode='reflect',
+    cval=0.0,
+    origin=0,
+):
+    """
+    Multidimensional morphological laplace.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the morphological laplace. Optional if ``footprint`` or
+            ``structure`` is provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for morphological laplace. Non-zero
+            values give the set of neighbors of the center over which opening
+            is chosen.
+        structure (array of ints): Structuring element used for the
+            morphological laplace. ``structure`` may be a non-flat
+            structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: The morphological laplace of the input.
+
+    .. seealso:: :func:`scipy.ndimage.morphological_laplace`
+    """
+    tmp1 = grey_dilation(
+        input, size, footprint, structure, None, mode, cval, origin
+    )
+    if isinstance(output, cupy.ndarray):
+        grey_erosion(
+            input, size, footprint, structure, output, mode, cval, origin
+        )
+        cupy.add(tmp1, output, output)
+        cupy.subtract(output, input, output)
+        return cupy.subtract(output, input, output)
+    else:
+        tmp2 = grey_erosion(
+            input, size, footprint, structure, None, mode, cval, origin
+        )
+        cupy.add(tmp1, tmp2, tmp2)
+        cupy.subtract(tmp2, input, tmp2)
+        cupy.subtract(tmp2, input, tmp2)
+        return tmp2
+
+
+def white_tophat(
+    input,
+    size=None,
+    footprint=None,
+    structure=None,
+    output=None,
+    mode='reflect',
+    cval=0.0,
+    origin=0,
+):
+    """
+    Multidimensional white tophat filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the white tophat. Optional if ``footprint`` or ``structure`` is
+            provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for the white tophat. Non-zero values
+            give the set of neighbors of the center over which opening is
+            chosen.
+        structure (array of ints): Structuring element used for the white
+            tophat. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarray: Result of the filter of ``input`` with ``structure``.
+
+    .. seealso:: :func:`scipy.ndimage.white_tophat`
+    """
+    if (size is not None) and (footprint is not None):
+        warnings.warn(
+            'ignoring size because footprint is set', UserWarning, stacklevel=2
+        )
+    tmp = grey_erosion(
+        input, size, footprint, structure, None, mode, cval, origin
+    )
+    tmp = grey_dilation(
+        tmp, size, footprint, structure, output, mode, cval, origin
+    )
+    if input.dtype == numpy.bool_ and tmp.dtype == numpy.bool_:
+        cupy.bitwise_xor(input, tmp, out=tmp)
+    else:
+        cupy.subtract(input, tmp, out=tmp)
+    return tmp
+
+
+def black_tophat(
+    input,
+    size=None,
+    footprint=None,
+    structure=None,
+    output=None,
+    mode='reflect',
+    cval=0.0,
+    origin=0,
+):
+    """
+    Multidimensional black tophat filter.
+
+    Args:
+        input (cupy.ndarray): The input array.
+        size (tuple of ints): Shape of a flat and full structuring element used
+            for the black tophat. Optional if ``footprint`` or ``structure`` is
+            provided.
+        footprint (array of ints): Positions of non-infinite elements of a flat
+            structuring element used for the black tophat. Non-zero values
+            give the set of neighbors of the center over which opening is
+            chosen.
+        structure (array of ints): Structuring element used for the black
+            tophat. ``structure`` may be a non-flat structuring element.
+        output (cupy.ndarray, dtype or None): The array in which to place the
+            output.
+        mode (str): The array borders are handled according to the given mode
+            (``'reflect'``, ``'constant'``, ``'nearest'``, ``'mirror'``,
+            ``'wrap'``). Default is ``'reflect'``.
+        cval (scalar): Value to fill past edges of input if mode is
+            ``constant``. Default is ``0.0``.
+        origin (scalar or tuple of scalar): The origin parameter controls the
+            placement of the filter, relative to the center of the current
+            element of the input. Default of 0 is equivalent to
+            ``(0,)*input.ndim``.
+
+    Returns:
+        cupy.ndarry : Result of the filter of ``input`` with ``structure``.
+
+    .. seealso:: :func:`scipy.ndimage.black_tophat`
+    """
+    if (size is not None) and (footprint is not None):
+        warnings.warn(
+            'ignoring size because footprint is set', UserWarning, stacklevel=2
+        )
+    tmp = grey_dilation(
+        input, size, footprint, structure, None, mode, cval, origin
+    )
+    tmp = grey_erosion(
+        tmp, size, footprint, structure, output, mode, cval, origin
+    )
+    if input.dtype == numpy.bool_ and tmp.dtype == numpy.bool_:
+        cupy.bitwise_xor(tmp, input, out=tmp)
+    else:
+        cupy.subtract(tmp, input, out=tmp)
+    return tmp
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_kernel_weights.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_kernel_weights.py
new file mode 100644
index 000000000..b2fc84449
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_kernel_weights.py
@@ -0,0 +1,73 @@
+"""Determination of spline kernel weights (adapted from SciPy)
+
+See more verbose comments for each case there:
+https://github.com/scipy/scipy/blob/eba29d69846ab1299976ff4af71c106188397ccc/scipy/ndimage/src/ni_splines.c#L7   # NOQA
+
+``spline_weights_inline`` is a dict where the key is the spline order and the
+value is the spline weight initialization code.
+"""
+
+spline_weights_inline = {}
+
+# Note: This order = 1 case is currently unused (order = 1 has a different code
+# path in _interp_kernels.py). I think that existing code is a bit more
+# efficient.
+spline_weights_inline[1] = '''
+wx = c_{j} - floor({order} & 1 ? c_{j} : c_{j} + 0.5);
+weights_{j}[0] = 1.0 - wx;
+weights_{j}[1] = wx;
+'''
+
+spline_weights_inline[2] = '''
+wx = c_{j} - floor({order} & 1 ? c_{j} : c_{j} + 0.5);
+weights_{j}[1] = 0.75 - wx * wx;
+wy = 0.5 - wx;
+weights_{j}[0] = 0.5 * wy * wy;
+weights_{j}[2] = 1.0 - weights_{j}[0] - weights_{j}[1];
+'''
+
+spline_weights_inline[3] = '''
+wx = c_{j} - floor({order} & 1 ? c_{j} : c_{j} + 0.5);
+wy = 1.0 - wx;
+weights_{j}[1] = (wx * wx * (wx - 2.0) * 3.0 + 4.0) / 6.0;
+weights_{j}[2] = (wy * wy * (wy - 2.0) * 3.0 + 4.0) / 6.0;
+weights_{j}[0] = wy * wy * wy / 6.0;
+weights_{j}[3] = 1.0 - weights_{j}[0] - weights_{j}[1] - weights_{j}[2];
+'''
+
+spline_weights_inline[4] = '''
+wx = c_{j} - floor({order} & 1 ? c_{j} : c_{j} + 0.5);
+wy = wx * wx;
+weights_{j}[2] = wy * (wy * 0.25 - 0.625) + 115.0 / 192.0;
+wy = 1.0 + wx;
+weights_{j}[1] = wy * (wy * (wy * (5.0 - wy) / 6.0 - 1.25) + 5.0 / 24.0) +
+             55.0 / 96.0;
+wy = 1.0 - wx;
+weights_{j}[3] = wy * (wy * (wy * (5.0 - wy) / 6.0 - 1.25) + 5.0 / 24.0) +
+             55.0 / 96.0;
+wy = 0.5 - wx;
+wy = wy * wy;
+weights_{j}[0] = wy * wy / 24.0;
+weights_{j}[4] = 1.0 - weights_{j}[0] - weights_{j}[1]
+                     - weights_{j}[2] - weights_{j}[3];
+'''
+
+spline_weights_inline[5] = '''
+wx = c_{j} - floor({order} & 1 ? c_{j} : c_{j} + 0.5);
+wy = wx * wx;
+weights_{j}[2] = wy * (wy * (0.25 - wx / 12.0) - 0.5) + 0.55;
+wy = 1.0 - wx;
+wy = wy * wy;
+weights_{j}[3] = wy * (wy * (0.25 - (1.0 - wx) / 12.0) - 0.5) + 0.55;
+wy = wx + 1.0;
+weights_{j}[1] = wy * (wy * (wy * (wy * (wy / 24.0 - 0.375) + 1.25) - 1.75)
+                   + 0.625) + 0.425;
+wy = 2.0 - wx;
+weights_{j}[4] = wy * (wy * (wy * (wy * (wy / 24.0 - 0.375) + 1.25) - 1.75)
+                  + 0.625) + 0.425;
+wy = 1.0 - wx;
+wy = wy * wy;
+weights_{j}[0] = (1.0 - wx) * wy * wy / 120.0;
+weights_{j}[5] = 1.0 - weights_{j}[0] - weights_{j}[1] - weights_{j}[2]
+                     - weights_{j}[3] - weights_{j}[4];
+'''
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_prefilter_core.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_prefilter_core.py
new file mode 100644
index 000000000..c44df836e
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_spline_prefilter_core.py
@@ -0,0 +1,256 @@
+"""
+Spline poles and boundary handling implemented as in SciPy
+
+https://github.com/scipy/scipy/blob/ee6ae72f83a0995aeb34929aed881d3f36fccfda/scipy/ndimage/src/ni_splines.c  # noqa
+"""
+import functools
+import math
+import operator
+import textwrap
+
+import cupy
+
+
+def get_poles(order):
+    if order == 2:
+        # sqrt(8.0) - 3.0
+        return (-0.171572875253809902396622551580603843,)
+    elif order == 3:
+        # sqrt(3.0) - 2.0
+        return (-0.267949192431122706472553658494127633,)
+    elif order == 4:
+        # sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0
+        # sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0
+        return (-0.361341225900220177092212841325675255,
+                -0.013725429297339121360331226939128204)
+    elif order == 5:
+        # sqrt(67.5 - sqrt(4436.25)) + sqrt(26.25) - 6.5
+        # sqrt(67.5 + sqrt(4436.25)) - sqrt(26.25) - 6.5
+        return (-0.430575347099973791851434783493520110,
+                -0.043096288203264653822712376822550182)
+    else:
+        raise ValueError('only order 2-5 supported')
+
+
+def get_gain(poles):
+    return functools.reduce(operator.mul,
+                            [(1.0 - z) * (1.0 - 1.0 / z) for z in poles])
+
+
+def _causal_init_code(mode):
+    """Code for causal initialization step of IIR filtering.
+
+    c is a 1d array of length n and z is a filter pole
+    """
+    code = f'''
+        // causal init for mode={mode}'''
+    if mode == 'mirror':
+        code += '''
+        z_i = z;
+        z_n_1 = pow(z, (P)(n - 1));
+
+        c[0] = c[0] + z_n_1 * c[(n - 1) * element_stride];
+        for (i = 1; i < min(n - 1, static_cast<idx_t>({n_boundary})); ++i) {{
+            c[0] += z_i * (c[i * element_stride] +
+                           z_n_1 * c[(n - 1 - i) * element_stride]);
+            z_i *= z;
+        }}
+        c[0] /= 1 - z_n_1 * z_n_1;'''
+    elif mode == 'grid-wrap':
+        code += '''
+        z_i = z;
+
+        for (i = 1; i < min(n, static_cast<idx_t>({n_boundary})); ++i) {{
+            c[0] += z_i * c[(n - i) * element_stride];
+            z_i *= z;
+        }}
+        c[0] /= 1 - z_i; /* z_i = pow(z, n) */'''
+    elif mode == 'reflect':
+        code += '''
+        z_i = z;
+        z_n = pow(z, (P)n);
+        c0 = c[0];
+
+        c[0] = c[0] + z_n * c[(n - 1) * element_stride];
+        for (i = 1; i < min(n, static_cast<idx_t>({n_boundary})); ++i) {{
+            c[0] += z_i * (c[i * element_stride] +
+                           z_n * c[(n - 1 - i) * element_stride]);
+            z_i *= z;
+        }}
+        c[0] *= z / (1 - z_n * z_n);
+        c[0] += c0;'''
+    else:
+        raise ValueError('invalid mode: {}'.format(mode))
+    return code
+
+
+def _anticausal_init_code(mode):
+    """Code for the anti-causal initialization step of IIR filtering.
+
+    c is a 1d array of length n and z is a filter pole
+    """
+    code = f'''
+        // anti-causal init for mode={mode}'''
+    if mode == 'mirror':
+        code += '''
+        c[(n - 1) * element_stride] = (
+            z * c[(n - 2) * element_stride] +
+            c[(n - 1) * element_stride]) * z / (z * z - 1);'''
+    elif mode == 'grid-wrap':
+        code += '''
+        z_i = z;
+
+        for (i = 0; i < min(n - 1, static_cast<idx_t>({n_boundary})); ++i) {{
+            c[(n - 1) * element_stride] += z_i * c[i * element_stride];
+            z_i *= z;
+        }}
+        c[(n - 1) * element_stride] *= z / (z_i - 1); /* z_i = pow(z, n) */'''
+    elif mode == 'reflect':
+        code += '''
+        c[(n - 1) * element_stride] *= z / (z - 1);'''
+    else:
+        raise ValueError('invalid mode: {}'.format(mode))
+    return code
+
+
+def _get_spline_mode(mode):
+    """spline boundary mode for interpolation with order >= 2."""
+    if mode in ['mirror', 'reflect', 'grid-wrap']:
+        # exact analytic boundary conditions exist for these modes.
+        return mode
+    elif mode == 'grid-mirror':
+        # grid-mirror is a synonym for 'reflect'
+        return 'reflect'
+    # No exact analytical spline boundary condition implemented. Reflect gives
+    # lower error than using mirror or wrap for mode 'nearest'. Otherwise, a
+    # mirror spline boundary condition is used.
+    return 'reflect' if mode == 'nearest' else 'mirror'
+
+
+def _get_spline1d_code(mode, poles, n_boundary):
+    """Generates the code required for IIR filtering of a single 1d signal.
+
+    Prefiltering is done by causal filtering followed by anti-causal filtering.
+    Multiple boundary conditions have been implemented.
+    """
+    code = ['''
+    __device__ void spline_prefilter1d(
+        T* __restrict__ c, idx_t signal_length, idx_t element_stride)
+    {{''']
+
+    # variables common to all boundary modes
+    code.append('''
+        idx_t i, n = signal_length;
+        P z, z_i;''')
+
+    # retrieve the spline boundary extension mode to use
+    mode = _get_spline_mode(mode)
+
+    if mode == 'mirror':
+        # variables specific to mirror boundary mode
+        code.append('''
+        P z_n_1;''')
+    elif mode == 'reflect':
+        # variables specific to reflect boundary mode
+        code.append('''
+        P z_n;
+        T c0;''')
+
+    for pole in poles:
+
+        code.append(f'''
+        // select the current pole
+        z = {pole};''')
+
+        # initialize and apply the causal filter
+        code.append(_causal_init_code(mode))
+        code.append('''
+        // apply the causal filter for the current pole
+        for (i = 1; i < n; ++i) {{
+            c[i * element_stride] += z * c[(i - 1) * element_stride];
+        }}''')
+        # initialize and apply the anti-causal filter
+        code.append(_anticausal_init_code(mode))
+        code.append('''
+        // apply the anti-causal filter for the current pole
+        for (i = n - 2; i >= 0; --i) {{
+            c[i * element_stride] = z * (c[(i + 1) * element_stride] -
+                                         c[i * element_stride]);
+        }}''')
+
+    code += ['''
+    }}''']
+    return textwrap.dedent('\n'.join(code)).format(n_boundary=n_boundary)
+
+
+_FILTER_GENERAL = '''
+#include "cupy/carray.cuh"
+#include "cupy/complex.cuh"
+typedef {data_type} T;
+typedef {pole_type} P;
+typedef {index_type} idx_t;
+template <typename T>
+__device__ T* row(
+        T* ptr, idx_t i, idx_t axis, idx_t ndim, const idx_t* shape) {{
+    idx_t index = 0, stride = 1;
+    for (idx_t a = ndim - 1; a > 0; --a) {{
+        if (a != axis) {{
+            index += (i % shape[a]) * stride;
+            i /= shape[a];
+        }}
+        stride *= shape[a];
+    }}
+    return ptr + index + stride * i;
+}}
+'''
+
+
+_batch_spline1d_strided_template = """
+extern "C" __global__
+__launch_bounds__({block_size})
+void {kernel_name}(T* __restrict__ y, const idx_t* __restrict__ info) {{
+    const idx_t n_signals = info[0], n_samples = info[1],
+        * __restrict__ shape = info+2;
+    idx_t y_elem_stride = 1;
+    for (int a = {ndim} - 1; a > {axis}; --a) {{ y_elem_stride *= shape[a]; }}
+    idx_t unraveled_idx = blockDim.x * blockIdx.x + threadIdx.x;
+    idx_t batch_idx = unraveled_idx;
+    if (batch_idx < n_signals)
+    {{
+        T* __restrict__ y_i = row(y, batch_idx, {axis}, {ndim}, shape);
+        spline_prefilter1d(y_i, n_samples, y_elem_stride);
+    }}
+}}
+"""
+
+
+@cupy.memoize(for_each_device=True)
+def get_raw_spline1d_kernel(axis, ndim, mode, order, index_type='int',
+                            data_type='double', pole_type='double',
+                            block_size=128):
+    """Generate a kernel for applying a spline prefilter along a given axis."""
+    poles = get_poles(order)
+
+    # determine number of samples for the boundary approximation
+    # (SciPy uses n_boundary = n_samples but this is excessive)
+    largest_pole = max([abs(p) for p in poles])
+    # tol < 1e-7 fails test cases comparing to SciPy at atol = rtol = 1e-5
+    tol = 1e-10 if pole_type == 'float' else 1e-18
+    n_boundary = math.ceil(math.log(tol, largest_pole))
+
+    # headers and general utility function for extracting rows of data
+    code = _FILTER_GENERAL.format(index_type=index_type,
+                                  data_type=data_type,
+                                  pole_type=pole_type)
+
+    # generate source for a 1d function for a given boundary mode and poles
+    code += _get_spline1d_code(mode, poles, n_boundary)
+
+    # generate code handling batch operation of the 1d filter
+    mode_str = mode.replace('-', '_')  # cannot have '-' in kernel name
+    kernel_name = (f'cupyx_scipy_ndimage_spline_filter_{ndim}d_ord{order}_'
+                   f'axis{axis}_{mode_str}')
+    code += _batch_spline1d_strided_template.format(ndim=ndim, axis=axis,
+                                                    block_size=block_size,
+                                                    kernel_name=kernel_name)
+    return cupy.RawKernel(code, kernel_name)
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
index 1d87b3256..1a6a830a5 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
@@ -5,6 +5,33 @@
 import numpy
 
 
+def _is_integer_output(output, input):
+    if output is None:
+        return input.dtype.kind in 'iu'
+    elif isinstance(output, cupy.ndarray):
+        return output.dtype.kind in 'iu'
+    return cupy.dtype(output).kind in 'iu'
+
+
+def _check_cval(mode, cval, integer_output):
+    if mode == 'constant' and integer_output and not cupy.isfinite(cval):
+        raise NotImplementedError("Non-finite cval is not supported for "
+                                  "outputs with integer dtype.")
+
+
+def _get_weights_dtype(input, weights, use_cucim_casting=False):
+    if weights.dtype.kind == "c" or input.dtype.kind == "c":
+        return cupy.promote_types(input.real.dtype, cupy.complex64)
+    elif weights.dtype.kind in 'iub':
+        if use_cucim_casting:
+            from cucim.skimage._shared.utils import _supported_float_type
+            return _supported_float_type(weights.dtype)
+        else:
+            # convert integer dtype weights to double as in SciPy
+            return cupy.float64
+    return cupy.promote_types(input.real.dtype, cupy.float32)
+
+
 def _get_output(output, input, shape=None, complex_output=False):
     shape = input.shape if shape is None else shape
     if output is None:
@@ -12,17 +39,17 @@ def _get_output(output, input, shape=None, complex_output=False):
             _dtype = cupy.promote_types(input.dtype, cupy.complex64)
         else:
             _dtype = input.dtype
-        output = cupy.zeros(shape, dtype=_dtype)
+        output = cupy.empty(shape, dtype=_dtype)
     elif isinstance(output, (type, cupy.dtype)):
         if complex_output and cupy.dtype(output).kind != 'c':
             warnings.warn("promoting specified output dtype to complex")
             output = cupy.promote_types(output, cupy.complex64)
-        output = cupy.zeros(shape, dtype=output)
+        output = cupy.empty(shape, dtype=output)
     elif isinstance(output, str):
         output = numpy.sctypeDict[output]
         if complex_output and cupy.dtype(output).kind != 'c':
             raise RuntimeError("output must have complex dtype")
-        output = cupy.zeros(shape, dtype=output)
+        output = cupy.empty(shape, dtype=output)
     elif output.shape != shape:
         raise RuntimeError("output shape not correct")
     elif complex_output and output.dtype.kind != 'c':
@@ -44,6 +71,13 @@ def _fix_sequence_arg(arg, ndim, name, conv=lambda x: x):
     return lst
 
 
+def _check_origin(origin, width):
+    origin = int(origin)
+    if (width // 2 + origin < 0) or (width // 2 + origin >= width):
+        raise ValueError('invalid origin')
+    return origin
+
+
 def _check_mode(mode):
     if mode not in ('reflect', 'constant', 'nearest', 'mirror', 'wrap',
                     'grid-mirror', 'grid-wrap', 'grid-reflect'):
@@ -52,13 +86,6 @@ def _check_mode(mode):
     return mode
 
 
-def _check_origin(origin, width):
-    origin = int(origin)
-    if (width // 2 + origin < 0) or (width // 2 + origin >= width):
-        raise ValueError('invalid origin')
-    return origin
-
-
 def _get_inttype(input):
     # The integer type to use for indices in the input array
     # The indices actually use byte positions and we can't just use
@@ -70,53 +97,118 @@ def _get_inttype(input):
 
 
 def _generate_boundary_condition_ops(mode, ix, xsize, int_t="int",
-                                     float_ix=False):
+                                     float_ix=False, separate=False):
+    """Generate boundary conditions
+
+    If separate = True, a pair of conditions for the (lower, upper) boundary
+    are provided instead of a single expression.
+    """
     min_func = "fmin" if float_ix else "min"
     max_func = "fmax" if float_ix else "max"
     if mode in ['reflect', 'grid-mirror']:
-        ops = '''
-        if ({ix} < 0) {{
-            {ix} = - 1 -{ix};
-        }}
-        {ix} %= {xsize} * 2;
-        {ix} = {min}({ix}, 2 * {xsize} - 1 - {ix});'''.format(
-            ix=ix, xsize=xsize, min=min_func)
-    elif mode == 'mirror':
-        ops = '''
-        if ({xsize} == 1) {{
-            {ix} = 0;
-        }} else {{
+        if separate:
+            ops_upper = f'''
+            {ix} %= {xsize} * 2;
+            {ix} = {min_func}({ix}, 2 * {xsize} - 1 - {ix});
+            '''
+            ops_lower = f'''
             if ({ix} < 0) {{
-                {ix} = -{ix};
+                {ix} = - 1 -{ix};
             }}
-            {ix} = 1 + ({ix} - 1) % (({xsize} - 1) * 2);
-            {ix} = {min}({ix}, 2 * {xsize} - 2 - {ix});
-        }}'''.format(ix=ix, xsize=xsize, min=min_func)
+            ''' + ops_upper
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''
+            if ({ix} < 0) {{
+                {ix} = - 1 -{ix};
+            }}
+            {ix} %= {xsize} * 2;
+            {ix} = {min_func}({ix}, 2 * {xsize} - 1 - {ix});'''
+    elif mode == 'mirror':
+        if separate:
+            temp1 = f'''
+            if ({xsize} == 1) {{
+                {ix} = 0;
+            }} else {{
+            '''
+            temp2 = f'''
+                if ({ix} < 0) {{
+                    {ix} = -{ix};
+                }}
+            '''
+            temp3 = f'''
+                {ix} = 1 + ({ix} - 1) % (({xsize} - 1) * 2);
+                {ix} = {min_func}({ix}, 2 * {xsize} - 2 - {ix});
+            }}'''
+            ops_lower = temp1 + temp2 + temp3
+            ops_upper = temp1 + temp3
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''
+            if ({xsize} == 1) {{
+                {ix} = 0;
+            }} else {{
+                if ({ix} < 0) {{
+                    {ix} = -{ix};
+                }}
+                {ix} = 1 + ({ix} - 1) % (({xsize} - 1) * 2);
+                {ix} = {min_func}({ix}, 2 * {xsize} - 2 - {ix});
+            }}'''
     elif mode == 'nearest':
-        ops = '''
-        {ix} = {min}({max}(({T}){ix}, ({T})0), ({T})({xsize} - 1));'''.format(
-            ix=ix, xsize=xsize, min=min_func, max=max_func,
-            # force using 64-bit signed integer for ptrdiff_t,
-            # see cupy/cupy#6048
-            T=('int' if int_t == 'int' else 'long long'))
+        T = 'int' if int_t == 'int' else 'long long'
+        if separate:
+            ops_lower = f'''{ix} = {max_func}(({T}){ix}, ({T})0);'''
+            ops_upper = f'''{ix} = {min_func}(({T}){ix}, ({T})({xsize} - 1));'''  # noqa
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''{ix} = {min_func}({max_func}(({T}){ix}, ({T})0), ({T})({xsize} - 1));'''  # noqa
     elif mode == 'grid-wrap':
-        ops = '''
-        {ix} %= {xsize};
-        if ({ix} < 0) {{
-            {ix} += {xsize};
-        }}'''.format(ix=ix, xsize=xsize)
+        if separate:
+            ops_upper = f'''
+            {ix} %= {xsize};
+            '''
+            ops_lower = ops_upper + f'''
+            if ({ix} < 0) {{
+                {ix} += {xsize};
+            }}'''
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''
+            {ix} %= {xsize};
+            if ({ix} < 0) {{
+                {ix} += {xsize};
+            }}'''
+
     elif mode == 'wrap':
-        ops = '''
-        if ({ix} < 0) {{
-            {ix} += ({sz} - 1) * (({int_t})(-{ix} / ({sz} - 1)) + 1);
-        }} else if ({ix} > ({sz} - 1)) {{
-            {ix} -= ({sz} - 1) * ({int_t})({ix} / ({sz} - 1));
-        }};'''.format(ix=ix, sz=xsize, int_t=int_t)
+        if separate:
+            ops_lower = f'''{ix} += ({xsize} - 1) * (({int_t})(-{ix} / ({xsize} - 1)) + 1);'''  # noqa
+            ops_upper = f'''{ix} -= ({xsize} - 1) * ({int_t})({ix} / ({xsize} - 1));'''  # noqa
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''
+            if ({ix} < 0) {{
+                {ix} += ({xsize} - 1) * (({int_t})(-{ix} / ({xsize} - 1)) + 1);
+            }} else if ({ix} > ({xsize} - 1)) {{
+                {ix} -= ({xsize} - 1) * ({int_t})({ix} / ({xsize} - 1));
+            }};'''
     elif mode in ['constant', 'grid-constant']:
-        ops = '''
-        if (({ix} < 0) || {ix} >= {xsize}) {{
-            {ix} = -1;
-        }}'''.format(ix=ix, xsize=xsize)
+        if separate:
+            ops_lower = f'''
+            if ({ix} < 0) {{
+                {ix} = -1;
+            }}'''
+            ops_upper = f'''
+            if ({ix} >= {xsize}) {{
+                {ix} = -1;
+            }}'''
+            ops = (ops_lower, ops_upper)
+        else:
+            ops = f'''
+            if (({ix} < 0) || {ix} >= {xsize}) {{
+                {ix} = -1;
+            }}'''
+        if separate:
+            ops = (ops, ops)
     return ops
 
 
diff --git a/python/cucim/src/cucim/skimage/_vendored/_texture.py b/python/cucim/src/cucim/skimage/_vendored/_texture.py
new file mode 100644
index 000000000..5c18bce71
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_texture.py
@@ -0,0 +1,194 @@
+import cupy
+from cupy import _core
+from cupy.cuda import runtime, texture
+
+_affine_transform_2d_array_kernel = _core.ElementwiseKernel(
+    'U texObj, raw float32 m, uint64 width', 'T transformed_image',
+    '''
+    float3 pixel = make_float3(
+        (float)(i / width),
+        (float)(i % width),
+        1.0f
+    );
+    float x = dot(pixel, make_float3(m[0],  m[1],  m[2])) + .5f;
+    float y = dot(pixel, make_float3(m[3],  m[4],  m[5])) + .5f;
+    transformed_image = tex2D<T>(texObj, y, x);
+    ''',
+    'cupyx_texture_affine_transformation_2d_array',
+    preamble='''
+    inline __host__ __device__ float dot(float3 a, float3 b)
+    {
+        return a.x * b.x + a.y * b.y + a.z * b.z;
+    }
+    ''')
+
+
+_affine_transform_3d_array_kernel = _core.ElementwiseKernel(
+    'U texObj, raw float32 m, uint64 height, uint64 width',
+    'T transformed_volume',
+    '''
+    float4 voxel = make_float4(
+        (float)(i / (width * height)),
+        (float)((i % (width * height)) / width),
+        (float)((i % (width * height)) % width),
+        1.0f
+    );
+    float x = dot(voxel, make_float4(m[0],  m[1],  m[2],  m[3])) + .5f;
+    float y = dot(voxel, make_float4(m[4],  m[5],  m[6],  m[7])) + .5f;
+    float z = dot(voxel, make_float4(m[8],  m[9],  m[10], m[11])) + .5f;
+    transformed_volume = tex3D<T>(texObj, z, y, x);
+    ''',
+    'cupyx_texture_affine_transformation_3d_array',
+    preamble='''
+    inline __host__ __device__ float dot(float4 a, float4 b)
+    {
+        return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
+    }
+    ''')
+
+
+def _create_texture_object(data,
+                           address_mode: str,
+                           filter_mode: str,
+                           read_mode: str,
+                           border_color=0):
+
+    if cupy.issubdtype(data.dtype, cupy.unsignedinteger):
+        fmt_kind = runtime.cudaChannelFormatKindUnsigned
+    elif cupy.issubdtype(data.dtype, cupy.integer):
+        fmt_kind = runtime.cudaChannelFormatKindSigned
+    elif cupy.issubdtype(data.dtype, cupy.floating):
+        fmt_kind = runtime.cudaChannelFormatKindFloat
+    else:
+        raise ValueError(f'Unsupported data type {data.dtype}')
+
+    if address_mode == 'nearest':
+        address_mode = runtime.cudaAddressModeClamp
+    elif address_mode == 'constant':
+        address_mode = runtime.cudaAddressModeBorder
+    else:
+        raise ValueError(
+            f'Unsupported address mode {address_mode} '
+            '(supported: constant, nearest)')
+
+    if filter_mode == 'nearest':
+        filter_mode = runtime.cudaFilterModePoint
+    elif filter_mode == 'linear':
+        filter_mode = runtime.cudaFilterModeLinear
+    else:
+        raise ValueError(
+            f'Unsupported filter mode {filter_mode} '
+            f'(supported: nearest, linear)')
+
+    if read_mode == 'element_type':
+        read_mode = runtime.cudaReadModeElementType
+    elif read_mode == 'normalized_float':
+        read_mode = runtime.cudaReadModeNormalizedFloat
+    else:
+        raise ValueError(
+            f'Unsupported read mode {read_mode} '
+            '(supported: element_type, normalized_float)')
+
+    texture_fmt = texture.ChannelFormatDescriptor(
+        data.itemsize * 8, 0, 0, 0, fmt_kind)
+    # CUDAArray: last dimension is the fastest changing dimension
+    array = texture.CUDAarray(texture_fmt, *data.shape[::-1])
+    res_desc = texture.ResourceDescriptor(
+        runtime.cudaResourceTypeArray, cuArr=array)
+    # TODO(the-lay): each dimension can have a different addressing mode
+    # TODO(the-lay): border color/value can be defined for up to 4 channels
+    tex_desc = texture.TextureDescriptor(
+        (address_mode, ) * data.ndim, filter_mode, read_mode,
+        borderColors=(border_color, ))
+    tex_obj = texture.TextureObject(res_desc, tex_desc)
+    array.copy_from(data)
+
+    return tex_obj
+
+
+def affine_transformation(data,
+                          transformation_matrix,
+                          output_shape=None,
+                          output=None,
+                          interpolation: str = 'linear',
+                          mode: str = 'constant',
+                          border_value=0):
+    """
+    Apply an affine transformation.
+
+    The method uses texture memory and supports only 2D and 3D float32 arrays
+    without channel dimension.
+
+    Args:
+        data (cupy.ndarray): The input array or texture object.
+        transformation_matrix (cupy.ndarray): Affine transformation matrix.
+            Must be a homogeneous and have shape ``(ndim + 1, ndim + 1)``.
+        output_shape (tuple of ints): Shape of output. If not specified,
+            the input array shape is used. Default is None.
+        output (cupy.ndarray or ~cupy.dtype): The array in which to place the
+            output, or the dtype of the returned array. If not specified,
+            creates the output array with shape of ``output_shape``. Default is
+            None.
+        interpolation (str): Specifies interpolation mode: ``'linear'`` or
+            ``'nearest'``. Default is ``'linear'``.
+        mode (str): Specifies addressing mode for points outside of the array:
+            (`'constant'``, ``'nearest'``). Default is ``'constant'``.
+        border_value: Specifies value to be used for coordinates outside
+            of the array for ``'constant'`` mode. Default is 0.
+
+    Returns:
+        cupy.ndarray:
+            The transformed input.
+
+    .. seealso:: :func:`cupyx.scipy.ndimage.affine_transform`
+    """
+
+    ndim = data.ndim
+    if (ndim < 2) or (ndim > 3):
+        raise ValueError(
+            'Texture memory affine transformation is defined only for '
+            '2D and 3D arrays without channel dimension.')
+
+    dtype = data.dtype
+    if dtype != cupy.float32:
+        raise ValueError(f'Texture memory affine transformation is available '
+                         f'only for float32 data type (not {dtype})')
+
+    if interpolation not in ['linear', 'nearest']:
+        raise ValueError(
+            f'Unsupported interpolation {interpolation} '
+            f'(supported: linear, nearest)')
+
+    if transformation_matrix.shape != (ndim + 1, ndim + 1):
+        raise ValueError('Matrix must be have shape (ndim + 1, ndim + 1)')
+
+    texture_object = _create_texture_object(data,
+                                            address_mode=mode,
+                                            filter_mode=interpolation,
+                                            read_mode='element_type',
+                                            border_color=border_value)
+
+    if ndim == 2:
+        kernel = _affine_transform_2d_array_kernel
+    else:
+        kernel = _affine_transform_3d_array_kernel
+
+    if output_shape is None:
+        output_shape = data.shape
+
+    if output is None:
+        output = cupy.zeros(output_shape, dtype=dtype)
+    elif isinstance(output, (type, cupy.dtype)):
+        if output != cupy.float32:
+            raise ValueError(f'Texture memory affine transformation is '
+                             f'available only for float32 data type (not '
+                             f'{output})')
+        output = cupy.zeros(output_shape, dtype=output)
+    elif isinstance(output, cupy.ndarray):
+        if output.shape != output_shape:
+            raise ValueError('Output shapes do not match')
+    else:
+        raise ValueError('Output must be None, cupy.ndarray or cupy.dtype')
+
+    kernel(texture_object, transformation_matrix, *output_shape[1:], output)
+    return output
diff --git a/python/cucim/src/cucim/skimage/_vendored/ndimage.py b/python/cucim/src/cucim/skimage/_vendored/ndimage.py
new file mode 100644
index 000000000..cd97578a8
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/ndimage.py
@@ -0,0 +1,82 @@
+# locally defined filters that are more efficient than in CuPy
+from cucim.skimage._vendored._ndimage_filters import correlate  # NOQA
+from cucim.skimage._vendored._ndimage_filters import convolve  # NOQA
+from cucim.skimage._vendored._ndimage_filters import correlate1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import convolve1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import uniform_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import uniform_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import gaussian_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import gaussian_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import prewitt  # NOQA
+from cucim.skimage._vendored._ndimage_filters import sobel  # NOQA
+from cucim.skimage._vendored._ndimage_filters import generic_laplace  # NOQA
+from cucim.skimage._vendored._ndimage_filters import laplace  # NOQA
+from cucim.skimage._vendored._ndimage_filters import gaussian_laplace  # NOQA
+from cucim.skimage._vendored._ndimage_filters import generic_gradient_magnitude  # NOQA
+from cucim.skimage._vendored._ndimage_filters import gaussian_gradient_magnitude  # NOQA
+from cucim.skimage._vendored._ndimage_filters import minimum_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import maximum_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import minimum_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import maximum_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import median_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import rank_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import percentile_filter  # NOQA
+
+# interpolation
+from cucim.skimage._vendored._ndimage_interpolation import affine_transform  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import map_coordinates  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import rotate  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import shift  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import spline_filter  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import spline_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_interpolation import zoom  # NOQA
+
+# morphology
+from cucim.skimage._vendored._ndimage_morphology import generate_binary_structure  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import iterate_structure  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_erosion  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_dilation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_opening  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_closing  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_hit_or_miss  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_fill_holes  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_propagation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_erosion  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_dilation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_closing  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_opening  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import morphological_gradient  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import morphological_laplace  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import white_tophat  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import black_tophat  # NOQA
+
+# Import the rest of the cupyx.scipy.ndimage API here
+
+# additional filters
+from cupyx.scipy.ndimage import generic_filter  # NOQA
+from cupyx.scipy.ndimage import generic_filter1d  # NOQA
+
+# fourier filters
+from cupyx.scipy.ndimage import fourier_ellipsoid  # NOQA
+from cupyx.scipy.ndimage import fourier_gaussian  # NOQA
+from cupyx.scipy.ndimage import fourier_shift  # NOQA
+from cupyx.scipy.ndimage import fourier_uniform  # NOQA
+
+# measurements
+from cupyx.scipy.ndimage import label  # NOQA
+try:
+    from cupyx.scipy.ndimage import sum_labels  # NOQA
+except ImportError:
+    from cupyx.scipy.ndimage import sum as sum_labels  # NOQA
+from cupyx.scipy.ndimage import mean  # NOQA
+from cupyx.scipy.ndimage import variance  # NOQA
+from cupyx.scipy.ndimage import standard_deviation  # NOQA
+from cupyx.scipy.ndimage import minimum  # NOQA
+from cupyx.scipy.ndimage import maximum  # NOQA
+from cupyx.scipy.ndimage import minimum_position  # NOQA
+from cupyx.scipy.ndimage import maximum_position  # NOQA
+from cupyx.scipy.ndimage import median  # NOQA
+from cupyx.scipy.ndimage import extrema  # NOQA
+from cupyx.scipy.ndimage import center_of_mass  # NOQA
+from cupyx.scipy.ndimage import histogram  # NOQA
+from cupyx.scipy.ndimage import labeled_comprehension  # NOQA
diff --git a/python/cucim/src/cucim/skimage/feature/_canny.py b/python/cucim/src/cucim/skimage/feature/_canny.py
index 99d85d465..56d0f5cab 100644
--- a/python/cucim/src/cucim/skimage/feature/_canny.py
+++ b/python/cucim/src/cucim/skimage/feature/_canny.py
@@ -12,7 +12,7 @@
 Original author: Lee Kamentsky
 """
 import cupy as cp
-import cupyx.scipy.ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from cucim.skimage.util import dtype_limits
 
diff --git a/python/cucim/src/cucim/skimage/feature/corner.py b/python/cucim/src/cucim/skimage/feature/corner.py
index 7a746eed4..f52f98ad7 100644
--- a/python/cucim/src/cucim/skimage/feature/corner.py
+++ b/python/cucim/src/cucim/skimage/feature/corner.py
@@ -3,7 +3,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 from scipy import spatial  # TODO: use RAPIDS cuSpatial?
 
 from cucim.skimage.util import img_as_float
diff --git a/python/cucim/src/cucim/skimage/feature/peak.py b/python/cucim/src/cucim/skimage/feature/peak.py
index eff5f6735..fbb2570f7 100644
--- a/python/cucim/src/cucim/skimage/feature/peak.py
+++ b/python/cucim/src/cucim/skimage/feature/peak.py
@@ -1,7 +1,7 @@
 from warnings import warn
 
 import cupy as cp
-import cupyx.scipy.ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 import numpy as np
 from scipy.ndimage import find_objects as cpu_find_objects
 
diff --git a/python/cucim/src/cucim/skimage/filters/_gabor.py b/python/cucim/src/cucim/skimage/filters/_gabor.py
index e55720357..df4c041c4 100644
--- a/python/cucim/src/cucim/skimage/filters/_gabor.py
+++ b/python/cucim/src/cucim/skimage/filters/_gabor.py
@@ -1,7 +1,8 @@
 import math
 
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
+
 
 from .._shared.utils import _supported_float_type, check_nD, warn
 
diff --git a/python/cucim/src/cucim/skimage/filters/_median.py b/python/cucim/src/cucim/skimage/filters/_median.py
index 821cb0149..de3f898c4 100644
--- a/python/cucim/src/cucim/skimage/filters/_median.py
+++ b/python/cucim/src/cucim/skimage/filters/_median.py
@@ -1,8 +1,8 @@
 from warnings import warn
 
 import numpy as np
-from cupyx.scipy import ndimage as ndi
 
+import cucim.skimage._vendored.ndimage as ndi
 from .._shared.utils import deprecate_kwarg
 
 
diff --git a/python/cucim/src/cucim/skimage/filters/_separable_filtering.py b/python/cucim/src/cucim/skimage/filters/_separable_filtering.py
new file mode 100644
index 000000000..eb8dad28a
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/filters/_separable_filtering.py
@@ -0,0 +1,990 @@
+import math
+
+import cupy as cp
+
+from cucim.skimage._vendored import _ndimage_util as util
+from cucim.skimage._vendored._internal import _normalize_axis_index, prod
+from cucim.skimage._vendored._ndimage_filters_core import (
+    _ndimage_CAST_FUNCTION, _ndimage_includes)
+
+
+def _get_constants(ndim, axis, kernel_size, anchor, patch_per_block=None):
+    if anchor is None:
+        anchor = kernel_size // 2
+    halo_pixels_needed = max(kernel_size - anchor, anchor)
+    if patch_per_block is None:
+        patch_per_block = 4
+
+    if ndim == 2:
+        # note, in 2d axis 0 = "y"
+        #             axis 1 = "x"
+        # for simplicity, keeping same halo size at both start and end
+        if axis == 1:
+            # as in OpenCV's column_filter.hpp
+            block_x = 16
+            block_y = 16
+            halo_size = math.ceil(halo_pixels_needed / block_x)
+        elif axis == 0:
+            # as in OpenCV's row_filter.hpp
+            block_x = 32  # 16 in CUDA example
+            block_y = 8  # 4 in CUDA example
+            halo_size = math.ceil(halo_pixels_needed / block_y)
+        # can have out of bounds access unless patch_per_block >= halo_size
+        patch_per_block = max(patch_per_block, halo_size)
+        block_z = 1
+    elif ndim == 3:
+        # note, in 3d axis 0 = "z"
+        #             axis 1 = "y"
+        #             axis 2 = "x"
+        # for simplicity, keeping same halo size at both start and end
+        if axis == 2:
+            # as in OpenCV's column_filter.hpp
+            block_x = 16
+            block_y = 4
+            block_z = 4
+            halo_size = math.ceil(halo_pixels_needed / block_x)
+        elif axis == 1:
+            # as in OpenCV's column_filter.hpp
+            block_x = 32
+            block_y = 4
+            block_z = 4
+            halo_size = math.ceil(halo_pixels_needed / block_y)
+        elif axis == 0:
+            # as in OpenCV's row_filter.hpp
+            block_x = 32
+            block_y = 4
+            block_z = 4
+            halo_size = math.ceil(halo_pixels_needed / block_z)
+        # can have out of bounds access unless patch_per_block >= halo_size
+        patch_per_block = max(patch_per_block, halo_size)
+    else:
+        raise NotImplementedError("Only 2D and 3D are currently supported")
+    block = (block_x, block_y, block_z)
+    return block, patch_per_block, halo_size
+
+
+def _get_smem_shape(ndim, axis, block, patch_per_block, halo_size, anchor=None,
+                    image_dtype=cp.float32):
+    bx, by, bz = block
+    if ndim == 2:
+        if axis == 0:
+            shape = ((patch_per_block + 2 * halo_size) * by, bx)
+        elif axis == 1:
+            shape = (by, (patch_per_block + 2 * halo_size) * bx)
+    elif ndim == 3:
+        if axis == 0:
+            shape = ((patch_per_block + 2 * halo_size) * bz, by, bx)
+        elif axis == 1:
+            shape = (bz, (patch_per_block + 2 * halo_size) * by, bx)
+        elif axis == 2:
+            shape = (bz, by, (patch_per_block + 2 * halo_size) * bx)
+    else:
+        raise NotImplementedError("TODO")
+    nbytes = cp.dtype(image_dtype).itemsize * prod(shape)
+    return shape, nbytes
+
+
+def _get_warp_size(device_id=None):
+    if device_id is None:
+        device_id = cp.cuda.runtime.getDevice()
+    device_props = cp.cuda.runtime.getDeviceProperties(device_id)
+    return device_props['warpSize']
+
+
+def _get_shmem_limits(device_id=None):
+    if device_id is None:
+        device_id = cp.cuda.runtime.getDevice()
+    device_props = cp.cuda.runtime.getDeviceProperties(device_id)
+    shared_mp = device_props.get('sharedMemPerMultiprocessor', None)
+    shared_block = device_props.get('sharedMemPerBlock', None)
+    shared_block_optin = device_props.get('sharedMemPerBlockOptin', None)
+    global_l1_cache_supported = device_props.get('globalL1CacheSupported', None)
+    local_l1_cache_supported = device_props.get('localL1CacheSupported', None)
+    l2_size = device_props.get('l2CacheSize', None)
+    warp_size = device_props.get('warpSize', None)
+    regs_per_block = device_props.get('regsPerBlock', None)
+    return {
+        'device_id': device_id,
+        'shared_mp': shared_mp,
+        'shared_block': shared_block,
+        'shared_block_optin': shared_block_optin,
+        'global_l1_cache_supported': global_l1_cache_supported,
+        'local_l1_cache_supported': local_l1_cache_supported,
+        'l2_size': l2_size,
+        'warp_size': warp_size,
+        'regs_per_block': regs_per_block,
+    }
+
+
+class ResourceLimitError(RuntimeError):
+    pass
+
+
+@cp.memoize(for_each_device=True)
+def _check_smem_availability(ndim, axis, kernel_size, anchor=None,
+                             patch_per_block=None, image_dtype=cp.float32,
+                             device_id=None):
+    block, patch_per_block, halo_size = _get_constants(
+        ndim, axis, kernel_size, anchor=anchor, patch_per_block=patch_per_block
+    )
+    shape, nbytes = _get_smem_shape(
+        ndim, axis, block, patch_per_block, halo_size, image_dtype
+    )
+    props = _get_shmem_limits(device_id=device_id)
+    if nbytes > props['shared_block']:
+        raise ResourceLimitError("inadequate shared memory available")
+
+
+_dtype_char_to_c_types = {
+    'e': 'float16',
+    'f': 'float',
+    'd': 'double',
+    'F': 'complex<float>',
+    'D': 'complex<double>',
+    '?': 'char',
+    'b': 'char',
+    'h': 'short',
+    'i': 'int',
+    'l': 'long long',
+    'B': 'unsigned char',
+    'H': 'unsigned short',
+    'I': 'unsigned int',
+    'L': 'unsigned long long',
+}
+
+
+def _get_code_stage1_shared_memory_load_2d(ndim, axis, mode, cval):
+    """Generates the first stage of the function body.
+
+    This involves just copying from the `src` array into the `smem` shared
+    memory array followed by a call to __syncthreads(). All boundary
+    handling also occurs within this function.
+    """
+
+    if ndim == 2 and axis == 0:
+        if mode not in ['constant', 'grid-constant']:
+            boundary_code_lower, boundary_code_upper = util._generate_boundary_condition_ops(mode, 'row', 'n_rows', separate=True)  # noqa
+
+        # as in OpenCV's column_filter.hpp
+        code = """
+        __shared__ T smem[(PATCH_PER_BLOCK + 2 * HALO_SIZE) * BLOCK_DIM_Y][BLOCK_DIM_X];
+        const int x = blockIdx.x * BLOCK_DIM_X + threadIdx.x;
+        if (x >= n_cols){
+            return;
+        }
+        const T* src_col = &src[x];
+        const int yStart = blockIdx.y * (BLOCK_DIM_Y * PATCH_PER_BLOCK) + threadIdx.y;
+
+        // memory is contiguous along last (columns) axis
+        const int row_stride = n_cols;  // stride (in elements) along axis 0
+        int row;
+
+        if (blockIdx.y > 0)
+        {
+            //Upper halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart - (HALO_SIZE - j) * BLOCK_DIM_Y) * row_stride];
+        }
+        else
+        {
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j) {
+                row = yStart - (HALO_SIZE - j) * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row < 0)
+                    smem[threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_lower
+        code += """
+                    smem[threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = src_col[row * row_stride];
+            }
+        }
+
+        if (blockIdx.y + 2 < gridDim.y)  // Note: +2 here assumes HALO_SIZE <= PATCH_PER_BLOCK so we ensure that elsewhere
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+                smem[threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart + j * BLOCK_DIM_Y) * row_stride];
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Y) * row_stride];
+        }
+        else
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j) {
+                row = yStart + j * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= n_rows)
+                    smem[threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[row * row_stride];
+            }
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+            {
+                row = yStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= n_rows)
+                    smem[threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[row * row_stride];
+            }
+        }
+        """  # noqa
+    elif ndim == 2 and axis == 1:
+        if mode not in ['constant', 'grid-constant']:
+            boundary_code_lower, boundary_code_upper = util._generate_boundary_condition_ops(mode, 'col', 'n_cols', separate=True)  # noqa
+
+        # as in OpenCV's row_filter.hpp
+        code = """
+        __shared__ T smem[BLOCK_DIM_Y][(PATCH_PER_BLOCK + 2 * HALO_SIZE) * BLOCK_DIM_X];
+        const int y = blockIdx.y * BLOCK_DIM_Y + threadIdx.y;
+        if (y >= n_rows) {
+            return;
+        }
+        const int row_stride = n_cols;  // stride (in elements) along axis 0
+        int col;
+        const T* src_row = &src[y * row_stride];
+        const int xStart = blockIdx.x * (PATCH_PER_BLOCK * BLOCK_DIM_X) + threadIdx.x;
+
+        if (blockIdx.x > 0)
+        {
+            //Load left halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = src_row[xStart - (HALO_SIZE - j) * BLOCK_DIM_X];
+        }
+        else
+        {
+            //Load left halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j){
+                col = xStart - (HALO_SIZE - j) * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col < 0)
+                    smem[threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_lower
+        code += """
+                    smem[threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = src_row[col];
+            }
+        }
+        if (blockIdx.x + 2 < gridDim.x)  // Note: +2 here assumes HALO_SIZE <= PATCH_PER_BLOCK so we ensure that elsewhere
+        {
+            //Load main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+                smem[threadIdx.y][threadIdx.x + (HALO_SIZE + j)* BLOCK_DIM_X] = src_row[xStart + j * BLOCK_DIM_X];
+
+            //Load right halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = src_row[xStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_X];
+        }
+        else
+        {
+            //Load main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j) {
+                col = xStart + j * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col >= n_cols)
+                    smem[threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X] = src_row[col];
+            }
+
+            //Load right halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j){
+                col = xStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col >= n_cols)
+                    smem[threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = src_row[col];
+            }
+        }
+        """  # noqa
+
+    code += """
+        __syncthreads();
+    """
+    return code
+
+
+def _get_code_stage1_shared_memory_load_3d(ndim, axis, mode, cval):
+    """Generates the first stage of the function body.
+
+    This involves just copying from the `src` array into the `smem` shared
+    memory array followed by a call to __syncthreads(). All boundary
+    handling also occurs within this function.
+    """
+
+    if ndim == 3 and axis == 0:
+        if mode not in ['constant', 'grid-constant']:
+            boundary_code_lower, boundary_code_upper = util._generate_boundary_condition_ops(mode, 'row', 's_0', separate=True)  # noqa
+
+        # as in OpenCV's column_filter.hpp
+        code = """
+        __shared__ T smem[(PATCH_PER_BLOCK + 2 * HALO_SIZE) * BLOCK_DIM_Z][BLOCK_DIM_Y][BLOCK_DIM_X];
+        const int x = blockIdx.x * BLOCK_DIM_X + threadIdx.x;
+        const int y = blockIdx.y * BLOCK_DIM_Y + threadIdx.y;
+        if ((x >= s_2) || (y >= s_1)) {
+            return;
+        }
+        // memory is contiguous along last (columns) axis
+        const int stride_0 = s_1 * s_2;  // stride (in elements) along axis 0
+        const int stride_1 = s_2;  // stride (in elements) along axis 1
+
+        const T* src_col = &src[x + stride_1 * y];
+        const int zStart = blockIdx.z * (BLOCK_DIM_Z * PATCH_PER_BLOCK) + threadIdx.z;
+
+        int row;
+
+        if (blockIdx.z > 0)
+        {
+            //Upper halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z + j * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[(zStart - (HALO_SIZE - j) * BLOCK_DIM_Z) * stride_0];
+        }
+        else
+        {
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j) {
+                row = zStart - (HALO_SIZE - j) * BLOCK_DIM_Z;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row < 0)
+                    smem[threadIdx.z + j * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_lower
+        code += """
+                    smem[threadIdx.z + j * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[row * stride_0];
+            }
+        }
+
+        if (blockIdx.z + 2 < gridDim.z)  // Note: +2 here assumes HALO_SIZE <= PATCH_PER_BLOCK so we ensure that elsewhere
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+                smem[threadIdx.z + (HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[(zStart + j * BLOCK_DIM_Z) * stride_0];
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[(zStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Z) * stride_0];
+        }
+        else
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j) {
+                row = zStart + j * BLOCK_DIM_Z;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= s_0)
+                    smem[threadIdx.z + (HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z + (HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[row * stride_0];
+            }
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+            {
+                row = zStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Z;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= s_0)
+                    smem[threadIdx.z + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Z][threadIdx.y][threadIdx.x] = src_col[row * stride_0];
+            }
+        }
+        """  # noqa
+    elif ndim == 3 and axis == 1:
+        if mode not in ['constant', 'grid-constant']:
+            boundary_code_lower, boundary_code_upper = util._generate_boundary_condition_ops(mode, 'row', 's_1', separate=True)  # noqa
+
+        # as in OpenCV's column_filter.hpp
+        code = """
+        __shared__ T smem[BLOCK_DIM_Z][(PATCH_PER_BLOCK + 2 * HALO_SIZE) * BLOCK_DIM_Y][BLOCK_DIM_X];
+        const int x = blockIdx.x * BLOCK_DIM_X + threadIdx.x;
+        const int z = blockIdx.z * BLOCK_DIM_Z + threadIdx.z;
+        if ((x >= s_2) || (z >= s_0)) {
+            return;
+        }
+        // memory is contiguous along last (columns) axis
+        const int stride_0 = s_1 * s_2;  // stride (in elements) along axis 0
+        const int stride_1 = s_2;  // stride (in elements) along axis 1
+
+        const T* src_col = &src[x + stride_0 * z];
+        const int yStart = blockIdx.y * (BLOCK_DIM_Y * PATCH_PER_BLOCK) + threadIdx.y;
+
+        int row;
+
+        if (blockIdx.y > 0)
+        {
+            //Upper halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z][threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart - (HALO_SIZE - j) * BLOCK_DIM_Y) * stride_1];
+        }
+        else
+        {
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j) {
+                row = yStart - (HALO_SIZE - j) * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row < 0)
+                    smem[threadIdx.z][threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_lower
+        code += """
+                    smem[threadIdx.z][threadIdx.y + j * BLOCK_DIM_Y][threadIdx.x] = src_col[row * stride_1];
+            }
+        }
+
+        if (blockIdx.y + 2 < gridDim.y)  // Note: +2 here assumes HALO_SIZE <= PATCH_PER_BLOCK so we ensure that elsewhere
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+                smem[threadIdx.z][threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart + j * BLOCK_DIM_Y) * stride_1];
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z][threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[(yStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Y) * stride_1];
+        }
+        else
+        {
+            //Main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j) {
+                row = yStart + j * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= s_1)
+                    smem[threadIdx.z][threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z][threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[row * stride_1];
+            }
+
+            //Lower halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+            {
+                row = yStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_Y;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (row >= s_1)
+                    smem[threadIdx.z][threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z][threadIdx.y + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_Y][threadIdx.x] = src_col[row * stride_1];
+            }
+        }
+        """  # noqa
+    elif ndim == 3 and axis == 2:
+        if mode not in ['constant', 'grid-constant']:
+            boundary_code_lower, boundary_code_upper = util._generate_boundary_condition_ops(mode, 'col', 's_2', separate=True)  # noqa
+
+        # as in OpenCV's row_filter.hpp
+        code = """
+        __shared__ T smem[BLOCK_DIM_Z][BLOCK_DIM_Y][(PATCH_PER_BLOCK + 2 * HALO_SIZE) * BLOCK_DIM_X];
+        const int y = blockIdx.y * BLOCK_DIM_Y + threadIdx.y;
+        const int z = blockIdx.z * BLOCK_DIM_Z + threadIdx.z;
+        if ((y >= s_1) || (z >= s_0)) {
+            return;
+        }
+        const int stride_0 = s_1 * s_2;  // stride (in elements) along axis 0
+        const int stride_1 = s_2;  // stride (in elements) along axis 1
+        int col;
+        const T* src_row = &src[z * stride_0 + y * stride_1];
+        const int xStart = blockIdx.x * (PATCH_PER_BLOCK * BLOCK_DIM_X) + threadIdx.x;
+
+        if (blockIdx.x > 0)
+        {
+            //Load left halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z][threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = src_row[xStart - (HALO_SIZE - j) * BLOCK_DIM_X];
+        }
+        else
+        {
+            //Load left halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j){
+                col = xStart - (HALO_SIZE - j) * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col < 0)
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_lower
+        code += """
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + j * BLOCK_DIM_X] = src_row[col];
+            }
+        }
+        if (blockIdx.x + 2 < gridDim.x)  // Note: +2 here assumes HALO_SIZE <= PATCH_PER_BLOCK so we ensure that elsewhere
+        {
+            //Load main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+                smem[threadIdx.z][threadIdx.y][threadIdx.x + (HALO_SIZE + j)* BLOCK_DIM_X] = src_row[xStart + j * BLOCK_DIM_X];
+
+            //Load right halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j)
+                smem[threadIdx.z][threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = src_row[xStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_X];
+        }
+        else
+        {
+            //Load main data
+            #pragma unroll
+            for (int j = 0; j < PATCH_PER_BLOCK; ++j) {
+                col = xStart + j * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col >= s_2)
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X] = src_row[col];
+            }
+
+            //Load right halo
+            #pragma unroll
+            for (int j = 0; j < HALO_SIZE; ++j){
+                col = xStart + (PATCH_PER_BLOCK + j) * BLOCK_DIM_X;
+        """  # noqa
+        if mode == 'constant':
+            code += f"""
+                if (col >= s_2)
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = static_cast<T>({cval});
+                else
+            """  # noqa
+        else:
+            code += boundary_code_upper
+        code += """
+                    smem[threadIdx.z][threadIdx.y][threadIdx.x + (PATCH_PER_BLOCK + HALO_SIZE + j) * BLOCK_DIM_X] = src_row[col];
+            }
+        }
+        """  # noqa
+
+    code += """
+        __syncthreads();
+    """
+    return code
+
+
+@cp.memoize(for_each_device=False)
+def _get_code_stage1_shared_memory_load(ndim, axis, mode, cval):
+    if ndim == 2:
+        return _get_code_stage1_shared_memory_load_2d(ndim, axis, mode, cval)
+    elif ndim == 3:
+        return _get_code_stage1_shared_memory_load_3d(ndim, axis, mode, cval)
+
+
+def _get_code_stage2_convolve_2d(ndim, axis, flip_kernel):
+    code = """
+    #pragma unroll
+    for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+    {
+    """
+    if flip_kernel:
+        kernel_idx = "KSIZE - 1 - k"
+    else:
+        kernel_idx = "k"
+
+    if ndim == 2 and axis == 0:
+        code += """
+        const int y = yStart + j * BLOCK_DIM_Y;
+
+        if (y < n_rows)
+        {
+        """
+        inner = f"""
+                sum = sum + static_cast<W>(smem[threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y - anchor + k][threadIdx.x]) * kernel[{kernel_idx}];
+        """  # noqa
+    elif ndim == 2 and axis == 1:
+        code += """
+        const int x = xStart + j * BLOCK_DIM_X;
+
+        if (x < n_cols)
+        {
+        """
+        inner = f"""
+                sum = sum + static_cast<W>(smem[threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X - anchor + k]) * kernel[{kernel_idx}];
+        """  # noqa
+    code += f"""
+            W sum = static_cast<W>(0);
+
+            #pragma unroll
+            for (int k = 0; k < KSIZE; ++k) {{
+               {inner}
+            }}
+            dst[y * row_stride + x] = cast<D>(sum);
+        }}
+    }}
+    """
+    return code
+
+
+def _get_code_stage2_convolve_3d(ndim, axis, flip_kernel):
+    code = """
+    #pragma unroll
+    for (int j = 0; j < PATCH_PER_BLOCK; ++j)
+    {
+    """
+    if flip_kernel:
+        kernel_idx = "KSIZE - 1 - k"
+    else:
+        kernel_idx = "k"
+
+    if ndim == 3 and axis == 0:
+        code += """
+        const int z = zStart + j * BLOCK_DIM_Z;
+
+        if (z < s_0)
+        {
+        """
+        inner = f"""
+                sum = sum + static_cast<W>(smem[threadIdx.z + (HALO_SIZE + j) * BLOCK_DIM_Z - anchor + k][threadIdx.y][threadIdx.x]) * kernel[{kernel_idx}];
+        """  # noqa
+    elif ndim == 3 and axis == 1:
+        code += """
+        const int y = yStart + j * BLOCK_DIM_Y;
+
+        if (y < s_1)
+        {
+        """
+        inner = f"""
+                sum = sum + static_cast<W>(smem[threadIdx.z][threadIdx.y + (HALO_SIZE + j) * BLOCK_DIM_Y - anchor + k][threadIdx.x]) * kernel[{kernel_idx}];
+        """  # noqa
+    elif ndim == 3 and axis == 2:
+        code += """
+        const int x = xStart + j * BLOCK_DIM_X;
+
+        if (x < s_2)
+        {
+        """
+        inner = f"""
+                sum = sum + static_cast<W>(smem[threadIdx.z][threadIdx.y][threadIdx.x + (HALO_SIZE + j) * BLOCK_DIM_X - anchor + k]) * kernel[{kernel_idx}];
+        """  # noqa
+    code += f"""
+            W sum = static_cast<W>(0);
+
+            #pragma unroll
+            for (int k = 0; k < KSIZE; ++k) {{
+               {inner}
+            }}
+            dst[z * stride_0 + y * stride_1 + x] = cast<D>(sum);
+        }}
+    }}
+    """
+    return code
+
+
+@cp.memoize(for_each_device=False)
+def _get_code_stage2_convolve(ndim, axis, flip_kernel):
+    if ndim == 2:
+        return _get_code_stage2_convolve_2d(ndim, axis, flip_kernel)
+    elif ndim == 3:
+        return _get_code_stage2_convolve_3d(ndim, axis, flip_kernel)
+
+
+@cp.memoize(for_each_device=True)
+def _get_separable_conv_kernel_src(
+    kernel_size, axis, ndim, anchor, image_c_type, kernel_c_type,
+    output_c_type, mode, cval, patch_per_block=None, flip_kernel=False
+):
+    blocks, patch_per_block, halo_size = _get_constants(
+        ndim, axis, kernel_size, anchor, patch_per_block
+    )
+    block_x, block_y, block_z = blocks
+
+    mode_str = mode
+    if 'constant' in mode_str:
+        mode_str += f'_{cval:0.2f}'.replace('.', '_')
+    mode_str = mode_str.replace('-', '_')
+    if flip_kernel:
+        func_name = f'convolve_s{kernel_size}_{ndim}d_ax{axis}_{mode_str}'
+    else:
+        func_name = f'correlate_s{kernel_size}_{ndim}d_ax{axis}_{mode_str}'
+    func_name += f"_T{image_c_type}_W{kernel_c_type}_D{output_c_type}".replace('complex<', 'c').replace('>', '').replace('long ', 'l').replace('unsigned ', 'u')  # noqa
+    func_name += f"_patch{patch_per_block}_halo{halo_size}"
+    # func_name += f"_bx{block_x}_by{block_y}" // these are fixed per axis
+
+    code = """
+    #include "cupy/carray.cuh"  // for float16
+    #include "cupy/complex.cuh"  // for complex<float>
+    """
+
+    # SciPy-style float -> unsigned integer casting for the output
+    # (use cast<D>(sum) instead of static_cast<D>(sum) for the output)
+    code += _ndimage_includes + _ndimage_CAST_FUNCTION
+
+    code += f"""
+    const int KSIZE = {kernel_size};
+    const int BLOCK_DIM_X = {block_x};
+    const int BLOCK_DIM_Y = {block_y};
+    const int BLOCK_DIM_Z = {block_z};
+    const int PATCH_PER_BLOCK = {patch_per_block};
+    const int HALO_SIZE = {halo_size};
+    typedef {image_c_type}  T;
+    typedef {output_c_type} D;
+    typedef {kernel_c_type} W;
+    """
+
+    if ndim == 2:
+        code += f"""
+        extern "C"{{
+        __global__ void {func_name}(const T *src, D *dst, const W* kernel, const int anchor, int n_rows, int n_cols)
+        {{
+        """  # noqa
+    elif ndim == 3:
+        code += f"""
+        extern "C"{{
+        __global__ void {func_name}(const T *src, D *dst, const W* kernel, const int anchor, int s_0, int s_1, int s_2)
+        {{
+        """  # noqa
+    code += _get_code_stage1_shared_memory_load(ndim, axis, mode, cval)
+    code += _get_code_stage2_convolve(ndim, axis, flip_kernel)
+    code += """
+    }  // end of function
+    }  // extern "C"
+    """
+    return func_name, blocks, patch_per_block, code
+
+
+@cp.memoize(for_each_device=True)
+def _get_separable_conv_kernel(kernel_size, axis, ndim, image_c_type,
+                               kernel_c_type, output_c_type, anchor=None,
+                               mode='nearest', cval=0,
+                               patch_per_block=None, flip_kernel=False):
+    func_name, block, patch_per_block, code = _get_separable_conv_kernel_src(
+        kernel_size=kernel_size,
+        axis=axis,
+        ndim=ndim,
+        image_c_type=image_c_type,
+        kernel_c_type=kernel_c_type,
+        output_c_type=output_c_type,
+        anchor=anchor,
+        mode=mode,
+        cval=cval,
+        patch_per_block=patch_per_block,
+        flip_kernel=flip_kernel,
+    )
+    options = ('--std=c++11', '-DCUPY_USE_JITIFY')
+    m = cp.RawModule(code=code, options=options)
+    return m.get_function(func_name), block, patch_per_block
+
+
+def _get_grid(shape, block, axis, patch_per_block):
+    """Determine grid size from image shape and block parameters"""
+    ndim = len(shape)
+    if ndim == 2:
+        if axis == 0:
+            # column filter
+            grid = (
+                math.ceil(shape[1] / block[0]),
+                math.ceil(shape[0] / (block[1] * patch_per_block)),
+                1,
+            )
+        elif axis == 1:
+            # row filter
+            grid = (
+                math.ceil(shape[1] / (block[0] * patch_per_block)),
+                math.ceil(shape[0] / block[1]),
+                1,
+            )
+        else:
+            raise ValueError(f"invalid axis: {axis}")
+    elif ndim == 3:
+        if axis == 0:
+            # column filter
+            grid = (
+                math.ceil(shape[2] / block[0]),
+                math.ceil(shape[1] / block[1]),
+                math.ceil(shape[0] / (block[2] * patch_per_block)),
+            )
+        elif axis == 1:
+            # row filter
+            grid = (
+                math.ceil(shape[2] / block[0]),
+                math.ceil(shape[1] / (block[1] * patch_per_block)),
+                math.ceil(shape[0] / block[2]),
+            )
+        elif axis == 2:
+            # row filter
+            grid = (
+                math.ceil(shape[2] / (block[0] * patch_per_block)),
+                math.ceil(shape[1] / block[1]),
+                math.ceil(shape[0] / block[2]),
+            )
+        else:
+            raise ValueError(f"invalid axis: {axis}")
+    else:
+        raise NotImplementedError(f"unsupported ndim: {ndim}")
+    return grid
+
+
+def _shmem_convolve1d(image, weights, axis=-1, output=None, mode="reflect",
+                      cval=0.0, origin=0, convolution=False):
+
+    ndim = image.ndim
+    if weights.ndim != 1:
+        raise ValueError("expected 1d weight array")
+    axis = _normalize_axis_index(axis, ndim)
+    origin = util._check_origin(origin, weights.size)
+    if weights.size == 0:
+        return cp.zeros_like(input)
+    util._check_mode(mode)
+
+    if convolution:
+        # use flip_kernel to avoid cp.ascontiguousarray(weights[::-1]))
+        origin = -origin
+        if weights.size % 2 == 0:
+            origin -= 1
+    elif weights.dtype.kind == "c":
+        # numpy.correlate conjugates weights rather than input.
+        weights = weights.conj()
+
+    anchor = weights.size // 2 + origin
+
+    if weights.size > 32:
+        # For large kernels, make sure we have adequate shared memory
+        _check_smem_availability(ndim, axis, weights.size, anchor=anchor,
+                                 patch_per_block=None, image_dtype=image.dtype,
+                                 device_id=None)
+
+    # CUDA kernels assume C-contiguous memory layout
+    if not image.flags.c_contiguous:
+        image = cp.ascontiguousarray(image)
+
+    complex_output = image.dtype.kind == 'c'
+    # Note: important to set use_cucim_casting=True for performance with
+    #       8 and 16-bit integer types. This causes the weights to get cast to
+    #       float32 rather than float64.
+    weights_dtype = util._get_weights_dtype(
+        image, weights, use_cucim_casting=True
+    )
+    if not weights.flags.c_contiguous or weights.dtype != weights_dtype:
+        weights = cp.ascontiguousarray(weights, weights_dtype)
+
+    # promote output to nearest complex dtype if necessary
+    complex_output = complex_output or weights.dtype.kind == 'c'
+    output = util._get_output(output, image, None, complex_output)
+
+    # handle potential overlap between input and output arrays
+    needs_temp = cp.shares_memory(output, image, 'MAY_SHARE_BOUNDS')
+    if needs_temp:
+        output, temp = util._get_output(output.dtype, input), output
+
+    # index_c_type = util._get_inttype(image)
+    image_c_type = _dtype_char_to_c_types[image.dtype.char]
+    weights_c_type = _dtype_char_to_c_types[weights.dtype.char]
+    output_c_type = _dtype_char_to_c_types[output.dtype.char]
+
+    conv_axis_kernel, block, patch_per_block = _get_separable_conv_kernel(
+        weights.size,
+        axis=axis,
+        ndim=ndim,
+        anchor=anchor,
+        image_c_type=image_c_type,
+        kernel_c_type=weights_c_type,
+        output_c_type=output_c_type,
+        mode=mode,
+        cval=cval,
+        patch_per_block=None,
+        flip_kernel=convolution,
+    )
+    grid = _get_grid(image.shape, block, axis, patch_per_block)
+    args = (image, output, weights, anchor) + image.shape[:ndim]
+    conv_axis_kernel(
+        grid,
+        block,
+        args,
+    )
+    if needs_temp:
+        output[:] = temp
+        output = temp
+    return output
diff --git a/python/cucim/src/cucim/skimage/filters/edges.py b/python/cucim/src/cucim/skimage/filters/edges.py
index 6b4dcff37..193da4f0d 100644
--- a/python/cucim/src/cucim/skimage/filters/edges.py
+++ b/python/cucim/src/cucim/skimage/filters/edges.py
@@ -13,7 +13,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import _supported_float_type, check_nD
 from ..restoration.uft import laplacian
diff --git a/python/cucim/src/cucim/skimage/filters/thresholding.py b/python/cucim/src/cucim/skimage/filters/thresholding.py
index 9f26ec982..14acc342e 100644
--- a/python/cucim/src/cucim/skimage/filters/thresholding.py
+++ b/python/cucim/src/cucim/skimage/filters/thresholding.py
@@ -6,7 +6,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 from skimage.filters import threshold_isodata as _threshold_isodata_cpu
 from skimage.filters import threshold_minimum as _threshold_minimum_cpu
 from skimage.filters import threshold_multiotsu as _threshold_multiotsu_cpu
diff --git a/python/cucim/src/cucim/skimage/measure/_blur_effect.py b/python/cucim/src/cucim/skimage/measure/_blur_effect.py
index 4e0fe6fe9..69421512e 100644
--- a/python/cucim/src/cucim/skimage/measure/_blur_effect.py
+++ b/python/cucim/src/cucim/skimage/measure/_blur_effect.py
@@ -1,5 +1,5 @@
 import cupy as cp
-import cupyx.scipy.ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from ..color import rgb2gray
 from ..util import img_as_float
diff --git a/python/cucim/src/cucim/skimage/measure/profile.py b/python/cucim/src/cucim/skimage/measure/profile.py
index 6643dc034..fc0b0cf1f 100644
--- a/python/cucim/src/cucim/skimage/measure/profile.py
+++ b/python/cucim/src/cucim/skimage/measure/profile.py
@@ -2,7 +2,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import _fix_ndimage_mode, _validate_interpolation_order
 
diff --git a/python/cucim/src/cucim/skimage/metrics/_structural_similarity.py b/python/cucim/src/cucim/skimage/metrics/_structural_similarity.py
index 6bfa07d31..44d0e28e5 100644
--- a/python/cucim/src/cucim/skimage/metrics/_structural_similarity.py
+++ b/python/cucim/src/cucim/skimage/metrics/_structural_similarity.py
@@ -1,7 +1,7 @@
 import functools
 
 import cupy as cp
-from cupyx.scipy.ndimage import uniform_filter
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared import utils
 from .._shared.filters import gaussian
@@ -187,7 +187,7 @@ def structural_similarity(im1, im2,
         filter_func = gaussian
         filter_args = {'sigma': sigma, 'truncate': truncate, 'mode': 'reflect'}
     else:
-        filter_func = uniform_filter
+        filter_func = ndi.uniform_filter
         filter_args = {'size': win_size}
 
     # ndimage filters need floating point data
diff --git a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
index 17ecd4a54..5038bedce 100644
--- a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
+++ b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
@@ -1,5 +1,5 @@
 import cupy as cp
-import cupyx.scipy.ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 import numpy as np
 
 from .._shared.utils import check_nD, deprecate_kwarg
diff --git a/python/cucim/src/cucim/skimage/morphology/binary.py b/python/cucim/src/cucim/skimage/morphology/binary.py
index 705c577cb..1df98fe00 100644
--- a/python/cucim/src/cucim/skimage/morphology/binary.py
+++ b/python/cucim/src/cucim/skimage/morphology/binary.py
@@ -4,7 +4,7 @@
 import functools
 
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import deprecate_kwarg
 from .footprints import _footprint_is_sequence
diff --git a/python/cucim/src/cucim/skimage/morphology/gray.py b/python/cucim/src/cucim/skimage/morphology/gray.py
index 96bb4c2bd..ce5daf759 100644
--- a/python/cucim/src/cucim/skimage/morphology/gray.py
+++ b/python/cucim/src/cucim/skimage/morphology/gray.py
@@ -4,7 +4,7 @@
 import functools
 
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import deprecate_kwarg
 from ..util import crop
diff --git a/python/cucim/src/cucim/skimage/morphology/misc.py b/python/cucim/src/cucim/skimage/morphology/misc.py
index d9b5304ea..bf3215a3f 100644
--- a/python/cucim/src/cucim/skimage/morphology/misc.py
+++ b/python/cucim/src/cucim/skimage/morphology/misc.py
@@ -2,7 +2,7 @@
 import functools
 
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import remove_arg, warn
 
diff --git a/python/cucim/src/cucim/skimage/restoration/j_invariant.py b/python/cucim/src/cucim/skimage/restoration/j_invariant.py
index 19f86a444..2add11e72 100644
--- a/python/cucim/src/cucim/skimage/restoration/j_invariant.py
+++ b/python/cucim/src/cucim/skimage/restoration/j_invariant.py
@@ -3,7 +3,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import _supported_float_type
 from ..metrics import mean_squared_error
diff --git a/python/cucim/src/cucim/skimage/segmentation/boundaries.py b/python/cucim/src/cucim/skimage/segmentation/boundaries.py
index 62724d9fc..d631c28e9 100644
--- a/python/cucim/src/cucim/skimage/segmentation/boundaries.py
+++ b/python/cucim/src/cucim/skimage/segmentation/boundaries.py
@@ -1,5 +1,5 @@
 import cupy as cp
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import _supported_float_type
 from ..color import gray2rgb
diff --git a/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py b/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
index d884df83d..dd3f60a82 100644
--- a/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
+++ b/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
@@ -4,7 +4,7 @@
 import cupy as cp
 import numpy as np
 from cupyx import rsqrt
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from cucim import _misc
 
diff --git a/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py b/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
index 8f2130be1..0ec9ed408 100644
--- a/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
+++ b/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
@@ -9,7 +9,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 from cupyx.scipy import sparse
 from cupyx.scipy.sparse.linalg import cg, spsolve
 
diff --git a/python/cucim/src/cucim/skimage/transform/_warps.py b/python/cucim/src/cucim/skimage/transform/_warps.py
index 3e8ae1af8..7f2f69b2b 100644
--- a/python/cucim/src/cucim/skimage/transform/_warps.py
+++ b/python/cucim/src/cucim/skimage/transform/_warps.py
@@ -2,7 +2,7 @@
 
 import cupy as cp
 import numpy as np
-from cupyx.scipy import ndimage as ndi
+import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import (_to_ndimage_mode, _validate_interpolation_order,
                              channel_as_last_axis, convert_to_float,


Test patch:
diff --git a/python/cucim/src/cucim/skimage/filters/tests/test_separable_filtering.py b/python/cucim/src/cucim/skimage/filters/tests/test_separable_filtering.py
new file mode 100644
index 000000000..0f2c5897f
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/filters/tests/test_separable_filtering.py
@@ -0,0 +1,333 @@
+import cupy as cp
+import pytest
+
+from cucim.skimage._vendored.ndimage import (convolve1d, correlate1d,
+                                             gaussian_filter, gaussian_filter1d,
+                                             gaussian_gradient_magnitude,
+                                             gaussian_laplace, laplace, prewitt,
+                                             sobel, uniform_filter,
+                                             uniform_filter1d)
+
+
+def _get_image(shape, dtype, seed=123):
+    rng = cp.random.default_rng(seed)
+    dtype = cp.dtype(dtype)
+    if dtype.kind == 'b':
+        image = rng.integers(0, 1, shape, dtype=cp.uint8).astype(bool)
+    elif dtype.kind in 'iu':
+        image = rng.integers(0, 128, shape, dtype=dtype)
+    elif dtype.kind in 'c':
+        real_dtype = cp.asarray([], dtype=dtype).real.dtype
+        image = rng.standard_normal(shape, dtype=real_dtype)
+        image = image + 1j * rng.standard_normal(shape, dtype=real_dtype)
+    else:
+        if dtype == cp.float16:
+            image = rng.standard_normal(shape).astype(dtype)
+        else:
+            image = rng.standard_normal(shape, dtype=dtype)
+    return image
+
+
+def _get_rtol_atol(dtype):
+    real_dtype = cp.array([], dtype=dtype).real.dtype
+    rtol = atol = 1e-5
+    if real_dtype == cp.float64:
+        rtol = atol = 1e-12
+    elif real_dtype == cp.float16:
+        rtol = atol = 1e-3
+    return rtol, atol
+
+
+def _compare_implementations(
+    shape, kernel_size, axis, dtype, mode, cval=0.0, origin=0,
+    output_dtype=None, kernel_dtype=None, output_preallocated=False,
+    function=convolve1d,
+):
+    dtype = cp.dtype(dtype)
+    if kernel_dtype is None:
+        kernel_dtype = dtype
+    image = _get_image(shape, dtype)
+    kernel = _get_image((kernel_size,), kernel_dtype)
+    rtol, atol = _get_rtol_atol(kernel.dtype)
+    kwargs = dict(axis=axis, mode=mode, cval=cval, origin=origin)
+    if output_dtype is not None:
+        output_dtype = cp.dtype(output_dtype)
+    if output_preallocated:
+        if output_dtype is None:
+            output_dtype = image.dtype
+        output1 = cp.empty(image.shape, dtype=output_dtype)
+        output2 = cp.empty(image.shape, dtype=output_dtype)
+        function(
+            image, kernel, output=output1, algorithm='elementwise', **kwargs
+        )
+        function(
+            image, kernel, output=output2, algorithm='shared_memory', **kwargs
+        )
+        cp.testing.assert_allclose(output1, output2, rtol=rtol, atol=atol)
+        return
+    output1 = function(
+        image, kernel, output=output_dtype, algorithm='elementwise', **kwargs
+    )
+    output2 = function(
+        image, kernel, output=output_dtype, algorithm='shared_memory', **kwargs
+    )
+    cp.testing.assert_allclose(output1, output2, rtol=rtol, atol=atol)
+    return
+
+
+def _compare_implementations_other(
+    shape, dtype, mode, cval=0.0,
+    output_dtype=None, kernel_dtype=None, output_preallocated=False,
+    function=convolve1d, func_kwargs={},
+):
+    dtype = cp.dtype(dtype)
+    image = _get_image(shape, dtype)
+    rtol, atol = _get_rtol_atol(image.dtype)
+    kwargs = dict(mode=mode, cval=cval)
+    if func_kwargs:
+        kwargs.update(func_kwargs)
+    if output_dtype is not None:
+        output_dtype = cp.dtype(output_dtype)
+    if output_preallocated:
+        if output_dtype is None:
+            output_dtype = image.dtype
+        output1 = cp.empty(image.shape, dtype=output_dtype)
+        output2 = cp.empty(image.shape, dtype=output_dtype)
+        function(image, output=output1, algorithm='elementwise', **kwargs)
+        function(image, output=output2, algorithm='shared_memory', **kwargs)
+        cp.testing.assert_allclose(output1, output2, rtol=rtol, atol=atol)
+        return
+    output1 = function(
+        image, output=output_dtype, algorithm='elementwise', **kwargs
+    )
+    output2 = function(
+        image, output=output_dtype, algorithm='shared_memory', **kwargs
+    )
+    cp.testing.assert_allclose(output1, output2, rtol=rtol, atol=atol)
+    return
+
+
+@pytest.mark.parametrize('shape', ((64, 57), (1000, 500)))
+@pytest.mark.parametrize('axis', (0, 1))
+@pytest.mark.parametrize('origin', ('min', 0, 'max'))
+@pytest.mark.parametrize('kernel_size', tuple(range(1, 17)))
+@pytest.mark.parametrize('function', [convolve1d, correlate1d])
+def test_separable_kernel_sizes_and_origins(
+    shape, axis, origin, kernel_size, function
+):
+    if kernel_size == 1:
+        origin = 0
+    elif origin == 'min':
+        origin = -(kernel_size // 2)
+    elif origin == 'max':
+        origin = kernel_size // 2
+        if kernel_size % 2 == 0:
+            origin -= 1
+    _compare_implementations(
+        shape,
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode='nearest',
+        origin=origin,
+        function=function,
+    )
+
+
+@pytest.mark.parametrize('shape', ((64, 57), (1000, 500)))
+@pytest.mark.parametrize('axis', (0, 1))
+@pytest.mark.parametrize(
+    'kernel_size',
+    tuple(range(17, 129, 11)) + tuple(range(145, 275, 41))
+)
+def test_separable_kernel_larger_sizes(shape, axis, kernel_size):
+    _compare_implementations(
+        shape,
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode='reflect',
+        origin=0,
+    )
+
+
+@pytest.mark.parametrize('shape', ((1000, 500),))
+@pytest.mark.parametrize('axis', (0, 1))
+def test_separable_elementwise_very_large_size_fallback(shape, axis):
+    """Very large kernel to make it likely shared memory will be exceeded."""
+    _compare_implementations(
+        shape,
+        kernel_size=901,
+        axis=axis,
+        dtype=cp.float64,
+        mode='nearest',
+        origin=0,
+    )
+
+
+@pytest.mark.parametrize('shape', ((4000, 2000), (1, 1), (5, 500), (1500, 5)))
+@pytest.mark.parametrize('axis', (-1, -2))
+@pytest.mark.parametrize('kernel_size', (1, 38, 129))
+@pytest.mark.parametrize(
+    'mode',
+    ('nearest', 'reflect', 'wrap', 'mirror', 'constant', ('constant', 1)),
+)
+def test_separable_image_shapes_and_modes(shape, axis, kernel_size, mode):
+
+    if isinstance(mode, tuple):
+        mode, cval = mode
+    else:
+        cval = 0
+
+    _compare_implementations(
+        shape,
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode=mode,
+        cval=cval,
+        origin=0,
+    )
+
+
+image_dtypes_tested = (
+    cp.float16, cp.float32, cp.float64, cp.complex64, cp.complex128, bool,
+    cp.int8, cp.uint8, cp.int16, cp.uint16, cp.int32, cp.uint32, cp.int64,
+    cp.uint64,
+)
+
+
+@pytest.mark.parametrize('axis', (0, 1))
+@pytest.mark.parametrize('image_dtype', image_dtypes_tested)
+@pytest.mark.parametrize(
+    'kernel_dtype', (None, cp.float32, cp.uint8, cp.complex64)
+)
+def test_separable_image_and_kernel_dtypes(axis, image_dtype, kernel_dtype):
+    """Test many kernel and image dtype combinations"""
+
+    _compare_implementations(
+        (64, 32),
+        kernel_size=3,
+        axis=axis,
+        dtype=image_dtype,
+        mode='nearest',
+        origin=0,
+        kernel_dtype=kernel_dtype,
+    )
+
+
+@pytest.mark.parametrize('axis', (0, 1))
+@pytest.mark.parametrize('image_dtype', image_dtypes_tested)
+@pytest.mark.parametrize(
+    'output_dtype', (None, cp.float32, cp.int32, cp.complex64)
+)
+@pytest.mark.parametrize('output_preallocated', (False, True))
+def test_separable_input_and_output_dtypes(
+    axis, image_dtype, output_dtype, output_preallocated
+):
+    """Test many kernel and image dtype combinations"""
+    if cp.dtype(image_dtype).kind == 'c' and output_dtype is not None:
+        if not cp.dtype(output_dtype).kind == 'c':
+            pytest.skip('cannot cast complex values to real')
+    _compare_implementations(
+        (64, 32),
+        kernel_size=3,
+        axis=axis,
+        dtype=image_dtype,
+        mode='nearest',
+        origin=0,
+        kernel_dtype=None,
+        output_dtype=output_dtype,
+        output_preallocated=output_preallocated,
+    )
+
+
+@pytest.mark.parametrize('shape', ((64, 57),))
+@pytest.mark.parametrize('axis', (0, 1))
+@pytest.mark.parametrize('origin', ('min', 0, 'max'))
+@pytest.mark.parametrize(
+    'function, func_kwargs',
+    [
+        (gaussian_filter, dict(sigma=1.5)),
+        (gaussian_filter1d, dict(sigma=1.5, axis=0)),
+        (gaussian_filter1d, dict(sigma=1.5, axis=1)),
+        (gaussian_gradient_magnitude, dict(sigma=3.5)),
+        (gaussian_laplace, dict(sigma=2.5)),
+        (laplace, {}),
+        (prewitt, {}),
+        (sobel, {}),
+        (uniform_filter, dict(size=7)),
+        (uniform_filter1d, dict(size=7, axis=0)),
+        (uniform_filter1d, dict(size=7, axis=1)),
+    ]
+)
+def test_separable_internal_kernel(
+    shape, axis, origin, function, func_kwargs
+):
+    """
+    Test case to make sure the 'algorithm' kwarg works for all other separable
+    ndimage filters as well.
+    """
+    _compare_implementations_other(
+        shape,
+        dtype=cp.float32,
+        mode='nearest',
+        function=function,
+        func_kwargs=func_kwargs,
+    )
+
+
+@pytest.mark.parametrize('shape', ((16, 24, 32), (192, 128, 160)))
+@pytest.mark.parametrize('axis', (0, 1, 2))
+@pytest.mark.parametrize('kernel_size', tuple(range(1, 17, 3)))
+@pytest.mark.parametrize('function', [convolve1d, correlate1d])
+def test_separable_kernel_sizes_3d(
+    shape, axis, kernel_size, function
+):
+    _compare_implementations(
+        shape,
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode='nearest',
+        origin=0,
+        function=function,
+    )
+
+
+@pytest.mark.parametrize('axis', (0, 1, 2))
+@pytest.mark.parametrize('kernel_size', (65, 129, 198))
+def test_separable_large_kernel_3d(axis, kernel_size):
+    _compare_implementations(
+        shape=(256, 128, 96),
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode='reflect',
+        origin=0,
+    )
+
+
+@pytest.mark.parametrize(
+    'shape', ((64, 5, 64), (5, 64, 64), (64, 64, 5), (32, 32, 32))
+)
+@pytest.mark.parametrize('axis', (-1, -2, -3))
+@pytest.mark.parametrize('kernel_size', (9,))
+@pytest.mark.parametrize(
+    'mode',
+    ('nearest', 'reflect', 'wrap', 'mirror', 'constant', ('constant', 1)),
+)
+def test_separable_image_shapes_and_modes_3d(shape, axis, kernel_size, mode):
+    if isinstance(mode, tuple):
+        mode, cval = mode
+    else:
+        cval = 0
+    _compare_implementations(
+        shape,
+        kernel_size=kernel_size,
+        axis=axis,
+        dtype=cp.float32,
+        mode=mode,
+        cval=cval,
+        origin=0,
+    )
diff --git a/python/cucim/src/cucim/skimage/measure/tests/test_blur_effect.py b/python/cucim/src/cucim/skimage/measure/tests/test_blur_effect.py
index ff5b2ac5e..3b0013005 100644
--- a/python/cucim/src/cucim/skimage/measure/tests/test_blur_effect.py
+++ b/python/cucim/src/cucim/skimage/measure/tests/test_blur_effect.py
@@ -1,6 +1,7 @@
 import cupy as cp
+import pytest
 from cupy.testing import assert_array_equal
-from skimage.data import astronaut, cells3d
+from skimage.data import astronaut
 
 from cucim.skimage.color import rgb2gray
 from cucim.skimage.filters import gaussian
@@ -47,6 +48,7 @@ def test_blur_effect_channel_axis():
 
 def test_blur_effect_3d():
     """Test that the blur metric works on a 3D image."""
+    cells3d = pytest.importorskip('skimage.data.cells3d')
     image_3d = cp.array(cells3d()[:, 1, :, :])  # grab just the nuclei
     B0 = blur_effect(image_3d)
     B1 = blur_effect(gaussian(image_3d, sigma=1))
diff --git a/python/cucim/src/cucim/skimage/registration/tests/test_masked_phase_cross_correlation.py b/python/cucim/src/cucim/skimage/registration/tests/test_masked_phase_cross_correlation.py
index 0bafa2c2f..404f64132 100644
--- a/python/cucim/src/cucim/skimage/registration/tests/test_masked_phase_cross_correlation.py
+++ b/python/cucim/src/cucim/skimage/registration/tests/test_masked_phase_cross_correlation.py
@@ -4,7 +4,7 @@
 from cupyx.scipy.ndimage import fourier_shift
 from cupyx.scipy.ndimage import shift as real_shift
 from numpy.testing import assert_almost_equal
-from skimage.data import brain, camera
+from skimage.data import camera
 from skimage.io import imread
 
 from cucim.skimage._shared.fft import fftmodule as fft
@@ -67,6 +67,7 @@ def test_masked_registration_random_masks():
 def test_masked_registration_3d_contiguous_mask():
     """masked_register_translation should be able to register translations
     between volumes with contiguous masks."""
+    brain = pytest.importorskip('skimage.data.brain')
     ref_vol = cp.array(brain()[:, ::2, ::2])
 
     offset = (1, -5, 10)