Incorporate bug fixes from skimage 0.19.3 (#312)

This PR incorporates the two bug fixes from skimage v0.19.3 (released on June 12th) that are relevant to cuCIM.

- bug fix to ensure all color channel share the same histogram bins in `cucim.skimage.exposure.histogram`
- clipping in warp functions should respect user-specified `cval` during clipping even if it is outside the original image range

There was also a fix to Canny edge detection, but we had already resolved that issue in cuCIM. I did go ahead and add the new test case for it to #310.

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

Approvers:
  - https://github.com/jakirkham

URL: #312

python/cucim/src/cucim/skimage/exposure/exposure.py

@@ -29,9 +29,7 @@ def _offset_array(arr, low_boundary, high_boundary):
             # prevent overflow errors when offsetting
             arr = arr.astype(offset_dtype)
         arr = arr - offset
-    else:
-        offset = 0
-    return arr, offset
+    return arr
 
 
 def _bincount_histogram_centers(image, source_range):
@@ -75,8 +73,10 @@ def _bincount_histogram(image, source_range, bin_centers=None):
     if bin_centers is None:
         bin_centers = _bincount_histogram_centers(image, source_range)
     image_min, image_max = bin_centers[0], bin_centers[-1]
-    image, offset = _offset_array(image, image_min.item(), image_max.item())  # synchronize  # noqa
-    hist = cp.bincount(image.ravel(), minlength=image_max - image_min + 1)
+    image = _offset_array(image, image_min.item(), image_max.item())  # synchronize  # noqa
+    hist = cp.bincount(
+        image.ravel(), minlength=image_max - min(image_min, 0) + 1
+    )
     if source_range == 'image':
         idx = max(image_min, 0)
         hist = hist[idx:]

python/cucim/src/cucim/skimage/exposure/tests/test_exposure.py

@@ -21,7 +21,9 @@
 def test_wrong_source_range():
     im = cp.array([-1, 100], dtype=cp.int8)
     with pytest.raises(ValueError):
-        frequencies, bin_centers = exposure.histogram(im, source_range="foobar")
+        frequencies, bin_centers = exposure.histogram(
+            im, source_range="foobar"
+        )
 
 
 def test_negative_overflow():
@@ -50,6 +52,15 @@ def test_int_range_image():
     assert bin_centers[-1] == 100
 
 
+def test_multichannel_int_range_image():
+    im = cp.array([[10, 5], [100, 102]], dtype=np.int8)
+    frequencies, bin_centers = exposure.histogram(im, channel_axis=-1)
+    for ch in range(im.shape[-1]):
+        assert len(frequencies[ch]) == len(bin_centers)
+    assert bin_centers[0] == 5
+    assert bin_centers[-1] == 102
+
+
 def test_peak_uint_range_dtype():
     im = cp.array([10, 100], dtype=cp.uint8)
     frequencies, bin_centers = exposure.histogram(im, source_range="dtype")

python/cucim/src/cucim/skimage/transform/_warps.py

@@ -766,19 +766,33 @@ def _clip_warp_output(input_image, output_image, mode, cval, clip):
     """
     if clip:
         min_val = input_image.min().item()
-        max_val = input_image.max().item()
-
-        preserve_cval = (mode == 'constant' and not
-                         (min_val <= cval <= max_val))
+        if np.isnan(min_val):
+            # NaNs detected, use NaN-safe min/max
+            min_func = cp.nanmin
+            max_func = cp.nanmax
+            min_val = min_func(input_image).item()
+        else:
+            min_func = cp.min
+            max_func = cp.max
+        max_val = max_func(input_image).item()
+
+        # Check if cval has been used such that it expands the effective input
+        # range
+        preserve_cval = (
+            mode == 'constant'
+            and not min_val <= cval <= max_val
+            and min_func(output_image) <= cval <= max_func(output_image)
+        )
 
+        # expand min/max range to account for cval
         if preserve_cval:
-            cval_mask = output_image == cval
+            # cast cval to the same dtype as the input image
+            cval = input_image.dtype.type(cval)
+            min_val = min(min_val, cval)
+            max_val = max(max_val, cval)
 
         cp.clip(output_image, min_val, max_val, out=output_image)
 
-        if preserve_cval:
-            output_image[cval_mask] = cval
-
 
 def warp(image, inverse_map, map_args={}, output_shape=None, order=None,
          mode='constant', cval=0., clip=True, preserve_range=False):

python/cucim/src/cucim/skimage/transform/tests/test_warps.py

@@ -118,6 +118,78 @@ def test_warp_clip():
     assert_almost_equal(float(outx.max()), 1)
 
 
+@pytest.mark.parametrize('order', [0, 1])
+def test_warp_clip_image_containing_nans(order):
+    # Test that clipping works as intended on an image with NaNs
+    # Orders >1 do not produce good output when the input image has
+    # NaNs, so those orders are not tested.
+
+    x = cp.ones((15, 15), dtype=cp.float64)
+    x[7, 7] = cp.nan
+
+    outx = rotate(x, 45, order=order, cval=2, resize=True, clip=True)
+
+    assert_almost_equal(cp.nanmin(outx).item(), 1)
+    assert_almost_equal(cp.nanmax(outx).item(), 2)
+
+
+@pytest.mark.parametrize('order', [0, 1])
+def test_warp_clip_cval_is_nan(order):
+    # Test that clipping works as intended when cval is NaN
+    # Orders > 1 do not produce good output when cval is NaN, so those
+    # orders are not tested.
+
+    x = cp.ones((15, 15), dtype=cp.float64)
+    x[5:-5, 5:-5] = 2
+
+    outx = rotate(x, 45, order=order, cval=cp.nan, resize=True, clip=True)
+
+    assert_almost_equal(cp.nanmin(outx).item(), 1)
+    assert_almost_equal(cp.nanmax(outx).item(), 2)
+
+
+@pytest.mark.parametrize('order', range(6))
+def test_warp_clip_cval_outside_input_range(order):
+    # Test that clipping behavior considers cval part of the input range
+
+    x = cp.ones((15, 15), dtype=cp.float64)
+
+    # Specify a cval that is outside the input range to check clipping
+    outx = rotate(x, 45, order=order, cval=2, resize=True, clip=True)
+
+    # The corners should be cval for all interpolation orders
+    outx = cp.asnumpy(outx)
+    assert_array_almost_equal([outx[0, 0], outx[0, -1],
+                               outx[-1, 0], outx[-1, -1]], 2)
+
+    # For all interpolation orders other than nearest-neighbor, the clipped
+    # output should have some pixels with values between the input (1) and
+    # cval (2) (i.e., clipping should not set them to 1)
+    if order > 0:
+        assert np.sum(np.less(1, outx) * np.less(outx, 2)) > 0
+
+
+@pytest.mark.parametrize('order', range(6))
+def test_warp_clip_cval_not_used(order):
+    # Test that clipping does not consider cval part of the input range if it
+    # is not used in the output image
+
+    x = cp.ones((15, 15), dtype=cp.float64)
+    x[5:-5, 5:-5] = 2
+
+    # Transform the image by stretching it out by one pixel on each side so
+    # that cval will not actually be used
+    scale = 15 / (15 + 2)
+    transform = AffineTransform(scale=scale, translation=(1, 1))
+    outx = warp(x, transform, mode='constant', order=order, cval=0, clip=True)
+
+    # At higher orders of interpolation, the transformed image has overshoots
+    # beyond the input range that should be clipped to the range 1 to 2.  Even
+    # though cval=0, the minimum value of the clipped output image should be
+    # 1 and not affected by the unused cval.
+    assert_array_almost_equal(outx.min(), 1)
+
+
 def test_homography():
     x = cp.zeros((5, 5), dtype=cp.double)
     x[1, 1] = 1