Repository: rapidsai/cucim (language: python)

Issue:
parameterization of `filters` and `features` benchmarks (v2)
Update of #278 for branch-22.08 with minor additional fixes. I have tested this locally and it seems to be operating as expected.


[DOC] audit See Also sections of docs and update deprecation warning strings
## Report incorrect documentation

### Location of incorrect documentation

See the ["See Also" section for `random_walker`](https://docs.rapids.ai/api/cucim/stable/api.html#module-cucim.skimage.segmentation). It mentions `skimage.morphology.watershed` but we do not currently have `watershed` in cuCIM. Similar cases can be found elsewhere as well.

### Describe the problems or issues found in the documentation

Some current **cucim.skimage** **See Also** sections (as copied from scikit-image docstrings) can be misleading when the referred function is not implemented in cuCIM. 

Similarly, docstrings notes regarding deprecations are sometime present and mention scikit-image version numbers (e.g. 0.19) while cuCIM uses a different version numbering scheme. 

### Suggested fix for documentation
We should audit all docstrings and remove references to unimplemented functions (or at least mark them as not implemented)

For deprecations, we should update the message to indicate that the argument is deprecated upstream in scikit-image and will eventually be removed here as well once the upstream deprecation is complete.

Add segmentation with the Chan-Vese active contours method
closes #335

This PR adds a function for segmentation based on active-contours. Test cases are directly from scikit-image, but the code itself was heavily refactored to make extensive use of kernel fusion via `cupy.fuse()`. I will add benchmark results here soon.

It should be reviewed after #318 is merged. The new commits only start from 502f41a



Discussion / hints:
Looks like there are some conflicts. Could you please take a look at these @grlee77? 🙂
@gpucibot merge
> Looks like there are some conflicts. Could you please take a look at these @grlee77? slightly_smiling_face

Resolved. should be good to go now. No point to wait on CI as it doesn't run these benchmarks.
Yeah it won't merge until CI passes. So we just need to wait unfortunately. On the bright side it should merge as soon as CI passes

This shows a huge speedup relative to scikit-image as image size increases (>200x):

shape | Duration (CPU) | Duration (GPU) | Acceleration
------|----------------|----------------|-------------
(256, 256) | 0.679 +/- 0.014 | 0.209 +/- 0.005 | 3.251
(512, 512) | 3.196 +/- 0.014 | 0.208 +/- 0.002 | 15.397
(1024, 1024) | 14.350 +/- 0.000 | 0.206 +/- 0.001 | 69.603
(2048, 2048) | 68.054 +/- 0.000 | 0.426 +/- 0.000 | 159.710
(3072, 3072) | 199.156 +/- 0.000 | 0.874 +/- 0.000 | 227.748

<details>
<summary>benchmark script</summary>

```Python
import cupy as cp
import matplotlib.pyplot as plt
from cupyx.profiler import benchmark
from skimage import data, segmentation

from cucim.skimage import img_as_float
from cucim.skimage.segmentation import chan_vese


_image = img_as_float(cp.array(data.camera()))

print("shape | Duration (CPU) | Duration (GPU) | Acceleration")
print("------|----------------|----------------|-------------")
for n in [0.5, 1, 2, 4, 6]:
    if n == 0.5:
        image = _image[128:-128, 128:-128]
    else:
        image = cp.tile(_image, (n, n))
    image_cpu = cp.asnumpy(image)
    kwargs = dict(mu=0.25, lambda1=1, lambda2=1, tol=1e-3, max_num_iter=200, dt=0.5, init_level_set="checkerboard", extended_output=True)
    perf_cpu = benchmark(segmentation.chan_vese, (image_cpu, ), kwargs=kwargs, n_warmup=0, n_repeat=10, max_duration=5)
    t_cpu = perf_cpu.gpu_times
    perf_gpu = benchmark(chan_vese, (image, ), kwargs=kwargs, n_warmup=4, n_repeat=100, max_duration=5)
    t_gpu = perf_gpu.gpu_times
    print(f"{image.shape} | {t_cpu.mean():0.3f} +/- {t_cpu.std():0.3f} | {t_gpu.mean():0.3f} +/- {t_gpu.std():0.3f} | {t_cpu.mean() / t_gpu.mean():0.3f}")
```

</details>


Fix patch:
diff --git a/3rdparty/LICENSE.pba+ b/3rdparty/LICENSE.pba+
new file mode 100644
index 000000000..9d0b4030a
--- /dev/null
+++ b/3rdparty/LICENSE.pba+
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2019 School of Computing, National University of Singapore
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8a7ca72dc..7eb977a03 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,7 @@ Please see https://github.com/rapidsai/cucim/releases/tag/v22.08.00a for the lat
 ## 🚀 New Features
 
 - add missing `cucim.skimage.segmentation.clear_border` function ([#267](https://github.com/rapidsai/cucim/pull/267)) [@grlee77](https://github.com/grlee77)
+- add `cucim.core.operations.color.stain_extraction_pca` and `cucim.core.operations.color.normalize_colors_pca` for digital pathology H&E stain extraction and normalization ([#273](https://github.com/rapidsai/cucim/pull/273)) [@grlee77](https://github.com/grlee77), [@drbeh](https://github.com/drbeh)
 
 ## 🛠️ Improvements
 
@@ -31,6 +32,7 @@ Please see https://github.com/rapidsai/cucim/releases/tag/v22.08.00a for the lat
 - Promote small integer types to single rather than double precision ([#278](https://github.com/rapidsai/cucim/pull/278)) [@grlee77](https://github.com/grlee77)
 - improve efficiency of histogram-based thresholding functions ([#276](https://github.com/rapidsai/cucim/pull/276)) [@grlee77](https://github.com/grlee77)
 - Remove unused dependencies in GPU tests job ([#268](https://github.com/rapidsai/cucim/pull/268)) [@Ethyling](https://github.com/Ethyling)
+- Enable footprint decomposition for morphology ([#274](https://github.com/rapidsai/cucim/pull/274)) [@grlee77](https://github.com/grlee77)
 - Use conda compilers ([#232](https://github.com/rapidsai/cucim/pull/232)) [@Ethyling](https://github.com/Ethyling)
 - Build packages using mambabuild ([#216](https://github.com/rapidsai/cucim/pull/216)) [@Ethyling](https://github.com/Ethyling)
 
diff --git a/LICENSE-3rdparty.md b/LICENSE-3rdparty.md
index 0a3789af7..984d0e952 100644
--- a/LICENSE-3rdparty.md
+++ b/LICENSE-3rdparty.md
@@ -281,3 +281,9 @@ StainTools
   - https://github.com/Peter554/StainTools/blob/master/LICENSE.txt
 - Copyright: Peter Byfield
 - Usage: reference for stain color normalization algorithm
+
+PBA+
+- License: MIT License
+  - https://github.com/orzzzjq/Parallel-Banding-Algorithm-plus/blob/master/LICENSE
+- Copyright: School of Computing, National University of Singapore
+- Usage: PBA+ is used to implement the Euclidean distance transform.
diff --git a/README.md b/README.md
index 3564282c9..612fda219 100644
--- a/README.md
+++ b/README.md
@@ -1,14 +1,43 @@
 # <div align="left"><img src="https://rapids.ai/assets/images/rapids_logo.png" width="90px"/>&nbsp;cuCIM</div>
 
+[RAPIDS](https://rapids.ai) cuCIM is an open-source, accelerated computer vision and image processing software library for multidimensional images used in biomedical, geospatial, material and life science, and remote sensing use cases.
 
-[RAPIDS](https://rapids.ai) cuCIM is an extensible toolkit designed to provide GPU accelerated I/O, computer vision & image processing primitives for N-Dimensional images with a focus on biomedical imaging.
+cuCIM offers:
+
+- Enhanced Image Processing Capabilities for large and n-dimensional tag image file format (TIFF) files
+- Accelerated performance through Graphics Processing Unit (GPU)-based image processing and computer vision primitives
+- A Straightforward Pythonic Interface with Matching Application Programming Interface (API) for Openslide
+
+cuCIM supports the following formats:
+
+- Aperio ScanScope Virtual Slide (SVS)
+- Philips TIFF
+- Generic Tiled, Multi-resolution RGB TIFF files with the following compression schemes:
+  - No Compression
+  - JPEG
+  - JPEG2000
+  - Lempel-Ziv-Welch (LZW)
+  - Deflate
 
 **NOTE:** For the latest stable [README.md](https://github.com/rapidsai/cucim/blob/main/README.md) ensure you are on the `main` branch.
 
-- [GTC 2021 cuCIM: A GPU Image I/O and Processing Toolkit [S32194]](https://www.nvidia.com/en-us/on-demand/search/?facet.mimetype[]=event%20session&layout=list&page=1&q=cucim&sort=date)
-  - [video](https://www.nvidia.com/en-us/on-demand/session/gtcspring21-s32194/)
+- [GTC 2022 Accelerating Storage IO to GPUs with Magnum IO [S41347]](https://events.rainfocus.com/widget/nvidia/gtcspring2022/sessioncatalog/session/1634960000577001Etxp)
+  - cuCIM's GDS API examples: <https://github.com/NVIDIA/MagnumIO/tree/main/gds/readers/cucim-gds>
 - [SciPy 2021 cuCIM - A GPU image I/O and processing library](https://www.scipy2021.scipy.org/)
   - [video](https://youtu.be/G46kOOM9xbQ)
+- [GTC 2021 cuCIM: A GPU Image I/O and Processing Toolkit [S32194]](https://www.nvidia.com/en-us/on-demand/search/?facet.mimetype[]=event%20session&layout=list&page=1&q=cucim&sort=date)
+  - [video](https://www.nvidia.com/en-us/on-demand/session/gtcspring21-s32194/)
+
+**Blogs**
+- [Enhanced Image Analysis with Multidimensional Image Processing](https://developer.nvidia.com/blog/enhanced-image-analysis-with-multidimensional-image-processing/)
+- [Accelerating Scikit-Image API with cuCIM: n-Dimensional Image Processing and IO on GPUs](https://developer.nvidia.com/blog/cucim-rapid-n-dimensional-image-processing-and-i-o-on-gpus/)
+- [Accelerating Digital Pathology Pipelines with NVIDIA Clara™ Deploy](https://developer.nvidia.com/blog/accelerating-digital-pathology-pipelines-with-nvidia-clara-deploy-2/)
+
+**Webinars**
+
+- [cuCIM: a GPU Image IO and Processing Library](https://www.youtube.com/watch?v=G46kOOM9xbQ)
+
+**[Documentation](https://docs.rapids.ai/api/cucim/stable)**
 
 **Release notes** are available on our [wiki page](https://github.com/rapidsai/cucim/wiki/Release-Notes).
 
diff --git a/benchmarks/skimage/cucim_color_bench.py b/benchmarks/skimage/cucim_color_bench.py
index d253eb0dc..87955dd0a 100644
--- a/benchmarks/skimage/cucim_color_bench.py
+++ b/benchmarks/skimage/cucim_color_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import os
 import pickle
 
@@ -14,6 +15,7 @@
 
 from _image_bench import ImageBench
 
+func_name_choices = ['convert_colorspace', 'rgb2hed', 'hed2rgb', 'lab2lch', 'lch2lab', 'xyz2lab', 'lab2xyz', 'rgba2rgb', 'label2rgb']
 
 class ColorBench(ImageBench):
     def set_args(self, dtype):
@@ -47,6 +49,7 @@ def __init__(
         index_str=None,
         module_cpu=scipy.ndimage,
         module_gpu=cupyx.scipy.ndimage,
+        run_cpu=True,
     ):
         self.contiguous_labels = contiguous_labels
         super().__init__(
@@ -58,6 +61,7 @@ def __init__(
             index_str=index_str,
             module_cpu=module_cpu,
             module_gpu=module_gpu,
+            run_cpu=run_cpu,
         )
 
     def set_args(self, dtype):
@@ -70,14 +74,12 @@ def set_args(self, dtype):
             ],
             dtype=int,
         )
-        tiling = tuple(s // a_s for s, a_s in zip(shape, a.shape))
+        tiling = tuple(s // a_s for s, a_s in zip(self.shape, a.shape))
         if self.contiguous_labels:
             label = np.kron(a, np.ones(tiling, dtype=a.dtype))
         else:
             label = np.tile(a, tiling)
         labeld = cp.asarray(label)
-        if self.shape[-1] != 3:
-            raise ValueError("shape must be 3 on the last axis")
         imaged = cupy.testing.shaped_random(labeld.shape, xp=cp, dtype=dtype, scale=1.0)
         image = cp.asnumpy(imaged)
         self.args_cpu = (
@@ -90,96 +92,125 @@ def set_args(self, dtype):
         )
 
 
-pfile = "cucim_color_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
-all_colorspaces = True
+def main(args):
+
+    pfile = "cucim_color_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    all_colorspaces = False
 
-for shape in [(256, 256, 3), (3840, 2160, 3), (192, 192, 192, 3)]:
     ndim = len(shape)
 
-    if all_colorspaces:
-        color_spaces = ["RGB", "HSV", "RGB CIE", "XYZ", "YUV", "YIQ", "YPbPr", "YCbCr", "YDbDr"]
-    else:
-        color_spaces = ["RGB", "HSV", "YUV", "XYZ"]
-    for fromspace in color_spaces:
-        for tospace in color_spaces:
-            if fromspace == tospace:
-                continue
+    for function_name in func_name_choices:
+
+        if function_name != args.func_name:
+            continue
+
+        if function_name == 'convert_colorspace':
+            if all_colorspaces:
+                color_spaces = ["RGB", "HSV", "RGB CIE", "XYZ", "YUV", "YIQ", "YPbPr", "YCbCr", "YDbDr"]
+            else:
+                color_spaces = ["RGB", "HSV", "YUV", "XYZ"]
+            for fromspace in color_spaces:
+                for tospace in color_spaces:
+                    if fromspace == tospace:
+                        continue
+
+                    B = ColorBench(
+                        function_name="convert_colorspace",
+                        shape=shape + (3,),
+                        dtypes=dtypes,
+                        fixed_kwargs=dict(fromspace=fromspace, tospace=tospace),
+                        var_kwargs={},
+                        index_str=f"{fromspace.lower()}2{tospace.lower()}",
+                        module_cpu=skimage.color,
+                        module_gpu=cucim.skimage.color,
+                        run_cpu=run_cpu,
+                    )
+                    results = B.run_benchmark(duration=args.duration)
+                    all_results = pd.concat([all_results, results["full"]])
+
+        elif function_name == 'rgba2rgb':
+            B = RGBABench(
+                function_name="rgba2rgb",
+                shape=shape[:-1] + (4,),
+                dtypes=dtypes,
+                fixed_kwargs={},
+                var_kwargs={},
+                module_cpu=skimage.color,
+                module_gpu=cucim.skimage.color,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
 
+        elif function_name == 'label2rgb':
+
+            for contiguous_labels in [True, False]:
+                if contiguous_labels:
+                    index_str = "contiguous"
+                else:
+                    index_str = None
+                B = LabelBench(
+                    function_name="label2rgb",
+                    shape=shape,
+                    dtypes=dtypes,
+                    contiguous_labels=contiguous_labels,
+                    index_str=index_str,
+                    fixed_kwargs=dict(bg_label=0),
+                    var_kwargs=dict(kind=["avg", "overlay"]),
+                    module_cpu=skimage.color,
+                    module_gpu=cucim.skimage.color,
+                    run_cpu=run_cpu,
+                )
+                results = B.run_benchmark(duration=args.duration)
+                all_results = pd.concat([all_results, results["full"]])
+
+        elif function_name in [
+            'rgb2hed', 'hed2rgb', 'lab2lch', 'lch2lab', 'xyz2lab', 'lab2xyz'
+        ]:
             B = ColorBench(
-                function_name="convert_colorspace",
-                shape=shape,
+                function_name=function_name,
+                shape=shape + (3,),
                 dtypes=dtypes,
-                fixed_kwargs=dict(fromspace=fromspace, tospace=tospace),
+                fixed_kwargs={},
                 var_kwargs={},
-                index_str=f"{fromspace.lower()}2{tospace.lower()}",
                 module_cpu=skimage.color,
                 module_gpu=cucim.skimage.color,
+                run_cpu=run_cpu,
             )
-            results = B.run_benchmark(duration=1)
+            results = B.run_benchmark(duration=args.duration)
             all_results = pd.concat([all_results, results["full"]])
 
-    # rgb2hed and hed2rgb test combine_stains and separate_stains and all other
-    # stains should have equivalent performance.
-    #
-    # Probably not necessary to benchmark combine_stains and separate_stains
-    # e.g.
-    #    ihc_hdx = separate_stains(ihc, hdx_from_rgb)
-    #    ihc = combine_stains(ihc_hdx, rgb_from_hdx)
-    #
-
-    for fname in ["rgb2hed", "hed2rgb", "lab2lch", "lch2lab", "xyz2lab",
-                  "lab2xyz"]:
-        B = ColorBench(
-            function_name=fname,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs={},
-            var_kwargs={},
-            module_cpu=skimage.color,
-            module_gpu=cucim.skimage.color,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-    B = RGBABench(
-        function_name="rgba2rgb",
-        shape=shape[:-1] + (4,),
-        dtypes=dtypes,
-        fixed_kwargs={},
-        var_kwargs={},
-        module_cpu=skimage.color,
-        module_gpu=cucim.skimage.color,
-    )
-    results = B.run_benchmark(duration=1)
-    all_results = pd.concat([all_results, results["full"]])
-
-    for contiguous_labels in [True, False]:
-        if contiguous_labels:
-            index_str = "contiguous"
-        else:
-            index_str = None
-        B = LabelBench(
-            function_name="label2rgb",
-            shape=shape,
-            dtypes=dtypes,
-            contiguous_labels=contiguous_labels,
-            index_str=index_str,
-            fixed_kwargs=dict(bg_label=0),
-            var_kwargs=dict(kind=["avg", "overlay"]),
-            module_cpu=skimage.color,
-            module_gpu=cucim.skimage.color,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM color conversion functions')
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_exposure_bench.py b/benchmarks/skimage/cucim_exposure_bench.py
index 7293afc18..d0fb4dbda 100644
--- a/benchmarks/skimage/cucim_exposure_bench.py
+++ b/benchmarks/skimage/cucim_exposure_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import os
 import pickle
 
@@ -27,6 +28,7 @@ def set_args(self, dtype):
 
 class MatchHistogramBench(ImageBench):
     def set_args(self, dtype):
+
         if np.dtype(dtype).kind in "iu":
             scale = 256
         else:
@@ -39,95 +41,101 @@ def set_args(self, dtype):
         self.args_gpu = (imaged, imaged2)
 
 
-pfile = "cucim_exposure_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.uint8, np.float32]
-
-exposure_config = {
-    "equalize_adapthist": dict(
-        fixed_kwargs=dict(clip_limit=0.01, nbins=256),
-        variable_kwargs=dict(),
-        color_required=False,
-        grayscale_only=False,
-        dtypes=None,
-        shapes=None,
-    ),
-    "histogram": dict(
-        fixed_kwargs=dict(source_range="image"),
-        variable_kwargs=dict(nbins=[16, 256], normalize=[True, False]),
-        color_required=False,
-        grayscale_only=True,
-        dtypes=None,
-        shapes=None,
-    ),
-}
-
-for function_name, fixed_kwargs, var_kwargs, allow_color in [
-    ("equalize_adapthist", dict(clip_limit=0.01, nbins=256), dict(), True),
-    (
-        "histogram",
-        dict(source_range="image"),
-        dict(nbins=[16, 256], normalize=[True, False]),
-        False,
-    ),
-    ("cumulative_distribution", dict(), dict(nbins=[16, 256]), False),
-    ("equalize_hist", dict(mask=None), dict(nbins=[16, 256]), False),
-    ("rescale_intensity", dict(in_range="image", out_range="dtype"), dict(), False),
-    ("adjust_gamma", dict(), dict(), False),
-    ("adjust_log", dict(), dict(), False),
-    ("adjust_sigmoid", dict(), dict(inv=[False, True]), False),
-    ("is_low_contrast", dict(), dict(), False),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
-        ndim = len(shape)
-        if shape[-1] == 3 and not allow_color:
+def main(args):
+
+    pfile = "cucim_exposure_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color in [
+        ("equalize_adapthist", dict(clip_limit=0.01, nbins=256), dict(), True),
+        (
+            "histogram",
+            dict(source_range="image"),
+            dict(nbins=[16, 256], normalize=[True, False]),
+            False,
+        ),
+        ("cumulative_distribution", dict(), dict(nbins=[16, 256]), False),
+        ("equalize_hist", dict(mask=None), dict(nbins=[16, 256]), False),
+        ("rescale_intensity", dict(in_range="image", out_range="dtype"), dict(), False),
+        ("adjust_gamma", dict(), dict(), False),
+        ("adjust_log", dict(), dict(), False),
+        ("adjust_sigmoid", dict(), dict(inv=[False, True]), False),
+        ("is_low_contrast", dict(), dict(), False),
+    ]:
+
+        if function_name != args.func_name:
             continue
 
-        if function_name == "equalize_adapthist":
-            #  TODO: fix equalize_adapthist for size (3840, 2160) and kernel_size = [16, 16]
-            size_factors = [4, 8, 16]
-            kernel_sizes = []
-            for size_factor in size_factors:
-                kernel_sizes.append([max(s // size_factor, 1) for s in shape if s != 3])
-            var_kwargs.update(dict(kernel_size=kernel_sizes))
-
-        B = ExposureBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.exposure,
-            module_gpu=cucim.skimage.exposure,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
-    ndim = len(shape)
-
-    channel_axis = -1 if shape[-1] in [3, 4] else None
-
-    B = MatchHistogramBench(
-        function_name="match_histograms",
-        shape=shape,
-        dtypes=dtypes,
-        fixed_kwargs=dict(channel_axis=channel_axis),
-        var_kwargs=dict(),
-        module_cpu=skimage.exposure,
-        module_gpu=cucim.skimage.exposure,
-    )
-    results = B.run_benchmark(duration=1)
-    all_results = pd.concat([all_results, results["full"]])
-
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+        ndim = len(shape)
+
+        if function_name == 'match_histograms':
+
+            channel_axis = -1 if shape[-1] in [3, 4] else None
+
+            B = MatchHistogramBench(
+                function_name="match_histograms",
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=dict(channel_axis=channel_axis),
+                var_kwargs=dict(),
+                module_cpu=skimage.exposure,
+                module_gpu=cucim.skimage.exposure,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+        else:
+
+            if shape[-1] == 3 and not allow_color:
+                continue
+
+            if function_name == "equalize_adapthist":
+                #  TODO: fix equalize_adapthist for size (3840, 2160) and kernel_size = [16, 16]
+                size_factors = [4, 8, 16]
+                kernel_sizes = []
+                for size_factor in size_factors:
+                    kernel_sizes.append([max(s // size_factor, 1) for s in shape if s != 3])
+                var_kwargs.update(dict(kernel_size=kernel_sizes))
+
+            B = ExposureBench(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.exposure,
+                module_gpu=cucim.skimage.exposure,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    with open(fbase + ".md", "wt") as f:
+        f.write(all_results.to_markdown())
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM exposure functions')
+    func_name_choices = ['equalize_adapthist', 'cumulative_distribution', 'equalize_hist', 'rescale_intensity', 'adjust_gamma', 'adjust_log', 'adjust_sigmoid', 'is_low_contrast', 'match_histograms']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices=dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices=func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_feature_bench.py b/benchmarks/skimage/cucim_feature_bench.py
index 396d3a8d3..05f61493f 100644
--- a/benchmarks/skimage/cucim_feature_bench.py
+++ b/benchmarks/skimage/cucim_feature_bench.py
@@ -127,13 +127,8 @@ def main(args):
     fbase = os.path.splitext(pfile)[0]
     all_results.to_csv(fbase + ".csv")
     all_results.to_pickle(pfile)
-    try:
-        import tabulate
-
-        with open(fbase + ".md", "wt") as f:
-            f.write(all_results.to_markdown())
-    except ImportError:
-        pass
+    with open(fbase + ".md", "wt") as f:
+        f.write(all_results.to_markdown())
 
 
 if __name__ == '__main__':
diff --git a/benchmarks/skimage/cucim_filters_bench.py b/benchmarks/skimage/cucim_filters_bench.py
index 0ee8c51a6..e12a9a9e8 100644
--- a/benchmarks/skimage/cucim_filters_bench.py
+++ b/benchmarks/skimage/cucim_filters_bench.py
@@ -161,19 +161,14 @@ def main(args):
             module_gpu=cucim.skimage.filters,
             run_cpu=not args.no_cpu,
         )
-        results = B.run_benchmark(duration=1)
+        results = B.run_benchmark(duration=args.duration)
         all_results = pd.concat([all_results, results["full"]])
 
     fbase = os.path.splitext(pfile)[0]
     all_results.to_csv(fbase + ".csv")
     all_results.to_pickle(pfile)
-    try:
-        import tabulate
-
-        with open(fbase + ".md", "wt") as f:
-            f.write(all_results.to_markdown())
-    except ImportError:
-        pass
+    with open(fbase + ".md", "wt") as f:
+        f.write(all_results.to_markdown())
 
 
 if __name__ == '__main__':
diff --git a/benchmarks/skimage/cucim_measure_bench.py b/benchmarks/skimage/cucim_measure_bench.py
index 5047798ff..c5f207b64 100644
--- a/benchmarks/skimage/cucim_measure_bench.py
+++ b/benchmarks/skimage/cucim_measure_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import math
 import os
 import pickle
@@ -25,6 +26,7 @@ def __init__(
         index_str=None,
         module_cpu=skimage.measure,
         module_gpu=cucim.skimage.measure,
+        run_cpu=True,
     ):
 
         self.contiguous_labels = contiguous_labels
@@ -38,6 +40,7 @@ def __init__(
             index_str=index_str,
             module_cpu=module_cpu,
             module_gpu=module_gpu,
+            run_cpu=run_cpu,
         )
 
     def set_args(self, dtype):
@@ -49,7 +52,7 @@ def set_args(self, dtype):
                 [0, 0, 0, 0, 0, 5, 0, 0],
             ]
         )
-        tiling = tuple(s // a_s for s, a_s in zip(shape, a.shape))
+        tiling = tuple(s // a_s for s, a_s in zip(self.shape, a.shape))
         if self.contiguous_labels:
             image = np.kron(a, np.ones(tiling, dtype=a.dtype))
         else:
@@ -71,6 +74,7 @@ def __init__(
         index_str=None,
         module_cpu=skimage.measure,
         module_gpu=cucim.skimage.measure,
+        run_cpu=True,
     ):
 
         self.contiguous_labels = contiguous_labels
@@ -84,6 +88,7 @@ def __init__(
             index_str=index_str,
             module_cpu=module_cpu,
             module_gpu=module_gpu,
+            run_cpu=run_cpu,
         )
 
     def set_args(self, dtype):
@@ -95,7 +100,7 @@ def set_args(self, dtype):
                 [0, 0, 0, 0, 0, 5, 0, 0],
             ]
         )
-        tiling = tuple(s // a_s for s, a_s in zip(shape, a.shape))
+        tiling = tuple(s // a_s for s, a_s in zip(self.shape, a.shape))
         if self.contiguous_labels:
             image = np.kron(a, np.ones(tiling, dtype=a.dtype))
         else:
@@ -127,32 +132,70 @@ def set_args(self, dtype):
         self.args_gpu = (imaged,)
 
 
-pfile = "cucim_measure_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _gaussian.py
-    (
-        "label",
-        dict(return_num=False, background=0),
-        dict(connectivity=[1, 2]),
-        False,
-        True,
-    ),
-    # regionprops.py
-    ("regionprops", dict(), dict(), False, True),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+def main(args):
+
+    pfile = "cucim_measure_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _gaussian.py
+        (
+            "label",
+            dict(return_num=False, background=0),
+            dict(connectivity=[1, 2]),
+            False,
+            True,
+        ),
+        # regionprops.py
+        ("regionprops", dict(), dict(), False, True),
+        # _moments.py
+        ("moments", dict(), dict(order=[1, 2, 3, 4]), False, False),
+        ("moments_central", dict(), dict(order=[1, 2, 3]), False, True),
+        # omited from benchmarks (only tiny arrays): moments_normalized, moments_hu
+        ("centroid", dict(), dict(), False, True),
+        ("inertia_tensor", dict(), dict(), False, True),
+        ("inertia_tensor_eigvals", dict(), dict(), False, True),
+        # _polygon.py
+        # TODO: approximate_polygon, subdivide_polygon
+        # block.py
+        (
+            "block_reduce",
+            dict(),
+            dict(
+                func=[
+                    cp.sum,
+                ]
+            ),
+            True,
+            True,
+        ),  # variable block_size configured below
+        # entropy.py
+        ("shannon_entropy", dict(base=2), dict(), True, True),
+        # profile.py
+        (
+            "profile_line",
+            dict(src=(5, 7)),
+            dict(reduce_func=[cp.mean], linewidth=[1, 2, 4], order=[1, 3]),
+            True,
+            False,
+        ),  # variable block_size configured below
+    ]:
+
+        if function_name != args.func_name:
+            continue
 
         ndim = len(shape)
         if not allow_nd:
-            if allow_color:
+            if not allow_color:
                 if ndim > 2:
                     continue
             else:
@@ -161,118 +204,89 @@ def set_args(self, dtype):
         if shape[-1] == 3 and not allow_color:
             continue
 
-        Tester = LabelBench if function_name == "label" else RegionpropsBench
+        if function_name in ['label', 'regionprops']:
+
+            Tester = LabelBench if function_name == "label" else RegionpropsBench
+
+            for contiguous_labels in [True, False]:
+                if contiguous_labels:
+                    index_str = f"contiguous"
+                else:
+                    index_str = None
+                B = Tester(
+                    function_name=function_name,
+                    shape=shape,
+                    dtypes=dtypes,
+                    contiguous_labels=contiguous_labels,
+                    index_str=index_str,
+                    fixed_kwargs=fixed_kwargs,
+                    var_kwargs=var_kwargs,
+                    module_cpu=skimage.measure,
+                    module_gpu=cucim.skimage.measure,
+                    run_cpu=run_cpu,
+                )
+                results = B.run_benchmark(duration=args.duration)
+                all_results = pd.concat([all_results, results["full"]])
+        else:
+
 
-        for contiguous_labels in [True, False]:
-            if contiguous_labels:
-                index_str = f"contiguous"
-            else:
-                index_str = None
-            B = Tester(
+            if function_name == "gabor" and np.prod(shape) > 1000000:
+                # avoid cases that are too slow on the CPU
+                var_kwargs["frequency"] = [f for f in var_kwargs["frequency"] if f >= 0.1]
+
+            if function_name == "block_reduce":
+                ndim = len(shape)
+                if shape[-1] == 3:
+                    block_sizes = [(b,) * (ndim - 1) + (3,) for b in (16, 32, 64)]
+                else:
+                    block_sizes = [(b,) * ndim for b in (16, 32, 64)]
+                var_kwargs["block_size"] = block_sizes
+
+            if function_name == "profile_line":
+                fixed_kwargs["dst"] = (shape[0] - 32, shape[1] + 9)
+
+            if function_name == "median":
+                footprints = []
+                ndim = len(shape)
+                footprint_sizes = [3, 5, 7, 9] if ndim == 2 else [3, 5, 7]
+                for footprint_size in [3, 5, 7, 9]:
+                    footprints.append(
+                        np.ones((footprint_sizes,) * ndim, dtype=bool)
+                    )
+                var_kwargs["footprint"] = footprints
+
+            if function_name in ["gaussian", "unsharp_mask"]:
+                fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
+
+            B = FiltersBench(
                 function_name=function_name,
                 shape=shape,
                 dtypes=dtypes,
-                contiguous_labels=contiguous_labels,
-                index_str=index_str,
                 fixed_kwargs=fixed_kwargs,
                 var_kwargs=var_kwargs,
                 module_cpu=skimage.measure,
                 module_gpu=cucim.skimage.measure,
+                run_cpu=run_cpu,
             )
-            results = B.run_benchmark(duration=1)
+            results = B.run_benchmark(duration=args.duration)
             all_results = pd.concat([all_results, results["full"]])
 
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _moments.py
-    ("moments", dict(), dict(order=[1, 2, 3, 4]), False, False),
-    ("moments_central", dict(), dict(order=[1, 2, 3]), False, True),
-    # omited from benchmarks (only tiny arrays): moments_normalized, moments_hu
-    ("centroid", dict(), dict(), False, True),
-    ("inertia_tensor", dict(), dict(), False, True),
-    ("inertia_tensor_eigvals", dict(), dict(), False, True),
-    # _polygon.py
-    # TODO: approximate_polygon, subdivide_polygon
-    # block.py
-    (
-        "block_reduce",
-        dict(),
-        dict(
-            func=[
-                cp.sum,
-            ]
-        ),
-        True,
-        True,
-    ),  # variable block_size configured below
-    # entropy.py
-    ("shannon_entropy", dict(base=2), dict(), True, True),
-    # profile.py
-    (
-        "profile_line",
-        dict(src=(5, 7)),
-        dict(reduce_func=[cp.mean], linewidth=[1, 2, 4], order=[1, 3]),
-        True,
-        False,
-    ),  # variable block_size configured below
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
-
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
-            else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
-
-        if function_name == "gabor" and np.prod(shape) > 1000000:
-            # avoid cases that are too slow on the CPU
-            var_kwargs["frequency"] = [f for f in var_kwargs["frequency"] if f >= 0.1]
-
-        if function_name == "block_reduce":
-            ndim = len(shape)
-            if shape[-1] == 3:
-                block_sizes = [(b,) * (ndim - 1) + (3,) for b in (16, 32, 64)]
-            else:
-                block_sizes = [(b,) * ndim for b in (16, 32, 64)]
-            var_kwargs["block_size"] = block_sizes
-
-        if function_name == "profile_line":
-            fixed_kwargs["dst"] = (shape[0] - 32, shape[1] + 9)
-
-        if function_name == "median":
-            footprints = []
-            ndim = len(shape)
-            footprint_sizes = [3, 5, 7, 9] if ndim == 2 else [3, 5, 7]
-            for footprint_size in [3, 5, 7, 9]:
-                footprints.append(
-                    np.ones((footprint_sizes,) * ndim, dtype=bool)
-                )
-            var_kwargs["footprint"] = footprints
-
-        if function_name in ["gaussian", "unsharp_mask"]:
-            fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
-
-        B = FiltersBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.measure,
-            module_gpu=cucim.skimage.measure,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    with open(fbase + ".md", "wt") as f:
+        f.write(all_results.to_markdown())
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM measure functions')
+    func_name_choices = ['label', 'regionprops', 'moments', 'moments_central', 'centroid', 'inertia_tensor', 'inertia_tensor_eigvals', 'block_reduce', 'shannon_entropy', 'profile_line']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices=dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices=func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_metrics_bench.py b/benchmarks/skimage/cucim_metrics_bench.py
index b567298b4..af6b94814 100644
--- a/benchmarks/skimage/cucim_metrics_bench.py
+++ b/benchmarks/skimage/cucim_metrics_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import os
 import pickle
 
@@ -17,44 +18,51 @@ def set_args(self, dtype):
         imaged = cp.testing.shaped_arange(self.shape, dtype=dtype)
         imaged2 = cp.testing.shaped_arange(self.shape, dtype=dtype)
         imaged2 = imaged2 + 0.05 * cp.random.standard_normal(self.shape)
-        imaged /= imaged.max()
-        imaged2 /= imaged2.max()
+        imaged = imaged / imaged.max()
+        imaged2 = imaged2 / imaged2.max()
         imaged2 = imaged2.clip(0, 1.0)
         self.args_cpu = (cp.asnumpy(imaged), cp.asnumpy(imaged2))
         self.args_gpu = (imaged, imaged2)
 
 
-pfile = "cucim_metrics_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
-
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _structural_similarity.py
-    (
-        "structural_similarity",
-        dict(data_range=1.0),
-        dict(gradient=[False, True], gaussian_weights=[False, True]),
-        True,
-        True,
-    ),
-    # simple_metrics.py
-    ("mean_squared_error", dict(), dict(), True, True),
-    (
-        "normalized_root_mse",
-        dict(),
-        dict(normalization=["euclidean", "min-max", "mean"]),
-        True,
-        True,
-    ),
-    ("peak_signal_noise_ratio", dict(data_range=1.0), dict(), True, True),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+def main(args):
+
+    pfile = "cucim_metrics_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _structural_similarity.py
+        (
+            "structural_similarity",
+            dict(data_range=1.0),
+            dict(gradient=[False, True], gaussian_weights=[False, True]),
+            True,
+            True,
+        ),
+        # simple_metrics.py
+        ("mean_squared_error", dict(), dict(), True, True),
+        (
+            "normalized_root_mse",
+            dict(),
+            dict(normalization=["euclidean", "min-max", "mean"]),
+            True,
+            True,
+        ),
+        ("peak_signal_noise_ratio", dict(data_range=1.0), dict(), True, True),
+        ("normalized_mutual_information", dict(bins=100), dict(), True, True),
+
+    ]:
+        if function_name != args.func_name:
+            continue
 
         ndim = len(shape)
         if not allow_nd:
@@ -78,13 +86,28 @@ def set_args(self, dtype):
             var_kwargs=var_kwargs,
             module_cpu=skimage.metrics,
             module_gpu=cucim.skimage.metrics,
+            run_cpu=run_cpu,
         )
-        results = B.run_benchmark(duration=1)
+        results = B.run_benchmark(duration=args.duration)
         all_results = pd.concat([all_results, results["full"]])
 
 
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    with open(fbase + ".md", "wt") as f:
+        f.write(all_results.to_markdown())
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM metrics functions')
+    func_name_choices = ['structural_similarity', 'mean_squared_error', 'normalized_root_mse', 'peak_signal_noise_ratio', 'normalized_mutual_information']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices=dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices=func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_morphology_bench.py b/benchmarks/skimage/cucim_morphology_bench.py
index 4776a13b6..8fe1d6a95 100644
--- a/benchmarks/skimage/cucim_morphology_bench.py
+++ b/benchmarks/skimage/cucim_morphology_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import copy
 import functools
 import math
@@ -30,6 +31,7 @@ def __init__(
         var_kwargs={},
         module_cpu=skimage.morphology,
         module_gpu=cucim.skimage.morphology,
+        run_cpu=True,
     ):
 
         array_kwargs = dict(footprint=footprint)
@@ -47,6 +49,7 @@ def __init__(
             index_str=index_str,
             module_cpu=module_cpu,
             module_gpu=module_gpu,
+            run_cpu=run_cpu,
         )
 
     def set_args(self, dtype):
@@ -104,169 +107,169 @@ def set_args(self, dtype):
         self.args_gpu = (a, 5)
 
 
-pfile = "cucim_morphology_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes_gray = [np.float32]
-
-
-for function_name, fixed_kwargs, var_kwargs, allow_nd in [
-    ("thin", dict(), dict(), True),
-]:
-
-    for shape in [(512, 512), (3840, 2160)]:
-
-        ndim = len(shape)
-        if ndim != 2:
-            raise ValueError("only 2d benchmark data has been implemented")
-
-        if not allow_nd and ndim > 2:
+def main(args):
+
+    pfile = "cucim_morphology_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # binary.py
+        ("binary_erosion", dict(), dict(), False, True),
+        ("binary_dilation", dict(), dict(), False, True),
+        ("binary_opening", dict(), dict(), False, True),
+        ("binary_closing", dict(), dict(), False, True),
+        # misc.py
+        ("remove_small_objects", dict(), dict(), False, True),
+        ("remove_small_holes", dict(), dict(), False, True),
+        # gray.py
+        ("erosion", dict(), dict(), False, True),
+        ("dilation", dict(), dict(), False, True),
+        ("opening", dict(), dict(), False, True),
+        ("closing", dict(), dict(), False, True),
+        ("white_tophat", dict(), dict(), False, True),
+        ("black_tophat", dict(), dict(), False, True),
+        # _skeletonize.py
+        ("thin", dict(), dict(), False, True),
+        # grayreconstruct.py
+        ("reconstruction", dict(), dict(), False, True),
+        # footprints.py
+        # OMIT the functions from this file (each creates a structuring element)
+    ]:
+
+        if function_name != args.func_name:
             continue
 
-
-        B = SkeletonizeBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=[bool],
-            fixed_kwargs={},
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.morphology,
-            module_gpu=cucim.skimage.morphology,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-for function_name, fixed_kwargs, var_kwargs, allow_nd in [
-    ("binary_erosion", dict(), dict(), True),
-    ("binary_dilation", dict(), dict(), True),
-    ("binary_opening", dict(), dict(), True),
-    ("binary_closing", dict(), dict(), True),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (192, 192, 192)]:
-
         ndim = len(shape)
-        if not allow_nd and ndim > 2:
-            continue
+        if function_name == 'thin':
+            if ndim != 2:
+                raise ValueError("only 2d benchmark data has been implemented")
 
-        for connectivity in range(1, ndim + 1):
-            index_str = f"conn={connectivity}"
-            footprint = ndi.generate_binary_structure(ndim, connectivity)
+            if not allow_nd and ndim > 2:
+                continue
 
-            B = BinaryMorphologyBench(
+            B = SkeletonizeBench(
                 function_name=function_name,
                 shape=shape,
                 dtypes=[bool],
-                footprint=footprint,
                 fixed_kwargs={},
                 var_kwargs=var_kwargs,
-                index_str=index_str,
                 module_cpu=skimage.morphology,
                 module_gpu=cucim.skimage.morphology,
+                run_cpu=run_cpu,
             )
-            results = B.run_benchmark(duration=1)
-            all_results = pd.concat([all_results, results["full"]])
 
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+        elif function_name.startswith('binary_'):
+
+            if not allow_nd and ndim > 2:
+                continue
+
+            for connectivity in range(1, ndim + 1):
+                index_str = f"conn={connectivity}"
+                footprint = ndi.generate_binary_structure(ndim, connectivity)
+
+                B = BinaryMorphologyBench(
+                    function_name=function_name,
+                    shape=shape,
+                    dtypes=[bool],
+                    footprint=footprint,
+                    fixed_kwargs={},
+                    var_kwargs=var_kwargs,
+                    index_str=index_str,
+                    module_cpu=skimage.morphology,
+                    module_gpu=cucim.skimage.morphology,
+                    run_cpu=run_cpu,
+                )
+                results = B.run_benchmark(duration=args.duration)
+                all_results = pd.concat([all_results, results["full"]])
 
-for function_name, fixed_kwargs, var_kwargs, allow_nd in [
-    # misc.py
-    ("remove_small_objects", dict(), dict(), True),
-    ("remove_small_holes", dict(), dict(), True),
-]:
 
-    for shape in [(512, 512), (3840, 2160), (192, 192, 192)]:
+        elif function_name in ['remove_small_holes', 'remove_small_objects']:
+            if not allow_nd and ndim > 2:
+                continue
 
-        ndim = len(shape)
-        if not allow_nd and ndim > 2:
-            continue
-
-        if function_name == "remove_small_objects":
-            TestClass = RemoveSmallObjectsBench
-        elif function_name == "remove_small_holes":
-            TestClass = RemoveSmallHolesBench
-        else:
-            raise ValueError(f"unknown function: {function_name}")
-        B = TestClass(
-            function_name=function_name,
-            shape=shape,
-            dtypes=[bool],
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.morphology,
-            module_gpu=cucim.skimage.morphology,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # grey.py
-    ("erosion", dict(), dict(), False, True),
-    ("dilation", dict(), dict(), False, True),
-    ("opening", dict(), dict(), False, True),
-    ("closing", dict(), dict(), False, True),
-    ("white_tophat", dict(), dict(), False, True),
-    ("black_tophat", dict(), dict(), False, True),
-    # greyreconstruct.py
-    ("reconstruction", dict(), dict(), False, True),
-    # footprints.py
-    # OMIT the functions from this file (each creates a structuring element)
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+            if function_name == "remove_small_objects":
+                TestClass = RemoveSmallObjectsBench
+            elif function_name == "remove_small_holes":
+                TestClass = RemoveSmallHolesBench
 
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
             else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
+                raise ValueError(f"unknown function: {function_name}")
+            B = TestClass(
+                function_name=function_name,
+                shape=shape,
+                dtypes=[bool],
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.morphology,
+                module_gpu=cucim.skimage.morphology,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
 
-        if function_name == "gabor" and np.prod(shape) > 1000000:
-            # avoid cases that are too slow on the CPU
-            var_kwargs["frequency"] = [f for f in var_kwargs["frequency"] if f >= 0.1]
-
-        if function_name == "median":
-            footprints = []
-            ndim = len(shape)
-            footprint_sizes = [3, 5, 7, 9] if ndim == 2 else [3, 5, 7]
-            for footprint_size in [3, 5, 7, 9]:
-                footprints.append(
-                    np.ones((footprint_sizes,) * ndim, dtype=bool)
-                )
-            var_kwargs["footprint"] = footprints
+        else:
 
-        if function_name in ["gaussian", "unsharp_mask"]:
-            fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
+            if not allow_nd:
+                if not allow_color:
+                    if ndim > 2:
+                        continue
+                else:
+                    if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
+                        continue
 
-        if function_name == "reconstruction":
-            TestClass = ReconstructionBench
-        else:
-            TestClass = ImageBench
+            if shape[-1] == 3 and not allow_color:
+                continue
 
-        B = TestClass(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes_gray,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.morphology,
-            module_gpu=cucim.skimage.morphology,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
+            if function_name == "reconstruction":
+                TestClass = ReconstructionBench
+            else:
+                TestClass = ImageBench
 
+            B = TestClass(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.morphology,
+                module_gpu=cucim.skimage.morphology,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
 
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM morphology functions')
+    func_name_choices = ['binary_erosion', 'binary_dilation', 'binary_opening', 'binary_closing', 'remove_small_objects', 'remove_small_holes', 'erosion', 'dilation', 'opening', 'closing', 'white_tophat', 'black_tophat', 'thin', 'reconstruction']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_registration_bench.py b/benchmarks/skimage/cucim_registration_bench.py
index df1380d7f..6d275dcec 100644
--- a/benchmarks/skimage/cucim_registration_bench.py
+++ b/benchmarks/skimage/cucim_registration_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import math
 import os
 import pickle
@@ -33,48 +34,93 @@ def set_args(self, dtype):
         self.args_gpu = (imaged, imaged2)
 
 
-pfile = "cucim_registration_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
-
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _phase_cross_correlation.py
-    ("phase_cross_correlation", dict(), dict(), False, True),
-]:
+def main(args):
+
+    pfile = "cucim_registration_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _phase_cross_correlation.py
+        ("phase_cross_correlation", dict(), dict(), False, True),
+        # optical flow functions
+        ("optical_flow_tvl1", dict(), dict(num_iter=[10], num_warp=[5]), False, True),
+        (
+            "optical_flow_ilk",
+            dict(),
+            dict(radius=[3, 7], num_warp=[10], gaussian=[False, True], prefilter=[False, True]),
+            False,
+            True,
+        ),
+
+    ]:
+
+        if function_name != args.func_name:
+            continue
 
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+        if function_name == 'phase_cross_correlation':
+
+            ndim = len(shape)
+            if not allow_nd:
+                if not allow_color:
+                    if ndim > 2:
+                        continue
+                else:
+                    if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
+                        continue
+            if shape[-1] == 3 and not allow_color:
+                continue
+
+            for masked in [True, False]:
+
+                index_str = f"masked={masked}"
+                if masked:
+                    moving_mask = cp.ones(shape, dtype=bool)
+                    moving_mask[20:-20, :] = 0
+                    moving_mask[:, 20:-20] = 0
+                    reference_mask = cp.ones(shape, dtype=bool)
+                    reference_mask[80:-80, :] = 0
+                    reference_mask[:, 80:-80] = 0
+                    fixed_kwargs["moving_mask"] = moving_mask
+                    fixed_kwargs["reference_mask"] = reference_mask
+                else:
+                    fixed_kwargs["moving_mask"] = None
+                    fixed_kwargs["reference_mask"] = None
+
+                B = RegistrationBench(
+                    function_name=function_name,
+                    shape=shape,
+                    dtypes=dtypes,
+                    fixed_kwargs=fixed_kwargs,
+                    var_kwargs=var_kwargs,
+                    index_str=index_str,
+                    module_cpu=skimage.registration,
+                    module_gpu=cucim.skimage.registration,
+                    run_cpu=run_cpu,
+                )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
 
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
-            else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
+        else:
 
-        for masked in [True, False]:
-
-            index_str = f"masked={masked}"
-            if masked:
-                moving_mask = cp.ones(shape, dtype=bool)
-                moving_mask[20:-20, :] = 0
-                moving_mask[:, 20:-20] = 0
-                reference_mask = cp.ones(shape, dtype=bool)
-                reference_mask[80:-80, :] = 0
-                reference_mask[:, 80:-80] = 0
-                fixed_kwargs["moving_mask"] = moving_mask
-                fixed_kwargs["reference_mask"] = reference_mask
-            else:
-                fixed_kwargs["moving_mask"] = None
-                fixed_kwargs["reference_mask"] = None
+            ndim = len(shape)
+            if not allow_nd:
+                if not allow_color:
+                    if ndim > 2:
+                        continue
+                else:
+                    if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
+                        continue
+            if shape[-1] == 3 and not allow_color:
+                continue
 
             B = RegistrationBench(
                 function_name=function_name,
@@ -82,54 +128,35 @@ def set_args(self, dtype):
                 dtypes=dtypes,
                 fixed_kwargs=fixed_kwargs,
                 var_kwargs=var_kwargs,
-                index_str=index_str,
                 module_cpu=skimage.registration,
                 module_gpu=cucim.skimage.registration,
+                run_cpu=run_cpu,
             )
-            results = B.run_benchmark(duration=1)
+            results = B.run_benchmark(duration=args.duration)
             all_results = pd.concat([all_results, results["full"]])
 
 
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _phase_cross_correlation.py
-    ("optical_flow_tvl1", dict(), dict(num_iter=[10], num_warp=[5]), False, True),
-    (
-        "optical_flow_ilk",
-        dict(),
-        dict(radius=[3, 7], num_warp=[10], gaussian=[False, True], prefilter=[False, True]),
-        False,
-        True,
-    ),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
-
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
-            else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM registration functions')
+    func_name_choices = ['phase_cross_correlation', 'optical_flow_tvl1', 'optical_flow_ilk']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
 
-        B = RegistrationBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.registration,
-            module_gpu=cucim.skimage.registration,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_restoration_bench.py b/benchmarks/skimage/cucim_restoration_bench.py
index 0de362d5a..a1523dbae 100644
--- a/benchmarks/skimage/cucim_restoration_bench.py
+++ b/benchmarks/skimage/cucim_restoration_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import math
 import os
 import pickle
@@ -97,65 +98,35 @@ def set_args(self, dtype):
         self.args_gpu = (imaged, psfd)
 
 
-pfile = "cucim_restoration_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
 
+def main(args):
 
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _denoise.py
-    ("denoise_tv_chambolle", dict(), dict(weight=[0.02]), True, True),
-    # j_invariant.py
-    ("calibrate_denoiser", dict(), dict(), False, True),
-]:
+    pfile = "cucim_restoration_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
 
-    for shape in [(512, 512), (1980, 1080), (1980, 1080, 3), (128, 128, 128)]:
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
 
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
-            else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
 
-        if function_name == "denoise_tv_chambolle":
-            fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _denoise.py
+        ("denoise_tv_chambolle", dict(), dict(weight=[0.02]), True, True),
+        # j_invariant.py
+        ("calibrate_denoiser", dict(), dict(), False, True),
+        # deconvolution.py
+        ("wiener", dict(balance=100.0), dict(), False, False),
+        ("unsupervised_wiener", dict(), dict(), False, False),
+        ("richardson_lucy", dict(), dict(num_iter=[5]), False, True),
+    ]:
 
-        if function_name == "calibrate_denoiser":
-            denoise_class = CalibratedDenoiseBench
-        else:
-            denoise_class = DenoiseBench
-
-        B = denoise_class(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.restoration,
-            module_gpu=cucim.skimage.restoration,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-# function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd = ('unsupervised_wiener', dict(), dict(), False, True)
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # deconvolution.py
-    ("wiener", dict(balance=100.0), dict(), False, False),
-    ("unsupervised_wiener", dict(), dict(), False, False),
-    ("richardson_lucy", dict(), dict(num_iter=[5]), False, True),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+        if function_name != args.func_name:
+            continue
 
         ndim = len(shape)
         if not allow_nd:
@@ -168,21 +139,65 @@ def set_args(self, dtype):
         if shape[-1] == 3 and not allow_color:
             continue
 
-        B = DeconvolutionBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.restoration,
-            module_gpu=cucim.skimage.restoration,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+        if function_name in ['denoise_tv_chambolle', 'calibrate_denoiser']:
+
+            if function_name == "denoise_tv_chambolle":
+                fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
+
+            if function_name == "calibrate_denoiser":
+                denoise_class = CalibratedDenoiseBench
+            else:
+                denoise_class = DenoiseBench
+
+            B = denoise_class(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.restoration,
+                module_gpu=cucim.skimage.restoration,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+        elif function_name in ['wiener', 'unsupervised_wiener', 'richardson_lucy']:
+
+            B = DeconvolutionBench(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.restoration,
+                module_gpu=cucim.skimage.restoration,
+                run_cpu=run_cpu,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM restoration functions')
+    func_name_choices = ['denoise_tv_chambolle', 'calibrate_denoiser', 'wiener', 'unsupervised_wiener', 'richardson_lucy']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
+
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_segmentation_bench.py b/benchmarks/skimage/cucim_segmentation_bench.py
index 47d7d0834..382349a9b 100644
--- a/benchmarks/skimage/cucim_segmentation_bench.py
+++ b/benchmarks/skimage/cucim_segmentation_bench.py
@@ -1,8 +1,11 @@
+import argparse
 import math
 import os
 import pickle
 
 import cucim.skimage
+import cucim.skimage.data
+import cucim.skimage.exposure
 import cucim.skimage.segmentation
 import cupy as cp
 import numpy as np
@@ -25,6 +28,7 @@ def __init__(
         index_str=None,
         module_cpu=skimage.measure,
         module_gpu=cucim.skimage.measure,
+        run_cpu=True,
     ):
 
         self.contiguous_labels = contiguous_labels
@@ -38,6 +42,7 @@ def __init__(
             index_str=index_str,
             module_cpu=module_cpu,
             module_gpu=module_gpu,
+            run_cpu=run_cpu,
         )
 
     def set_args(self, dtype):
@@ -50,14 +55,38 @@ def set_args(self, dtype):
             ],
             dtype=dtype,
         )
-        tiling = tuple(s // a_s for s, a_s in zip(shape, a.shape))
+        tiling = tuple(s // a_s for s, a_s in zip(self.shape, a.shape))
         if self.contiguous_labels:
-            image = np.kron(a, np.ones(tiling, dtype=a.dtype))
+            labels = np.kron(a, np.ones(tiling, dtype=a.dtype))
         else:
-            image = np.tile(a, tiling)
-        imaged = cp.asarray(image)
-        self.args_cpu = (image,)
-        self.args_gpu = (imaged,)
+            labels = np.tile(a, tiling)
+        labels_d = cp.asarray(labels)
+        self.args_cpu = (labels,)
+        self.args_gpu = (labels_d,)
+
+
+class LabelAndImageBench(LabelBench):
+
+    def set_args(self, dtype):
+        a = np.array(
+            [
+                [0, 0, 1, 1, 0, 0, 0, 0],
+                [0, 0, 0, 1, 0, 0, 4, 0],
+                [2, 2, 0, 0, 3, 0, 4, 4],
+                [0, 0, 0, 0, 0, 5, 0, 0],
+            ],
+            dtype=dtype,
+        )
+        tiling = tuple(s // a_s for s, a_s in zip(self.shape, a.shape))
+        if self.contiguous_labels:
+            labels = np.kron(a, np.ones(tiling, dtype=a.dtype))
+        else:
+            labels = np.tile(a, tiling)
+        labels_d = cp.asarray(labels)
+        image_d = cp.random.standard_normal(labels.shape).astype(np.float32)
+        image = cp.asnumpy(image_d)
+        self.args_cpu = (image, labels)
+        self.args_gpu = (image_d, labels_d)
 
 
 class MorphGeodesicBench(ImageBench):
@@ -75,7 +104,9 @@ def set_args(self, dtype):
         imaged = cp.tile(im1, n_tile)[slices]
 
         # need this preprocessing for morphological_geodesic_active_contour
-        imaged = skimage.segmentation.inverse_gaussian_gradient(imaged)
+        imaged = cp.array(
+            skimage.segmentation.inverse_gaussian_gradient(cp.asnumpy(imaged))
+        )
 
         image = cp.asnumpy(imaged)
         assert imaged.dtype == dtype
@@ -84,29 +115,118 @@ def set_args(self, dtype):
         self.args_gpu = (imaged,)
 
 
-pfile = "cucim_segmentation_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.int32]
+class RandomWalkerBench(ImageBench):
+
+
+    def set_args(self, dtype):
+        # Note: dtype only used for merkers array, data is hard-coded as float32
+
+        if np.dtype(dtype).kind not in 'iu':
+            raise ValueError("random_walker markers require integer dtype")
+
+        n_dim = len(self.shape)
+        data = cucim.skimage.img_as_float(
+            cucim.skimage.data.binary_blobs(
+                length=max(self.shape), n_dim=n_dim, seed=1
+            )
+        )
+        data = data[tuple(slice(s) for s in self.shape)]
+        sigma = 0.35
+        rng = np.random.default_rng(5)
+        data += cp.array(rng.normal(loc=0, scale=sigma, size=data.shape))
+        data = cucim.skimage.exposure.rescale_intensity(
+            data, in_range=(-sigma, 1 + sigma), out_range=(-1, 1)
+        )
+        data = data.astype(cp.float32)
+        data_cpu = cp.asnumpy(data)
+
+        # The range of the binary image spans over (-1, 1).
+        # We choose the hottest and the coldest pixels as markers.
+        markers = cp.zeros(data.shape, dtype=dtype)
+        markers[data < -0.95] = 1
+        markers[data > 0.95] = 2
+        markers_cpu = cp.asnumpy(markers)
+        self.args_cpu = (data_cpu, markers_cpu)
+        self.args_gpu = (data, markers)
+
+
+
+def main(args):
+
+    pfile = "cucim_segmentation_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
 
+    dtypes = [np.dtype(args.dtype)]
+    dtypes_label = [np.dtype(args.dtype_label)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
 
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _denoise.py
-    (
-        "find_boundaries",
-        dict(),
-        dict(connectivity=[1], mode=["thick", "inner", "outer", "subpixel"]),
-        False,
-        True,
-    ),
-]:
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _clear_border.py
+        (
+            "clear_border",
+            dict(),
+            dict(),
+            False,
+            True,
+        ),
+        # _join.py
+        (
+            "relabel_sequential",
+            dict(offset=5),
+            dict(),
+            False,
+            True,
+        ),
+        # boundaries.py
+        (
+            "find_boundaries",
+            dict(),
+            dict(connectivity=[1], mode=["thick", "inner", "outer", "subpixel"]),
+            False,
+            True,
+        ),
+        (
+            "mark_boundaries",
+            dict(),
+            dict(),
+            False,
+            True,
+        ),
+        (
+            "random_walker",
+            dict(beta=4, tol=1.e-4, prob_tol=1.e-2),
+            dict(mode=['cg', 'cg_j']),
+            False,
+            True,
+        ),
+        # morphsnakes.py
+        ("inverse_gaussian_gradient", dict(), dict(), False, True),
+        (
+            "morphological_geodesic_active_contour",
+            dict(),
+            dict(num_iter=[16], init_level_set=["checkerboard", "disk"]),
+            False,
+            False,
+        ),
+        (
+            "morphological_chan_vese",
+            dict(),
+            dict(num_iter=[16], init_level_set=["checkerboard", "disk"]),
+            False,
+            False,
+        ),
+        # omit: disk_level_set (simple array generation function)
+        # omit: checkerboard_level_set (simple array generation function)
+    ]:
 
-    for shape in [
-        (64, 64),
-    ]:  # (512, 512), (1980, 1080), (1980, 1080, 3), (128, 128, 128)]:
+        if function_name != args.func_name:
+            continue
 
         ndim = len(shape)
         if not allow_nd:
@@ -119,73 +239,72 @@ def set_args(self, dtype):
         if shape[-1] == 3 and not allow_color:
             continue
 
-        B = LabelBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.segmentation,
-            module_gpu=cucim.skimage.segmentation,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
-
-
-dtypes = [np.float32]
-# function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd = ('unsupervised_wiener', dict(), dict(), False, True)
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # morphsnakes.py
-    ("inverse_gaussian_gradient", dict(), dict(), False, True),
-    (
-        "morphological_geodesic_active_contour",
-        dict(),
-        dict(num_iter=[16], init_level_set=["checkerboard", "disk"]),
-        False,
-        False,
-    ),
-    (
-        "morphological_chan_vese",
-        dict(),
-        dict(num_iter=[16], init_level_set=["checkerboard", "disk"]),
-        False,
-        False,
-    ),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+        if function_name in ["clear_border", "relabel_sequential", "find_boundaries", "mark_boundaries", "random_walker"]:
+            if function_name == 'random_walker':
+                fixed_kwargs['channel_axis'] = -1 if shape[-1] == 3 else None
 
-        ndim = len(shape)
-        if not allow_nd:
-            if not allow_color:
-                if ndim > 2:
-                    continue
+            if function_name == 'mark_boundaries':
+                bench_func = LabelAndImageBench
+            elif function_name == 'random_walker':
+                bench_func = RandomWalkerBench
             else:
-                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                    continue
-        if shape[-1] == 3 and not allow_color:
-            continue
+                bench_func = LabelBench
 
-        if function_name == "morphological_geodesic_active_contour":
-            bench_class = MorphGeodesicBench
-        else:
-            bench_class = ImageBench
+            B = bench_func(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes_label,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.segmentation,
+                module_gpu=cucim.skimage.segmentation,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+
+        elif function_name in ["inverse_gaussian_gradient", "morphological_geodesic_active_contour", "morphological_chan_vese"]:
+
+            if function_name == "morphological_geodesic_active_contour":
+                bench_class = MorphGeodesicBench
+            else:
+                bench_class = ImageBench
+
+            B = bench_class(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.segmentation,
+                module_gpu=cucim.skimage.segmentation,
+            )
+            results = B.run_benchmark(duration=args.duration)
+            all_results = pd.concat([all_results, results["full"]])
+
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
 
-        B = ImageBench(
-            function_name=function_name,
-            shape=shape,
-            dtypes=dtypes,
-            fixed_kwargs=fixed_kwargs,
-            var_kwargs=var_kwargs,
-            module_cpu=skimage.segmentation,
-            module_gpu=cucim.skimage.segmentation,
-        )
-        results = B.run_benchmark(duration=1)
-        all_results = pd.concat([all_results, results["full"]])
 
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM segmentation functions')
+    func_name_choices = ["clear_border", "relabel_sequential", "find_boundaries", "mark_boundaries", "random_walker", "inverse_gaussian_gradient", "morphological_geodesic_active_contour", "morphological_chan_vese"]
+    label_dtype_choices = ['int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('--dtype_label', type=str, help='Dtype of input image', choices = label_dtype_choices, required=False, default='uint8')
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
 
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/cucim_transform_bench.py b/benchmarks/skimage/cucim_transform_bench.py
index 91e2fb27f..b6aa315f9 100644
--- a/benchmarks/skimage/cucim_transform_bench.py
+++ b/benchmarks/skimage/cucim_transform_bench.py
@@ -1,3 +1,4 @@
+import argparse
 import os
 import pickle
 
@@ -10,75 +11,83 @@
 
 from _image_bench import ImageBench
 
-pfile = "cucim_transform_results.pickle"
-if os.path.exists(pfile):
-    with open(pfile, "rb") as f:
-        all_results = pickle.load(f)
-else:
-    all_results = pd.DataFrame()
-dtypes = [np.float32]
-
-for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
-    # _warps.py
-    (
-        "resize",
-        dict(preserve_range=True),
-        dict(order=[0, 1, 3], mode=["reflect"], anti_aliasing=[True]),
-        True,
-        True,
-    ),  # scale handled in loop below
-    (
-        "resize_local_mean",
-        dict(preserve_range=True),
-        {},
-        True,
-        True,
-    ),  # scale handled in loop below
-    (
-        "rescale",
-        dict(preserve_range=True),
-        dict(order=[0, 1, 3], mode=["reflect"], anti_aliasing=[True]),
-        True,
-        True,
-    ),  # output_shape handled in loop below
-    (
-        "rotate",
-        dict(angle=15, preserve_range=True),
-        dict(order=[0, 1, 3], mode=["reflect"], resize=[False, True]),
-        False,
-        False,
-    ),
-    ("downscale_local_mean", dict(), dict(), True, True),  # factors handled in loop below
-    (
-        "swirl",
-        dict(strength=1, preserve_range=True),
-        dict(order=[0, 1, 3], mode=["reflect"]),
-        False,
-        False,
-    ),
-    # TODO : warp? already indirectly benchmarked via swirl, etc
-    ("warp_polar", dict(), dict(scaling=["linear", "log"]), True, False),
-    # integral.py
-    ("integral_image", dict(), dict(), False, True),
-    # TODO: integrate
-    # pyramids.py
-    (
-        "pyramid_gaussian",
-        dict(max_layer=6, downscale=2, preserve_range=True),
-        dict(order=[0, 1, 3]),
-        True,
-        True,
-    ),
-    (
-        "pyramid_laplacian",
-        dict(max_layer=6, downscale=2, preserve_range=True),
-        dict(order=[0, 1, 3]),
-        True,
-        True,
-    ),
-]:
-
-    for shape in [(512, 512), (3840, 2160), (3840, 2160, 3), (192, 192, 192)]:
+
+def main(args):
+
+    pfile = "cucim_transform_results.pickle"
+    if os.path.exists(pfile):
+        with open(pfile, "rb") as f:
+            all_results = pickle.load(f)
+    else:
+        all_results = pd.DataFrame()
+
+    dtypes = [np.dtype(args.dtype)]
+    # image sizes/shapes
+    shape = tuple(list(map(int,(args.img_size.split(',')))))
+    run_cpu = not args.no_cpu
+
+    for function_name, fixed_kwargs, var_kwargs, allow_color, allow_nd in [
+        # _warps.py
+        (
+            "resize",
+            dict(preserve_range=True),
+            dict(order=[0, 1, 3], mode=["reflect"], anti_aliasing=[True]),
+            True,
+            True,
+        ),  # scale handled in loop below
+        (
+            "resize_local_mean",
+            dict(preserve_range=True),
+            {},
+            True,
+            True,
+        ),  # scale handled in loop below
+        (
+            "rescale",
+            dict(preserve_range=True),
+            dict(order=[0, 1, 3], mode=["reflect"], anti_aliasing=[True]),
+            True,
+            True,
+        ),  # output_shape handled in loop below
+        (
+            "rotate",
+            dict(angle=15, preserve_range=True),
+            dict(order=[0, 1, 3], mode=["reflect"], resize=[False, True]),
+            False,
+            False,
+        ),
+        ("downscale_local_mean", dict(), dict(), True, True),  # factors handled in loop below
+        (
+            "swirl",
+            dict(strength=1, preserve_range=True),
+            dict(order=[0, 1, 3], mode=["reflect"]),
+            False,
+            False,
+        ),
+        # TODO : warp? already indirectly benchmarked via swirl, etc
+        ("warp_polar", dict(), dict(scaling=["linear", "log"]), True, False),
+        # integral.py
+        ("integral_image", dict(), dict(), False, True),
+        # TODO: integrate
+        # pyramids.py
+        (
+            "pyramid_gaussian",
+            dict(max_layer=6, downscale=2, preserve_range=True),
+            dict(order=[0, 1, 3]),
+            True,
+            True,
+        ),
+        (
+            "pyramid_laplacian",
+            dict(max_layer=6, downscale=2, preserve_range=True),
+            dict(order=[0, 1, 3]),
+            True,
+            True,
+        ),
+    ]:
+
+        if function_name != args.func_name:
+            continue
 
         ndim = len(shape)
         if not allow_nd:
@@ -131,12 +140,30 @@
             module_gpu=cucim.skimage.transform,
             function_is_generator=function_is_generator,
         )
-        results = B.run_benchmark(duration=1)
+        results = B.run_benchmark(duration=args.duration)
         all_results = pd.concat([all_results, results["full"]])
 
+    fbase = os.path.splitext(pfile)[0]
+    all_results.to_csv(fbase + ".csv")
+    all_results.to_pickle(pfile)
+    try:
+        import tabular
+
+        with open(fbase + ".md", "wt") as f:
+            f.write(all_results.to_markdown())
+    except ImportError:
+        pass
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Benchmarking cuCIM transform functions')
+    func_name_choices = ['resize', 'resize_local_mean', 'rescale', 'rotate', 'downscale_local_mean', 'warp_polar', 'integral_image', 'pyramid_gaussian', 'pyramid_laplacian']
+    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
+    parser.add_argument('-i','--img_size', type=str, help='Size of input image (omit color channel, it will be appended as needed)', required=True)
+    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices = dtype_choices, required=True)
+    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices = func_name_choices, required=True)
+    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
+    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
 
-fbase = os.path.splitext(pfile)[0]
-all_results.to_csv(fbase + ".csv")
-all_results.to_pickle(pfile)
-with open(fbase + ".md", "wt") as f:
-    f.write(all_results.to_markdown())
+    args = parser.parse_args()
+    main(args)
diff --git a/benchmarks/skimage/requirements-bench.txt b/benchmarks/skimage/requirements-bench.txt
new file mode 100644
index 000000000..2cc6ac30e
--- /dev/null
+++ b/benchmarks/skimage/requirements-bench.txt
@@ -0,0 +1,2 @@
+pandas>=1.0
+tabulate>=0.8.7
diff --git a/benchmarks/skimage/run-nv-bench-color.sh b/benchmarks/skimage/run-nv-bench-color.sh
new file mode 100755
index 000000000..14a0809ec
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-color.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 192,192,192)
+param_filt=(convert_colorspace rgb2hed hed2rgb lab2lch lch2lab xyz2lab lab2xyz rgba2rgb label2rgb)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_color_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-exposure.sh b/benchmarks/skimage/run-nv-bench-exposure.sh
new file mode 100755
index 000000000..4ed9e0748
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-exposure.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(equalize_adapthist cumulative_distribution equalize_hist rescale_intensity adjust_gamma adjust_log adjust_sigmoid is_low_contrast match_histograms)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_exposure_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-feature.sh b/benchmarks/skimage/run-nv-bench-feature.sh
index c777bf40a..fbab913c6 100755
--- a/benchmarks/skimage/run-nv-bench-feature.sh
+++ b/benchmarks/skimage/run-nv-bench-feature.sh
@@ -1,7 +1,7 @@
 #!/bin/bash
 param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
 param_filt=(multiscale_basic_features canny daisy structure_tensor hessian_matrix hessian_matrix_det shape_index corner_kitchen_rosenfeld corner_harris corner_shi_tomasi corner_foerstner corner_peaks match_template)
-param_dt=('float64' 'float32' 'float16')
+param_dt=(float64 float32 float16)
 for shape in "${param_shape[@]}"; do
     for filt in "${param_filt[@]}"; do
         for dt in "${param_dt[@]}"; do
diff --git a/benchmarks/skimage/run-nv-bench-filters.sh b/benchmarks/skimage/run-nv-bench-filters.sh
index 52f94e0a6..37388e9e9 100755
--- a/benchmarks/skimage/run-nv-bench-filters.sh
+++ b/benchmarks/skimage/run-nv-bench-filters.sh
@@ -2,7 +2,7 @@
 param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
 param_filt=(gabor gaussian median rank_order unsharp_mask sobel prewitt scharr roberts roberts_pos_diag roberts_neg_diag farid laplace meijering sato frangi hessian threshold_isodata threshold_otsu threshold_yen threshold_local threshold_li threshold_minimum threshold_mean threshold_triangle threshold_niblack threshold_sauvola apply_hysteresis_threshold threshold_multiotsu)
 # param_filt=(rank_order )
-param_dt=('float64' 'float32' 'float16')
+param_dt=(float64 float32 float16)
 for shape in "${param_shape[@]}"; do
     for filt in "${param_filt[@]}"; do
         for dt in "${param_dt[@]}"; do
diff --git a/benchmarks/skimage/run-nv-bench-measure.sh b/benchmarks/skimage/run-nv-bench-measure.sh
new file mode 100755
index 000000000..5c2b4f413
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-measure.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(label regionprops moments moments_central centroid inertia_tensor inertia_tensor_eigvals block_reduce shannon_entropy profile_line)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_measure_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-metrics.sh b/benchmarks/skimage/run-nv-bench-metrics.sh
new file mode 100755
index 000000000..b24ba775b
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-metrics.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(structural_similarity mean_squared_error normalized_root_mse peak_signal_noise_ratio normalized_mutual_information)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_metrics_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-morphology.sh b/benchmarks/skimage/run-nv-bench-morphology.sh
new file mode 100755
index 000000000..4981a58ad
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-morphology.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(binary_erosion binary_dilation binary_opening binary_closing remove_small_objects remove_small_holes erosion dilation opening closing white_tophat black_tophat thin reconstruction)
+# Note: user-specified dtype ignored for binary_* functions and thin (these only accept binary input)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_morphology_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-registration.sh b/benchmarks/skimage/run-nv-bench-registration.sh
new file mode 100755
index 000000000..18f572b53
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-registration.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(phase_cross_correlation optical_flow_tvl1 optical_flow_ilk)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_registration_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-restoration.sh b/benchmarks/skimage/run-nv-bench-restoration.sh
new file mode 100755
index 000000000..df7c9141f
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-restoration.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(denoise_tv_chambolle calibrate_denoiser wiener unsupervised_wiener richardson_lucy)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_restoration_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-segmentation.sh b/benchmarks/skimage/run-nv-bench-segmentation.sh
new file mode 100755
index 000000000..4723f0c60
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-segmentation.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+
+# these require an integer-valued label image
+param_filt=(clear_border relabel_sequential find_boundaries mark_boundaries random_walker)
+param_dt=(float32)
+param_dt_label=(uint8 uint32)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            for dt_label in "${param_dt_label[@]}"; do
+                python cucim_segmentation_bench.py -f $filt -i $shape -d $dt --dtype_label $dt_label -t 10
+            done
+        done
+    done
+done
+
+# these do not require an integer-valued input image
+param_filt=(inverse_gaussian_gradient morphological_geodesic_active_contour morphological_chan_vese)
+param_dt=(float32)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_segmentation_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/benchmarks/skimage/run-nv-bench-transform.sh b/benchmarks/skimage/run-nv-bench-transform.sh
new file mode 100755
index 000000000..c880cb699
--- /dev/null
+++ b/benchmarks/skimage/run-nv-bench-transform.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
+param_filt=(resize resize_local_mean rescale rotate downscale_local_mean warp_polar integral_image pyramid_gaussian pyramid_laplacian)
+param_dt=(float32 uint8)
+for shape in "${param_shape[@]}"; do
+    for filt in "${param_filt[@]}"; do
+        for dt in "${param_dt[@]}"; do
+            python cucim_transform_bench.py -f $filt -i $shape -d $dt -t 10
+        done
+    done
+done
diff --git a/conda/recipes/cucim/meta.yaml b/conda/recipes/cucim/meta.yaml
index fa9377b93..a38d1bdbc 100644
--- a/conda/recipes/cucim/meta.yaml
+++ b/conda/recipes/cucim/meta.yaml
@@ -29,7 +29,7 @@ requirements:
     - python {{ python_version }}.*
     - libcucim {{ version }}.*
     - click
-    - cupy >=9,<11.0.0a0
+    - cupy >=9,<12.0.0a0
     - numpy 1.19
     - scipy
     - scikit-image >=0.18.1,<0.20.0a0
@@ -38,7 +38,7 @@ requirements:
     - python {{ python_version }}.*
     - libcucim {{ version }}.*
     - click
-    - cupy >=9,<11.0.0a0
+    - cupy >=9,<12.0.0a0
     - {{ pin_compatible('numpy') }}
     - scipy
     - scikit-image >=0.18.1,<0.20.0a0
diff --git a/docs/source/_static/params.css b/docs/source/_static/params.css
deleted file mode 100644
index c080d3669..000000000
--- a/docs/source/_static/params.css
+++ /dev/null
@@ -1,30 +0,0 @@
-/* Mirrors the change in:
- *   https://github.com/sphinx-doc/sphinx/pull/5976
- * which is not showing up in our theme.
- */
-.classifier:before {
-    font-style: normal;
-    margin: 0.5em;
-    content: ":";
-}
-
-:root {
-
-  --pst-color-active-navigation: 114, 83, 237;
-  --pst-color-navbar-link: 77, 77, 77;
-  --pst-color-navbar-link-hover: var(--pst-color-active-navigation);
-  --pst-color-navbar-link-active: var(--pst-color-active-navigation);
-  --pst-color-sidebar-link: 77, 77, 77;
-  --pst-color-sidebar-link-hover: var(--pst-color-active-navigation);
-  --pst-color-sidebar-link-active: var(--pst-color-active-navigation);
-  --pst-color-sidebar-expander-background-hover: 244, 244, 244;
-  --pst-color-sidebar-caption: 77, 77, 77;
-  --pst-color-toc-link: 119, 117, 122;
-  --pst-color-toc-link-hover: var(--pst-color-active-navigation);
-  --pst-color-toc-link-active: var(--pst-color-active-navigation);
-
-}
-
-.special-table td, .special-table th {
-    border: 1px solid #dee2e6;
-}
diff --git a/docs/source/conf.py b/docs/source/conf.py
index e51d398cd..3efe01e03 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -182,5 +182,5 @@
 
 
 def setup(app):
-    app.add_css_file('params.css')
     app.add_css_file("https://docs.rapids.ai/assets/css/custom.css")
+    app.add_js_file("https://docs.rapids.ai/assets/js/custom.js")
diff --git a/python/cucim/CHANGELOG.md b/python/cucim/CHANGELOG.md
index d29d31b41..81d223e82 100644
--- a/python/cucim/CHANGELOG.md
+++ b/python/cucim/CHANGELOG.md
@@ -1,6 +1,30 @@
 
 # Changelog (See [Release Notes](https://github.com/rapidsai/cucim/wiki/Release-Notes))
 
+## [22.02.06](https://github.com/rapidsai/cucim/wiki/release_notes_v22.06.00)
+
+- [Update/Breaking] Promote small integer types to single rather than double precision ([#278](https://github.com/rapidsai/cucim/pull/278)) [@grlee77](https://github.com/grlee77)
+- [Bug] Populate correct channel names for RGBA image ([#294](https://github.com/rapidsai/cucim/pull/294)) [@gigony](https://github.com/gigony)
+- [Bug] Merge branch-22.04 into branch-22.06 ([#258](https://github.com/rapidsai/cucim/pull/258)) [@jakirkham](https://github.com/jakirkham)
+- [New] add missing `cucim.skimage.segmentation.clear_border` function ([#267](https://github.com/rapidsai/cucim/pull/267)) [@grlee77](https://github.com/grlee77)
+- [New] add `cucim.core.operations.color.stain_extraction_pca` and `cucim.core.operations.color.normalize_colors_pca` for digital pathology H&E stain extraction and normalization ([#273](https://github.com/rapidsai/cucim/pull/273)) [@grlee77](https://github.com/grlee77), [@drbeh](https://github.com/drbeh)
+- [Update] Update to use DLPack v0.6 ([#295](https://github.com/rapidsai/cucim/pull/295)) [@gigony](https://github.com/gigony)
+- [Update] Remove plugin-related messages temporarily ([#291](https://github.com/rapidsai/cucim/pull/291)) [@gigony](https://github.com/gigony)
+- [Update] Simplify recipes ([#286](https://github.com/rapidsai/cucim/pull/286)) [@Ethyling](https://github.com/Ethyling)
+- [Update] Use cupy.fuse to improve efficiency hessian_matrix_eigvals ([#280](https://github.com/rapidsai/cucim/pull/280)) [@grlee77](https://github.com/grlee77)
+- [Update] improve efficiency of histogram-based thresholding functions ([#276](https://github.com/rapidsai/cucim/pull/276)) [@grlee77](https://github.com/grlee77)
+- [Update] Remove unused dependencies in GPU tests job ([#268](https://github.com/rapidsai/cucim/pull/268)) [@Ethyling](https://github.com/Ethyling)
+- [Update] Enable footprint decomposition for morphology ([#274](https://github.com/rapidsai/cucim/pull/274)) [@grlee77](https://github.com/grlee77)
+- [Update] Use conda compilers ([#232](https://github.com/rapidsai/cucim/pull/232)) [@Ethyling](https://github.com/Ethyling)
+- [Update] Build packages using mambabuild ([#216](https://github.com/rapidsai/cucim/pull/216)) [@Ethyling](https://github.com/Ethyling)
+- [Doc] update outdated links to example data ([#289](https://github.com/rapidsai/cucim/pull/289)) [@grlee77](https://github.com/grlee77)
+- [Doc] Add missing API docs ([#275](https://github.com/rapidsai/cucim/pull/275)) [@grlee77](https://github.com/grlee77)
+
+## [22.02.04](https://github.com/rapidsai/cucim/wiki/release_notes_v22.04.00)
+
+- [Bug] Fix ImportError from vendored code ([#252](https://github.com/rapidsai/cucim/pull/252)) [@grlee77](https://github.com/grlee77)
+- [Bug] Fix wrong dimension in metadata ([#248](https://github.com/rapidsai/cucim/pull/248)) [@gigony](https://github.com/gigony)
+
 ## [22.02.01](https://github.com/rapidsai/cucim/wiki/release_notes_v22.02.01)
 
 - [Bug] Check nullptr of handler in CuFileDriver::close() ([#229](https://github.com/rapidsai/cucim/pull/229)) [@gigony](https://github.com/gigony)
diff --git a/python/cucim/README.md b/python/cucim/README.md
index 300e4c5bf..8f1413137 100644
--- a/python/cucim/README.md
+++ b/python/cucim/README.md
@@ -6,10 +6,13 @@
 
 **NOTE:** For the latest stable [README.md](https://github.com/rapidsai/cucim/blob/main/README.md) ensure you are on the `main` branch.
 
-- [GTC 2021 cuCIM: A GPU Image I/O and Processing Toolkit [S32194]](https://www.nvidia.com/en-us/on-demand/search/?facet.mimetype[]=event%20session&layout=list&page=1&q=cucim&sort=date)
-  - [video](https://www.nvidia.com/en-us/on-demand/session/gtcspring21-s32194/)
+- [GTC 2022 Accelerating Storage IO to GPUs with Magnum IO [S41347]](https://events.rainfocus.com/widget/nvidia/gtcspring2022/sessioncatalog/session/1634960000577001Etxp)
+  - cuCIM's GDS API examples: <https://github.com/NVIDIA/MagnumIO/tree/main/gds/readers/cucim-gds>
 - [SciPy 2021 cuCIM - A GPU image I/O and processing library](https://www.scipy2021.scipy.org/)
   - [video](https://youtu.be/G46kOOM9xbQ)
+- [GTC 2021 cuCIM: A GPU Image I/O and Processing Toolkit [S32194]](https://www.nvidia.com/en-us/on-demand/search/?facet.mimetype[]=event%20session&layout=list&page=1&q=cucim&sort=date)
+  - [video](https://www.nvidia.com/en-us/on-demand/session/gtcspring21-s32194/)
+
 
 ## Quick Start
 
@@ -178,4 +181,4 @@ is used in this project.
 
 Apache-2.0 License (see `LICENSE` file).
 
-Copyright (c) 2020-2021, NVIDIA CORPORATION.
+Copyright (c) 2020-2022, NVIDIA CORPORATION.
diff --git a/python/cucim/src/cucim/core/__init__.py b/python/cucim/src/cucim/core/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/python/cucim/src/cucim/core/operations/color/stain_normalizer.py b/python/cucim/src/cucim/core/operations/color/stain_normalizer.py
index f5141ede2..35dcb6e0e 100644
--- a/python/cucim/src/cucim/core/operations/color/stain_normalizer.py
+++ b/python/cucim/src/cucim/core/operations/color/stain_normalizer.py
@@ -314,7 +314,7 @@ def stain_extraction_pca(image, source_intensity=240, alpha=1, beta=0.345,
     )
 
     # remove transparent pixels
-    absorbance = absorbance[:, cp.all(absorbance > beta, axis=0)]
+    absorbance = absorbance[:, cp.any(absorbance > beta, axis=0)]
     if absorbance.size == 0:
         raise ValueError(
             "All pixels of the input image are below the threshold."
@@ -383,7 +383,10 @@ def _get_raw_concentrations(src_stain_coeff, absorbance):
 def _normalized_from_concentrations(conc_raw, max_percentile, ref_stain_coeff,
                                     ref_max_conc, source_intensity,
                                     original_shape, channel_axis):
-    """Determine normalized image from concentrations."""
+    """Determine normalized image from concentrations.
+
+    Note: This function will also modify conc_raw in-place.
+    """
 
     # verify conc_raw is shape (2, n_pixels)
     if conc_raw.ndim != 2 or conc_raw.shape[0] != 2:
@@ -410,10 +413,10 @@ def _normalized_from_concentrations(conc_raw, max_percentile, ref_stain_coeff,
          for ch_raw in conc_raw]
     )
     normalization_factors = ref_max_conc / max_conc
-    conc_norm = conc_raw * normalization_factors[:, cp.newaxis]
+    conc_raw *= normalization_factors[:, cp.newaxis]
 
     # reconstruct the image based on the reference stain matrix
-    absorbance_norm = ref_stain_coeff.dot(conc_norm)
+    absorbance_norm = ref_stain_coeff.dot(conc_raw)
     image_norm = absorbance_to_image(
         absorbance_norm, source_intensity=source_intensity, dtype=np.uint8
     )
@@ -527,8 +530,8 @@ def normalize_colors_pca(
     image_norm = _normalized_from_concentrations(
         conc_raw=conc_raw,
         max_percentile=100 - alpha,
-        ref_max_conc=cp.asarray(ref_max_conc),
-        ref_stain_coeff=cp.asarray(ref_stain_coeff),
+        ref_max_conc=cp.asarray(ref_max_conc, dtype=conc_raw.dtype),
+        ref_stain_coeff=cp.asarray(ref_stain_coeff, dtype=conc_raw.dtype),
         source_intensity=source_intensity,
         channel_axis=channel_axis,
         original_shape=image.shape,
diff --git a/python/cucim/src/cucim/core/operations/morphology/__init__.py b/python/cucim/src/cucim/core/operations/morphology/__init__.py
new file mode 100644
index 000000000..a9c676edb
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/__init__.py
@@ -0,0 +1,3 @@
+from ._distance_transform import distance_transform_edt
+
+__all__ = ["distance_transform_edt"]
diff --git a/python/cucim/src/cucim/core/operations/morphology/_distance_transform.py b/python/cucim/src/cucim/core/operations/morphology/_distance_transform.py
new file mode 100644
index 000000000..e781938b2
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/_distance_transform.py
@@ -0,0 +1,182 @@
+import numpy as np
+
+from ._pba_2d import _pba_2d
+from ._pba_3d import _pba_3d
+
+# TODO: support sampling distances
+#       support the distances and indices output arguments
+#       support chamfer, chessboard and l1/manhattan distances too?
+
+
+def distance_transform_edt(image, sampling=None, return_distances=True,
+                           return_indices=False, distances=None, indices=None,
+                           *, block_params=None, float64_distances=False):
+    """Exact Euclidean distance transform.
+
+    This function calculates the distance transform of the `input`, by
+    replacing each foreground (non-zero) element, with its shortest distance to
+    the background (any zero-valued element).
+
+    In addition to the distance transform, the feature transform can be
+    calculated. In this case the index of the closest background element to
+    each foreground element is returned in a separate array.
+
+    Parameters
+    ----------
+    image : array_like
+        Input data to transform. Can be any type but will be converted into
+        binary: 1 wherever image equates to True, 0 elsewhere.
+    sampling : float, or sequence of float, optional
+        Spacing of elements along each dimension. If a sequence, must be of
+        length equal to the image rank; if a single number, this is used for
+        all axes. If not specified, a grid spacing of unity is implied.
+    return_distances : bool, optional
+        Whether to calculate the distance transform.
+    return_indices : bool, optional
+        Whether to calculate the feature transform.
+    distances : float32 cupy.ndarray, optional
+        An output array to store the calculated distance transform, instead of
+        returning it. `return_distances` must be True. It must be the same
+        shape as `image`.
+    indices : int32 cupy.ndarray, optional
+        An output array to store the calculated feature transform, instead of
+        returning it. `return_indicies` must be True. Its shape must be
+        `(image.ndim,) + image.shape`.
+
+    Other Parameters
+    ----------------
+    block_params : 3-tuple of int
+        The m1, m2, m3 algorithm parameters as described in [2]_. If None,
+        suitable defaults will be chosen. Note: This parameter is specific to
+        cuCIM and does not exist in SciPy.
+    float64_distances : bool, optional
+        If True, use double precision in the distance computation (to match
+        SciPy behavior). Otherwise, single precision will be used for
+        efficiency. Note: This parameter is specific to cuCIM and does not
+        exist in SciPy.
+
+    Returns
+    -------
+    distances : float64 ndarray, optional
+        The calculated distance transform. Returned only when
+        `return_distances` is True and `distances` is not supplied. It will
+        have the same shape as `image`.
+    indices : int32 ndarray, optional
+        The calculated feature transform. It has an image-shaped array for each
+        dimension of the image. See example below. Returned only when
+        `return_indices` is True and `indices` is not supplied.
+
+    Notes
+    -----
+    The Euclidean distance transform gives values of the Euclidean distance::
+
+                    n
+      y_i = sqrt(sum (x[i]-b[i])**2)
+                    i
+
+    where b[i] is the background point (value 0) with the smallest Euclidean
+    distance to input points x[i], and n is the number of dimensions.
+
+    Note that the `indices` output may differ from the one given by
+    `scipy.ndimage.distance_transform_edt` in the case of input pixels that are
+    equidistant from multiple background points.
+
+    The parallel banding algorithm implemented here was originally described in
+    [1]_. The kernels used here correspond to the revised PBA+ implementation
+    that is described on the author's website [2]_. The source code of the
+    author's PBA+ implementation is available at [3]_.
+
+    References
+    ----------
+    ..[1] Thanh-Tung Cao, Ke Tang, Anis Mohamed, and Tiow-Seng Tan. 2010.
+        Parallel Banding Algorithm to compute exact distance transform with the
+        GPU. In Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive
+        3D Graphics and Games (I3D ’10). Association for Computing Machinery,
+        New York, NY, USA, 83–90.
+        DOI:https://doi.org/10.1145/1730804.1730818
+    .. [2] https://www.comp.nus.edu.sg/~tants/pba.html
+    .. [3] https://github.com/orzzzjq/Parallel-Banding-Algorithm-plus
+
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from cucim.core.operations import morphology
+    >>> a = cp.array(([0,1,1,1,1],
+    ...               [0,0,1,1,1],
+    ...               [0,1,1,1,1],
+    ...               [0,1,1,1,0],
+    ...               [0,1,1,0,0]))
+    >>> morphology.distance_transform_edt(a)
+    array([[ 0.    ,  1.    ,  1.4142,  2.2361,  3.    ],
+           [ 0.    ,  0.    ,  1.    ,  2.    ,  2.    ],
+           [ 0.    ,  1.    ,  1.4142,  1.4142,  1.    ],
+           [ 0.    ,  1.    ,  1.4142,  1.    ,  0.    ],
+           [ 0.    ,  1.    ,  1.    ,  0.    ,  0.    ]])
+
+    With a sampling of 2 units along x, 1 along y:
+
+    >>> morphology.distance_transform_edt(a, sampling=[2,1])
+    array([[ 0.    ,  1.    ,  2.    ,  2.8284,  3.6056],
+           [ 0.    ,  0.    ,  1.    ,  2.    ,  3.    ],
+           [ 0.    ,  1.    ,  2.    ,  2.2361,  2.    ],
+           [ 0.    ,  1.    ,  2.    ,  1.    ,  0.    ],
+           [ 0.    ,  1.    ,  1.    ,  0.    ,  0.    ]])
+
+    Asking for indices as well:
+
+    >>> edt, inds = morphology.distance_transform_edt(a, return_indices=True)
+    >>> inds
+    array([[[0, 0, 1, 1, 3],
+            [1, 1, 1, 1, 3],
+            [2, 2, 1, 3, 3],
+            [3, 3, 4, 4, 3],
+            [4, 4, 4, 4, 4]],
+           [[0, 0, 1, 1, 4],
+            [0, 1, 1, 1, 4],
+            [0, 0, 1, 4, 4],
+            [0, 0, 3, 3, 4],
+            [0, 0, 3, 3, 4]]])
+
+    """
+    if distances is not None:
+        raise NotImplementedError(
+            "preallocated distances image is not supported"
+        )
+    if indices is not None:
+        raise NotImplementedError(
+            "preallocated indices image is not supported"
+        )
+    scalar_sampling = None
+    if sampling is not None:
+        sampling = np.unique(np.atleast_1d(sampling))
+        if len(sampling) == 1:
+            scalar_sampling = float(sampling)
+            sampling = None
+        else:
+            raise NotImplementedError(
+                "non-uniform values in sampling is not currently supported"
+            )
+
+    if image.ndim == 3:
+        pba_func = _pba_3d
+    elif image.ndim == 2:
+        pba_func = _pba_2d
+    else:
+        raise NotImplementedError(
+            "Only 2D and 3D distance transforms are supported.")
+
+    vals = pba_func(
+        image,
+        sampling=sampling,
+        return_distances=return_distances,
+        return_indices=return_indices,
+        block_params=block_params
+    )
+
+    if return_distances and scalar_sampling is not None:
+        vals = (vals[0] * scalar_sampling,) + vals[1:]
+
+    if len(vals) == 1:
+        vals = vals[0]
+
+    return vals
diff --git a/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py b/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py
new file mode 100644
index 000000000..edb10a983
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/_pba_2d.py
@@ -0,0 +1,300 @@
+import math
+import os
+
+import cupy
+
+pba2d_defines_template = """
+
+// MARKER is used to mark blank pixels in the texture.
+// Any uncolored pixels will have x = MARKER.
+// Input texture should have x = MARKER for all pixels other than sites
+#define MARKER      {marker}
+#define BLOCKSIZE   {block_size_2d}
+#define pixel_int2_t {pixel_int2_t}                // typically short2 (int2 for images with > 32k pixels per side)
+#define make_pixel(x, y)  {make_pixel_func}(x, y)  // typically make_short2 (make_int2 images with > 32k pixels per side
+
+"""  # noqa
+
+
+def _init_marker(int_dtype):
+    """use a minimum value that is appropriate to the integer dtype"""
+    if int_dtype == cupy.int16:
+        # marker = cupy.iinfo(int_dtype).min
+        marker = -32768
+    elif int_dtype == cupy.int32:
+        # divide by two so we don't have to promote other intermediate int
+        # variables to 64-bit int
+        marker = -2147483648 // 2
+    else:
+        raise ValueError(
+            "expected int_dtype to be either cupy.int16 or cupy.int32"
+        )
+    return marker
+
+
+@cupy.memoize(True)
+def get_pba2d_src(block_size_2d=64, marker=-32768, pixel_int2_t='short2'):
+    make_pixel_func = 'make_' + pixel_int2_t
+
+    pba2d_code = pba2d_defines_template.format(
+        block_size_2d=block_size_2d,
+        marker=marker,
+        pixel_int2_t=pixel_int2_t,
+        make_pixel_func=make_pixel_func
+    )
+    kernel_directory = os.path.join(os.path.dirname(__file__), 'cuda')
+    with open(os.path.join(kernel_directory, 'pba_kernels_2d.h'), 'rt') as f:
+        pba2d_kernels = '\n'.join(f.readlines())
+
+    pba2d_code += pba2d_kernels
+    return pba2d_code
+
+
+def _get_block_size(check_warp_size=False):
+    if check_warp_size:
+        dev = cupy.cuda.runtime.getDevice()
+        device_properties = cupy.cuda.runtime.getDeviceProperties(dev)
+        return int(device_properties['warpSize'])
+    else:
+        return 32
+
+
+def _pack_int2(arr, marker=-32768, int_dtype=cupy.int16):
+    if arr.ndim != 2:
+        raise ValueError("only 2d arr suppported")
+    input_x = cupy.zeros(arr.shape, dtype=int_dtype)
+    input_y = cupy.zeros(arr.shape, dtype=int_dtype)
+    # TODO: create custom kernel for setting values in input_x, input_y
+    cond = arr == 0
+    y, x = cupy.where(cond)
+    input_x[cond] = x
+    mask = arr != 0
+    input_x[mask] = marker  # 1 << 32
+    input_y[cond] = y
+    input_y[mask] = marker  # 1 << 32
+    int2_dtype = cupy.dtype({'names': ['x', 'y'], 'formats': [int_dtype] * 2})
+    # in C++ code x is the contiguous axis and corresponds to width
+    #             y is the non-contiguous axis and corresponds to height
+    # given that, store input_x as the last axis here
+    return cupy.squeeze(
+        cupy.stack((input_x, input_y), axis=-1).view(int2_dtype)
+    )
+
+
+def _unpack_int2(img, make_copy=False, int_dtype=cupy.int16):
+    temp = img.view(int_dtype).reshape(img.shape + (2,))
+    if make_copy:
+        temp = temp.copy()
+    return temp
+
+
+def _determine_padding(shape, padded_size, block_size):
+    # all kernels assume equal size along both axes, so pad up to equal size if
+    # shape is not isotropic
+    orig_sy, orig_sx = shape
+    if orig_sx != padded_size or orig_sy != padded_size:
+        padding_width = (
+            (0, padded_size - orig_sy), (0, padded_size - orig_sx)
+        )
+    else:
+        padding_width = None
+    return padding_width
+
+
+def _pba_2d(arr, sampling=None, return_distances=True, return_indices=False,
+            block_params=None, check_warp_size=False, *,
+            float64_distances=False):
+
+    # input_arr: a 2D image
+    #    For each site at (x, y), the pixel at coordinate (x, y) should contain
+    #    the pair (x, y). Pixels that are not sites should contain the pair
+    #    (MARKER, MARKER)
+
+    # Note: could query warp size here, but for now just assume 32 to avoid
+    #       overhead of querying properties
+    block_size = _get_block_size(check_warp_size)
+
+    if sampling is not None:
+        raise NotImplementedError("sampling not yet supported")
+        # if len(sampling) != 2:
+        #     raise ValueError("sampling must be a sequence of two values.")
+
+    padded_size = math.ceil(max(arr.shape) / block_size) * block_size
+    if block_params is None:
+        # should be <= size / 64. sy must be a multiple of m1
+        m1 = max(1, min(padded_size // block_size, 32))
+
+        # size must be a multiple of m2
+        m2 = max(1, min(padded_size // block_size, 32))
+        # m2 must also be a power of two
+        m2 = 2**math.floor(math.log2(m2))
+        if padded_size % m2 != 0:
+            raise RuntimeError("error in setting default m2")
+
+        # should be <= 64. texture size must be a multiple of m3
+        m3 = min(min(m1, m2), 2)
+    else:
+        m1, m2, m3 = block_params
+
+    if m1 > padded_size // block_size:
+        raise ValueError("m1 too large. must be <= arr.shape[0] // 32")
+    if m2 > padded_size // block_size:
+        raise ValueError("m2 too large. must be <= arr.shape[1] // 32")
+    for m in (m1, m2, m3):
+        if padded_size % m != 0:
+            raise ValueError(
+                f"Largest dimension of image ({padded_size}) must be evenly "
+                f"disivible by each element of block_params: {(m1, m2, m3)}."
+            )
+
+    shape_max = max(arr.shape)
+    if shape_max <= 32768:
+        int_dtype = cupy.int16
+        pixel_int2_type = 'short2'
+    else:
+        if shape_max > (1 << 24):
+            # limit to coordinate range to 2**24 due to use of __mul24 in
+            # coordinate TOID macro
+            raise ValueError(
+                f"maximum axis size of {1 << 24} exceeded, for image with "
+                f"shape {arr.shape}"
+            )
+        int_dtype = cupy.int32
+        pixel_int2_type = 'int2'
+
+    marker = _init_marker(int_dtype)
+
+    orig_sy, orig_sx = arr.shape
+    padding_width = _determine_padding(arr.shape, padded_size, block_size)
+    if padding_width is not None:
+        arr = cupy.pad(arr, padding_width, mode='constant', constant_values=1)
+    size = arr.shape[0]
+
+    input_arr = _pack_int2(arr, marker=marker, int_dtype=int_dtype)
+    output = cupy.zeros_like(input_arr)
+
+    int2_dtype = cupy.dtype({'names': ['x', 'y'], 'formats': [int_dtype] * 2})
+    margin = cupy.empty((2 * m1 * size,), dtype=int2_dtype)
+
+    # phase 1 of PBA. m1 must divide texture size and be <= 64
+    pba2d = cupy.RawModule(
+        code=get_pba2d_src(
+            block_size_2d=block_size,
+            marker=marker,
+            pixel_int2_t=pixel_int2_type,
+        )
+    )
+    kernelFloodDown = pba2d.get_function('kernelFloodDown')
+    kernelFloodUp = pba2d.get_function('kernelFloodUp')
+    kernelPropagateInterband = pba2d.get_function('kernelPropagateInterband')
+    kernelUpdateVertical = pba2d.get_function('kernelUpdateVertical')
+    kernelProximatePoints = pba2d.get_function('kernelProximatePoints')
+    kernelCreateForwardPointers = pba2d.get_function(
+        'kernelCreateForwardPointers'
+    )
+    kernelMergeBands = pba2d.get_function('kernelMergeBands')
+    kernelDoubleToSingleList = pba2d.get_function('kernelDoubleToSingleList')
+    kernelColor = pba2d.get_function('kernelColor')
+
+    block = (block_size, 1, 1)
+    grid = (math.ceil(size / block[0]), m1, 1)
+    bandSize1 = size // m1
+    # kernelFloodDown modifies input_arr in-place
+    kernelFloodDown(
+        grid,
+        block,
+        (input_arr, input_arr, size, bandSize1),
+    )
+    # kernelFloodUp modifies input_arr in-place
+    kernelFloodUp(
+        grid,
+        block,
+        (input_arr, input_arr, size, bandSize1),
+    )
+    # kernelFloodUp fills values into margin
+    kernelPropagateInterband(
+        grid,
+        block,
+        (input_arr, margin, size, bandSize1),
+    )
+    # kernelUpdateVertical stores output into an intermediate array of
+    # transposed shape
+    kernelUpdateVertical(
+        grid,
+        block,
+        (input_arr, margin, output, size, bandSize1),
+    )
+
+    # phase 2
+    block = (block_size, 1, 1)
+    grid = (math.ceil(size / block[0]), m2, 1)
+    bandSize2 = size // m2
+    kernelProximatePoints(
+        grid,
+        block,
+        (output, input_arr, size, bandSize2),
+    )
+    kernelCreateForwardPointers(
+        grid,
+        block,
+        (input_arr, input_arr, size, bandSize2),
+    )
+    # Repeatly merging two bands into one
+    noBand = m2
+    while noBand > 1:
+        grid = (math.ceil(size / block[0]), noBand // 2)
+        kernelMergeBands(
+            grid,
+            block,
+            (output, input_arr, input_arr, size, size // noBand),
+        )
+        noBand //= 2
+    # Replace the forward link with the X coordinate of the seed to remove
+    # the need of looking at the other texture. We need it for coloring.
+    grid = (math.ceil(size / block[0]), size)
+    kernelDoubleToSingleList(
+        grid,
+        block,
+        (output, input_arr, input_arr, size),
+    )
+
+    # Phase 3 of PBA
+    block = (block_size, m3, 1)
+    grid = (math.ceil(size / block[0]), 1, 1)
+    kernelColor(
+        grid,
+        block,
+        (input_arr, output, size),
+    )
+
+    output = _unpack_int2(output, make_copy=False, int_dtype=int_dtype)
+    # make sure to crop any padding that was added here!
+    x = output[:orig_sy, :orig_sx, 0]
+    y = output[:orig_sy, :orig_sx, 1]
+
+    # raise NotImplementedError("TODO")
+    vals = ()
+    if return_distances:
+        # TODO: custom kernel for more efficient distance computation
+        y0, x0 = cupy.meshgrid(
+            *(
+                cupy.arange(s, dtype=cupy.int32)
+                for s in (orig_sy, orig_sx)
+            ),
+            indexing='ij',
+            sparse=True,
+        )
+        tmp = (x - x0)
+        dist = tmp * tmp
+        tmp = (y - y0)
+        dist += tmp * tmp
+        if float64_distances:
+            dist = cupy.sqrt(dist)
+        else:
+            dist = dist.astype(cupy.float32)
+            cupy.sqrt(dist, out=dist)
+        vals = vals + (dist,)
+    if return_indices:
+        indices = cupy.stack((y, x), axis=0)
+        vals = vals + (indices,)
+    return vals
diff --git a/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py b/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py
new file mode 100644
index 000000000..dab484e90
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/_pba_3d.py
@@ -0,0 +1,311 @@
+import functools
+import math
+import numbers
+import os
+
+import cupy
+import numpy as np
+
+from ._pba_2d import _get_block_size
+
+try:
+    # math.lcm was introduced in Python 3.9
+    from math import lcm
+except ImportError:
+
+    """Fallback implementation of least common multiple (lcm)
+
+    TODO: remove once minimum Python requirement is >= 3.9
+    """
+
+    def _lcm(a, b):
+        return abs(b * (a // math.gcd(a, b)))
+
+    @functools.lru_cache()
+    def lcm(*args):
+        nargs = len(args)
+        if not all(isinstance(a, numbers.Integral) for a in args):
+            raise TypeError("all arguments must be integers")
+        if nargs == 0:
+            return 1
+        res = int(args[0])
+        if nargs == 1:
+            return abs(res)
+        for i in range(1, nargs):
+            x = int(args[i])
+            res = _lcm(res, x)
+        return res
+
+
+pba3d_defines_template = """
+
+#define MARKER     {marker}
+#define MAX_INT    {max_int}
+#define BLOCKSIZE  {block_size_3d}
+
+"""
+
+# For efficiency, the original PBA+ packs three 10-bit integers and two binary
+# flags into a single 32-bit integer. The defines in
+# `pba3d_defines_encode_32bit` handle this format.
+pba3d_defines_encode_32bit = """
+// Sites     : ENCODE(x, y, z, 0, 0)
+// Not sites : ENCODE(0, 0, 0, 1, 0) or MARKER
+#define ENCODED_INT_TYPE int
+#define ZERO 0
+#define ONE 1
+#define ENCODE(x, y, z, a, b)  (((x) << 20) | ((y) << 10) | (z) | ((a) << 31) | ((b) << 30))
+#define DECODE(value, x, y, z) \
+    x = ((value) >> 20) & 0x3ff; \
+    y = ((value) >> 10) & 0x3ff; \
+    z = (value) & 0x3ff
+
+#define NOTSITE(value)  (((value) >> 31) & 1)
+#define HASNEXT(value)  (((value) >> 30) & 1)
+
+#define GET_X(value)    (((value) >> 20) & 0x3ff)
+#define GET_Y(value)    (((value) >> 10) & 0x3ff)
+#define GET_Z(value)    ((NOTSITE((value))) ? MAX_INT : ((value) & 0x3ff))
+
+""" # noqa
+
+
+# 64bit version of ENCODE/DECODE to allow a 20-bit integer per coordinate axis.
+pba3d_defines_encode_64bit = """
+// Sites     : ENCODE(x, y, z, 0, 0)
+// Not sites : ENCODE(0, 0, 0, 1, 0) or MARKER
+#define ENCODED_INT_TYPE long long
+#define ZERO 0L
+#define ONE 1L
+#define ENCODE(x, y, z, a, b)  (((x) << 40) | ((y) << 20) | (z) | ((a) << 61) | ((b) << 60))
+#define DECODE(value, x, y, z) \
+    x = ((value) >> 40) & 0xfffff; \
+    y = ((value) >> 20) & 0xfffff; \
+    z = (value) & 0xfffff
+
+#define NOTSITE(value)  (((value) >> 61) & 1)
+#define HASNEXT(value)  (((value) >> 60) & 1)
+
+#define GET_X(value)    (((value) >> 40) & 0xfffff)
+#define GET_Y(value)    (((value) >> 20) & 0xfffff)
+#define GET_Z(value)    ((NOTSITE((value))) ? MAX_INT : ((value) & 0xfffff))
+
+""" # noqa
+
+
+@cupy.memoize(True)
+def get_pba3d_src(block_size_3d=32, marker=-2147483648, max_int=2147483647,
+                  size_max=1024):
+    pba3d_code = pba3d_defines_template.format(
+        block_size_3d=block_size_3d, marker=marker, max_int=max_int
+    )
+    if size_max > 1024:
+        pba3d_code += pba3d_defines_encode_64bit
+    else:
+        pba3d_code += pba3d_defines_encode_32bit
+    kernel_directory = os.path.join(os.path.dirname(__file__), 'cuda')
+    with open(os.path.join(kernel_directory, 'pba_kernels_3d.h'), 'rt') as f:
+        pba3d_kernels = '\n'.join(f.readlines())
+    pba3d_code += pba3d_kernels
+    return pba3d_code
+
+
+# TODO: custom kernel for encode3d
+def encode3d(arr, marker=-2147483648, bit_depth=32, size_max=1024):
+    if arr.ndim != 3:
+        raise ValueError("only 3d arr suppported")
+    if bit_depth not in [32, 64]:
+        raise ValueError("only bit_depth of 32 or 64 is supported")
+    if size_max > 1024:
+        dtype = np.int64
+    else:
+        dtype = np.int32
+    image = cupy.zeros(arr.shape, dtype=dtype, order='C')
+    cond = arr == 0
+    z, y, x = cupy.where(cond)
+    # z, y, x so that x is the contiguous axis
+    # (must match TOID macro in the C++ code!)
+    if size_max > 1024:
+        image[cond] = (((x) << 40) | ((y) << 20) | (z))
+    else:
+        image[cond] = (((x) << 20) | ((y) << 10) | (z))
+    image[arr != 0] = marker  # 1 << 32
+    return image
+
+
+# TODO: custom kernel for decode3d
+def decode3d(output, size_max=1024):
+    if size_max > 1024:
+        x = (output >> 40) & 0xfffff
+        y = (output >> 20) & 0xfffff
+        z = output & 0xfffff
+    else:
+        x = (output >> 20) & 0x3ff
+        y = (output >> 10) & 0x3ff
+        z = output & 0x3ff
+    return (x, y, z)
+
+
+def _determine_padding(shape, block_size, m1, m2, m3, blockx, blocky):
+    # TODO: can possibly revise to consider only particular factors for LCM on
+    #       a given axis
+    LCM = lcm(block_size, m1, m2, m3, blockx, blocky)
+    orig_sz, orig_sy, orig_sx = shape
+    round_up = False
+    if orig_sx % LCM != 0:
+        # round up size to a multiple of the band size
+        round_up = True
+        sx = LCM * math.ceil(orig_sx / LCM)
+    else:
+        sx = orig_sx
+    if orig_sy % LCM != 0:
+        # round up size to a multiple of the band size
+        round_up = True
+        sy = LCM * math.ceil(orig_sy / LCM)
+    else:
+        sy = orig_sy
+    if orig_sz % LCM != 0:
+        # round up size to a multiple of the band size
+        round_up = True
+        sz = LCM * math.ceil(orig_sz / LCM)
+    else:
+        sz = orig_sz
+
+    aniso = not (sx == sy == sz)
+    if aniso or round_up:
+        smax = max(sz, sy, sx)
+        padding_width = (
+            (0, smax - orig_sz), (0, smax - orig_sy), (0, smax - orig_sx)
+        )
+    else:
+        padding_width = None
+    return padding_width
+
+
+def _pba_3d(arr, sampling=None, return_distances=True, return_indices=False,
+            block_params=None, check_warp_size=False, *,
+            float64_distances=False):
+    if arr.ndim != 3:
+        raise ValueError(f"expected a 3D array, got {arr.ndim}D")
+
+    if sampling is not None:
+        raise NotImplementedError("sampling not yet supported")
+        # if len(sampling) != 3:
+        #     raise ValueError("sampling must be a sequence of three values.")
+
+    if block_params is None:
+        m1 = 1
+        m2 = 1
+        m3 = 2
+    else:
+        m1, m2, m3 = block_params
+
+    # reduce blockx for small inputs
+    s_min = min(arr.shape)
+    if s_min <= 4:
+        blockx = 4
+    elif s_min <= 8:
+        blockx = 8
+    elif s_min <= 16:
+        blockx = 16
+    else:
+        blockx = 32
+    blocky = 4
+
+    block_size = _get_block_size(check_warp_size)
+
+    orig_sz, orig_sy, orig_sx = arr.shape
+    padding_width = _determine_padding(
+        arr.shape, block_size, m1, m2, m3, blockx, blocky
+    )
+    if padding_width is not None:
+        arr = cupy.pad(arr, padding_width, mode='constant', constant_values=1)
+    size = arr.shape[0]
+
+    # pba algorithm was implemented to use 32-bit integer to store compressed
+    # coordinates. input_arr will be C-contiguous, int32
+    size_max = max(arr.shape)
+    input_arr = encode3d(arr, size_max=size_max)
+    buffer_idx = 0
+    output = cupy.zeros_like(input_arr)
+    pba_images = [input_arr, output]
+
+    block = (blockx, blocky, 1)
+    grid = (size // block[0], size // block[1], 1)
+    pba3d = cupy.RawModule(
+        code=get_pba3d_src(block_size_3d=block_size, size_max=size_max)
+    )
+
+    kernelFloodZ = pba3d.get_function('kernelFloodZ')
+    kernelMaurerAxis = pba3d.get_function('kernelMaurerAxis')
+    kernelColorAxis = pba3d.get_function('kernelColorAxis')
+
+    kernelFloodZ(
+        grid,
+        block,
+        (pba_images[buffer_idx], pba_images[1 - buffer_idx], size)
+    )
+    buffer_idx = 1 - buffer_idx
+
+    block = (blockx, blocky, 1)
+    grid = (size // block[0], size // block[1], 1)
+    kernelMaurerAxis(
+        grid,
+        block,
+        (pba_images[buffer_idx], pba_images[1 - buffer_idx], size),
+    )
+
+    block = (block_size, m3, 1)
+    grid = (size // block[0], size, 1)
+    kernelColorAxis(
+        grid,
+        block,
+        (pba_images[1 - buffer_idx], pba_images[buffer_idx], size),
+    )
+
+    block = (blockx, blocky, 1)
+    grid = (size // block[0], size // block[1], 1)
+    kernelMaurerAxis(
+        grid,
+        block,
+        (pba_images[buffer_idx], pba_images[1 - buffer_idx], size),
+    )
+
+    block = (block_size, m3, 1)
+    grid = (size // block[0], size, 1)
+    kernelColorAxis(
+        grid,
+        block,
+        (pba_images[1 - buffer_idx], pba_images[buffer_idx], size),
+    )
+
+    output = pba_images[buffer_idx]
+    if return_distances or return_indices:
+        x, y, z = decode3d(output[:orig_sz, :orig_sy, :orig_sx],
+                           size_max=size_max)
+
+    vals = ()
+    if return_distances:
+        # TODO: custom kernel for more efficient distance computation
+        orig_shape = (orig_sz, orig_sy, orig_sx)
+        z0, y0, x0 = cupy.meshgrid(
+            *(cupy.arange(s, dtype=cupy.int32) for s in orig_shape),
+            indexing='ij',
+            sparse=True
+        )
+        tmp = (x - x0)
+        dist = tmp * tmp
+        tmp = (y - y0)
+        dist += tmp * tmp
+        tmp = (z - z0)
+        dist += tmp * tmp
+        if float64_distances:
+            dist = cupy.sqrt(dist)
+        else:
+            dist = dist.astype(cupy.float32)
+            cupy.sqrt(dist, out=dist)
+        vals = vals + (dist,)
+    if return_indices:
+        indices = cupy.stack((z, y, x), axis=0)
+        vals = vals + (indices,)
+    return vals
diff --git a/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_2d.h b/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_2d.h
new file mode 100644
index 000000000..61677c682
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_2d.h
@@ -0,0 +1,451 @@
+// Euclidean Distance Transform
+// 
+// Kernels for the 2D version of the Parallel Banding Algorithm (PBA+). 
+// 
+// MIT license: see 3rdparty/LICENSE.pba+
+// Copyright: (c) 2019 School of Computing, National University of Singapore
+//
+// Modifications by Gregory Lee (2022) (NVIDIA)
+// - add user-defined pixel_int2_t to enable
+// - replace __mul24 operations with standard multiplication operator
+
+
+// START OF DEFINITIONS OVERRIDDEN BY THE PYTHON SCRIPT
+
+// The values included in this header file are those defined in the original
+// PBA+ implementation
+
+// However, the Python code generation can potentially generate a different
+// ENCODE/DECODE that use 20 bits per coordinates instead of 10 bits per
+// coordinate with ENCODED_INT_TYPE as `long long`.
+
+#ifndef MARKER
+#define MARKER -32768
+#endif
+
+#ifndef BLOCKSIZE
+#define BLOCKSIZE 32
+#endif
+
+#ifndef pixel_int2_t
+#define pixel_int2_t short2
+#define make_pixel(x, y)  make_short2(x, y)
+#endif
+
+// END OF DEFINITIONS OVERRIDDEN BY THE PYTHON SCRIPT
+
+
+#define TOID(x, y, size)  ((y) * (size) + (x))
+
+#define LL long long
+__device__ bool dominate(LL x1, LL y1, LL x2, LL y2, LL x3, LL y3, LL x0)
+{
+    LL k1 = y2 - y1, k2 = y3 - y2;
+    return (k1 * (y1 + y2) + (x2 - x1) * ((x1 + x2) - (x0 << 1))) * k2 > \
+            (k2 * (y2 + y3) + (x3 - x2) * ((x2 + x3) - (x0 << 1))) * k1;
+}
+#undef LL
+
+
+extern "C"{
+
+__global__ void kernelFloodDown(pixel_int2_t *input, pixel_int2_t *output, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * bandSize;
+    int id = TOID(tx, ty, size);
+
+    pixel_int2_t pixel1, pixel2;
+
+    pixel1 = make_pixel(MARKER, MARKER);
+
+    for (int i = 0; i < bandSize; i++, id += size) {
+        pixel2 = input[id];
+
+        if (pixel2.x != MARKER)
+            pixel1 = pixel2;
+
+        output[id] = pixel1;
+    }
+}
+
+__global__ void kernelFloodUp(pixel_int2_t *input, pixel_int2_t *output, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = (blockIdx.y+1) * bandSize - 1;
+    int id = TOID(tx, ty, size);
+
+    pixel_int2_t pixel1, pixel2;
+    int dist1, dist2;
+
+    pixel1 = make_pixel(MARKER, MARKER);
+
+    for (int i = 0; i < bandSize; i++, id -= size) {
+        dist1 = abs(pixel1.y - ty + i);
+
+        pixel2 = input[id];
+        dist2 = abs(pixel2.y - ty + i);
+
+        if (dist2 < dist1)
+            pixel1 = pixel2;
+
+        output[id] = pixel1;
+    }
+}
+
+__global__ void kernelPropagateInterband(pixel_int2_t *input, pixel_int2_t *margin_out, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int inc = bandSize * size;
+    int ny, nid, nDist;
+    pixel_int2_t pixel;
+
+    // Top row, look backward
+    int ty = blockIdx.y * bandSize;
+    int topId = TOID(tx, ty, size);
+    int bottomId = TOID(tx, ty + bandSize - 1, size);
+    int tid = blockIdx.y * size + tx;
+    int bid = tid + (size * size / bandSize);
+
+    pixel = input[topId];
+    int myDist = abs(pixel.y - ty);
+    margin_out[tid] = pixel;
+
+    for (nid = bottomId - inc; nid >= 0; nid -= inc) {
+        pixel = input[nid];
+
+        if (pixel.x != MARKER) {
+            nDist = abs(pixel.y - ty);
+
+            if (nDist < myDist)
+                margin_out[tid] = pixel;
+
+            break;
+        }
+    }
+
+    // Last row, look downward
+    ty = ty + bandSize - 1;
+    pixel = input[bottomId];
+    myDist = abs(pixel.y - ty);
+    margin_out[bid] = pixel;
+
+    for (ny = ty + 1, nid = topId + inc; ny < size; ny += bandSize, nid += inc) {
+        pixel = input[nid];
+
+        if (pixel.x != MARKER) {
+            nDist = abs(pixel.y - ty);
+
+            if (nDist < myDist)
+                margin_out[bid] = pixel;
+
+            break;
+        }
+    }
+}
+
+__global__ void kernelUpdateVertical(pixel_int2_t *color, pixel_int2_t *margin, pixel_int2_t *output, int size, int bandSize)
+{
+    __shared__ pixel_int2_t block[BLOCKSIZE][BLOCKSIZE];
+
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * bandSize;
+
+    pixel_int2_t top = margin[blockIdx.y * size + tx];
+    pixel_int2_t bottom = margin[(blockIdx.y + size / bandSize) * size + tx];
+    pixel_int2_t pixel;
+
+    int dist, myDist;
+
+    int id = TOID(tx, ty, size);
+
+    int n_step = bandSize / blockDim.x;
+    for(int step = 0; step < n_step; ++step) {
+        int y_start = blockIdx.y * bandSize + step * blockDim.x;
+        int y_end = y_start + blockDim.x;
+
+        for (ty = y_start; ty < y_end; ++ty, id += size) {
+            pixel = color[id];
+            myDist = abs(pixel.y - ty);
+
+            dist = abs(top.y - ty);
+            if (dist < myDist) { myDist = dist; pixel = top; }
+
+            dist = abs(bottom.y - ty);
+            if (dist < myDist) pixel = bottom;
+
+            // temporary result is stored in block
+            block[threadIdx.x][ty - y_start] = make_pixel(pixel.y, pixel.x);
+        }
+
+        __syncthreads();
+
+        // block is written to a transposed location in the output
+
+        int tid = TOID(blockIdx.y * bandSize + step * blockDim.x + threadIdx.x, \
+                        blockIdx.x * blockDim.x, size);
+
+        for(int i = 0; i < blockDim.x; ++i, tid += size) {
+            output[tid] = block[i][threadIdx.x];
+        }
+
+        __syncthreads();
+    }
+}
+
+__global__ void kernelProximatePoints(pixel_int2_t *input, pixel_int2_t *stack, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * bandSize;
+    int id = TOID(tx, ty, size);
+    int lasty = -1;
+    pixel_int2_t last1, last2, current;
+
+    last1.y = -1; last2.y = -1;
+
+    for (int i = 0; i < bandSize; i++, id += size) {
+        current = input[id];
+
+        if (current.x != MARKER) {
+            while (last2.y >= 0) {
+                if (!dominate(last1.x, last2.y, last2.x, \
+                    lasty, current.x, current.y, tx))
+                    break;
+
+                lasty = last2.y; last2 = last1;
+
+                if (last1.y >= 0)
+                    last1 = stack[TOID(tx, last1.y, size)];
+            }
+
+            last1 = last2; last2 = make_pixel(current.x, lasty); lasty = current.y;
+
+            stack[id] = last2;
+        }
+    }
+
+    // Store the pointer to the tail at the last pixel of this band
+    if (lasty != ty + bandSize - 1)
+        stack[TOID(tx, ty + bandSize - 1, size)] = make_pixel(MARKER, lasty);
+}
+
+__global__ void kernelCreateForwardPointers(pixel_int2_t *input, pixel_int2_t *output, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = (blockIdx.y+1) * bandSize - 1;
+    int id = TOID(tx, ty, size);
+    int lasty = -1, nexty;
+    pixel_int2_t current;
+
+    // Get the tail pointer
+    current = input[id];
+
+    if (current.x == MARKER)
+        nexty = current.y;
+    else
+        nexty = ty;
+
+    for (int i = 0; i < bandSize; i++, id -= size)
+        if (ty - i == nexty) {
+            current = make_pixel(lasty, input[id].y);
+            output[id] = current;
+
+            lasty = nexty;
+            nexty = current.y;
+        }
+
+    // Store the pointer to the head at the first pixel of this band
+    if (lasty != ty - bandSize + 1)
+        output[id + size] = make_pixel(lasty, MARKER);
+}
+
+__global__ void kernelMergeBands(pixel_int2_t *color, pixel_int2_t *link, pixel_int2_t *output, int size, int bandSize)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int band1 = blockIdx.y * 2;
+    int band2 = band1 + 1;
+    int firsty, lasty;
+    pixel_int2_t last1, last2, current;
+    // last1 and last2: x component store the x coordinate of the site,
+    // y component store the backward pointer
+    // current: y component store the x coordinate of the site,
+    // x component store the forward pointer
+
+    // Get the two last items of the first list
+    lasty = band2 * bandSize - 1;
+    last2 = make_pixel(color[TOID(tx, lasty, size)].x,
+        link[TOID(tx, lasty, size)].y);
+
+    if (last2.x == MARKER) {
+        lasty = last2.y;
+
+        if (lasty >= 0)
+            last2 = make_pixel(color[TOID(tx, lasty, size)].x,
+            link[TOID(tx, lasty, size)].y);
+        else
+            last2 = make_pixel(MARKER, MARKER);
+    }
+
+    if (last2.y >= 0) {
+        // Second item at the top of the stack
+        last1 = make_pixel(color[TOID(tx, last2.y, size)].x,
+            link[TOID(tx, last2.y, size)].y);
+    }
+
+    // Get the first item of the second band
+    firsty = band2 * bandSize;
+    current = make_pixel(link[TOID(tx, firsty, size)].x,
+        color[TOID(tx, firsty, size)].x);
+
+    if (current.y == MARKER) {
+        firsty = current.x;
+
+        if (firsty >= 0)
+            current = make_pixel(link[TOID(tx, firsty, size)].x,
+            color[TOID(tx, firsty, size)].x);
+        else
+            current = make_pixel(MARKER, MARKER);
+    }
+
+    // Count the number of item in the second band that survive so far.
+    // Once it reaches 2, we can stop.
+    int top = 0;
+
+    while (top < 2 && current.y >= 0) {
+        // While there's still something on the left
+        while (last2.y >= 0) {
+
+            if (!dominate(last1.x, last2.y, last2.x, \
+                lasty, current.y, firsty, tx))
+                break;
+
+            lasty = last2.y; last2 = last1;
+            top--;
+
+            if (last1.y >= 0)
+                last1 = make_pixel(color[TOID(tx, last1.y, size)].x,
+                link[TOID(tx, last1.y, size)].y);
+        }
+
+        // Update the current pointer
+        output[TOID(tx, firsty, size)] = make_pixel(current.x, lasty);
+
+        if (lasty >= 0)
+            output[TOID(tx, lasty, size)] = make_pixel(firsty, last2.y);
+
+        last1 = last2; last2 = make_pixel(current.y, lasty); lasty = firsty;
+        firsty = current.x;
+
+        top = max(1, top + 1);
+
+        // Advance the current pointer to the next one
+        if (firsty >= 0)
+            current = make_pixel(link[TOID(tx, firsty, size)].x,
+            color[TOID(tx, firsty, size)].x);
+        else
+            current = make_pixel(MARKER, MARKER);
+    }
+
+    // Update the head and tail pointer.
+    firsty = band1 * bandSize;
+    lasty = band2 * bandSize;
+    current = link[TOID(tx, firsty, size)];
+
+    if (current.y == MARKER && current.x < 0) { // No head?
+        last1 = link[TOID(tx, lasty, size)];
+
+        if (last1.y == MARKER)
+            current.x = last1.x;
+        else
+            current.x = lasty;
+
+        output[TOID(tx, firsty, size)] = current;
+    }
+
+    firsty = band1 * bandSize + bandSize - 1;
+    lasty = band2 * bandSize + bandSize - 1;
+    current = link[TOID(tx, lasty, size)];
+
+    if (current.x == MARKER && current.y < 0) { // No tail?
+        last1 = link[TOID(tx, firsty, size)];
+
+        if (last1.x == MARKER)
+            current.y = last1.y;
+        else
+            current.y = firsty;
+
+        output[TOID(tx, lasty, size)] = current;
+    }
+}
+
+__global__ void kernelDoubleToSingleList(pixel_int2_t *color, pixel_int2_t *link, pixel_int2_t *output, int size)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y;
+    int id = TOID(tx, ty, size);
+
+    output[id] = make_pixel(color[id].x, link[id].y);
+}
+
+__global__ void kernelColor(pixel_int2_t *input, pixel_int2_t *output, int size)
+{
+    __shared__ pixel_int2_t block[BLOCKSIZE][BLOCKSIZE];
+
+    int col = threadIdx.x;
+    int tid = threadIdx.y;
+    int tx = blockIdx.x * blockDim.x + col;
+    int dx, dy, lasty;
+    unsigned int best, dist;
+    pixel_int2_t last1, last2;
+
+    lasty = size - 1;
+
+    last2 = input[TOID(tx, lasty, size)];
+
+    if (last2.x == MARKER) {
+        lasty = last2.y;
+        last2 = input[TOID(tx, lasty, size)];
+    }
+
+    if (last2.y >= 0)
+        last1 = input[TOID(tx, last2.y, size)];
+
+    int y_start, y_end, n_step = size / blockDim.x;
+    for(int step = 0; step < n_step; ++step) {
+        y_start = size - step * blockDim.x - 1;
+        y_end = size - (step + 1) * blockDim.x;
+
+        for (int ty = y_start - tid; ty >= y_end; ty -= blockDim.y) {
+            dx = last2.x - tx; dy = lasty - ty;
+            best = dist = dx * dx + dy * dy;
+
+            while (last2.y >= 0) {
+                dx = last1.x - tx; dy = last2.y - ty;
+                dist = dx * dx + dy * dy;
+
+                if (dist > best)
+                    break;
+
+                best = dist; lasty = last2.y; last2 = last1;
+
+                if (last2.y >= 0)
+                    last1 = input[TOID(tx, last2.y, size)];
+            }
+
+            block[threadIdx.x][ty - y_end] = make_pixel(lasty, last2.x);
+        }
+
+        __syncthreads();
+
+        // note: transposes back to original shape here
+        if(!threadIdx.y) {
+            int id = TOID(y_end + threadIdx.x, blockIdx.x * blockDim.x, size);
+            for(int i = 0; i < blockDim.x; ++i, id+=size) {
+                output[id] = block[i][threadIdx.x];
+            }
+        }
+
+        __syncthreads();
+    }
+}
+} // extern C
diff --git a/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_3d.h b/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_3d.h
new file mode 100644
index 000000000..c09f4b51f
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/cuda/pba_kernels_3d.h
@@ -0,0 +1,237 @@
+// Euclidean Distance Transform
+// 
+// Kernels for the 3D version of the Parallel Banding Algorithm (PBA+). 
+// 
+// MIT license: see 3rdparty/LICENSE.pba+
+//
+// Modifications by Gregory Lee (2022) (NVIDIA)
+// - allow user-defined ENCODED_INT_TYPE, ENCODE, DECODE
+
+
+// START OF DEFINITIONS OVERRIDDEN BY THE PYTHON SCRIPT
+
+// The values included in this header file are those defined in the original
+// PBA+ implementation
+
+// However, the Python code generation can potentially generate a different
+// ENCODE/DECODE that use 20 bits per coordinates instead of 10 bits per
+// coordinate with ENCODED_INT_TYPE as `long long`.
+
+
+#ifndef MARKER
+#define MARKER     -2147483648
+#endif  // MARKER
+
+#ifndef MAX_INT
+#define MAX_INT    2147483647
+#endif
+
+#ifndef BLOCKSIZE
+#define BLOCKSIZE  32
+#endif
+
+#ifndef ENCODE
+
+// Sites     : ENCODE(x, y, z, 0, 0)
+// Not sites : ENCODE(0, 0, 0, 1, 0) or MARKER
+#define ENCODED_INT_TYPE int
+#define ZERO 0
+#define ONE 1
+#define ENCODE(x, y, z, a, b)  (((x) << 20) | ((y) << 10) | (z) | ((a) << 31) | ((b) << 30))
+#define DECODE(value, x, y, z) \
+    x = ((value) >> 20) & 0x3ff; \
+    y = ((value) >> 10) & 0x3ff; \
+    z = (value) & 0x3ff
+
+#define NOTSITE(value)  (((value) >> 31) & 1)
+#define HASNEXT(value)  (((value) >> 30) & 1)
+
+#define GET_X(value)    (((value) >> 20) & 0x3ff)
+#define GET_Y(value)    (((value) >> 10) & 0x3ff)
+#define GET_Z(value)    ((NOTSITE((value))) ? MAX_INT : ((value) & 0x3ff))
+
+#endif // ENCODE
+
+// END OF DEFINITIONS DEFINED IN THE PYTHON SCRIPT
+
+
+#define LL long long
+__device__ bool dominate(LL x_1, LL y_1, LL z_1, LL x_2, LL y_2, LL z_2, LL x_3, LL y_3, LL z_3, LL x_0, LL z_0)
+{
+    LL k_1 = y_2 - y_1, k_2 = y_3 - y_2;
+
+    return (((y_1 + y_2) * k_1 + ((x_2 - x_1) * (x_1 + x_2 - (x_0 << 1)) + (z_2 - z_1) * (z_1 + z_2 - (z_0 << 1)))) * k_2 > \
+            ((y_2 + y_3) * k_2 + ((x_3 - x_2) * (x_2 + x_3 - (x_0 << 1)) + (z_3 - z_2) * (z_2 + z_3 - (z_0 << 1)))) * k_1);
+}
+#undef LL
+
+#define TOID(x, y, z, size)    ((((z) * (size)) + (y)) * (size) + (x))
+
+
+extern "C"{
+
+__global__ void kernelFloodZ(ENCODED_INT_TYPE *input, ENCODED_INT_TYPE *output, int size)
+{
+
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+    int tz = 0;
+
+    int plane = size * size;
+    int id = TOID(tx, ty, tz, size);
+    ENCODED_INT_TYPE pixel1, pixel2;
+
+    pixel1 = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO);
+
+    // Sweep down
+    for (int i = 0; i < size; i++, id += plane) {
+        pixel2 = input[id];
+
+        if (!NOTSITE(pixel2))
+            pixel1 = pixel2;
+
+        output[id] = pixel1;
+    }
+
+    ENCODED_INT_TYPE dist1, dist2, nz;
+
+    id -= plane + plane;
+
+    // Sweep up
+    for (int i = size - 2; i >= 0; i--, id -= plane) {
+        nz = GET_Z(pixel1);
+        dist1 = abs(nz - (tz + i));
+
+        pixel2 = output[id];
+        nz = GET_Z(pixel2);
+        dist2 = abs(nz - (tz + i));
+
+        if (dist2 < dist1)
+            pixel1 = pixel2;
+
+        output[id] = pixel1;
+    }
+}
+
+
+__global__ void kernelMaurerAxis(ENCODED_INT_TYPE *input, ENCODED_INT_TYPE *stack, int size)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int tz = blockIdx.y * blockDim.y + threadIdx.y;
+    int ty = 0;
+
+    int id = TOID(tx, ty, tz, size);
+
+    ENCODED_INT_TYPE lasty = 0;
+    ENCODED_INT_TYPE x1, y1, z1, x2, y2, z2, nx, ny, nz;
+    ENCODED_INT_TYPE p = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO), s1 = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO), s2 = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO);
+    ENCODED_INT_TYPE flag = 0;
+
+    for (ty = 0; ty < size; ++ty, id += size) {
+        p = input[id];
+
+        if (!NOTSITE(p)) {
+
+            while (HASNEXT(s2)) {
+                DECODE(s1, x1, y1, z1);
+                DECODE(s2, x2, y2, z2);
+                DECODE(p, nx, ny, nz);
+
+                if (!dominate(x1, y2, z1, x2, lasty, z2, nx, ty, nz, tx, tz))
+                    break;
+
+                lasty = y2; s2 = s1; y2 = y1;
+
+                if (HASNEXT(s2))
+                    s1 = stack[TOID(tx, y2, tz, size)];
+            }
+
+            DECODE(p, nx, ny, nz);
+            s1 = s2;
+            s2 = ENCODE(nx, lasty, nz, ZERO, flag);
+            y2 = lasty;
+            lasty = ty;
+
+            stack[id] = s2;
+
+            flag = ONE;
+        }
+    }
+
+    if (NOTSITE(p))
+        stack[TOID(tx, ty - 1, tz, size)] = ENCODE(ZERO, lasty, ZERO, ONE, flag);
+}
+
+__global__ void kernelColorAxis(ENCODED_INT_TYPE *input, ENCODED_INT_TYPE *output, int size)
+{
+    __shared__ ENCODED_INT_TYPE block[BLOCKSIZE][BLOCKSIZE];
+
+    int col = threadIdx.x;
+    int tid = threadIdx.y;
+    int tx = blockIdx.x * blockDim.x + col;
+    int tz = blockIdx.y;
+
+    ENCODED_INT_TYPE x1, y1, z1, x2, y2, z2;
+    ENCODED_INT_TYPE last1 = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO), last2 = ENCODE(ZERO,ZERO,ZERO,ONE,ZERO), lasty;
+    long long dx, dy, dz, best, dist;
+
+    lasty = size - 1;
+
+    last2 = input[TOID(tx, lasty, tz, size)];
+    DECODE(last2, x2, y2, z2);
+
+    if (NOTSITE(last2)) {
+        lasty = y2;
+        if(HASNEXT(last2)) {
+            last2 = input[TOID(tx, lasty, tz, size)];
+            DECODE(last2, x2, y2, z2);
+        }
+    }
+
+    if (HASNEXT(last2)) {
+        last1 = input[TOID(tx, y2, tz, size)];
+        DECODE(last1, x1, y1, z1);
+    }
+
+    int y_start, y_end, n_step = size / blockDim.x;
+    for(int step = 0; step < n_step; ++step) {
+        y_start = size - step * blockDim.x - 1;
+        y_end = size - (step + 1) * blockDim.x;
+
+        for (int ty = y_start - tid; ty >= y_end; ty -= blockDim.y) {
+            dx = x2 - tx; dy = lasty - ty; dz = z2 - tz;
+            best = dx * dx + dy * dy + dz * dz;
+
+            while (HASNEXT(last2)) {
+                dx = x1 - tx; dy = y2 - ty; dz = z1 - tz;
+                dist = dx * dx + dy * dy + dz * dz;
+
+                if(dist > best) break;
+
+                best = dist; lasty = y2; last2 = last1;
+                DECODE(last2, x2, y2, z2);
+
+                if (HASNEXT(last2)) {
+                    last1 = input[TOID(tx, y2, tz, size)];
+                    DECODE(last1, x1, y1, z1);
+                }
+            }
+
+            block[threadIdx.x][ty - y_end] = ENCODE(lasty, x2, z2, NOTSITE(last2), ZERO);
+        }
+
+        __syncthreads();
+
+        if(!threadIdx.y) {
+            int id = TOID(y_end + threadIdx.x, blockIdx.x * blockDim.x, tz, size);
+            for(int i = 0; i < blockDim.x; i++, id+=size) {
+                output[id] = block[i][threadIdx.x];
+            }
+        }
+
+        __syncthreads();
+    }
+}
+
+
+} // extern C
diff --git a/python/cucim/src/cucim/skimage/_shared/utils.py b/python/cucim/src/cucim/skimage/_shared/utils.py
index 6798628c7..1e014ef3a 100644
--- a/python/cucim/src/cucim/skimage/_shared/utils.py
+++ b/python/cucim/src/cucim/skimage/_shared/utils.py
@@ -95,10 +95,11 @@ def __call__(self, func):
         parameters = inspect.signature(func).parameters
         arg_idx = list(parameters.keys()).index(self.arg_name)
         warning_msg = (
-            f"{self.arg_name} argument is deprecated and will be removed "
-            f"in version {self.changed_version}. To avoid this warning, "
-            f"please do not use the {self.arg_name} argument. Please "
-            f"see {func.__name__} documentation for more details.")
+            f"{self.arg_name} argument is deprecated in upstream scikit-image "
+            f"and will be removed in cuCIM {self.changed_version}. To avoid "
+            f"this warning, please do not use the {self.arg_name} argument. "
+            f"Please see {func.__name__} documentation for more details."
+        )
 
         if self.help_msg is not None:
             warning_msg += f" {self.help_msg}"
@@ -206,7 +207,7 @@ def __init__(self, kwarg_mapping, deprecated_version, warning_msg=None,
             self.warning_msg = ("`{old_arg}` is a deprecated argument name "
                                 "for `{func_name}`. ")
             if removed_version is not None:
-                self.warning_msg += (f'It will be removed in '
+                self.warning_msg += (f'It will be removed in cuCIM '
                                      f'version {removed_version}.')
             self.warning_msg += "Please use `{new_arg}` instead."
         else:
@@ -248,10 +249,11 @@ class deprecate_multichannel_kwarg(deprecate_kwarg):
 
     """
 
-    def __init__(self, removed_version='1.0', multichannel_position=None):
+    def __init__(self, removed_version='2023.02.00',
+                 multichannel_position=None):
         super().__init__(
             kwarg_mapping={'multichannel': 'channel_axis'},
-            deprecated_version='0.19',
+            deprecated_version='22.02.00',
             warning_msg=None,
             removed_version=removed_version)
         self.position = multichannel_position
@@ -293,7 +295,7 @@ def fixed_func(*args, **kwargs):
 
         if func.__doc__ is not None:
             newdoc = docstring_add_deprecated(
-                func, {'multichannel': 'channel_axis'}, '0.19')
+                func, {'multichannel': 'channel_axis'}, '22.02.00')
             fixed_func.__doc__ = newdoc
         return fixed_func
 
diff --git a/python/cucim/src/cucim/skimage/_vendored/signaltools.py b/python/cucim/src/cucim/skimage/_vendored/signaltools.py
index c54415e99..3184055fa 100644
--- a/python/cucim/src/cucim/skimage/_vendored/signaltools.py
+++ b/python/cucim/src/cucim/skimage/_vendored/signaltools.py
@@ -424,7 +424,7 @@ def choose_conv_method(in1, in2, mode="full", measure=False):
     returns `direct` (e.g., to protect against floating point integer
     precision).
 
-    .. versionadded:: 0.19
+    .. versionadded:: 22.02.00
 
     Examples
     --------
diff --git a/python/cucim/src/cucim/skimage/filters/_median.py b/python/cucim/src/cucim/skimage/filters/_median.py
index d0da048c9..821cb0149 100644
--- a/python/cucim/src/cucim/skimage/filters/_median.py
+++ b/python/cucim/src/cucim/skimage/filters/_median.py
@@ -6,8 +6,9 @@
 from .._shared.utils import deprecate_kwarg
 
 
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def median(image, footprint=None, out=None, mode='nearest', cval=0.0,
            behavior='ndimage'):
     """Return local median of an image.
diff --git a/python/cucim/src/cucim/skimage/filters/thresholding.py b/python/cucim/src/cucim/skimage/filters/thresholding.py
index 8ca274209..9f26ec982 100644
--- a/python/cucim/src/cucim/skimage/filters/thresholding.py
+++ b/python/cucim/src/cucim/skimage/filters/thresholding.py
@@ -735,8 +735,8 @@ def threshold_li(image, *, tolerance=None, initial_guess=None,
     return threshold
 
 
-@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def threshold_minimum(image=None, nbins=256, max_num_iter=10000, *, hist=None):
     """Return threshold value based on minimum method.
 
diff --git a/python/cucim/src/cucim/skimage/measure/_label.py b/python/cucim/src/cucim/skimage/measure/_label.py
index 9dc2fc2bd..c9969cdff 100644
--- a/python/cucim/src/cucim/skimage/measure/_label.py
+++ b/python/cucim/src/cucim/skimage/measure/_label.py
@@ -16,8 +16,8 @@ def _get_structure(ndim, connectivity):
 
 # TODO: currently uses int32 for the labels. should add int64 option as well
 @deprecate_kwarg({'input': 'label_image'},
-                 deprecated_version='0.19',
-                 removed_version='1.0')
+                 deprecated_version='22.02.00',
+                 removed_version='23.02.00')
 def label(label_image, background=None, return_num=False, connectivity=None):
     r"""Label connected regions of an integer array.
 
diff --git a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
index 4a93622de..17ecd4a54 100644
--- a/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
+++ b/python/cucim/src/cucim/skimage/morphology/_skeletonize.py
@@ -39,8 +39,8 @@
                        0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=bool)
 
 
-@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def thin(image, max_num_iter=None):
     """
     Perform morphological thinning of a binary image.
@@ -62,7 +62,7 @@ def thin(image, max_num_iter=None):
 
     See Also
     --------
-    skeletonize, medial_axis
+    medial_axis
 
     Notes
     -----
diff --git a/python/cucim/src/cucim/skimage/morphology/binary.py b/python/cucim/src/cucim/skimage/morphology/binary.py
index 56c5b6654..705c577cb 100644
--- a/python/cucim/src/cucim/skimage/morphology/binary.py
+++ b/python/cucim/src/cucim/skimage/morphology/binary.py
@@ -37,8 +37,9 @@ def _iterate_binary_func(binary_func, image, footprint, out):
 # default with the same dimension as the input image and size 3 along each
 # axis.
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def binary_erosion(image, footprint=None, out=None):
     """Return fast binary morphological erosion of an image.
 
@@ -94,8 +95,9 @@ def binary_erosion(image, footprint=None, out=None):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def binary_dilation(image, footprint=None, out=None):
     """Return fast binary morphological dilation of an image.
 
@@ -149,8 +151,9 @@ def binary_dilation(image, footprint=None, out=None):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def binary_opening(image, footprint=None, out=None):
     """Return fast binary morphological opening of an image.
 
@@ -199,8 +202,9 @@ def binary_opening(image, footprint=None, out=None):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def binary_closing(image, footprint=None, out=None):
     """Return fast binary morphological closing of an image.
 
diff --git a/python/cucim/src/cucim/skimage/morphology/footprints.py b/python/cucim/src/cucim/skimage/morphology/footprints.py
index 744f11535..70e8fda7c 100644
--- a/python/cucim/src/cucim/skimage/morphology/footprints.py
+++ b/python/cucim/src/cucim/skimage/morphology/footprints.py
@@ -175,8 +175,8 @@ def _decompose_size(size, kernel_size=3):
 
 
 @deprecate_kwarg({'height': 'ncols', 'width': 'nrows'},
-                 deprecated_version='0.18.0',
-                 removed_version='0.20.0')
+                 deprecated_version='21.06.00',
+                 removed_version='22.02.00')
 def rectangle(nrows, ncols, dtype=cp.uint8, *, decomposition=None):
     """Generates a flat, rectangular-shaped footprint.
 
diff --git a/python/cucim/src/cucim/skimage/morphology/gray.py b/python/cucim/src/cucim/skimage/morphology/gray.py
index ca6fa2845..96bb4c2bd 100644
--- a/python/cucim/src/cucim/skimage/morphology/gray.py
+++ b/python/cucim/src/cucim/skimage/morphology/gray.py
@@ -157,8 +157,8 @@ def func_out(image, footprint, out=None, *args, **kwargs):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 def erosion(image, footprint=None, out=None, shift_x=False, shift_y=False):
     """Return grayscale morphological erosion of an image.
 
@@ -236,8 +236,8 @@ def erosion(image, footprint=None, out=None, shift_x=False, shift_y=False):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 def dilation(image, footprint=None, out=None, shift_x=False, shift_y=False):
     """Return grayscale morphological dilation of an image.
 
@@ -324,8 +324,8 @@ def dilation(image, footprint=None, out=None, shift_x=False, shift_y=False):
     return out
 
 
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 @default_footprint
 @pad_for_eccentric_footprints
 def opening(image, footprint=None, out=None):
@@ -390,8 +390,8 @@ def opening(image, footprint=None, out=None):
     return out
 
 
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 @default_footprint
 @pad_for_eccentric_footprints
 def closing(image, footprint=None, out=None):
@@ -474,8 +474,8 @@ def _white_tophat_seqence(image, footprints, out):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 def white_tophat(image, footprint=None, out=None):
     """Return white top hat of an image.
 
@@ -565,8 +565,8 @@ def white_tophat(image, footprint=None, out=None):
 
 
 @default_footprint
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 def black_tophat(image, footprint=None, out=None):
     """Return black top hat of an image.
 
diff --git a/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py b/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py
index 8847aaaf0..646503f13 100644
--- a/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py
+++ b/python/cucim/src/cucim/skimage/morphology/grayreconstruct.py
@@ -15,8 +15,8 @@
 from .._shared.utils import deprecate_kwarg
 
 
-@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg(kwarg_mapping={'selem': 'footprint'},
+                 removed_version="23.02.00", deprecated_version="22.02.00")
 def reconstruction(seed, mask, method='dilation', footprint=None, offset=None):
     """Perform a morphological reconstruction of an image.
 
diff --git a/python/cucim/src/cucim/skimage/morphology/misc.py b/python/cucim/src/cucim/skimage/morphology/misc.py
index 62ff9bfd0..d9b5304ea 100644
--- a/python/cucim/src/cucim/skimage/morphology/misc.py
+++ b/python/cucim/src/cucim/skimage/morphology/misc.py
@@ -73,7 +73,7 @@ def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False,
         labelling if `ar` is bool.
     in_place : bool, optional (default: False)
         If ``True``, remove the objects in the input array itself.
-        Otherwise, make a copy. Deprecated since version 0.19. Please
+        Otherwise, make a copy. Deprecated since version 22.02.00. Please
         use `out` instead.
     out : ndarray
         Array of the same shape as `ar`, into which the output is
@@ -172,7 +172,7 @@ def remove_small_holes(ar, area_threshold=64, connectivity=1, in_place=False,
         The connectivity defining the neighborhood of a pixel.
     in_place : bool, optional (default: False)
         If `True`, remove the connected components in the input array
-        itself. Otherwise, make a copy. Deprecated since version 0.19.
+        itself. Otherwise, make a copy. Deprecated since version 22.02.00.
         Please use `out` instead.
     out : ndarray
         Array of the same shape as `ar` and bool dtype, into which the
diff --git a/python/cucim/src/cucim/skimage/restoration/_denoise.py b/python/cucim/src/cucim/skimage/restoration/_denoise.py
index 614c5e7b3..392fc4ad2 100644
--- a/python/cucim/src/cucim/skimage/restoration/_denoise.py
+++ b/python/cucim/src/cucim/skimage/restoration/_denoise.py
@@ -91,8 +91,9 @@ def _denoise_tv_chambolle_nd(image, weight=0.1, eps=2.0e-4, max_num_iter=200):
     return out
 
 
-@utils.deprecate_kwarg({'n_iter_max': 'max_num_iter'}, removed_version="1.0",
-                       deprecated_version="0.19.2")
+@utils.deprecate_kwarg({'n_iter_max': 'max_num_iter'},
+                       removed_version="23.02.00",
+                       deprecated_version="22.06.00")
 @utils.deprecate_multichannel_kwarg(multichannel_position=4)
 def denoise_tv_chambolle(image, weight=0.1, eps=2.0e-4, max_num_iter=200,
                          multichannel=False, *, channel_axis=None):
diff --git a/python/cucim/src/cucim/skimage/restoration/deconvolution.py b/python/cucim/src/cucim/skimage/restoration/deconvolution.py
index 9809637eb..2a1cfa67c 100644
--- a/python/cucim/src/cucim/skimage/restoration/deconvolution.py
+++ b/python/cucim/src/cucim/skimage/restoration/deconvolution.py
@@ -383,8 +383,8 @@ 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="1.0",
-                 deprecated_version="0.19")
+@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/segmentation/morphsnakes.py b/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
index 9938cb9d1..4afad4ae5 100644
--- a/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
+++ b/python/cucim/src/cucim/skimage/segmentation/morphsnakes.py
@@ -219,8 +219,9 @@ def inverse_gaussian_gradient(image, alpha=100.0, sigma=5.0):
     return 1.0 / cp.sqrt(1.0 + alpha * gradnorm)
 
 
-@deprecate_kwarg({'iterations': 'num_iter'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg({'iterations': 'num_iter'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def morphological_chan_vese(image, num_iter, init_level_set='checkerboard',
                             smoothing=1, lambda1=1, lambda2=1,
                             iter_callback=lambda x: None):
@@ -326,8 +327,9 @@ def morphological_chan_vese(image, num_iter, init_level_set='checkerboard',
     return u
 
 
-@deprecate_kwarg({'iterations': 'num_iter'}, removed_version="1.0",
-                 deprecated_version="0.19")
+@deprecate_kwarg({'iterations': 'num_iter'},
+                 removed_version="23.02.00",
+                 deprecated_version="22.02.00")
 def morphological_geodesic_active_contour(gimage, num_iter,
                                           init_level_set='disk', smoothing=1,
                                           threshold='auto', balloon=0,
diff --git a/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py b/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
index 2acaaf8ba..8f2130be1 100644
--- a/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
+++ b/python/cucim/src/cucim/skimage/segmentation/random_walker_segmentation.py
@@ -343,12 +343,6 @@ def random_walker(data, labels, beta=130, mode='cg_j', tol=1.e-3, copy=True,
           probability that label `label_nb` reaches the pixel `(i, j)`
           first.
 
-    See Also
-    --------
-    skimage.morphology.watershed : watershed segmentation
-        A segmentation algorithm based on mathematical morphology
-        and "flooding" of regions from markers.
-
     Notes
     -----
     Multichannel inputs are scaled with all channel data combined. Ensure all
diff --git a/python/cucim/src/cucim/skimage/transform/_warps.py b/python/cucim/src/cucim/skimage/transform/_warps.py
index 7e3628ed5..3e8ae1af8 100644
--- a/python/cucim/src/cucim/skimage/transform/_warps.py
+++ b/python/cucim/src/cucim/skimage/transform/_warps.py
@@ -260,8 +260,8 @@ def rescale(image, scale, order=None, mode='reflect', cval=0, clip=True,
         Otherwise, this parameter indicates which axis of the array corresponds
         to channels.
 
-        .. versionadded:: 0.19
-           ``channel_axis`` was added in 0.19.
+        .. versionadded:: 22.02.00
+           ``channel_axis`` was added in 22.02.00.
 
     Notes
     -----
@@ -1094,8 +1094,8 @@ def warp_polar(image, center=None, *, radius=None, output_shape=None,
         Otherwise, this parameter indicates which axis of the array corresponds
         to channels.
 
-        .. versionadded:: 0.19
-           ``channel_axis`` was added in 0.19.
+        .. versionadded:: 22.02.00
+           ``channel_axis`` was added in 22.02.00.
     **kwargs : keyword arguments
         Passed to `transform.warp`.
 


Test patch:
diff --git a/python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py b/python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py
new file mode 100644
index 000000000..ab46ae1fb
--- /dev/null
+++ b/python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py
@@ -0,0 +1,145 @@
+from copy import copy
+
+import cupy as cp
+import numpy as np
+import pytest
+import scipy.ndimage as ndi_cpu
+
+from cucim.core.operations.morphology import distance_transform_edt
+
+
+def binary_image(shape, pct_true=50):
+    rng = cp.random.default_rng(123)
+    x = rng.integers(0, 100, size=shape, dtype=cp.uint8)
+    return x >= pct_true
+
+
+def assert_percentile_equal(arr1, arr2, pct=95):
+    """Assert that at least pct% of the entries in arr1 and arr2 are equal."""
+    pct_mismatch = (100 - pct) / 100
+    arr1 = cp.asnumpy(arr1)
+    arr2 = cp.asnumpy(arr2)
+    mismatch = np.sum(arr1 != arr2) / arr1.size
+    assert mismatch < pct_mismatch
+
+
+@pytest.mark.parametrize('return_indices', [False, True])
+@pytest.mark.parametrize('return_distances', [False, True])
+@pytest.mark.parametrize(
+    'shape, sampling',
+    [
+        ((256, 128), None),
+        ((384, 256), (1.5, 1.5)),
+        ((14, 32, 50), None),
+        ((50, 32, 24), (2, 2, 2)),
+    ]
+)
+@pytest.mark.parametrize('density', [5, 50, 95])
+@pytest.mark.parametrize('block_params', [None, (1, 1, 1)])
+def test_distance_transform_edt(
+    shape, sampling, return_distances, return_indices, density, block_params
+):
+
+    if not (return_indices or return_distances):
+        return
+
+    kwargs_scipy = dict(
+        sampling=sampling,
+        return_distances=return_distances,
+        return_indices=return_indices,
+    )
+    kwargs_cucim = copy(kwargs_scipy)
+    kwargs_cucim['block_params'] = block_params
+    img = binary_image(shape, pct_true=density)
+    out = distance_transform_edt(img, **kwargs_cucim)
+    expected = ndi_cpu.distance_transform_edt(cp.asnumpy(img), **kwargs_scipy)
+    if return_indices and return_distances:
+        assert len(out) == 2
+        cp.testing.assert_allclose(out[0], expected[0])
+        # May differ at a small % of coordinates where multiple points were
+        # equidistant.
+        assert_percentile_equal(out[1], expected[1], pct=95)
+    elif return_distances:
+        cp.testing.assert_allclose(out, expected)
+    elif return_indices:
+        assert_percentile_equal(out, expected, pct=95)
+
+
+@pytest.mark.parametrize('return_indices', [False, True])
+@pytest.mark.parametrize('return_distances', [False, True])
+@pytest.mark.parametrize(
+    'shape, sampling',
+    [
+        ((384, 256), (1, 3)),
+        ((50, 32, 24), (1, 2, 4)),
+    ]
+)
+@pytest.mark.parametrize('density', [5, 50, 95])
+def test_distance_transform_edt_nonuniform_sampling(
+    shape, sampling, return_distances, return_indices, density
+):
+
+    if not (return_indices or return_distances):
+        return
+
+    kwargs_scipy = dict(
+        sampling=sampling,
+        return_distances=return_distances,
+        return_indices=return_indices,
+    )
+    kwargs_cucim = copy(kwargs_scipy)
+    img = binary_image(shape, pct_true=density)
+    if sampling is not None and len(np.unique(sampling)) != 1:
+        with pytest.raises(NotImplementedError):
+            distance_transform_edt(img, **kwargs_cucim)
+        return
+
+
+@pytest.mark.parametrize('value', [0, 1, 3])
+@pytest.mark.parametrize('ndim', [2, 3])
+def test_distance_transform_edt_uniform_valued(value, ndim):
+    """ensure default block_params is robust to anisotropic shape."""
+    img = cp.full((48, ) * ndim, value, dtype=cp.uint8)
+    # ensure there is at least 1 pixel at background intensity
+    img[(slice(24, 25),) * ndim] = 0
+    out = distance_transform_edt(img)
+    expected = ndi_cpu.distance_transform_edt(cp.asnumpy(img))
+    cp.testing.assert_allclose(out, expected)
+
+
+@pytest.mark.parametrize('sx', list(range(16)))
+@pytest.mark.parametrize('sy', list(range(16)))
+def test_distance_transform_edt_2d_aniso(sx, sy):
+    """ensure default block_params is robust to anisotropic shape."""
+    shape = (128 + sy, 128 + sx)
+    img = binary_image(shape, pct_true=80)
+    out = distance_transform_edt(img)
+    expected = ndi_cpu.distance_transform_edt(cp.asnumpy(img))
+    cp.testing.assert_allclose(out, expected)
+
+
+@pytest.mark.parametrize('sx', list(range(4)))
+@pytest.mark.parametrize('sy', list(range(4)))
+@pytest.mark.parametrize('sz', list(range(4)))
+def test_distance_transform_edt_3d_aniso(sx, sy, sz):
+    """ensure default block_params is robust to anisotropic shape."""
+    shape = (16 + sz, 32 + sy, 48 + sx)
+    img = binary_image(shape, pct_true=80)
+    out = distance_transform_edt(img)
+    expected = ndi_cpu.distance_transform_edt(cp.asnumpy(img))
+    cp.testing.assert_allclose(out, expected)
+
+
+@pytest.mark.parametrize('ndim', [1, 4, 5])
+def test_distance_transform_edt_unsupported_ndim(ndim):
+    with pytest.raises(NotImplementedError):
+        distance_transform_edt(cp.zeros((8,) * ndim))
+
+
+@pytest.mark.skip(reason="excessive memory requirement")
+def test_distance_transform_edt_3d_int64():
+    shape = (1280, 1280, 1280)
+    img = binary_image(shape, pct_true=80)
+    distance_transform_edt(img)
+    # Note: no validation vs. scipy.ndimage due to excessive run time
+    return
diff --git a/python/cucim/src/cucim/skimage/_shared/tests/test_utils.py b/python/cucim/src/cucim/skimage/_shared/tests/test_utils.py
index f7135a0e3..c7d817b84 100644
--- a/python/cucim/src/cucim/skimage/_shared/tests/test_utils.py
+++ b/python/cucim/src/cucim/skimage/_shared/tests/test_utils.py
@@ -66,13 +66,13 @@ def bar(arg0, arg1=0, arg2=1):
 
 def test_deprecate_kwarg():
 
-    @deprecate_kwarg({'old_arg1': 'new_arg1'}, '0.19')
+    @deprecate_kwarg({'old_arg1': 'new_arg1'}, '22.02.00')
     def foo(arg0, new_arg1=1, arg2=None):
         """Expected docstring"""
         return arg0, new_arg1, arg2
 
     @deprecate_kwarg({'old_arg1': 'new_arg1'},
-                     deprecated_version='0.19',
+                     deprecated_version='22.02.00',
                      warning_msg="Custom warning message")
     def bar(arg0, new_arg1=1, arg2=None):
         """Expected docstring"""
@@ -116,7 +116,7 @@ def bar(arg0, new_arg1=1, arg2=None):
     old_arg1 : DEPRECATED
         Deprecated in favor of `new_arg1`.
 
-        .. deprecated:: 0.19
+        .. deprecated:: 22.02.00
 """
 
     assert len(recorded) == 0
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 1c53458f1..c0f6792ff 100644
--- a/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
+++ b/python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py
@@ -297,11 +297,7 @@ def test_rescale_in_range_clip():
 def test_rescale_out_range(dtype):
     """Check that output range is correct.
 
-    .. versionchanged:: 0.17
-        This function used to return dtype matching the input dtype. It now
-        matches the output.
-
-    .. versionchanged:: 0.19
+    .. versionchanged:: 22.02.00
         float16 and float32 inputs now result in float32 output. Formerly they
         would give float64 outputs.
     """
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 b6a697789..58f3bc0a1 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
@@ -25,7 +25,7 @@
 # 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':
+    if 'channel_axis' in kwargs and skimage.__version__ < '22.02.00':
         # 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
@@ -62,7 +62,7 @@ def test_invariant_denoise_color(dtype):
 
 def test_invariant_denoise_color_deprecated():
 
-    # can remove \A\Z here if only testing with scikit-image 0.19
+    # can remove \A\Z here if only testing with cuCIM 22.02.00
     with expected_warnings([r"`multichannel` is a deprecated argument|\A\Z"]):
         denoised_img_color = _invariant_denoise(
             noisy_img_color, _denoise_wavelet,