Repository: rapidsai/cucim (language: python)

Issue:
[FEA] add cucim.skimage.segmentation.chan_vese
Once #318 is merged, we can add the `chan_vese` function for the Chan-Vese segmentation algorithm based on active contours. I will make a PR for this soon.

Add Euclidean distance transform for images/volumes
### 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](https://github.com/InsightSoftwareConsortium/itk_cucim/issues/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 ( https://github.com/cupy/cupy/issues/6919 )



Discussion / hints:

The failure here is due to use of `math.lcm` in a helper function. Apparently this function was only recently introduced in Python 3.9, so will need to introduce a workaround for that.
I am pretty confident that this is working as expected based on the tests here as well as use in other pending feature PRs. 

It would be good to review and merge this soon so I can submit a series of smaller new feature PRs that rely on this functionality (issues #334, #335, #336).

> I don't understand detailed algorithms so I feel I am not capable of reviewing this patch in detail.
> Overall, it looks good to me.

Thanks @gigony. I think the kernels are pretty well tested both here and in the original PBA+ repo, so probably doesn't need to much review in the details of the kernels in the `.h` files. I can hopefully revisit adding support for non-uniform values for `sampling` in a follow-up PR (for a future release).

This PR is a bit unique for cuCIM in that I kept the bulk of the kernel code in the form of CUDA C++ header files so that it is easier to view than if it was stored within a Python string. The file gets loaded into a string for purposes of runtime compilation, but keeping it in a separate file makes it easier to view with proper code highlighting either on GitHub or in an IDE.
Made an edit to the OP to xref the related upstream CuPy issue. Hope that is ok 🙂
@gpucibot merge

Fix patch:
diff --git a/python/cucim/src/cucim/skimage/segmentation/__init__.py b/python/cucim/src/cucim/skimage/segmentation/__init__.py
index 154467ee7..55b65058f 100644
--- a/python/cucim/src/cucim/skimage/segmentation/__init__.py
+++ b/python/cucim/src/cucim/skimage/segmentation/__init__.py
@@ -1,3 +1,4 @@
+from ._chan_vese import chan_vese
 from ._clear_border import clear_border
 from ._expand_labels import expand_labels
 from ._join import join_segmentations, relabel_sequential
@@ -15,6 +16,7 @@
     "clear_border",
     "join_segmentations",
     "relabel_sequential",
+    "chan_vese",
     "morphological_geodesic_active_contour",
     "morphological_chan_vese",
     "inverse_gaussian_gradient",
diff --git a/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
new file mode 100644
index 000000000..0401c1f07
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/segmentation/_chan_vese.py
@@ -0,0 +1,442 @@
+import cupy as cp
+import numpy as np
+from cupyx import rsqrt  # reciprocal sqrt
+
+from cucim.core.operations.morphology import distance_transform_edt
+
+from .._shared.utils import _supported_float_type, deprecate_kwarg
+
+
+@cp.fuse()
+def _fused_curvature(phi, x_start, x_end, y_start, y_end, ul, ur, ll, lr):
+    fy = (y_end - y_start) / 2.0
+    fx = (x_end - x_start) / 2.0
+    fyy = y_end + y_start - 2 * phi
+    fxx = x_end + x_start - 2 * phi
+    fxy = .25 * (lr + ul - ur - ll)
+    grad2 = fx**2 + fy**2
+    K = (fxx * fy**2 - 2 * fxy * fx * fy + fyy * fx**2)
+    K /= (grad2 * cp.sqrt(grad2) + 1e-8)
+    return K
+
+
+def _cv_curvature(phi):
+    """Returns the 'curvature' of a level set 'phi'.
+    """
+    P = cp.pad(phi, 1, mode='edge')
+    y_start = P[:-2, 1:-1]
+    y_end = P[2:, 1:-1]
+    x_start = P[1:-1, :-2]
+    x_end = P[1:-1, 2:]
+
+    lower_right = P[2:, 2:]
+    lower_left = P[2:, :-2]
+    upper_right = P[:-2, 2:]
+    upper_left = P[:-2, :-2]
+    K = _fused_curvature(phi, x_start, x_end, y_start, y_end, upper_left,
+                         upper_right, lower_left, lower_right)
+    return K
+
+
+@cp.fuse()
+def _fused_variance_kernel1(eta, x_start, x_mid, x_end, y_start, y_mid, y_end):
+    phixp = x_end - x_mid
+    phixn = x_mid - x_start
+    phix0 = x_end - x_start
+    phix0 /= 2.0
+    phixp *= phixp
+    phixn *= phixn
+    phix0 *= phix0
+
+    phiyp = y_end - y_mid
+    phiyn = y_mid - y_start
+    phiy0 = y_end - y_start
+    phiy0 /= 2.0
+    phiyp *= phiyp
+    phiyn *= phiyn
+    phiy0 *= phiy0
+
+    C1 = rsqrt(eta + phixp + phiy0)
+    C2 = rsqrt(eta + phixn + phiy0)
+    C3 = rsqrt(eta + phix0 + phiyp)
+    C4 = rsqrt(eta + phix0 + phiyn)
+
+    K = x_end * C1
+    K += x_start * C2
+    K += y_end * C3
+    K += y_start * C4
+
+    Csum = C1
+    Csum += C2
+    Csum += C3
+    Csum += C4
+
+    return K, Csum
+
+
+@cp.fuse()
+def _fused_hphi_hinv(phi):
+    Hphi = (phi > 0).astype(phi.dtype)
+    Hinv = 1.0 - Hphi
+    return Hphi, Hinv
+
+
+@cp.fuse()
+def _fused_variance_kernel2(
+    image, c1, c2, lam1, lam2, phi, K, dt, mu, delta_phi, Csum
+):
+    difference_term = image - c1
+    difference_term *= difference_term
+    difference_term *= -lam1
+
+    term2 = image - c2
+    term2 *= term2
+    term2 *= lam2
+    difference_term += term2
+
+    new_phi = phi + (dt * delta_phi) * (mu * K + difference_term)
+    out = new_phi / (1 + mu * dt * delta_phi * Csum)
+    return out
+
+
+def _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt):
+    """Returns the variation of level set 'phi' based on algorithm parameters.
+    """
+    eta = 1e-16
+    P = cp.pad(phi, 1, mode='edge')
+
+    x_end = P[1:-1, 2:]
+    x_mid = P[1:-1, 1:-1]
+    x_start = P[1:-1, :-2]
+
+    y_end = P[2:, 1:-1]
+    y_mid = P[1:-1, 1:-1]
+    y_start = P[:-2, 1:-1]
+
+    K, Csum = _fused_variance_kernel1(
+        eta, x_start, x_mid, x_end, y_start, y_mid, y_end
+    )
+    Hphi, Hinv = _fused_hphi_hinv(phi)
+    c1, c2 = _cv_calculate_averages(image, Hphi, Hinv)
+    delta_phi = _cv_delta(phi)
+    out = _fused_variance_kernel2(
+        image, c1, c2, lambda1, lambda2, phi, K, dt, mu, delta_phi, Csum
+    )
+    return out
+
+
+@cp.fuse()
+def _cv_heavyside(x, eps=1.):
+    """Returns the result of a regularised heavyside function of the
+    input value(s).
+    """
+    return 0.5 * (1. + (2. / cp.pi) * cp.arctan(x / eps))
+
+
+@cp.fuse()
+def _cv_delta(x, eps=1.):
+    """Returns the result of a regularised dirac function of the
+    input value(s).
+    """
+    return eps / (eps * eps + x * x)
+
+
+@cp.fuse()
+def _fused_inplace_eps_div(num, denom, eps):
+    denom += eps
+    num /= denom
+    return
+
+
+def _cv_calculate_averages(image, H, Hinv):
+    """Returns the average values 'inside' and 'outside'.
+    """
+    Hsum = cp.sum(H)
+    Hinvsum = cp.sum(Hinv)
+    avg_inside = cp.sum(image * H)
+    avg_oustide = cp.sum(image * Hinv)
+
+    eps = 10 * cp.finfo(image.dtype).eps
+    _fused_inplace_eps_div(avg_inside, Hsum, eps)
+    _fused_inplace_eps_div(avg_oustide, Hinvsum, eps)
+    return (avg_inside, avg_oustide)
+
+
+@cp.fuse()
+def _fused_difference_op1(image, c, h, lam):
+    out = image - c
+    out *= out
+    out *= h
+    out *= lam
+    return out
+
+
+def _cv_difference_from_average_term(image, Hphi, lambda_pos, lambda_neg):
+    """Returns the 'energy' contribution due to the difference from
+    the average value within a region at each point.
+    """
+    Hinv = 1. - Hphi
+    (c1, c2) = _cv_calculate_averages(image, Hphi, Hinv)
+    out = _fused_difference_op1(image, c1, Hphi, lambda_pos)
+    out += _fused_difference_op1(image, c2, Hinv, lambda_neg)
+    return out
+
+
+def _cv_edge_length_term(phi, mu):
+    """Returns the 'energy' contribution due to the length of the
+    edge between regions at each point, multiplied by a factor 'mu'.
+    """
+    e = _cv_curvature(phi)
+    e *= mu
+    return e
+
+
+def _cv_energy(image, phi, mu, lambda1, lambda2):
+    """Returns the total 'energy' of the current level set function.
+    """
+    H = _cv_heavyside(phi)
+    avgenergy = _cv_difference_from_average_term(image, H, lambda1, lambda2)
+    lenenergy = _cv_edge_length_term(phi, mu)
+    return cp.sum(avgenergy) + cp.sum(lenenergy)
+
+
+def _cv_checkerboard(image_size, square_size, dtype=cp.float64):
+    """Generates a checkerboard level set function.
+
+    According to Pascal Getreuer, such a level set function has fast
+    convergence.
+    """
+    yv = cp.arange(image_size[0], dtype=dtype)[:, np.newaxis]
+    xv = cp.arange(image_size[1], dtype=dtype)[np.newaxis, :]
+    sf = cp.pi / square_size
+    xv *= sf
+    yv *= sf
+    cp.sin(xv, out=xv)
+    cp.sin(yv, out=yv)
+    return xv * yv
+
+
+def _cv_large_disk(image_size):
+    """Generates a disk level set function.
+
+    The disk covers the whole image along its smallest dimension.
+    """
+    res = cp.ones(image_size, dtype=bool)
+    centerY = int((image_size[0] - 1) / 2)
+    centerX = int((image_size[1] - 1) / 2)
+    res[centerY, centerX] = 0.
+    radius = float(min(centerX, centerY))
+    out = radius - distance_transform_edt(res)
+    out /= radius
+    return out
+
+
+def _cv_small_disk(image_size):
+    """Generates a disk level set function.
+
+    The disk covers half of the image along its smallest dimension.
+    """
+    res = cp.ones(image_size, dtype=bool)
+    centerY = int((image_size[0] - 1) / 2)
+    centerX = int((image_size[1] - 1) / 2)
+    res[centerY, centerX] = 0.
+    radius = float(min(centerX, centerY)) / 2.0
+    out = radius - distance_transform_edt(res)
+    out /= radius * 3
+    return out
+
+
+def _cv_init_level_set(init_level_set, image_shape, dtype=cp.float64):
+    """Generates an initial level set function conditional on input arguments.
+    """
+    if type(init_level_set) == str:
+        if init_level_set == 'checkerboard':
+            res = _cv_checkerboard(image_shape, 5, dtype)
+        elif init_level_set == 'disk':
+            res = _cv_large_disk(image_shape)
+        elif init_level_set == 'small disk':
+            res = _cv_small_disk(image_shape)
+        else:
+            raise ValueError("Incorrect name for starting level set preset.")
+    else:
+        res = init_level_set
+    return res.astype(dtype, copy=False)
+
+
+@deprecate_kwarg({'max_iter': 'max_num_iter'}, removed_version="1.0",
+                 deprecated_version="0.19")
+def chan_vese(image, mu=0.25, lambda1=1.0, lambda2=1.0, tol=1e-3,
+              max_num_iter=500, dt=0.5, init_level_set='checkerboard',
+              extended_output=False):
+    """Chan-Vese segmentation algorithm.
+
+    Active contour model by evolving a level set. Can be used to
+    segment objects without clearly defined boundaries.
+
+    Parameters
+    ----------
+    image : (M, N) ndarray
+        Grayscale image to be segmented.
+    mu : float, optional
+        'edge length' weight parameter. Higher `mu` values will
+        produce a 'round' edge, while values closer to zero will
+        detect smaller objects.
+    lambda1 : float, optional
+        'difference from average' weight parameter for the output
+        region with value 'True'. If it is lower than `lambda2`, this
+        region will have a larger range of values than the other.
+    lambda2 : float, optional
+        'difference from average' weight parameter for the output
+        region with value 'False'. If it is lower than `lambda1`, this
+        region will have a larger range of values than the other.
+    tol : float, positive, optional
+        Level set variation tolerance between iterations. If the
+        L2 norm difference between the level sets of successive
+        iterations normalized by the area of the image is below this
+        value, the algorithm will assume that the solution was
+        reached.
+    max_num_iter : uint, optional
+        Maximum number of iterations allowed before the algorithm
+        interrupts itself.
+    dt : float, optional
+        A multiplication factor applied at calculations for each step,
+        serves to accelerate the algorithm. While higher values may
+        speed up the algorithm, they may also lead to convergence
+        problems.
+    init_level_set : str or (M, N) ndarray, optional
+        Defines the starting level set used by the algorithm.
+        If a string is inputted, a level set that matches the image
+        size will automatically be generated. Alternatively, it is
+        possible to define a custom level set, which should be an
+        array of float values, with the same shape as 'image'.
+        Accepted string values are as follows.
+
+        'checkerboard'
+            the starting level set is defined as
+            sin(x/5*pi)*sin(y/5*pi), where x and y are pixel
+            coordinates. This level set has fast convergence, but may
+            fail to detect implicit edges.
+        'disk'
+            the starting level set is defined as the opposite
+            of the distance from the center of the image minus half of
+            the minimum value between image width and image height.
+            This is somewhat slower, but is more likely to properly
+            detect implicit edges.
+        'small disk'
+            the starting level set is defined as the
+            opposite of the distance from the center of the image
+            minus a quarter of the minimum value between image width
+            and image height.
+    extended_output : bool, optional
+        If set to True, the return value will be a tuple containing
+        the three return values (see below). If set to False which
+        is the default value, only the 'segmentation' array will be
+        returned.
+
+    Returns
+    -------
+    segmentation : (M, N) ndarray, bool
+        Segmentation produced by the algorithm.
+    phi : (M, N) ndarray of floats
+        Final level set computed by the algorithm.
+    energies : list of floats
+        Shows the evolution of the 'energy' for each step of the
+        algorithm. This should allow to check whether the algorithm
+        converged.
+
+    Notes
+    -----
+    The Chan-Vese Algorithm is designed to segment objects without
+    clearly defined boundaries. This algorithm is based on level sets
+    that are evolved iteratively to minimize an energy, which is
+    defined by weighted values corresponding to the sum of differences
+    intensity from the average value outside the segmented region, the
+    sum of differences from the average value inside the segmented
+    region, and a term which is dependent on the length of the
+    boundary of the segmented region.
+
+    This algorithm was first proposed by Tony Chan and Luminita Vese,
+    in a publication entitled "An Active Contour Model Without Edges"
+    [1]_.
+
+    This implementation of the algorithm is somewhat simplified in the
+    sense that the area factor 'nu' described in the original paper is
+    not implemented, and is only suitable for grayscale images.
+
+    Typical values for `lambda1` and `lambda2` are 1. If the
+    'background' is very different from the segmented object in terms
+    of distribution (for example, a uniform black image with figures
+    of varying intensity), then these values should be different from
+    each other.
+
+    Typical values for mu are between 0 and 1, though higher values
+    can be used when dealing with shapes with very ill-defined
+    contours.
+
+    The 'energy' which this algorithm tries to minimize is defined
+    as the sum of the differences from the average within the region
+    squared and weighed by the 'lambda' factors to which is added the
+    length of the contour multiplied by the 'mu' factor.
+
+    Supports 2D grayscale images only, and does not implement the area
+    term described in the original article.
+
+    References
+    ----------
+    .. [1] An Active Contour Model without Edges, Tony Chan and
+           Luminita Vese, Scale-Space Theories in Computer Vision,
+           1999, :DOI:`10.1007/3-540-48236-9_13`
+    .. [2] Chan-Vese Segmentation, Pascal Getreuer Image Processing On
+           Line, 2 (2012), pp. 214-224,
+           :DOI:`10.5201/ipol.2012.g-cv`
+    .. [3] The Chan-Vese Algorithm - Project Report, Rami Cohen, 2011
+           :arXiv:`1107.2782`
+    """
+    if len(image.shape) != 2:
+        raise ValueError("Input image should be a 2D array.")
+
+    float_dtype = _supported_float_type(image.dtype)
+    phi = _cv_init_level_set(init_level_set, image.shape, dtype=float_dtype)
+    if type(phi) != cp.ndarray or phi.shape != image.shape:
+        raise ValueError("The dimensions of initial level set do not "
+                         "match the dimensions of image.")
+
+    image = image.astype(float_dtype, copy=False)
+    image = image - cp.min(image)
+    if cp.max(image) != 0:
+        image = image / cp.max(image)
+
+    i = 0
+    if extended_output:
+        old_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
+        energies = []
+    phivar = tol + 1
+
+    while phivar > tol and i < max_num_iter:
+        # Save old level set values
+        oldphi = phi
+
+        # Calculate new level set
+        phi = _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt)
+        phivar = phi - oldphi
+        phivar *= phivar
+        phivar = cp.sqrt(phivar.mean())
+
+        if extended_output:
+            # Extract energy
+            new_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
+
+            # Could compare energy to the previous level set to see if
+            # continuing is necessary
+
+            # Save old energy values
+            energies.append(old_energy)
+            old_energy = new_energy
+        i += 1
+
+    segmentation = phi > 0
+
+    if extended_output:
+        return (segmentation, phi, energies)
+    else:
+        return segmentation


Test patch:
diff --git a/python/cucim/src/cucim/skimage/segmentation/tests/test_chan_vese.py b/python/cucim/src/cucim/skimage/segmentation/tests/test_chan_vese.py
new file mode 100644
index 000000000..15b991883
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/segmentation/tests/test_chan_vese.py
@@ -0,0 +1,105 @@
+import cupy as cp
+import pytest
+from cupy.testing import assert_array_equal
+
+from cucim.skimage._shared._warnings import expected_warnings
+from cucim.skimage._shared.utils import _supported_float_type
+from cucim.skimage.segmentation import chan_vese
+
+
+@pytest.mark.parametrize('dtype', [cp.float32, cp.float64])
+def test_chan_vese_flat_level_set(dtype):
+    # because the algorithm evolves the level set around the
+    # zero-level, it the level-set has no zero level, the algorithm
+    # will not produce results in theory. However, since a continuous
+    # approximation of the delta function is used, the algorithm
+    # still affects the entirety of the level-set. Therefore with
+    # infinite time, the segmentation will still converge.
+    img = cp.zeros((10, 10), dtype=dtype)
+    img[3:6, 3:6] = 1
+    ls = cp.full((10, 10), 1000, dtype=dtype)
+    result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set=ls)
+    assert_array_equal(result.astype(float), cp.ones((10, 10)))
+    result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set=-ls)
+    assert_array_equal(result.astype(float), cp.zeros((10, 10)))
+
+
+def test_chan_vese_small_disk_level_set():
+    img = cp.zeros((10, 10))
+    img[3:6, 3:6] = 1
+    result = chan_vese(img, mu=0.0, tol=1e-3, init_level_set="small disk")
+    assert_array_equal(result.astype(float), img)
+
+
+def test_chan_vese_simple_shape():
+    img = cp.zeros((10, 10))
+    img[3:6, 3:6] = 1
+    result = chan_vese(img, mu=0.0, tol=1e-8).astype(float)
+    assert_array_equal(result, img)
+
+
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+def test_chan_vese_extended_output(dtype):
+    img = cp.zeros((10, 10), dtype=dtype)
+    img[3:6, 3:6] = 1
+    result = chan_vese(img, mu=0.0, tol=1e-8, extended_output=True)
+    float_dtype = _supported_float_type(dtype)
+    assert result[1].dtype == float_dtype
+    assert all(arr.dtype == float_dtype for arr in result[2])
+    assert_array_equal(len(result), 3)
+
+
+def test_chan_vese_remove_noise():
+    ref = cp.zeros((10, 10))
+    ref[1:6, 1:6] = cp.array([[0, 1, 1, 1, 0],
+                              [1, 1, 1, 1, 1],
+                              [1, 1, 1, 1, 1],
+                              [1, 1, 1, 1, 1],
+                              [0, 1, 1, 1, 0]])
+    img = ref.copy()
+    img[8, 3] = 1
+    result = chan_vese(img, mu=0.3, tol=1e-3, max_num_iter=100, dt=10,
+                       init_level_set="disk").astype(float)
+    assert_array_equal(result, ref)
+
+
+def test_chan_vese_incorrect_image_type():
+    img = cp.zeros((10, 10, 3))
+    ls = cp.zeros((10, 9))
+    with pytest.raises(ValueError):
+        chan_vese(img, mu=0.0, init_level_set=ls)
+
+
+def test_chan_vese_gap_closing():
+    ref = cp.zeros((20, 20))
+    ref[8:15, :] = cp.ones((7, 20))
+    img = ref.copy()
+    img[:, 6] = cp.zeros((20))
+    result = chan_vese(img, mu=0.7, tol=1e-3, max_num_iter=1000, dt=1000,
+                       init_level_set="disk").astype(float)
+    assert_array_equal(result, ref)
+
+
+def test_chan_vese_max_iter_deprecation():
+    img = cp.zeros((20, 20))
+    with expected_warnings(["`max_iter` is a deprecated argument"]):
+        chan_vese(img, max_iter=10)
+
+
+def test_chan_vese_incorrect_level_set():
+    img = cp.zeros((10, 10))
+    ls = cp.zeros((10, 9))
+    with pytest.raises(ValueError):
+        chan_vese(img, mu=0.0, init_level_set=ls)
+    with pytest.raises(ValueError):
+        chan_vese(img, mu=0.0, init_level_set="a")
+
+
+def test_chan_vese_blank_image():
+    img = cp.zeros((10, 10))
+    level_set = cp.random.rand(10, 10)
+    ref = level_set > 0
+    result = chan_vese(img, mu=0.0, tol=0.0, init_level_set=level_set)
+    assert_array_equal(result, ref)