diff --git a/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py b/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py
index 8327aa6cb..85aacd06e 100644
--- a/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py
+++ b/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py
@@ -5,6 +5,7 @@
 
 import cupy
 
+from cucim.skimage._vendored import pad
 from cucim.skimage._vendored._ndimage_util import _get_inttype
 
 try:
@@ -352,7 +353,7 @@ def _pba_2d(arr, sampling=None, return_distances=True, return_indices=False,
     orig_sy, orig_sx = arr.shape
     padding_width = _determine_padding(arr.shape, padded_size, block_size)
     if padding_width is not None:
-        arr = cupy.pad(arr, padding_width, mode="constant", constant_values=1)
+        arr = pad(arr, padding_width, mode="constant", constant_values=1)
     size = arr.shape[0]
 
     input_arr = _pack_int2(arr, marker=marker, int_dtype=int_dtype)
diff --git a/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py b/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py
index b62706546..f7c572bfe 100644
--- a/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py
+++ b/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py
@@ -4,6 +4,7 @@
 import cupy
 import numpy as np
 
+from cucim.skimage._vendored import pad
 from cucim.skimage._vendored._ndimage_util import _get_inttype
 
 from ._pba_2d import (_check_distances, _check_indices,
@@ -366,7 +367,7 @@ def _pba_3d(arr, sampling=None, return_distances=True, return_indices=False,
         arr.shape, block_size, m1, m2, m3, blockx, blocky
     )
     if padding_width is not None:
-        arr = cupy.pad(arr, padding_width, mode="constant", constant_values=1)
+        arr = pad(arr, padding_width, mode="constant", constant_values=1)
     size = arr.shape[0]
 
     # pba algorithm was implemented to use 32-bit integer to store compressed
diff --git a/python/cucim/src/cucim/skimage/_vendored/__init__.py b/python/cucim/src/cucim/skimage/_vendored/__init__.py
index 1f7c85383..8b74f933d 100644
--- a/python/cucim/src/cucim/skimage/_vendored/__init__.py
+++ b/python/cucim/src/cucim/skimage/_vendored/__init__.py
@@ -4,4 +4,5 @@
 """
 
 from cucim.skimage._vendored._pearsonr import pearsonr
+from cucim.skimage._vendored.pad import pad
 from cucim.skimage._vendored.signaltools import *  # noqa
diff --git a/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py
index 18275213b..ab396a0e5 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_interpolation.py
@@ -11,6 +11,7 @@
 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 import pad
 from cucim.skimage._vendored._internal import _normalize_axis_index, prod
 
 
@@ -223,7 +224,7 @@ def _prepad_for_spline_filter(input, mode, cval):
             kwargs = dict(mode='constant', constant_values=cval)
         else:
             kwargs = dict(mode='edge')
-        padded = cupy.pad(input, npad, **kwargs)
+        padded = pad(input, npad, **kwargs)
     else:
         npad = 0
         padded = input
@@ -236,7 +237,7 @@ def _filter_input(image, prefilter, mode, cval, order):
     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.
+    prepadding of the input with pad is used to maintain accuracy.
     ``npad`` is an integer corresponding to the amount of padding at each edge
     of the array.
     """
@@ -289,8 +290,8 @@ def map_coordinates(input, coordinates, output=None, order=3,
     _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)
+        input = pad(input, [(1, 1)] * input.ndim, 'constant',
+                    constant_values=cval)
         coordinates = cupy.add(coordinates, 1)
         mode = 'constant'
 
diff --git a/python/cucim/src/cucim/skimage/_vendored/pad.py b/python/cucim/src/cucim/skimage/_vendored/pad.py
new file mode 100644
index 000000000..3d744b431
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/pad.py
@@ -0,0 +1,813 @@
+"""version of cupy.pad that dispatches to elementwise kernels for some boundary
+modes: {'edge', 'wrap', 'symmetric', 'reflect'}.
+
+New utility _use_elementwise_kernel determines when to use the new elementwise
+kernels, otherwise the existing implementations as in cupy.pad are used.
+"""
+
+import numbers
+
+import cupy
+import numpy
+
+###############################################################################
+# Private utility functions.
+
+
+def _round_if_needed(arr, dtype):
+    """Rounds arr inplace if the destination dtype is an integer.
+    """
+    if cupy.issubdtype(dtype, cupy.integer):
+        arr.round(out=arr)  # bug in round so use rint (cupy/cupy#2330)
+
+
+def _slice_at_axis(sl, axis):
+    """Constructs a tuple of slices to slice an array in the given dimension.
+
+    Args:
+      sl(slice): The slice for the given dimension.
+      axis(int): The axis to which `sl` is applied. All other dimensions are
+          left "unsliced".
+
+    Returns:
+      tuple of slices: A tuple with slices matching `shape` in length.
+    """
+    return (slice(None),) * axis + (sl,) + (Ellipsis,)
+
+
+def _view_roi(array, original_area_slice, axis):
+    """Gets a view of the current region of interest during iterative padding.
+
+    When padding multiple dimensions iteratively corner values are
+    unnecessarily overwritten multiple times. This function reduces the
+    working area for the first dimensions so that corners are excluded.
+
+    Args:
+      array(cupy.ndarray): The array with the region of interest.
+      original_area_slice(tuple of slices): Denotes the area with original
+          values of the unpadded array.
+      axis(int): The currently padded dimension assuming that `axis` is padded
+          before `axis` + 1.
+
+    Returns:
+    """
+    axis += 1
+    sl = (slice(None),) * axis + original_area_slice[axis:]
+    return array[sl]
+
+
+def _pad_simple(array, pad_width, fill_value=None):
+    """Pads an array on all sides with either a constant or undefined values.
+
+    Args:
+      array(cupy.ndarray): Array to grow.
+      pad_width(sequence of tuple[int, int]): Pad width on both sides for each
+          dimension in `arr`.
+      fill_value(scalar, optional): If provided the padded area is
+          filled with this value, otherwise the pad area left undefined.
+          (Default value = None)
+    """
+    # Allocate grown array
+    new_shape = tuple(
+        left + size + right
+        for size, (left, right) in zip(array.shape, pad_width)
+    )
+    order = 'F' if array.flags.fnc else 'C'  # Fortran and not also C-order
+    padded = cupy.empty(new_shape, dtype=array.dtype, order=order)
+
+    if fill_value is not None:
+        padded.fill(fill_value)
+
+    # Copy old array into correct space
+    original_area_slice = tuple(
+        slice(left, left + size)
+        for size, (left, right) in zip(array.shape, pad_width)
+    )
+    padded[original_area_slice] = array
+
+    return padded, original_area_slice
+
+
+def _set_pad_area(padded, axis, width_pair, value_pair):
+    """Set an empty-padded area in given dimension.
+    """
+    left_slice = _slice_at_axis(slice(None, width_pair[0]), axis)
+    padded[left_slice] = value_pair[0]
+
+    right_slice = _slice_at_axis(
+        slice(padded.shape[axis] - width_pair[1], None), axis
+    )
+    padded[right_slice] = value_pair[1]
+
+
+def _get_edges(padded, axis, width_pair):
+    """Retrieves edge values from an empty-padded array along a given axis.
+
+    Args:
+      padded(cupy.ndarray): Empty-padded array.
+      axis(int): Dimension in which the edges are considered.
+      width_pair((int, int)): Pair of widths that mark the pad area on both
+          sides in the given dimension.
+    """
+    left_index = width_pair[0]
+    left_slice = _slice_at_axis(slice(left_index, left_index + 1), axis)
+    left_edge = padded[left_slice]
+
+    right_index = padded.shape[axis] - width_pair[1]
+    right_slice = _slice_at_axis(slice(right_index - 1, right_index), axis)
+    right_edge = padded[right_slice]
+
+    return left_edge, right_edge
+
+
+def _get_linear_ramps(padded, axis, width_pair, end_value_pair):
+    """Constructs linear ramps for an empty-padded array along a given axis.
+
+    Args:
+      padded(cupy.ndarray): Empty-padded array.
+      axis(int): Dimension in which the ramps are constructed.
+      width_pair((int, int)): Pair of widths that mark the pad area on both
+          sides in the given dimension.
+      end_value_pair((scalar, scalar)): End values for the linear ramps which
+          form the edge of the fully padded array. These values are included in
+          the linear ramps.
+    """
+    edge_pair = _get_edges(padded, axis, width_pair)
+
+    left_ramp = cupy.linspace(
+        start=end_value_pair[0],
+        # squeeze axis replaced by linspace
+        stop=edge_pair[0].squeeze(axis),
+        num=width_pair[0],
+        endpoint=False,
+        dtype=padded.dtype,
+        axis=axis,
+    )
+
+    right_ramp = cupy.linspace(
+        start=end_value_pair[1],
+        # squeeze axis replaced by linspace
+        stop=edge_pair[1].squeeze(axis),
+        num=width_pair[1],
+        endpoint=False,
+        dtype=padded.dtype,
+        axis=axis,
+    )
+    # Reverse linear space in appropriate dimension
+    right_ramp = right_ramp[_slice_at_axis(slice(None, None, -1), axis)]
+
+    return left_ramp, right_ramp
+
+
+def _get_stats(padded, axis, width_pair, length_pair, stat_func):
+    """Calculates a statistic for an empty-padded array along a given axis.
+
+    Args:
+      padded(cupy.ndarray): Empty-padded array.
+      axis(int): Dimension in which the statistic is calculated.
+      width_pair((int, int)): Pair of widths that mark the pad area on both
+          sides in the given dimension.
+      length_pair(2-element sequence of None or int): Gives the number of
+          values in valid area from each side that is taken into account when
+          calculating the statistic. If None the entire valid area in `padded`
+          is considered.
+      stat_func(function): Function to compute statistic. The expected
+          signature is
+          ``stat_func(x: ndarray, axis: int, keepdims: bool) -> ndarray``.
+    """
+    # Calculate indices of the edges of the area with original values
+    left_index = width_pair[0]
+    right_index = padded.shape[axis] - width_pair[1]
+    # as well as its length
+    max_length = right_index - left_index
+
+    # Limit stat_lengths to max_length
+    left_length, right_length = length_pair
+    if left_length is None or max_length < left_length:
+        left_length = max_length
+    if right_length is None or max_length < right_length:
+        right_length = max_length
+
+    # Calculate statistic for the left side
+    left_slice = _slice_at_axis(
+        slice(left_index, left_index + left_length), axis
+    )
+    left_chunk = padded[left_slice]
+    left_stat = stat_func(left_chunk, axis=axis, keepdims=True)
+    _round_if_needed(left_stat, padded.dtype)
+
+    if left_length == right_length == max_length:
+        # return early as right_stat must be identical to left_stat
+        return left_stat, left_stat
+
+    # Calculate statistic for the right side
+    right_slice = _slice_at_axis(
+        slice(right_index - right_length, right_index), axis
+    )
+    right_chunk = padded[right_slice]
+    right_stat = stat_func(right_chunk, axis=axis, keepdims=True)
+    _round_if_needed(right_stat, padded.dtype)
+    return left_stat, right_stat
+
+
+def _set_reflect_both(padded, axis, width_pair, method, include_edge=False):
+    """Pads an `axis` of `arr` using reflection.
+
+    Args:
+      padded(cupy.ndarray): Input array of arbitrary shape.
+      axis(int): Axis along which to pad `arr`.
+      width_pair((int, int)): Pair of widths that mark the pad area on both
+          sides in the given dimension.
+      method(str): Controls method of reflection; options are 'even' or 'odd'.
+      include_edge(bool, optional): If true, edge value is included in
+          reflection, otherwise the edge value forms the symmetric axis to the
+          reflection. (Default value = False)
+    """
+    left_pad, right_pad = width_pair
+    old_length = padded.shape[axis] - right_pad - left_pad
+
+    if include_edge:
+        # Edge is included, we need to offset the pad amount by 1
+        edge_offset = 1
+    else:
+        edge_offset = 0  # Edge is not included, no need to offset pad amount
+        old_length -= 1  # but must be omitted from the chunk
+
+    if left_pad > 0:
+        # Pad with reflected values on left side:
+        # First limit chunk size which can't be larger than pad area
+        chunk_length = min(old_length, left_pad)
+        # Slice right to left, stop on or next to edge, start relative to stop
+        stop = left_pad - edge_offset
+        start = stop + chunk_length
+        left_slice = _slice_at_axis(slice(start, stop, -1), axis)
+        left_chunk = padded[left_slice]
+
+        if method == 'odd':
+            # Negate chunk and align with edge
+            edge_slice = _slice_at_axis(slice(left_pad, left_pad + 1), axis)
+            left_chunk = 2 * padded[edge_slice] - left_chunk
+
+        # Insert chunk into padded area
+        start = left_pad - chunk_length
+        stop = left_pad
+        pad_area = _slice_at_axis(slice(start, stop), axis)
+        padded[pad_area] = left_chunk
+        # Adjust pointer to left edge for next iteration
+        left_pad -= chunk_length
+
+    if right_pad > 0:
+        # Pad with reflected values on right side:
+        # First limit chunk size which can't be larger than pad area
+        chunk_length = min(old_length, right_pad)
+        # Slice right to left, start on or next to edge, stop relative to start
+        start = -right_pad + edge_offset - 2
+        stop = start - chunk_length
+        right_slice = _slice_at_axis(slice(start, stop, -1), axis)
+        right_chunk = padded[right_slice]
+
+        if method == 'odd':
+            # Negate chunk and align with edge
+            edge_slice = _slice_at_axis(
+                slice(-right_pad - 1, -right_pad), axis
+            )
+            right_chunk = 2 * padded[edge_slice] - right_chunk
+
+        # Insert chunk into padded area
+        start = padded.shape[axis] - right_pad
+        stop = start + chunk_length
+        pad_area = _slice_at_axis(slice(start, stop), axis)
+        padded[pad_area] = right_chunk
+        # Adjust pointer to right edge for next iteration
+        right_pad -= chunk_length
+
+    return left_pad, right_pad
+
+
+def _set_wrap_both(padded, axis, width_pair):
+    """Pads an `axis` of `arr` with wrapped values.
+
+    Args:
+      padded(cupy.ndarray): Input array of arbitrary shape.
+      axis(int): Axis along which to pad `arr`.
+      width_pair((int, int)): Pair of widths that mark the pad area on both
+          sides in the given dimension.
+    """
+    left_pad, right_pad = width_pair
+    period = padded.shape[axis] - right_pad - left_pad
+
+    # If the current dimension of `arr` doesn't contain enough valid values
+    # (not part of the undefined pad area) we need to pad multiple times.
+    # Each time the pad area shrinks on both sides which is communicated with
+    # these variables.
+    new_left_pad = 0
+    new_right_pad = 0
+
+    if left_pad > 0:
+        # Pad with wrapped values on left side
+        # First slice chunk from right side of the non-pad area.
+        # Use min(period, left_pad) to ensure that chunk is not larger than
+        # pad area
+        right_slice = _slice_at_axis(
+            slice(
+                -right_pad - min(period, left_pad),
+                -right_pad if right_pad != 0 else None,
+            ),
+            axis,
+        )
+        right_chunk = padded[right_slice]
+
+        if left_pad > period:
+            # Chunk is smaller than pad area
+            pad_area = _slice_at_axis(slice(left_pad - period, left_pad), axis)
+            new_left_pad = left_pad - period
+        else:
+            # Chunk matches pad area
+            pad_area = _slice_at_axis(slice(None, left_pad), axis)
+        padded[pad_area] = right_chunk
+
+    if right_pad > 0:
+        # Pad with wrapped values on right side
+        # First slice chunk from left side of the non-pad area.
+        # Use min(period, right_pad) to ensure that chunk is not larger than
+        # pad area
+        left_slice = _slice_at_axis(
+            slice(left_pad, left_pad + min(period, right_pad)), axis
+        )
+        left_chunk = padded[left_slice]
+
+        if right_pad > period:
+            # Chunk is smaller than pad area
+            pad_area = _slice_at_axis(
+                slice(-right_pad, -right_pad + period), axis
+            )
+            new_right_pad = right_pad - period
+        else:
+            # Chunk matches pad area
+            pad_area = _slice_at_axis(slice(-right_pad, None), axis)
+        padded[pad_area] = left_chunk
+
+    return new_left_pad, new_right_pad
+
+
+def _as_pairs(x, ndim, as_index=False):
+    """Broadcasts `x` to an array with shape (`ndim`, 2).
+
+    A helper function for `pad` that prepares and validates arguments like
+    `pad_width` for iteration in pairs.
+
+    Args:
+      x(scalar or array-like, optional): The object to broadcast to the shape
+          (`ndim`, 2).
+      ndim(int): Number of pairs the broadcasted `x` will have.
+      as_index(bool, optional): If `x` is not None, try to round each
+          element of `x` to an integer (dtype `cupy.intp`) and ensure every
+          element is positive. (Default value = False)
+
+    Returns:
+      nested iterables, shape (`ndim`, 2): The broadcasted version of `x`.
+    """
+    if x is None:
+        # Pass through None as a special case, otherwise cupy.round(x) fails
+        # with an AttributeError
+        return ((None, None),) * ndim
+    elif isinstance(x, numbers.Number):
+        if as_index:
+            x = round(x)
+        return ((x, x),) * ndim
+
+    x = numpy.array(x)
+    if as_index:
+        x = numpy.asarray(numpy.round(x), dtype=numpy.intp)
+
+    if x.ndim < 3:
+        # Optimization: Possibly use faster paths for cases where `x` has
+        # only 1 or 2 elements. `numpy.broadcast_to` could handle these as well
+        # but is currently slower
+
+        if x.size == 1:
+            # x was supplied as a single value
+            x = x.ravel()  # Ensure x[0] works for x.ndim == 0, 1, 2
+            if as_index and x < 0:
+                raise ValueError("index can't contain negative values")
+            return ((x[0], x[0]),) * ndim
+
+        if x.size == 2 and x.shape != (2, 1):
+            # x was supplied with a single value for each side
+            # but except case when each dimension has a single value
+            # which should be broadcasted to a pair,
+            # e.g. [[1], [2]] -> [[1, 1], [2, 2]] not [[1, 2], [1, 2]]
+            x = x.ravel()  # Ensure x[0], x[1] works
+            if as_index and (x[0] < 0 or x[1] < 0):
+                raise ValueError("index can't contain negative values")
+            return ((x[0], x[1]),) * ndim
+
+    if as_index and x.min() < 0:
+        raise ValueError("index can't contain negative values")
+
+    # Converting the array with `tolist` seems to improve performance
+    # when iterating and indexing the result (see usage in `pad`)
+    x_view = x.view()
+    x_view.shape = (ndim, 2)
+    return x_view.tolist()
+
+
+def _use_elementwise_kernel(arr, mode, kwargs):
+    """Determine if we can use an ElementwiseKernel from pad_elementwise.py"""
+    use_elementwise = False
+    if arr.ndim == 0 or arr.size == 0:
+        return False
+    if mode in ('edge', 'wrap'):
+        use_elementwise = True
+    # elif mode == 'constant':
+    #     # Only a uniform constant is supported in the Elementwise kernel.
+    #     # A per-axis constant is not currently supported.
+    #     return isinstance(kwargs.get('constant_values', 0), numbers.Number)
+    elif mode in ('symmetric', 'reflect'):
+        # only the default 'even' reflect type is supported
+        use_elementwise = kwargs.get('reflect_type', 'even') == 'even'
+    if use_elementwise:
+        if (arr.ndim > 2 and (arr.flags.fnc and arr.nbytes > 5_000_000)):
+            # Empirically found slower performance for large Fortran-ordered
+            # arrays with ndim > 2.
+            return False
+        return True
+    return False
+
+
+###############################################################################
+# Public functions
+
+
+# @array_function_dispatch(_pad_dispatcher, module='numpy')
+def pad(array, pad_width, mode='constant', **kwargs):
+    """Pads an array with specified widths and values.
+
+    Args:
+      array(cupy.ndarray): The array to pad.
+      pad_width(sequence, array_like or int): Number of values padded to the
+          edges of each axis. ((before_1, after_1), ... (before_N, after_N))
+          unique pad widths for each axis. ((before, after),) yields same
+          before and after pad for each axis. (pad,) or int is a shortcut for
+          before = after = pad width for all axes. You cannot specify
+          ``cupy.ndarray``.
+      mode(str or function, optional): One of the following string values or a
+          user supplied function
+
+          'constant' (default)
+              Pads with a constant value.
+          'edge'
+              Pads with the edge values of array.
+          'linear_ramp'
+              Pads with the linear ramp between end_value and the array edge
+              value.
+          'maximum'
+              Pads with the maximum value of all or part of the vector along
+              each axis.
+          'mean'
+              Pads with the mean value of all or part of the vector along each
+              axis.
+          'median'
+              Pads with the median value of all or part of the vector along
+              each axis. (Not Implemented)
+          'minimum'
+              Pads with the minimum value of all or part of the vector along
+              each axis.
+          'reflect'
+              Pads with the reflection of the vector mirrored on the first and
+              last values of the vector along each axis.
+          'symmetric'
+               Pads with the reflection of the vector mirrored along the edge
+               of the array.
+          'wrap'
+              Pads with the wrap of the vector along the axis. The first
+              values are used to pad the end and the end values are used to
+              pad the beginning.
+          'empty'
+              Pads with undefined values.
+          <function>
+              Padding function, see Notes.
+      stat_length(sequence or int, optional): Used in 'maximum', 'mean',
+          'median', and 'minimum'.  Number of values at edge of each axis used
+          to calculate the statistic value.
+          ((before_1, after_1), ... (before_N, after_N)) unique statistic
+          lengths for each axis. ((before, after),) yields same before and
+          after statistic lengths for each axis. (stat_length,) or int is a
+          shortcut for before = after = statistic length for all axes.
+          Default is ``None``, to use the entire axis. You cannot specify
+          ``cupy.ndarray``.
+      constant_values(sequence or scalar, optional): Used in 'constant'. The
+          values to set the padded values for each axis.
+          ((before_1, after_1), ... (before_N, after_N)) unique pad constants
+          for each axis.
+          ((before, after),) yields same before and after constants for each
+          axis.
+          (constant,) or constant is a shortcut for before = after = constant
+          for all axes.
+          Default is 0. You cannot specify ``cupy.ndarray``.
+      end_values(sequence or scalar, optional): Used in 'linear_ramp'. The
+          values used for the ending value of the linear_ramp and that will
+          form the edge of the padded array.
+          ((before_1, after_1), ... (before_N, after_N)) unique end values
+          for each axis.
+          ((before, after),) yields same before and after end
+          values for each axis.
+          (constant,) or constant is a shortcut for before = after = constant
+          for all axes.
+          Default is 0. You cannot specify ``cupy.ndarray``.
+      reflect_type({'even', 'odd'}, optional): Used in 'reflect', and
+          'symmetric'.  The 'even' style is the default with an unaltered
+          reflection around the edge value.  For the 'odd' style, the extended
+          part of the array is created by subtracting the reflected values from
+          two times the edge value.
+
+    Returns:
+      cupy.ndarray: Padded array with shape extended by ``pad_width``.
+
+    .. note::
+        For an array with rank greater than 1, some of the padding of later
+        axes is calculated from padding of previous axes.  This is easiest to
+        think about with a rank 2 array where the corners of the padded array
+        are calculated by using padded values from the first axis.
+
+        The padding function, if used, should modify a rank 1 array in-place.
+        It has the following signature:
+
+        ``padding_func(vector, iaxis_pad_width, iaxis, kwargs)``
+
+        where
+
+        vector (cupy.ndarray)
+            A rank 1 array already padded with zeros.  Padded values are
+            ``vector[:iaxis_pad_width[0]]`` and
+            ``vector[-iaxis_pad_width[1]:]``.
+        iaxis_pad_width (tuple)
+            A 2-tuple of ints, ``iaxis_pad_width[0]`` represents the number of
+            values padded at the beginning of vector where
+            ``iaxis_pad_width[1]`` represents the number of values padded at
+            the end of vector.
+        iaxis (int)
+            The axis currently being calculated.
+        kwargs (dict)
+            Any keyword arguments the function requires.
+
+    Examples
+    --------
+    >>> a = cupy.array([1, 2, 3, 4, 5])
+    >>> cupy.pad(a, (2, 3), 'constant', constant_values=(4, 6))
+    array([4, 4, 1, ..., 6, 6, 6])
+
+    >>> cupy.pad(a, (2, 3), 'edge')
+    array([1, 1, 1, ..., 5, 5, 5])
+
+    >>> cupy.pad(a, (2, 3), 'linear_ramp', end_values=(5, -4))
+    array([ 5,  3,  1,  2,  3,  4,  5,  2, -1, -4])
+
+    >>> cupy.pad(a, (2,), 'maximum')
+    array([5, 5, 1, 2, 3, 4, 5, 5, 5])
+
+    >>> cupy.pad(a, (2,), 'mean')
+    array([3, 3, 1, 2, 3, 4, 5, 3, 3])
+
+    >>> a = cupy.array([[1, 2], [3, 4]])
+    >>> cupy.pad(a, ((3, 2), (2, 3)), 'minimum')
+    array([[1, 1, 1, 2, 1, 1, 1],
+           [1, 1, 1, 2, 1, 1, 1],
+           [1, 1, 1, 2, 1, 1, 1],
+           [1, 1, 1, 2, 1, 1, 1],
+           [3, 3, 3, 4, 3, 3, 3],
+           [1, 1, 1, 2, 1, 1, 1],
+           [1, 1, 1, 2, 1, 1, 1]])
+
+    >>> a = cupy.array([1, 2, 3, 4, 5])
+    >>> cupy.pad(a, (2, 3), 'reflect')
+    array([3, 2, 1, 2, 3, 4, 5, 4, 3, 2])
+
+    >>> cupy.pad(a, (2, 3), 'reflect', reflect_type='odd')
+    array([-1,  0,  1,  2,  3,  4,  5,  6,  7,  8])
+
+    >>> cupy.pad(a, (2, 3), 'symmetric')
+    array([2, 1, 1, 2, 3, 4, 5, 5, 4, 3])
+
+    >>> cupy.pad(a, (2, 3), 'symmetric', reflect_type='odd')
+    array([0, 1, 1, 2, 3, 4, 5, 5, 6, 7])
+
+    >>> cupy.pad(a, (2, 3), 'wrap')
+    array([4, 5, 1, 2, 3, 4, 5, 1, 2, 3])
+
+    >>> def pad_with(vector, pad_width, iaxis, kwargs):
+    ...     pad_value = kwargs.get('padder', 10)
+    ...     vector[:pad_width[0]] = pad_value
+    ...     vector[-pad_width[1]:] = pad_value
+    >>> a = cupy.arange(6)
+    >>> a = a.reshape((2, 3))
+    >>> cupy.pad(a, 2, pad_with)
+    array([[10, 10, 10, 10, 10, 10, 10],
+           [10, 10, 10, 10, 10, 10, 10],
+           [10, 10,  0,  1,  2, 10, 10],
+           [10, 10,  3,  4,  5, 10, 10],
+           [10, 10, 10, 10, 10, 10, 10],
+           [10, 10, 10, 10, 10, 10, 10]])
+    >>> cupy.pad(a, 2, pad_with, padder=100)
+    array([[100, 100, 100, 100, 100, 100, 100],
+           [100, 100, 100, 100, 100, 100, 100],
+           [100, 100,   0,   1,   2, 100, 100],
+           [100, 100,   3,   4,   5, 100, 100],
+           [100, 100, 100, 100, 100, 100, 100],
+           [100, 100, 100, 100, 100, 100, 100]])
+    """
+    if isinstance(pad_width, numbers.Integral):
+        pad_width = ((pad_width, pad_width),) * array.ndim
+    else:
+        pad_width = numpy.asarray(pad_width)
+
+        if not pad_width.dtype.kind == 'i':
+            raise TypeError('`pad_width` must be of integral type.')
+
+        # Broadcast to shape (array.ndim, 2)
+        pad_width = _as_pairs(pad_width, array.ndim, as_index=True)
+
+    if callable(mode):
+        # Old behavior: Use user-supplied function with numpy.apply_along_axis
+        function = mode
+        # Create a new zero padded array
+        padded, _ = _pad_simple(array, pad_width, fill_value=0)
+        # And apply along each axis
+
+        for axis in range(padded.ndim):
+            # Iterate using ndindex as in apply_along_axis, but assuming that
+            # function operates inplace on the padded array.
+
+            # view with the iteration axis at the end
+            view = cupy.moveaxis(padded, axis, -1)
+
+            # compute indices for the iteration axes, and append a trailing
+            # ellipsis to prevent 0d arrays decaying to scalars (gh-8642)
+            inds = numpy.ndindex(view.shape[:-1])
+            inds = (ind + (Ellipsis,) for ind in inds)
+            for ind in inds:
+                function(view[ind], pad_width[axis], axis, kwargs)
+
+        return padded
+
+    # Make sure that no unsupported keywords were passed for the current mode
+    allowed_kwargs = {
+        'empty': [],
+        'edge': [],
+        'wrap': [],
+        'constant': ['constant_values'],
+        'linear_ramp': ['end_values'],
+        'maximum': ['stat_length'],
+        'mean': ['stat_length'],
+        # 'median': ['stat_length'],
+        'minimum': ['stat_length'],
+        'reflect': ['reflect_type'],
+        'symmetric': ['reflect_type'],
+    }
+    try:
+        unsupported_kwargs = set(kwargs) - set(allowed_kwargs[mode])
+    except KeyError:
+        raise ValueError("mode '{}' is not supported".format(mode))
+    if unsupported_kwargs:
+        raise ValueError(
+            "unsupported keyword arguments for mode '{}': {}".format(
+                mode, unsupported_kwargs
+            )
+        )
+
+    if _use_elementwise_kernel(array, mode, kwargs):
+        # import here to avoid circular import
+        from cucim.skimage._vendored.pad_elementwise import _get_pad_kernel
+
+        if mode == 'reflect' and min(array.shape) > 1:
+            mode = 'reflect_no_singleton_dim'
+
+        if not array.flags.forc:
+            # make non-contiguous input C-contiguous
+            array = cupy.ascontiguousarray(array)
+
+        # Allocate grown array
+        new_shape = tuple(
+            left + size + right
+            for size, (left, right) in zip(array.shape, pad_width)
+        )
+        order = 'F' if array.flags.fnc else 'C'  # Fortran and not also C-order
+        padded = cupy.empty(new_shape, dtype=array.dtype, order=order)
+
+        (int_type, np_type) = (('int', cupy.int32) if padded.size < (1 << 31)
+                               else ('ptrdiff_t', cupy.intp))
+        kern = _get_pad_kernel(
+            pad_starts=tuple(p[0] for p in pad_width),
+            mode=mode,
+            int_type=int_type,
+            order=order,
+        )
+        # pad_width must be C-contiguous
+        if mode == 'constant':
+            # `_use_elementwise_kernel` excludes cases with non-scalar cval
+            cval = float(kwargs.get('constant_values', 0))
+            kern(array, cval, padded, size=padded.size)
+        else:
+            kern(array, padded, size=padded.size)
+        return padded
+
+    if mode == 'constant':
+        values = kwargs.get('constant_values', 0)
+        if isinstance(values, numbers.Number) and values == 0 and (
+                array.ndim == 1 or array.size < 4e6):
+            # faster path for 1d arrays or small n-dimensional arrays
+            return _pad_simple(array, pad_width, 0)[0]
+
+    stat_functions = {
+        'maximum': cupy.max,
+        'minimum': cupy.min,
+        'mean': cupy.mean,
+        # 'median': cupy.median,
+    }
+
+    # Create array with final shape and original values
+    # (padded area is undefined)
+    padded, original_area_slice = _pad_simple(array, pad_width)
+    # And prepare iteration over all dimensions
+    # (zipping may be more readable than using enumerate)
+    axes = range(padded.ndim)
+
+    if mode == 'constant':
+        values = _as_pairs(values, padded.ndim)
+        for axis, width_pair, value_pair in zip(axes, pad_width, values):
+            roi = _view_roi(padded, original_area_slice, axis)
+            _set_pad_area(roi, axis, width_pair, value_pair)
+
+    elif mode == 'empty':
+        pass  # Do nothing as _pad_simple already returned the correct result
+
+    elif array.size == 0:
+        # Only modes 'constant' and 'empty' can extend empty axes, all other
+        # modes depend on `array` not being empty
+        # -> ensure every empty axis is only 'padded with 0'
+        for axis, width_pair in zip(axes, pad_width):
+            if array.shape[axis] == 0 and any(width_pair):
+                raise ValueError(
+                    "can't extend empty axis {} using modes other than "
+                    "'constant' or 'empty'".format(axis)
+                )
+        # passed, don't need to do anything more as _pad_simple already
+        # returned the correct result
+
+    elif mode == 'edge':
+        for axis, width_pair in zip(axes, pad_width):
+            roi = _view_roi(padded, original_area_slice, axis)
+            edge_pair = _get_edges(roi, axis, width_pair)
+            _set_pad_area(roi, axis, width_pair, edge_pair)
+
+    elif mode == 'linear_ramp':
+        end_values = kwargs.get('end_values', 0)
+        end_values = _as_pairs(end_values, padded.ndim)
+        for axis, width_pair, value_pair in zip(axes, pad_width, end_values):
+            roi = _view_roi(padded, original_area_slice, axis)
+            ramp_pair = _get_linear_ramps(roi, axis, width_pair, value_pair)
+            _set_pad_area(roi, axis, width_pair, ramp_pair)
+
+    elif mode in stat_functions:
+        func = stat_functions[mode]
+        length = kwargs.get('stat_length', None)
+        length = _as_pairs(length, padded.ndim, as_index=True)
+        for axis, width_pair, length_pair in zip(axes, pad_width, length):
+            roi = _view_roi(padded, original_area_slice, axis)
+            stat_pair = _get_stats(roi, axis, width_pair, length_pair, func)
+            _set_pad_area(roi, axis, width_pair, stat_pair)
+
+    elif mode in {'reflect', 'symmetric'}:
+        method = kwargs.get('reflect_type', 'even')
+        include_edge = True if mode == 'symmetric' else False
+        for axis, (left_index, right_index) in zip(axes, pad_width):
+            if array.shape[axis] == 1 and (left_index > 0 or right_index > 0):
+                # Extending singleton dimension for 'reflect' is legacy
+                # behavior; it really should raise an error.
+                edge_pair = _get_edges(padded, axis, (left_index, right_index))
+                _set_pad_area(
+                    padded, axis, (left_index, right_index), edge_pair
+                )
+                continue
+
+            roi = _view_roi(padded, original_area_slice, axis)
+            while left_index > 0 or right_index > 0:
+                # Iteratively pad until dimension is filled with reflected
+                # values. This is necessary if the pad area is larger than
+                # the length of the original values in the current dimension.
+                left_index, right_index = _set_reflect_both(
+                    roi, axis, (left_index, right_index), method, include_edge
+                )
+
+    elif mode == 'wrap':
+        for axis, (left_index, right_index) in zip(axes, pad_width):
+            roi = _view_roi(padded, original_area_slice, axis)
+            while left_index > 0 or right_index > 0:
+                # Iteratively pad until dimension is filled with wrapped
+                # values. This is necessary if the pad area is larger than
+                # the length of the original values in the current dimension.
+                left_index, right_index = _set_wrap_both(
+                    roi, axis, (left_index, right_index)
+                )
+
+    return padded
diff --git a/python/cucim/src/cucim/skimage/_vendored/pad_elementwise.py b/python/cucim/src/cucim/skimage/_vendored/pad_elementwise.py
new file mode 100644
index 000000000..9b58c2fe0
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/pad_elementwise.py
@@ -0,0 +1,230 @@
+import cupy
+
+
+def _pad_boundary_ops(mode, var_name, size, int_t="int", no_singleton=False):
+    T = 'int' if int_t == 'int' else 'long long'
+    min_func = 'min'
+    max_func = 'max'
+    if mode == 'constant':
+        ops = f'''
+        if (({var_name} < 0) || {var_name} >= {size}) {{
+            {var_name} = -1;
+        }}'''
+    elif mode == 'symmetric':
+        ops = f'''
+            if ({var_name} < 0) {{
+                {var_name} = - 1 -{var_name};
+            }}
+            {var_name} %= {size} * 2;
+            {var_name} = {min_func}({var_name}, 2 * {size} - 1 - {var_name});
+        '''
+    elif mode == 'reflect':
+        ops = f'''
+        if ({size} == 1) {{
+            {var_name} = 0;
+        }} else {{
+            if ({var_name} < 0) {{
+                {var_name} = -{var_name};
+            }}
+            if ({var_name} >= {size}) {{
+                {var_name} = 1 + ({var_name} - 1) % (({size} - 1) * 2);
+                {var_name} = {min_func}({var_name},
+                                        2 * {size} - 2 - {var_name});
+            }}
+        }}'''  # noqa
+    elif mode == 'reflect_no_singleton_dim':
+        # the same as reflect, but without the extra `{size} == 1` check
+        ops = f'''
+        if ({var_name} < 0) {{
+            {var_name} = -{var_name};
+        }}
+        if ({var_name} >= {size}) {{
+            {var_name} = 1 + ({var_name} - 1) % (({size} - 1) * 2);
+            {var_name} = {min_func}({var_name}, 2 * {size} - 2 - {var_name});
+        }}
+        '''
+    elif mode == 'edge':
+        ops = f'''
+        {var_name} = {min_func}(
+            {max_func}(static_cast<{T}>({var_name}), static_cast<{T}>(0)),
+            static_cast<{T}>({size} - 1));
+        '''
+    elif mode == 'wrap':
+        ops = f'''
+        {var_name} %= {size};
+        if ({var_name} < 0) {{
+            {var_name} += {size};
+        }}
+        '''
+    return ops + "\n"
+
+
+def _generate_size_vars(
+    ndim, arr_name='arr', size_prefix='size', int_type='int'
+):
+    """Store shape of a raw array into individual variables.
+
+    Examples
+    --------
+    >>> print(_generate_size_vars(3, 'arr', 'size', 'int'))
+    int size_0 = arr.shape()[0];
+    int size_1 = arr.shape()[1];
+    int size_2 = arr.shape()[2];
+    """
+    set_size_vars = [f'{int_type} {size_prefix}_{i} = {arr_name}.shape()[{i}];'
+                     for i in range(ndim)]
+    return '\n'.join(set_size_vars) + '\n'
+
+
+def _generate_stride_vars(
+    ndim, arr_name='arr', size_prefix='stride', int_type='int'
+):
+    """Store stride (in bytes) of a raw array into individual variables.
+
+    Examples
+    --------
+    >>> print(_generate_size_vars(3, 'arr', 'size', 'int'))
+    int stride_0 = arr.strides()[0];
+    int stride_1 = arr.strides()[1];
+    int stride_2 = arr.strides()[2];
+    """
+    set_size_vars = [
+        f'{int_type} {size_prefix}_{i} = {arr_name}.strides()[{i}];'
+        for i in range(ndim)
+    ]
+    return '\n'.join(set_size_vars) + '\n'
+
+
+def _generate_indices_ops(
+    ndim, size_prefix='size', int_type='int', index_prefix='ind', order='C',
+):
+    """Generate indices based existing variables.
+
+    Assumes variables f'{size_prefix}_{i}' has the size along axis, i.
+
+    Examples
+    --------
+    >>> print(_generate_indices_ops(3, 'size', 'int', 'ind', 'C'))
+    int _i = i;
+    int ind_2 = _i % size_2; _i /= size_2;
+    int ind_1 = _i % size_1; _i /= size_1;
+    int ind_0 = _i;
+    """
+    if order == 'C':
+        _range = range(ndim - 1, 0, -1)
+        idx_largest_stride = 0
+    elif order == 'F':
+        _range = range(ndim - 1)
+        idx_largest_stride = ndim - 1
+    else:
+        raise ValueError(f"Unknown order: {order}. Must be one of {'C', 'F'}.")
+    body = [f'{int_type} {index_prefix}_{j} = _i % {size_prefix}_{j}; _i /= {size_prefix}_{j};'  # noqa
+            for j in _range]
+    body = '\n'.join(body)
+    code = f'{int_type} _i = i;\n'
+    code += body + '\n'
+    code += f'{int_type} {index_prefix}_{idx_largest_stride} = _i;\n'
+    return code
+
+
+def _gen_raveled(ndim, stride_prefix='stride', index_prefix='i', order=None):
+    """Generate raveled index for c-ordered memory layout
+
+    For index_prefix='i', the indices are (i_0, i_1, ....)
+    For stride_prefix='stride', the stride is (stride_0, stride_1, ....)
+    """
+    return ' + '.join(
+        f'{stride_prefix}_{j} * {index_prefix}_{j}' for j in range(ndim)
+    )
+
+
+def _get_pad_kernel_code(pad_starts, int_type='int', mode='edge', order='C'):
+    # variables storing shape of the output array
+    ndim = len(pad_starts)
+    out_size_prefix = 'shape'
+    operation = _generate_size_vars(
+        ndim, arr_name='out', size_prefix=out_size_prefix, int_type=int_type
+    )
+
+    # variables storing shape of the input array
+    in_size_prefix = 'ishape'
+    in_stride_prefix = 'istride'
+    operation += _generate_size_vars(
+        ndim, arr_name='arr', size_prefix=in_size_prefix, int_type=int_type
+    )
+    operation += _generate_stride_vars(
+        ndim, arr_name='arr', size_prefix=in_stride_prefix, int_type=int_type
+    )
+
+    # unraveled indices into the output array
+    out_index_prefix = 'oi'
+    # Note: Regardless of actual memory layout, need order='C' here to match
+    #       the behavior of the index raveling used by ElementwiseKernel.
+    operation += _generate_indices_ops(
+        ndim, size_prefix=out_size_prefix, int_type=int_type,
+        index_prefix=out_index_prefix, order='C'
+    )
+
+    # compute unraveled indices into the input array
+    # (i_0, i_1, ...)
+    in_index_prefix = 'i'
+    operation += '\n'.join(
+        [f'{int_type} {in_index_prefix}_{j} = {out_index_prefix}_{j} - {pad_starts[j]};'  # noqa
+         for j in range(ndim)]
+    )
+    operation += '\n'
+    input_indices = tuple(f'{in_index_prefix}_{j}' for j in range(ndim))
+
+    # impose boundary condition
+
+    range_cond = " || ".join(
+        f"({coord} < 0) || ({coord} >= {in_size_prefix}_{j})"
+        for j, coord in enumerate(input_indices)
+    )
+    operation += f"bool range_cond = {range_cond};"
+    operation += "if (range_cond) {\n"
+    if mode == "constant":
+        for j, coord in enumerate(input_indices):
+            operation += _pad_boundary_ops(
+                mode, coord, f"{in_size_prefix}_{j}", int_type
+            )
+            operation += f"""
+                if ({coord} == -1) {{
+                    out[i] = static_cast<F>(cval);
+                    return;
+                }}
+            """
+    else:
+        for j, coord in enumerate(input_indices):
+            operation += _pad_boundary_ops(
+                mode, coord, f"{in_size_prefix}_{j}", int_type
+            )
+    operation += "}\n"
+
+    raveled_idx = _gen_raveled(
+        ndim, stride_prefix=in_stride_prefix, index_prefix=in_index_prefix,
+        order=order,
+    )
+    operation += f"""
+    // set output based on raveled index into the input array
+    const char* char_arr = reinterpret_cast<const char*>(&arr[0]);
+    out[i] = *reinterpret_cast<const F*>(char_arr + {raveled_idx});
+    """
+    return operation
+
+
+@cupy._util.memoize(for_each_device=True)
+def _get_pad_kernel(pad_starts, int_type='int', mode='edge', order='C'):
+    in_params = "raw F arr"
+    if mode == 'constant':
+        in_params += ", float64 cval"
+
+    kernel_name = f"pad_{len(pad_starts)}d_order{order}_{mode}"
+    if int_type != "int":
+        kernel_name += f"_{int_type.replace(' ', '_')}_idx"
+
+    return cupy.ElementwiseKernel(
+        in_params=in_params,
+        out_params="raw F out",
+        operation=_get_pad_kernel_code(pad_starts, int_type, mode, order),
+        name=kernel_name)
diff --git a/python/cucim/src/cucim/skimage/exposure/_adapthist.py b/python/cucim/src/cucim/skimage/exposure/_adapthist.py
index d9dbf10c2..4d42e2ae3 100644
--- a/python/cucim/src/cucim/skimage/exposure/_adapthist.py
+++ b/python/cucim/src/cucim/skimage/exposure/_adapthist.py
@@ -27,6 +27,7 @@
 from cucim.skimage.exposure.exposure import rescale_intensity
 
 from .._shared.utils import _supported_float_type
+from .._vendored import pad
 from ..color.adapt_rgb import adapt_rgb, hsv_value
 from ..util import img_as_uint
 
@@ -141,9 +142,11 @@ def _clahe(image, kernel_size, clip_limit, nbins):
     pad_end_per_dim = [(k - s % k) % k + math.ceil(k / 2.)
                        for k, s in zip(kernel_size, image.shape)]
 
-    image = cp.pad(image, [[p_i, p_f] for p_i, p_f in
-                           zip(pad_start_per_dim, pad_end_per_dim)],
-                   mode='reflect')
+    image = pad(
+        image,
+        [[p_i, p_f] for p_i, p_f in zip(pad_start_per_dim, pad_end_per_dim)],
+        mode='reflect',
+    )
 
     # determine gray value bins
     bin_size = 1 + NR_OF_GRAY // nbins
@@ -194,9 +197,9 @@ def _clahe(image, kernel_size, clip_limit, nbins):
     hist = hist.reshape(hist_block_assembled_shape[:ndim] + (-1,))
 
     # duplicate leading mappings in each dim
-    map_array = cp.pad(hist,
-                       [(1, 1) for _ in range(ndim)] + [(0, 0)],
-                       mode='edge')
+    map_array = pad(
+        hist, [(1, 1) for _ in range(ndim)] + [(0, 0)], mode='edge'
+    )
 
     # Perform multilinear interpolation of graylevel mappings
     # using the convention described here:
diff --git a/python/cucim/src/cucim/skimage/feature/template.py b/python/cucim/src/cucim/skimage/feature/template.py
index ce3443c27..3eca3e266 100644
--- a/python/cucim/src/cucim/skimage/feature/template.py
+++ b/python/cucim/src/cucim/skimage/feature/template.py
@@ -3,6 +3,7 @@
 from cucim import _misc
 
 from .._shared.utils import _supported_float_type, check_nD
+from .._vendored import pad
 
 # from cupyx.scipy import signal
 
@@ -141,10 +142,10 @@ def match_template(image, template, pad_input=False, mode='constant',
 
     pad_width = tuple((width, width) for width in template.shape)
     if mode == 'constant':
-        image = cp.pad(image, pad_width=pad_width, mode=mode,
-                       constant_values=constant_values)
+        image = pad(image, pad_width=pad_width, mode=mode,
+                    constant_values=constant_values)
     else:
-        image = cp.pad(image, pad_width=pad_width, mode=mode)
+        image = pad(image, pad_width=pad_width, mode=mode)
 
     # Use special case for 2-D images for much better performance in
     # computation of integral images
diff --git a/python/cucim/src/cucim/skimage/feature/tests/test_blob.py b/python/cucim/src/cucim/skimage/feature/tests/test_blob.py
index 09346f6ad..91683111f 100644
--- a/python/cucim/src/cucim/skimage/feature/tests/test_blob.py
+++ b/python/cucim/src/cucim/skimage/feature/tests/test_blob.py
@@ -6,6 +6,7 @@
 from skimage.draw.draw3d import ellipsoid
 
 from cucim.skimage import feature
+from cucim.skimage._vendored import pad
 from cucim.skimage.feature import blob_dog, blob_doh, blob_log
 
 
@@ -77,9 +78,9 @@ def radius(x):
 def test_blob_dog_3d(dtype, threshold_type):
     # Testing 3D
     r = 10
-    pad = 10
+    pad_width = 10
     im3 = cp.asarray(ellipsoid(r, r, r))
-    im3 = cp.pad(im3, pad, mode='constant')
+    im3 = pad(im3, pad_width, mode='constant')
 
     if threshold_type == 'absolute':
         threshold = 0.001
@@ -99,9 +100,9 @@ def test_blob_dog_3d(dtype, threshold_type):
     b = blobs[0]
 
     assert b.shape == (4,)
-    assert b[0] == r + pad + 1
-    assert b[1] == r + pad + 1
-    assert b[2] == r + pad + 1
+    assert b[0] == r + pad_width + 1
+    assert b[1] == r + pad_width + 1
+    assert b[2] == r + pad_width + 1
     assert abs(math.sqrt(3) * b[3] - r) < 1.1
 
 
@@ -112,9 +113,9 @@ def test_blob_dog_3d(dtype, threshold_type):
 def test_blob_dog_3d_anisotropic(dtype, threshold_type):
     # Testing 3D anisotropic
     r = 10
-    pad = 10
+    pad_width = 10
     im3 = cp.asarray(ellipsoid(r / 2, r, r))
-    im3 = cp.pad(im3, pad, mode='constant')
+    im3 = pad(im3, pad_width, mode='constant')
 
     if threshold_type == 'absolute':
         threshold = 0.001
@@ -134,9 +135,9 @@ def test_blob_dog_3d_anisotropic(dtype, threshold_type):
     b = blobs[0]
 
     assert b.shape == (6,)
-    assert b[0] == r / 2 + pad + 1
-    assert b[1] == r + pad + 1
-    assert b[2] == r + pad + 1
+    assert b[0] == r / 2 + pad_width + 1
+    assert b[1] == r + pad_width + 1
+    assert b[2] == r + pad_width + 1
     assert abs(math.sqrt(3) * b[3] - r / 2) < 1.1
     assert abs(math.sqrt(3) * b[4] - r) < 1.1
     assert abs(math.sqrt(3) * b[5] - r) < 1.1
@@ -294,26 +295,26 @@ def test_blob_log_no_warnings():
 def test_blob_log_3d():
     # Testing 3D
     r = 6
-    pad = 10
+    pad_width = 10
     im3 = cp.asarray(ellipsoid(r, r, r))
-    im3 = cp.pad(im3, pad, mode='constant')
+    im3 = pad(im3, pad_width, mode='constant')
 
     blobs = blob_log(im3, min_sigma=3, max_sigma=10)
     b = blobs[0]
 
     assert b.shape == (4,)
-    assert b[0] == r + pad + 1
-    assert b[1] == r + pad + 1
-    assert b[2] == r + pad + 1
+    assert b[0] == r + pad_width + 1
+    assert b[1] == r + pad_width + 1
+    assert b[2] == r + pad_width + 1
     assert abs(math.sqrt(3) * b[3] - r) < 1
 
 
 def test_blob_log_3d_anisotropic():
     # Testing 3D anisotropic
     r = 6
-    pad = 10
+    pad_width = 10
     im3 = cp.asarray(ellipsoid(r / 2, r, r))
-    im3 = cp.pad(im3, pad, mode='constant')
+    im3 = pad(im3, pad_width, mode='constant')
 
     blobs = blob_log(
         im3,
@@ -323,9 +324,9 @@ def test_blob_log_3d_anisotropic():
 
     b = blobs[0]
     assert b.shape == (6,)
-    assert b[0] == r / 2 + pad + 1
-    assert b[1] == r + pad + 1
-    assert b[2] == r + pad + 1
+    assert b[0] == r / 2 + pad_width + 1
+    assert b[1] == r + pad_width + 1
+    assert b[2] == r + pad_width + 1
     assert abs(math.sqrt(3) * b[3] - r / 2) < 1
     assert abs(math.sqrt(3) * b[4] - r) < 1
     assert abs(math.sqrt(3) * b[5] - r) < 1
@@ -501,11 +502,10 @@ def test_blob_log_overlap_3d():
     r1, r2 = 7, 6
     pad1, pad2 = 11, 12
     blob1 = cp.asarray(ellipsoid(r1, r1, r1))
-    blob1 = cp.pad(blob1, pad1, mode='constant')
+    blob1 = pad(blob1, pad1, mode='constant')
     blob2 = cp.asarray(ellipsoid(r2, r2, r2))
-    blob2 = cp.pad(blob2, [(pad2, pad2), (pad2 - 9, pad2 + 9),
-                           (pad2, pad2)],
-                   mode='constant')
+    blob2 = pad(blob2, [(pad2, pad2), (pad2 - 9, pad2 + 9), (pad2, pad2)],
+                mode='constant')
     im3 = cp.logical_or(blob1, blob2)
 
     blobs = blob_log(im3, min_sigma=2, max_sigma=10, overlap=0.1)
diff --git a/python/cucim/src/cucim/skimage/feature/tests/test_corner.py b/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
index f2437b0d7..dc752ecd7 100644
--- a/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
+++ b/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
@@ -9,6 +9,7 @@
 from cucim.skimage import img_as_float
 from cucim.skimage._shared._warnings import expected_warnings
 from cucim.skimage._shared.utils import _supported_float_type
+from cucim.skimage._vendored import pad
 from cucim.skimage.color import rgb2gray
 from cucim.skimage.feature import (corner_foerstner, corner_harris,
                                    corner_kitchen_rosenfeld, corner_peaks,
@@ -24,9 +25,9 @@
 @pytest.fixture
 def im3d():
     r = 10
-    pad = 10
+    pad_width = 10
     im3 = draw.ellipsoid(r, r, r)
-    im3 = np.pad(im3, pad, mode='constant').astype(np.uint8)
+    im3 = np.pad(im3, pad_width, mode='constant').astype(np.uint8)
     return cp.asarray(im3)
 
 
@@ -226,9 +227,9 @@ def test_structure_tensor_eigenvalues(dtype):
 
 
 def test_structure_tensor_eigenvalues_3d():
-    image = cp.pad(cube(9), 5, mode='constant') * 1000
-    boundary = (cp.pad(cube(9), 5, mode='constant')
-                - cp.pad(cube(7), 6, mode='constant')).astype(bool)
+    image = pad(cube(9), 5, mode='constant') * 1000
+    boundary = (pad(cube(9), 5, mode='constant')
+                - pad(cube(7), 6, mode='constant')).astype(bool)
     A_elems = structure_tensor(image, sigma=0.1)
     e0, e1, e2 = structure_tensor_eigenvalues(A_elems)
     # e0 should detect facets
diff --git a/python/cucim/src/cucim/skimage/filters/_fft_based.py b/python/cucim/src/cucim/skimage/filters/_fft_based.py
index d54777610..5e1b76ba0 100644
--- a/python/cucim/src/cucim/skimage/filters/_fft_based.py
+++ b/python/cucim/src/cucim/skimage/filters/_fft_based.py
@@ -5,6 +5,7 @@
 import numpy as np
 
 from .._shared.utils import _supported_float_type
+from .._vendored import pad
 
 
 def _get_nd_butterworth_filter(shape, factor, order, high_pass, real,
@@ -157,7 +158,7 @@ def butterworth(
         raise ValueError("npad must be >= 0")
     elif npad > 0:
         center_slice = tuple(slice(npad, s + npad) for s in image.shape)
-        image = cp.pad(image, npad, mode='edge')
+        image = pad(image, npad, mode='edge')
     fft_shape = (image.shape if channel_axis is None
                  else np.delete(image.shape, channel_axis))
     is_real = cp.isrealobj(image)
diff --git a/python/cucim/src/cucim/skimage/filters/_median_hist.py b/python/cucim/src/cucim/skimage/filters/_median_hist.py
index 03e10c688..1dfbd8d09 100644
--- a/python/cucim/src/cucim/skimage/filters/_median_hist.py
+++ b/python/cucim/src/cucim/skimage/filters/_median_hist.py
@@ -7,6 +7,7 @@
 import numpy as np
 
 from .._shared.utils import _to_np_mode
+from .._vendored import pad
 
 if hasattr(math, 'prod'):
     prod = math.prod
@@ -495,7 +496,7 @@ def _median_hist(image, footprint, output=None, mode='mirror', cval=0,
             pad_kwargs = dict(mode=mode, constant_values=cval)
         else:
             pad_kwargs = dict(mode=mode)
-        image = cp.pad(image, npad, **pad_kwargs)
+        image = pad(image, npad, **pad_kwargs)
         # must update n_rows, n_cols after padding!
         n_rows, n_cols = image.shape[:2]
 
diff --git a/python/cucim/src/cucim/skimage/filters/_sparse.py b/python/cucim/src/cucim/skimage/filters/_sparse.py
index 6c893d4d9..d202ce3e9 100644
--- a/python/cucim/src/cucim/skimage/filters/_sparse.py
+++ b/python/cucim/src/cucim/skimage/filters/_sparse.py
@@ -2,6 +2,7 @@
 import numpy as np
 
 from .._shared.utils import _supported_float_type, _to_np_mode
+from .._vendored import pad
 
 
 def _validate_window_size(axis_sizes):
@@ -118,7 +119,7 @@ def correlate_sparse(image, kernel, mode='reflect'):
     else:
         np_mode = _to_np_mode(mode)
         _validate_window_size(kernel.shape)
-        padded_image = cp.pad(
+        padded_image = pad(
             image,
             [(w // 2, w // 2) for w in kernel.shape],
             mode=np_mode,
diff --git a/python/cucim/src/cucim/skimage/filters/tests/test_gaussian.py b/python/cucim/src/cucim/skimage/filters/tests/test_gaussian.py
index 3f79d774b..15cc7d1c3 100644
--- a/python/cucim/src/cucim/skimage/filters/tests/test_gaussian.py
+++ b/python/cucim/src/cucim/skimage/filters/tests/test_gaussian.py
@@ -3,6 +3,7 @@
 import pytest
 
 from cucim.skimage._shared._warnings import expected_warnings
+from cucim.skimage._vendored import pad
 from cucim.skimage.filters._gaussian import difference_of_gaussians, gaussian
 
 
@@ -182,7 +183,7 @@ def test_shared_mem_check_fix(dtype, sigma):
 def test_deprecated_automatic_channel_detection():
     rgb = cp.zeros((5, 5, 3))
     rgb[1, 1] = cp.arange(1, 4)
-    gray = cp.pad(rgb, pad_width=((0, 0), (0, 0), (1, 0)))
+    gray = pad(rgb, pad_width=((0, 0), (0, 0), (1, 0)))
 
     # Warning is raised if channel_axis is not set and shape is (M, N, 3)
     with pytest.warns(
diff --git a/python/cucim/src/cucim/skimage/filters/thresholding.py b/python/cucim/src/cucim/skimage/filters/thresholding.py
index e0ffc803c..367eb35d8 100644
--- a/python/cucim/src/cucim/skimage/filters/thresholding.py
+++ b/python/cucim/src/cucim/skimage/filters/thresholding.py
@@ -18,6 +18,7 @@
 from .._shared.filters import gaussian
 from .._shared.utils import _supported_float_type, deprecate_kwarg, warn
 from .._shared.version_requirements import require
+from .._vendored import pad
 from ..exposure import histogram
 from ..transform import integral_image
 from ..util import dtype_limits
@@ -955,8 +956,8 @@ def _mean_std(image, w):
 
     float_dtype = _supported_float_type(image.dtype)
     pad_width = tuple((k // 2 + 1, k // 2) for k in w)
-    padded = cp.pad(image.astype(float_dtype, copy=False), pad_width,
-                    mode='reflect')
+    padded = pad(image.astype(float_dtype, copy=False), pad_width,
+                 mode='reflect')
 
     # Note: keep float64 integral images for accuracy. Outputs of
     # _correlate_sparse can later be safely cast to float_dtype
diff --git a/python/cucim/src/cucim/skimage/measure/_regionprops.py b/python/cucim/src/cucim/skimage/measure/_regionprops.py
index 941b14c32..daf9de458 100644
--- a/python/cucim/src/cucim/skimage/measure/_regionprops.py
+++ b/python/cucim/src/cucim/skimage/measure/_regionprops.py
@@ -9,6 +9,8 @@
 from cupyx.scipy import ndimage as ndi
 from scipy.ndimage import find_objects as cpu_find_objects
 
+from cucim.skimage._vendored import pad
+
 from . import _moments
 from ._regionprops_utils import euler_number, perimeter, perimeter_crofton
 
@@ -479,7 +481,7 @@ def feret_diameter_max(self):
 
         # TODO: implement marching cubes, etc.
         warn("feret diameter_max currently not implemented on GPU.")
-        identity_convex_hull = cp.pad(
+        identity_convex_hull = pad(
             self.image_convex, 2, mode="constant", constant_values=0
         )
         identity_convex_hull = cp.asnumpy(identity_convex_hull)
diff --git a/python/cucim/src/cucim/skimage/measure/_regionprops_utils.py b/python/cucim/src/cucim/skimage/measure/_regionprops_utils.py
index 358d6f7da..54812a11b 100644
--- a/python/cucim/src/cucim/skimage/measure/_regionprops_utils.py
+++ b/python/cucim/src/cucim/skimage/measure/_regionprops_utils.py
@@ -4,6 +4,8 @@
 import cupyx.scipy.ndimage as ndi
 import numpy as np
 
+from cucim.skimage._vendored import pad
+
 from .._shared.utils import deprecate_kwarg
 
 # Don't allocate STREL_* on GPU as we don't know in advance which device
@@ -146,7 +148,7 @@ def euler_number(image, connectivity=None):
 
     # as image can be a label image, transform it to binary
     image = (image > 0).astype(int)
-    image = cp.pad(image, pad_width=1, mode='constant')
+    image = pad(image, pad_width=1, mode='constant')
 
     # check connectivity
     if connectivity is None:
@@ -323,7 +325,7 @@ def perimeter_crofton(image, directions=4):
 
     # as image could be a label image, transform it to binary image
     image = (image > 0).astype(cp.uint8)
-    image = cp.pad(image, pad_width=1, mode="constant")
+    image = pad(image, pad_width=1, mode="constant")
     XF = ndi.convolve(image, cp.array([[0, 0, 0], [0, 1, 4], [0, 2, 8]]),
                       mode='constant', cval=0)
 
diff --git a/python/cucim/src/cucim/skimage/measure/block.py b/python/cucim/src/cucim/skimage/measure/block.py
index f21c543c9..3db2ea742 100644
--- a/python/cucim/src/cucim/skimage/measure/block.py
+++ b/python/cucim/src/cucim/skimage/measure/block.py
@@ -1,6 +1,7 @@
 import cupy as cp
 import numpy as np
 
+from .._vendored import pad
 from ..util import view_as_blocks
 
 
@@ -84,8 +85,8 @@ def block_reduce(image, block_size=2, func=cp.sum, cval=0, func_kwargs=None):
             after_width = 0
         pad_width.append((0, after_width))
 
-    image = cp.pad(image, pad_width=pad_width, mode='constant',
-                   constant_values=cval)
+    image = pad(image, pad_width=pad_width, mode='constant',
+                constant_values=cval)
 
     blocked = view_as_blocks(image, block_size)
 
diff --git a/python/cucim/src/cucim/skimage/measure/tests/test_regionprops.py b/python/cucim/src/cucim/skimage/measure/tests/test_regionprops.py
index f6be66c4a..67763b8eb 100644
--- a/python/cucim/src/cucim/skimage/measure/tests/test_regionprops.py
+++ b/python/cucim/src/cucim/skimage/measure/tests/test_regionprops.py
@@ -11,6 +11,7 @@
 
 from cucim.skimage import transform
 from cucim.skimage._shared._warnings import expected_warnings
+from cucim.skimage._vendored import pad
 from cucim.skimage.measure import (euler_number, perimeter, perimeter_crofton,
                                    regionprops, regionprops_table)
 from cucim.skimage.measure._regionprops import \
@@ -239,7 +240,7 @@ def test_bbox():
 
 
 def test_area_bbox():
-    padded = cp.pad(SAMPLE, 5, mode='constant')
+    padded = pad(SAMPLE, 5, mode='constant')
     bbox_area = regionprops(padded)[0].area_bbox
     assert_array_almost_equal(bbox_area, SAMPLE.size)
 
@@ -391,7 +392,7 @@ def test_coordinates_scaled():
 
 
 def test_slice():
-    padded = cp.pad(SAMPLE, ((2, 4), (5, 2)), mode="constant")
+    padded = pad(SAMPLE, ((2, 4), (5, 2)), mode="constant")
     nrow, ncol = SAMPLE.shape
     result = regionprops(padded)[0].slice
     expected = (slice(2, 2 + nrow), slice(5, 5 + ncol))
diff --git a/python/cucim/src/cucim/skimage/metrics/simple_metrics.py b/python/cucim/src/cucim/skimage/metrics/simple_metrics.py
index 1b29e602f..3e1a7bdbb 100644
--- a/python/cucim/src/cucim/skimage/metrics/simple_metrics.py
+++ b/python/cucim/src/cucim/skimage/metrics/simple_metrics.py
@@ -2,6 +2,7 @@
 from cupyx.scipy.stats import entropy
 
 from .._shared.utils import _supported_float_type, check_shape_equality, warn
+from .._vendored import pad
 from ..util.dtype import dtype_range
 
 __all__ = ['mean_squared_error',
@@ -188,7 +189,7 @@ def _pad_to(arr, shape):
         raise ValueError(f'Target shape {shape} cannot be smaller than input'
                          f'shape {arr.shape} along any axis.')
     padding = [(0, s - i) for s, i in zip(shape, arr.shape)]
-    return cp.pad(arr, pad_width=padding, mode='constant', constant_values=0)
+    return pad(arr, pad_width=padding, mode='constant', constant_values=0)
 
 
 def normalized_mutual_information(image0, image1, *, bins=100):
diff --git a/python/cucim/src/cucim/skimage/morphology/gray.py b/python/cucim/src/cucim/skimage/morphology/gray.py
index f900c1607..48f87a40e 100644
--- a/python/cucim/src/cucim/skimage/morphology/gray.py
+++ b/python/cucim/src/cucim/skimage/morphology/gray.py
@@ -8,6 +8,7 @@
 import cucim.skimage._vendored.ndimage as ndi
 
 from .._shared.utils import deprecate_kwarg
+from .._vendored import pad
 from ..util import crop
 from .footprints import _footprint_is_sequence, _shape_from_sequence
 from .misc import default_footprint
@@ -144,7 +145,7 @@ def func_out(image, footprint, out=None, *args, **kwargs):
                 axis_pad_width = 0
             pad_widths.append((axis_pad_width,) * 2)
         if padding:
-            image = cp.pad(image, pad_widths, mode='edge')
+            image = pad(image, pad_widths, mode='edge')
             out_temp = cp.empty_like(image)
         else:
             out_temp = out
diff --git a/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
index e63521d6b..bab2f1887 100644
--- a/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
+++ b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
@@ -5,6 +5,7 @@
 from cucim.core.operations.morphology import distance_transform_edt
 
 from .._shared.utils import _supported_float_type, deprecate_kwarg
+from .._vendored import pad
 
 
 @cp.fuse()
@@ -23,7 +24,7 @@ def _fused_curvature(phi, x_start, x_end, y_start, y_end, ul, ur, ll, lr):
 def _cv_curvature(phi):
     """Returns the 'curvature' of a level set 'phi'.
     """
-    P = cp.pad(phi, 1, mode='edge')
+    P = pad(phi, 1, mode='edge')
     y_start = P[:-2, 1:-1]
     y_end = P[2:, 1:-1]
     x_start = P[1:-1, :-2]
@@ -103,7 +104,7 @@ def _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt):
     """Returns the variation of level set 'phi' based on algorithm parameters.
     """
     eta = 1e-16
-    P = cp.pad(phi, 1, mode='edge')
+    P = pad(phi, 1, mode='edge')
 
     x_end = P[1:-1, 2:]
     x_mid = P[1:-1, 1:-1]
diff --git a/python/cucim/src/cucim/skimage/transform/_warps.py b/python/cucim/src/cucim/skimage/transform/_warps.py
index 995961436..dfc617152 100644
--- a/python/cucim/src/cucim/skimage/transform/_warps.py
+++ b/python/cucim/src/cucim/skimage/transform/_warps.py
@@ -8,6 +8,7 @@
 from .._shared.utils import (_to_ndimage_mode, _validate_interpolation_order,
                              channel_as_last_axis, convert_to_float,
                              safe_as_int, warn)
+from .._vendored import pad
 from ..measure import block_reduce
 from ._geometric import (AffineTransform, ProjectiveTransform,
                          SimilarityTransform)
@@ -1203,14 +1204,14 @@ def _local_mean_weights(old_size, new_size, grid_mode, dtype):
         new_breaks = cp.linspace(0, old_size, num=new_size + 1, dtype=dtype)
     else:
         old, new = old_size - 1, new_size - 1
-        old_breaks = cp.pad(cp.linspace(0.5, old - 0.5, old, dtype=dtype),
-                            1, 'constant', constant_values=(0, old))
+        old_breaks = pad(cp.linspace(0.5, old - 0.5, old, dtype=dtype),
+                         1, 'constant', constant_values=(0, old))
         if new == 0:
             val = np.inf
         else:
             val = 0.5 * old / new
-        new_breaks = cp.pad(cp.linspace(val, old - val, new, dtype=dtype),
-                            1, 'constant', constant_values=(0, old))
+        new_breaks = pad(cp.linspace(val, old - val, new, dtype=dtype),
+                         1, 'constant', constant_values=(0, old))
 
     upper = cp.minimum(new_breaks[1:, np.newaxis], old_breaks[np.newaxis, 1:])
     lower = cp.maximum(new_breaks[:-1, np.newaxis],