diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 3df5863d..1b4153f6 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -5,7 +5,7 @@ from .radon_transform import radon, iradon, iradon_sart from .finite_radon_transform import frt2, ifrt2 from .integral import integral_image, integrate from ._geometric import (warp, warp_coords, estimate_transform, - matrix_transform, + matrix_transform, EuclideanTransform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform, PiecewiseAffineTransform) @@ -31,6 +31,7 @@ __all__ = ['hough_circle', 'warp_coords', 'estimate_transform', 'matrix_transform', + 'EuclideanTransform', 'SimilarityTransform', 'AffineTransform', 'ProjectiveTransform', diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 3793c9f5..e73decdc 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -1,3 +1,4 @@ +from __future__ import division import six import math import numpy as np @@ -67,8 +68,83 @@ def _center_and_normalize_points(points): return matrix, new_points +def _umeyama(src, dst, estimate_scale): + """Estimate N-D similarity transformation with or without scaling. + + Parameters + ---------- + src : (M, N) array + Source coordinates. + dst : (M, N) array + Destination coordinates. + estimate_scale : bool + Whether to estimate scaling factor. + + Returns + ------- + T : (N + 1, N + 1) + The homogeneous similarity transformation matrix. The matrix contains + NaN values only if the problem is not well-conditioned. + + References + ---------- + .. [1] "Least-squares estimation of transformation parameters between two + point patterns", Shinji Umeyama, PAMI 1991, DOI: 10.1109/34.88573 + + """ + + num = src.shape[0] + dim = src.shape[1] + + # Compute mean of src and dst. + src_mean = src.mean(axis=0) + dst_mean = dst.mean(axis=0) + + # Subtract mean from src and dst. + src_demean = src - src_mean + dst_demean = dst - dst_mean + + # Eq. (38). + A = np.dot(dst_demean.T, src_demean) / num + + # Eq. (39). + d = np.ones((dim,), dtype=np.double) + if np.linalg.det(A) < 0: + d[dim - 1] = -1 + + T = np.eye(dim + 1, dtype=np.double) + + U, S, V = np.linalg.svd(A) + + # Eq. (40) and (43). + rank = np.linalg.matrix_rank(A) + if rank == 0: + return np.nan * T + elif rank == dim - 1: + if np.linalg.det(U) * np.linalg.det(V) > 0: + T[:dim, :dim] = np.dot(U, V) + else: + s = d[dim - 1] + d[dim - 1] = -1 + T[:dim, :dim] = np.dot(U, np.dot(np.diag(d), V)) + d[dim - 1] = s + else: + T[:dim, :dim] = np.dot(U, np.dot(np.diag(d), V.T)) + + if estimate_scale: + # Eq. (41) and (42). + scale = 1.0 / src_demean.var(axis=0).sum() * np.dot(S, d) + else: + scale = 1.0 + + T[:dim, dim] = dst_mean - scale * np.dot(T[:dim, :dim], src_mean.T) + T[:dim, :dim] *= scale + + return T + + class GeometricTransform(object): - """Perform geometric transformations on a set of coordinates. + """Base class for geometric transformations. """ def __call__(self, coords): @@ -216,8 +292,7 @@ class ProjectiveTransform(GeometricTransform): return self._apply_mat(coords, self._inv_matrix) def estimate(self, src, dst): - """Set the transformation matrix with the explicit transformation - parameters. + """Estimate the transformation from a set of corresponding points. You can determine the over-, well- and under-determined parameters with the total least-squares method. @@ -341,7 +416,6 @@ class ProjectiveTransform(GeometricTransform): class AffineTransform(ProjectiveTransform): - """2D affine transformation of the form: X = a0*x + a1*y + a2 = @@ -350,7 +424,7 @@ class AffineTransform(ProjectiveTransform): Y = b0*x + b1*y + b2 = = sx*x*sin(rotation) + sy*y*cos(rotation + shear) + b2 - where ``sx`` and ``sy`` are zoom factors in the x and y directions, + where ``sx`` and ``sy`` are scale factors in the x and y directions, and the homogeneous transformation matrix is:: [[a0 a1 a2] @@ -433,7 +507,6 @@ class AffineTransform(ProjectiveTransform): class PiecewiseAffineTransform(GeometricTransform): - """2D piecewise affine transformation. Control points are used to define the mapping. The transform is based on @@ -456,7 +529,7 @@ class PiecewiseAffineTransform(GeometricTransform): self.inverse_affines = None def estimate(self, src, dst): - """Set the control points with which to perform the piecewise mapping. + """Estimate the transformation from a set of corresponding points. Number of source and destination coordinates must match. @@ -567,21 +640,122 @@ class PiecewiseAffineTransform(GeometricTransform): return out -class SimilarityTransform(ProjectiveTransform): - """2D similarity transformation of the form: +class EuclideanTransform(ProjectiveTransform): + """2D Euclidean transformation of the form: X = a0 * x - b0 * y + a1 = - = m * x * cos(rotation) - m * y * sin(rotation) + a1 + = x * cos(rotation) - y * sin(rotation) + a1 Y = b0 * x + a0 * y + b1 = - = m * x * sin(rotation) + m * y * cos(rotation) + b1 + = x * sin(rotation) + y * cos(rotation) + b1 - where ``m`` is a zoom factor and the homogeneous transformation matrix is:: + where the homogeneous transformation matrix is:: [[a0 b0 a1] [b0 a0 b1] [0 0 1]] + The Euclidean transformation is a rigid transformation with rotation and + translation parameters. The similarity transformation extends the Euclidean + transformation with a single scaling factor. + + Parameters + ---------- + matrix : (3, 3) array, optional + Homogeneous transformation matrix. + rotation : float, optional + Rotation angle in counter-clockwise direction as radians. + translation : (tx, ty) as array, list or tuple, optional + x, y translation parameters. + + Attributes + ---------- + params : (3, 3) array + Homogeneous transformation matrix. + + """ + + def __init__(self, matrix=None, rotation=None, translation=None): + params = any(param is not None + for param in (rotation, translation)) + + if params and matrix is not None: + raise ValueError("You cannot specify the transformation matrix and" + " the implicit parameters at the same time.") + elif matrix is not None: + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + self.params = matrix + elif params: + if rotation is None: + rotation = 0 + if translation is None: + translation = (0, 0) + + self.params = np.array([ + [math.cos(rotation), - math.sin(rotation), 0], + [math.sin(rotation), math.cos(rotation), 0], + [ 0, 0, 1] + ]) + self.params[0:2, 2] = translation + else: + # default to an identity transform + self.params = np.eye(3) + + def estimate(self, src, dst): + """Estimate the transformation from a set of corresponding points. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + self.params = _umeyama(src, dst, False) + + return True + + @property + def rotation(self): + return math.atan2(self.params[1, 0], self.params[1, 1]) + + @property + def translation(self): + return self.params[0:2, 2] + + +class SimilarityTransform(EuclideanTransform): + """2D similarity transformation of the form: + + X = a0 * x - b0 * y + a1 = + = s * x * cos(rotation) - s * y * sin(rotation) + a1 + + Y = b0 * x + a0 * y + b1 = + = s * x * sin(rotation) + s * y * cos(rotation) + b1 + + where ``s`` is a scale factor and the homogeneous transformation matrix is:: + + [[a0 b0 a1] + [b0 a0 b1] + [0 0 1]] + + The similarity transformation extends the Euclidean transformation with a + single scaling factor in addition to the rotation and translation + parameters. + Parameters ---------- matrix : (3, 3) array, optional @@ -632,38 +806,13 @@ class SimilarityTransform(ProjectiveTransform): self.params = np.eye(3) def estimate(self, src, dst): - """Set the transformation matrix with the explicit parameters. + """Estimate the transformation from a set of corresponding points. You can determine the over-, well- and under-determined parameters with the total least-squares method. Number of source and destination coordinates must match. - The transformation is defined as:: - - X = a0 * x - b0 * y + a1 - Y = b0 * x + a0 * y + b1 - - These equations can be transformed to the following form:: - - 0 = a0 * x - b0 * y + a1 - X - 0 = b0 * x + a0 * y + b1 - Y - - which exist for each set of corresponding points, so we have a set of - N * 2 equations. The coefficients appear linearly so we can write - A x = 0, where:: - - A = [[x 1 -y 0 -X] - [y 0 x 1 -Y] - ... - ... - ] - x.T = [a0 a1 b0 b1 c3] - - In case of total least-squares the solution of this homogeneous system - of equations is the right singular vector of A which corresponds to the - smallest singular value normed by the coefficient c3. - Parameters ---------- src : (N, 2) array @@ -678,44 +827,7 @@ class SimilarityTransform(ProjectiveTransform): """ - 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 False - - xs = src[:, 0] - ys = src[:, 1] - xd = dst[:, 0] - yd = dst[:, 1] - rows = src.shape[0] - - # params: a0, a1, b0, b1 - A = np.zeros((rows * 2, 5)) - A[:rows, 0] = xs - A[:rows, 2] = - ys - A[:rows, 1] = 1 - A[rows:, 2] = xs - A[rows:, 0] = ys - A[rows:, 3] = 1 - A[:rows, 4] = xd - A[rows:, 4] = yd - - _, _, V = np.linalg.svd(A) - - # solution is right singular vector that corresponds to smallest - # singular value - a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -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 + self.params = _umeyama(src, dst, True) return True @@ -728,17 +840,9 @@ class SimilarityTransform(ProjectiveTransform): scale = self.params[0, 0] / math.cos(self.rotation) return scale - @property - def rotation(self): - return math.atan2(self.params[1, 0], self.params[1, 1]) - - @property - def translation(self): - return self.params[0:2, 2] - class PolynomialTransform(GeometricTransform): - """2D transformation of the form: + """2D polynomial transformation of the form: X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) @@ -766,8 +870,7 @@ class PolynomialTransform(GeometricTransform): self.params = params def estimate(self, src, dst, order=2): - """Set the transformation matrix with the explicit transformation - parameters. + """Estimate the transformation from a set of corresponding points. You can determine the over-, well- and under-determined parameters with the total least-squares method. @@ -885,6 +988,7 @@ class PolynomialTransform(GeometricTransform): TRANSFORMS = { + 'euclidean': EuclideanTransform, 'similarity': SimilarityTransform, 'affine': AffineTransform, 'piecewise-affine': PiecewiseAffineTransform, @@ -909,13 +1013,14 @@ def estimate_transform(ttype, src, dst, **kwargs): Parameters ---------- - ttype : {'similarity', 'affine', 'piecewise-affine', 'projective', \ - 'polynomial'} + ttype : {'euclidean', similarity', 'affine', 'piecewise-affine', \ + 'projective', 'polynomial'} Type of transform. kwargs : array or int Function parameters (src, dst, n, angle):: NAME / TTYPE FUNCTION PARAMETERS + 'euclidean' `src, `dst` 'similarity' `src, `dst` 'affine' `src, `dst` 'piecewise-affine' `src, `dst` diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index a12eb50b..98bb7524 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -4,9 +4,9 @@ from numpy.testing import (assert_equal, assert_almost_equal, from skimage.transform._geometric import _stackcopy from skimage.transform._geometric import GeometricTransform from skimage.transform import (estimate_transform, matrix_transform, - SimilarityTransform, AffineTransform, - ProjectiveTransform, PolynomialTransform, - PiecewiseAffineTransform) + EuclideanTransform, SimilarityTransform, + AffineTransform, ProjectiveTransform, + PolynomialTransform, PiecewiseAffineTransform) from skimage._shared._warnings import expected_warnings @@ -42,7 +42,8 @@ def test_stackcopy(): def test_estimate_transform(): - for tform in ('similarity', 'affine', 'projective', 'polynomial'): + for tform in ('euclidean', 'similarity', 'affine', 'projective', + 'polynomial'): estimate_transform(tform, SRC[:2, :], DST[:2, :]) assert_raises(ValueError, estimate_transform, 'foobar', SRC[:2, :], DST[:2, :]) @@ -53,18 +54,65 @@ def test_matrix_transform(): assert_equal(tform(SRC), matrix_transform(SRC, tform.params)) +def test_euclidean_estimation(): + # exact solution + tform = estimate_transform('euclidean', SRC[:2, :], SRC[:2, :] + 10) + assert_almost_equal(tform(SRC[:2, :]), SRC[:2, :] + 10) + assert_almost_equal(tform.params[0, 0], tform.params[1, 1]) + assert_almost_equal(tform.params[0, 1], - tform.params[1, 0]) + + # over-determined + tform2 = estimate_transform('euclidean', SRC, DST) + assert_almost_equal(tform2.inverse(tform2(SRC)), SRC) + assert_almost_equal(tform2.params[0, 0], tform2.params[1, 1]) + assert_almost_equal(tform2.params[0, 1], - tform2.params[1, 0]) + + # via estimate method + tform3 = EuclideanTransform() + tform3.estimate(SRC, DST) + assert_almost_equal(tform3.params, tform2.params) + + +def test_euclidean_init(): + # init with implicit parameters + rotation = 1 + translation = (1, 1) + tform = EuclideanTransform(rotation=rotation, translation=translation) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) + + # init with transformation matrix + tform2 = EuclideanTransform(tform.params) + assert_almost_equal(tform2.rotation, rotation) + assert_almost_equal(tform2.translation, translation) + + # test special case for scale if rotation=0 + rotation = 0 + translation = (1, 1) + tform = EuclideanTransform(rotation=rotation, translation=translation) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) + + # test special case for scale if rotation=90deg + rotation = np.pi / 2 + translation = (1, 1) + tform = EuclideanTransform(rotation=rotation, translation=translation) + assert_almost_equal(tform.rotation, rotation) + assert_almost_equal(tform.translation, translation) + + def test_similarity_estimation(): # exact solution tform = estimate_transform('similarity', 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]) + assert_almost_equal(tform.params[0, 0], tform.params[1, 1]) + assert_almost_equal(tform.params[0, 1], - tform.params[1, 0]) # over-determined tform2 = estimate_transform('similarity', SRC, DST) 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]) + assert_almost_equal(tform2.params[0, 0], tform2.params[1, 1]) + assert_almost_equal(tform2.params[0, 1], - tform2.params[1, 0]) # via estimate method tform3 = SimilarityTransform() @@ -240,11 +288,14 @@ def test_invalid_input(): assert_raises(ValueError, ProjectiveTransform, np.zeros((2, 3))) assert_raises(ValueError, AffineTransform, np.zeros((2, 3))) assert_raises(ValueError, SimilarityTransform, np.zeros((2, 3))) + assert_raises(ValueError, EuclideanTransform, np.zeros((2, 3))) assert_raises(ValueError, AffineTransform, matrix=np.zeros((2, 3)), scale=1) assert_raises(ValueError, SimilarityTransform, matrix=np.zeros((2, 3)), scale=1) + assert_raises(ValueError, EuclideanTransform, + matrix=np.zeros((2, 3)), translation=(0, 0)) assert_raises(ValueError, PolynomialTransform, np.zeros((3, 3)))