From eb1e71114ca161615d7bce080ae9a2ccc0bf465a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 8 Aug 2012 19:13:45 +0200 Subject: [PATCH] fix and improve estimation of geometric transformation parameters Design matrix was not composed correctly as functional model was incorrect. Additionally estimation is now based on total least-squares method. --- skimage/transform/_geometric.py | 49 ++++++++++++++++++----- skimage/transform/tests/test_geometric.py | 6 +-- 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c1c39b2b..0c767cd7 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -133,6 +133,11 @@ class ProjectiveTransform(GeometricTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -148,7 +153,7 @@ class ProjectiveTransform(GeometricTransform): rows = src.shape[0] #: params: a0, a1, a2, b0, b1, b2, c0, c1 - A = np.zeros((rows * 2, 8)) + A = np.zeros((rows * 2, 9)) A[:rows, 0] = xs A[:rows, 1] = ys A[:rows, 2] = 1 @@ -159,14 +164,18 @@ class ProjectiveTransform(GeometricTransform): A[rows:, 5] = 1 A[rows:, 6] = - yd * xs A[rows:, 7] = - yd * ys + A[:rows, 8] = xd + A[rows:, 8] = yd # Select relevant columns, depending on coeffs - A = A[:, self._coefs] + A = A[:, self._coefs + [8]] - b = np.hstack([xd, yd]) + _, _, V = np.linalg.svd(A) H = np.zeros((3, 3)) - H.flat[self._coefs] = np.linalg.lstsq(A, b)[0] + # solution is eigen vector that corresponds to smallest eigen value + # and normed by c3 + H.flat[self._coefs + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 self._matrix = H @@ -294,6 +303,11 @@ class SimilarityTransform(ProjectiveTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -309,17 +323,20 @@ class SimilarityTransform(ProjectiveTransform): rows = src.shape[0] #: params: a0, a1, b0, b1 - A = np.zeros((rows * 2, 4)) + 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 - b = np.hstack([xd, yd]) + _, _, V = np.linalg.svd(A) + + a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1] - a0, a1, b0, b1 = np.linalg.lstsq(A, b)[0] self._matrix = np.array([[a0, -b0, a1], [b0, a0, b1], [ 0, 0, 1]]) @@ -392,6 +409,11 @@ class PolynomialTransform(GeometricTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -411,7 +433,7 @@ class PolynomialTransform(GeometricTransform): # number of unknown polynomial coefficients u = (order + 1) * (order + 2) - A = np.zeros((rows * 2, u)) + A = np.zeros((rows * 2, u + 1)) pidx = 0 for j in xrange(order + 1): for i in xrange(j + 1): @@ -419,9 +441,14 @@ class PolynomialTransform(GeometricTransform): A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i pidx += 1 - b = np.hstack([xd, yd]) + A[:rows, -1] = xd + A[rows:, -1] = yd - self._coeffs = np.linalg.lstsq(A, b)[0].reshape((2, u / 2)) + _, _, V = np.linalg.svd(A) + + coeffs = - V[-1, :-1] / V[-1, -1] + + self._coeffs = coeffs.reshape((2, u / 2)) def __call__(self, coords): """Apply forward transformation. @@ -473,7 +500,7 @@ def estimate_transform(ttype, src, dst, **kwargs): """Estimate 2D geometric transformation parameters. You can determine the over-, well- and under-determined parameters - with the least-squares method. + with the total least-squares method. Number of source must match number of destination coordinates. diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index f99bab7a..c7305956 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -42,7 +42,6 @@ def test_similarity_estimation(): #: exact solution tform = estimate_transform('similarity', SRC[:2, :], DST[:2, :]) assert_array_almost_equal(tform(SRC[:2, :]), DST[:2, :]) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) @@ -68,7 +67,6 @@ def test_affine_estimation(): #: exact solution tform = estimate_transform('affine', SRC[:3, :], DST[:3, :]) assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :]) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) #: over-determined tform = estimate_transform('affine', SRC, DST) @@ -91,10 +89,10 @@ def test_affine_implicit(): def test_projective(): #: exact solution tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) + assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :]) #: over-determined - tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) + tform = estimate_transform('projective', SRC, DST) assert_array_almost_equal(tform.inverse(tform(SRC)), SRC)