Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Merged
merged 2 commits into from
Mar 25, 2022
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
be2ede5
apply fixes to geometric transform estimation
grlee77 Jan 26, 2022
2c0270f
Apply upstream fix to PiecewiseAffineTransform.estimate
grlee77 Jan 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 29 additions & 28 deletions python/cucim/src/cucim/skimage/transform/_geometric.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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))
Comment on lines +1062 to +1063
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we open an issue about this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could, it is probably not very easy to address though. We would need a Delaunay triangulation class with a vertices attribute and find_simplex method for the GPU (There are several other attributes/methods on this class in SciPy, but they are not used here). On the SciPy side the implementation relies on the QHull C library.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sure there are GPU implementations of Delaunay triangulation, but hadn't looked into whether there is a compatibly licensed one we could reuse.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we did find one, ideally it would be available via upstream cupyx.scipy.spatial

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Idk maybe it could happen. There has been work in the past like CudaHull


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.
Expand All @@ -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))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we coerce to a NumPy array here?

Note: This comes up again below, but only commenting here for simplicity.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self._tesselation is of class Delaunay from SciPy so this find_simplex method is wrapping some CPU implementation from QHull under the hood.

simplex = xp.asarray(simplex)

# coordinates outside of mesh
out[simplex == -1, :] = -1
Expand All @@ -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):
Expand All @@ -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
Expand All @@ -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


Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down
50 changes: 40 additions & 10 deletions python/cucim/src/cucim/skimage/transform/tests/test_geometric.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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
Expand All @@ -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):
Expand Down Expand Up @@ -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)
Expand Down