Add Euclidean distance transform for images/volumes (#318)

### Overview

closes #319

This PR adds an implementation of `scipy.ndimage.distance_transform_edt`. For each foreground pixel in a binary image, this function computes the minimal Euclidean distance to reach a background pixel. This is a SciPy function rather than a scikit-image one so I have put it under `cucim.core` instead of `cucim.skimage`. Longer term we could move this upstream to CuPy, but I think we want to have a copy here so that we can use it in the near term to implement some of the missing scikit-image functionality. This function is used by the following functions in the scikit-image API, so this PR will help us implement these in future PRs:

`skimage.morphology.medial_axis` (**update**: see #342)
`skimage.segmentation.expand_labels` (**update**: see #341)
`skimage.segmentation.chan_vese` (**update**: see #343)
It is also used in examples related to `skimage.segmentation.watershed` and is needed for implementation of ITK's `SignedMaurerDistanceMapImageFilter` for [itk-cucim](InsightSoftwareConsortium/itk_cucim#13).

The algorithm used here is adapted from the C++ code in the [PBA+ repository](https://github.com/orzzzjq/Parallel-Banding-Algorithm-plus) (MIT-licensed)

### extensions made to PBA+ kernel source
- larger sizes can be supported due to runtime code generation of the integer datatypes and encode/decode functions.

### current known limitations
- 2D and 3D only (this should cover the most common use cases)
- The `sampling` argument is not fully supported. This can likely be done with a bit of effort, but will require modifications to the CUDA kernels
- User-specified output arguments are not currently supported. We can potentially add this in the future
- The indices returned by `return_indices=True` are equally valid, but not always identical to those returned by SciPy. This is because there can be multiple indices with an identical distance, so which one gets chosen in that case is implementation-dependent.

### initial benchmarks relative to SciPy

Here `% true` is the percentage of the image that corresponds to foreground pixels. Even for fairly small images there is some benefit, with the benefit becoming two orders of magnitude at larger sizes.

shape      | % true |  cuCIM |  SciPy  | acceleration
-----------------|--------|---------|---------|--------------
(256, 256) | 5 | 0.00108 | 0.00353 | 3.25
(512, 512) | 5 | 0.00108 | 0.01552 | 14.42
(2048, 2048) | 5 | 0.00252 | 0.34434 | 136.86
(4096, 4096) | 5 | 0.00765 | 1.58948 | 207.88
(32, 32, 32) | 5 | 0.00103 | 0.00305 | 2.98
(128, 128, 128) | 5 | 0.00153 | 0.30103 | 196.26
(256, 256, 256) | 5 | 0.00763 | 3.17872 | 416.37
(384, 384, 384) | 5 | 0.02460 | 14.28779 | 580.89
(256, 256) | 50 | 0.00107 | 0.00430 | 4.01
(512, 512) | 50 | 0.00109 | 0.01878 | 17.30
(2048, 2048) | 50 | 0.00299 | 0.39304 | 131.60
(4096, 4096) | 50 | 0.00896 | 1.84686 | 206.19
(32, 32, 32) | 50 | 0.00102 | 0.00361 | 3.53
(128, 128, 128) | 50 | 0.00163 | 0.31657 | 194.66
(256, 256, 256) | 50 | 0.00914 | 3.35914 | 367.49
(384, 384, 384) | 50 | 0.03005 | 13.83219 | 460.30
(256, 256) | 95 | 0.00107 | 0.00344 | 3.22
(512, 512) | 95 | 0.00108 | 0.01638 | 15.17
(2048, 2048) | 95 | 0.00314 | 0.36996 | 117.81
(4096, 4096) | 95 | 0.00943 | 1.90475 | 202.05
(32, 32, 32) | 95 | 0.00102 | 0.00314 | 3.07
(128, 128, 128) | 95 | 0.00180 | 0.28843 | 159.80
(256, 256, 256) | 95 | 0.01073 | 3.23450 | 301.41
(384, 384, 384) | 95 | 0.03577 | 12.40526 | 346.82

### other comments
- can likely reduce memory overhead and improve performance a bit by refactoring some of the pre/post-processing code into elementwise kernels. (e.g. `encode3d`/`decode3d`, etc.)

- (JK) This may be able to leverage CuPy in the future ( cupy/cupy#6919 )

Authors:
  - Gregory Lee (https://github.com/grlee77)

Approvers:
  - Gigon Bae (https://github.com/gigony)
  - https://github.com/jakirkham

URL: #318

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.

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.

python/cucim/src/cucim/core/operations/morphology/__init__.py

@@ -0,0 +1,3 @@
+from ._distance_transform import distance_transform_edt
+
+__all__ = ["distance_transform_edt"]

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