Repository: rapidsai/cucim (language: python)

Issue:
Forward-merge branch-23.06 to branch-23.08
Forward-merge triggered by push to `branch-23.06` that creates a PR to keep `branch-23.08` up-to-date. If this PR is unable to be immediately merged due to conflicts, it will remain open for the team to manually merge.
[FEA] Missing IntelliSense capability in VSCode
**Is your feature request related to a problem? Please describe.**
After installing cucim, unfortunately I cannot get IntelliSense features like auto-complete in Visual Studio Code. According to https://github.com/microsoft/pylance-release/issues/242 this is due to the fact that the Pylance extension does not support IntelliSense for compiled libraries.

**Describe the solution you'd like**
As described in https://github.com/microsoft/pylance-release/issues/242 it would be necessary to provide a type stub ("pyi") file for the library interface.
I would very much appreciate if such a type stub file could be provided.

**Describe alternatives you've considered**
Interestingly, IntelliSense features are working for me in PyCharm. Haven't figured out the reason why PyCharm can auto-generate the type stubs and VSCode can't..

Disable cuFile / GDS on SBSA ( `aarch64` )
As there are no cuFile / GDS binaries for SSBA ( `aarch64` ) currently, it would be good to disable downloading cuFile (or needing to build with it) in the SSBA ( `aarch64` ) case.
Pin `libwebp-base`
This PR pins `libwebp-base`, since newer versions conflict with other RAPIDS packages.

The screenshot below shows the solver issues.

![image](https://user-images.githubusercontent.com/7400326/231845892-02c19e96-4671-4568-86b8-c96c095583b4.png)

It appears that `cucim` and `cuxfilter` are currently incompatible, likely due to the changes outlined in the issue/commit below:

- https://github.com/conda-forge/conda-forge.github.io/issues/673
- https://github.com/conda-forge/libjpeg-turbo-feedstock/commit/6cee8324c6109060155cfe31c4175a62f7628eef
Flaky test: `test_read_region_cuda_memleak`
Previously we saw ( https://github.com/rapidsai/cucim/pull/593#issuecomment-1662504461 ) this failure [on CI]( https://github.com/rapidsai/cucim/actions/runs/5740435207/job/15558989875?pr=593#step:7:20114 )

```python
=================================== FAILURES ===================================
________________________ test_read_region_cuda_memleak _________________________
[gw1] linux -- Python 3.10.12 /opt/conda/envs/test/bin/python3.10
tests/performance/clara/test_read_region_memory_usage.py:44: in test_read_region_cuda_memleak
    assert mem_usage_history[4] - mem_usage_history[1] < 180.0
E   assert (14623.0 - 14257.0) < 180.0
----------------------------- Captured stdout call -----------------------------
[12999.0, 14257.0, 14569.0, 14617.0, 14623.0, 14621.0]
```

In discussion offline it sounds like this test can be a bit flaky. Wondering if there are ways we can make it less flaky or if we should just `skip` it


Discussion / hints:

Thanks, @alexkrz. The use of `.pyi` files was introduced in recent months to scikit-image and should appear in the 0.20 release. 

Let me update to use .pyi files here as well (port changes from https://github.com/scikit-image/scikit-image/pull/6577). I will target this for the 2023.04.00 milestone.

CC: @AjayThorve for cuxfilter
Thanks all! 🙏


Fix patch:
diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
index 9b4e8a7ac..f4349bb89 100644
--- a/.github/workflows/build.yaml
+++ b/.github/workflows/build.yaml
@@ -28,7 +28,7 @@ concurrency:
 jobs:
   cpp-build:
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.08
     with:
       build_type: ${{ inputs.build_type || 'branch' }}
       branch: ${{ inputs.branch }}
@@ -38,7 +38,7 @@ jobs:
     if: github.ref_type == 'branch'
     needs: [python-build]
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/custom-job.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/custom-job.yaml@branch-23.08
     with:
       arch: "amd64"
       branch: ${{ inputs.branch }}
@@ -51,7 +51,7 @@ jobs:
   python-build:
     needs: [cpp-build]
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-build.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-build.yaml@branch-23.08
     with:
       build_type: ${{ inputs.build_type || 'branch' }}
       branch: ${{ inputs.branch }}
@@ -60,7 +60,7 @@ jobs:
   upload-conda:
     needs: [cpp-build, python-build]
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/conda-upload-packages.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-upload-packages.yaml@branch-23.08
     with:
       build_type: ${{ inputs.build_type || 'branch' }}
       branch: ${{ inputs.branch }}
diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml
index bd614e619..150debb86 100644
--- a/.github/workflows/pr.yaml
+++ b/.github/workflows/pr.yaml
@@ -18,32 +18,32 @@ jobs:
       - conda-python-tests
       - docs-build
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/pr-builder.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/pr-builder.yaml@branch-23.08
   checks:
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/checks.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/checks.yaml@branch-23.08
   conda-cpp-build:
     needs: checks
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-cpp-build.yaml@branch-23.08
     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.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-build.yaml@branch-23.08
     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.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-tests.yaml@branch-23.08
     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.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/custom-job.yaml@branch-23.08
     with:
       build_type: pull-request
       node_type: "gpu-v100-latest-1"
diff --git a/CHANGELOG.md b/CHANGELOG.md
index a2cf4758b..b29bc2e00 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,30 @@
+# cuCIM 23.08.00 (9 Aug 2023)
+
+## 🚨 Breaking Changes
+
+- Sync cuCIM API with scikit-image 0.21 ([#573](https://github.com/rapidsai/cucim/pull/573)) [@grlee77](https://github.com/grlee77)
+
+## 🐛 Bug Fixes
+
+- Remove libwebp-base upper bound. ([#599](https://github.com/rapidsai/cucim/pull/599)) [@bdice](https://github.com/bdice)
+- Fix ignore_run_exports_from. ([#596](https://github.com/rapidsai/cucim/pull/596)) [@bdice](https://github.com/bdice)
+- Add ignore_run_exports for CUDA 11 ([#593](https://github.com/rapidsai/cucim/pull/593)) [@raydouglass](https://github.com/raydouglass)
+- Use linalg &amp; inline `_get_manders_overlap_coeff` ([#578](https://github.com/rapidsai/cucim/pull/578)) [@jakirkham](https://github.com/jakirkham)
+- Fix canny and butterworth (recent CuPy and NumPy compatibility) ([#574](https://github.com/rapidsai/cucim/pull/574)) [@grlee77](https://github.com/grlee77)
+
+## 🛠️ Improvements
+
+- Postponed the assessment of GPU memory for testing ([#601](https://github.com/rapidsai/cucim/pull/601)) [@gigony](https://github.com/gigony)
+- Do not use x86_64 GDS binaries for aarch64 ([#590](https://github.com/rapidsai/cucim/pull/590)) [@gigony](https://github.com/gigony)
+- remove checks for versions of scikit-image that are no longer supported ([#587](https://github.com/rapidsai/cucim/pull/587)) [@grlee77](https://github.com/grlee77)
+- Allow scikit-image 0.21.0 ([#580](https://github.com/rapidsai/cucim/pull/580)) [@jakirkham](https://github.com/jakirkham)
+- Drop unneeded selector on `libwebp-base` ([#579](https://github.com/rapidsai/cucim/pull/579)) [@jakirkham](https://github.com/jakirkham)
+- switch from bundled lazy loading code to the public lazy_loader package ([#575](https://github.com/rapidsai/cucim/pull/575)) [@grlee77](https://github.com/grlee77)
+- Sync cuCIM API with scikit-image 0.21 ([#573](https://github.com/rapidsai/cucim/pull/573)) [@grlee77](https://github.com/grlee77)
+- cuCIM: Build CUDA 12 packages ([#572](https://github.com/rapidsai/cucim/pull/572)) [@jakirkham](https://github.com/jakirkham)
+- use rapids-upload-docs script ([#570](https://github.com/rapidsai/cucim/pull/570)) [@AyodeAwe](https://github.com/AyodeAwe)
+- Remove documentation build scripts for Jenkins ([#567](https://github.com/rapidsai/cucim/pull/567)) [@ajschmidt8](https://github.com/ajschmidt8)
+
 # cuCIM 23.06.00 (7 Jun 2023)
 
 ## 🚨 Breaking Changes
diff --git a/README.md b/README.md
index 081702200..1249e06dd 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-23.06/notebooks/Welcome.ipynb))
+Please check out our [Welcome](notebooks/Welcome.ipynb) notebook ([NBViewer](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-23.08/notebooks/Welcome.ipynb))
 
 #### Downloading sample images
 
diff --git a/VERSION b/VERSION
index 4f3f0be7c..e7ab23867 100644
--- a/VERSION
+++ b/VERSION
@@ -1,1 +1,1 @@
-23.06.00
+23.08.00
diff --git a/ci/build_docs.sh b/ci/build_docs.sh
index 33d14a1dc..9b228b4eb 100755
--- a/ci/build_docs.sh
+++ b/ci/build_docs.sh
@@ -18,22 +18,22 @@ 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="23.06"
 
 rapids-mamba-retry install \
     --channel "${CPP_CHANNEL}" \
     --channel "${PYTHON_CHANNEL}" \
     cucim libcucim
 
-# Build Python docs
+export RAPIDS_VERSION_NUMBER="23.08"
+export RAPIDS_DOCS_DIR="$(mktemp -d)"
+
 rapids-logger "Build Python docs"
 pushd docs
 sphinx-build -b dirhtml ./source _html
 sphinx-build -b text ./source _text
+mkdir -p "${RAPIDS_DOCS_DIR}/cucim/"{html,txt}
+mv _html/* "${RAPIDS_DOCS_DIR}/cucim/html"
+mv _text/* "${RAPIDS_DOCS_DIR}/cucim/txt"
 popd
 
-if [[ "${RAPIDS_BUILD_TYPE}" != "pull-request" ]]; 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
+rapids-upload-docs
diff --git a/ci/docs/build.sh b/ci/docs/build.sh
deleted file mode 100644
index ed70f8a7d..000000000
--- a/ci/docs/build.sh
+++ /dev/null
@@ -1,59 +0,0 @@
-#!/bin/bash
-# Copyright (c) 2020, NVIDIA CORPORATION.
-#################################
-# cuCIM Docs build script for CI #
-#################################
-
-if [ -z "$PROJECT_WORKSPACE" ]; then
-    echo ">>>> ERROR: Could not detect PROJECT_WORKSPACE in environment"
-    echo ">>>> WARNING: This script contains git commands meant for automated building, do not run locally"
-    exit 1
-fi
-
-export DOCS_WORKSPACE=$WORKSPACE/docs
-export PATH=/opt/conda/bin:/usr/local/cuda/bin:$PATH
-export HOME=$WORKSPACE
-export PROJECT_WORKSPACE=/rapids/cucim
-export PROJECTS=(cucim)
-
-gpuci_logger "Check environment"
-env
-
-gpuci_logger "Check GPU usage"
-nvidia-smi
-
-gpuci_logger "Activate conda env"
-. /opt/conda/etc/profile.d/conda.sh
-conda activate rapids
-
-gpuci_logger "Installing cuCIM / Deps / Docs into new env"
-gpuci_conda_retry install -y -c rapidsai-nightly \
-    rapids-doc-env \
-    cucim
-
-gpuci_logger "Check versions"
-python --version
-$CC --version
-$CXX --version
-
-gpuci_logger "Check conda environment"
-conda info
-conda config --show-sources
-conda list --show-channel-urls
-
-# Build Python docs
-gpuci_logger "Build Sphinx docs"
-cd $PROJECT_WORKSPACE/docs
-make html
-
-#Commit to Website
-cd $DOCS_WORKSPACE
-
-for PROJECT in ${PROJECTS[@]}; do
-    if [ ! -d "api/$PROJECT/$BRANCH_VERSION" ]; then
-        mkdir -p api/$PROJECT/$BRANCH_VERSION
-    fi
-    rm -rf $DOCS_WORKSPACE/api/$PROJECT/$BRANCH_VERSION/*
-done
-
-mv $PROJECT_WORKSPACE/docs/build/html/* $DOCS_WORKSPACE/api/cucim/$BRANCH_VERSION
diff --git a/ci/release/update-version.sh b/ci/release/update-version.sh
index ed13f2f12..e14f389e2 100755
--- a/ci/release/update-version.sh
+++ b/ci/release/update-version.sh
@@ -48,4 +48,4 @@ sed_runner "s#branch-${CURRENT_MAJOR}.${CURRENT_MINOR}#branch-${NEXT_MAJOR}.${NE
 for FILE in .github/workflows/*.yaml; do
   sed_runner "/shared-action-workflows/ s/@.*/@branch-${NEXT_SHORT_TAG}/g" "${FILE}"
 done
-sed_runner "s/VERSION_NUMBER=\".*/VERSION_NUMBER=\"${NEXT_SHORT_TAG}\"/g" ci/build_docs.sh
+sed_runner "s/RAPIDS_VERSION_NUMBER=\".*/RAPIDS_VERSION_NUMBER=\"${NEXT_SHORT_TAG}\"/g" ci/build_docs.sh
diff --git a/conda/environments/all_cuda-118_arch-x86_64.yaml b/conda/environments/all_cuda-118_arch-x86_64.yaml
index 4b06f29aa..b9a4de842 100644
--- a/conda/environments/all_cuda-118_arch-x86_64.yaml
+++ b/conda/environments/all_cuda-118_arch-x86_64.yaml
@@ -10,13 +10,19 @@ dependencies:
 - c-compiler
 - click
 - cmake>=3.23.1,!=3.25.0
-- cudatoolkit=11.8
+- cuda-version=11.8
+- cudatoolkit
 - cupy >=12.0.0
 - cxx-compiler
 - gcc_linux-64=11.*
 - imagecodecs>=2021.6.8
 - ipython
 - jbig
+- lazy_loader >=0.1
+- libcufile-dev=1.4.0.31
+- libcufile=1.4.0.31
+- libnvjpeg-dev=11.6.0.55
+- libnvjpeg=11.6.0.55
 - libwebp-base
 - nbsphinx
 - ninja
@@ -34,7 +40,7 @@ dependencies:
 - pytest>=6.2.4
 - python>=3.8,<3.11
 - recommonmark
-- scikit-image >=0.19.0,<0.21.0a0
+- scikit-image >=0.19.0,<0.22.0a0
 - scipy
 - sphinx<6
 - sysroot_linux-64==2.17
diff --git a/conda/environments/all_cuda-120_arch-x86_64.yaml b/conda/environments/all_cuda-120_arch-x86_64.yaml
new file mode 100644
index 000000000..a4cb89dcd
--- /dev/null
+++ b/conda/environments/all_cuda-120_arch-x86_64.yaml
@@ -0,0 +1,53 @@
+# 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
+- cuda-cudart-dev
+- cuda-nvcc
+- cuda-version=12.0
+- cupy >=12.0.0
+- cxx-compiler
+- gcc_linux-64=11.*
+- imagecodecs>=2021.6.8
+- ipython
+- jbig
+- lazy_loader >=0.1
+- libcufile-dev
+- libnvjpeg-dev
+- libnvjpeg-static
+- libwebp-base
+- nbsphinx
+- ninja
+- numpy >=1.21.3
+- numpydoc
+- 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.22.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-120_arch-x86_64
diff --git a/conda/recipes/cucim/conda_build_config.yaml b/conda/recipes/cucim/conda_build_config.yaml
index ccc49851a..fb04ce5c8 100644
--- a/conda/recipes/cucim/conda_build_config.yaml
+++ b/conda/recipes/cucim/conda_build_config.yaml
@@ -4,5 +4,11 @@ c_compiler_version:
 cxx_compiler_version:
   - 11
 
+cuda_compiler:
+  - cuda-nvcc
+
+cuda11_compiler:
+  - nvcc
+
 sysroot_version:
   - "2.17"
diff --git a/conda/recipes/cucim/meta.yaml b/conda/recipes/cucim/meta.yaml
index afbf439bc..609ca09c5 100644
--- a/conda/recipes/cucim/meta.yaml
+++ b/conda/recipes/cucim/meta.yaml
@@ -1,12 +1,12 @@
 # Copyright (c) 2021-2023, NVIDIA CORPORATION.
 
 {% set version = environ.get('GIT_DESCRIBE_TAG', '0.0.0.dev').lstrip('v') %}
-{% 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] %}
 {% set date_string = environ['RAPIDS_DATE_STRING'] %}
 
+
 package:
   name: cucim
   version: {{ version }}
@@ -17,6 +17,10 @@ source:
 build:
   number: {{ GIT_DESCRIBE_NUMBER }}
   string: cuda{{ cuda_major }}_py{{ py_version }}_{{ date_string }}_{{ GIT_DESCRIBE_HASH }}_{{ GIT_DESCRIBE_NUMBER }}
+  ignore_run_exports_from:
+    {% if cuda_major == "11" %}
+    - {{ compiler('cuda11') }}
+    {% endif %}
   script_env:
     - AWS_ACCESS_KEY_ID
     - AWS_SECRET_ACCESS_KEY
@@ -38,33 +42,40 @@ requirements:
   build:
     - {{ compiler("c") }}
     - {{ compiler("cxx") }}
+    {% if cuda_major == "11" %}
+    - {{ compiler('cuda11') }} ={{ cuda_version }}
+    {% else %}
+    - {{ compiler('cuda') }}
+    {% endif %}
+    - cuda-version ={{ cuda_version }}
     - cmake >=3.23.1,!=3.25.0
     - make
     - ninja
     - sysroot_{{ target_platform }} {{ sysroot_version }}
   host:
     - click
-    - cudatoolkit ={{ cuda_version }}
+    - cuda-version ={{ cuda_version }}
     - cupy >=12.0.0
     - libcucim ={{ version }}
     - numpy 1.21
     - python
-    - scikit-image >=0.19.0,<0.21.0a0
+    - scikit-image >=0.19.0,<0.22.0a0
     - scipy
   run:
-    - {{ pin_compatible('cudatoolkit', max_pin='x', min_pin='x') }}
+    - {{ pin_compatible('cuda-version', max_pin='x', min_pin='x') }}
     - {{ pin_compatible('numpy') }}
     - click
     - cupy >=12.0.0
+    - lazy_loader >=0.1
     - libcucim ={{ version }}
     # - openslide # skipping here but benchmark binary would need openslide library
     - python
-    - scikit-image >=0.19.0,<0.21.0a0
+    - scikit-image >=0.19.0,<0.22.0a0
     - scipy
 
 tests:
   requirements:
-    - cudatoolkit ={{ cuda_version }}
+    - cuda-version ={{ cuda_version }}
   imports:
     - cucim
 
diff --git a/conda/recipes/libcucim/build.sh b/conda/recipes/libcucim/build.sh
index 7d21aebe3..d9da2e18d 100644
--- a/conda/recipes/libcucim/build.sh
+++ b/conda/recipes/libcucim/build.sh
@@ -4,7 +4,6 @@ CUCIM_BUILD_TYPE=${CUCIM_BUILD_TYPE:-release}
 
 echo "CC          : ${CC}"
 echo "CXX         : ${CXX}"
-echo "CUDA        : ${CUDA}"
 
 # CUDA needs to include $PREFIX/include as system include path
 export CUDAFLAGS="-isystem $BUILD_PREFIX/include -isystem $PREFIX/include "
diff --git a/conda/recipes/libcucim/conda_build_config.yaml b/conda/recipes/libcucim/conda_build_config.yaml
index 6c2368ed8..6e58b05e2 100644
--- a/conda/recipes/libcucim/conda_build_config.yaml
+++ b/conda/recipes/libcucim/conda_build_config.yaml
@@ -4,8 +4,24 @@ c_compiler_version:
 cxx_compiler_version:
   - 11
 
+cuda_compiler:
+  - cuda-nvcc
+
+cuda11_compiler:
+  - nvcc
+
 sysroot_version:
   - "2.17"
 
-libwebp_base_version:
-  - "<1.3a0"
+# The CTK libraries below are missing from the conda-forge::cudatoolkit package
+# for CUDA 11. The "*_host_*" version specifiers correspond to `11.8` packages
+# and the "*_run_*" version specifiers correspond to `11.x` packages.
+
+cuda11_libcufile_host_version:
+  - "1.4.0.31"
+
+cuda11_libcufile_run_version:
+  - ">=1.0.0.82,<=1.4.0.31"
+
+cuda11_libnvjpeg_host_version:
+  - "11.6.0.55"
diff --git a/conda/recipes/libcucim/meta.yaml b/conda/recipes/libcucim/meta.yaml
index 11486062d..4140bbf1d 100644
--- a/conda/recipes/libcucim/meta.yaml
+++ b/conda/recipes/libcucim/meta.yaml
@@ -1,11 +1,11 @@
 # Copyright (c) 2021-2023, NVIDIA CORPORATION.
 
 {% set version = environ.get('GIT_DESCRIBE_TAG', '0.0.0.dev').lstrip('v') %}
-{% 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] %}
 {% set date_string = environ['RAPIDS_DATE_STRING'] %}
 
+
 package:
   name: libcucim
   version: {{ version }}
@@ -18,6 +18,10 @@ build:
   string: cuda{{ cuda_major }}_{{ date_string }}_{{ GIT_DESCRIBE_HASH }}_{{ GIT_DESCRIBE_NUMBER }}
   ignore_run_exports:
     - libwebp-base
+  ignore_run_exports_from:
+    {% if cuda_major == "11" %}
+    - {{ compiler('cuda11') }}
+    {% endif %}
   script_env:
     - AWS_ACCESS_KEY_ID
     - AWS_SECRET_ACCESS_KEY
@@ -39,6 +43,12 @@ requirements:
   build:
     - {{ compiler("c") }}
     - {{ compiler("cxx") }}
+    {% if cuda_major == "11" %}
+    - {{ compiler('cuda11') }} ={{ cuda_version }}
+    {% else %}
+    - {{ compiler('cuda') }}
+    {% endif %}
+    - cuda-version ={{ cuda_version }}
     - binutils
     - cmake >=3.23.1,!=3.25.0
     - make
@@ -46,15 +56,30 @@ requirements:
     - sysroot_{{ target_platform }} {{ sysroot_version }}
     - yasm # [x86_64]
   host:
-    - cudatoolkit ={{ cuda_version }}
+    - cuda-version ={{ cuda_version }}
+    {% if cuda_major == "11" %}
+    - libcufile {{ cuda11_libcufile_host_version }}         # [linux64]
+    - libcufile-dev {{ cuda11_libcufile_host_version }}     # [linux64]
+    - libnvjpeg {{ cuda11_libnvjpeg_host_version }}
+    - libnvjpeg-dev {{ cuda11_libnvjpeg_host_version }}
+    {% else %}
+    - libcufile-dev     # [linux64]
+    - libnvjpeg-dev
+    - libnvjpeg-static
+    {% endif %}
     - jbig
-    - libwebp-base {{ libwebp_base_version }}  # [linux or osx]
+    - libwebp-base
+    - nvtx-c >=3.1.0
     - openslide
     - xz
     - zlib
     - zstd
   run:
-    - {{ pin_compatible('cudatoolkit', max_pin='x', min_pin='x') }}
+    - {{ pin_compatible('cuda-version', max_pin='x', min_pin='x') }}
+    {% if cuda_major == "11" %}
+    - cudatoolkit
+    - libcufile {{ cuda11_libcufile_run_version }}      # [linux64]
+    {% endif %}
     - {{ pin_compatible('libwebp-base', max_pin='x.x') }}
     - jbig
     # - openslide # skipping here but benchmark binary would need openslide library
diff --git a/cpp/plugins/cucim.kit.cumed/VERSION b/cpp/plugins/cucim.kit.cumed/VERSION
index 4f3f0be7c..e7ab23867 100644
--- a/cpp/plugins/cucim.kit.cumed/VERSION
+++ b/cpp/plugins/cucim.kit.cumed/VERSION
@@ -1,1 +1,1 @@
-23.06.00
+23.08.00
diff --git a/cpp/plugins/cucim.kit.cuslide/VERSION b/cpp/plugins/cucim.kit.cuslide/VERSION
index 4f3f0be7c..e7ab23867 100644
--- a/cpp/plugins/cucim.kit.cuslide/VERSION
+++ b/cpp/plugins/cucim.kit.cuslide/VERSION
@@ -1,1 +1,1 @@
-23.06.00
+23.08.00
diff --git a/cucim.code-workspace b/cucim.code-workspace
index eca71fd21..ed540b80b 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@23.06.00.so"
+					"CUCIM_TEST_PLUGIN_PATH": "cucim.kit.cuslide@23.08.00.so"
 				},
 				"cwd": "${workspaceDirectory}",
 				"catch2": {
@@ -226,7 +226,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@23.06.00.so"
+						"value": "cucim.kit.cuslide@23.08.00.so"
 					}
 				],
 				"console": "externalTerminal",
@@ -254,7 +254,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@23.06.00.so"
+						"value": "cucim.kit.cuslide@23.08.00.so"
 					}
 				],
 				"console": "externalTerminal",
@@ -286,7 +286,7 @@
 					},
 					{
 						"name": "CUCIM_TEST_PLUGIN_PATH",
-						"value": "cucim.kit.cuslide@23.06.00.so"
+						"value": "cucim.kit.cuslide@23.08.00.so"
 					}
 				],
 				"console": "externalTerminal",
diff --git a/dependencies.yaml b/dependencies.yaml
index d314b2e19..cdf60e9cf 100644
--- a/dependencies.yaml
+++ b/dependencies.yaml
@@ -3,7 +3,7 @@ files:
   all:
     output: conda
     matrix:
-      cuda: ["11.8"]
+      cuda: ["11.8", "12.0"]
       arch: [x86_64]
     includes:
       - build
@@ -66,12 +66,28 @@ dependencies:
               arch: x86_64
               cuda: "11.8"
             packages:
+              - cuda-version=11.8
               - nvcc_linux-64=11.8
+              - libcufile=1.4.0.31
+              - libcufile-dev=1.4.0.31
+              - libnvjpeg=11.6.0.55
+              - libnvjpeg-dev=11.6.0.55
           - matrix:
               arch: aarch64
               cuda: "11.8"
             packages:
+              - cuda-version=11.8
               - nvcc_linux-aarch64=11.8
+              - libnvjpeg=11.6.0.55
+              - libnvjpeg-dev=11.6.0.55
+          - matrix:
+              cuda: "12.0"
+            packages:
+              - cuda-version=12.0
+              - cuda-nvcc
+              - libcufile-dev
+              - libnvjpeg-dev
+              - libnvjpeg-static
   checks:
     common:
       - output_types: [conda, requirements]
@@ -84,19 +100,30 @@ dependencies:
           - matrix:
               cuda: "11.2"
             packages:
-              - cudatoolkit=11.2
+              - cuda-version=11.2
+              - cudatoolkit
           - matrix:
               cuda: "11.4"
             packages:
-              - cudatoolkit=11.4
+              - cuda-version=11.4
+              - cudatoolkit
           - matrix:
               cuda: "11.5"
             packages:
-              - cudatoolkit=11.5
+              - cuda-version=11.5
+              - cudatoolkit
           - matrix:
               cuda: "11.8"
             packages:
-              - cudatoolkit=11.8
+              - cuda-version=11.8
+              - cudatoolkit
+          - matrix:
+              cuda: "12.0"
+            packages:
+              - cuda-version=12.0
+              - cuda-cudart-dev
+              - libnvjpeg-dev
+              - libcufile-dev
   docs:
     common:
       - output_types: [conda, requirements]
@@ -135,9 +162,10 @@ dependencies:
           - click
           - cupy >=12.0.0
           - jbig
+          - lazy_loader >=0.1
           - libwebp-base
           - numpy >=1.21.3
-          - scikit-image >=0.19.0,<0.21.0a0
+          - scikit-image >=0.19.0,<0.22.0a0
           - scipy
           - xz
           - zlib
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 44423d8d7..3fea3380b 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -72,9 +72,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '23.06'
+version = '23.08'
 # The full version, including alpha/beta/rc tags.
-release = '23.06.00'
+release = '23.08.00'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/python/cucim/MANIFEST.in b/python/cucim/MANIFEST.in
index 4244cebec..3c2e31645 100644
--- a/python/cucim/MANIFEST.in
+++ b/python/cucim/MANIFEST.in
@@ -22,4 +22,6 @@ include tox.ini .travis.yml .appveyor.yml .readthedocs.yml
 include versioneer.py
 include src/cucim/clara/*.so*
 
+recursive-include src/cucim *.py *.pyi *.cu *.h *.npy *.txt *.md
+
 global-exclude *.py[cod] __pycache__/*
diff --git a/python/cucim/README.md b/python/cucim/README.md
index 081702200..1249e06dd 100755
--- a/python/cucim/README.md
+++ b/python/cucim/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-23.06/notebooks/Welcome.ipynb))
+Please check out our [Welcome](notebooks/Welcome.ipynb) notebook ([NBViewer](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-23.08/notebooks/Welcome.ipynb))
 
 #### Downloading sample images
 
diff --git a/python/cucim/VERSION b/python/cucim/VERSION
index 4f3f0be7c..e7ab23867 100644
--- a/python/cucim/VERSION
+++ b/python/cucim/VERSION
@@ -1,1 +1,1 @@
-23.06.00
+23.08.00
diff --git a/python/cucim/ci/requirements.txt b/python/cucim/ci/requirements.txt
index d7f5177e6..a0f53e12a 100644
--- a/python/cucim/ci/requirements.txt
+++ b/python/cucim/ci/requirements.txt
@@ -2,3 +2,4 @@ virtualenv>=16.6.0
 pip>=19.1.1
 setuptools>=18.0.1
 six>=1.14.0
+lazy_loader>=0.1
diff --git a/python/cucim/docs/getting_started/index.md b/python/cucim/docs/getting_started/index.md
index 09c345372..a0228a17e 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-v23.06.00-linux.tar.gz`).
+Please download the latest SDK package (`cuCIM-v23.08.00-linux.tar.gz`).
 
 Untar the downloaded file.
 
 ```bash
-mkdir -p cuCIM-v23.06.00
-tar -xzvf cuCIM-v23.06.00-linux.tar.gz -C cuCIM-v23.06.00
+mkdir -p cuCIM-v23.08.00
+tar -xzvf cuCIM-v23.08.00-linux.tar.gz -C cuCIM-v23.08.00
 
-cd cuCIM-v23.06.00
+cd cuCIM-v23.08.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@23.06.00.so
+[Plugin: cucim.kit.cuslide] Loading the dynamic library from: cucim.kit.cuslide@23.08.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 ea88d7af4..e5000f865 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 23.06.00](release_notes/v23.06.00.md).
+Current latest version is [Version 23.08.00](release_notes/v23.08.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 1b7179b6b..214a44f25 100755
--- a/python/cucim/setup.py
+++ b/python/cucim/setup.py
@@ -39,6 +39,7 @@ def read(*names, **kwargs):
     url='https://developer.nvidia.com/multidimensional-image-processing',
     packages=find_packages('src'),
     package_dir={'cucim': 'src/cucim'},
+    package_data={"": ["*.pyi", "*.h", "*.cu"]},
     include_package_data=True,
     zip_safe=False,
     classifiers=[
@@ -56,11 +57,10 @@ def read(*names, **kwargs):
         'Programming Language :: C++',
         'Programming Language :: Python',
         'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.6',
-        'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',
         'Programming Language :: Python :: 3.9',
         'Programming Language :: Python :: 3.10',
+        'Programming Language :: Python :: 3.11',
         # 'Operating System :: OS Independent',
         # 'Operating System :: Unix',
         # 'Operating System :: POSIX',
@@ -87,7 +87,7 @@ def read(*names, **kwargs):
     setup_requires=SETUP_REQUIRES,
     install_requires=[
         # TODO: Check cupy dependency based on cuda version
-        'click', 'numpy',  # 'scipy', 'scikit-image'
+        'click', 'numpy', "lazy_loader>=0.1",  # 'scipy', 'scikit-image'
         # eg: 'aspectlib==1.1.1', 'six>=1.7',
     ],
     extras_require={
diff --git a/python/cucim/src/cucim/__init__.py b/python/cucim/src/cucim/__init__.py
index d51c2cf2e..37bfa2cb3 100644
--- a/python/cucim/src/cucim/__init__.py
+++ b/python/cucim/src/cucim/__init__.py
@@ -59,13 +59,9 @@
     del _version
 
 
-from .skimage._shared import lazy
+import lazy_loader as lazy
 
-__getattr__, __lazy_dir__, _ = lazy.attach(
-    __name__,
-    submodules,
-    submod_attrs,
-)
+__getattr__, __lazy_dir__, _ = lazy.attach_stub(__name__, __file__)
 
 
 def __dir__():
diff --git a/python/cucim/src/cucim/__init__.pyi b/python/cucim/src/cucim/__init__.pyi
new file mode 100644
index 000000000..985551d7e
--- /dev/null
+++ b/python/cucim/src/cucim/__init__.pyi
@@ -0,0 +1,35 @@
+#
+# Copyright (c) 2020-2021, NVIDIA CORPORATION.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+submodules = []
+
+try:
+    import cupy
+    _is_cupy_available = True
+    submodules += ['core', 'skimage']
+except ImportError:
+    pass
+
+try:
+    from .clara import CuImage, __version__, cli
+    _is_clara_available = True
+    submodules += ['clara']
+except ImportError:
+    from ._version import get_versions
+    __version__ = get_versions()['version']
+    del get_versions
+    del _version
+
+__all__ = submodules + ['__version__', 'is_available']
diff --git a/python/cucim/src/cucim/skimage/__init__.py b/python/cucim/src/cucim/skimage/__init__.py
index b48ceace8..c6cf7d4eb 100644
--- a/python/cucim/src/cucim/skimage/__init__.py
+++ b/python/cucim/src/cucim/skimage/__init__.py
@@ -60,35 +60,17 @@
 
 """
 
-from ._shared import lazy
+import lazy_loader as lazy
 
-submodules = [
-    'color',
-    'data',
-    'exposure',
-    'feature',
-    'filters',
-    'measure',
-    'metrics',
-    'morphology',
-    'registration',
-    'restoration',
-    'segmentation',
-    'transform',
-    'util',
-]
-
-
-__getattr__, __lazy_dir__, _ = lazy.attach(
-    __name__,
-    submodules,
-    {'util.dtype': ['dtype_limits', 'img_as_bool', 'img_as_float',
-                    'img_as_float32', 'img_as_float64', 'img_as_int',
-                    'img_as_ubyte', 'img_as_uint'],
-     'util.lookfor': ['lookfor'],
-    }
-)
+__getattr__, __lazy_dir__, _ = lazy.attach_stub(__name__, __file__)
 
 
 def __dir__():
     return __lazy_dir__()
+
+
+# Legacy imports into the root namespace; not advertised in __all__
+from .util.dtype import (dtype_limits, img_as_bool, img_as_float,
+                         img_as_float32, img_as_float64, img_as_int,
+                         img_as_ubyte, img_as_uint)
+from .util.lookfor import lookfor
diff --git a/python/cucim/src/cucim/skimage/__init__.pyi b/python/cucim/src/cucim/skimage/__init__.pyi
new file mode 100644
index 000000000..1c3c07ea6
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/__init__.pyi
@@ -0,0 +1,23 @@
+import lazy_loader as lazy
+
+submodules = [
+    'color',
+    'data',
+    'exposure',
+    'feature',
+    'filters',
+    'measure',
+    'metrics',
+    'morphology',
+    'registration',
+    'restoration',
+    'segmentation',
+    'transform',
+    'util',
+]
+
+__all__ = submodules
+
+from . import (color, data, exposure, feature, filters, measure, metrics,
+               morphology, registration, restoration, segmentation, transform,
+               util)
diff --git a/python/cucim/src/cucim/skimage/_shared/lazy.py b/python/cucim/src/cucim/skimage/_shared/lazy.py
deleted file mode 100644
index 2b7270c66..000000000
--- a/python/cucim/src/cucim/skimage/_shared/lazy.py
+++ /dev/null
@@ -1,139 +0,0 @@
-import importlib
-import importlib.util
-import os
-import sys
-
-
-def attach(package_name, submodules=None, submod_attrs=None):
-    """Attach lazily loaded submodules, functions, or other attributes.
-
-    Typically, modules import submodules and attributes as follows::
-
-      import mysubmodule
-      import anothersubmodule
-
-      from .foo import someattr
-
-    The idea is to replace a package's `__getattr__`, `__dir__`, and
-    `__all__`, such that all imports work exactly the way they did
-    before, except that they are only imported when used.
-
-    The typical way to call this function, replacing the above imports, is::
-
-      __getattr__, __lazy_dir__, __all__ = lazy.attach(
-        __name__,
-        ['mysubmodule', 'anothersubmodule'],
-        {'foo': 'someattr'}
-      )
-
-    This functionality requires Python 3.7 or higher.
-
-    Parameters
-    ----------
-    package_name : str
-        Typically use ``__name__``.
-    submodules : set
-        List of submodules to attach.
-    submod_attrs : dict
-        Dictionary of submodule -> list of attributes / functions.
-        These attributes are imported as they are used.
-
-    Returns
-    -------
-    __getattr__, __dir__, __all__
-
-    """
-    if submod_attrs is None:
-        submod_attrs = {}
-
-    if submodules is None:
-        submodules = set()
-    else:
-        submodules = set(submodules)
-
-    attr_to_modules = {
-        attr: mod for mod, attrs in submod_attrs.items() for attr in attrs
-    }
-
-    __all__ = list(submodules | attr_to_modules.keys())
-
-    def __getattr__(name):
-        if name in submodules:
-            return importlib.import_module(f'{package_name}.{name}')
-        elif name in attr_to_modules:
-            submod = importlib.import_module(
-                f'{package_name}.{attr_to_modules[name]}'
-            )
-            return getattr(submod, name)
-        else:
-            raise AttributeError(f'No {package_name} attribute {name}')
-
-    def __dir__():
-        return __all__
-
-    eager_import = os.environ.get('EAGER_IMPORT', '')
-    if eager_import not in ['', '0', 'false']:
-        for attr in set(attr_to_modules.keys()) | submodules:
-            __getattr__(attr)
-
-    return __getattr__, __dir__, list(__all__)
-
-
-def load(fullname):
-    """Return a lazily imported proxy for a module.
-
-    We often see the following pattern::
-
-      def myfunc():
-          import scipy as sp
-          sp.argmin(...)
-          ....
-
-    This is to prevent a module, in this case `scipy`, from being
-    imported at function definition time, since that can be slow.
-
-    This function provides a proxy module that, upon access, imports
-    the actual module.  So the idiom equivalent to the above example is::
-
-      sp = lazy.load("scipy")
-
-      def myfunc():
-          sp.argmin(...)
-          ....
-
-    The initial import time is fast because the actual import is delayed
-    until the first attribute is requested. The overall import time may
-    decrease as well for users that don't make use of large portions
-    of the library.
-
-    Parameters
-    ----------
-    fullname : str
-        The full name of the module or submodule to import.  For example::
-
-          sp = lazy.load('scipy')  # import scipy as sp
-          spla = lazy.load('scipy.linalg')  # import scipy.linalg as spla
-
-    Returns
-    -------
-    pm : importlib.util._LazyModule
-        Proxy module.  Can be used like any regularly imported module.
-        Actual loading of the module occurs upon first attribute request.
-
-    """
-    try:
-        return sys.modules[fullname]
-    except KeyError:
-        pass
-
-    spec = importlib.util.find_spec(fullname)
-    if spec is None:
-        raise ModuleNotFoundError(f"No module name '{fullname}'")
-
-    module = importlib.util.module_from_spec(spec)
-    sys.modules[fullname] = module
-
-    loader = importlib.util.LazyLoader(spec.loader)
-    loader.exec_module(module)
-
-    return module
diff --git a/python/cucim/src/cucim/skimage/_shared/utils.py b/python/cucim/src/cucim/skimage/_shared/utils.py
index b4371ff1b..51b25e198 100644
--- a/python/cucim/src/cucim/skimage/_shared/utils.py
+++ b/python/cucim/src/cucim/skimage/_shared/utils.py
@@ -9,17 +9,48 @@
 
 from ._warnings import all_warnings, warn  # noqa
 
+__all__ = ['deprecate_func', 'get_bound_method_class', 'all_warnings',
+           'safe_as_int', 'check_shape_equality', 'check_nD', 'warn',
+           'reshape_nd', 'identity', 'slice_at_axis']
 
-class skimage_deprecation(Warning):
-    """Create our own deprecation class, since Python >= 2.7
-    silences deprecations by default.
 
+def _get_stack_rank(func):
+    """Return function rank in the call stack."""
+    if _is_wrapped(func):
+        return 1 + _get_stack_rank(func.__wrapped__)
+    else:
+        return 0
+
+
+def _is_wrapped(func):
+    return "__wrapped__" in dir(func)
+
+
+def _get_stack_length(func):
+    """Return function call stack length."""
+    return _get_stack_rank(func.__globals__.get(func.__name__, func))
+
+
+class _DecoratorBaseClass:
+    """Used to manage decorators' warnings stacklevel.
+
+    The `_stack_length` class variable is used to store the number of
+    times a function is wrapped by a decorator.
+
+    Let `stack_length` be the total number of times a decorated
+    function is wrapped, and `stack_rank` be the rank of the decorator
+    in the decorators stack. The stacklevel of a warning is then
+    `stacklevel = 1 + stack_length - stack_rank`.
     """
 
-    pass
+    _stack_length = {}
+
+    def get_stack_length(self, func):
+        return self._stack_length.get(func.__name__,
+                                      _get_stack_length(func))
 
 
-class change_default_value:
+class change_default_value(_DecoratorBaseClass):
     """Decorator for changing the default value of an argument.
 
     Parameters
@@ -48,6 +79,7 @@ def __call__(self, func):
         arg_idx = list(parameters.keys()).index(self.arg_name)
         old_value = parameters[self.arg_name].default
 
+        stack_rank = _get_stack_rank(func)
         if self.warning_msg is None:
             self.warning_msg = (
                 f"The new recommended value for {self.arg_name} is "
@@ -59,15 +91,17 @@ def __call__(self, func):
 
         @functools.wraps(func)
         def fixed_func(*args, **kwargs):
+            stacklevel = 1 + self.get_stack_length(func) - stack_rank
             if len(args) < arg_idx + 1 and self.arg_name not in kwargs.keys():
                 # warn that arg_name default value changed:
-                warnings.warn(self.warning_msg, FutureWarning, stacklevel=2)
+                warnings.warn(self.warning_msg, FutureWarning,
+                              stacklevel=stacklevel)
             return func(*args, **kwargs)
 
         return fixed_func
 
 
-class remove_arg:
+class remove_arg(_DecoratorBaseClass):
     """Decorator to remove an argument from function's signature.
 
     Parameters
@@ -100,21 +134,25 @@ def __call__(self, func):
         if self.help_msg is not None:
             warning_msg += f" {self.help_msg}"
 
+        stack_rank = _get_stack_rank(func)
+
         @functools.wraps(func)
         def fixed_func(*args, **kwargs):
+            stacklevel = 1 + self.get_stack_length(func) - stack_rank
             if len(args) > arg_idx or self.arg_name in kwargs.keys():
                 # warn that arg_name is deprecated
-                warnings.warn(warning_msg, FutureWarning, stacklevel=2)
+                warnings.warn(warning_msg, FutureWarning,
+                              stacklevel=stacklevel)
             return func(*args, **kwargs)
 
         return fixed_func
 
 
-def docstring_add_deprecated(func, kwarg_mapping, deprecated_version):
+def _docstring_add_deprecated(func, kwarg_mapping, deprecated_version):
     """Add deprecated kwarg(s) to the "Other Params" section of a docstring.
 
     Parameters
-    ---------
+    ----------
     func : function
         The function whose docstring we wish to update.
     kwarg_mapping : dict
@@ -176,7 +214,7 @@ def docstring_add_deprecated(func, kwarg_mapping, deprecated_version):
     return final_docstring
 
 
-class deprecate_kwarg:
+class deprecate_kwarg(_DecoratorBaseClass):
     """Decorator ensuring backward compatibility when argument names are
     modified in a function definition.
 
@@ -213,14 +251,18 @@ def __init__(self, kwarg_mapping, deprecated_version, warning_msg=None,
 
     def __call__(self, func):
 
+        stack_rank = _get_stack_rank(func)
+
         @functools.wraps(func)
         def fixed_func(*args, **kwargs):
+            stacklevel = 1 + self.get_stack_length(func) - stack_rank
             for old_arg, new_arg in self.kwarg_mapping.items():
                 if old_arg in kwargs:
                     #  warn that the function interface has changed:
                     warnings.warn(self.warning_msg.format(
                         old_arg=old_arg, func_name=func.__name__,
-                        new_arg=new_arg), FutureWarning, stacklevel=2)
+                        new_arg=new_arg), FutureWarning,
+                        stacklevel=stacklevel)
                     # Substitute new_arg to old_arg
                     kwargs[new_arg] = kwargs.pop(old_arg)
 
@@ -228,13 +270,13 @@ def fixed_func(*args, **kwargs):
             return func(*args, **kwargs)
 
         if func.__doc__ is not None:
-            newdoc = docstring_add_deprecated(func, self.kwarg_mapping,
-                                              self.deprecated_version)
+            newdoc = _docstring_add_deprecated(func, self.kwarg_mapping,
+                                               self.deprecated_version)
             fixed_func.__doc__ = newdoc
         return fixed_func
 
 
-class channel_as_last_axis():
+class channel_as_last_axis:
     """Decorator for automatically making channels axis last for all arrays.
 
     This decorator reorders axes for compatibility with functions that only
@@ -311,53 +353,66 @@ def fixed_func(*args, **kwargs):
         return fixed_func
 
 
-class deprecated(object):
-    """Decorator to mark deprecated functions with warning.
+class deprecate_func(_DecoratorBaseClass):
+    """Decorate a deprecated function and warn when it is called.
 
     Adapted from <http://wiki.python.org/moin/PythonDecoratorLibrary>.
 
     Parameters
     ----------
-    alt_func : str
-        If given, tell user what function to use instead.
-    behavior : {'warn', 'raise'}
-        Behavior during call to deprecated function: 'warn' = warn user that
-        function is deprecated; 'raise' = raise error.
+    deprecated_version : str
+        The package version when the deprecation was introduced.
     removed_version : str
         The package version in which the deprecated function will be removed.
+    hint : str, optional
+        A hint on how to address this deprecation,
+        e.g., "Use `skimage.submodule.alternative_func` instead."
+    Examples
+    --------
+    >>> @deprecate_func(
+    ...     deprecated_version="1.0.0",
+    ...     removed_version="1.2.0",
+    ...     hint="Use `bar` instead."
+    ... )
+    ... def foo():
+    ...     pass
+    Calling ``foo`` will warn with::
+        FutureWarning: `foo` is deprecated since version 1.0.0
+        and will be removed in version 1.2.0. Use `bar` instead.
     """
 
-    def __init__(self, alt_func=None, behavior='warn', removed_version=None):
-        self.alt_func = alt_func
-        self.behavior = behavior
+    def __init__(self, *, deprecated_version, removed_version=None, hint=None):
+        self.deprecated_version = deprecated_version
         self.removed_version = removed_version
+        self.hint = hint
 
     def __call__(self, func):
 
-        alt_msg = ''
-        if self.alt_func is not None:
-            alt_msg = f' Use ``{self.alt_func}`` instead.'
-        rmv_msg = ''
-        if self.removed_version is not None:
-            rmv_msg = f' and will be removed in version {self.removed_version}'
+        message = (
+            f"`{func.__name__}` is deprecated since version "
+            f"{self.deprecated_version}"
+        )
+        if self.removed_version:
+            message += (
+                f" and will be removed in version {self.removed_version}."
+            )
+        if self.hint:
+            message += f" {self.hint.rstrip('.')}."
 
-        msg = f'Function ``{func.__name__}`` is deprecated{rmv_msg}.{alt_msg}'
+        stack_rank = _get_stack_rank(func)
 
         @functools.wraps(func)
         def wrapped(*args, **kwargs):
-            if self.behavior == 'warn':
-                func_code = func.__code__
-                warnings.simplefilter('always', skimage_deprecation)
-                warnings.warn_explicit(msg,
-                                       category=skimage_deprecation,
-                                       filename=func_code.co_filename,
-                                       lineno=func_code.co_firstlineno + 1)
-            elif self.behavior == 'raise':
-                raise skimage_deprecation(msg)
+            stacklevel = 1 + self.get_stack_length(func) - stack_rank
+            warnings.warn(
+                message,
+                category=FutureWarning,
+                stacklevel=stacklevel
+            )
             return func(*args, **kwargs)
 
         # modify doc string to display deprecation warning
-        doc = '**Deprecated function**.' + alt_msg
+        doc = f'**Deprecated:** {message}'
         if wrapped.__doc__ is None:
             wrapped.__doc__ = doc
         else:
diff --git a/python/cucim/src/cucim/skimage/color/__init__.py b/python/cucim/src/cucim/skimage/color/__init__.py
index 390668712..ef9049f34 100644
--- a/python/cucim/src/cucim/skimage/color/__init__.py
+++ b/python/cucim/src/cucim/skimage/color/__init__.py
@@ -9,12 +9,14 @@
                         rgb_from_bex, rgb_from_bpx, rgb_from_bro, rgb_from_fgx,
                         rgb_from_gdx, rgb_from_hax, rgb_from_hdx, rgb_from_hed,
                         rgb_from_hpx, rgb_from_rbd, rgba2rgb, rgbcie2rgb,
-                        separate_stains, xyz2lab, xyz2luv, xyz2rgb, ycbcr2rgb,
-                        ydbdr2rgb, yiq2rgb, ypbpr2rgb, yuv2rgb)
+                        separate_stains, xyz2lab, xyz2luv, xyz2rgb,
+                        xyz_tristimulus_values, ycbcr2rgb, ydbdr2rgb, yiq2rgb,
+                        ypbpr2rgb, yuv2rgb)
 from .colorlabel import color_dict, label2rgb
 from .delta_e import deltaE_cie76, deltaE_ciede94, deltaE_ciede2000, deltaE_cmc
 
 __all__ = ['convert_colorspace',
+           'xyz_tristimulus_values',
            'rgba2rgb',
            'rgb2hsv',
            'hsv2rgb',
diff --git a/python/cucim/src/cucim/skimage/color/colorconv.py b/python/cucim/src/cucim/skimage/color/colorconv.py
index ece630492..f7a1e3678 100644
--- a/python/cucim/src/cucim/skimage/color/colorconv.py
+++ b/python/cucim/src/cucim/skimage/color/colorconv.py
@@ -56,7 +56,7 @@
 from scipy import linalg
 
 from .._shared.utils import (_supported_float_type, channel_as_last_axis,
-                             identity)
+                             deprecate_func, identity)
 from ..util import dtype, dtype_limits
 
 
@@ -546,7 +546,7 @@ def hsv2rgb(hsv, *, channel_axis=-1):
 #    ----------
 #    .. [1] https://en.wikipedia.org/wiki/Standard_illuminant
 
-illuminants = \
+_illuminants = \
     {"A": {'2': (1.098466069456375, 1, 0.3558228003436005),
            '10': (1.111420406956693, 1, 0.3519978321919493),
            'R': (1.098466069456375, 1, 0.3558228003436005)},
@@ -573,7 +573,83 @@ def hsv2rgb(hsv, *, channel_axis=-1):
            'R': (1.0, 1.0, 1.0)}}
 
 
-def get_xyz_coords(illuminant, observer):
+def xyz_tristimulus_values(*, illuminant, observer, dtype=None):
+    """Get the CIE XYZ tristimulus values.
+
+    Given an illuminant and observer, this function returns the CIE XYZ
+    tristimulus values [2]_ scaled such that :math:`Y = 1`.
+
+    Parameters
+    ----------
+    illuminant : {"A", "B", "C", "D50", "D55", "D65", "D75", "E"}
+        The name of the illuminant (the function is NOT case sensitive).
+    observer : {"2", "10", "R"}
+        One of: 2-degree observer, 10-degree observer, or 'R' observer as in
+        R function ``grDevices::convertColor`` [3]_.
+    dtype : np.dtype, optional
+        This argument is ignored in the cuCIM implementation of
+        `xyz_tristimulus_values` since an array is not returned. The output is
+        always a 3-tuple of float.
+
+    Returns
+    -------
+    values : 3-tuple of float
+        Three elements :math:`X, Y, Z` containing the CIE XYZ tristimulus values
+        of the given illuminant.
+
+    Raises
+    ------
+    ValueError
+        If either the illuminant or the observer angle are not supported or
+        unknown.
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Standard_illuminant#White_points_of_standard_illuminants
+    .. [2] https://en.wikipedia.org/wiki/CIE_1931_color_space#Meaning_of_X,_Y_and_Z
+    .. [3] https://www.rdocumentation.org/packages/grDevices/versions/3.6.2/topics/convertColor
+
+    Notes
+    -----
+    The return type of this function differs from the one in scikit-image as it
+    always returns a 3-tuple of float rather than an array with a
+    user-specified dtype.
+
+    The CIE XYZ tristimulus values are calculated from :math:`x, y` [1]_, using the
+    formula
+
+    .. math:: X = x / y
+
+    .. math:: Y = 1
+
+    .. math:: Z = (1 - x - y) / y
+
+    The only exception is the illuminant "D65" with aperture angle 2° for
+    backward-compatibility reasons.
+
+    Examples
+    --------
+    Get the CIE XYZ tristimulus values for a "D65" illuminant for a 10 degree
+    field of view
+
+    >>> xyz_tristimulus_values(illuminant="D65", observer="10")
+    array([0.94809668, 1.        , 1.07305136])
+    """  # noqa
+    illuminant = illuminant.upper()
+    observer = observer.upper()
+    try:
+        return _illuminants[illuminant][observer]
+    except KeyError:
+        raise ValueError(f'Unknown illuminant/observer combination '
+                         f'(`{illuminant}`, `{observer}`)')
+
+
+@deprecate_func(
+    hint="Use `skimage.color.xyz_tristimulus_values` instead.",
+    deprecated_version="23.08",
+    removed_version="24.06",
+)
+def get_xyz_coords(illuminant, observer, dtype=float):
     """Get the XYZ coordinates of the given illuminant and observer [1]_.
 
     Parameters
@@ -602,13 +678,7 @@ def get_xyz_coords(illuminant, observer):
     ----------
     .. [1] https://en.wikipedia.org/wiki/Standard_illuminant
     """
-    illuminant = illuminant.upper()
-    observer = observer.upper()
-    try:
-        return illuminants[illuminant][observer]
-    except KeyError:
-        raise ValueError(f'Unknown illuminant/observer combination '
-                         f'(`{illuminant}`, `{observer}`)')
+    return xyz_tristimulus_values(illuminant=illuminant, observer=observer)
 
 
 # Haematoxylin-Eosin-DAB colorspace
@@ -625,7 +695,6 @@ def get_xyz_coords(illuminant, observer):
                          [0.27, 0.57, 0.78]])
 hed_from_rgb = linalg.inv(rgb_from_hed)
 
-
 # Following matrices are adapted form the Java code written by G.Landini.
 # The original code is available at:
 # https://web.archive.org/web/20160624145052/http://www.mecourse.com/landinig/software/cdeconv/cdeconv.html
@@ -1071,9 +1140,12 @@ def gray2rgb(image, *, channel_axis=-1):
 @cp.memoize(for_each_device=True)
 def _get_xyz_to_lab_kernel(xyz_ref_white, name='xyz2lab'):
     _xyz_to_lab = f"""
+        // scale by CIE XYZ tristimulus values of the reference white point
         arr[3*i] /= {xyz_ref_white[0]};
         arr[3*i + 1] /= {xyz_ref_white[1]};
         arr[3*i + 2] /= {xyz_ref_white[2]};
+
+        // Nonlinear distortion and linear transformation
         for (int ch=0; ch < 3; ch++)
         {{
             if (arr[3*i + ch] > 0.008856)
@@ -1083,6 +1155,8 @@ def _get_xyz_to_lab_kernel(xyz_ref_white, name='xyz2lab'):
                 arr[3*i + ch] = 7.787 * arr[3*i + ch] + 16.0 / 116.0;
             }}
         }}
+
+        // Vector scaling
         lab[3*i] = (116. * arr[3*i + 1]) - 16.0;
         lab[3*i + 1] = 500.0 * (arr[3*i] - arr[3*i + 1]);
         lab[3*i + 2] = 200.0 * (arr[3*i + 1] - arr[3*i + 2]);
@@ -1130,13 +1204,13 @@ def xyz2lab(xyz, 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.
+    x_ref=95.047, y_ref=100., z_ref=108.883. See function
+    :func:`~.xyz_tristimulus_values` 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
+    .. [1] http://www.easyrgb.com/en/math.php
+    .. [2] https://en.wikipedia.org/wiki/CIELAB_color_space
 
     Examples
     --------
@@ -1150,8 +1224,9 @@ def xyz2lab(xyz, illuminant="D65", observer="2", *, channel_axis=-1):
     xyz = _prepare_colorarray(xyz, force_copy=True, force_c_contiguous=True,
                               channel_axis=channel_axis)
 
-    xyz_ref_white = get_xyz_coords(illuminant, observer)
-
+    xyz_ref_white = xyz_tristimulus_values(
+        illuminant=illuminant, observer=observer
+    )
     name = f'xyz2lab_{xyz.dtype.char}'
     kern = _get_xyz_to_lab_kernel(xyz_ref_white, name=name)
     lab = cp.empty_like(xyz)
@@ -1222,7 +1297,7 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     Returns
     -------
     out : (..., 3, ...) ndarray
-        The image in XYZ color space. Same dimensions as input.
+        The image in XYZ color space, of same shape as input.
 
     Raises
     ------
@@ -1237,18 +1312,49 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     Notes
     -----
     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
+    z_ref = 108.883. See function :func:`~.xyz_tristimulus_values` for a list of
     supported illuminants.
 
+    See Also
+    --------
+    xyz2lab
+
     References
     ----------
-    .. [1] http://www.easyrgb.com/index.php?X=MATH&H=07
+    .. [1] http://www.easyrgb.com/en/math.php
     .. [2] https://en.wikipedia.org/wiki/CIELAB_color_space
     """
+    xyz, n_invalid = _lab2xyz(lab, illuminant, observer, channel_axis)
+    if n_invalid > 0:
+        warn(
+            "Conversion from CIE-LAB to XYZ color space resulted in "
+            f"{n_invalid} negative Z values that have been clipped to zero",
+            stacklevel=3,
+        )
+
+    return xyz
+
+
+def _lab2xyz(lab, illuminant, observer, channel_axis):
+    """Convert CIE-LAB to XYZ color space.
+
+    Internal function for :func:`~.lab2xyz` and others. In addition to the
+    converted image, return the number of invalid pixels in the Z channel for
+    correct warning propagation.
+
+    Returns
+    -------
+    out : (..., 3, ...) ndarray
+        The image in XYZ format. Same dimensions as input.
+    n_invalid : int
+        Number of invalid pixels in the Z channel after conversion.
+    """
     lab = _prepare_colorarray(lab, force_c_contiguous=True,
                               channel_axis=channel_axis)
 
-    xyz_ref_white = get_xyz_coords(illuminant, observer)
+    xyz_ref_white = xyz_tristimulus_values(
+        illuminant=illuminant, observer=observer
+    )
 
     name = f'lab2xyz_{lab.dtype.char}'
     kern = _get_lab_to_xyz_kernel(xyz_ref_white, name=name)
@@ -1258,12 +1364,8 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     #       operations?
     warnings = cp.zeros(lab.shape[:-1], dtype=np.int32)
     kern(lab, xyz, warnings, size=lab.size // 3)
-
-    nwarn = int(cp.count_nonzero(warnings))
-    if nwarn > 0:  # synchronize!
-        warn(f'Color data out of range: Z < 0 in {nwarn} pixels',
-             stacklevel=3)
-    return xyz
+    n_invalid = int(cp.count_nonzero(warnings))  # synchronize!
+    return xyz, n_invalid
 
 
 @channel_as_last_axis()
@@ -1302,8 +1404,8 @@ def rgb2lab(rgb, illuminant="D65", observer="2", *, channel_axis=-1):
 
     This function uses rgb2xyz and xyz2lab.
     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.
+    x_ref=95.047, y_ref=100., z_ref=108.883. See function
+    :func:`~.xyz_tristimulus_values` for a list of supported illuminants.
 
     References
     ----------
@@ -1335,7 +1437,7 @@ def lab2rgb(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     Returns
     -------
     out : (..., 3, ...) ndarray
-        The image in RGB format. Same dimensions as input.
+        The image in sRGB color space, of same shape as input.
 
     Raises
     ------
@@ -1344,17 +1446,28 @@ def lab2rgb(lab, illuminant="D65", observer="2", *, channel_axis=-1):
 
     Notes
     -----
-    This function uses lab2xyz and xyz2rgb.
+    This function uses :func:`~.lab2xyz` and :func:`~.xyz2rgb`.
     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
+    z_ref = 108.883. See function :func:`~.xyz_tristimulus_values` for a list of
     supported illuminants.
 
+    See Also
+    --------
+    rgb2lab
+
     References
     ----------
     .. [1] https://en.wikipedia.org/wiki/Standard_illuminant
     .. [2] https://en.wikipedia.org/wiki/CIELAB_color_space
     """
-    return xyz2rgb(lab2xyz(lab, illuminant, observer))
+    xyz, n_invalid = _lab2xyz(lab, illuminant, observer, channel_axis)
+    if n_invalid != 0:
+        warn(
+            "Conversion from CIE-LAB, via XYZ to sRGB color space resulted in "
+            f"{n_invalid} negative Z values that have been clipped to zero",
+            stacklevel=3,
+        )
+    return xyz2rgb(xyz, channel_axis=channel_axis)
 
 
 @cp.memoize(for_each_device=True)
@@ -1442,12 +1555,12 @@ def xyz2luv(xyz, illuminant="D65", observer="2", *, channel_axis=-1):
     -----
     By default XYZ conversion weights use observer=2A. Reference whitepoint
     for D65 Illuminant, with XYZ tristimulus values of ``(95.047, 100.,
-    108.883)``. See function 'get_xyz_coords' for a list of supported
-    illuminants.
+    108.883)``. See function :func:`~.xyz_tristimulus_values` for a list of
+    supported illuminants.
 
     References
     ----------
-    .. [1] http://www.easyrgb.com/index.php?X=MATH&H=16#text16
+    .. [1] http://www.easyrgb.com/en/math.php
     .. [2] https://en.wikipedia.org/wiki/CIELUV
 
     Examples
@@ -1466,7 +1579,9 @@ def xyz2luv(xyz, illuminant="D65", observer="2", *, channel_axis=-1):
     xyz = _prepare_colorarray(xyz, force_c_contiguous=True,
                               channel_axis=channel_axis)
 
-    xyz_ref_white = get_xyz_coords(illuminant, observer)
+    xyz_ref_white = xyz_tristimulus_values(
+        illuminant=illuminant, observer=observer
+    )
     kern = _get_xyz_to_luv_kernel(xyz_ref_white, xyz.dtype)
     luv = cp.empty_like(xyz)
     kern(xyz, luv, size=xyz.size // 3)
@@ -1546,16 +1661,19 @@ def luv2xyz(luv, illuminant="D65", observer="2", *, channel_axis=-1):
     -----
     XYZ conversion weights use observer=2A. Reference whitepoint for D65
     Illuminant, with XYZ tristimulus values of ``(95.047, 100., 108.883)``. See
-    function 'get_xyz_coords' for a list of supported illuminants.
+    function :func:`~.xyz_tristimulus_values` for a list of supported
+    illuminants.
 
     References
     ----------
-    .. [1] http://www.easyrgb.com/index.php?X=MATH&H=16#text16
+    .. [1] http://www.easyrgb.com/en/math.php
     .. [2] https://en.wikipedia.org/wiki/CIELUV
     """
     luv = _prepare_colorarray(luv, force_c_contiguous=True,
                               channel_axis=channel_axis)
-    xyz_ref_white = get_xyz_coords(illuminant, observer)
+    xyz_ref_white = xyz_tristimulus_values(
+        illuminant=illuminant, observer=observer
+    )
     kern = _get_luv_to_xyz_kernel(xyz_ref_white, luv.dtype)
     xyz = cp.empty_like(luv)
     kern(luv, xyz, size=luv.size // 3)
@@ -1591,9 +1709,8 @@ def rgb2luv(rgb, *, channel_axis=-1):
 
     References
     ----------
-    .. [1] http://www.easyrgb.com/index.php?X=MATH&H=16#text16
-    .. [2] http://www.easyrgb.com/index.php?X=MATH&H=02#text2
-    .. [3] https://en.wikipedia.org/wiki/CIELUV
+    .. [1] http://www.easyrgb.com/en/math.php
+    .. [2] https://en.wikipedia.org/wiki/CIELUV
     """
     return xyz2luv(rgb2xyz(rgb))
 
diff --git a/python/cucim/src/cucim/skimage/color/delta_e.py b/python/cucim/src/cucim/skimage/color/delta_e.py
index 14ed971db..bf07162a8 100644
--- a/python/cucim/src/cucim/skimage/color/delta_e.py
+++ b/python/cucim/src/cucim/skimage/color/delta_e.py
@@ -29,7 +29,7 @@
 
 def _float_inputs(lab1, lab2, allow_float32=True):
     if allow_float32:
-        float_dtype = _supported_float_type([lab1.dtype, lab2.dtype])
+        float_dtype = _supported_float_type((lab1.dtype, lab2.dtype))
     else:
         float_dtype = cp.float64
     lab1 = lab1.astype(float_dtype, copy=False)
@@ -375,7 +375,7 @@ def get_dH2(lab1, lab2, *, channel_axis=-1):
     """
     # This function needs double precision internally for accuracy
     input_is_float_32 = _supported_float_type(
-        [lab1.dtype, lab2.dtype]
+        (lab1.dtype, lab2.dtype)
     ) == cp.float32
     lab1, lab2 = _float_inputs(lab1, lab2, allow_float32=False)
     a1, b1 = cp.moveaxis(lab1, source=channel_axis, destination=0)[1:3]
diff --git a/python/cucim/src/cucim/skimage/data/__init__.py b/python/cucim/src/cucim/skimage/data/__init__.py
index 00c7daddc..62d86f776 100644
--- a/python/cucim/src/cucim/skimage/data/__init__.py
+++ b/python/cucim/src/cucim/skimage/data/__init__.py
@@ -1,9 +1,3 @@
-from .._shared import lazy
+import lazy_loader as lazy
 
-__getattr__, __dir__, __all__ = lazy.attach(
-    __name__,
-    submodules={},
-    submod_attrs={
-        '_binary_blobs': ['binary_blobs'],
-    }
-)
+__getattr__, __dir__, __all__ = lazy.attach_stub(__name__, __file__)
diff --git a/python/cucim/src/cucim/skimage/data/__init__.pyi b/python/cucim/src/cucim/skimage/data/__init__.pyi
new file mode 100644
index 000000000..88deb1e76
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/data/__init__.pyi
@@ -0,0 +1,5 @@
+__all__ = [
+    'binary_blobs',
+]
+
+from ._binary_blobs import binary_blobs
diff --git a/python/cucim/src/cucim/skimage/exposure/_adapthist.py b/python/cucim/src/cucim/skimage/exposure/_adapthist.py
index 4d42e2ae3..939601973 100644
--- a/python/cucim/src/cucim/skimage/exposure/_adapthist.py
+++ b/python/cucim/src/cucim/skimage/exposure/_adapthist.py
@@ -4,14 +4,7 @@
 
 http://tog.acm.org/resources/GraphicsGems/
 
-The Graphics Gems code is copyright-protected.  In other words, you cannot
-claim the text of the code as your own and resell it. Using the code is
-permitted in any program, product, or library, non-commercial or commercial.
-Giving credit is not required, though is a nice gesture.  The code comes as-is,
-and if there are any flaws or problems with any Gems code, nobody involved with
-Gems - authors, editors, publishers, or webmasters - are to be held
-responsible.  Basically, don't be a jerk, and remember that anything free
-comes with no guarantee.
+Relicensed with permission of the author under the Modified BSD license.
 """
 import functools
 import itertools
@@ -97,7 +90,7 @@ def equalize_adapthist(image, kernel_size=None,
     elif isinstance(kernel_size, numbers.Number):
         kernel_size = (kernel_size,) * image.ndim
     elif len(kernel_size) != image.ndim:
-        ValueError(f'Incorrect value of `kernel_size`: {kernel_size}')
+        raise ValueError(f'Incorrect value of `kernel_size`: {kernel_size}')
 
     kernel_size = [int(k) for k in kernel_size]
 
diff --git a/python/cucim/src/cucim/skimage/exposure/exposure.py b/python/cucim/src/cucim/skimage/exposure/exposure.py
index 4230f36c9..c6038cae4 100644
--- a/python/cucim/src/cucim/skimage/exposure/exposure.py
+++ b/python/cucim/src/cucim/skimage/exposure/exposure.py
@@ -72,8 +72,14 @@ def _bincount_histogram(image, source_range, bin_centers=None):
     """
     if bin_centers is None:
         bin_centers = _bincount_histogram_centers(image, source_range)
-    image_min, image_max = bin_centers[0], bin_centers[-1]
-    image = _offset_array(image, image_min.item(), image_max.item())  # synchronize  # noqa
+    image_min, image_max = bin_centers[0].item(), bin_centers[-1].item()
+    image = _offset_array(image, image_min, image_max)  # synchronize  # noqa
+
+    # Casting back to unsigned dtype seems necessary to avoid incorrect
+    # results for larger integer ranges with CUDA 12.x.
+    unsigned_dtype_char = image.dtype.char.upper()
+    image = image.astype(unsigned_dtype_char, copy=False)
+
     hist = cp.bincount(
         image.ravel(), minlength=image_max - min(image_min, 0) + 1
     )
@@ -178,7 +184,9 @@ def _get_numpy_hist_range(image, source_range):
     elif source_range == 'dtype':
         hist_range = dtype_limits(image, clip_negative=False)
     else:
-        ValueError('Wrong value for the `source_range` argument')
+        raise ValueError(
+            f'Incorrect value for `source_range` argument: {source_range}'
+        )
     return hist_range
 
 
diff --git a/python/cucim/src/cucim/skimage/feature/__init__.py b/python/cucim/src/cucim/skimage/feature/__init__.py
index d26dfb347..fcfa62c03 100644
--- a/python/cucim/src/cucim/skimage/feature/__init__.py
+++ b/python/cucim/src/cucim/skimage/feature/__init__.py
@@ -1,4 +1,3 @@
-from .._shared.utils import deprecated
 from ._basic_features import multiscale_basic_features
 from ._canny import canny
 from ._daisy import daisy
diff --git a/python/cucim/src/cucim/skimage/feature/_canny.py b/python/cucim/src/cucim/skimage/feature/_canny.py
index ae61cb603..1d60dbf86 100644
--- a/python/cucim/src/cucim/skimage/feature/_canny.py
+++ b/python/cucim/src/cucim/skimage/feature/_canny.py
@@ -235,6 +235,13 @@ def _nonmaximum_suppression_bilinear(
     kernel = _get_nonmax_kernel(large_int=large_int)
 
     out = cp.empty_like(magnitude)
+
+    if isinstance(low_threshold, cp.ndarray):
+        # if array scalar was provided, make sure dtype matches other arrays
+        if low_threshold.ndim > 0:
+            raise ValueError("expected scalar low_treshold")
+        low_threshold = float(low_threshold)
+
     kernel(isobel, jsobel, magnitude, eroded_mask, low_threshold, out)
     return out
 
@@ -351,6 +358,9 @@ def canny(image, sigma=1., low_threshold=None, high_threshold=None, mask=None,
     # mask by one and then mask the output. We also mask out the border points
     # because who knows what lies beyond the edge of the image?
 
+    if (image.dtype.kind in 'iu' and image.dtype.itemsize >= 8):
+        raise ValueError("64-bit or larger integer images are not supported")
+
     check_nD(image, 2)
     dtype_max = dtype_limits(image, clip_negative=False)[1]
 
diff --git a/python/cucim/src/cucim/skimage/filters/__init__.py b/python/cucim/src/cucim/skimage/filters/__init__.py
index 29ebabac1..62d86f776 100644
--- a/python/cucim/src/cucim/skimage/filters/__init__.py
+++ b/python/cucim/src/cucim/skimage/filters/__init__.py
@@ -1,30 +1,3 @@
-from .._shared import lazy
+import lazy_loader as lazy
 
-__getattr__, __dir__, __all__ = lazy.attach(
-    __name__,
-    # submodules={'rank'},
-    submod_attrs={
-        'lpi_filter': ['filter_inverse', 'wiener', 'LPIFilter2D'],
-        '_gaussian': ['gaussian', 'difference_of_gaussians'],
-        'edges': ['sobel', 'sobel_h', 'sobel_v',
-                  'scharr', 'scharr_h', 'scharr_v',
-                  'prewitt', 'prewitt_h', 'prewitt_v',
-                  'roberts', 'roberts_pos_diag', 'roberts_neg_diag',
-                  'laplace',
-                  'farid', 'farid_h', 'farid_v'],
-        '_rank_order': ['rank_order'],
-        '_gabor': ['gabor_kernel', 'gabor'],
-        'thresholding': ['threshold_local', 'threshold_otsu', 'threshold_yen',
-                         'threshold_isodata', 'threshold_li', 'threshold_minimum',
-                         'threshold_mean', 'threshold_triangle',
-                         'threshold_niblack', 'threshold_sauvola',
-                         'threshold_multiotsu', 'try_all_threshold',
-                         'apply_hysteresis_threshold'],
-        'ridges': ['meijering', 'sato', 'frangi', 'hessian'],
-        '_median': ['median'],
-        '_sparse': ['correlate_sparse'],
-        '_unsharp_mask': ['unsharp_mask'],
-        '_window': ['window'],
-        '_fft_based': ['butterworth']
-    }
-)
+__getattr__, __dir__, __all__ = lazy.attach_stub(__name__, __file__)
diff --git a/python/cucim/src/cucim/skimage/filters/__init__.pyi b/python/cucim/src/cucim/skimage/filters/__init__.pyi
new file mode 100644
index 000000000..c7703e933
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/filters/__init__.pyi
@@ -0,0 +1,68 @@
+__all__ = [
+'LPIFilter2D',
+'apply_hysteresis_threshold',
+'butterworth',
+'correlate_sparse',
+'difference_of_gaussians',
+'farid',
+'farid_h',
+'farid_v',
+'filter_forward',
+'filter_inverse',
+'frangi',
+'gabor',
+'gabor_kernel',
+'gaussian',
+'hessian',
+'laplace',
+'median',
+'meijering',
+'prewitt',
+'prewitt_h',
+'prewitt_v',
+'rank_order',
+'roberts',
+'roberts_neg_diag',
+'roberts_pos_diag',
+'sato',
+'scharr',
+'scharr_h',
+'scharr_v',
+'sobel',
+'sobel_h',
+'sobel_v',
+'threshold_isodata',
+'threshold_li',
+'threshold_local',
+'threshold_mean',
+'threshold_minimum',
+'threshold_multiotsu',
+'threshold_niblack',
+'threshold_otsu',
+'threshold_sauvola',
+'threshold_triangle',
+'threshold_yen',
+'try_all_threshold',
+'unsharp_mask',
+'wiener',
+'window'
+]
+
+from ._fft_based import butterworth
+from ._gabor import gabor, gabor_kernel
+from ._gaussian import difference_of_gaussians, gaussian
+from ._median import median
+from ._rank_order import rank_order
+from ._sparse import correlate_sparse
+from ._unsharp_mask import unsharp_mask
+from ._window import window
+from .edges import (farid, farid_h, farid_v, laplace, prewitt, prewitt_h,
+                    prewitt_v, roberts, roberts_neg_diag, roberts_pos_diag,
+                    scharr, scharr_h, scharr_v, sobel, sobel_h, sobel_v)
+from .lpi_filter import LPIFilter2D, filter_forward, filter_inverse, wiener
+from .ridges import frangi, hessian, meijering, sato
+from .thresholding import (apply_hysteresis_threshold, threshold_isodata,
+                           threshold_li, threshold_local, threshold_mean,
+                           threshold_minimum, threshold_multiotsu,
+                           threshold_niblack, threshold_otsu, threshold_sauvola,
+                           threshold_triangle, threshold_yen, try_all_threshold)
diff --git a/python/cucim/src/cucim/skimage/filters/_fft_based.py b/python/cucim/src/cucim/skimage/filters/_fft_based.py
index 5e1b76ba0..9a68ea7b1 100644
--- a/python/cucim/src/cucim/skimage/filters/_fft_based.py
+++ b/python/cucim/src/cucim/skimage/filters/_fft_based.py
@@ -159,8 +159,8 @@ def butterworth(
     elif npad > 0:
         center_slice = tuple(slice(npad, s + npad) for s in image.shape)
         image = pad(image, npad, mode='edge')
-    fft_shape = (image.shape if channel_axis is None
-                 else np.delete(image.shape, channel_axis))
+    fft_shape = tuple(image.shape if channel_axis is None
+                      else np.delete(image.shape, channel_axis))
     is_real = cp.isrealobj(image)
     float_dtype = _supported_float_type(image.dtype, allow_complex=True)
     image = image.astype(float_dtype, copy=False)
diff --git a/python/cucim/src/cucim/skimage/filters/lpi_filter.py b/python/cucim/src/cucim/skimage/filters/lpi_filter.py
index d6cd59c0b..e443a0723 100644
--- a/python/cucim/src/cucim/skimage/filters/lpi_filter.py
+++ b/python/cucim/src/cucim/skimage/filters/lpi_filter.py
@@ -7,7 +7,7 @@
 import numpy as np
 from cupyx.scipy import fft
 
-from .._shared.utils import _supported_float_type, check_nD, deprecated
+from .._shared.utils import _supported_float_type, check_nD, deprecate_func
 
 eps = np.finfo(float).eps
 
@@ -19,7 +19,7 @@ def _min_limit(x, val=eps):
 
 def _center(x, oshape):
     """Return an array of shape ``oshape`` from the center of array ``x``."""
-    start = (np.array(x.shape) - np.array(oshape)) // 2 + 1
+    start = (np.array(x.shape) - np.array(oshape)) // 2
     out = x[tuple(slice(s, s + n) for s, n in zip(start, oshape))]
     return out
 
@@ -65,11 +65,10 @@ def __init__(self, impulse_response, **filter_params):
 
         Examples
         --------
-        Gaussian filter: Use a 1-D gaussian in each direction without
-        normalization coefficients.
+        Gaussian filter without normalization of coefficients:
 
         >>> def filt_func(r, c, sigma = 1):
-        ...     return cp.exp(-cp.hypot(r, c)/sigma)
+        ...     return cp.exp(-(r**2 + c**2)/(2 * sigma**2))
         >>> filter = LPIFilter2D(filt_func)
 
         """
@@ -83,14 +82,22 @@ def __init__(self, impulse_response, **filter_params):
     def _prepare(self, data):
         """Calculate filter and data FFT in preparation for filtering."""
         dshape = np.array(data.shape)
-        dshape += dshape % 2 == 0  # all filter dimensions must be uneven
-        oshape = np.array(data.shape) * 2 - 1
+        # all filter dimensions must be uneven
+        even_offset = tuple(int(s % 2 == 0) for s in data.shape)
+        dshape = tuple(
+            s + offset for s, offset in zip(data.shape, even_offset)
+        )
+
+        oshape = tuple(s * 2 - 1 for s in data.shape)
 
         float_dtype = _supported_float_type(data.dtype)
         data = data.astype(float_dtype, copy=False)
 
         if self._cache is None or np.any(self._cache.shape != oshape):
-            coords = cp.mgrid[[slice(0, float(n)) for n in dshape]]
+            coords = cp.mgrid[
+                [slice(0 + offset, float(n + offset))
+                 for (n, offset) in zip(dshape, even_offset)]
+            ]
             # this steps over two sets of coordinates,
             # not over the coordinates individually
             for k, coord in enumerate(coords):
@@ -127,7 +134,7 @@ def __call__(self, data):
         return out
 
 
-def filter_forward(data, impulse_response=None, filter_params={},
+def filter_forward(data, impulse_response=None, filter_params=None,
                    predefined_filter=None):
     """Apply the given filter to data.
 
@@ -137,7 +144,7 @@ def filter_forward(data, impulse_response=None, filter_params={},
         Input data.
     impulse_response : callable `f(r, c, **filter_params)`
         Impulse response of the filter.  See LPIFilter2D.__init__.
-    filter_params : dict
+    filter_params : dict, optional
         Additional keyword parameters to the impulse_response function.
 
     Other Parameters
@@ -148,32 +155,35 @@ def filter_forward(data, impulse_response=None, filter_params={},
 
     Examples
     --------
+    Gaussian filter without normalization:
 
-    Gaussian filter:
-
-    >>> import cupy as cp
-    >>> def filt_func(r, c):
-    ...     return cp.exp(-cp.hypot(r, c)/1)
+    >>> def filt_func(r, c, sigma=1):
+    ...     return cp.exp(-(r**2 + c**2)/(2 * sigma**2))
     >>>
     >>> from skimage import data
     >>> filtered = filter_forward(cp.array(data.coins()), filt_func)
 
     """
+    if filter_params is None:
+        filter_params = {}
     check_nD(data, 2, 'data')
     if predefined_filter is None:
         predefined_filter = LPIFilter2D(impulse_response, **filter_params)
     return predefined_filter(data)
 
 
-@deprecated(alt_func='cucim.skimage.filters.lpi_filter.filter_inverse',
-            removed_version='2023.06.00')
-def inverse(data, impulse_response=None, filter_params={}, max_gain=2,
+@deprecate_func(
+    hint="use cucim.skimage.filters.lpi_filter.filter_inverse instead",
+    deprecated_version="",
+    removed_version="2023.12.00",
+)
+def inverse(data, impulse_response=None, filter_params=None, max_gain=2,
             predefined_filter=None):
     return filter_inverse(data, impulse_response, filter_params,
                           max_gain, predefined_filter)
 
 
-def filter_inverse(data, impulse_response=None, filter_params={}, max_gain=2,
+def filter_inverse(data, impulse_response=None, filter_params=None, max_gain=2,
                    predefined_filter=None):
     """Apply the filter in reverse to the given data.
 
@@ -182,21 +192,24 @@ def filter_inverse(data, impulse_response=None, filter_params={}, max_gain=2,
     data : (M, N) ndarray
         Input data.
     impulse_response : callable `f(r, c, **filter_params)`
-        Impulse response of the filter.  See LPIFilter2D.__init__.
-    filter_params : dict
+        Impulse response of the filter.  See :class:`~.LPIFilter2D`. This is a
+        required argument unless a `predifined_filter` is provided.
+    filter_params : dict, optional
         Additional keyword parameters to the impulse_response function.
-    max_gain : float
+    max_gain : float, optional
         Limit the filter gain.  Often, the filter contains zeros, which would
         cause the inverse filter to have infinite gain.  High gain causes
         amplification of artefacts, so a conservative limit is recommended.
 
     Other Parameters
     ----------------
-    predefined_filter : LPIFilter2D
+    predefined_filter : LPIFilter2D, optional
         If you need to apply the same filter multiple times over different
         images, construct the LPIFilter2D and specify it here.
 
     """
+    if filter_params is None:
+        filter_params = {}
     check_nD(data, 2, 'data')
     if predefined_filter is None:
         filt = LPIFilter2D(impulse_response, **filter_params)
@@ -213,7 +226,7 @@ def filter_inverse(data, impulse_response=None, filter_params={}, max_gain=2,
     return _center(cp.abs(fft.ifftshift(fft.ifftn(G * F))), data.shape)
 
 
-def wiener(data, impulse_response=None, filter_params={}, K=0.25,
+def wiener(data, impulse_response=None, filter_params=None, K=0.25,
            predefined_filter=None):
     """Minimum Mean Square Error (Wiener) inverse filter.
 
@@ -226,7 +239,7 @@ def wiener(data, impulse_response=None, filter_params={}, K=0.25,
         image.
     impulse_response : callable `f(r, c, **filter_params)`
         Impulse response of the filter.  See LPIFilter2D.__init__.
-    filter_params : dict
+    filter_params : dict, optional
         Additional keyword parameters to the impulse_response function.
 
     Other Parameters
@@ -236,6 +249,8 @@ def wiener(data, impulse_response=None, filter_params={}, K=0.25,
         images, construct the LPIFilter2D and specify it here.
 
     """
+    if filter_params is None:
+        filter_params = {}
     check_nD(data, 2, 'data')
 
     if not isinstance(K, float):
diff --git a/python/cucim/src/cucim/skimage/filters/ridges.py b/python/cucim/src/cucim/skimage/filters/ridges.py
index 5d861b8f7..3615278bc 100644
--- a/python/cucim/src/cucim/skimage/filters/ridges.py
+++ b/python/cucim/src/cucim/skimage/filters/ridges.py
@@ -14,13 +14,13 @@ class of ridge filters relies on the eigenvalues of the Hessian matrix of
 import cupy as cp
 import numpy as np
 
-from .._shared.utils import _supported_float_type, check_nD, deprecated
+from .._shared.utils import _supported_float_type, check_nD, deprecate_func
 from ..feature.corner import (_symmetric_compute_eigenvalues, hessian_matrix,
                               hessian_matrix_eigvals)
 from ..util import img_as_float
 
 
-@deprecated(removed_version="2023.06.01")
+@deprecate_func(deprecated_version="", removed_version="2023.06.01")
 def compute_hessian_eigenvalues(image, sigma, sorting='none',
                                 mode='constant', cval=0,
                                 use_gaussian_derivatives=False):
diff --git a/python/cucim/src/cucim/skimage/filters/thresholding.py b/python/cucim/src/cucim/skimage/filters/thresholding.py
index 367eb35d8..c2dedb2fa 100644
--- a/python/cucim/src/cucim/skimage/filters/thresholding.py
+++ b/python/cucim/src/cucim/skimage/filters/thresholding.py
@@ -704,6 +704,10 @@ def threshold_li(image, *, tolerance=None, initial_guess=None,
                                    weights=hist[foreground])
             mean_back = np.average(bin_centers[background],
                                    weights=hist[background])
+
+            if mean_back == 0:
+                break
+
             eps = 100 * np.finfo(float).eps
             mean_back = float(mean_back)
             mean_fore = float(mean_fore)
@@ -726,6 +730,9 @@ def threshold_li(image, *, tolerance=None, initial_guess=None,
             mean_fore = float(cp.mean(image[foreground]))
             mean_back = float(cp.mean(image[~foreground]))
 
+            if mean_back == 0:
+                break
+
             t_next = ((mean_back - mean_fore) /
                       (math.log(mean_back + eps) - math.log(mean_fore + eps)))
 
diff --git a/python/cucim/src/cucim/skimage/measure/_colocalization.py b/python/cucim/src/cucim/skimage/measure/_colocalization.py
index f79c43a44..3587a331b 100644
--- a/python/cucim/src/cucim/skimage/measure/_colocalization.py
+++ b/python/cucim/src/cucim/skimage/measure/_colocalization.py
@@ -176,12 +176,6 @@ def manders_coloc_coeff(image0, image1_mask, mask=None):
     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
 
@@ -261,7 +255,8 @@ def manders_overlap_coeff(image0, image1, mask=None):
     if image1.min() < 0:
         raise ValueError("image1 contains negative values")
 
-    return _get_manders_overlap_coeff(image0, image1)
+    denom = cp.linalg.norm(image0) * cp.linalg.norm(image1)
+    return cp.vdot(image0, image1) / denom
 
 
 def intersection_coeff(image0_mask, image1_mask, mask=None):
diff --git a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
index d4ca12fde..a58ba8681 100644
--- a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
+++ b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
@@ -153,10 +153,10 @@ def thin(image, max_num_iter=None):
 # --------- Skeletonization by medial axis transform --------
 
 
-def _get_tiebreaker(n, random_seed):
+def _get_tiebreaker(n, seed):
     # CuPy generator doesn't currently have the permutation method, so
     # fall back to cp.random.permutation instead.
-    cp.random.seed(random_seed)
+    cp.random.seed(seed)
     if n < 2 << 31:
         dtype = np.int32
     else:
@@ -165,7 +165,12 @@ def _get_tiebreaker(n, random_seed):
     return tiebreaker
 
 
-def medial_axis(image, mask=None, return_distance=False, *, random_state=None):
+@deprecate_kwarg(
+    {'random_state': 'seed'},
+    deprecated_version='23.08',
+    removed_version='24.06'
+)
+def medial_axis(image, mask=None, return_distance=False, *, seed=None):
     """Compute the medial axis transform of a binary image.
 
     Parameters
@@ -177,13 +182,12 @@ def medial_axis(image, mask=None, return_distance=False, *, random_state=None):
         value in `mask` are used for computing the medial axis.
     return_distance : bool, optional
         If true, the distance transform is returned as well as the skeleton.
-    random_state : {None, int, `numpy.random.Generator`}, optional
-        If `random_state` is None the `numpy.random.Generator` singleton is
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None, the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used, seeded with
+        `seed`.
+        If `seed` is already a ``Generator`` instance, then that instance is
         used.
-        If `random_state` is an int, a new ``Generator`` instance is used,
-        seeded with `random_state`.
-        If `random_state` is already a ``Generator`` instance then that
-        instance is used.
 
         .. versionadded:: 0.19
 
@@ -295,7 +299,7 @@ def medial_axis(image, mask=None, return_distance=False, *, random_state=None):
     # We use a random # for tiebreaking. Assign each pixel in the image a
     # predictable, random # so that masking doesn't affect arbitrary choices
     # of skeletons
-    tiebreaker = _get_tiebreaker(n=distance.size, random_seed=random_state)
+    tiebreaker = _get_tiebreaker(n=distance.size, seed=seed)
     order = cp.lexsort(
         cp.stack(
             (tiebreaker, corner_score[masked_image], distance),
diff --git a/python/cucim/src/cucim/skimage/restoration/__init__.py b/python/cucim/src/cucim/skimage/restoration/__init__.py
index aa0b63ecb..07eaeb6b1 100644
--- a/python/cucim/src/cucim/skimage/restoration/__init__.py
+++ b/python/cucim/src/cucim/skimage/restoration/__init__.py
@@ -1,6 +1,6 @@
 from ._denoise import denoise_tv_chambolle
 from .deconvolution import richardson_lucy, unsupervised_wiener, wiener
-from .j_invariant import calibrate_denoiser
+from .j_invariant import calibrate_denoiser, denoise_invariant
 
 __all__ = [
     "wiener",
@@ -8,4 +8,5 @@
     "richardson_lucy",
     "denoise_tv_chambolle",
     "calibrate_denoiser",
+    "denoise_invariant",
 ]
diff --git a/python/cucim/src/cucim/skimage/restoration/deconvolution.py b/python/cucim/src/cucim/skimage/restoration/deconvolution.py
index 923505a4a..f9dda2e30 100644
--- a/python/cucim/src/cucim/skimage/restoration/deconvolution.py
+++ b/python/cucim/src/cucim/skimage/restoration/deconvolution.py
@@ -146,8 +146,10 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
     return deconv
 
 
+@deprecate_kwarg({'random_state': 'seed'}, removed_version="23.08.00",
+                 deprecated_version="24.06.00")
 def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
-                        clip=True, *, random_state=None):
+                        clip=True, *, seed=None):
     """Unsupervised Wiener-Hunt deconvolution.
 
     Return the deconvolution with a Wiener-Hunt approach, where the
@@ -172,13 +174,12 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
     clip : boolean, optional
        True by default. If true, pixel values of the result above 1 or
        under -1 are thresholded for skimage pipeline compatibility.
-    random_state : {None, int, `cupy.random.Generator`}, optional
-        If `random_state` is None the `cupy.random.Generator` singleton is
+    seed : {None, int, `numpy.random.Generator`}, optional
+        If `seed` is None, the `numpy.random.Generator` singleton is used.
+        If `seed` is an int, a new ``Generator`` instance is used, seeded with
+        `seed`.
+        If `seed` is already a ``Generator`` instance, then that instance is
         used.
-        If `random_state` is an int, a new ``Generator`` instance is used,
-        seeded with `random_state`.
-        If `random_state` is already a ``Generator`` instance then that
-        instance is used.
 
     Returns
     -------
@@ -304,10 +305,10 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         data_spectrum = uft.ufft2(image)
 
     try:
-        rng = cp.random.default_rng(random_state)
+        rng = cp.random.default_rng(seed)
     except AttributeError:
         # older CuPy without default_rng
-        rng = cp.random.RandomState(random_state)
+        rng = cp.random.RandomState(seed)
 
     # Gibbs sampling
     for iteration in range(params["max_num_iter"]):
@@ -383,8 +384,6 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
     return (x_postmean, {'noise': gn_chain, 'prior': gx_chain})
 
 
-@deprecate_kwarg({'iterations': 'num_iter'}, removed_version="23.02.00",
-                 deprecated_version="22.02.00")
 def richardson_lucy(image, psf, num_iter=50, clip=True, filter_epsilon=None):
     """Richardson-Lucy deconvolution.
 
diff --git a/python/cucim/src/cucim/skimage/restoration/j_invariant.py b/python/cucim/src/cucim/skimage/restoration/j_invariant.py
index e8a58dab2..1e19a92ca 100644
--- a/python/cucim/src/cucim/skimage/restoration/j_invariant.py
+++ b/python/cucim/src/cucim/skimage/restoration/j_invariant.py
@@ -89,14 +89,16 @@ def _generate_grid_slice(shape, *, offset, stride=3):
     return mask
 
 
-def _invariant_denoise(image, denoise_function, *, stride=4,
-                       masks=None, denoiser_kwargs=None):
+def denoise_invariant(image, denoise_function, *, stride=4, masks=None,
+                      denoiser_kwargs=None):
     """Apply a J-invariant version of `denoise_function`.
 
     Parameters
     ----------
-    image : ndarray
-        Input data to be denoised (converted using `img_as_float`).
+    image : ndarray ([M[, N[, ...P]][, C]) of ints, uints or floats
+        Input data to be denoised. `image` can be of any numeric type,
+        but it is cast into a ndarray of floats (using `img_as_float`) for the
+        computation of the denoised image.
     denoise_function : function
         Original denoising function.
     stride : int, optional
@@ -112,6 +114,32 @@ def _invariant_denoise(image, denoise_function, *, stride=4,
     -------
     output : ndarray
         Denoised image, of same shape as `image`.
+
+    Notes
+    -----
+    A denoising function is J-invariant if the prediction it makes for each
+    pixel does not depend on the value of that pixel in the original image.
+    The prediction for each pixel may instead use all the relevant information
+    contained in the rest of the image, which is typically quite significant.
+    Any function can be converted into a J-invariant one using a simple masking
+    procedure, as described in [1].
+
+    The pixel-wise error of a J-invariant denoiser is uncorrelated to the noise,
+    so long as the noise in each pixel is independent. Consequently, the average
+    difference between the denoised image and the oisy image, the
+    *self-supervised loss*, is the same as the difference between the denoised
+    image and the original clean image, the *ground-truth loss* (up to a
+    constant).
+
+    This means that the best J-invariant denoiser for a given image can be found
+    using the noisy data alone, by selecting the denoiser minimizing the self-
+    supervised loss.
+
+    References
+    ----------
+    .. [1] J. Batson & L. Royer. Noise2Self: Blind Denoising by
+       Self-Supervision, International Conference on Machine Learning,
+       p. 524-533 (2019).
     """
     image = img_as_float(image)
 
@@ -252,7 +280,7 @@ def calibrate_denoiser(image, denoise_function, denoise_parameters, *,
     best_parameters = parameters_tested[idx]
 
     best_denoise_function = functools.partial(
-        _invariant_denoise,
+        denoise_invariant,
         denoise_function=denoise_function,
         stride=stride,
         denoiser_kwargs=best_parameters,
@@ -303,7 +331,7 @@ def _calibrate_denoiser_search(image, denoise_function, denoise_parameters, *,
             multichannel = \
                 denoiser_kwargs.get('channel_axis', None) is not None
         if not approximate_loss:
-            denoised = _invariant_denoise(
+            denoised = denoise_invariant(
                 image, denoise_function,
                 stride=stride,
                 denoiser_kwargs=denoiser_kwargs
@@ -315,7 +343,7 @@ def _calibrate_denoiser_search(image, denoise_function, denoise_parameters, *,
             mask = _generate_grid_slice(image.shape[:spatialdims],
                                         offset=n_masks // 2, stride=stride)
 
-            masked_denoised = _invariant_denoise(
+            masked_denoised = denoise_invariant(
                 image, denoise_function,
                 masks=[mask],
                 denoiser_kwargs=denoiser_kwargs
diff --git a/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
index bab2f1887..84a7feb2b 100644
--- a/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
+++ b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
@@ -8,37 +8,6 @@
 from .._vendored import pad
 
 
-@cp.fuse()
-def _fused_curvature(phi, x_start, x_end, y_start, y_end, ul, ur, ll, lr):
-    fy = (y_end - y_start) / 2.0
-    fx = (x_end - x_start) / 2.0
-    fyy = y_end + y_start - 2 * phi
-    fxx = x_end + x_start - 2 * phi
-    fxy = .25 * (lr + ul - ur - ll)
-    grad2 = fx**2 + fy**2
-    K = (fxx * fy**2 - 2 * fxy * fx * fy + fyy * fx**2)
-    K /= (grad2 * cp.sqrt(grad2) + 1e-8)
-    return K
-
-
-def _cv_curvature(phi):
-    """Returns the 'curvature' of a level set 'phi'.
-    """
-    P = pad(phi, 1, mode='edge')
-    y_start = P[:-2, 1:-1]
-    y_end = P[2:, 1:-1]
-    x_start = P[1:-1, :-2]
-    x_end = P[1:-1, 2:]
-
-    lower_right = P[2:, 2:]
-    lower_left = P[2:, :-2]
-    upper_right = P[:-2, 2:]
-    upper_left = P[:-2, :-2]
-    K = _fused_curvature(phi, x_start, x_end, y_start, y_end, upper_left,
-                         upper_right, lower_left, lower_right)
-    return K
-
-
 @cp.fuse()
 def _fused_variance_kernel1(eta, x_start, x_mid, x_end, y_start, y_mid, y_end):
     phixp = x_end - x_mid
@@ -102,6 +71,21 @@ def _fused_variance_kernel2(
 
 def _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt):
     """Returns the variation of level set 'phi' based on algorithm parameters.
+
+    This corresponds to equation (22) of the paper by Pascal Getreuer,
+    which computes the next iteration of the level set based on a current
+    level set.
+
+    A full explanation regarding all the terms is beyond the scope of the
+    present description, but there is one difference of particular import.
+    In the original algorithm, convergence is accelerated, and required
+    memory is reduced, by using a single array. This array, therefore, is a
+    combination of non-updated and updated values. If this were to be
+    implemented in python, this would require a double loop, where the
+    benefits of having fewer iterations would be outweided by massively
+    increasing the time required to perform each individual iteration. A
+    similar approach is used by Rami Cohen, and it is from there that the
+    C1-4 notation is taken.
     """
     eta = 1e-16
     P = pad(phi, 1, mode='edge')
@@ -183,17 +167,44 @@ def _cv_difference_from_average_term(image, Hphi, lambda_pos, lambda_neg):
     return out
 
 
+@cp.fuse()
+def _fused_edge_length(mu, phi, x_start, x_end, y_start, y_end):
+    fy = (y_end - y_start) / 2.0
+    fx = (x_end - x_start) / 2.0
+    grad2 = fx**2 + fy**2
+
+    del_phi = 1.0 / (1.0 + phi * phi)  # _cv_delta(phi)
+    K = mu * del_phi * cp.sqrt(grad2)
+    return K
+
+
 def _cv_edge_length_term(phi, mu):
-    """Returns the 'energy' contribution due to the length of the
-    edge between regions at each point, multiplied by a factor 'mu'.
+    """Returns the 'curvature' of a level set 'phi'.
     """
-    e = _cv_curvature(phi)
-    e *= mu
-    return e
+    P = pad(phi, 1, mode='edge')
+    y_start = P[:-2, 1:-1]
+    y_end = P[2:, 1:-1]
+    x_start = P[1:-1, :-2]
+    x_end = P[1:-1, 2:]
+
+    K = _fused_edge_length(mu, phi, x_start, x_end, y_start, y_end)
+    return K
 
 
 def _cv_energy(image, phi, mu, lambda1, lambda2):
     """Returns the total 'energy' of the current level set function.
+
+    This corresponds to equation (7) of the paper by Pascal Getreuer,
+    which is the weighted sum of the following:
+    (A) the length of the contour produced by the zero values of the
+    level set,
+    (B) the area of the "foreground" (area of the image where the
+    level set is positive),
+    (C) the variance of the image inside the foreground,
+    (D) the variance of the image outside of the foreground
+
+    Each value is computed for each pixel, and then summed. The weight
+    of (B) is set to 0 in this implementation.
     """
     H = _cv_heavyside(phi)
     avgenergy = _cv_difference_from_average_term(image, H, lambda1, lambda2)
diff --git a/python/cucim/src/cucim/skimage/segmentation/_join.py b/python/cucim/src/cucim/skimage/segmentation/_join.py
index 3744de969..6227a2fea 100644
--- a/python/cucim/src/cucim/skimage/segmentation/_join.py
+++ b/python/cucim/src/cucim/skimage/segmentation/_join.py
@@ -3,7 +3,7 @@
 from ..util._map_array import ArrayMap, map_array
 
 
-def join_segmentations(s1, s2):
+def join_segmentations(s1, s2, return_mapping: bool = False):
     """Return the join of the two input segmentations.
 
     The join J of S1 and S2 is defined as the segmentation in which two
@@ -14,11 +14,18 @@ def join_segmentations(s1, s2):
     ----------
     s1, s2 : numpy arrays
         s1 and s2 are label fields of the same shape.
+    return_mapping : bool, optional
+        If true, return mappings for joined segmentation labels to the original
+        labels.
 
     Returns
     -------
     j : numpy array
         The join segmentation of s1 and s2.
+    map_j_to_s1 : ArrayMap, optional
+        Mapping from labels of the joined segmentation j to labels of s1.
+    map_j_to_s2 : ArrayMap, optional
+        Mapping from labels of the joined segmentation j to labels of s2.
 
     Examples
     --------
@@ -38,11 +45,21 @@ def join_segmentations(s1, s2):
     if s1.shape != s2.shape:
         raise ValueError("Cannot join segmentations of different shape. "
                          f"s1.shape: {s1.shape}, s2.shape: {s2.shape}")
-    s1 = relabel_sequential(s1)[0]
-    s2 = relabel_sequential(s2)[0]
-    j = (s2.max() + 1) * s1 + s2
-    j = relabel_sequential(j)[0]
-    return j
+    s1_relabeled, _, backward_map1 = relabel_sequential(s1)
+    s2_relabeled, _, backward_map2 = relabel_sequential(s2)
+    factor = s2.max() + 1
+    j_initial = factor * s1_relabeled + s2_relabeled
+    j, _, map_j_to_j_initial = relabel_sequential(j_initial)
+    if not return_mapping:
+        return j
+    # Determine label mapping
+    labels_j = cp.unique(j_initial)
+    labels_s1_relabeled, labels_s2_relabeled = cp.divmod(labels_j, factor)
+    map_j_to_s1 = ArrayMap(map_j_to_j_initial.in_values,
+                           backward_map1[labels_s1_relabeled])
+    map_j_to_s2 = ArrayMap(map_j_to_j_initial.in_values,
+                           backward_map2[labels_s2_relabeled])
+    return j, map_j_to_s1, map_j_to_s2
 
 
 def relabel_sequential(label_field, offset=1):
diff --git a/python/cucim/src/cucim/skimage/transform/_geometric.py b/python/cucim/src/cucim/skimage/transform/_geometric.py
index bb6bde4bf..4da5b11f4 100644
--- a/python/cucim/src/cucim/skimage/transform/_geometric.py
+++ b/python/cucim/src/cucim/skimage/transform/_geometric.py
@@ -487,6 +487,39 @@ class EssentialMatrixTransform(FundamentalMatrixTransform):
     params : (3, 3) ndarray
         Essential matrix.
 
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from cucim.skimage import transform
+    >>>
+    >>> tform_matrix = transform.EssentialMatrixTransform(
+    ...     rotation=cp.eye(3), translation=cp.array([0, 0, 1])
+    ... )
+    >>> tform_matrix.params
+    array([[ 0., -1.,  0.],
+           [ 1.,  0.,  0.],
+           [ 0.,  0.,  0.]])
+    >>> src = cp.array([[ 1.839035, 1.924743],
+    ...                 [ 0.543582, 0.375221],
+    ...                 [ 0.47324 , 0.142522],
+    ...                 [ 0.96491 , 0.598376],
+    ...                 [ 0.102388, 0.140092],
+    ...                 [15.994343, 9.622164],
+    ...                 [ 0.285901, 0.430055],
+    ...                 [ 0.09115 , 0.254594]])
+    >>> dst = cp.array([[1.002114, 1.129644],
+    ...                 [1.521742, 1.846002],
+    ...                 [1.084332, 0.275134],
+    ...                 [0.293328, 0.588992],
+    ...                 [0.839509, 0.08729 ],
+    ...                 [1.779735, 1.116857],
+    ...                 [0.878616, 0.602447],
+    ...                 [0.642616, 1.028681]])
+    >>> tform_matrix.estimate(src, dst)
+    True
+    >>> tform_matrix.residuals(src, dst)
+    array([0.42455187, 0.01460448, 0.13847034, 0.12140951, 0.27759346,
+           0.32453118, 0.00210776, 0.26512283])
     """
 
     # CuPy Backend: if matrix is None cannot infer array module from it
@@ -874,14 +907,26 @@ class AffineTransform(ProjectiveTransform):
 
     Has the following form::
 
-        X = a0*x + a1*y + a2 =
-          = sx*x*cos(rotation) - sy*y*sin(rotation + shear) + a2
+        X = a0 * x + a1 * y + a2
+          =   sx * x * [cos(rotation) + tan(shear_y) * sin(rotation)]
+            - sy * y * [tan(shear_x) * cos(rotation) + sin(rotation)]
+            + translation_x
+
+        Y = b0 * x + b1 * y + b2
+          =   sx * x * [sin(rotation) - tan(shear_y) * cos(rotation)]
+            - sy * y * [tan(shear_x) * sin(rotation) - cos(rotation)]
+            + translation_y
+
+    where ``sx`` and ``sy`` are scale factors in the x and y directions.
 
-        Y = b0*x + b1*y + b2 =
-          = sx*x*sin(rotation) + sy*y*cos(rotation + shear) + b2
+    This is equivalent to applying the operations in the following order:
 
-    where ``sx`` and ``sy`` are scale factors in the x and y directions,
-    and the homogeneous transformation matrix is::
+    1. Scale
+    2. Shear
+    3. Rotate
+    4. Translate
+
+    The homogeneous transformation matrix is::
 
         [[a0  a1  a2]
          [b0  b1  b2]
@@ -909,11 +954,12 @@ class AffineTransform(ProjectiveTransform):
         .. versionadded:: 0.17
            Added support for supplying a single scalar value.
     rotation : float, optional
-        Rotation angle in counter-clockwise direction as radians. Only
-        available for 2D.
-    shear : float, optional
-        Shear angle in counter-clockwise direction as radians. Only available
-        for 2D.
+        Rotation angle, clockwise, as radians. Only available for 2D.
+    shear : float or 2-tuple of float, optional
+        The x and y shear angles, clockwise, by which these axes are
+        rotated around the origin [2].
+        If a single value is given, take that to be the x shear angle, with
+        the y angle remaining 0. Only available in 2D.
     translation : (tx, ty) as ndarray, list or tuple, optional
         Translation parameters. Only available for 2D.
     dimensionality : int, optional
@@ -929,6 +975,39 @@ class AffineTransform(ProjectiveTransform):
     ------
     ValueError
         If both ``matrix`` and any of the other parameters are provided.
+
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from cucim.skimage import transform
+    >>> from skimage import data
+    >>> img = cp.array(data.astronaut())
+
+    Define source and destination points:
+
+    >>> src = cp.array([[150, 150],
+    ...                 [250, 100],
+    ...                 [150, 200]])
+    >>> dst = cp.array([[200, 200],
+    ...                 [300, 150],
+    ...                 [150, 400]])
+
+    Estimate the transformation matrix:
+
+    >>> tform = transform.AffineTransform()
+    >>> tform.estimate(src, dst)
+    True
+
+    Apply the transformation:
+
+    >>> warped = transform.warp(img, inverse_map=tform.inverse)
+
+    References
+    ----------
+    .. [1] Wikipedia, "Affine transformation",
+           https://en.wikipedia.org/wiki/Affine_transformation#Image_transformation
+    .. [2] Wikipedia, "Shear mapping",
+           https://en.wikipedia.org/wiki/Shear_mapping
     """
 
     def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
@@ -968,17 +1047,27 @@ def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
             else:
                 sx, sy = scale
 
-            # fmt: off
-            self.params = np.array(
-                [
-                    [sx * _cos(rotation), -sy * _sin(rotation + shear), 0],  # NOQA
-                    [sx * _sin(rotation),  sy * _cos(rotation + shear), 0],  # NOQA
-                    [                  0,                            0, 1],  # NOQA
-                ]
+            if np.isscalar(shear):
+                shear_x, shear_y = (shear, 0)
+            else:
+                shear_x, shear_y = shear
+
+            a0 = sx * (
+                math.cos(rotation) + math.tan(shear_y) * math.sin(rotation)
             )
-            # fmt: on
-            self.params[0:2, 2] = translation
-            self.params = xp.asarray(self.params)
+            a1 = -sy * (
+                math.tan(shear_x) * math.cos(rotation) + math.sin(rotation)
+            )
+            a2 = translation[0]
+
+            b0 = sx * (
+                math.sin(rotation) - math.tan(shear_y) * math.cos(rotation)
+            )
+            b1 = -sy * (
+                math.tan(shear_x) * math.sin(rotation) - math.cos(rotation)
+            )
+            b2 = translation[1]
+            self.params = xp.array([[a0, a1, a2], [b0, b1, b2], [0, 0, 1]])
         else:
             # default to an identity transform
             self.params = xp.eye(dimensionality + 1)
@@ -986,9 +1075,14 @@ def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
     @property
     def scale(self):
         xp = cp.get_array_module(self.params)
-        return xp.sqrt(xp.sum(self.params * self.params, axis=0))[
-            : self.dimensionality
-        ]
+        if self.dimensionality != 2:
+            return xp.sqrt(xp.sum(self.params * self.params, axis=0))[
+                : self.dimensionality
+            ]
+        else:
+            ss = xp.sum(self.params * self.params, axis=0)
+            ss[1] = ss[1] / (math.tan(self.shear) ** 2 + 1)
+            return xp.sqrt(ss)[:self.dimensionality]
 
     @property
     def rotation(self):
@@ -1066,7 +1160,13 @@ def estimate(self, src, dst):
 
         # find affine mapping from source positions to destination
         self.affines = []
-        for tri in xp.asarray(self._tesselation.vertices):
+
+        try:
+            tesselation_simplices = self._tesselation.simplices
+        except AttributeError:
+            # vertices is deprecated and scheduled for removal in SciPy 1.11
+            tesselation_simplices = self._tesselation.vertices
+        for tri in xp.asarray(tesselation_simplices):
             affine = AffineTransform(dimensionality=ndim)
             ok &= affine.estimate(src[tri, :], dst[tri, :])
             self.affines.append(affine)
@@ -1077,7 +1177,12 @@ def estimate(self, src, dst):
         self._inverse_tesselation = spatial.Delaunay(cp.asnumpy(dst))
         # find affine mapping from source positions to destination
         self.inverse_affines = []
-        for tri in xp.asarray(self._inverse_tesselation.vertices):
+        try:
+            inv_tesselation_simplices = self._inverse_tesselation.simplices
+        except AttributeError:
+            # vertices is deprecated and scheduled for removal in SciPy 1.11
+            inv_tesselation_simplices = self._inverse_tesselation.vertices
+        for tri in xp.asarray(inv_tesselation_simplices):
             affine = AffineTransform(dimensionality=ndim)
             ok &= affine.estimate(dst[tri, :], src[tri, :])
             self.inverse_affines.append(affine)
@@ -1112,7 +1217,12 @@ def __call__(self, coords):
         # coordinates outside of mesh
         out[simplex == -1, :] = -1
 
-        for index in range(len(self._tesselation.vertices)):
+        try:
+            tesselation_simplices = self._tesselation.simplices
+        except AttributeError:
+            # vertices is deprecated and scheduled for removal in SciPy 1.11
+            tesselation_simplices = self._tesselation.vertices
+        for index in range(len(tesselation_simplices)):
             # affine transform for triangle
             affine = self.affines[index]
             # all coordinates within triangle
@@ -1150,7 +1260,12 @@ def inverse(self, coords):
         # coordinates outside of mesh
         out[simplex == -1, :] = -1
 
-        for index in range(len(self._inverse_tesselation.vertices)):
+        try:
+            inv_tesselation_simplices = self._inverse_tesselation.simplices
+        except AttributeError:
+            # vertices is deprecated and scheduled for removal in SciPy 1.11
+            inv_tesselation_simplices = self._inverse_tesselation.vertices
+        for index in range(len(inv_tesselation_simplices)):
             # affine transform for triangle
             affine = self.inverse_affines[index]
             # all coordinates within triangle
@@ -1244,7 +1359,7 @@ class EuclideanTransform(ProjectiveTransform):
     matrix : (D+1, D+1) ndarray, optional
         Homogeneous transformation matrix.
     rotation : float or sequence of float, optional
-        Rotation angle in counter-clockwise direction as radians. If given as
+        Rotation angle, clockwise, as radians. If given as
         a vector, it is interpreted as Euler rotation angles [1]_. Only 2D
         (single rotation) and 3D (Euler rotations) values are supported. For
         higher dimensions, you must provide or estimate the transformation
@@ -1390,7 +1505,7 @@ class SimilarityTransform(EuclideanTransform):
     scale : float, optional
         Scale factor. Implemented only for 2D and 3D.
     rotation : float, optional
-        Rotation angle in counter-clockwise direction as radians.
+        Rotation angle, clockwise, as radians.
         Implemented only for 2D and 3D. For 3D, this is given in XZX Euler
         angles.
     translation : (dim,) ndarray-like, optional
diff --git a/python/cucim/src/cucim/skimage/transform/_warps.py b/python/cucim/src/cucim/skimage/transform/_warps.py
index 88a395591..343f1f0f7 100644
--- a/python/cucim/src/cucim/skimage/transform/_warps.py
+++ b/python/cucim/src/cucim/skimage/transform/_warps.py
@@ -166,9 +166,6 @@ def resize(image, output_shape, order=None, mode='reflect', cval=0, clip=None,
     if clip is None:
         clip = True if order > 1 else False
 
-    # Save input value range for clip
-    img_bounds = (image.min(), image.max()) if clip else None
-
     # Translate modes used by np.pad to those used by scipy.ndimage
     ndi_mode = _to_ndimage_mode(mode)
     if anti_aliasing:
@@ -192,13 +189,15 @@ def resize(image, output_shape, order=None, mode='reflect', cval=0, clip=None,
         _ndi_mode = {'grid-constant': 'constant', 'grid-wrap':'wrap'}.get(ndi_mode, ndi_mode)  # noqa
         # keep ndi.gaussian_filter rather than cucim.skimage.filters.gaussian
         # to avoid undesired dtype coercion
-        image = ndi.gaussian_filter(image, anti_aliasing_sigma, cval=cval,
-                                    mode=_ndi_mode)
+        filtered = ndi.gaussian_filter(image, anti_aliasing_sigma, cval=cval,
+                                       mode=_ndi_mode)
+    else:
+        filtered = image
 
     zoom_factors = [1 / f for f in factors]
-    out = ndi.zoom(image, zoom_factors, order=order, mode=ndi_mode,
+    out = ndi.zoom(filtered, zoom_factors, order=order, mode=ndi_mode,
                    cval=cval, grid_mode=True)
-    _clip_warp_output(img_bounds, out, mode, cval, order, clip)
+    _clip_warp_output(image, out, mode, cval, order, clip)
     return out
 
 
@@ -821,7 +820,7 @@ def _clip_warp_output(input_image, output_image, mode, cval, order, clip):
     cp.clip(output_image, min_val, max_val, out=output_image)
 
 
-def warp(image, inverse_map, map_args={}, output_shape=None, order=None,
+def warp(image, inverse_map, map_args=None, output_shape=None, order=None,
          mode='constant', cval=0., clip=None, preserve_range=False):
     """Warp an image according to a given coordinate transformation.
 
@@ -963,6 +962,9 @@ def warp(image, inverse_map, map_args={}, output_shape=None, order=None,
     >>> warped = warp(cube, coords)
 
     """  # noqa
+    if map_args is None:
+        map_args = {}
+
     if image.size == 0:
         raise ValueError("Cannot warp empty image with dimensions",
                          image.shape)
diff --git a/run b/run
index 35a2e6674..a9e666350 100755
--- a/run
+++ b/run
@@ -234,6 +234,30 @@ retry() {
     return 0
 }
 
+get_arch() {
+    local platform="${1:-$(uname -p)}"
+    local platform_str
+
+    case "${platform}" in
+        amd64|x86_64|x86|linux/amd64)
+            platform_str="x86_64"
+            ;;
+        arm64|aarch64|arm|linux/arm64)
+            platform_str="aarch64"
+            ;;
+    esac
+
+    echo -n "${platform_str}"
+}
+
+is_aarch64() {
+    get_arch | grep -q "aarch64" && return 0 || return 1
+}
+
+is_x86_64() {
+    get_arch | grep -q "x86_64" && return 0 || return 1
+}
+
 #==================================================================================
 # Section: Build
 #==================================================================================
@@ -432,7 +456,22 @@ build_local() {
 
 copy_gds_files_() {
     local root_folder=${1:-${TOP}}
-    local cufile_search="${root_folder}/temp/cuda/include:${root_folder}/temp/cuda/lib64"  #"/usr/local/cuda/include:/usr/local/cuda/lib64 ${PREFIX:-}/include:${PREFIX:-}/lib ${CONDA_PREFIX:-}/include:${CONDA_PREFIX:-}/lib ${root_folder}/temp/cuda/include:${root_folder}/temp/cuda/lib64"
+    local cufile_search="${root_folder}/temp/cuda/include:${root_folder}/temp/cuda/lib64"
+    local targets_dir=""
+    if is_x86_64; then
+      targets_dir="targets/x86_64-linux"
+    elif is_aarch64; then
+      targets_dir="targets/sbsa-linux"
+    else
+      echo "Unknown architecture"
+      exit 1
+    fi
+    if [ ! -z "${CONDA_PREFIX+x}" ]; then
+      cufile_search="${CONDA_PREFIX:-}/include:${CONDA_PREFIX:-}/lib ${CONDA_PREFIX:-}/${targets_dir}/include:${CONDA_PREFIX:-}/${targets_dir}/lib ${cufile_search}"
+    fi
+    if [ ! -z "${PREFIX+x}" ]; then
+      cufile_search="${PREFIX:-}/include:${PREFIX:-}/lib ${PREFIX:-}/${targets_dir}/include:${PREFIX:-}/${targets_dir}/lib ${cufile_search}"
+    fi
     local gds_version=1.0.0
     local candidate
     local cufile_include
@@ -459,8 +498,11 @@ copy_gds_files_() {
         pushd ${temp_tgz_dir}
         run_command wget https://developer.download.nvidia.com/gds/redist/rel-${gds_version}/gds-redistrib-${gds_version}.tgz
         run_command tar xzvf gds-redistrib-${gds_version}.tgz
+        # There are no GDS binaries for aarch64 SBSA yet, but we need cufile.h for building cufile_stub
         run_command cp -P gds-redistrib-${gds_version}/targets/x86_64-linux/include/cufile.h ${root_folder}/temp/cuda/include/
-        run_command cp -P gds-redistrib-${gds_version}/targets/x86_64-linux/lib/* ${root_folder}/temp/cuda/lib64/
+        if is_x86_64; then
+            run_command cp -P gds-redistrib-${gds_version}/targets/x86_64-linux/lib/* ${root_folder}/temp/cuda/lib64/
+        fi
         popd > /dev/null
         run_command rm -r ${temp_tgz_dir}
     else
@@ -486,6 +528,21 @@ copy_nvjpeg_files_() {
     local root_folder=${1:-${TOP}}
     local arch_name="$(get_arch_name_)"
     local nvjpeg_search="${root_folder}/temp/cuda/include:${root_folder}/temp/cuda/lib64"
+    local targets_dir=""
+    if is_x86_64; then
+      targets_dir="targets/x86_64-linux"
+    elif is_aarch64; then
+      targets_dir="targets/sbsa-linux"
+    else
+      echo "Unknown architecture"
+      exit 1
+    fi
+    if [ ! -z "${CONDA_PREFIX+x}" ]; then
+      nvjpeg_search="${CONDA_PREFIX:-}/include:${CONDA_PREFIX:-}/lib ${CONDA_PREFIX:-}/${targets_dir}/include:${CONDA_PREFIX:-}/${targets_dir}/lib ${nvjpeg_search}"
+    fi
+    if [ ! -z "${PREFIX+x}" ]; then
+      nvjpeg_search="${PREFIX:-}/include:${PREFIX:-}/lib ${PREFIX:-}/${targets_dir}/include:${PREFIX:-}/${targets_dir}/lib ${nvjpeg_search}"
+    fi
     local cuda_version="11.6"
     local nvjpeg_version="11-6_11.6.0.55-1"
     local candidate


Test patch:
diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index d5c918a2f..39b8ac83b 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -16,7 +16,7 @@ on:
 jobs:
   conda-python-tests:
     secrets: inherit
-    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-tests.yaml@branch-23.06
+    uses: rapidsai/shared-action-workflows/.github/workflows/conda-python-tests.yaml@branch-23.08
     with:
       build_type: nightly
       branch: ${{ inputs.branch }}
diff --git a/python/cucim/src/cucim/skimage/_shared/testing.py b/python/cucim/src/cucim/skimage/_shared/testing.py
index 8488c100a..2442db08d 100644
--- a/python/cucim/src/cucim/skimage/_shared/testing.py
+++ b/python/cucim/src/cucim/skimage/_shared/testing.py
@@ -15,14 +15,11 @@
 
 have_fetch = True
 try:
-    # scikit-image 0.19
+    # scikit-image >=0.19
     from skimage.data._fetchers import _fetch
 except ImportError:
-    # scikit-image 0.18
-    try:
-        from skimage.data import _fetch
-    except ImportError:
-        have_fetch = False
+    # skip this test if private API changed on scikit-image end
+    have_fetch = False
 
 
 def fetch(data_filename):
diff --git a/python/cucim/src/cucim/skimage/color/tests/test_colorconv.py b/python/cucim/src/cucim/skimage/color/tests/test_colorconv.py
index 0d50cf119..f3bd45a1e 100644
--- a/python/cucim/src/cucim/skimage/color/tests/test_colorconv.py
+++ b/python/cucim/src/cucim/skimage/color/tests/test_colorconv.py
@@ -666,9 +666,13 @@ def test_lab_full_gamut(self):
         a, b = cp.meshgrid(cp.arange(-100, 100), cp.arange(-100, 100))
         L = cp.ones(a.shape)
         lab = cp.dstack((L, a, b))
+        regex = (
+            "Conversion from CIE-LAB to XYZ color space resulted in "
+            "\\d+ negative Z values that have been clipped to zero"
+        )
         for value in [0, 10, 20]:
             lab[:, :, 0] = value
-            with expected_warnings(['Color data out of range']):
+            with pytest.warns(UserWarning, match=regex):
                 lab2xyz(lab)
 
     @pytest.mark.parametrize("channel_axis", [0, 1, -1, -2])
@@ -800,6 +804,17 @@ def test_rgb2yiq_conversion(self):
                         )
         assert_array_almost_equal(yiq, gt, decimal=2)
 
+    @pytest.mark.parametrize("func", [lab2rgb, lab2xyz])
+    def test_warning_stacklevel(self, func):
+        regex = (
+            "Conversion from CIE-LAB.* XYZ.*color space resulted in "
+            "1 negative Z values that have been clipped to zero"
+        )
+        with pytest.warns(UserWarning, match=regex) as messages:
+            func(lab=cp.array([[[0, 0, 300.]]]))
+        assert len(messages) == 1
+        assert messages[0].filename == __file__, "warning points at wrong file"
+
 
 def test_gray2rgb():
     x = cp.asarray([0, 0.5, 1])
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 7087611f0..11f6fb673 100644
--- a/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
+++ b/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
@@ -25,7 +25,8 @@
 
 def test_wrong_source_range():
     im = cp.array([-1, 100], dtype=cp.int8)
-    with pytest.raises(ValueError):
+    match = "Incorrect value for `source_range` argument"
+    with pytest.raises(ValueError, match=match):
         frequencies, bin_centers = exposure.histogram(
             im, source_range="foobar"
         )
@@ -592,6 +593,12 @@ def norm_brightness_err(img1, img2):
     return nbe
 
 
+def test_adapthist_incorrect_kernel_size():
+    img = cp.ones((8, 8), dtype=float)
+    with pytest.raises(ValueError, match="Incorrect value of `kernel_size`"):
+        exposure.equalize_adapthist(img, (3, 3, 3))
+
+
 # Test Gamma Correction
 # =====================
 
diff --git a/python/cucim/src/cucim/skimage/feature/tests/test_canny.py b/python/cucim/src/cucim/skimage/feature/tests/test_canny.py
index 6e06114fc..c7880bb60 100644
--- a/python/cucim/src/cucim/skimage/feature/tests/test_canny.py
+++ b/python/cucim/src/cucim/skimage/feature/tests/test_canny.py
@@ -76,8 +76,7 @@ def test_mask_none(self):
         result2 = feature.canny(cp.zeros((20, 20)), 4, 0, 0)
         assert cp.all(result1 == result2)
 
-    @cp.testing.with_requires("scikit-image>=0.18")
-    @pytest.mark.parametrize('image_dtype', [cp.uint8, cp.int64, cp.float32,
+    @pytest.mark.parametrize('image_dtype', [cp.uint8, cp.int32, cp.float32,
                                              cp.float64])
     def test_use_quantiles(self, image_dtype):
         dtype = cp.dtype(image_dtype)
@@ -146,3 +145,24 @@ def test_dtype(self):
             feature.canny(image_float, 1.0, low, high),
             feature.canny(image_uint8, 1.0, 255 * low, 255 * high)
         )
+
+    def test_full_mask_matches_no_mask(self):
+        """The masked and unmasked algorithms should return the same result.
+
+        """
+        image = cp.array(data.camera())
+
+        for mode in ('constant', 'nearest', 'reflect'):
+            cp.testing.assert_array_equal(
+                feature.canny(image, mode=mode),
+                feature.canny(image, mode=mode,
+                              mask=cp.ones_like(image, dtype=bool))
+            )
+
+    @pytest.mark.parametrize('dtype', (cp.int64, cp.uint64))
+    def test_unsupported_int64(self, dtype):
+        image = cp.zeros((10, 10), dtype=dtype)
+        image[3, 3] = cp.iinfo(dtype).max
+        match = "64-bit or larger integer images are not supported"
+        with pytest.raises(ValueError, match=match):
+            feature.canny(image)
diff --git a/python/cucim/src/cucim/skimage/filters/tests/test_thresholding.py b/python/cucim/src/cucim/skimage/filters/tests/test_thresholding.py
index 50db541a9..5e393a06a 100644
--- a/python/cucim/src/cucim/skimage/filters/tests/test_thresholding.py
+++ b/python/cucim/src/cucim/skimage/filters/tests/test_thresholding.py
@@ -282,20 +282,17 @@ def test_threshold_sauvola_iterable_window_size(self):
         assert_array_equal(ref, out)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_otsu_camera_image():
     camera = util.img_as_ubyte(camerad)
     assert 101 < threshold_otsu(camera) < 103
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_otsu_camera_image_histogram():
     camera = util.img_as_ubyte(camerad)
     hist = histogram(camera.ravel(), 256, source_range="image")
     assert 101 < threshold_otsu(hist=hist) < 103
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_otsu_camera_image_counts():
     camera = util.img_as_ubyte(camerad)
     counts, bin_centers = histogram(camera.ravel(), 256, source_range="image")
@@ -341,7 +338,6 @@ def test_otsu_one_color_image_3d():
     assert threshold_otsu(img) == 1
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_li_camera_image():
     image = util.img_as_ubyte(camerad)
     threshold = threshold_li(image)
@@ -431,20 +427,17 @@ def test_li_pathological_arrays():
     assert cp.all(cp.isfinite(thresholds))
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_yen_camera_image():
     camera = util.img_as_ubyte(camerad)
     assert 145 < threshold_yen(camera) < 147
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_yen_camera_image_histogram():
     camera = util.img_as_ubyte(camerad)
     hist = histogram(camera.ravel(), 256, source_range="image")
     assert 145 < threshold_yen(hist=hist) < 147
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_yen_camera_image_counts():
     camera = util.img_as_ubyte(camerad)
     counts, bin_centers = histogram(camera.ravel(), 256, source_range='image')
@@ -467,7 +460,6 @@ def test_local_even_block_size_error():
         threshold_local(img, block_size=4)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_isodata_camera_image():
     camera = util.img_as_ubyte(camerad)
 
@@ -479,7 +471,6 @@ def test_isodata_camera_image():
     assert_array_equal(threshold_isodata(camera, return_all=True), [102, 103])
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_isodata_camera_image_histogram():
     camera = util.img_as_ubyte(camerad)
     hist = histogram(camera.ravel(), 256, source_range='image')
@@ -487,7 +478,6 @@ def test_isodata_camera_image_histogram():
     assert threshold == 102
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_isodata_camera_image_counts():
     camera = util.img_as_ubyte(camerad)
     counts, bin_centers = histogram(camera.ravel(), 256, source_range='image')
@@ -551,7 +541,6 @@ def test_isodata_moon_image_negative_float():
     # fmt: on
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_threshold_minimum():
     camera = util.img_as_ubyte(camerad)
 
@@ -563,7 +552,6 @@ def test_threshold_minimum():
     assert_array_equal(threshold, 114)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_threshold_minimum_histogram():
     camera = util.img_as_ubyte(camerad)
     hist = histogram(camera.ravel(), 256, source_range='image')
@@ -578,7 +566,6 @@ def test_threshold_minimum_deprecated_max_iter_kwarg():
         threshold_minimum(hist=hist, max_iter=5000)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_threshold_minimum_counts():
     camera = util.img_as_ubyte(camerad)
     counts, bin_centers = histogram(camera.ravel(), 256, source_range='image')
@@ -774,7 +761,6 @@ def test_multiotsu_more_classes_then_values():
         threshold_multiotsu(img, classes=4)
 
 
-# @testing.with_requires("scikit-image>=0.18")
 # @pytest.mark.parametrize(
 #     "thresholding, lower, upper",
 #     [
diff --git a/python/cucim/src/cucim/skimage/metrics/tests/test_simple_metrics.py b/python/cucim/src/cucim/skimage/metrics/tests/test_simple_metrics.py
index 3f9deee0f..8c062a2dc 100644
--- a/python/cucim/src/cucim/skimage/metrics/tests/test_simple_metrics.py
+++ b/python/cucim/src/cucim/skimage/metrics/tests/test_simple_metrics.py
@@ -24,7 +24,6 @@
 
 
 @pytest.mark.parametrize('dtype', [cp.uint8, cp.float32, cp.float64])
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_PSNR_vs_IPOL(dtype):
     """Tests vs. imdiff result from the following IPOL article and code:
     https://www.ipol.im/pub/art/2011/g_lmii/.
diff --git a/python/cucim/src/cucim/skimage/metrics/tests/test_structural_similarity.py b/python/cucim/src/cucim/skimage/metrics/tests/test_structural_similarity.py
index 19049b51a..c11569df9 100644
--- a/python/cucim/src/cucim/skimage/metrics/tests/test_structural_similarity.py
+++ b/python/cucim/src/cucim/skimage/metrics/tests/test_structural_similarity.py
@@ -188,7 +188,6 @@ def test_structural_similarity_multichannel_chelsea():
     assert_equal(structural_similarity(Xc, Xc, channel_axis=-1), 1.0)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_gaussian_structural_similarity_vs_IPOL():
     """Tests vs. imdiff result from the following IPOL article and code:
     https://www.ipol.im/pub/art/2011/g_lmii/.
@@ -209,7 +208,6 @@ def test_gaussian_structural_similarity_vs_IPOL():
 @pytest.mark.parametrize(
     'dtype', [np.uint8, np.int32, np.float16, np.float32, np.float64]
 )
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_mssim_vs_legacy(dtype):
     # check that ssim with default options matches skimage 0.11 result
     mssim_skimage_0pt17 = 0.3674518327910367
diff --git a/python/cucim/src/cucim/skimage/restoration/tests/test_j_invariant.py b/python/cucim/src/cucim/skimage/restoration/tests/test_j_invariant.py
index 88d814cbb..f42bc252f 100644
--- a/python/cucim/src/cucim/skimage/restoration/tests/test_j_invariant.py
+++ b/python/cucim/src/cucim/skimage/restoration/tests/test_j_invariant.py
@@ -1,7 +1,6 @@
 import cupy as cp
 import numpy as np
 import pytest
-import skimage
 from skimage.data import camera, chelsea
 # from cucim.skimage.restoration import denoise_wavelet
 from skimage.restoration import denoise_wavelet
@@ -10,7 +9,7 @@
 from cucim.skimage.data import binary_blobs
 from cucim.skimage.metrics import mean_squared_error as mse
 from cucim.skimage.restoration import calibrate_denoiser, denoise_tv_chambolle
-from cucim.skimage.restoration.j_invariant import _invariant_denoise
+from cucim.skimage.restoration.j_invariant import denoise_invariant
 from cucim.skimage.util import img_as_float, random_noise
 
 test_img = img_as_float(cp.asarray(camera()))
@@ -23,12 +22,6 @@
 
 # TODO: replace with CuPy version once completed
 def _denoise_wavelet(image, rescale_sigma=True, **kwargs):
-
-    if 'channel_axis' in kwargs and skimage.__version__ < '0.19.0':
-        # convert channel_axis to older multichannel kwarg for skimage 0.18.x
-        channel_axis = kwargs.pop('channel_axis')
-        kwargs['multichannel'] = False if channel_axis is None else True
-
     return cp.asarray(
         denoise_wavelet(
             cp.asnumpy(image), rescale_sigma=rescale_sigma, **kwargs
@@ -36,9 +29,9 @@ def _denoise_wavelet(image, rescale_sigma=True, **kwargs):
     )
 
 
-def test_invariant_denoise():
-    # denoised_img = _invariant_denoise(noisy_img, _denoise_wavelet)
-    denoised_img = _invariant_denoise(noisy_img, denoise_tv_chambolle)
+def test_denoise_invariant():
+    # denoised_img = denoise_invariant(noisy_img, _denoise_wavelet)
+    denoised_img = denoise_invariant(noisy_img, denoise_tv_chambolle)
 
     denoised_mse = mse(denoised_img, test_img)
     original_mse = mse(noisy_img, test_img)
@@ -46,8 +39,8 @@ def test_invariant_denoise():
 
 
 @pytest.mark.parametrize('dtype', [cp.float16, cp.float32, cp.float64])
-def test_invariant_denoise_color(dtype):
-    denoised_img_color = _invariant_denoise(
+def test_denoise_invariant_color(dtype):
+    denoised_img_color = denoise_invariant(
         noisy_img_color.astype(dtype),
         _denoise_wavelet,
         denoiser_kwargs=dict(channel_axis=-1),
@@ -59,8 +52,8 @@ def test_invariant_denoise_color(dtype):
     assert denoised_img_color.dtype == _supported_float_type(dtype)
 
 
-def test_invariant_denoise_3d():
-    denoised_img_3d = _invariant_denoise(noisy_img_3d, _denoise_wavelet)
+def test_denoise_invariant_3d():
+    denoised_img_3d = denoise_invariant(noisy_img_3d, _denoise_wavelet)
 
     denoised_mse = mse(denoised_img_3d, test_img_3d)
     original_mse = mse(noisy_img_3d, test_img_3d)
@@ -76,8 +69,8 @@ def test_calibrate_denoiser_extra_output():
         extra_output=True
     )
 
-    all_denoised = [_invariant_denoise(noisy_img, _denoise_wavelet,
-                                       denoiser_kwargs=denoiser_kwargs)
+    all_denoised = [denoise_invariant(noisy_img, _denoise_wavelet,
+                                      denoiser_kwargs=denoiser_kwargs)
                     for denoiser_kwargs in parameters_tested]
 
     ground_truth_losses = [float(mse(img, test_img)) for img in all_denoised]
diff --git a/python/cucim/src/cucim/skimage/restoration/tests/test_restoration.py b/python/cucim/src/cucim/skimage/restoration/tests/test_restoration.py
index 5d5ccd5a9..adbf0833b 100644
--- a/python/cucim/src/cucim/skimage/restoration/tests/test_restoration.py
+++ b/python/cucim/src/cucim/skimage/restoration/tests/test_restoration.py
@@ -1,7 +1,6 @@
 import cupy as cp
 import numpy as np
 import pytest
-from cupy import testing
 from cupyx.scipy import ndimage as ndi
 from scipy import signal
 
@@ -82,8 +81,7 @@ def test_unsupervised_wiener(dtype):
 
     psf = cp.asarray(psf, dtype=dtype)
     data = cp.asarray(data, dtype=dtype)
-    deconvolved, _ = restoration.unsupervised_wiener(data, psf,
-                                                     random_state=seed)
+    deconvolved, _ = restoration.unsupervised_wiener(data, psf, seed=seed)
     float_type = _supported_float_type(dtype)
     assert deconvolved.dtype == float_type
 
@@ -109,7 +107,7 @@ def test_unsupervised_wiener(dtype):
         user_params={"callback": lambda x: None,
                      "max_num_iter": 200,
                      "min_num_iter": 30},
-        random_state=seed,
+        seed=seed,
     )[0]
     assert deconvolved2.real.dtype == float_type
 
@@ -129,14 +127,14 @@ def test_unsupervised_wiener_deprecated_user_param():
     otf = uft.ir2tf(psf, data.shape, is_real=False)
     _, laplacian = uft.laplacian(2, data.shape)
     with expected_warnings(["`max_iter` is a deprecated key",
-                            "`min_iter` is a deprecated key"]):
+                            "`min_iter` is a deprecated key",
+                            "`random_state` is a deprecated argument name"]):
         restoration.unsupervised_wiener(
             data, otf, reg=laplacian, is_real=False,
             user_params={"max_iter": 200, "min_iter": 30}, random_state=5
         )
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_image_shape():
     """Test that shape of output image in deconvolution is same as input.
 
@@ -175,20 +173,8 @@ def test_richardson_lucy():
     cp.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-4)
 
 
-def test_richardson_lucy_deprecated_iterations_kwarg():
-    psf = np.ones((5, 5)) / 25
-    data = signal.convolve2d(cp.asnumpy(test_img), psf, 'same')
-    np.random.seed(0)
-    data += 0.1 * data.std() * np.random.standard_normal(data.shape)
-    data = cp.array(data)
-    psf = cp.array(psf)
-    with expected_warnings(["`iterations` is a deprecated argument"]):
-        restoration.richardson_lucy(data, psf, iterations=5)
-
-
 @pytest.mark.parametrize('dtype_image', [cp.float16, cp.float32, cp.float64])
 @pytest.mark.parametrize('dtype_psf', [cp.float32, cp.float64])
-@testing.with_requires("scikit-image>=0.18")
 def test_richardson_lucy_filtered(dtype_image, dtype_psf):
     if dtype_image == cp.float64:
         atol = 1e-8
diff --git a/python/cucim/src/cucim/skimage/segmentation/tests/test_join.py b/python/cucim/src/cucim/skimage/segmentation/tests/test_join.py
index af528073c..e68df7b8e 100644
--- a/python/cucim/src/cucim/skimage/segmentation/tests/test_join.py
+++ b/python/cucim/src/cucim/skimage/segmentation/tests/test_join.py
@@ -23,6 +23,12 @@ def test_join_segmentations():
                       [0, 5, 3, 2],
                       [4, 5, 5, 3]])
     assert_array_equal(j, j_ref)
+
+    # test correct mapping
+    j, m1, m2 = join_segmentations(s1, s2, return_mapping=True)
+    assert_array_equal(m1[j], s1)
+    assert_array_equal(m2[j], s2)
+
     # fmt: on
     # test correct exception when arrays are different shapes
     s3 = cp.array([[0, 0, 1, 1], [0, 2, 2, 1]])
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 2b2891542..d9aee75db 100644
--- a/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
+++ b/python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
@@ -269,6 +269,38 @@ def test_affine_init():
                               AffineTransform(scale=(0.5, 0.5)).scale)
 
 
+@pytest.mark.parametrize("xp", [np, cp])
+def test_affine_shear(xp):
+    shear = 0.1
+    # expected horizontal shear transform
+    cx = -np.tan(shear)
+    # fmt: off
+    expected = xp.array([
+        [1, cx, 0],
+        [0,  1, 0],  # noqa: E241
+        [0,  0, 1],  # noqa: E241
+    ])
+    # fmt: on
+
+    tform = AffineTransform(shear=shear, xp=xp)
+    xp.testing.assert_array_almost_equal(tform.params, expected)
+
+    shear = (1.2, 0.8)
+    # expected x, y shear transform
+    cx = -np.tan(shear[0])
+    cy = -np.tan(shear[1])
+    # fmt: off
+    expected = xp.array([
+        [ 1, cx, 0],  # noqa: E201
+        [cy,  1, 0],  # noqa: E241
+        [ 0,  0, 1],  # noqa: E241, E201
+    ])
+    # fmt: on
+
+    tform = AffineTransform(shear=shear, xp=xp)
+    xp.testing.assert_array_almost_equal(tform.params, expected)
+
+
 def test_piecewise_affine():
     tform = PiecewiseAffineTransform()
     assert tform.estimate(SRC, DST)
diff --git a/python/cucim/src/cucim/skimage/transform/tests/test_warps.py b/python/cucim/src/cucim/skimage/transform/tests/test_warps.py
index de7bbb06c..dbb47738e 100644
--- a/python/cucim/src/cucim/skimage/transform/tests/test_warps.py
+++ b/python/cucim/src/cucim/skimage/transform/tests/test_warps.py
@@ -472,10 +472,12 @@ def test_resize_clip(order, preserve_range, anti_aliasing, dtype):
     x = cp.ones((5, 5), dtype=dtype)
     if dtype == cp.uint8:
         x *= 255
+    else:
+        x[0, 0] = cp.NaN
     resized = resize(x, (3, 3), order=order, preserve_range=preserve_range,
                      anti_aliasing=anti_aliasing)
 
-    assert abs(float(resized.max()) - expected_max) < 1e-14
+    assert abs(float(cp.nanmax(resized)) - expected_max) < 1e-14
 
 
 @pytest.mark.parametrize('dtype', [cp.float16, cp.float32, cp.float64])
diff --git a/python/cucim/src/cucim/skimage/util/tests/test_random_noise.py b/python/cucim/src/cucim/skimage/util/tests/test_random_noise.py
index 0bbbf6a30..ae601e4bb 100644
--- a/python/cucim/src/cucim/skimage/util/tests/test_random_noise.py
+++ b/python/cucim/src/cucim/skimage/util/tests/test_random_noise.py
@@ -171,7 +171,6 @@ def test_clip_poisson():
     assert (cam_poisson2.max() > 1.3) and (cam_poisson2.min() == -1.0)
 
 
-@cp.testing.with_requires("scikit-image>=0.18")
 def test_clip_gaussian():
     seed = 42
     data = camerad  # 512x512 grayscale uint8
diff --git a/python/cucim/tests/performance/clara/test_read_region_memory_usage.py b/python/cucim/tests/performance/clara/test_read_region_memory_usage.py
index 676168937..0ae507ed5 100644
--- a/python/cucim/tests/performance/clara/test_read_region_memory_usage.py
+++ b/python/cucim/tests/performance/clara/test_read_region_memory_usage.py
@@ -30,7 +30,7 @@ def test_read_region_cuda_memleak(testimg_tiff_stripe_4096x4096_256_jpeg):
     gpu = gpus[0]
     mem_usage_history = [gpu.memoryUsed]
 
-    for i in range(5):
+    for i in range(10):
         _ = img.read_region(device='cuda')
         gpus = GPUtil.getGPUs()
         gpu = gpus[0]
@@ -40,8 +40,10 @@ def test_read_region_cuda_memleak(testimg_tiff_stripe_4096x4096_256_jpeg):
 
     # The difference in memory usage should be less than 180MB.
     # Note: Since we cannot measure GPU memory usage for a process,
-    #       we use a rough number (experimentally measured).
-    assert mem_usage_history[4] - mem_usage_history[1] < 180.0
+    #       we use a rough number.
+    #       (experimentally measured, assuming that each image load
+    #        consumes around 50MB of GPU memory).
+    assert mem_usage_history[5] - mem_usage_history[9] < 180.0
 
 
 def test_read_region_cpu_memleak(testimg_tiff_stripe_4096x4096_256):