diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index c2a71f71..1b960d4b 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -239,16 +239,12 @@ def test_ransac_dynamic_max_trials(): def test_ransac_invalid_input(): - assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=-1, - residual_threshold=0) assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2, - residual_threshold=0, max_trials=-1) + residual_threshold=0, max_trials=-1) assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2, - residual_threshold=0, - stop_probability=-1) + residual_threshold=0, stop_probability=-1) assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2, - residual_threshold=0, - stop_probability=1.01) + residual_threshold=0, stop_probability=1.01) def test_deprecated_params_attribute(): diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index e658d23e..3a2972c0 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -10,6 +10,54 @@ from ..exposure import rescale_intensity from ._warps_cy import _warp_fast +def _center_and_normalize_points(points): + """Center and normalize image points. + + The points are transformed in a two-step procedure that is expressed + as a transformation matrix. The matrix of the resulting points is usually + better conditioned than the matrix of the original points. + + Center the image points, such that the new coordinate system has its + origin at the centroid of the image points. + + Normalize the image points, such that the mean distance from the points + to the origin of the coordinate system is sqrt(2). + + Parameters + ---------- + points : (N, 2) array + The coordinates of the image points. + + Returns + ------- + matrix : (3, 3) array + The transformation matrix to obtain the new points. + new_points : (N, 2) array + The transformed image points. + + """ + + centroid = np.mean(points, axis=0) + + rms = math.sqrt(np.sum((points - centroid) ** 2) / points.shape[0]) + + norm_factor = math.sqrt(2) / rms + + matrix = np.array([[norm_factor, 0, -norm_factor * centroid[0]], + [0, norm_factor, -norm_factor * centroid[1]], + [0, 0, 1]]) + + pointsh = np.row_stack([points.T, np.ones((points.shape[0]),)]) + + new_pointsh = np.dot(matrix, pointsh).T + + new_points = new_pointsh[:, :2] + new_points[:, 0] /= new_pointsh[:, 2] + new_points[:, 1] /= new_pointsh[:, 2] + + return matrix, new_points + + class GeometricTransform(object): """Perform geometric transformations on a set of coordinates. @@ -216,6 +264,14 @@ class ProjectiveTransform(GeometricTransform): Destination coordinates. """ + + try: + src_matrix, src = _center_and_normalize_points(src) + dst_matrix, dst = _center_and_normalize_points(dst) + except ZeroDivisionError: + self.params = np.nan * np.empty((3, 3)) + return + xs = src[:, 0] ys = src[:, 1] xd = dst[:, 0] @@ -248,6 +304,9 @@ class ProjectiveTransform(GeometricTransform): H.flat[list(self._coeffs) + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 + # De-center and de-normalize + H = np.dot(np.linalg.inv(dst_matrix), np.dot(H, src_matrix)) + self.params = H def __add__(self, other): @@ -596,6 +655,14 @@ class SimilarityTransform(ProjectiveTransform): Destination coordinates. """ + + try: + src_matrix, src = _center_and_normalize_points(src) + dst_matrix, dst = _center_and_normalize_points(dst) + except ZeroDivisionError: + self.params = np.nan * np.empty((3, 3)) + return + xs = src[:, 0] ys = src[:, 1] xd = dst[:, 0] @@ -619,9 +686,15 @@ class SimilarityTransform(ProjectiveTransform): # singular value a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1] - self.params = np.array([[a0, -b0, a1], - [b0, a0, b1], - [ 0, 0, 1]]) + S = np.array([[a0, -b0, a1], + [b0, a0, b1], + [ 0, 0, 1]]) + + # De-center and de-normalize + S = np.dot(np.linalg.inv(dst_matrix), np.dot(S, src_matrix)) + + self.params = S + @property def scale(self): diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index a27c1d42..5a9c8aa2 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -88,18 +88,23 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True, else: # 2-dimensional interpolation - # 3 control points necessary to estimate exact AffineTransform - src_corners = np.array([[1, 1], [1, rows], [cols, rows]]) - 1 - dst_corners = np.zeros(src_corners.shape, dtype=np.double) - # take into account that 0th pixel is at position (0.5, 0.5) - dst_corners[:, 0] = col_scale * (src_corners[:, 0] + 0.5) - 0.5 - dst_corners[:, 1] = row_scale * (src_corners[:, 1] + 0.5) - 0.5 + if rows == 1 and cols == 1: + tform = AffineTransform(translation=(orig_cols / 2.0 - 0.5, + orig_rows / 2.0 - 0.5)) + else: + # 3 control points necessary to estimate exact AffineTransform + src_corners = np.array([[1, 1], [1, rows], [cols, rows]]) - 1 + dst_corners = np.zeros(src_corners.shape, dtype=np.double) + # take into account that 0th pixel is at position (0.5, 0.5) + dst_corners[:, 0] = col_scale * (src_corners[:, 0] + 0.5) - 0.5 + dst_corners[:, 1] = row_scale * (src_corners[:, 1] + 0.5) - 0.5 - tform = AffineTransform() - tform.estimate(src_corners, dst_corners) + tform = AffineTransform() + tform.estimate(src_corners, dst_corners) out = warp(image, tform, output_shape=output_shape, order=order, - mode=mode, cval=cval, clip=clip, preserve_range=preserve_range) + mode=mode, cval=cval, clip=clip, + preserve_range=preserve_range) return out diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 125faf6c..d34e5e8f 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -1,5 +1,5 @@ import numpy as np -from numpy.testing import (assert_equal, assert_array_almost_equal, +from numpy.testing import (assert_equal, assert_almost_equal, assert_raises) from skimage.transform._geometric import _stackcopy from skimage.transform._geometric import GeometricTransform @@ -38,7 +38,7 @@ def test_stackcopy(): y = np.eye(3, 3) _stackcopy(x, y) for i in range(layers): - assert_array_almost_equal(x[..., i], y) + assert_almost_equal(x[..., i], y) def test_estimate_transform(): @@ -56,20 +56,20 @@ def test_matrix_transform(): def test_similarity_estimation(): # exact solution tform = estimate_transform('similarity', SRC[:2, :], DST[:2, :]) - assert_array_almost_equal(tform(SRC[:2, :]), DST[:2, :]) + assert_almost_equal(tform(SRC[:2, :]), DST[:2, :]) assert_equal(tform.params[0, 0], tform.params[1, 1]) assert_equal(tform.params[0, 1], - tform.params[1, 0]) # over-determined tform2 = estimate_transform('similarity', SRC, DST) - assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + assert_almost_equal(tform2.inverse(tform2(SRC)), SRC) assert_equal(tform2.params[0, 0], tform2.params[1, 1]) assert_equal(tform2.params[0, 1], - tform2.params[1, 0]) # via estimate method tform3 = SimilarityTransform() tform3.estimate(SRC, DST) - assert_array_almost_equal(tform3.params, tform2.params) + assert_almost_equal(tform3.params, tform2.params) def test_similarity_init(): @@ -79,15 +79,15 @@ def test_similarity_init(): translation = (1, 1) tform = SimilarityTransform(scale=scale, rotation=rotation, translation=translation) - assert_array_almost_equal(tform.scale, scale) - assert_array_almost_equal(tform.rotation, rotation) - assert_array_almost_equal(tform.translation, translation) + assert_almost_equal(tform.scale, scale) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) # init with transformation matrix tform2 = SimilarityTransform(tform.params) - assert_array_almost_equal(tform2.scale, scale) - assert_array_almost_equal(tform2.rotation, rotation) - assert_array_almost_equal(tform2.translation, translation) + assert_almost_equal(tform2.scale, scale) + assert_almost_equal(tform2.rotation, rotation) + assert_almost_equal(tform2.translation, translation) # test special case for scale if rotation=0 scale = 0.1 @@ -95,9 +95,9 @@ def test_similarity_init(): translation = (1, 1) tform = SimilarityTransform(scale=scale, rotation=rotation, translation=translation) - assert_array_almost_equal(tform.scale, scale) - assert_array_almost_equal(tform.rotation, rotation) - assert_array_almost_equal(tform.translation, translation) + assert_almost_equal(tform.scale, scale) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) # test special case for scale if rotation=90deg @@ -106,24 +106,24 @@ def test_similarity_init(): translation = (1, 1) tform = SimilarityTransform(scale=scale, rotation=rotation, translation=translation) - assert_array_almost_equal(tform.scale, scale) - assert_array_almost_equal(tform.rotation, rotation) - assert_array_almost_equal(tform.translation, translation) + assert_almost_equal(tform.scale, scale) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) def test_affine_estimation(): # exact solution tform = estimate_transform('affine', SRC[:3, :], DST[:3, :]) - assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :]) + assert_almost_equal(tform(SRC[:3, :]), DST[:3, :]) # over-determined tform2 = estimate_transform('affine', SRC, DST) - assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + assert_almost_equal(tform2.inverse(tform2(SRC)), SRC) # via estimate method tform3 = AffineTransform() tform3.estimate(SRC, DST) - assert_array_almost_equal(tform3.params, tform2.params) + assert_almost_equal(tform3.params, tform2.params) def test_affine_init(): @@ -134,71 +134,71 @@ def test_affine_init(): translation = (1, 1) tform = AffineTransform(scale=scale, rotation=rotation, shear=shear, translation=translation) - assert_array_almost_equal(tform.scale, scale) - assert_array_almost_equal(tform.rotation, rotation) - assert_array_almost_equal(tform.shear, shear) - assert_array_almost_equal(tform.translation, translation) + assert_almost_equal(tform.scale, scale) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.shear, shear) + assert_almost_equal(tform.translation, translation) # init with transformation matrix tform2 = AffineTransform(tform.params) - assert_array_almost_equal(tform2.scale, scale) - assert_array_almost_equal(tform2.rotation, rotation) - assert_array_almost_equal(tform2.shear, shear) - assert_array_almost_equal(tform2.translation, translation) + assert_almost_equal(tform2.scale, scale) + assert_almost_equal(tform2.rotation, rotation) + assert_almost_equal(tform2.shear, shear) + assert_almost_equal(tform2.translation, translation) def test_piecewise_affine(): tform = PiecewiseAffineTransform() 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) + assert_almost_equal(tform(SRC), DST) + assert_almost_equal(tform.inverse(DST), SRC) def test_projective_estimation(): # exact solution tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) - assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :]) + assert_almost_equal(tform(SRC[:4, :]), DST[:4, :]) # over-determined tform2 = estimate_transform('projective', SRC, DST) - assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + assert_almost_equal(tform2.inverse(tform2(SRC)), SRC) # via estimate method tform3 = ProjectiveTransform() tform3.estimate(SRC, DST) - assert_array_almost_equal(tform3.params, tform2.params) + assert_almost_equal(tform3.params, tform2.params) def test_projective_init(): tform = estimate_transform('projective', SRC, DST) # init with transformation matrix tform2 = ProjectiveTransform(tform.params) - assert_array_almost_equal(tform2.params, tform.params) + assert_almost_equal(tform2.params, tform.params) def test_polynomial_estimation(): # over-determined tform = estimate_transform('polynomial', SRC, DST, order=10) - assert_array_almost_equal(tform(SRC), DST, 6) + assert_almost_equal(tform(SRC), DST, 6) # via estimate method tform2 = PolynomialTransform() tform2.estimate(SRC, DST, order=10) - assert_array_almost_equal(tform2.params, tform.params) + assert_almost_equal(tform2.params, tform.params) def test_polynomial_init(): tform = estimate_transform('polynomial', SRC, DST, order=10) # init with transformation parameters tform2 = PolynomialTransform(tform.params) - assert_array_almost_equal(tform2.params, tform.params) + assert_almost_equal(tform2.params, tform.params) def test_polynomial_default_order(): tform = estimate_transform('polynomial', SRC, DST) tform2 = estimate_transform('polynomial', SRC, DST, order=2) - assert_array_almost_equal(tform2.params, tform.params) + assert_almost_equal(tform2.params, tform.params) def test_polynomial_inverse(): @@ -210,17 +210,17 @@ def test_union(): tform2 = SimilarityTransform(scale=0.1, rotation=0.9) tform3 = SimilarityTransform(scale=0.1 ** 2, rotation=0.3 + 0.9) tform = tform1 + tform2 - assert_array_almost_equal(tform.params, tform3.params) + assert_almost_equal(tform.params, tform3.params) tform1 = AffineTransform(scale=(0.1, 0.1), rotation=0.3) tform2 = SimilarityTransform(scale=0.1, rotation=0.9) tform3 = SimilarityTransform(scale=0.1 ** 2, rotation=0.3 + 0.9) tform = tform1 + tform2 - assert_array_almost_equal(tform.params, tform3.params) + assert_almost_equal(tform.params, tform3.params) assert tform.__class__ == ProjectiveTransform tform = AffineTransform(scale=(0.1, 0.1), rotation=0.3) - assert_array_almost_equal((tform + tform.inverse).params, np.eye(3)) + assert_almost_equal((tform + tform.inverse).params, np.eye(3)) def test_union_differing_types(): @@ -260,6 +260,22 @@ def test_deprecated_params_attributes(): assert_equal(tform._params, tform.params) +def test_degenerate(): + src = dst = np.zeros((10, 2)) + + tform = SimilarityTransform() + tform.estimate(src, dst) + assert np.all(np.isnan(tform.params)) + + tform = AffineTransform() + tform.estimate(src, dst) + assert np.all(np.isnan(tform.params)) + + tform = ProjectiveTransform() + tform.estimate(src, dst) + assert np.all(np.isnan(tform.params)) + + if __name__ == "__main__": from numpy.testing import run_module_suite run_module_suite()