wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-25 13:30:51 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2044 from ahojnnes/umeyama

Browse Source 
Improved SimilarityTransform, new EuclideanTransform class
This commit is contained in:
Juan Nunez-Iglesias
2016-06-21 18:19:15 -04:00
committed by GitHub
parent 378e14fe8e bd3de29a8d
commit 9bc03189f8
3 changed files with 255 additions and 98 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/transform/__init__.py
+2 -1
View File
@@ -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',
skimage/transform/_geometric.py
+194 -89
View File
@@ -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`
skimage/transform/tests/test_geometric.py
+59 -8
View File
@@ -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)))