mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-30 18:51:53 +08:00
Merge pull request #1319 from ahojnnes/geom-est
Improve conditioning of geometric transformations
This commit is contained in:
Stefan van der Walt
@@ -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():
|
||||
|
||||
@@ -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):
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -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()
|
||||
|
||||
Reference in New Issue
Block a user