Repository: rapidsai/cucim (language: python)

Issue:
Fix conflicts in "Forward-merge branch-22.12 to branch-23.02" [skip gpuci]
xref: https://github.com/rapidsai/cucim/pull/467
[BUG] Hue color jitter on an all gray or all white image returns a black image
## Describe the bug

Hue jitter, via `cucim.core.operations.color.color_jitter`, on images with a lot (or all) gray seems to overflow and give either black pixels, or an all-black image.

I noticed this when using augmenting images of planes in, e.g., CIFAR: `color_jitter` with a non-zero hue parameter would render some sky pixels black. I found this issue only comes up when the hue parameter is non-zero -- other jitter parameters didn't seem to affect the behavior. Further experimenting showed that inputting an image of all white or all one shade of gray resulted in `color_jitter` outputting an image of all zeros. 

More specifically:
```python
import cupy as cp
import cucim.core.operations.color as ccl

x_cupy = cp.ones((3, 64, 64)).astype(cp.uint8)*255
x_jitter = ccl.color_jitter(x_cupy,0,0,0,0.1)

#Outputs all zeros, i.e.
(x_jitter == 0).all() # is True
```
The same behavior occurs for any gray image, e.g. if I make `x_cupy` a tensor of all intensity 128.

Perhaps this is some rounding/overflow issue since `color_jitter` seems to convert all inputs to `uint8`? 

## Expected behavior 

I would expect the output to not be all black but, if the image is all gray then hue jitter should have no effect -- it should output the same gray image. The torch implementation, for instance, behaves how I would expect:
```python
import torch
import numpy as np 
import torchvision
x_np = np.ones((3, 64, 64)).astype(np.uint8)*255
x_torch = torch.tensor(x_np)
x_torch_jitter = torchvision.transforms.ColorJitter(0, 0, 0, 0.1)(x_torch)

#Jittered output is same as input image
(x_torch_jitter == 255).all() # is True
```

## Environment

Running Ubuntu 20.04. cudatoolkit 11.6.2.

Installed cuCIM with pip. 
python==3.8.13
cupy==11.2.0
cucim==22.10.00

(torch version for the above example == 1.13.0)

## Other notes

* I came across this issue when using MONAI's `RandCuCIM` transform for data augmentation, as it calls cucim's `color_jitter`.
* I've tried a number of `hue` parameters -- anything non-zero seems to create the issue. Though the docstring seems to suggest anything between [-0.5, 0.5] is accepted, including negative values, negative values in fact return `ValueError: If hue is a single number, it must be non negative.` But I suppose that's a separate issue.
[BUG] Hue color jitter on an all gray or all white image returns a black image
## Describe the bug

Hue jitter, via `cucim.core.operations.color.color_jitter`, on images with a lot (or all) gray seems to overflow and give either black pixels, or an all-black image.

I noticed this when using augmenting images of planes in, e.g., CIFAR: `color_jitter` with a non-zero hue parameter would render some sky pixels black. I found this issue only comes up when the hue parameter is non-zero -- other jitter parameters didn't seem to affect the behavior. Further experimenting showed that inputting an image of all white or all one shade of gray resulted in `color_jitter` outputting an image of all zeros. 

More specifically:
```python
import cupy as cp
import cucim.core.operations.color as ccl

x_cupy = cp.ones((3, 64, 64)).astype(cp.uint8)*255
x_jitter = ccl.color_jitter(x_cupy,0,0,0,0.1)

#Outputs all zeros, i.e.
(x_jitter == 0).all() # is True
```
The same behavior occurs for any gray image, e.g. if I make `x_cupy` a tensor of all intensity 128.

Perhaps this is some rounding/overflow issue since `color_jitter` seems to convert all inputs to `uint8`? 

## Expected behavior 

I would expect the output to not be all black but, if the image is all gray then hue jitter should have no effect -- it should output the same gray image. The torch implementation, for instance, behaves how I would expect:
```python
import torch
import numpy as np 
import torchvision
x_np = np.ones((3, 64, 64)).astype(np.uint8)*255
x_torch = torch.tensor(x_np)
x_torch_jitter = torchvision.transforms.ColorJitter(0, 0, 0, 0.1)(x_torch)

#Jittered output is same as input image
(x_torch_jitter == 255).all() # is True
```

## Environment

Running Ubuntu 20.04. cudatoolkit 11.6.2.

Installed cuCIM with pip. 
python==3.8.13
cupy==11.2.0
cucim==22.10.00

(torch version for the above example == 1.13.0)

## Other notes

* I came across this issue when using MONAI's `RandCuCIM` transform for data augmentation, as it calls cucim's `color_jitter`.
* I've tried a number of `hue` parameters -- anything non-zero seems to create the issue. Though the docstring seems to suggest anything between [-0.5, 0.5] is accepted, including negative values, negative values in fact return `ValueError: If hue is a single number, it must be non negative.` But I suppose that's a separate issue.


Discussion / hints:

Fix patch:
diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS
index 204ecac8b..745f4a5e4 100644
--- a/.github/CODEOWNERS
+++ b/.github/CODEOWNERS
@@ -10,9 +10,10 @@ python/            @rapidsai/cucim-python-codeowners
 **/cmake/          @rapidsai/cucim-cmake-codeowners
 
 #build/ops code owners
-.github/           @rapidsai/ops-codeowners 
+.github/           @rapidsai/ops-codeowners
 ci/                @rapidsai/ops-codeowners
 conda/             @rapidsai/ops-codeowners
 **/Dockerfile      @rapidsai/ops-codeowners
 **/.dockerignore   @rapidsai/ops-codeowners
 docker/            @rapidsai/ops-codeowners
+dependencies.yaml  @rapidsai/ops-codeowners
diff --git a/.github/ops-bot.yaml b/.github/ops-bot.yaml
index 0a52b6795..4d53902ba 100644
--- a/.github/ops-bot.yaml
+++ b/.github/ops-bot.yaml
@@ -6,3 +6,5 @@ branch_checker: true
 label_checker: true
 release_drafter: true
 external_contributors: false
+copy_prs: true
+recently_updated: true
diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
new file mode 100644
index 000000000..c64bfc4a7
--- /dev/null
+++ b/.github/workflows/build.yaml
@@ -0,0 +1,65 @@
+name: build
+
+on:
+  push:
+    branches:
+      - "branch-*"
+    tags:
+      - v[0-9][0-9].[0-9][0-9].[0-9][0-9]
+  workflow_dispatch:
+    inputs:
+      branch:
+        required: true
+        type: string
+      date:
+        required: true
+        type: string
+      sha:
+        required: true
+        type: string
+      build_type:
+        type: string
+        default: nightly
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  cpp-build:
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.02
+    with:
+      build_type: ${{ inputs.build_type || 'branch' }}
+      branch: ${{ inputs.branch }}
+      date: ${{ inputs.date }}
+      sha: ${{ inputs.sha }}
+  docs-build:
+    if: ${{ startsWith(github.ref, 'refs/heads/branch-') }}
+    needs: [python-build]
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/custom-job.yaml@branch-23.02
+    with:
+      build_type: branch
+      node_type: "gpu-latest-1"
+      arch: "amd64"
+      container_image: "rapidsai/ci:latest"
+      run_script: "ci/build_docs.sh"
+  python-build:
+    needs: [cpp-build]
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-build.yaml@branch-23.02
+    with:
+      build_type: ${{ inputs.build_type || 'branch' }}
+      branch: ${{ inputs.branch }}
+      date: ${{ inputs.date }}
+      sha: ${{ inputs.sha }}
+  upload-conda:
+    needs: [cpp-build, python-build]
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-upload-packages.yaml@branch-23.02
+    with:
+      build_type: ${{ inputs.build_type || 'branch' }}
+      branch: ${{ inputs.branch }}
+      date: ${{ inputs.date }}
+      sha: ${{ inputs.sha }}
diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml
new file mode 100644
index 000000000..a0433495f
--- /dev/null
+++ b/.github/workflows/pr.yaml
@@ -0,0 +1,52 @@
+name: pr
+
+on:
+  push:
+    branches:
+      - "pull-request/[0-9]+"
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  pr-builder:
+    needs:
+      - checks
+      - conda-cpp-build
+      - conda-python-build
+      - conda-python-tests
+      - docs-build
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/pr-builder.yaml@branch-23.02
+  checks:
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/checks.yaml@branch-23.02
+  conda-cpp-build:
+    needs: checks
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.02
+    with:
+      build_type: pull-request
+  conda-python-build:
+    needs: conda-cpp-build
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-build.yaml@branch-23.02
+    with:
+      build_type: pull-request
+  conda-python-tests:
+    needs: conda-python-build
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-tests.yaml@branch-23.02
+    with:
+      build_type: pull-request
+  docs-build:
+    needs: conda-python-build
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/custom-job.yaml@branch-23.02
+    with:
+      build_type: pull-request
+      node_type: "gpu-latest-1"
+      arch: "amd64"
+      container_image: "rapidsai/ci:latest"
+      run_script: "ci/build_docs.sh"
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 020613427..5f2a03745 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,20 +4,23 @@
 #   pre-commit install-hooks
 exclude: '^(\.tox|ci/templates|\.bumpversion\.cfg)(/|$)'
 repos:
-  - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: master
-    hooks:
-      - id: trailing-whitespace
-      - id: end-of-file-fixer
-      - id: debug-statements
-  - repo: https://github.com/pycqa/isort
-    rev: 5.6.4
+  # TODO: re-enable this after GitHub Actions migration, it currently makes a
+  # large number of changes that shouldn't go in the GitHub Actions PR.
+  #- repo: https://github.com/pre-commit/pre-commit-hooks
+  #  rev: v4.4.0
+  #  hooks:
+  #    - id: trailing-whitespace
+  #    - id: end-of-file-fixer
+  #    - id: debug-statements
+  - repo: https://github.com/PyCQA/isort
+    rev: 5.12.0
     hooks:
       - id: isort
         files: ^python/cucim/src/.*
         args: ["--settings-path=python/cucim/setup.cfg"]
-  - repo: https://gitlab.com/pycqa/flake8
-    rev: master
+  - repo: https://github.com/PyCQA/flake8
+    rev: 6.0.0
     hooks:
       - id: flake8
+        args: ["--config=python/cucim/setup.cfg"]
         files: ^python/cucim/.*
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2917b901b..80db73ccc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,37 @@
+# cuCIM 23.02.00 (9 Feb 2023)
+
+## 🚨 Breaking Changes
+
+- Add disambiguation option to phase_cross_correlation (skimage 0.20 feature) ([#486](https://github.com/rapidsai/cucim/pull/486)) [@grlee77](https://github.com/grlee77)
+
+## 🐛 Bug Fixes
+
+- apply bug fix to vendored ndimage code ([#494](https://github.com/rapidsai/cucim/pull/494)) [@grlee77](https://github.com/grlee77)
+- Closes #490 -- fixes bug in hue jitter ([#491](https://github.com/rapidsai/cucim/pull/491)) [@benlansdell](https://github.com/benlansdell)
+- Fix random seed used in test_3d_similarity_estimation ([#472](https://github.com/rapidsai/cucim/pull/472)) [@grlee77](https://github.com/grlee77)
+
+## 📖 Documentation
+
+- Fix documentation author ([#475](https://github.com/rapidsai/cucim/pull/475)) [@bdice](https://github.com/bdice)
+
+## 🚀 New Features
+
+- Add colocalization measures ([#488](https://github.com/rapidsai/cucim/pull/488)) [@grlee77](https://github.com/grlee77)
+- Add disambiguation option to phase_cross_correlation (skimage 0.20 feature) ([#486](https://github.com/rapidsai/cucim/pull/486)) [@grlee77](https://github.com/grlee77)
+
+## 🛠️ Improvements
+
+- Update shared workflow branches ([#501](https://github.com/rapidsai/cucim/pull/501)) [@ajschmidt8](https://github.com/ajschmidt8)
+- Update `isort` version to 5.12.0 ([#492](https://github.com/rapidsai/cucim/pull/492)) [@ajschmidt8](https://github.com/ajschmidt8)
+- Improve rank filtering performance by removing use of footprint kernel when possible ([#485](https://github.com/rapidsai/cucim/pull/485)) [@grlee77](https://github.com/grlee77)
+- use vendored version of cupy.pad with added performance optimizations ([#482](https://github.com/rapidsai/cucim/pull/482)) [@grlee77](https://github.com/grlee77)
+- add docs builds to Github Actions ([#481](https://github.com/rapidsai/cucim/pull/481)) [@AjayThorve](https://github.com/AjayThorve)
+- Update `numpy` version specifier ([#480](https://github.com/rapidsai/cucim/pull/480)) [@ajschmidt8](https://github.com/ajschmidt8)
+- Build CUDA `11.8` and Python `3.10` Packages ([#476](https://github.com/rapidsai/cucim/pull/476)) [@ajschmidt8](https://github.com/ajschmidt8)
+- Add GitHub Actions Workflows. ([#471](https://github.com/rapidsai/cucim/pull/471)) [@bdice](https://github.com/bdice)
+- Fix conflicts in &quot;Forward-merge branch-22.12 to branch-23.02&quot; ([#468](https://github.com/rapidsai/cucim/pull/468)) [@jakirkham](https://github.com/jakirkham)
+- Enable copy_prs. ([#465](https://github.com/rapidsai/cucim/pull/465)) [@bdice](https://github.com/bdice)
+
 # cuCIM 22.12.00 (8 Dec 2022)
 
 ## 🚨 Breaking Changes
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index a778ff1a0..f89d1007e 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -107,18 +107,19 @@ cd $CUCIM_HOME
 ```
 ## Local Development using Conda Environment (for gcc 9.x and nvcc 11.0.x)
 
-Conda can be used to setup GCC 9.x, CUDA Toolkit (including nvcc) 11.0.x, and other dependent libraries (as shown in `./conda/environments/env.yml`) for building cuCIM.
+Conda can be used to setup an environment which includes all of the necessary dependencies (as shown in `./conda/environments/all_cuda-118_arch-x86_64.yaml`) for building cuCIM.
 
 Otherwise, you may need to install dependencies (such as zlib, xz, yasm) through your OS's package manager (`apt`, `yum`, and so on).
 
 
 ### Creating the Conda Development Environment `cucim`
 
-Note that `./conda/environments/env.yml` is currently set to use gcc 9.x (gxx_linux-64) and CUDA 11.0.x (cudatoolkit & cudatoolkit-dev). If you want to change the version of gcc or CUDA toolkit
-package, please update `./conda/environments/env.yml` before executing the following commands.
+Note that `./conda/environments/all_cuda-118_arch-x86_64.yaml` is currently set to use specific versions of gcc (gxx_linux-64) and CUDA (cudatoolkit & cudatoolkit-dev).
+
+If you want to change the version of gcc or CUDA toolkit package, please update `./conda/environments/all_cuda-118_arch-x86_64.yaml` before executing the following commands.
 
 ```bash
-conda env create -f ./conda/environments/env.yml
+conda env create -n cucim -f ./conda/environments/all_cuda-118_arch-x86_64.yaml
 # activate the environment
 conda activate cucim
 ```
diff --git a/README.md b/README.md
index efaaa560f..9a92b8de5 100644
--- a/README.md
+++ b/README.md
@@ -70,7 +70,7 @@ pip install scipy scikit-image cupy-cuda110
 
 ### Notebooks
 
-Please check out our [Welcome](notebooks/Welcome.ipynb) notebook ([NBViewer](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-22.12/notebooks/Welcome.ipynb))
+Please check out our [Welcome](notebooks/Welcome.ipynb) notebook ([NBViewer](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-23.02/notebooks/Welcome.ipynb))
 
 #### Downloading sample images
 
diff --git a/VERSION b/VERSION
index ed6a80946..7c8431de5 100644
--- a/VERSION
+++ b/VERSION
@@ -1,1 +1,1 @@
-22.12.00
+23.02.00
diff --git a/benchmarks/skimage/cucim_measure_bench.py b/benchmarks/skimage/cucim_measure_bench.py
index c5f207b64..6a8f4d1cc 100644
--- a/benchmarks/skimage/cucim_measure_bench.py
+++ b/benchmarks/skimage/cucim_measure_bench.py
@@ -12,6 +12,7 @@
 import skimage.measure
 
 from _image_bench import ImageBench
+from cucim_metrics_bench import MetricsBench
 
 
 class LabelBench(ImageBench):
@@ -132,6 +133,32 @@ def set_args(self, dtype):
         self.args_gpu = (imaged,)
 
 
+class BinaryImagePairBench(ImageBench):
+    def set_args(self, dtype):
+        rng = cp.random.default_rng(seed=123)
+        img0d = rng.integers(0, 2, self.shape, dtype=cp.uint8).view(bool)
+        img1d = rng.integers(0, 2, self.shape, dtype=cp.uint8).view(bool)
+        img0 = cp.asnumpy(img0d)
+        img1 = cp.asnumpy(img1d)
+
+        self.args_cpu = (img0, img1)
+        self.args_gpu = (img0d, img1d)
+
+
+class MandersColocBench(ImageBench):
+    def set_args(self, dtype):
+        # image
+        imaged = cp.testing.shaped_arange(self.shape, dtype=dtype)
+        imaged = imaged / imaged.max()
+
+        # binary mask
+        rng = cp.random.default_rng(seed=123)
+        maskd = rng.integers(0, 2, self.shape, dtype=cp.uint8).view(bool)
+
+        self.args_cpu = (cp.asnumpy(imaged), cp.asnumpy(maskd))
+        self.args_gpu = (imaged, maskd)
+
+
 def main(args):
 
     pfile = "cucim_measure_results.pickle"
@@ -188,6 +215,12 @@ def main(args):
             True,
             False,
         ),  # variable block_size configured below
+
+        # binary image overlap measures
+        ("intersection_coeff", dict(mask=None), dict(), False, True),
+        ("manders_coloc_coeff", dict(mask=None), dict(), False, True),
+        ("manders_overlap_coeff", dict(mask=None), dict(), False, True),
+        ("pearson_corr_coeff", dict(mask=None), dict(), False, True),
     ]:
 
         if function_name != args.func_name:
@@ -225,8 +258,29 @@ def main(args):
                     module_gpu=cucim.skimage.measure,
                     run_cpu=run_cpu,
                 )
-                results = B.run_benchmark(duration=args.duration)
-                all_results = pd.concat([all_results, results["full"]])
+        elif function_name in ['intersection_coeff', 'manders_coloc_coeff',
+                               'manders_overlap_coeff', 'pearson_corr_coeff']:
+
+            if function_name in ["pearson_corr_coeff", "manders_overlap_coeff"]:
+                # arguments are two images of matching dtype
+                Tester = MetricsBench
+            elif function_name == "manders_coloc_coeff":
+                # arguments are one image + binary mask
+                Tester = MandersColocBench
+            else:
+                # arguments are two binary images
+                Tester = BinaryImagePairBench
+
+            B = Tester(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.measure,
+                module_gpu=cucim.skimage.measure,
+                run_cpu=run_cpu,
+            )
         else:
 
 
@@ -268,8 +322,8 @@ def main(args):
                 module_gpu=cucim.skimage.measure,
                 run_cpu=run_cpu,
             )
-            results = B.run_benchmark(duration=args.duration)
-            all_results = pd.concat([all_results, results["full"]])
+        results = B.run_benchmark(duration=args.duration)
+        all_results = pd.concat([all_results, results["full"]])
 
     fbase = os.path.splitext(pfile)[0]
     all_results.to_csv(fbase + ".csv")
@@ -280,8 +334,8 @@ def main(args):
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Benchmarking cuCIM measure functions')
-    func_name_choices = ['label', 'regionprops', 'moments', 'moments_central', 'centroid', 'inertia_tensor', 'inertia_tensor_eigvals', 'block_reduce', 'shannon_entropy', 'profile_line']
-    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    func_name_choices = ['label', 'regionprops', 'moments', 'moments_central', 'centroid', 'inertia_tensor', 'inertia_tensor_eigvals', 'block_reduce', 'shannon_entropy', 'profile_line', 'intersection_coeff', 'manders_coloc_coeff', 'manders_overlap_coeff', 'pearson_corr_coeff']
+    dtype_choices = ['bool', 'float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
     parser.add_argument('-i','--img_size', type=str, help='Size of input image', required=True)
     parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices=dtype_choices, required=True)
     parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices=func_name_choices, required=True)
diff --git a/benchmarks/skimage/run-nv-bench-measure.sh b/benchmarks/skimage/run-nv-bench-measure.sh
index 5c2b4f413..613307b6f 100755
--- a/benchmarks/skimage/run-nv-bench-measure.sh
+++ b/benchmarks/skimage/run-nv-bench-measure.sh
@@ -9,3 +9,28 @@ for shape in "${param_shape[@]}"; do
         done
     done
 done
+
+# # commenting out colocalization metrics below this point
+# # (scikit-image 0.20 is not yet officially released)
+# param_shape=(512,512 3840,2160 192,192,192)
+# param_filt=(manders_coloc_coeff manders_overlap_coeff pearson_corr_coeff)
+# param_dt=(float32 uint8)
+# for shape in "${param_shape[@]}"; do
+#     for filt in "${param_filt[@]}"; do
+#         for dt in "${param_dt[@]}"; do
+#             python cucim_measure_bench.py -f $filt -i $shape -d $dt -t 10
+#         done
+#     done
+# done
+
+# # only supports binary-valued images
+# param_shape=(512,512 3840,2160 192,192,192)
+# param_filt=(intersection_coeff)
+# param_dt=(bool)
+# for shape in "${param_shape[@]}"; do
+#     for filt in "${param_filt[@]}"; do
+#         for dt in "${param_dt[@]}"; do
+#             python cucim_measure_bench.py -f $filt -i $shape -d $dt -t 10
+#         done
+#     done
+# done
diff --git a/ci/build_cpp.sh b/ci/build_cpp.sh
new file mode 100755
index 000000000..269a0be40
--- /dev/null
+++ b/ci/build_cpp.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+# Copyright (c) 2022, NVIDIA CORPORATION.
+
+set -euo pipefail
+
+source rapids-env-update
+
+export CMAKE_GENERATOR=Ninja
+
+rapids-print-env
+
+rapids-logger "Begin cpp build"
+
+rapids-mamba-retry mambabuild conda/recipes/libcucim
+
+rapids-upload-conda-to-s3 cpp
diff --git a/ci/build_docs.sh b/ci/build_docs.sh
new file mode 100755
index 000000000..6dc8e9841
--- /dev/null
+++ b/ci/build_docs.sh
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+set -euo pipefail
+
+rapids-logger "Create test conda environment"
+. /opt/conda/etc/profile.d/conda.sh
+
+rapids-dependency-file-generator \
+    --output conda \
+    --file_key docs \
+    --matrix "cuda=${RAPIDS_CUDA_VERSION%.*};arch=$(arch);py=${RAPIDS_PY_VERSION}" | tee env.yaml
+
+rapids-mamba-retry env create --force -f env.yaml -n docs
+conda activate docs
+
+rapids-print-env
+
+rapids-logger "Downloading artifacts from previous jobs"
+CPP_CHANNEL=$(rapids-download-conda-from-s3 cpp)
+PYTHON_CHANNEL=$(rapids-download-conda-from-s3 python)
+VERSION_NUMBER=$(rapids-get-rapids-version-from-git)
+
+rapids-mamba-retry install \
+    --channel "${CPP_CHANNEL}" \
+    --channel "${PYTHON_CHANNEL}" \
+    cucim libcucim
+
+# Build Python docs
+rapids-logger "Build Python docs"
+pushd docs
+sphinx-build -b dirhtml ./source _html
+sphinx-build -b text ./source _text
+popd
+
+if [[ "${RAPIDS_BUILD_TYPE}" == "branch" ]]; then
+  rapids-logger "Upload Docs to S3"
+  aws s3 sync --no-progress --delete docs/_html "s3://rapidsai-docs/cucim/${VERSION_NUMBER}/html"
+  aws s3 sync --no-progress --delete docs/_text "s3://rapidsai-docs/cucim/${VERSION_NUMBER}/txt"
+fi
diff --git a/ci/build_python.sh b/ci/build_python.sh
new file mode 100755
index 000000000..d12cd4208
--- /dev/null
+++ b/ci/build_python.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+# Copyright (c) 2022, NVIDIA CORPORATION.
+
+set -euo pipefail
+
+source rapids-env-update
+
+export CMAKE_GENERATOR=Ninja
+
+rapids-print-env
+
+rapids-logger "Begin py build"
+
+CPP_CHANNEL=$(rapids-download-conda-from-s3 cpp)
+
+# TODO: Remove `--no-test` flag once importing on a CPU
+# node works correctly
+rapids-mamba-retry mambabuild \
+  --no-test \
+  --channel "${CPP_CHANNEL}" \
+  conda/recipes/cucim
+
+rapids-upload-conda-to-s3 python
diff --git a/ci/check_style.sh b/ci/check_style.sh
new file mode 100755
index 000000000..be3ac3f4b
--- /dev/null
+++ b/ci/check_style.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+# Copyright (c) 2020-2022, NVIDIA CORPORATION.
+
+set -euo pipefail
+
+rapids-logger "Create checks conda environment"
+. /opt/conda/etc/profile.d/conda.sh
+
+rapids-dependency-file-generator \
+  --output conda \
+  --file_key checks \
+  --matrix "cuda=${RAPIDS_CUDA_VERSION%.*};arch=$(arch);py=${RAPIDS_PY_VERSION}" | tee env.yaml
+
+rapids-mamba-retry env create --force -f env.yaml -n checks
+conda activate checks
+
+# Run pre-commit checks
+pre-commit run --hook-stage manual --all-files --show-diff-on-failure
diff --git a/ci/cpu/build.sh b/ci/cpu/build.sh
index d86582d90..7e995d6f5 100755
--- a/ci/cpu/build.sh
+++ b/ci/cpu/build.sh
@@ -36,6 +36,10 @@ echo "PYTHON_VER  : ${PYTHON_VER}"
 export GPUCI_CONDA_RETRY_MAX=1
 export GPUCI_CONDA_RETRY_SLEEP=30
 
+# Workaround to keep Jenkins builds working
+# until we migrate fully to GitHub Actions
+export RAPIDS_CUDA_VERSION="${CUDA}"
+
 export CONDA_BLD_DIR="${WORKSPACE}/.conda-bld"
 
 ################################################################################
diff --git a/ci/gpu/build.sh b/ci/gpu/build.sh
index eb38e9640..ffcef97d4 100755
--- a/ci/gpu/build.sh
+++ b/ci/gpu/build.sh
@@ -38,6 +38,10 @@ echo "CUDA_VERSION: ${CUDA_VERSION}"
 echo "CUDA_VER    : ${CUDA_VER}"
 echo "PYTHON_VER  : ${PYTHON_VER}"
 
+# Workaround to keep Jenkins builds working
+# until we migrate fully to GitHub Actions
+export RAPIDS_CUDA_VERSION="${CUDA}"
+
 ################################################################################
 # SETUP - Check environment
 ################################################################################
diff --git a/ci/release/update-version.sh b/ci/release/update-version.sh
index 771d9a0a8..d05d36b99 100755
--- a/ci/release/update-version.sh
+++ b/ci/release/update-version.sh
@@ -44,3 +44,7 @@ sed_runner "s#cucim.kit.cuslide@${CURRENT_LONG_TAG}.so#cucim.kit.cuslide@${NEXT_
 sed_runner "s#cucim.kit.cuslide@${CURRENT_LONG_TAG}.so#cucim.kit.cuslide@${NEXT_FULL_TAG}.so#g" cucim.code-workspace
 sed_runner "s#branch-${CURRENT_MAJOR}.${CURRENT_MINOR}#branch-${NEXT_MAJOR}.${NEXT_MINOR}#g" README.md
 sed_runner "s#branch-${CURRENT_MAJOR}.${CURRENT_MINOR}#branch-${NEXT_MAJOR}.${NEXT_MINOR}#g" python/cucim/README.md
+
+for FILE in .github/workflows/*.yaml; do
+  sed_runner "/shared-action-workflows/ s/@.*/@branch-${NEXT_SHORT_TAG}/g" "${FILE}"
+done
diff --git a/conda/environments/all_cuda-118_arch-x86_64.yaml b/conda/environments/all_cuda-118_arch-x86_64.yaml
new file mode 100644
index 000000000..5d46021d7
--- /dev/null
+++ b/conda/environments/all_cuda-118_arch-x86_64.yaml
@@ -0,0 +1,49 @@
+# This file is generated by `rapids-dependency-file-generator`.
+# To make changes, edit ../../dependencies.yaml and run `rapids-dependency-file-generator`.
+channels:
+- rapidsai
+- rapidsai-nightly
+- conda-forge
+- nvidia
+dependencies:
+- GPUtil>=1.4.0
+- c-compiler
+- click
+- cmake>=3.23.1,!=3.25.0
+- cudatoolkit=11.8
+- cupy >=10,<12.0.0a0
+- cxx-compiler
+- gcc_linux-64=9.*
+- imagecodecs>=2021.6.8
+- ipython
+- jbig
+- jpeg
+- libwebp-base
+- nbsphinx
+- ninja
+- numpy >=1.21.3
+- numpydoc
+- nvcc_linux-64=11.8
+- openslide-python>=1.1.2
+- pip
+- pre-commit
+- psutil>=5.8.0
+- pydata-sphinx-theme
+- pytest-cov>=2.12.1
+- pytest-lazy-fixture>=0.6.3
+- pytest-xdist
+- pytest>=6.2.4
+- python>=3.8,<3.11
+- recommonmark
+- scikit-image >=0.19.0,<0.20.0a0
+- scipy
+- sphinx<6
+- sysroot_linux-64==2.17
+- tifffile>=2022.7.28
+- xz
+- yasm
+- zlib
+- zstd
+- pip:
+  - opencv-python-headless>=4.6
+name: all_cuda-118_arch-x86_64
diff --git a/conda/environments/env.yml b/conda/environments/env.yml
deleted file mode 100644
index 0cee99491..000000000
--- a/conda/environments/env.yml
+++ /dev/null
@@ -1,28 +0,0 @@
-name: cucim
-channels:
-  - conda-forge
-dependencies:
-  - cupy>=9
-  - scikit-image>=0.19.0
-  - openslide
-  - zlib
-  - jpeg
-  - jbig
-  - xz
-  - zstd
-  - libwebp-base  # [linux or osx]
-  - numpy
-  - xorg-libxcb
-  - scipy
-  - python=3.8
-  - cudatoolkit=11.0
-  - cudatoolkit-dev=11.0
-  - cxx-compiler
-  - c-compiler
-  - cmake
-  - automake
-  - make
-  - click
-  - yasm # [x86_64]
-  - isort
-  - flake8
diff --git a/conda/recipes/cucim/meta.yaml b/conda/recipes/cucim/meta.yaml
index fd691d636..252e705f6 100644
--- a/conda/recipes/cucim/meta.yaml
+++ b/conda/recipes/cucim/meta.yaml
@@ -1,11 +1,10 @@
-# Copyright (c) 2021, NVIDIA CORPORATION.
+# Copyright (c) 2021-2022, NVIDIA CORPORATION.
 
 {% set version = environ.get('GIT_DESCRIBE_TAG', '0.0.0.dev').lstrip('v') + environ.get('VERSION_SUFFIX', '') %}
-{% set minor_version =  version.split('.')[0] + '.' + version.split('.')[2] %}
-{% set py_version=environ.get('CONDA_PY', 37) %}
-{% set python_version=environ.get('PYTHON_VER', '3.7') %}
-{% set cuda_version='.'.join(environ.get('CUDA', '11.0').split('.')[:2]) %}
-{% set cuda_major=cuda_version.split('.')[0] %}
+{% set minor_version = version.split('.')[0] + '.' + version.split('.')[2] %}
+{% set py_version = environ['CONDA_PY'] %}
+{% set cuda_version = '.'.join(environ['RAPIDS_CUDA_VERSION'].split('.')[:2]) %}
+{% set cuda_major = cuda_version.split('.')[0] %}
 
 package:
   name: cucim
@@ -17,36 +16,52 @@ source:
 build:
   number: {{ GIT_DESCRIBE_NUMBER }}
   string: cuda_{{ cuda_major }}_py{{ py_version }}_{{ GIT_DESCRIBE_HASH }}_{{ GIT_DESCRIBE_NUMBER }}
+  script_env:
+    - AWS_ACCESS_KEY_ID
+    - AWS_SECRET_ACCESS_KEY
+    - CMAKE_C_COMPILER_LAUNCHER
+    - CMAKE_CUDA_COMPILER_LAUNCHER
+    - CMAKE_CXX_COMPILER_LAUNCHER
+    - CMAKE_GENERATOR
+    - PARALLEL_LEVEL
+    - SCCACHE_BUCKET
+    - SCCACHE_IDLE_TIMEOUT
+    - SCCACHE_REGION
+    - SCCACHE_S3_KEY_PREFIX=cucim-aarch64 # [aarch64]
+    - SCCACHE_S3_KEY_PREFIX=cucim-linux64 # [linux64]
 
 requirements:
   build:
-    - cmake>=3.23.1,!=3.25.0
     - {{ compiler("c") }}
     - {{ compiler("cxx") }}
+    - cmake >=3.23.1,!=3.25.0
+    - make
+    - ninja
     - sysroot_{{ target_platform }} {{ sysroot_version }}
   host:
-    - cudatoolkit {{ cuda_version }}.*
-    - python {{ python_version }}.*
-    - libcucim {{ version }}.*
     - click
-    - cupy >=9,<12.0.0a0
-    - numpy 1.19
-    - scipy
+    - cudatoolkit ={{ cuda_version }}
+    - cupy >=10,<12.0.0a0
+    - libcucim ={{ version }}
+    - numpy 1.21
+    - python
     - scikit-image >=0.19.0,<0.20.0a0
+    - scipy
   run:
     - {{ pin_compatible('cudatoolkit', max_pin='x', min_pin='x') }}
-    - python {{ python_version }}.*
-    - libcucim {{ version }}.*
-    - click
-    - cupy >=9,<12.0.0a0
     - {{ pin_compatible('numpy') }}
-    - scipy
+    - click
+    - cupy >=10,<12.0.0a0
+    - libcucim ={{ version }}
+    # - openslide # skipping here but benchmark binary would need openslide library
+    - python
     - scikit-image >=0.19.0,<0.20.0a0
-    # - openslide # skipping here but benchmark binary would needs openslide library
+    - scipy
 
+# TODO: Remove the linux64 tags on tests after disabling gpuCI / Jenkins
 tests:                                 # [linux64]
   requirements:                        # [linux64]
-    - cudatoolkit {{ cuda_version }}.* # [linux64]
+    - cudatoolkit ={{ cuda_version }}  # [linux64]
   imports:                             # [linux64]
     - cucim                            # [linux64]
 
diff --git a/conda/recipes/libcucim/meta.yaml b/conda/recipes/libcucim/meta.yaml
index ae4c1d3c8..92a33b5b1 100644
--- a/conda/recipes/libcucim/meta.yaml
+++ b/conda/recipes/libcucim/meta.yaml
@@ -1,10 +1,9 @@
-# Copyright (c) 2021, NVIDIA CORPORATION.
+# Copyright (c) 2021-2022, NVIDIA CORPORATION.
 
 {% set version = environ.get('GIT_DESCRIBE_TAG', '0.0.0.dev').lstrip('v') + environ.get('VERSION_SUFFIX', '') %}
-{% set minor_version =  version.split('.')[0] + '.' + version.split('.')[1] %}
-{% set python_version=environ.get('PYTHON_VER', '3.7') %}
-{% set cuda_version='.'.join(environ.get('CUDA', '11.0').split('.')[:2]) %}
-{% set cuda_major=cuda_version.split('.')[0] %}
+{% set minor_version = version.split('.')[0] + '.' + version.split('.')[1] %}
+{% set cuda_version = '.'.join(environ['RAPIDS_CUDA_VERSION'].split('.')[:2]) %}
+{% set cuda_major = cuda_version.split('.')[0] %}
 
 package:
   name: libcucim
@@ -16,32 +15,48 @@ source:
 build:
   number: {{ GIT_DESCRIBE_NUMBER }}
   string: cuda{{ cuda_major }}_{{ GIT_DESCRIBE_HASH }}_{{ GIT_DESCRIBE_NUMBER }}
+  script_env:
+    - AWS_ACCESS_KEY_ID
+    - AWS_SECRET_ACCESS_KEY
+    - CMAKE_C_COMPILER_LAUNCHER
+    - CMAKE_CUDA_COMPILER_LAUNCHER
+    - CMAKE_CXX_COMPILER_LAUNCHER
+    - CMAKE_GENERATOR
+    - PARALLEL_LEVEL
+    - SCCACHE_BUCKET
+    - SCCACHE_IDLE_TIMEOUT
+    - SCCACHE_REGION
+    - SCCACHE_S3_KEY_PREFIX=libcucim-aarch64 # [aarch64]
+    - SCCACHE_S3_KEY_PREFIX=libcucim-linux64 # [linux64]
 
 requirements:
   build:
-    - cmake>=3.23.1,!=3.25.0
     - {{ compiler("c") }}
     - {{ compiler("cxx") }}
+    - binutils
+    - cmake >=3.23.1,!=3.25.0
+    - make
+    - ninja
     - sysroot_{{ target_platform }} {{ sysroot_version }}
     - yasm # [x86_64]
   host:
-    - cudatoolkit {{ cuda_version }}.*
-    - openslide
-    - zlib
-    - jpeg
+    - cudatoolkit ={{ cuda_version }}
     - jbig
+    - jpeg
+    - libwebp-base  # [linux or osx]
+    - openslide
     - xz
+    - zlib
     - zstd
-    - libwebp-base  # [linux or osx]
   run:
     - {{ pin_compatible('cudatoolkit', max_pin='x', min_pin='x') }}
-    # - openslide # skipping here but benchmark binary would needs openslide library
-    - zlib
-    - jpeg
     - jbig
+    - jpeg
+    - libwebp-base  # [linux or osx]
+    # - openslide # skipping here but benchmark binary would need openslide library
     - xz
+    - zlib
     - zstd
-    - libwebp-base  # [linux or osx]
 
 about:
   home: https://developer.nvidia.com/multidimensional-image-processing
diff --git a/cpp/plugins/cucim.kit.cumed/VERSION b/cpp/plugins/cucim.kit.cumed/VERSION
index ed6a80946..7c8431de5 100644
--- a/cpp/plugins/cucim.kit.cumed/VERSION
+++ b/cpp/plugins/cucim.kit.cumed/VERSION
@@ -1,1 +1,1 @@
-22.12.00
+23.02.00
diff --git a/cpp/plugins/cucim.kit.cuslide/VERSION b/cpp/plugins/cucim.kit.cuslide/VERSION
index ed6a80946..7c8431de5 100644
--- a/cpp/plugins/cucim.kit.cuslide/VERSION
+++ b/cpp/plugins/cucim.kit.cuslide/VERSION
@@ -1,1 +1,1 @@
-22.12.00
+23.02.00
diff --git a/cucim.code-workspace b/cucim.code-workspace
index cbd28718d..d53fa8936 100644
--- a/cucim.code-workspace
+++ b/cucim.code-workspace
@@ -33,7 +33,7 @@
 					"CUCIM_TESTDATA_FOLDER": "${workspaceDirectory}/test_data",
 					// Add cuslide plugin's library path to LD_LIBRARY_PATH
 					"LD_LIBRARY_PATH": "${workspaceDirectory}/build-debug/lib:${workspaceDirectory}/cpp/plugins/cucim.kit.cuslide/build-debug/lib:${workspaceDirectory}/temp/cuda/lib64:${os_env:LD_LIBRARY_PATH}",
-					"CUCIM_TEST_PLUGIN_PATH": "cucim.kit.cuslide@22.12.00.so"
+					"CUCIM_TEST_PLUGIN_PATH": "cucim.kit.cuslide@23.02.00.so"
 				},
 				"cwd": "${workspaceDirectory}",
 				"catch2": {
@@ -226,7 +226,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@22.12.00.so"
+						"value": "cucim.kit.cuslide@23.02.00.so"
 					}
 				],
 				"console": "externalTerminal",
@@ -254,7 +254,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@22.12.00.so"
+						"value": "cucim.kit.cuslide@23.02.00.so"
 					}
 				],
 				"console": "externalTerminal",
@@ -286,7 +286,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@22.12.00.so"
+						"value": "cucim.kit.cuslide@23.02.00.so"
 					}
 				],
 				"console": "externalTerminal",
diff --git a/dependencies.yaml b/dependencies.yaml
new file mode 100644
index 000000000..1fbff3985
--- /dev/null
+++ b/dependencies.yaml
@@ -0,0 +1,169 @@
+# Dependency list for https://github.com/rapidsai/dependency-file-generator
+files:
+  all:
+    output: conda
+    matrix:
+      cuda: ["11.8"]
+      arch: [x86_64]
+    includes:
+      - build
+      - cudatoolkit
+      - checks
+      - docs
+      - py_version
+      - run
+      - test_python
+  test_python:
+    output: none
+    includes:
+      - cudatoolkit
+      - py_version
+      - test_python
+  checks:
+    output: none
+    includes:
+      - checks
+      - py_version
+  docs:
+    output: none
+    includes:
+      - cudatoolkit
+      - docs
+      - py_version
+channels:
+  - rapidsai
+  - rapidsai-nightly
+  - conda-forge
+  - nvidia
+dependencies:
+  build:
+    common:
+      - output_types: [conda, requirements]
+        packages:
+          - cmake>=3.23.1,!=3.25.0
+          - ninja
+      - output_types: conda
+        packages:
+          - c-compiler
+          - cxx-compiler
+    specific:
+      - output_types: conda
+        matrices:
+          - matrix:
+              arch: x86_64
+            packages:
+              - gcc_linux-64=9.*
+              - sysroot_linux-64==2.17
+              - yasm
+          - matrix:
+              arch: aarch64
+            packages:
+              - gcc_linux-aarch64=9.*
+              - sysroot_linux-aarch64==2.17
+      - output_types: conda
+        matrices:
+          - matrix:
+              arch: x86_64
+              cuda: "11.8"
+            packages:
+              - nvcc_linux-64=11.8
+          - matrix:
+              arch: aarch64
+              cuda: "11.8"
+            packages:
+              - nvcc_linux-aarch64=11.8
+  checks:
+    common:
+      - output_types: [conda, requirements]
+        packages:
+          - pre-commit
+  cudatoolkit:
+    specific:
+      - output_types: conda
+        matrices:
+          - matrix:
+              cuda: "11.2"
+            packages:
+              - cudatoolkit=11.2
+          - matrix:
+              cuda: "11.4"
+            packages:
+              - cudatoolkit=11.4
+          - matrix:
+              cuda: "11.5"
+            packages:
+              - cudatoolkit=11.5
+          - matrix:
+              cuda: "11.8"
+            packages:
+              - cudatoolkit=11.8
+  docs:
+    common:
+      - output_types: [conda, requirements]
+        packages:
+          - ipython
+          - nbsphinx
+          - numpydoc
+          - pydata-sphinx-theme
+          - recommonmark
+          # TODO: sphinx upper version limit can likely be removed after pydata 0.13
+          # is released. See https://github.com/pydata/pydata-sphinx-theme/issues/1094
+          - sphinx<6
+  py_version:
+    specific:
+      - output_types: conda
+        matrices:
+          - matrix:
+              py: "3.8"
+            packages:
+              - python=3.8
+          - matrix:
+              py: "3.9"
+            packages:
+              - python=3.9
+          - matrix:
+              py: "3.10"
+            packages:
+              - python=3.10
+          - matrix:
+            packages:
+              - python>=3.8,<3.11
+  run:
+    common:
+      - output_types: conda
+        packages:
+          - click
+          - cupy >=10,<12.0.0a0
+          - jbig
+          - jpeg
+          - libwebp-base
+          - numpy >=1.21.3
+          - scikit-image >=0.19.0,<0.20.0a0
+          - scipy
+          - xz
+          - zlib
+          - zstd
+          # Not sure where these go, if anywhere:
+          # - openslide
+          # - xorg-libxcb
+  test_python:
+    common:
+      - output_types: [conda, requirements]
+        packages:
+          - GPUtil>=1.4.0
+          - imagecodecs>=2021.6.8
+          - openslide-python>=1.1.2
+          - psutil>=5.8.0
+          - pytest-cov>=2.12.1
+          - pytest-lazy-fixture>=0.6.3
+          - pytest-xdist
+          - pytest>=6.2.4
+          - tifffile>=2022.7.28
+      - output_types: [conda]
+        packages:
+          - pip
+          - pip:
+            - opencv-python-headless>=4.6
+      - output_types: [requirements]
+        packages:
+          - opencv-python-headless>=4.6
diff --git a/docs/source/conf.py b/docs/source/conf.py
index fbe3fb514..1c057a067 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -72,9 +72,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '22.12'
+version = '23.02'
 # The full version, including alpha/beta/rc tags.
-release = '22.12.00'
+release = '23.02.00'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -147,7 +147,7 @@
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
     (master_doc, 'cucim.tex', 'cucim Documentation',
-     'Continuum Analytics', 'manual'),
+     'NVIDIA Corporation', 'manual'),
 ]
 
 
diff --git a/python/cucim/.pre-commit-config.yaml b/python/cucim/.pre-commit-config.yaml
deleted file mode 100644
index 6e974cba9..000000000
--- a/python/cucim/.pre-commit-config.yaml
+++ /dev/null
@@ -1,20 +0,0 @@
-# To install the git pre-commit hook run:
-#   pre-commit install
-# To update the pre-commit hooks run:
-#   pre-commit install-hooks
-exclude: '^(\.tox|ci/templates|\.bumpversion\.cfg)(/|$)'
-repos:
-  - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: master
-    hooks:
-      - id: trailing-whitespace
-      - id: end-of-file-fixer
-      - id: debug-statements
-  - repo: https://github.com/timothycrosley/isort
-    rev: master
-    hooks:
-      - id: isort
-  - repo: https://gitlab.com/pycqa/flake8
-    rev: master
-    hooks:
-      - id: flake8
diff --git a/python/cucim/MANIFEST.in b/python/cucim/MANIFEST.in
index 462e3037f..4244cebec 100644
--- a/python/cucim/MANIFEST.in
+++ b/python/cucim/MANIFEST.in
@@ -17,7 +17,7 @@ include ../../LICENSE
 include ../../LICENSE-3rdparty.md
 include README.md
 
-include tox.ini .travis.yml .appveyor.yml .readthedocs.yml .pre-commit-config.yaml
+include tox.ini .travis.yml .appveyor.yml .readthedocs.yml
 
 include versioneer.py
 include src/cucim/clara/*.so*
diff --git a/python/cucim/README.md b/python/cucim/README.md
deleted file mode 120000
index fe8400541..000000000
--- a/python/cucim/README.md
+++ /dev/null
@@ -1,1 +0,0 @@
-../../README.md
\ No newline at end of file
diff --git a/python/cucim/README.md b/python/cucim/README.md
new file mode 100755
index 000000000..9a92b8de5
--- /dev/null
+++ b/python/cucim/README.md
@@ -0,0 +1,113 @@
+# <div align="left"><img src="https://rapids.ai/assets/images/rapids_logo.png" width="90px"/>&nbsp;cuCIM</div>
+
+[RAPIDS](https://rapids.ai) cuCIM is an open-source, accelerated computer vision and image processing software library for multidimensional images used in biomedical, geospatial, material and life science, and remote sensing use cases.
+
+cuCIM offers:
+
+- Enhanced Image Processing Capabilities for large and n-dimensional tag image file format (TIFF) files
+- Accelerated performance through Graphics Processing Unit (GPU)-based image processing and computer vision primitives
+- A Straightforward Pythonic Interface with Matching Application Programming Interface (API) for Openslide
+
+cuCIM supports the following formats:
+
+- Aperio ScanScope Virtual Slide (SVS)
+- Philips TIFF
+- Generic Tiled, Multi-resolution RGB TIFF files with the following compression schemes:
+  - No Compression
+  - JPEG
+  - JPEG2000
+  - Lempel-Ziv-Welch (LZW)
+  - Deflate
+
+**NOTE:** For the latest stable [README.md](https://github.com/rapidsai/cucim/blob/main/README.md) ensure you are on the `main` branch.
+
+- [GTC 2022 Accelerating Storage IO to GPUs with Magnum IO [S41347]](https://events.rainfocus.com/widget/nvidia/gtcspring2022/sessioncatalog/session/1634960000577001Etxp)
+  - cuCIM's GDS API examples: <https://github.com/NVIDIA/MagnumIO/tree/main/gds/readers/cucim-gds>
+- [SciPy 2021 cuCIM - A GPU image I/O and processing library](https://www.scipy2021.scipy.org/)
+  - [video](https://youtu.be/G46kOOM9xbQ)
+- [GTC 2021 cuCIM: A GPU Image I/O and Processing Toolkit [S32194]](https://www.nvidia.com/en-us/on-demand/search/?facet.mimetype[]=event%20session&layout=list&page=1&q=cucim&sort=date)
+  - [video](https://www.nvidia.com/en-us/on-demand/session/gtcspring21-s32194/)
+
+**[Developer Page](https://developer.nvidia.com/multidimensional-image-processing)**
+
+**Blogs**
+- [Enhanced Image Analysis with Multidimensional Image Processing](https://developer.nvidia.com/blog/enhanced-image-analysis-with-multidimensional-image-processing/)
+- [Accelerating Scikit-Image API with cuCIM: n-Dimensional Image Processing and IO on GPUs](https://developer.nvidia.com/blog/cucim-rapid-n-dimensional-image-processing-and-i-o-on-gpus/)
+- [Accelerating Digital Pathology Pipelines with NVIDIA Clara™ Deploy](https://developer.nvidia.com/blog/accelerating-digital-pathology-pipelines-with-nvidia-clara-deploy-2/)
+
+**Webinars**
+
+- [cuCIM: a GPU Image IO and Processing Library](https://www.youtube.com/watch?v=G46kOOM9xbQ)
+
+**[Documentation](https://docs.rapids.ai/api/cucim/stable)**
+
+**Release notes** are available on our [wiki page](https://github.com/rapidsai/cucim/wiki/Release-Notes).
+
+## Install cuCIM
+
+### Conda
+
+#### [Conda (stable)](https://anaconda.org/rapidsai/cucim)
+
+> conda create -n cucim -c rapidsai -c conda-forge cucim cudatoolkit=`<CUDA version>`
+
+`<CUDA version>` should be 11.0+ (e.g., `11.0`, `11.2`, etc.)
+
+#### [Conda (nightlies)](https://anaconda.org/rapidsai-nightly/cucim)
+
+> conda create -n cucim -c rapidsai-nightly -c conda-forge cucim cudatoolkit=`<CUDA version>`
+
+`<CUDA version>` should be 11.0+ (e.g., `11.0`, `11.2`, etc)
+
+### [PyPI](https://pypi.org/project/cucim/)
+
+```bash
+pip install cucim
+
+# Install dependencies for `cucim.skimage` (assuming that CUDA 11.0 is used for CuPy)
+pip install scipy scikit-image cupy-cuda110
+```
+
+### Notebooks
+
+Please check out our [Welcome](notebooks/Welcome.ipynb) notebook ([NBViewer](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-23.02/notebooks/Welcome.ipynb))
+
+#### Downloading sample images
+
+To download images used in the notebooks, please execute the following commands from the repository root folder to copy sample input images into `notebooks/input` folder:
+
+(You will need [Docker](https://www.docker.com/) installed in your system)
+
+```bash
+./run download_testdata
+```
+or
+
+```bash
+mkdir -p notebooks/input
+tmp_id=$(docker create gigony/svs-testdata:little-big)
+docker cp $tmp_id:/input notebooks
+docker rm -v ${tmp_id}
+```
+
+## Build/Install from Source
+
+See build [instructions](CONTRIBUTING.md#setting-up-your-build-environment).
+
+## Contributing Guide
+
+Contributions to cuCIM are more than welcome!
+Please review the [CONTRIBUTING.md](https://github.com/rapidsai/cucim/blob/main/CONTRIBUTING.md) file for information on how to contribute code and issues to the project.
+
+## Acknowledgments
+
+Without awesome third-party open source software, this project wouldn't exist.
+
+Please find [LICENSE-3rdparty.md](LICENSE-3rdparty.md) to see which third-party open source software
+is used in this project.
+
+## License
+
+Apache-2.0 License (see [LICENSE](LICENSE) file).
+
+Copyright (c) 2020-2022, NVIDIA CORPORATION.
diff --git a/python/cucim/VERSION b/python/cucim/VERSION
index ed6a80946..7c8431de5 100644
--- a/python/cucim/VERSION
+++ b/python/cucim/VERSION
@@ -1,1 +1,1 @@
-22.12.00
+23.02.00
diff --git a/python/cucim/docs/getting_started/index.md b/python/cucim/docs/getting_started/index.md
index 941d719b7..f03db8a11 100644
--- a/python/cucim/docs/getting_started/index.md
+++ b/python/cucim/docs/getting_started/index.md
@@ -14,15 +14,15 @@
 
 ## Installation
 
-Please download the latest SDK package (`cuCIM-v22.12.00-linux.tar.gz`).
+Please download the latest SDK package (`cuCIM-v23.02.00-linux.tar.gz`).
 
 Untar the downloaded file.
 
 ```bash
-mkdir -p cuCIM-v22.12.00
-tar -xzvf cuCIM-v22.12.00-linux.tar.gz -C cuCIM-v22.12.00
+mkdir -p cuCIM-v23.02.00
+tar -xzvf cuCIM-v23.02.00-linux.tar.gz -C cuCIM-v23.02.00
 
-cd cuCIM-v22.12.00
+cd cuCIM-v23.02.00
 ```
 
 ## Run command
@@ -147,7 +147,7 @@ Its execution would show some metadata information and create two files -- `outp
 ```
 $ ./bin/tiff_image notebooks/input/image.tif .
 [Plugin: cucim.kit.cuslide] Loading...
-[Plugin: cucim.kit.cuslide] Loading the dynamic library from: cucim.kit.cuslide@22.12.00.so
+[Plugin: cucim.kit.cuslide] Loading the dynamic library from: cucim.kit.cuslide@23.02.00.so
 [Plugin: cucim.kit.cuslide] loaded successfully. Version: 0
 Initializing plugin: cucim.kit.cuslide (interfaces: [cucim::io::IImageFormat v0.1]) (impl: cucim.kit.cuslide)
 is_loaded: true
diff --git a/python/cucim/docs/index.md b/python/cucim/docs/index.md
index 614784bc3..6fc64c1cd 100644
--- a/python/cucim/docs/index.md
+++ b/python/cucim/docs/index.md
@@ -18,7 +18,7 @@ development/index -->
 
 # cuCIM Documentation
 
-Current latest version is [Version 22.12.00](release_notes/v22.12.00.md).
+Current latest version is [Version 23.02.00](release_notes/v23.02.00.md).
 
 **cuCIM** a toolkit to provide GPU accelerated I/O, image  processing & computer vision primitives for N-Dimensional images with a focus on biomedical imaging.
 
diff --git a/python/cucim/setup.py b/python/cucim/setup.py
index 57e390ad9..1b7179b6b 100755
--- a/python/cucim/setup.py
+++ b/python/cucim/setup.py
@@ -60,6 +60,7 @@ def read(*names, **kwargs):
         'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',
         'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
         # 'Operating System :: OS Independent',
         # 'Operating System :: Unix',
         # 'Operating System :: POSIX',
diff --git a/python/cucim/src/cucim/core/operations/color/jitter.py b/python/cucim/src/cucim/core/operations/color/jitter.py
index 6e36d8369..fbae3dcbc 100755
--- a/python/cucim/src/cucim/core/operations/color/jitter.py
+++ b/python/cucim/src/cucim/core/operations/color/jitter.py
@@ -161,7 +161,7 @@ def _adjust_saturation(input_arr, saturation):
 
 # hue jitter
 def _adjust_hue(input_arr, hue):
-    if not(-0.5 <= hue <= 0.5):
+    if not (-0.5 <= hue <= 0.5):
         raise ValueError('hue factor({}) is not in [-0.5, 0.5].'.
                          format(hue))
 
diff --git a/python/cucim/src/cucim/core/operations/color/kernel/cuda_kernel_source.py b/python/cucim/src/cucim/core/operations/color/kernel/cuda_kernel_source.py
index 563b53e4c..ac60de508 100644
--- a/python/cucim/src/cucim/core/operations/color/kernel/cuda_kernel_source.py
+++ b/python/cucim/src/cucim/core/operations/color/kernel/cuda_kernel_source.py
@@ -205,7 +205,7 @@
     float s = __uint2float_rn((unsigned int)s_char);
     float v = __uint2float_rn((unsigned int)v_char);
 
-    if (s == 0) {
+    if (v == 0) {
       // write zero and out
       output_rgb[idx] = 0;
       output_rgb[idx+pitch] = 0;
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 8a2eb72f4..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,12 +4,13 @@
 import cupy
 import numpy as np
 
-from ._pba_2d import (_check_distances, _check_indices,
-                      _distance_tranform_arg_check, _get_block_size,
-                      _generate_shape, _generate_indices_ops,
-                      lcm)
+from cucim.skimage._vendored import pad
 from cucim.skimage._vendored._ndimage_util import _get_inttype
 
+from ._pba_2d import (_check_distances, _check_indices,
+                      _distance_tranform_arg_check, _generate_indices_ops,
+                      _generate_shape, _get_block_size, lcm)
+
 pba3d_defines_template = """
 
 #define MARKER     {marker}
@@ -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/_shared/utils.py b/python/cucim/src/cucim/skimage/_shared/utils.py
index cb1b2ce3f..b4371ff1b 100644
--- a/python/cucim/src/cucim/skimage/_shared/utils.py
+++ b/python/cucim/src/cucim/skimage/_shared/utils.py
@@ -7,11 +7,7 @@
 import cupy as cp
 import numpy as np
 
-from ._warnings import all_warnings, warn
-
-__all__ = ['deprecated', 'get_bound_method_class', 'all_warnings',
-           'safe_as_int', 'check_shape_equality', 'check_nD', 'warn',
-           'reshape_nd', 'identity', 'slice_at_axis']
+from ._warnings import all_warnings, warn  # noqa
 
 
 class skimage_deprecation(Warning):
@@ -444,10 +440,11 @@ def safe_as_int(val, atol=1e-3):
     return np.round(val).astype(np.int64)
 
 
-def check_shape_equality(im1, im2):
-    """Raise an error if the shape do not match."""
-    if not im1.shape == im2.shape:
-        raise ValueError("Input images must have the same dimensions.")
+def check_shape_equality(*images):
+    """Check that all images have the same shape"""
+    image0 = images[0]
+    if not all(image0.shape == image.shape for image in images[1:]):
+        raise ValueError('Input images must have the same dimensions.')
     return
 
 
@@ -721,3 +718,23 @@ def _supported_float_type(input_dtype, allow_complex=False):
 def identity(image, *args, **kwargs):
     """Returns the first argument unmodified."""
     return image
+
+
+def as_binary_ndarray(array, *, variable_name):
+    """Return `array` as a numpy.ndarray of dtype bool.
+
+    Raises
+    ------
+    ValueError:
+        An error including the given `variable_name` if `array` can not be
+        safely cast to a boolean array.
+    """
+    array = cp.asarray(array)
+    if array.dtype != bool:
+        if cp.any((array != 1) & (array != 0)):
+            raise ValueError(
+                f"{variable_name} array is not of dtype boolean or "
+                f"contains values other than 0 and 1 so cannot be "
+                f"safely cast to boolean array."
+            )
+    return cp.asarray(array, dtype=bool)
diff --git a/python/cucim/src/cucim/skimage/_vendored/__init__.py b/python/cucim/src/cucim/skimage/_vendored/__init__.py
index 6c7d90791..8b74f933d 100644
--- a/python/cucim/src/cucim/skimage/_vendored/__init__.py
+++ b/python/cucim/src/cucim/skimage/_vendored/__init__.py
@@ -3,4 +3,6 @@
 not yet been merged to CuPy master.
 """
 
+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_filters.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
index 60f173435..633371887 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters.py
@@ -848,7 +848,8 @@ def _min_or_max_filter(input, size, ftprnt, structure, output, mode, cval,
             input, sizes, output, mode, cval, origin)
 
     origins, int_type = _filters_core._check_nd_args(input, ftprnt,
-                                                     mode, origin, 'footprint')
+                                                     mode, origin, 'footprint',
+                                                     sizes=sizes)
     if structure is not None and structure.ndim != input.ndim:
         raise RuntimeError('structure array has incorrect shape')
 
@@ -1081,28 +1082,51 @@ def get_rank(fs):
 
 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)
+    sizes, footprint, _ = _filters_core._check_size_footprint_structure(
+        input.ndim, size, footprint, None, force_footprint=False)
     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:
+                                                     mode, origin, 'footprint',
+                                                     sizes=sizes)
+    has_weights = True
+    if sizes is not None:
+        has_weights = False
+        filter_size = internal.prod(sizes)
+        if filter_size == 0:
+            return cupy.zeros_like(input)
+        footprint_shape = tuple(sizes)
+    elif footprint.size == 0:
         return cupy.zeros_like(input)
-    filter_size = int(footprint.sum())
+    else:
+        footprint_shape = footprint.shape
+        filter_size = int(footprint.sum())
+        if filter_size == footprint.size:
+            # can omit passing the footprint if it is all ones
+            sizes = footprint.shape
+            has_weights = False
+
     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,
+        min_max_op = 'min'
+    elif rank == filter_size - 1:
+        min_max_op = 'max'
+    else:
+        min_max_op = None
+    if min_max_op is not None:
+        if sizes is not None:
+            return _min_or_max_filter(input, sizes[0], None, None, output,
+                                      mode, cval, origins, min_max_op)
+        else:
+            return _min_or_max_filter(input, None, footprint, None, output,
+                                      mode, cval, origins, min_max_op)
+    offsets = _filters_core._origins_to_offsets(origins, footprint_shape)
+    kernel = _get_rank_kernel(filter_size, rank, mode, footprint_shape,
+                              offsets, float(cval), int_type,
+                              has_weights=has_weights)
+    return _filters_core._call_kernel(kernel, input, None, output,
                                       weights_dtype=bool)
 
 
@@ -1134,7 +1158,7 @@ def _get_shell_gap(filter_size):
 
 @cupy._util.memoize(for_each_device=True)
 def _get_rank_kernel(filter_size, rank, mode, w_shape, offsets, cval,
-                     int_type):
+                     int_type, has_weights):
     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,
@@ -1184,4 +1208,5 @@ def _get_rank_kernel(filter_size, rank, mode, w_shape, offsets, cval,
         '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)
+        mode, w_shape, int_type, offsets, cval, has_weights=has_weights,
+        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 469dff9d9..a09e22306 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_filters_core.py
@@ -56,14 +56,24 @@ def _convert_1d_args(ndim, weights, origin, axis):
     return weights, tuple(origins)
 
 
-def _check_nd_args(input, weights, mode, origin, wghts_name='filter weights'):
+def _check_nd_args(input, weights, mode, origin, wghts_name='filter weights',
+                   sizes=None):
     _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))
+    if weights is not None:
+        # 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)
+            )
+    elif sizes is None:
+        raise ValueError("must specify either weights array or sizes")
+    else:
+        weight_dims = sizes
     origins = _util._fix_sequence_arg(origin, len(weight_dims), 'origin', int)
     for origin, width in zip(origins, weight_dims):
         _util._check_origin(origin, width)
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/_ndimage_util.py b/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
index 1a6a830a5..0b0b0be6d 100644
--- a/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
+++ b/python/cucim/src/cucim/skimage/_vendored/_ndimage_util.py
@@ -168,7 +168,7 @@ def _generate_boundary_condition_ops(mode, ix, xsize, int_t="int",
             {ix} %= {xsize};
             '''
             ops_lower = ops_upper + f'''
-            if ({ix} < 0) {{
+            while ({ix} < 0) {{
                 {ix} += {xsize};
             }}'''
             ops = (ops_lower, ops_upper)
diff --git a/python/cucim/src/cucim/skimage/_vendored/_pearsonr.py b/python/cucim/src/cucim/skimage/_vendored/_pearsonr.py
new file mode 100644
index 000000000..0da8541f7
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/_vendored/_pearsonr.py
@@ -0,0 +1,268 @@
+"""Implementation of cupyx.scipy.stats.pearsonr (currently missing in CuPy)
+
+Simple port of SciPy's pearsonr in scipy/stats/_stats_py.py.
+
+Note that this is based on the simpler implementation from SciPy<1.9 where the
+return type is just a 2-tuple of (statistic, p-value). In SciPy>=1.9, the
+return type changed to a PearsonRResult object that also has a
+`confidence_interval` method. Since we do not need that method in the cuCIM
+API, a simple tuple-based return is used
+"""
+import warnings
+
+import cupy as cp
+import numpy as np
+from scipy import special
+
+
+class PearsonRConstantInputWarning(RuntimeWarning):
+    """Warning generated by `pearsonr` when an input is constant."""
+
+    def __init__(self, msg=None):
+        if msg is None:
+            msg = ("An input array is constant; the correlation coefficient "
+                   "is not defined.")
+        self.args = (msg,)
+
+
+class PearsonRNearConstantInputWarning(RuntimeWarning):
+    """Warning generated by `pearsonr` when an input is nearly constant."""
+
+    def __init__(self, msg=None):
+        if msg is None:
+            msg = ("An input array is nearly constant; the computed "
+                   "correlation coefficient may be inaccurate.")
+        self.args = (msg,)
+
+
+# Note: adapted from scipy.stats._stats_py.pearsonr in SciPy 1.8.1
+def pearsonr(x, y, *, disable_checks=False):
+    r"""
+    Pearson correlation coefficient and p-value for testing non-correlation.
+
+    The Pearson correlation coefficient [1]_ measures the linear relationship
+    between two datasets.  The calculation of the p-value relies on the
+    assumption that each dataset is normally distributed.  (See Kowalski [3]_
+    for a discussion of the effects of non-normality of the input on the
+    distribution of the correlation coefficient.)  Like other correlation
+    coefficients, this one varies between -1 and +1 with 0 implying no
+    correlation. Correlations of -1 or +1 imply an exact linear relationship.
+
+    Parameters
+    ----------
+    x : (N,) array_like
+        Input array.
+    y : (N,) array_like
+        Input array.
+
+    Returns
+    -------
+    r : float
+        Pearson's correlation coefficient.
+    p-value : float
+        Two-tailed p-value.
+
+    Warns
+    -----
+    PearsonRConstantInputWarning
+        Raised if an input is a constant array.  The correlation coefficient
+        is not defined in this case, so ``cp.nan`` is returned.
+
+    PearsonRNearConstantInputWarning
+        Raised if an input is "nearly" constant.  The array ``x`` is considered
+        nearly constant if ``norm(x - mean(x)) < 1e-13 * abs(mean(x))``.
+        Numerical errors in the calculation ``x - mean(x)`` in this case might
+        result in an inaccurate calculation of r.
+
+    See Also
+    --------
+    spearmanr : Spearman rank-order correlation coefficient.
+    kendalltau : Kendall's tau, a correlation measure for ordinal data.
+
+    Notes
+    -----
+    The correlation coefficient is calculated as follows:
+
+    .. math::
+
+        r = \frac{\sum (x - m_x) (y - m_y)}
+                 {\sqrt{\sum (x - m_x)^2 \sum (y - m_y)^2}}
+
+    where :math:`m_x` is the mean of the vector x and :math:`m_y` is
+    the mean of the vector y.
+
+    Under the assumption that x and y are drawn from
+    independent normal distributions (so the population correlation coefficient
+    is 0), the probability density function of the sample correlation
+    coefficient r is ([1]_, [2]_):
+
+    .. math::
+
+        f(r) = \frac{{(1-r^2)}^{n/2-2}}{\mathrm{B}(\frac{1}{2},\frac{n}{2}-1)}
+
+    where n is the number of samples, and B is the beta function.  This
+    is sometimes referred to as the exact distribution of r.  This is
+    the distribution that is used in `pearsonr` to compute the p-value.
+    The distribution is a beta distribution on the interval [-1, 1],
+    with equal shape parameters a = b = n/2 - 1.  In terms of SciPy's
+    implementation of the beta distribution, the distribution of r is::
+
+        dist = scipy.stats.beta(n/2 - 1, n/2 - 1, loc=-1, scale=2)
+
+    The p-value returned by `pearsonr` is a two-sided p-value. The p-value
+    roughly indicates the probability of an uncorrelated system
+    producing datasets that have a Pearson correlation at least as extreme
+    as the one computed from these datasets. More precisely, for a
+    given sample with correlation coefficient r, the p-value is
+    the probability that abs(r') of a random sample x' and y' drawn from
+    the population with zero correlation would be greater than or equal
+    to abs(r). In terms of the object ``dist`` shown above, the p-value
+    for a given r and length n can be computed as::
+
+        p = 2*dist.cdf(-abs(r))
+
+    When n is 2, the above continuous distribution is not well-defined.
+    One can interpret the limit of the beta distribution as the shape
+    parameters a and b approach a = b = 0 as a discrete distribution with
+    equal probability masses at r = 1 and r = -1.  More directly, one
+    can observe that, given the data x = [x1, x2] and y = [y1, y2], and
+    assuming x1 != x2 and y1 != y2, the only possible values for r are 1
+    and -1.  Because abs(r') for any sample x' and y' with length 2 will
+    be 1, the two-sided p-value for a sample of length 2 is always 1.
+
+    References
+    ----------
+    .. [1] "Pearson correlation coefficient", Wikipedia,
+           https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
+    .. [2] Student, "Probable error of a correlation coefficient",
+           Biometrika, Volume 6, Issue 2-3, 1 September 1908, pp. 302-310.
+    .. [3] C. J. Kowalski, "On the Effects of Non-Normality on the Distribution
+           of the Sample Product-Moment Correlation Coefficient"
+           Journal of the Royal Statistical Society. Series C (Applied
+           Statistics), Vol. 21, No. 1 (1972), pp. 1-12.
+
+    Examples
+    --------
+    >>> from scipy import stats
+    >>> stats.pearsonr([1, 2, 3, 4, 5], [10, 9, 2.5, 6, 4])
+    (-0.7426106572325057, 0.1505558088534455)
+
+    There is a linear dependence between x and y if y = a + b*x + e, where
+    a,b are constants and e is a random error term, assumed to be independent
+    of x. For simplicity, assume that x is standard normal, a=0, b=1 and let
+    e follow a normal distribution with mean zero and standard deviation s>0.
+
+    >>> s = 0.5
+    >>> x = stats.norm.rvs(size=500)
+    >>> e = stats.norm.rvs(scale=s, size=500)
+    >>> y = x + e
+    >>> stats.pearsonr(x, y)
+    (0.9029601878969703, 8.428978827629898e-185) # may vary
+
+    This should be close to the exact value given by
+
+    >>> 1/cp.sqrt(1 + s**2)
+    0.8944271909999159
+
+    For s=0.5, we observe a high level of correlation. In general, a large
+    variance of the noise reduces the correlation, while the correlation
+    approaches one as the variance of the error goes to zero.
+
+    It is important to keep in mind that no correlation does not imply
+    independence unless (x, y) is jointly normal. Correlation can even be zero
+    when there is a very simple dependence structure: if X follows a
+    standard normal distribution, let y = abs(x). Note that the correlation
+    between x and y is zero. Indeed, since the expectation of x is zero,
+    cov(x, y) = E[x*y]. By definition, this equals E[x*abs(x)] which is zero
+    by symmetry. The following lines of code illustrate this observation:
+
+    >>> y = cp.abs(x)
+    >>> stats.pearsonr(x, y)
+    (-0.016172891856853524, 0.7182823678751942) # may vary
+
+    A non-zero correlation coefficient can be misleading. For example, if X has
+    a standard normal distribution, define y = x if x < 0 and y = 0 otherwise.
+    A simple calculation shows that corr(x, y) = sqrt(2/Pi) = 0.797...,
+    implying a high level of correlation:
+
+    >>> y = cp.where(x < 0, x, 0)
+    >>> stats.pearsonr(x, y)
+    (0.8537091583771509, 3.183461621422181e-143) # may vary
+
+    This is unintuitive since there is no dependence of x and y if x is larger
+    than zero which happens in about half of the cases if we sample x and y.
+    """
+    # inputs must be 1D
+    n = len(x)
+    if n != len(y):
+        raise ValueError('x and y must have the same length.')
+
+    if n < 2:
+        raise ValueError('x and y must have length at least 2.')
+
+    if not disable_checks:
+        # If an input is constant, the correlation coefficient is not defined.
+        if (x == x[0]).all() or (y == y[0]).all():
+            warnings.warn(PearsonRConstantInputWarning())
+            return cp.nan, cp.nan
+
+    # dtype is the data type for the calculations.  This expression ensures
+    # that the data type is at least 64 bit floating point.  It might have
+    # more precision if the input is, for example, cp.longdouble.
+    dtype = cp.result_type(x.dtype, y.dtype, float)
+
+    if n == 2:
+        # faster on host for such a small inputs
+        x = cp.asnumpy(x)
+        y = cp.asnumpy(y)
+        return dtype(np.sign(x[1] - x[0]) * np.sign(y[1] - y[0])), 1.0
+
+    xmean = x.mean(dtype=dtype)
+    ymean = y.mean(dtype=dtype)
+
+    # By using `astype(dtype)`, we ensure that the intermediate calculations
+    # use at least 64 bit floating point.
+    xm = x.astype(dtype) - xmean
+    ym = y.astype(dtype) - ymean
+
+    # TODO: use cupyx.scipy.linalg.norm from CuPy once available
+    # Unlike cp.linalg.norm or the expression sqrt((xm*xm).sum()),
+    # scipy.linalg.norm(xm) does not overflow if xm is, for example,
+    # [-5e210, 5e210, 3e200, -3e200]
+    normxm = cp.linalg.norm(xm)
+    normym = cp.linalg.norm(ym)
+
+    if not disable_checks:
+        threshold = 1e-13
+        if normxm < threshold * abs(xmean) or normym < threshold * abs(ymean):
+            # If all the values in x (likewise y) are very close to the mean,
+            # the loss of precision that occurs in the subtraction
+            # xm = x - xmean might result in large errors in r.
+            warnings.warn(PearsonRNearConstantInputWarning())
+
+    r = float(cp.dot(xm / normxm, ym / normym))
+
+    # Presumably, if abs(r) > 1, then it is only some small artifact of
+    # floating point arithmetic.
+    r = max(min(r, 1.0), -1.0)
+
+    # As explained in the docstring, the p-value can be computed as
+    #     p = 2*dist.cdf(-abs(r))
+    # where dist is the beta distribution on [-1, 1] with shape parameters
+    # a = b = n/2 - 1.  `special.btdtr` is the CDF for the beta distribution
+    # on [0, 1].  To use it, we make the transformation  x = (r + 1)/2; the
+    # shape parameters do not change.  Then -abs(r) used in `cdf(-abs(r))`
+    # becomes x = (-abs(r) + 1)/2 = 0.5*(1 - abs(r)).  (r is cast to float64
+    # to avoid a TypeError raised by btdtr when r is higher precision.)
+    ab = n / 2 - 1
+    # scalar valued, so use special.btdtr from SciPy, not CuPy
+    prob = 2 * special.btdtr(ab, ab, 0.5 * (1.0 - abs(r)))
+    if disable_checks:
+        # warn only based on output values to avoid overhead of host/device
+        # synchronization needed for the disabled checks above.
+        if np.isnan(r) or np.isnan(prob):
+            warnings.warn(
+                "NaN encountered during Pearson R calculation. This may occur "
+                "in same cases such as a nearly constant-valued input."
+            )
+    return r, prob
diff --git a/python/cucim/src/cucim/skimage/_vendored/ndimage.py b/python/cucim/src/cucim/skimage/_vendored/ndimage.py
index cd97578a8..31d8b20d0 100644
--- a/python/cucim/src/cucim/skimage/_vendored/ndimage.py
+++ b/python/cucim/src/cucim/skimage/_vendored/ndimage.py
@@ -1,82 +1,93 @@
 # locally defined filters that are more efficient than in CuPy
-from cucim.skimage._vendored._ndimage_filters import correlate  # NOQA
+# measurements
+# fourier filters
+# additional 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
+from cupyx.scipy.ndimage import generic_filter  # NOQA
+from cupyx.scipy.ndimage import generic_filter1d  # NOQA
+from cupyx.scipy.ndimage import label  # 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 correlate  # NOQA
+from cucim.skimage._vendored._ndimage_filters import correlate1d  # 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 gaussian_filter1d  # NOQA
+from cucim.skimage._vendored._ndimage_filters import \
+    gaussian_gradient_magnitude  # 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 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 minimum_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import minimum_filter1d  # NOQA
 from cucim.skimage._vendored._ndimage_filters import percentile_filter  # NOQA
-
+from cucim.skimage._vendored._ndimage_filters import prewitt  # NOQA
+from cucim.skimage._vendored._ndimage_filters import rank_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import sobel  # NOQA
+from cucim.skimage._vendored._ndimage_filters import uniform_filter  # NOQA
+from cucim.skimage._vendored._ndimage_filters import uniform_filter1d  # 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 \
+    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 \
+    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_closing  # NOQA
 from cucim.skimage._vendored._ndimage_morphology import binary_dilation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import binary_erosion  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import \
+    binary_fill_holes  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import \
+    binary_hit_or_miss  # 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 \
+    binary_propagation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import black_tophat  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import \
+    generate_binary_structure  # NOQA
 from cucim.skimage._vendored._ndimage_morphology import grey_closing  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_dilation  # NOQA
+from cucim.skimage._vendored._ndimage_morphology import grey_erosion  # 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 \
+    iterate_structure  # 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 extrema  # NOQA
 from cupyx.scipy.ndimage import histogram  # NOQA
 from cupyx.scipy.ndimage import labeled_comprehension  # NOQA
+from cupyx.scipy.ndimage import maximum  # NOQA
+from cupyx.scipy.ndimage import maximum_position  # NOQA
+from cupyx.scipy.ndimage import mean  # NOQA
+from cupyx.scipy.ndimage import median  # NOQA
+from cupyx.scipy.ndimage import minimum  # NOQA
+from cupyx.scipy.ndimage import minimum_position  # NOQA
+from cupyx.scipy.ndimage import standard_deviation  # NOQA
+from cupyx.scipy.ndimage import variance  # NOQA
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/color/colorconv.py b/python/cucim/src/cucim/skimage/color/colorconv.py
index 2d3b42176..ece630492 100644
--- a/python/cucim/src/cucim/skimage/color/colorconv.py
+++ b/python/cucim/src/cucim/skimage/color/colorconv.py
@@ -1201,13 +1201,16 @@ def _get_lab_to_xyz_kernel(xyz_ref_white, name='lab2xyz'):
 
 @channel_as_last_axis()
 def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
-    """CIE-LAB to XYZcolor space conversion.
+    """Convert image in CIE-LAB to XYZ color space.
 
     Parameters
     ----------
     lab : (..., 3, ...) array_like
-        The image in Lab format. By default, the final dimension denotes
+        The input image in CIE-LAB color space.
+        Unless `channel_axis` is set, the final dimension denotes the CIE-LAB
         channels.
+        The L* values range from 0 to 100;
+        the a* and b* values range from -128 to 127.
     illuminant : {"A", "B", "C", "D50", "D55", "D65", "D75", "E"}, optional
         The name of the illuminant (the function is NOT case sensitive).
     observer : {"2", "10", "R"}, optional
@@ -1219,7 +1222,7 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     Returns
     -------
     out : (..., 3, ...) ndarray
-        The image in XYZ format. Same dimensions as input.
+        The image in XYZ color space. Same dimensions as input.
 
     Raises
     ------
@@ -1233,14 +1236,14 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
 
     Notes
     -----
-    By default Observer="2", Illuminant="D65". CIE XYZ tristimulus values x_ref
-    = 95.047, y_ref = 100., z_ref = 108.883. See function 'get_xyz_coords' for
-    a list of supported illuminants.
+    The CIE XYZ tristimulus values are x_ref = 95.047, y_ref = 100., and
+    z_ref = 108.883. See function :func:`~.get_xyz_coords` for a list of
+    supported illuminants.
 
     References
     ----------
     .. [1] http://www.easyrgb.com/index.php?X=MATH&H=07
-    .. [2] https://en.wikipedia.org/wiki/Lab_color_space
+    .. [2] https://en.wikipedia.org/wiki/CIELAB_color_space
     """
     lab = _prepare_colorarray(lab, force_c_contiguous=True,
                               channel_axis=channel_axis)
@@ -1311,13 +1314,16 @@ def rgb2lab(rgb, illuminant="D65", observer="2", *, channel_axis=-1):
 
 @channel_as_last_axis()
 def lab2rgb(lab, illuminant="D65", observer="2", *, channel_axis=-1):
-    """Lab to RGB color space conversion.
+    """Convert image in CIE-LAB to sRGB color space.
 
     Parameters
     ----------
     lab : (..., 3, ...) array_like
-        The image in Lab format. By default, the final dimension denotes
+        The input image in CIE-LAB color space.
+        Unless `channel_axis` is set, the final dimension denotes the CIE-LAB
         channels.
+        The L* values range from 0 to 100;
+        the a* and b* values range from -128 to 127.
     illuminant : {"A", "B", "C", "D50", "D55", "D65", "D75", "E"}, optional
         The name of the illuminant (the function is NOT case sensitive).
     observer : {"2", "10", "R"}, optional
@@ -1339,13 +1345,14 @@ def lab2rgb(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     Notes
     -----
     This function uses lab2xyz and xyz2rgb.
-    By default Observer="2", Illuminant="D65". CIE XYZ tristimulus values
-    x_ref=95.047, y_ref=100., z_ref=108.883. See function `get_xyz_coords` for
-    a list of supported illuminants.
+    The CIE XYZ tristimulus values are x_ref = 95.047, y_ref = 100., and
+    z_ref = 108.883. See function :func:`~.get_xyz_coords` for a list of
+    supported illuminants.
 
     References
     ----------
     .. [1] https://en.wikipedia.org/wiki/Standard_illuminant
+    .. [2] https://en.wikipedia.org/wiki/CIELAB_color_space
     """
     return xyz2rgb(lab2xyz(lab, illuminant, observer))
 
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/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.py b/python/cucim/src/cucim/skimage/filters/_median.py
index be93dcc61..e5d89d2c2 100644
--- a/python/cucim/src/cucim/skimage/filters/_median.py
+++ b/python/cucim/src/cucim/skimage/filters/_median.py
@@ -29,12 +29,16 @@ def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
     ----------
     image : array-like
         Input image.
-    footprint : ndarray, optional
-        If ``behavior=='rank'``, ``footprint`` is a 2-D array of 1's and 0's.
-        If ``behavior=='ndimage'``, ``footprint`` is a N-D array of 1's and 0's
-        with the same number of dimension than ``image``.
-        If None, ``footprint`` will be a N-D array with 3 elements for each
-        dimension (e.g., vector, square, cube, etc.)
+    footprint : ndarray, tuple of int, or None
+        If ``None``, ``footprint`` will be a N-D array with 3 elements for each
+        dimension (e.g., vector, square, cube, etc.). If `footprint` is a
+        tuple of integers, it will be an array of ones with the given shape.
+        Otherwise, if ``behavior=='rank'``, ``footprint`` is a 2-D array of 1's
+        and 0's. If ``behavior=='ndimage'``, ``footprint`` is a N-D array of
+        1's and 0's with the same number of dimension as ``image``.
+        Note that upstream scikit-image currently does not support supplying
+        a tuple for `footprint`. It is added here to avoid overhead of
+        generating a small weights array in cases where it is not needed.
     out : ndarray, (same dtype as image), optional
         If None, a new array is allocated.
     mode : {'reflect', 'constant', 'nearest', 'mirror','‘wrap'}, optional
@@ -121,12 +125,21 @@ def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
         # return generic.median(image, selem=selem, out=out)
 
     if footprint is None:
-        footprint = ndi.generate_binary_structure(image.ndim, image.ndim)
+        footprint_shape = (3,) * image.ndim
+    elif isinstance(footprint, tuple):
+        if len(footprint) != image.ndim:
+            raise ValueError("tuple footprint must have ndim matching image")
+        footprint_shape = footprint
+        footprint = None
+    else:
+        footprint_shape = footprint.shape
 
     if algorithm == 'sorting':
         can_use_histogram = False
     elif algorithm in ['auto', 'histogram']:
-        can_use_histogram, reason = _can_use_histogram(image, footprint)
+        can_use_histogram, reason = _can_use_histogram(
+            image, footprint, footprint_shape
+        )
     else:
         raise ValueError(f"unknown algorithm: {algorithm}")
 
@@ -142,7 +155,7 @@ def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
     use_histogram = can_use_histogram
     if algorithm == 'auto':
         # prefer sorting-based algorithm if footprint shape is small
-        use_histogram = use_histogram and prod(footprint.shape) > 150
+        use_histogram = use_histogram and prod(footprint_shape) > 150
 
     if use_histogram:
         try:
@@ -162,8 +175,16 @@ def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
 
             # TODO: Can't currently pass an output array into _median_hist as a
             #       new array currently needs to be created during padding.
-            temp = _median_hist(image, footprint, mode=mode, cval=cval,
-                                **algorithm_kwargs)
+
+            # pass shape if explicit footprint isn't needed
+            # (use new variable name in case KernelResourceError occurs)
+            temp = _median_hist(
+                image,
+                footprint_shape if footprint is None else footprint,
+                mode=mode,
+                cval=cval,
+                **algorithm_kwargs
+            )
             if output_array_provided:
                 out[:] = temp
             else:
@@ -175,12 +196,13 @@ def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
             # Fall back to sorting-based implementation if we encounter a
             # resource limit (e.g. insufficient shared memory per block).
             warn("Kernel resource error encountered in histogram-based "
-                 f"median kerne: {e}\n"
+                 f"median kernel: {e}\n"
                  "Falling back to sorting-based median instead.")
 
     if algorithm_kwargs:
         warn(f"algorithm_kwargs={algorithm_kwargs} ignored for sorting-based "
              f"algorithm")
 
-    return ndi.median_filter(image, footprint=footprint, output=out, mode=mode,
-                             cval=cval)
+    size = footprint_shape if footprint is None else None
+    return ndi.median_filter(image, size=size, footprint=footprint, output=out,
+                             mode=mode, cval=cval)
diff --git a/python/cucim/src/cucim/skimage/filters/_median_hist.py b/python/cucim/src/cucim/skimage/filters/_median_hist.py
index 0700238c0..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
@@ -230,7 +231,7 @@ def _check_global_scratch_space_size(
     return (n_fine + n_coarse) * cp.dtype(hist_dtype).itemsize
 
 
-def _can_use_histogram(image, footprint):
+def _can_use_histogram(image, footprint, footprint_shape=None):
     """Validate compatibility with histogram-based median.
 
     Parameters
@@ -255,23 +256,32 @@ def _can_use_histogram(image, footprint):
         return False, "Only 8 and 16-bit integer image types (signed or "
         "unsigned)."
 
+    if footprint is None:
+        if footprint_shape is None:
+            raise ValueError(
+                "must provide either footprint or footprint_shape"
+            )
+    else:
+        footprint_shape = footprint.shape
+
     # only odd-sized footprints are supported
-    if not all(s % 2 == 1 for s in footprint.shape):
+    if not all(s % 2 == 1 for s in footprint_shape):
         return False, "footprint must have odd size on both axes"
 
-    if any(s == 1 for s in footprint.shape):
+    if any(s == 1 for s in footprint_shape):
         return False, "footprint must have size >= 3"
 
     # footprint radius can't be larger than the image
     # TODO: need to check if we need this exact restriction
     #       (may be specific to OpenCV's boundary handling)
-    radii = tuple(s // 2 for s in footprint.shape)
+    radii = tuple(s // 2 for s in footprint_shape)
     if any(r > s for r, s in zip(radii, image.shape)):
         return False, "footprint half-width cannot exceed the image extent"
 
-    # only fully populated footprint is supported
-    if not np.all(footprint):  # synchronizes!
-        return False, "footprint must be 1 everywhere"
+    if footprint is not None:
+        # only fully populated footprint is supported
+        if not np.all(footprint):  # synchronizes!
+            return False, "footprint must be 1 everywhere"
 
     return True, None
 
@@ -439,7 +449,19 @@ def _median_hist(image, footprint, output=None, mode='mirror', cval=0,
             "Use of a user-defined output array has not been implemented"
         )
 
-    compatible_image, reason = _can_use_histogram(image, footprint)
+    if footprint is None:
+        footprint_shape = (3,) * image.ndim
+        med_pos = prod(footprint_shape) // 2
+    elif isinstance(footprint, tuple):
+        footprint_shape = footprint
+        footprint = None
+        med_pos = prod(footprint_shape) // 2
+    else:
+        footprint_shape = footprint.shape
+        med_pos = footprint.size // 2
+    compatible_image, reason = _can_use_histogram(
+        image, footprint, footprint_shape
+    )
     if not compatible_image:
         raise ValueError(reason)
 
@@ -451,13 +473,12 @@ def _median_hist(image, footprint, output=None, mode='mirror', cval=0,
     if image.dtype.kind not in 'iu':
         raise ValueError("only integer-type images are accepted")
 
-    radii = tuple(s // 2 for s in footprint.shape)
+    radii = tuple(s // 2 for s in footprint_shape)
     # med_pos is the index corresponding to the median
     # (calculation here assumes all elements of the footprint are True)
-    med_pos = footprint.size // 2
 
     params = _get_kernel_params(
-        image, footprint.shape, value_range, partitions
+        image, footprint_shape, value_range, partitions
     )
 
     # pad as necessary to avoid boundary artifacts
@@ -475,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/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/__init__.py b/python/cucim/src/cucim/skimage/measure/__init__.py
index bd3f16bfc..f24751e14 100644
--- a/python/cucim/src/cucim/skimage/measure/__init__.py
+++ b/python/cucim/src/cucim/skimage/measure/__init__.py
@@ -1,4 +1,6 @@
 from ._blur_effect import blur_effect
+from ._colocalization import (intersection_coeff, manders_coloc_coeff,
+                              manders_overlap_coeff, pearson_corr_coeff)
 from ._label import label
 from ._moments import (centroid, inertia_tensor, inertia_tensor_eigvals,
                        moments, moments_central, moments_coords,
@@ -30,4 +32,8 @@
     "profile_line",
     "label",
     "shannon_entropy",
+    "intersection_coeff",
+    "manders_coloc_coeff",
+    "manders_overlap_coeff",
+    "pearson_corr_coeff",
 ]
diff --git a/python/cucim/src/cucim/skimage/measure/_colocalization.py b/python/cucim/src/cucim/skimage/measure/_colocalization.py
new file mode 100644
index 000000000..f79c43a44
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/measure/_colocalization.py
@@ -0,0 +1,304 @@
+import cupy as cp
+
+from .._shared.utils import as_binary_ndarray, check_shape_equality
+from .._vendored import pearsonr
+
+__all__ = [
+    "pearson_corr_coeff",
+    "manders_coloc_coeff",
+    "manders_overlap_coeff",
+    "intersection_coeff",
+]
+
+
+def pearson_corr_coeff(image0, image1, mask=None):
+    r"""Calculate Pearson's Correlation Coefficient between pixel intensities
+    in channels.
+
+    Parameters
+    ----------
+    image0 : (M, N) ndarray
+        Image of channel A.
+    image1 : (M, N) ndarray
+        Image of channel 2 to be correlated with channel B.
+        Must have same dimensions as `image0`.
+    mask : (M, N) ndarray of dtype bool, optional
+        Only `image0` and `image1` pixels within this region of interest mask
+        are included in the calculation. Must have same dimensions as `image0`.
+
+    Returns
+    -------
+    pcc : float
+        Pearson's correlation coefficient of the pixel intensities between
+        the two images, within the mask if provided.
+    p-value : float
+        Two-tailed p-value.
+
+    Notes
+    -----
+    Pearson's Correlation Coefficient (PCC) measures the linear correlation
+    between the pixel intensities of the two images. Its value ranges from -1
+    for perfect linear anti-correlation to +1 for perfect linear correlation.
+    The calculation of the p-value assumes that the intensities of pixels in
+    each input image are normally distributed.
+
+    Scipy's implementation of Pearson's correlation coefficient is used. Please
+    refer to it for further information and caveats [1]_.
+
+    .. math::
+        r = \frac{\sum (A_i - m_A_i) (B_i - m_B_i)}
+        {\sqrt{\sum (A_i - m_A_i)^2 \sum (B_i - m_B_i)^2}}
+
+    where
+        :math:`A_i` is the value of the :math:`i^{th}` pixel in `image0`
+        :math:`B_i` is the value of the :math:`i^{th}` pixel in `image1`,
+        :math:`m_A_i` is the mean of the pixel values in `image0`
+        :math:`m_B_i` is the mean of the pixel values in `image1`
+
+    A low PCC value does not necessarily mean that there is no correlation
+    between the two channel intensities, just that there is no linear
+    correlation. You may wish to plot the pixel intensities of each of the two
+    channels in a 2D scatterplot and use Spearman's rank correlation if a
+    non-linear correlation is visually identified [2]_. Also consider if you
+    are interested in correlation or co-occurence, in which case a method
+    involving segmentation masks (e.g. MCC or intersection coefficient) may be
+    more suitable [3]_ [4]_.
+
+    Providing the mask of only relevant sections of the image (e.g., cells, or
+    particular cellular compartments) and removing noise is important as the
+    PCC is sensitive to these measures [3]_ [4]_.
+
+    References
+    ----------
+    .. [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html  # noqa
+    .. [2] https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html  # noqa
+    .. [3] Dunn, K. W., Kamocka, M. M., & McDonald, J. H. (2011). A practical
+           guide to evaluating colocalization in biological microscopy.
+           American journal of physiology. Cell physiology, 300(4), C723–C742.
+           https://doi.org/10.1152/ajpcell.00462.2010
+    .. [4] Bolte, S. and Cordelières, F.P. (2006), A guided tour into
+           subcellular colocalization analysis in light microscopy. Journal of
+           Microscopy, 224: 213-232.
+           https://doi.org/10.1111/j.1365-2818.2006.01706.x
+    """
+    if mask is not None:
+        mask = as_binary_ndarray(mask, variable_name="mask")
+        check_shape_equality(image0, image1, mask)
+        image0 = image0[mask]
+        image1 = image1[mask]
+    else:
+        check_shape_equality(image0, image1)
+        # scipy pearsonr function only takes flattened arrays
+        image0 = image0.reshape(-1)
+        image1 = image1.reshape(-1)
+
+    return pearsonr(image0, image1, disable_checks=True)
+
+
+def manders_coloc_coeff(image0, image1_mask, mask=None):
+    r"""Manders' colocalization coefficient between two channels.
+
+    Parameters
+    ----------
+    image0 : (M, N) ndarray
+        Image of channel A. All pixel values should be non-negative.
+    image1_mask : (M, N) ndarray of dtype bool
+        Binary mask with segmented regions of interest in channel B.
+        Must have same dimensions as `image0`.
+    mask : (M, N) ndarray of dtype bool, optional
+        Only `image0` pixel values within this region of interest mask are
+        included in the calculation.
+        Must have same dimensions as `image0`.
+
+    Returns
+    -------
+    mcc : float
+        Manders' colocalization coefficient.
+
+    Notes
+    -----
+    Manders' Colocalization Coefficient (MCC) is the fraction of total
+    intensity of a certain channel (channel A) that is within the segmented
+    region of a second channel (channel B) [1]_. It ranges from 0 for no
+    colocalisation to 1 for complete colocalization. It is also referred to
+    as M1 and M2.
+
+    MCC is commonly used to measure the colocalization of a particular protein
+    in a subceullar compartment. Typically a segmentation mask for channel B
+    is generated by setting a threshold that the pixel values must be above
+    to be included in the MCC calculation. In this implementation,
+    the channel B mask is provided as the argument `image1_mask`, allowing
+    the exact segmentation method to be decided by the user beforehand.
+
+    The implemented equation is:
+
+    .. math::
+        r = \frac{\sum A_{i,coloc}}{\sum A_i}
+
+    where
+        :math:`A_i` is the value of the :math:`i^{th}` pixel in `image0`
+        :math:`A_{i,coloc} = A_i` if :math:`Bmask_i > 0`
+        :math:`Bmask_i` is the value of the :math:`i^{th}` pixel in
+        `mask`
+
+    MCC is sensitive to noise, with diffuse signal in the first channel
+    inflating its value. Images should be processed to remove out of focus and
+    background light before the MCC is calculated [2]_.
+
+    References
+    ----------
+    .. [1] Manders, E.M.M., Verbeek, F.J. and Aten, J.A. (1993), Measurement of
+           co-localization of objects in dual-colour confocal images. Journal
+           of Microscopy, 169: 375-382.
+           https://doi.org/10.1111/j.1365-2818.1993.tb03313.x
+           https://imagej.net/media/manders.pdf
+    .. [2] Dunn, K. W., Kamocka, M. M., & McDonald, J. H. (2011). A practical
+           guide to evaluating colocalization in biological microscopy.
+           American journal of physiology. Cell physiology, 300(4), C723–C742.
+           https://doi.org/10.1152/ajpcell.00462.2010
+
+    """
+    image1_mask = as_binary_ndarray(image1_mask, variable_name="image1_mask")
+    if mask is not None:
+        mask = as_binary_ndarray(mask, variable_name="mask")
+        check_shape_equality(image0, image1_mask, mask)
+        image0 = image0[mask]
+        image1_mask = image1_mask[mask]
+    else:
+        check_shape_equality(image0, image1_mask)
+    # check non-negative image
+    if image0.min() < 0:
+        raise ValueError("image contains negative values")
+
+    img_sum = cp.sum(image0)
+    if (img_sum == 0):
+        return 0
+    return cp.sum(image0 * image1_mask) / img_sum
+
+
+@cp.fuse()
+def _get_manders_overlap_coeff(image0, image1):
+    denom = (cp.sum(cp.square(image0)) * (cp.sum(cp.square(image1)))) ** 0.5
+    return cp.sum(cp.multiply(image0, image1)) / denom
+
+
+def manders_overlap_coeff(image0, image1, mask=None):
+    r"""Manders' overlap coefficient
+
+    Parameters
+    ----------
+    image0 : (M, N) ndarray
+        Image of channel A. All pixel values should be non-negative.
+    image1 : (M, N) ndarray
+        Image of channel B. All pixel values should be non-negative.
+        Must have same dimensions as `image0`
+    mask : (M, N) ndarray of dtype bool, optional
+        Only `image0` and `image1` pixel values within this region of interest
+        mask are included in the calculation.
+        Must have same dimensions as `image0`.
+
+    Returns
+    -------
+    moc: float
+        Manders' Overlap Coefficient of pixel intensities between the two
+        images.
+
+    Notes
+    -----
+    Manders' Overlap Coefficient (MOC) is given by the equation [1]_:
+
+    .. math::
+        r = \frac{\sum A_i B_i}{\sqrt{\sum A_i^2 \sum B_i^2}}
+
+    where
+        :math:`A_i` is the value of the :math:`i^{th}` pixel in `image0`
+        :math:`B_i` is the value of the :math:`i^{th}` pixel in `image1`
+
+    It ranges between 0 for no colocalization and 1 for complete colocalization
+    of all pixels.
+
+    MOC does not take into account pixel intensities, just the fraction of
+    pixels that have positive values for both channels[2]_ [3]_. Its usefulness
+    has been criticized as it changes in response to differences in both
+    co-occurence and correlation and so a particular MOC value could indicate
+    a wide range of colocalization patterns [4]_ [5]_.
+
+    References
+    ----------
+    .. [1] Manders, E.M.M., Verbeek, F.J. and Aten, J.A. (1993), Measurement of
+           co-localization of objects in dual-colour confocal images. Journal
+           of Microscopy, 169: 375-382.
+           https://doi.org/10.1111/j.1365-2818.1993.tb03313.x
+           https://imagej.net/media/manders.pdf
+    .. [2] Dunn, K. W., Kamocka, M. M., & McDonald, J. H. (2011). A practical
+           guide to evaluating colocalization in biological microscopy.
+           American journal of physiology. Cell physiology, 300(4), C723–C742.
+           https://doi.org/10.1152/ajpcell.00462.2010
+    .. [3] Bolte, S. and Cordelières, F.P. (2006), A guided tour into
+           subcellular colocalization analysis in light microscopy. Journal of
+           Microscopy, 224: 213-232.
+           https://doi.org/10.1111/j.1365-2818.2006.01
+    .. [4] Adler J, Parmryd I. (2010), Quantifying colocalization by
+           correlation: the Pearson correlation coefficient is
+           superior to the Mander's overlap coefficient. Cytometry A.
+           Aug;77(8):733-42.https://doi.org/10.1002/cyto.a.20896
+    .. [5] Adler, J, Parmryd, I. Quantifying colocalization: The case for
+           discarding the Manders overlap coefficient. Cytometry. 2021; 99:
+           910– 920. https://doi.org/10.1002/cyto.a.24336
+
+    """
+    if mask is not None:
+        mask = as_binary_ndarray(mask, variable_name="mask")
+        check_shape_equality(image0, image1, mask)
+        image0 = image0[mask]
+        image1 = image1[mask]
+    else:
+        check_shape_equality(image0, image1)
+
+    # check non-negative image
+    if image0.min() < 0:
+        raise ValueError("image0 contains negative values")
+    if image1.min() < 0:
+        raise ValueError("image1 contains negative values")
+
+    return _get_manders_overlap_coeff(image0, image1)
+
+
+def intersection_coeff(image0_mask, image1_mask, mask=None):
+    r"""Fraction of a channel's segmented binary mask that overlaps with a
+    second channel's segmented binary mask.
+
+    Parameters
+    ----------
+    image0_mask : (M, N) ndarray of dtype bool
+        Image mask of channel A.
+    image1_mask : (M, N) ndarray of dtype bool
+        Image mask of channel B.
+        Must have same dimensions as `image0_mask`.
+    mask : (M, N) ndarray of dtype bool, optional
+        Only `image0_mask` and `image1_mask` pixels within this region of
+        interest
+        mask are included in the calculation.
+        Must have same dimensions as `image0_mask`.
+
+    Returns
+    -------
+    Intersection coefficient, float
+        Fraction of `image0_mask` that overlaps with `image1_mask`.
+
+    """
+    image0_mask = as_binary_ndarray(image0_mask, variable_name="image0_mask")
+    image1_mask = as_binary_ndarray(image1_mask, variable_name="image1_mask")
+    if mask is not None:
+        mask = as_binary_ndarray(mask, variable_name="mask")
+        check_shape_equality(image0_mask, image1_mask, mask)
+        image0_mask = image0_mask[mask]
+        image1_mask = image1_mask[mask]
+    else:
+        check_shape_equality(image0_mask, image1_mask)
+
+    nonzero_image0 = cp.count_nonzero(image0_mask)
+    if nonzero_image0 == 0:
+        return 0
+    nonzero_joint = cp.count_nonzero(cp.logical_and(image0_mask, image1_mask))
+    return nonzero_joint / nonzero_image0
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/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/__init__.py b/python/cucim/src/cucim/skimage/morphology/__init__.py
index 91f85d1eb..b4648e0a9 100644
--- a/python/cucim/src/cucim/skimage/morphology/__init__.py
+++ b/python/cucim/src/cucim/skimage/morphology/__init__.py
@@ -6,8 +6,8 @@
 from .gray import (black_tophat, closing, dilation, erosion, opening,
                    white_tophat)
 from .grayreconstruct import reconstruction
-from .isotropic import (isotropic_dilation, isotropic_erosion,
-                        isotropic_opening, isotropic_closing)
+from .isotropic import (isotropic_closing, isotropic_dilation,
+                        isotropic_erosion, isotropic_opening)
 from .misc import remove_small_holes, remove_small_objects
 
 __all__ = [
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/morphology/isotropic.py b/python/cucim/src/cucim/skimage/morphology/isotropic.py
index 33dea9f19..4e47b1aaa 100644
--- a/python/cucim/src/cucim/skimage/morphology/isotropic.py
+++ b/python/cucim/src/cucim/skimage/morphology/isotropic.py
@@ -2,6 +2,7 @@
 Binary morphological operations
 """
 import cupy as cp
+
 from cucim.core.operations.morphology import distance_transform_edt
 
 
diff --git a/python/cucim/src/cucim/skimage/registration/_phase_cross_correlation.py b/python/cucim/src/cucim/skimage/registration/_phase_cross_correlation.py
index 98cdd4f19..0937125d0 100644
--- a/python/cucim/src/cucim/skimage/registration/_phase_cross_correlation.py
+++ b/python/cucim/src/cucim/skimage/registration/_phase_cross_correlation.py
@@ -3,9 +3,11 @@
 http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation
 """
 
+import itertools
 import math
 
 import cupy as cp
+import cupyx.scipy.ndimage as ndi
 import numpy as np
 
 from .._shared.fft import fftmodule as fft
@@ -113,8 +115,75 @@ def _compute_error(cross_correlation_max, src_amp, target_amp):
     return cp.sqrt(cp.abs(error))
 
 
+def _disambiguate_shift(reference_image, moving_image, shift):
+    """Determine the correct real-space shift based on periodic shift.
+
+    When determining a translation shift from phase cross-correlation in
+    Fourier space, the shift is only correct to within a period of the image
+    size along each axis, resulting in $2^n$ possible shifts, where $n$ is the
+    number of dimensions of the image. This function checks the
+    cross-correlation in real space for each of those shifts, and returns the
+    one with the highest cross-correlation.
+
+    The strategy we use is to perform the shift on the moving image *using the
+    'grid-wrap' mode* in `scipy.ndimage`. The moving image's original borders
+    then define $2^n$ quadrants, which we cross-correlate with the reference
+    image in turn using slicing. The entire operation is thus $O(2^n + m)$,
+    where $m$ is the number of pixels in the image (and typically dominates).
+
+    Parameters
+    ----------
+    reference_image : numpy array
+        The reference (non-moving) image.
+    moving_image : numpy array
+        The moving image: applying the shift to this image overlays it on the
+        reference image. Must be the same shape as the reference image.
+    shift : tuple of float
+        The shift to apply to each axis of the moving image, *modulo* image
+        size. The length of ``shift`` must be equal to ``moving_image.ndim``.
+
+    Returns
+    -------
+    real_shift : tuple of float
+        The shift disambiguated in real space.
+    """
+    shape = reference_image.shape
+    positive_shift = [shift_i % s for shift_i, s in zip(shift, shape)]
+    negative_shift = [shift_i - s
+                      for shift_i, s in zip(positive_shift, shape)]
+    subpixel = any(s % 1 != 0 for s in shift)
+    interp_order = 3 if subpixel else 0
+    shifted = ndi.shift(
+        moving_image, shift, mode='grid-wrap', order=interp_order
+    )
+    indices = tuple(round(s) for s in positive_shift)
+    splits_per_dim = [(slice(0, i), slice(i, None)) for i in indices]
+    max_corr = -1.0
+    max_slice = None
+    for test_slice in itertools.product(*splits_per_dim):
+        reference_tile = cp.reshape(reference_image[test_slice], -1)
+        moving_tile = cp.reshape(shifted[test_slice], -1)
+        corr = -1.0
+        if reference_tile.size > 2:
+            corr = float(cp.corrcoef(reference_tile, moving_tile)[0, 1])
+        if corr > max_corr:
+            max_corr = corr
+            max_slice = test_slice
+    real_shift_acc = []
+    for sl, pos_shift, neg_shift in zip(
+        max_slice, positive_shift, negative_shift
+    ):
+        real_shift_acc.append(pos_shift if sl.stop is None else neg_shift)
+    if not subpixel:
+        real_shift = tuple(map(int, real_shift_acc))
+    else:
+        real_shift = tuple(real_shift_acc)
+    return real_shift
+
+
 def phase_cross_correlation(reference_image, moving_image, *,
                             upsample_factor=1, space="real",
+                            disambiguate=False,
                             return_error=True, reference_mask=None,
                             moving_mask=None, overlap_ratio=0.3,
                             normalization="phase"):
@@ -145,10 +214,16 @@ def phase_cross_correlation(reference_image, moving_image, *,
         data will bypass FFT of input data. Case insensitive. Not
         used if any of ``reference_mask`` or ``moving_mask`` is not
         None.
+    disambiguate : bool
+        The shift returned by this function is only accurate *modulo* the
+        image shape, due to the periodic nature of the Fourier transform. If
+        this parameter is set to ``True``, the *real* space cross-correlation
+        is computed for each possible shift, and the shift with the highest
+        cross-correlation within the overlapping area is returned.
     return_error : bool, optional
-        Returns error and phase difference if on, otherwise only
-        shifts are returned. Has noeffect if any of ``reference_mask`` or
-        ``moving_mask`` is not None. In this case only shifts is returned.
+        Returns error and phase difference if "always" is given. If False, or
+        either ``reference_mask`` or ``moving_mask`` are given, only the shift
+        is returned.
     reference_mask : ndarray
         Boolean mask for ``reference_image``. The mask should evaluate
         to ``True`` (or 1) on valid pixels. ``reference_mask`` should
@@ -172,10 +247,10 @@ def phase_cross_correlation(reference_image, moving_image, *,
 
     Returns
     -------
-    shifts : ndarray
+    shift : tuple
         Shift vector (in pixels) required to register ``moving_image``
         with ``reference_image``. Axis ordering is consistent with
-        numpy (e.g. Z, Y, X)
+        the axis order of the input array.
     error : float
         Translation invariant normalized RMS error between
         ``reference_image`` and ``moving_image``.
@@ -250,14 +325,12 @@ def phase_cross_correlation(reference_image, moving_image, *,
     cross_correlation = fft.ifftn(image_product)
 
     # Locate maximum
-    maxima = cp.unravel_index(
-        cp.argmax(cp.abs(cross_correlation)), cross_correlation.shape
+    maxima = np.unravel_index(
+        int(cp.argmax(cp.abs(cross_correlation))), cross_correlation.shape
     )
-    midpoints = cp.asarray([np.fix(axis_size / 2) for axis_size in shape])
-
-    float_dtype = image_product.real.dtype
-    shifts = cp.stack([m.astype(float_dtype, copy=False) for m in maxima])
-    shifts[shifts > midpoints] -= cp.asarray(shape)[shifts > midpoints]
+    midpoint = tuple(float(axis_size // 2) for axis_size in shape)
+    shift = tuple(_max - axis_size if _max > mid else _max
+                  for _max, mid, axis_size in zip(maxima, midpoint, shape))
 
     if upsample_factor == 1:
         if return_error:
@@ -271,29 +344,29 @@ def phase_cross_correlation(reference_image, moving_image, *,
     # If upsampling > 1, then refine estimate with matrix multiply DFT
     else:
         # Initial shift estimate in upsampled grid
-        shifts = cp.around(shifts * upsample_factor) / upsample_factor
+        # shift = cp.around(shift * upsample_factor) / upsample_factor
+        upsample_factor = float(upsample_factor)
+        shift = tuple(round(s * upsample_factor) / upsample_factor
+                      for s in shift)
         upsampled_region_size = math.ceil(upsample_factor * 1.5)
         # Center of output array at dftshift + 1
-        dftshift = np.fix(upsampled_region_size / 2.0)
-        upsample_factor = float(upsample_factor)
+        dftshift = float(upsampled_region_size // 2)
         # Matrix multiply DFT around the current shift estimate
-        sample_region_offset = dftshift - shifts * upsample_factor
+        sample_region_offset = tuple(
+            dftshift - s * upsample_factor for s in shift
+        )
         cross_correlation = _upsampled_dft(image_product.conj(),
                                            upsampled_region_size,
                                            upsample_factor,
                                            sample_region_offset).conj()
 
         # Locate maximum and map back to original pixel grid
-        maxima = cp.unravel_index(cp.argmax(cp.abs(cross_correlation)),
+        maxima = np.unravel_index(int(cp.argmax(cp.abs(cross_correlation))),
                                   cross_correlation.shape)
         CCmax = cross_correlation[maxima]
 
-        maxima = (
-            cp.stack([m.astype(float_dtype, copy=False) for m in maxima])
-            - dftshift
-        )
-
-        shifts += maxima / upsample_factor
+        maxima = tuple(float(m) - dftshift for m in maxima)
+        shift = tuple(s + m / upsample_factor for s, m in zip(shift, maxima))
 
         if return_error:
             src_amp = cp.abs(src_freq)
@@ -305,9 +378,15 @@ def phase_cross_correlation(reference_image, moving_image, *,
 
     # If its only one row or column the shift along that dimension has no
     # effect. We set to zero.
-    for dim in range(src_freq.ndim):
-        if shape[dim] == 1:
-            shifts[dim] = 0
+    shift = tuple(
+        s if axis_size != 1 else 0 for s, axis_size in zip(shift, shape)
+    )
+
+    if disambiguate:
+        if space.lower() != 'real':
+            reference_image = fft.ifftn(reference_image)
+            moving_image = fft.ifftn(moving_image)
+        shift = _disambiguate_shift(reference_image, moving_image, shift)
 
     if return_error:
         # Redirect user to masked_phase_cross_correlation if NaNs are observed
@@ -320,7 +399,7 @@ def phase_cross_correlation(reference_image, moving_image, *,
                 "reference_mask=~np.isnan(reference_image), "
                 "moving_mask=~np.isnan(moving_image))")
 
-        return shifts, _compute_error(CCmax, src_amp, target_amp),\
+        return shift, _compute_error(CCmax, src_amp, target_amp),\
             _compute_phasediff(CCmax)
     else:
-        return shifts
+        return shift
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/segmentation/_clear_border.py b/python/cucim/src/cucim/skimage/segmentation/_clear_border.py
index 6f273a498..07ad03842 100644
--- a/python/cucim/src/cucim/skimage/segmentation/_clear_border.py
+++ b/python/cucim/src/cucim/skimage/segmentation/_clear_border.py
@@ -89,7 +89,7 @@ def clear_border(labels, buffer_size=0, bgval=0, mask=None, *, out=None):
         err_msg = (f'labels and mask should have the same shape but '
                    f'are {out.shape} and {mask.shape}')
         if out.shape != mask.shape:
-            raise(ValueError, err_msg)
+            raise ValueError(err_msg)
         if mask.dtype != bool:
             raise TypeError("mask should be of type bool.")
         borders = ~mask
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],


Test patch:
diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
new file mode 100644
index 000000000..5c18a0b1c
--- /dev/null
+++ b/.github/workflows/test.yaml
@@ -0,0 +1,24 @@
+name: test
+
+on:
+  workflow_dispatch:
+    inputs:
+      branch:
+        required: true
+        type: string
+      date:
+        required: true
+        type: string
+      sha:
+        required: true
+        type: string
+
+jobs:
+  conda-python-tests:
+    secrets: inherit
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-tests.yaml@branch-23.02
+    with:
+      build_type: nightly
+      branch: ${{ inputs.branch }}
+      date: ${{ inputs.date }}
+      sha: ${{ inputs.sha }}
diff --git a/ci/test_python.sh b/ci/test_python.sh
new file mode 100755
index 000000000..1b848fc80
--- /dev/null
+++ b/ci/test_python.sh
@@ -0,0 +1,66 @@
+#!/bin/bash
+# Copyright (c) 2022, NVIDIA CORPORATION.
+
+# Common setup steps shared by Python test jobs
+
+set -euo pipefail
+
+. /opt/conda/etc/profile.d/conda.sh
+
+rapids-logger "Generate Python testing dependencies"
+rapids-dependency-file-generator \
+  --output conda \
+  --file_key test_python \
+  --matrix "cuda=${RAPIDS_CUDA_VERSION%.*};arch=$(arch);py=${RAPIDS_PY_VERSION}" | tee env.yaml
+
+rapids-mamba-retry env create --force -f env.yaml -n test
+
+# Temporarily allow unbound variables for conda activation.
+set +u
+conda activate test
+set -u
+
+rapids-logger "Downloading artifacts from previous jobs"
+CPP_CHANNEL=$(rapids-download-conda-from-s3 cpp)
+PYTHON_CHANNEL=$(rapids-download-conda-from-s3 python)
+
+RAPIDS_TESTS_DIR=${RAPIDS_TESTS_DIR:-"${PWD}/test-results"}
+RAPIDS_COVERAGE_DIR=${RAPIDS_COVERAGE_DIR:-"${PWD}/coverage-results"}
+mkdir -p "${RAPIDS_TESTS_DIR}" "${RAPIDS_COVERAGE_DIR}"
+SUITEERROR=0
+
+rapids-print-env
+
+rapids-mamba-retry install \
+  --channel "${CPP_CHANNEL}" \
+  --channel "${PYTHON_CHANNEL}" \
+  libcucim cucim
+
+rapids-logger "Check GPU usage"
+nvidia-smi
+
+set +e
+
+rapids-logger "pytest cucim"
+pushd python/cucim
+pytest \
+  --cache-clear \
+  --junitxml="${RAPIDS_TESTS_DIR}/junit-cucim.xml" \
+  --numprocesses=8 \
+  --dist=loadscope \
+  --cov-config=.coveragerc \
+  --cov=cucim \
+  --cov-report=xml:"${RAPIDS_COVERAGE_DIR}/cucim-coverage.xml" \
+  --cov-report=term \
+  src \
+  tests/unit \
+  tests/performance
+exitcode=$?
+
+if (( ${exitcode} != 0 )); then
+    SUITEERROR=${exitcode}
+    echo "FAILED: 1 or more tests in cucim"
+fi
+popd
+
+exit ${SUITEERROR}
diff --git a/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py b/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
index 14277b12f..7087611f0 100644
--- a/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
+++ b/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
@@ -1,3 +1,4 @@
+import platform
 import warnings
 
 import cupy as cp
@@ -14,6 +15,10 @@
 from cucim.skimage.exposure.exposure import intensity_range
 from cucim.skimage.util.dtype import dtype_range
 
+# TODO: Some tests fail unexpectedly on ARM.
+ON_AARCH64 = platform.machine() == "aarch64"
+ON_AARCH64_REASON = "TODO: Test fails unexpectedly on ARM."
+
 # Test integer histograms
 # =======================
 
@@ -26,6 +31,7 @@ def test_wrong_source_range():
         )
 
 
+@pytest.mark.xfail(ON_AARCH64, reason=ON_AARCH64_REASON)
 def test_negative_overflow():
     im = cp.array([-1, 100], dtype=cp.int8)
     frequencies, bin_centers = exposure.histogram(im)
@@ -294,6 +300,7 @@ def test_rescale_in_range_clip():
 
 @pytest.mark.parametrize('dtype', [cp.int8, cp.int32, cp.float16, cp.float32,
                                    cp.float64])
+@pytest.mark.xfail(ON_AARCH64, reason=ON_AARCH64_REASON)
 def test_rescale_out_range(dtype):
     """Check that output range is correct.
 
@@ -383,6 +390,7 @@ def test_rescale_output_dtype(out_range, out_dtype):
     assert output_image.dtype == out_dtype
 
 
+@pytest.mark.xfail(ON_AARCH64, reason=ON_AARCH64_REASON)
 def test_rescale_no_overflow():
     image = cp.array([-128, 0, 127], dtype=cp.int8)
     output_image = exposure.rescale_intensity(image, out_range=cp.uint8)
@@ -390,6 +398,7 @@ def test_rescale_no_overflow():
     assert output_image.dtype == cp.uint8
 
 
+@pytest.mark.xfail(ON_AARCH64, reason=ON_AARCH64_REASON)
 def test_rescale_float_output():
     image = cp.array([-128, 0, 127], dtype=cp.int8)
     output_image = exposure.rescale_intensity(image, out_range=(0, 255))
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/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/tests/test_median.py b/python/cucim/src/cucim/skimage/filters/tests/test_median.py
index 6b0d8e002..c510aba20 100644
--- a/python/cucim/src/cucim/skimage/filters/tests/test_median.py
+++ b/python/cucim/src/cucim/skimage/filters/tests/test_median.py
@@ -69,16 +69,22 @@ def test_selem_kwarg_deprecation(image):
         (3, 3), (5, 5), (9, 15), (2, 2), (1, 1), (2, 7), (23, 23), (15, 35),
     ]
 )
+@pytest.mark.parametrize('footprint_tuple', (False, True))
 @pytest.mark.parametrize('out', [None, cp.uint8, cp.float32, 'array'])
-def test_median_behavior(camera, behavior, func, mode, footprint_shape, out):
-    footprint = cp.ones(footprint_shape, dtype=bool)
+def test_median_behavior(
+    camera, behavior, func, mode, footprint_shape, footprint_tuple, out
+):
+    if footprint_tuple:
+        footprint = footprint_shape
+    else:
+        footprint = cp.ones(footprint_shape, dtype=bool)
     cam2 = camera[:, :177]  # use anisotropic size
     assert cam2.dtype == cp.uint8
     if out == 'array':
         out = cp.zeros_like(cam2)
     assert_allclose(
         median(cam2, footprint, mode=mode, behavior=behavior, out=out),
-        func(cam2, size=footprint.shape, mode=mode, output=out),
+        func(cam2, size=footprint_shape, mode=mode, output=out),
     )
 
 
diff --git a/python/cucim/src/cucim/skimage/measure/tests/test_colocalization.py b/python/cucim/src/cucim/skimage/measure/tests/test_colocalization.py
new file mode 100644
index 000000000..ab1c05949
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/measure/tests/test_colocalization.py
@@ -0,0 +1,114 @@
+import cupy as cp
+import numpy as np
+import pytest
+
+from cucim.skimage.measure import (intersection_coeff, manders_coloc_coeff,
+                                   manders_overlap_coeff, pearson_corr_coeff)
+
+
+def test_invalid_input():
+    # images are not same size
+    img1 = cp.array([[i + j for j in range(4)] for i in range(4)])
+    img2 = cp.ones((3, 5, 6))
+    mask = cp.array([[i <= 1 for i in range(5)] for _ in range(5)])
+    non_binary_mask = cp.array([[2 for __ in range(4)] for _ in range(4)])
+
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        pearson_corr_coeff(img1, img1, mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        pearson_corr_coeff(img1, img2)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        pearson_corr_coeff(img1, img1, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        pearson_corr_coeff(img1, img1, non_binary_mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        manders_coloc_coeff(img1, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        manders_coloc_coeff(img1, non_binary_mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        manders_coloc_coeff(img1, img1 > 0, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        manders_coloc_coeff(img1, img1 > 0, non_binary_mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        manders_overlap_coeff(img1, img1, mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        manders_overlap_coeff(img1, img2)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        manders_overlap_coeff(img1, img1, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        manders_overlap_coeff(img1, img1, non_binary_mask)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        intersection_coeff(img1 > 2, img2 > 1, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        intersection_coeff(img1, img2)
+    with pytest.raises(ValueError, match=". must have the same dimensions"):
+        intersection_coeff(img1 > 2, img1 > 1, mask)
+    with pytest.raises(ValueError, match=". array is not of dtype boolean"):
+        intersection_coeff(img1 > 2, img1 > 1, non_binary_mask)
+
+
+def test_pcc():
+    # simple example
+    img1 = cp.array([[i + j for j in range(4)] for i in range(4)])
+    np.testing.assert_allclose(
+        pearson_corr_coeff(img1, img1),
+        (1.0, 0.0),
+        rtol=1e-12,
+    )
+
+    img2 = cp.where(img1 <= 2, 0, img1)
+    np.testing.assert_allclose(
+        pearson_corr_coeff(img1, img2),
+        (0.944911182523068, 3.5667540654536515e-08),
+        rtol=1e-12,
+    )
+
+    # change background of roi and see if values are same
+    roi = cp.where(img1 <= 2, 0, 1)
+    np.testing.assert_allclose(
+        pearson_corr_coeff(img1, img1, roi),
+        pearson_corr_coeff(img1, img2, roi),
+        rtol=1e-12,
+    )
+
+
+def test_mcc():
+    img1 = cp.array([[j for j in range(4)] for i in range(4)])
+    mask = cp.array([[i <= 1 for j in range(4)]for i in range(4)])
+    assert manders_coloc_coeff(img1, mask) == 0.5
+
+    # test negative values
+    img_negativeint = cp.where(img1 == 1, -1, img1)
+    img_negativefloat = img_negativeint / 2.0
+    with pytest.raises(ValueError):
+        manders_coloc_coeff(img_negativeint, mask)
+    with pytest.raises(ValueError):
+        manders_coloc_coeff(img_negativefloat, mask)
+
+
+def test_moc():
+    img1 = cp.ones((4, 4))
+    img2 = 2 * cp.ones((4, 4))
+    assert manders_overlap_coeff(img1, img2) == 1
+
+    # test negative values
+    img_negativeint = cp.where(img1 == 1, -1, img1)
+    img_negativefloat = img_negativeint / 2.0
+    with pytest.raises(ValueError):
+        manders_overlap_coeff(img_negativeint, img2)
+    with pytest.raises(ValueError):
+        manders_overlap_coeff(img1, img_negativeint)
+    with pytest.raises(ValueError):
+        manders_overlap_coeff(img_negativefloat, img2)
+    with pytest.raises(ValueError):
+        manders_overlap_coeff(img1, img_negativefloat)
+    with pytest.raises(ValueError):
+        manders_overlap_coeff(img_negativefloat, img_negativefloat)
+
+
+def test_intersection_coefficient():
+    img1_mask = cp.array([[j <= 1 for j in range(4)] for i in range(4)])
+    img2_mask = cp.array([[i <= 1 for j in range(4)] for i in range(4)])
+    img3_mask = cp.array([[1 for j in range(4)] for i in range(4)])
+    assert intersection_coeff(img1_mask, img2_mask) == 0.5
+    assert intersection_coeff(img1_mask, img3_mask) == 1
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/morphology/tests/test_isotropic.py b/python/cucim/src/cucim/skimage/morphology/tests/test_isotropic.py
index eb70196a7..3ca66550e 100644
--- a/python/cucim/src/cucim/skimage/morphology/tests/test_isotropic.py
+++ b/python/cucim/src/cucim/skimage/morphology/tests/test_isotropic.py
@@ -6,7 +6,6 @@
 from cucim.skimage import color, morphology
 from cucim.skimage.util import img_as_bool
 
-
 img = color.rgb2gray(cp.asarray(data.astronaut()))
 bw_img = img > 100 / 255.
 
diff --git a/python/cucim/src/cucim/skimage/registration/tests/test_phase_cross_correlation.py b/python/cucim/src/cucim/skimage/registration/tests/test_phase_cross_correlation.py
index 20696552a..5561b66c2 100644
--- a/python/cucim/src/cucim/skimage/registration/tests/test_phase_cross_correlation.py
+++ b/python/cucim/src/cucim/skimage/registration/tests/test_phase_cross_correlation.py
@@ -1,13 +1,15 @@
+import itertools
+
 import cupy as cp
+import numpy as np
 import pytest
 from cupy.testing import assert_allclose
 from cupyx.scipy.ndimage import fourier_shift
-from skimage.data import camera
+from skimage.data import camera, eagle
 
 from cucim.skimage import img_as_float
 from cucim.skimage._shared._warnings import expected_warnings
 from cucim.skimage._shared.fft import fftmodule as fft
-from cucim.skimage._shared.utils import _supported_float_type
 from cucim.skimage.data import binary_blobs
 from cucim.skimage.registration._phase_cross_correlation import (
     _upsampled_dft, phase_cross_correlation)
@@ -67,7 +69,8 @@ def test_real_input(dtype):
     result, error, diffphase = phase_cross_correlation(reference_image,
                                                        shifted_image,
                                                        upsample_factor=100)
-    assert result.dtype == _supported_float_type(dtype)
+    assert isinstance(result, tuple)
+    assert all(isinstance(s, float) for s in result)
     assert_allclose(result[:2], -cp.array(subpixel_shift), atol=0.05)
 
 
@@ -169,3 +172,40 @@ def test_mismatch_offsets_size():
     with pytest.raises(ValueError):
         _upsampled_dft(cp.ones((4, 4)), 3,
                        axis_offsets=[3, 2, 1, 4])
+
+
+@pytest.mark.parametrize(
+    ('shift0', 'shift1'),
+    itertools.product((100, -100, 350, -350), (100, -100, 350, -350)),
+)
+def test_disambiguate_2d(shift0, shift1):
+    image = cp.array(eagle()[500:, 900:])  # use a highly textured image region
+    shift = (shift0, shift1)
+    origin0 = []
+    for s in shift:
+        if s > 0:
+            origin0.append(0)
+        else:
+            origin0.append(-s)
+    origin1 = np.array(origin0) + shift
+    slice0 = tuple(slice(o, o + 450) for o in origin0)
+    slice1 = tuple(slice(o, o + 450) for o in origin1)
+    reference = image[slice0]
+    moving = image[slice1]
+    computed_shift, _, _ = phase_cross_correlation(
+        reference, moving, disambiguate=True, return_error='always'
+    )
+    np.testing.assert_equal(shift, computed_shift)
+
+
+def test_disambiguate_zero_shift():
+    """When the shift is 0, disambiguation becomes degenerate.
+    Some quadrants become size 0, which prevents computation of
+    cross-correlation. This test ensures that nothing bad happens in that
+    scenario.
+    """
+    image = cp.array(camera())
+    computed_shift, _, _ = phase_cross_correlation(
+        image, image, disambiguate=True, return_error='always'
+    )
+    assert computed_shift == (0, 0)
diff --git a/python/cucim/src/cucim/skimage/segmentation/tests/test_random_walker.py b/python/cucim/src/cucim/skimage/segmentation/tests/test_random_walker.py
index 1b0d0b21c..f0860182a 100644
--- a/python/cucim/src/cucim/skimage/segmentation/tests/test_random_walker.py
+++ b/python/cucim/src/cucim/skimage/segmentation/tests/test_random_walker.py
@@ -1,5 +1,8 @@
+import platform
+
 import cupy as cp
 import numpy as np
+import pytest
 
 from cucim.skimage._shared import testing
 from cucim.skimage._shared._warnings import expected_warnings
@@ -10,6 +13,10 @@
 # SciPy >= 1.8 is installed.
 cupy_warning = r"Please use `spmatrix` from the `scipy.sparse` |\A\Z"
 
+# TODO: Some tests fail unexpectedly on ARM.
+ON_AARCH64 = platform.machine() == "aarch64"
+ON_AARCH64_REASON = "TODO: Test fails unexpectedly on ARM."
+
 
 def make_2d_syntheticdata(lx, ly=None):
     if ly is None:
@@ -325,6 +332,7 @@ def test_spacing_1():
     assert (labels_aniso2[26:34, 13:17, 13:17] == 2).all()
 
 
+@pytest.mark.xfail(ON_AARCH64, reason=ON_AARCH64_REASON)
 def test_trivial_cases():
     # When all voxels are labeled
     img = cp.ones((10, 10))
diff --git a/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py b/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
index 8048f7c0f..2b2891542 100644
--- a/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
+++ b/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
@@ -147,13 +147,17 @@ def test_similarity_estimation():
 
 
 def test_3d_similarity_estimation():
-    src_points = cp.random.rand(1000, 3)
+    # Note: some choices of seed will produce poorly conditioned values,
+    #       leading to NaN in the estimate. We fix the seed here to avoid this.
+    seed = 10
+    rng = cp.random.default_rng(seed)
+    src_points = rng.random(size=(1000, 3))
 
     # Random transformation for testing
-    angles = cp.random.random((3,)) * 2 * np.pi - np.pi
-    scale = cp.random.randint(0, 20)
+    angles = rng.random(size=(3,)) * 2 * np.pi - np.pi
+    scale = rng.integers(0, 20, size=(1,))[0]
     rotation_matrix = cp.array(_euler_rotation_matrix(angles)) * scale
-    translation_vector = cp.random.random((3,))
+    translation_vector = rng.random((3,))
     dst_points = []
     for pt in src_points:
         pt_r = pt.reshape(3, 1)