Apply fixes to skimage.transform scheduled for scikit-image 0.19.2 (#208

)

Closes #199

Corresponds to:
scikit-image/scikit-image#6207
scikit-image/scikit-image#6211
scikit-image/scikit-image#6214

Additionally, the affines contained within the `PiecewiseAffineTransform` will have parameters in CuPy arrays when the inputs to `estimate` are CuPy arrays. This is one to make it consistent with the other transform classes.

@chrisroat, can you confirm if this fixes the issue for you?

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

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

URL: #208

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

@@ -763,7 +763,7 @@ def estimate(self, src, dst, weights=None):
         src_matrix, src, has_nan1 = _center_and_normalize_points(src)
         dst_matrix, dst, has_nan2 = _center_and_normalize_points(dst)
         if has_nan1 or has_nan2:
-            self.params = xp.full((d, d), np.nan)
+            self.params = xp.full((d + 1, d + 1), xp.nan)
             return False
         # params: a0, a1, a2, b0, b1, b2, c0, c1
         A = xp.zeros((n * d, (d + 1) ** 2))
@@ -792,6 +792,7 @@ def estimate(self, src, dst, weights=None):
         # because it is a rank-defective transform, which would map points
         # to a line rather than a plane.
         if xp.isclose(V[-1, -1], 0):
+            self.params = xp.full((d + 1, d + 1), xp.nan)
             return False
 
         H = np.zeros(
@@ -1050,39 +1051,38 @@ def estimate(self, src, dst):
         Returns
         -------
         success : bool
-            True, if model estimation succeeds.
+            True, if all pieces of the model are successfully estimated.
 
         """
 
         ndim = src.shape[1]
         # forward piecewise affine
         # triangulate input positions into mesh
         xp = cp.get_array_module(src)
-        if xp is cp:
-            # TODO: grlee77 :update if spatial.Delaunay is implemented for GPU
-            # transfer to CPU for use of spatial.Delaunay
-            src = cp.asnumpy(src)
-            dst = cp.asnumpy(dst)
+        # TODO: update if spatial.Delaunay become available in CuPy
+        self._tesselation = spatial.Delaunay(cp.asnumpy(src))
+
+        ok = True
 
-        self._tesselation = spatial.Delaunay(src)
         # find affine mapping from source positions to destination
         self.affines = []
-        for tri in self._tesselation.vertices:
+        for tri in xp.asarray(self._tesselation.vertices):
             affine = AffineTransform(dimensionality=ndim)
-            affine.estimate(src[tri, :], dst[tri, :])
+            ok &= affine.estimate(src[tri, :], dst[tri, :])
             self.affines.append(affine)
 
         # inverse piecewise affine
         # triangulate input positions into mesh
-        self._inverse_tesselation = spatial.Delaunay(dst)
+        # TODO: update if spatial.Delaunay become available in CuPy
+        self._inverse_tesselation = spatial.Delaunay(cp.asnumpy(dst))
         # find affine mapping from source positions to destination
         self.inverse_affines = []
-        for tri in self._inverse_tesselation.vertices:
+        for tri in xp.asarray(self._inverse_tesselation.vertices):
             affine = AffineTransform(dimensionality=ndim)
-            affine.estimate(dst[tri, :], src[tri, :])
+            ok &= affine.estimate(dst[tri, :], src[tri, :])
             self.inverse_affines.append(affine)
 
-        return True
+        return ok
 
     def __call__(self, coords):
         """Apply forward transformation.
@@ -1102,12 +1102,12 @@ def __call__(self, coords):
         """
 
         xp = cp.get_array_module(coords)
-        if xp == cp:
-            coords = coords.get()
-        out = np.empty_like(coords, np.double)
+        out = xp.empty_like(coords, xp.double)
 
         # determine triangle index for each coordinate
-        simplex = self._tesselation.find_simplex(coords)
+        # coords must be on host for calls to _tesselation methods
+        simplex = self._tesselation.find_simplex(cp.asnumpy(coords))
+        simplex = xp.asarray(simplex)
 
         # coordinates outside of mesh
         out[simplex == -1, :] = -1
@@ -1120,8 +1120,6 @@ def __call__(self, coords):
 
             out[index_mask, :] = affine(coords[index_mask, :])
 
-        if xp == cp:
-            out = xp.asarray(out)
         return out
 
     def inverse(self, coords):
@@ -1142,12 +1140,12 @@ def inverse(self, coords):
         """
 
         xp = cp.get_array_module(coords)
-        if xp == cp:
-            coords = coords.get()
-        out = np.empty_like(coords, np.double)
+        out = xp.empty_like(coords, xp.double)
 
         # determine triangle index for each coordinate
-        simplex = self._inverse_tesselation.find_simplex(coords)
+        # coords must be on host for calls to _tesselation methods
+        simplex = self._inverse_tesselation.find_simplex(cp.asnumpy(coords))
+        simplex = xp.asarray(simplex)
 
         # coordinates outside of mesh
         out[simplex == -1, :] = -1
@@ -1159,8 +1157,7 @@ def inverse(self, coords):
             index_mask = simplex == index
 
             out[index_mask, :] = affine(coords[index_mask, :])
-        if xp == cp:
-            out = xp.asarray(out)
+
         return out
 
 
@@ -1344,7 +1341,9 @@ def estimate(self, src, dst):
 
         self.params = _umeyama(src, dst, False)
 
-        return True
+        # _umeyama will return nan if the problem is not well-conditioned.
+        xp = cp.get_array_module(self.params)
+        return not xp.any(xp.isnan(self.params))
 
     @property
     def rotation(self):
@@ -1460,7 +1459,9 @@ def estimate(self, src, dst):
 
         self.params = _umeyama(src, dst, estimate_scale=True)
 
-        return True
+        # _umeyama will return nan if the problem is not well-conditioned.
+        xp = cp.get_array_module(self.params)
+        return not xp.any(xp.isnan(self.params))
 
     @property
     def scale(self):

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

@@ -70,7 +70,7 @@ def test_euclidean_estimation():
 
     # via estimate method
     tform3 = EuclideanTransform()
-    tform3.estimate(SRC, DST)
+    assert tform3.estimate(SRC, DST)
     assert_array_almost_equal(tform3.params, tform2.params)
 
 
@@ -117,7 +117,7 @@ def test_similarity_estimation():
 
     # via estimate method
     tform3 = SimilarityTransform()
-    tform3.estimate(SRC, DST)
+    assert tform3.estimate(SRC, DST)
     assert_array_almost_equal(tform3.params, tform2.params)
 
 
@@ -183,7 +183,7 @@ def test_affine_estimation():
 
     # via estimate method
     tform3 = AffineTransform()
-    tform3.estimate(SRC, DST)
+    assert tform3.estimate(SRC, DST)
     assert_array_almost_equal(tform3.params, tform2.params)
 
 
@@ -214,7 +214,7 @@ def test_affine_init():
 
 def test_piecewise_affine():
     tform = PiecewiseAffineTransform()
-    tform.estimate(SRC, DST)
+    assert tform.estimate(SRC, DST)
     # make sure each single affine transform is exactly estimated
     assert_array_almost_equal(tform(SRC), DST)
     assert_array_almost_equal(tform.inverse(DST), SRC)
@@ -357,7 +357,7 @@ def test_projective_estimation():
 
     # via estimate method
     tform3 = ProjectiveTransform()
-    tform3.estimate(SRC, DST)
+    assert tform3.estimate(SRC, DST)
     assert_array_almost_equal(tform3.params, tform2.params)
 
 
@@ -401,7 +401,7 @@ def test_polynomial_estimation():
 
     # via estimate method
     tform2 = PolynomialTransform()
-    tform2.estimate(SRC, DST, order=10)
+    assert tform2.estimate(SRC, DST, order=10)
     assert_array_almost_equal(tform2.params, tform.params)
 
 
@@ -557,15 +557,19 @@ def test_degenerate(xp=cp):
     src = dst = xp.zeros((10, 2))
 
     tform = SimilarityTransform()
-    tform.estimate(src, dst)
+    assert not tform.estimate(src, dst)
+    assert xp.all(xp.isnan(tform.params))
+
+    tform = EuclideanTransform()
+    assert not tform.estimate(src, dst)
     assert xp.all(xp.isnan(tform.params))
 
     tform = AffineTransform()
-    tform.estimate(src, dst)
+    assert not tform.estimate(src, dst)
     assert xp.all(xp.isnan(tform.params))
 
     tform = ProjectiveTransform()
-    tform.estimate(src, dst)
+    assert not tform.estimate(src, dst)
     assert xp.all(xp.isnan(tform.params))
 
     # See gh-3926 for discussion details
@@ -581,6 +585,32 @@ def test_degenerate(xp=cp):
         # a transform could be returned with nan values.
         assert not tform.estimate(src, dst) or xp.isfinite(tform.params).all()
 
+    src = xp.array([[0, 2, 0], [0, 2, 0], [0, 4, 0]])
+    dst = xp.array([[0, 1, 0], [0, 1, 0], [0, 3, 0]])
+    tform = AffineTransform()
+    assert not tform.estimate(src, dst)
+    assert xp.all(xp.isnan(tform.params))
+
+    # The tesselation on the following points produces one degenerate affine
+    # warp within PiecewiseAffineTransform.
+    src = xp.asarray([
+        [0, 192, 256], [0, 256, 256], [5, 0, 192], [5, 64, 0], [5, 64, 64],
+        [5, 64, 256], [5, 192, 192], [5, 256, 256], [0, 192, 256],
+    ])
+
+    dst = xp.asarray([
+        [0, 142, 206], [0, 206, 206], [5, -50, 142], [5, 14, 0], [5, 14, 64],
+        [5, 14, 206], [5, 142, 142], [5, 206, 206], [0, 142, 206],
+    ])
+    tform = PiecewiseAffineTransform()
+    assert not tform.estimate(src, dst)
+    assert np.all(np.isnan(tform.affines[4].params))  # degenerate affine
+    for idx, affine in enumerate(tform.affines):
+        if idx != 4:
+            assert not xp.all(xp.isnan(affine.params))
+    for affine in tform.inverse_affines:
+        assert not xp.all(xp.isnan(affine.params))
+
 
 @pytest.mark.parametrize('xp', [np, cp])
 def test_normalize_degenerate_points(xp):
@@ -659,7 +689,7 @@ def test_estimate_affine_3d(xp):
     dst = tf(src)
     dst_noisy = dst + xp.asarray(np.random.random((25, ndim)))
     tf2 = AffineTransform(dimensionality=ndim)
-    tf2.estimate(src, dst_noisy)
+    assert tf2.estimate(src, dst_noisy)
     # we check rot/scale/etc more tightly than translation because translation
     # estimation is on the 1 pixel scale
     assert_array_almost_equal(tf2.params[:, :-1], matrix[:, :-1], decimal=2)