additional minor updates (skimage 0.20)

python/cucim/src/cucim/skimage/color/colorconv.py

@@ -1259,7 +1259,7 @@ def lab2xyz(lab, illuminant="D65", observer="2", *, channel_axis=-1):
     nwarn = int(cp.count_nonzero(warnings))
     if nwarn > 0:  # synchronize!
         warn(f'Color data out of range: Z < 0 in {nwarn} pixels',
-             stacklevel=2)
+             stacklevel=3)
     return xyz
 
 

python/cucim/src/cucim/skimage/feature/corner.py

@@ -171,10 +171,11 @@ def _hessian_matrix_with_gaussian(image, sigma=1, mode='reflect', cval=0,
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
     order : {'rc', 'xy'}, optional
-        This parameter allows for the use of reverse or forward order of
-        the image axes in gradient computation. 'rc' indicates the use of
-        the first axis initially (Hrr, Hrc, Hcc), whilst 'xy' indicates the
-        usage of the last axis initially (Hxx, Hxy, Hyy)
+        NOTE: 'xy' is only an option for 2D images, higher dimensions must
+        always use 'rc' order. This parameter allows for the use of reverse or
+        forward order of the image axes in gradient computation. 'rc' indicates
+        the use of the first axis initially (Hrr, Hrc, Hcc), whilst 'xy'
+        indicates the usage of the last axis initially (Hxx, Hxy, Hyy).
 
     Returns
     -------
@@ -187,6 +188,10 @@ def _hessian_matrix_with_gaussian(image, sigma=1, mode='reflect', cval=0,
     image = img_as_float(image)
     float_dtype = _supported_float_type(image.dtype)
     image = image.astype(float_dtype, copy=False)
+    if image.ndim > 2 and order == "xy":
+        raise ValueError("Order 'xy' is only supported for 2D images.")
+    if order not in ["rc", "xy"]:
+        raise ValueError(f"unrecognized order: {order}")
 
     if np.isscalar(sigma):
         sigma = (sigma,) * image.ndim
@@ -223,7 +228,7 @@ def _hessian_matrix_with_gaussian(image, sigma=1, mode='reflect', cval=0,
 
     # 2.) apply the derivative along another axis as well
     axes = range(ndim)
-    if order == 'rc':
+    if order == 'xy':
         axes = reversed(axes)
     H_elems = [gaussian_(gradients[ax0], order=orders[ax1])
                for ax0, ax1 in combinations_with_replacement(axes, 2)]
@@ -257,10 +262,11 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc',
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
     order : {'rc', 'xy'}, optional
-        This parameter allows for the use of reverse or forward order of
-        the image axes in gradient computation. 'rc' indicates the use of
-        the first axis initially (Hrr, Hrc, Hcc), whilst 'xy' indicates the
-        usage of the last axis initially (Hxx, Hxy, Hyy)
+        NOTE: 'xy' is only an option for 2D images, higher dimensions must
+        always use 'rc' order. This parameter allows for the use of reverse or
+        forward order of the image axes in gradient computation. 'rc' indicates
+        the use of the first axis initially (Hrr, Hrc, Hcc), whilst 'xy'
+        indicates the usage of the last axis initially (Hxx, Hxy, Hyy).
     use_gaussian_derivatives : boolean, optional
         Indicates whether the Hessian is computed by convolving with Gaussian
         derivatives, or by a simple finite-difference operation.
@@ -310,6 +316,10 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc',
     image = img_as_float(image)
     float_dtype = _supported_float_type(image.dtype)
     image = image.astype(float_dtype, copy=False)
+    if image.ndim > 2 and order == "xy":
+        raise ValueError("Order 'xy' is only supported for 2D images.")
+    if order not in ["rc", "xy"]:
+        raise ValueError(f"unrecognized order: {order}")
 
     if use_gaussian_derivatives is None:
         use_gaussian_derivatives = False
@@ -333,7 +343,7 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc',
     gradients = gradient(gaussian_filtered)
     axes = range(image.ndim)
 
-    if order == "rc":
+    if order == "xy":
         axes = reversed(axes)
 
     H_elems = [
@@ -344,6 +354,31 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0, order='rc',
     return H_elems
 
 
+@cp.memoize()
+def _get_real_symmetric_2x2_det_kernel():
+    return cp.ElementwiseKernel(
+        in_params="F M00, F M01, F M11",
+        out_params="F det",
+        operation="det = M00 * M11 - M01 * M01;",
+        name="cucim_skimage_symmetric_det22_kernel")
+
+
+@cp.memoize()
+def _get_real_symmetric_3x3_det_kernel():
+
+    operation = """
+    det = M00 * (M11 * M22 - M12 * M12);
+    det -= M01 * (M01 * M22 - M12 * M02);
+    det += M02 * (M01 * M12 - M11 * M02);
+    """
+
+    return cp.ElementwiseKernel(
+        in_params="F M00, F M01, F M02, F M11, F M12, F M22",
+        out_params="F det",
+        operation=operation,
+        name="cucim_skimage_symmetric_det33_kernel")
+
+
 def hessian_matrix_det(image, sigma=1, approximate=True):
     """Compute the approximate Hessian Determinant over an image.
 
@@ -385,19 +420,31 @@ def hessian_matrix_det(image, sigma=1, approximate=True):
     image = image.astype(float_dtype, copy=False)
     if image.ndim == 2 and approximate:
         integral = integral_image(image)
-        return cp.asarray(_hessian_matrix_det(integral, sigma))
+        # integral image will promote to float64 for accuracy, but we can
+        # cast back to float_dtype here.
+        integral = integral.astype(float_dtype, copy=False)
+        return _hessian_matrix_det(integral, sigma)
     else:  # slower brute-force implementation for nD images
         if image.ndim in [2, 3]:
-            # Compute determinant as the product of the eigenvalues.
-            # This avoids the huge memory overhead of forming
-            # `_symmetric_image` as in the code below.
-            # Could optimize further by computing the determinant directly
-            # using ElementwiseKernels rather than reusing the eigenvalue ones.
+            det = cp.empty_like(image)
+            if image.ndim == 2:
+                kernel = _get_real_symmetric_2x2_det_kernel()
+            else:
+                kernel = _get_real_symmetric_3x3_det_kernel()
             H = hessian_matrix(image, sigma)
-            evs = hessian_matrix_eigvals(H)
-            return cp.prod(evs, axis=0)
-        hessian_mat_array = _symmetric_image(hessian_matrix(image, sigma))
-        return cp.linalg.det(hessian_mat_array)
+            kernel(*H, det)
+        else:
+            # # Compute determinant as the product of the eigenvalues.
+            # # This avoids the huge memory overhead of forming
+            # # `_symmetric_image` as in the code below.
+            # # Could optimize further by computing the determinant directly
+            # # using ElementwiseKernels rather than reusing the eigenvalue ones.
+            # H = hessian_matrix(image, sigma)
+            # evs = hessian_matrix_eigvals(H)
+            # return cp.prod(evs, axis=0)
+            hessian_mat_array = _symmetric_image(hessian_matrix(image, sigma))
+            det = cp.linalg.det(hessian_mat_array)
+        return det
 
 
 @cp.memoize()

python/cucim/src/cucim/skimage/feature/tests/test_corner.py

@@ -4,6 +4,7 @@
 from cupy.testing import assert_array_almost_equal, assert_array_equal
 from numpy.testing import assert_equal
 from skimage import data, draw
+from skimage.feature import hessian_matrix_det as hessian_matrix_det_cpu
 
 from cucim.skimage import img_as_float
 from cucim.skimage._shared._warnings import expected_warnings
@@ -129,23 +130,23 @@ def test_hessian_matrix(dtype):
     out_dtype = _supported_float_type(dtype)
     assert all(a.dtype == out_dtype for a in (Hrr, Hrc, Hcc))
     # fmt: off
-    assert_array_almost_equal(Hrr, cp.asarray([[0, 0,  0, 0, 0],    # noqa
+    assert_array_almost_equal(Hrr, cp.asarray([[0, 0,  2, 0, 0],    # noqa
                                                [0, 0,  0, 0, 0],    # noqa
-                                               [2, 0, -2, 0, 2],    # noqa
+                                               [0, 0, -2, 0, 0],    # noqa
                                                [0, 0,  0, 0, 0],    # noqa
-                                               [0, 0,  0, 0, 0]]))  # noqa
+                                               [0, 0,  2, 0, 0]]))  # noqa
 
     assert_array_almost_equal(Hrc, cp.asarray([[0,  0, 0,  0, 0],    # noqa
                                                [0,  1, 0, -1, 0],    # noqa
                                                [0,  0, 0,  0, 0],    # noqa
                                                [0, -1, 0,  1, 0],    # noqa
                                                [0,  0, 0,  0, 0]]))  # noqa
 
-    assert_array_almost_equal(Hcc, cp.asarray([[0, 0,  2, 0, 0],    # noqa
+    assert_array_almost_equal(Hcc, cp.asarray([[0, 0,  0, 0, 0],    # noqa
                                                [0, 0,  0, 0, 0],    # noqa
-                                               [0, 0, -2, 0, 0],    # noqa
+                                               [2, 0, -2, 0, 2],    # noqa
                                                [0, 0,  0, 0, 0],    # noqa
-                                               [0, 0,  2, 0, 0]]))  # noqa
+                                               [0, 0,  0, 0, 0]]))  # noqa
     # fmt: on
 
     with expected_warnings(["use_gaussian_derivatives currently defaults"]):
@@ -154,6 +155,25 @@ def test_hessian_matrix(dtype):
         hessian_matrix(square, sigma=0.1, order="rc")
 
 
+@pytest.mark.parametrize('use_gaussian_derivatives', [False, True])
+def test_hessian_matrix_order(use_gaussian_derivatives):
+    square = cp.zeros((5, 5), dtype=float)
+    square[2, 2] = 4
+
+    Hxx, Hxy, Hyy = hessian_matrix(
+        square, sigma=0.1, order="xy",
+        use_gaussian_derivatives=use_gaussian_derivatives)
+
+    Hrr, Hrc, Hcc = hessian_matrix(
+        square, sigma=0.1, order="rc",
+        use_gaussian_derivatives=use_gaussian_derivatives)
+
+    # verify results are equivalent, just reversed in order
+    cp.testing.assert_allclose(Hxx, Hcc, atol=1e-30)
+    cp.testing.assert_allclose(Hxy, Hrc, atol=1e-30)
+    cp.testing.assert_allclose(Hyy, Hrr, atol=1e-30)
+
+
 def test_hessian_matrix_3d():
     cube = cp.zeros((5, 5, 5))
     cube[2, 2, 2] = 4
@@ -170,6 +190,21 @@ def test_hessian_matrix_3d():
     # fmt: on
 
 
+@pytest.mark.parametrize('use_gaussian_derivatives', [False, True])
+def test_hessian_matrix_3d_xy(use_gaussian_derivatives):
+
+    img = cp.ones((5, 5, 5))
+
+    # order="xy" is only permitted for 2D
+    with pytest.raises(ValueError):
+        hessian_matrix(img, sigma=0.1, order="xy",
+                       use_gaussian_derivatives=use_gaussian_derivatives)
+
+    with pytest.raises(ValueError):
+        hessian_matrix(img, sigma=0.1, order='nonexistant',
+                       use_gaussian_derivatives=use_gaussian_derivatives)
+
+
 @pytest.mark.parametrize('dtype', [cp.float16, cp.float32, cp.float64])
 def test_structure_tensor_eigenvalues(dtype):
     square = cp.zeros((5, 5), dtype=dtype)
@@ -279,11 +314,37 @@ def test_custom_eigvals_kernels_vs_linalg_eigvalsh(shape, dtype):
 def test_hessian_matrix_det(approximate):
     image = cp.zeros((5, 5))
     image[2, 2] = 1
-    # TODO: approximate=True case not implemented
     det = hessian_matrix_det(image, 5, approximate=approximate)
     assert_array_almost_equal(det, 0, decimal=3)
 
 
+@pytest.mark.parametrize('approximate', [False, True])
+@pytest.mark.parametrize('ndim', [2, 3])
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+def test_hessian_matrix_det_vs_skimage(approximate, ndim, dtype):
+    if approximate and ndim != 2:
+        pytest.skip(reason="approximate only implemented for 2D images")
+    if approximate:
+        sigma = 3
+    else:
+        sigma = 1.5
+    rng = cp.random.default_rng(5)
+    if np.dtype(dtype).kind in 'iu':
+        image = rng.integers(0, 256, (16,) * ndim, dtype=dtype)
+    else:
+        image = rng.standard_normal((16,) * ndim).astype(dtype=dtype)
+    float_type = _supported_float_type(image.dtype)
+    expected = hessian_matrix_det_cpu(
+        cp.asnumpy(image), sigma, approximate=approximate
+    )
+    det = hessian_matrix_det(image, sigma, approximate=approximate)
+    assert det.dtype == float_type
+    tol = 1e-12 if det.dtype == np.float64 else 1e-6
+    cp.testing.assert_allclose(det, expected, rtol=tol, atol=tol)
+
+
 @pytest.mark.parametrize('dtype', [cp.float16, cp.float32, cp.float64])
 def test_hessian_matrix_det_3d(im3d, dtype):
     im3d = im3d.astype(dtype, copy=False)

python/cucim/src/cucim/skimage/metrics/_structural_similarity.py

@@ -36,7 +36,7 @@
 def _get_ssim_kernel():
 
     return cp.ElementwiseKernel(
-        in_params='F cov_norm, F ux, F uy, F uxx, F uyy, F uxy, float64 data_range, float64 K1, float64 K2',  # noqa
+        in_params='float64 cov_norm, F ux, F uy, F uxx, F uyy, F uxy, float64 data_range, float64 K1, float64 K2',  # noqa
         out_params='F ssim',
         operation=_ssim_operation,
         name='cucim_ssim'
@@ -47,7 +47,7 @@ def _get_ssim_kernel():
 def _get_ssim_grad_kernel():
 
     return cp.ElementwiseKernel(
-        in_params='F cov_norm, F ux, F uy, F uxx, F uyy, F uxy, float64 data_range, float64 K1, float64 K2',  # noqa
+        in_params='float64 cov_norm, F ux, F uy, F uxx, F uyy, F uxy, float64 data_range, float64 K1, float64 K2',  # noqa
         out_params='F ssim, F grad_temp1, F grad_temp2, F grad_temp3',
         operation=_ssim_operation + """
             grad_temp1 = A1 / D;
@@ -134,8 +134,9 @@ def structural_similarity(im1, im2,
     of Y), so one cannot calculate a channel-averaged SSIM with a single call
     to this function, as identical ranges are assumed for each channel.
 
-    To match the implementation of Wang et. al. [1]_, set `gaussian_weights`
-    to True, `sigma` to 1.5, and `use_sample_covariance` to False.
+    To match the implementation of Wang et al. [1]_, set `gaussian_weights`
+    to True, `sigma` to 1.5, `use_sample_covariance` to False, and
+    specify the `data_range` argument.
 
     .. versionchanged:: 0.16
         This function was renamed from ``skimage.measure.compare_ssim`` to
@@ -159,6 +160,12 @@ def structural_similarity(im1, im2,
     check_shape_equality(im1, im2)
     float_type = _supported_float_type(im1.dtype)
 
+    if isinstance(data_range, cp.ndarray):
+        if data_range.ndim != 0:
+            raise ValueError("data_range must be a scalar")
+        # need a host scalar
+        data_range = float(data_range)
+
     if channel_axis is not None:
         # loop over channels
         args = dict(win_size=win_size,
@@ -236,11 +243,22 @@ def structural_similarity(im1, im2,
         raise ValueError('Window size must be odd.')
 
     if data_range is None:
+        if (cp.issubdtype(im1.dtype, cp.floating) or
+            cp.issubdtype(im2.dtype, cp.floating)):
+            raise ValueError(
+                'Since image dtype is floating point, you must specify '
+                'the data_range parameter. Please read the documentation '
+                'carefully (including the note). It is recommended that '
+                'you always specify the data_range anyway.')
         if im1.dtype != im2.dtype:
-            warn("Inputs have mismatched dtype.  Setting data_range based on "
+            warn("Inputs have mismatched dtypes.  Setting data_range based on "
                  "im1.dtype.", stacklevel=2)
         dmin, dmax = dtype_range[im1.dtype.type]
-        data_range = dmax - dmin
+        data_range = float(dmax - dmin)
+        if cp.issubdtype(im1.dtype, cp.integer) and (im1.dtype != cp.uint8):
+            warn("Setting data_range based on im1.dtype. " +
+                 ("data_range = %.0f. " % data_range) +
+                 "Please specify data_range explicitly to avoid mistakes.", stacklevel=2)
 
     ndim = im1.ndim