mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-27 21:53:30 +08:00
Johannes Schönberger
1493 lines
46 KiB
Python
1493 lines
46 KiB
Python
from __future__ import division
|
|
import six
|
|
import math
|
|
import numpy as np
|
|
from scipy import spatial
|
|
from scipy import ndimage as ndi
|
|
|
|
from .._shared.utils import (get_bound_method_class, safe_as_int, warn)
|
|
from ..util import img_as_float
|
|
|
|
from ._warps_cy import _warp_fast
|
|
|
|
|
|
def _to_ndimage_mode(mode):
|
|
"""Convert from `numpy.pad` mode name to the corresponding ndimage mode."""
|
|
mode_translation_dict = dict(edge='nearest', symmetric='reflect',
|
|
reflect='mirror')
|
|
if mode in mode_translation_dict:
|
|
mode = mode_translation_dict[mode]
|
|
return mode
|
|
|
|
|
|
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
|
|
|
|
|
|
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):
|
|
"""Base class for geometric transformations.
|
|
|
|
"""
|
|
def __call__(self, coords):
|
|
"""Apply forward transformation.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
Source coordinates.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
raise NotImplementedError()
|
|
|
|
def inverse(self, coords):
|
|
"""Apply inverse transformation.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
Source coordinates.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
raise NotImplementedError()
|
|
|
|
def residuals(self, src, dst):
|
|
"""Determine residuals of transformed destination coordinates.
|
|
|
|
For each transformed source coordinate the euclidean distance to the
|
|
respective destination coordinate is determined.
|
|
|
|
Parameters
|
|
----------
|
|
src : (N, 2) array
|
|
Source coordinates.
|
|
dst : (N, 2) array
|
|
Destination coordinates.
|
|
|
|
Returns
|
|
-------
|
|
residuals : (N, ) array
|
|
Residual for coordinate.
|
|
|
|
"""
|
|
|
|
return np.sqrt(np.sum((self(src) - dst)**2, axis=1))
|
|
|
|
def __add__(self, other):
|
|
"""Combine this transformation with another.
|
|
|
|
"""
|
|
raise NotImplementedError()
|
|
|
|
|
|
class ProjectiveTransform(GeometricTransform):
|
|
"""Matrix transformation.
|
|
|
|
Apply a projective transformation (homography) on coordinates.
|
|
|
|
For each homogeneous coordinate :math:`\mathbf{x} = [x, y, 1]^T`, its
|
|
target position is calculated by multiplying with the given matrix,
|
|
:math:`H`, to give :math:`H \mathbf{x}`::
|
|
|
|
[[a0 a1 a2]
|
|
[b0 b1 b2]
|
|
[c0 c1 1 ]].
|
|
|
|
E.g., to rotate by theta degrees clockwise, the matrix should be::
|
|
|
|
[[cos(theta) -sin(theta) 0]
|
|
[sin(theta) cos(theta) 0]
|
|
[0 0 1]]
|
|
|
|
or, to translate x by 10 and y by 20::
|
|
|
|
[[1 0 10]
|
|
[0 1 20]
|
|
[0 0 1 ]].
|
|
|
|
Parameters
|
|
----------
|
|
matrix : (3, 3) array, optional
|
|
Homogeneous transformation matrix.
|
|
|
|
Attributes
|
|
----------
|
|
params : (3, 3) array
|
|
Homogeneous transformation matrix.
|
|
|
|
"""
|
|
|
|
_coeffs = range(8)
|
|
|
|
def __init__(self, matrix=None):
|
|
if matrix is None:
|
|
# default to an identity transform
|
|
matrix = np.eye(3)
|
|
if matrix.shape != (3, 3):
|
|
raise ValueError("invalid shape of transformation matrix")
|
|
self.params = matrix
|
|
|
|
@property
|
|
def _inv_matrix(self):
|
|
return np.linalg.inv(self.params)
|
|
|
|
def _apply_mat(self, coords, matrix):
|
|
coords = np.array(coords, copy=False, ndmin=2)
|
|
|
|
x, y = np.transpose(coords)
|
|
src = np.vstack((x, y, np.ones_like(x)))
|
|
dst = np.dot(src.transpose(), matrix.transpose())
|
|
|
|
# rescale to homogeneous coordinates
|
|
dst[:, 0] /= dst[:, 2]
|
|
dst[:, 1] /= dst[:, 2]
|
|
|
|
return dst[:, :2]
|
|
|
|
def __call__(self, coords):
|
|
return self._apply_mat(coords, self.params)
|
|
|
|
def inverse(self, coords):
|
|
"""Apply inverse transformation.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
Source coordinates.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
return self._apply_mat(coords, self._inv_matrix)
|
|
|
|
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.
|
|
|
|
The transformation is defined as::
|
|
|
|
X = (a0*x + a1*y + a2) / (c0*x + c1*y + 1)
|
|
Y = (b0*x + b1*y + b2) / (c0*x + c1*y + 1)
|
|
|
|
These equations can be transformed to the following form::
|
|
|
|
0 = a0*x + a1*y + a2 - c0*x*X - c1*y*X - X
|
|
0 = b0*x + b1*y + b2 - c0*x*Y - c1*y*Y - 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 y 1 0 0 0 -x*X -y*X -X]
|
|
[0 0 0 x y 1 -x*Y -y*Y -Y]
|
|
...
|
|
...
|
|
]
|
|
x.T = [a0 a1 a2 b0 b1 b2 c0 c1 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.
|
|
|
|
In case of the affine transformation the coefficients c0 and c1 are 0.
|
|
Thus the system of equations is::
|
|
|
|
A = [[x y 1 0 0 0 -X]
|
|
[0 0 0 x y 1 -Y]
|
|
...
|
|
...
|
|
]
|
|
x.T = [a0 a1 a2 b0 b1 b2 c3]
|
|
|
|
Parameters
|
|
----------
|
|
src : (N, 2) array
|
|
Source coordinates.
|
|
dst : (N, 2) array
|
|
Destination coordinates.
|
|
|
|
Returns
|
|
-------
|
|
success : bool
|
|
True, if model estimation succeeds.
|
|
|
|
"""
|
|
|
|
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, a2, b0, b1, b2, c0, c1
|
|
A = np.zeros((rows * 2, 9))
|
|
A[:rows, 0] = xs
|
|
A[:rows, 1] = ys
|
|
A[:rows, 2] = 1
|
|
A[:rows, 6] = - xd * xs
|
|
A[:rows, 7] = - xd * ys
|
|
A[rows:, 3] = xs
|
|
A[rows:, 4] = ys
|
|
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 params
|
|
A = A[:, list(self._coeffs) + [8]]
|
|
|
|
_, _, V = np.linalg.svd(A)
|
|
|
|
H = np.zeros((3, 3))
|
|
# solution is right singular vector that corresponds to smallest
|
|
# singular value
|
|
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
|
|
|
|
return True
|
|
|
|
def __add__(self, other):
|
|
"""Combine this transformation with another.
|
|
|
|
"""
|
|
if isinstance(other, ProjectiveTransform):
|
|
# combination of the same types result in a transformation of this
|
|
# type again, otherwise use general projective transformation
|
|
if type(self) == type(other):
|
|
tform = self.__class__
|
|
else:
|
|
tform = ProjectiveTransform
|
|
return tform(other.params.dot(self.params))
|
|
elif (hasattr(other, '__name__')
|
|
and other.__name__ == 'inverse'
|
|
and hasattr(get_bound_method_class(other), '_inv_matrix')):
|
|
return ProjectiveTransform(self._inv_matrix.dot(self.params))
|
|
else:
|
|
raise TypeError("Cannot combine transformations of differing "
|
|
"types.")
|
|
|
|
|
|
class AffineTransform(ProjectiveTransform):
|
|
"""2D affine transformation of the form:
|
|
|
|
X = a0*x + a1*y + a2 =
|
|
= sx*x*cos(rotation) - sy*y*sin(rotation + shear) + a2
|
|
|
|
Y = b0*x + b1*y + b2 =
|
|
= sx*x*sin(rotation) + sy*y*cos(rotation + shear) + b2
|
|
|
|
where ``sx`` and ``sy`` are scale factors in the x and y directions,
|
|
and the homogeneous transformation matrix is::
|
|
|
|
[[a0 a1 a2]
|
|
[b0 b1 b2]
|
|
[0 0 1]]
|
|
|
|
Parameters
|
|
----------
|
|
matrix : (3, 3) array, optional
|
|
Homogeneous transformation matrix.
|
|
scale : (sx, sy) as array, list or tuple, optional
|
|
Scale factors.
|
|
rotation : float, optional
|
|
Rotation angle in counter-clockwise direction as radians.
|
|
shear : float, optional
|
|
Shear angle in counter-clockwise direction as radians.
|
|
translation : (tx, ty) as array, list or tuple, optional
|
|
Translation parameters.
|
|
|
|
Attributes
|
|
----------
|
|
params : (3, 3) array
|
|
Homogeneous transformation matrix.
|
|
|
|
"""
|
|
|
|
_coeffs = range(6)
|
|
|
|
def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
|
|
translation=None):
|
|
params = any(param is not None
|
|
for param in (scale, rotation, shear, 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 scale is None:
|
|
scale = (1, 1)
|
|
if rotation is None:
|
|
rotation = 0
|
|
if shear is None:
|
|
shear = 0
|
|
if translation is None:
|
|
translation = (0, 0)
|
|
|
|
sx, sy = scale
|
|
self.params = np.array([
|
|
[sx * math.cos(rotation), -sy * math.sin(rotation + shear), 0],
|
|
[sx * math.sin(rotation), sy * math.cos(rotation + shear), 0],
|
|
[ 0, 0, 1]
|
|
])
|
|
self.params[0:2, 2] = translation
|
|
else:
|
|
# default to an identity transform
|
|
self.params = np.eye(3)
|
|
|
|
@property
|
|
def scale(self):
|
|
sx = math.sqrt(self.params[0, 0] ** 2 + self.params[1, 0] ** 2)
|
|
sy = math.sqrt(self.params[0, 1] ** 2 + self.params[1, 1] ** 2)
|
|
return sx, sy
|
|
|
|
@property
|
|
def rotation(self):
|
|
return math.atan2(self.params[1, 0], self.params[0, 0])
|
|
|
|
@property
|
|
def shear(self):
|
|
beta = math.atan2(- self.params[0, 1], self.params[1, 1])
|
|
return beta - self.rotation
|
|
|
|
@property
|
|
def translation(self):
|
|
return self.params[0:2, 2]
|
|
|
|
|
|
class PiecewiseAffineTransform(GeometricTransform):
|
|
"""2D piecewise affine transformation.
|
|
|
|
Control points are used to define the mapping. The transform is based on
|
|
a Delaunay triangulation of the points to form a mesh. Each triangle is
|
|
used to find a local affine transform.
|
|
|
|
Attributes
|
|
----------
|
|
affines : list of AffineTransform objects
|
|
Affine transformations for each triangle in the mesh.
|
|
inverse_affines : list of AffineTransform objects
|
|
Inverse affine transformations for each triangle in the mesh.
|
|
|
|
"""
|
|
|
|
def __init__(self):
|
|
self._tesselation = None
|
|
self._inverse_tesselation = None
|
|
self.affines = None
|
|
self.inverse_affines = None
|
|
|
|
def estimate(self, src, dst):
|
|
"""Estimate the transformation from a set of corresponding points.
|
|
|
|
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.
|
|
|
|
"""
|
|
|
|
# forward piecewise affine
|
|
# triangulate input positions into mesh
|
|
self._tesselation = spatial.Delaunay(src)
|
|
# find affine mapping from source positions to destination
|
|
self.affines = []
|
|
for tri in self._tesselation.vertices:
|
|
affine = AffineTransform()
|
|
affine.estimate(src[tri, :], dst[tri, :])
|
|
self.affines.append(affine)
|
|
|
|
# inverse piecewise affine
|
|
# triangulate input positions into mesh
|
|
self._inverse_tesselation = spatial.Delaunay(dst)
|
|
# find affine mapping from source positions to destination
|
|
self.inverse_affines = []
|
|
for tri in self._inverse_tesselation.vertices:
|
|
affine = AffineTransform()
|
|
affine.estimate(dst[tri, :], src[tri, :])
|
|
self.inverse_affines.append(affine)
|
|
|
|
return True
|
|
|
|
def __call__(self, coords):
|
|
"""Apply forward transformation.
|
|
|
|
Coordinates outside of the mesh will be set to `- 1`.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
Source coordinates.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
|
|
out = np.empty_like(coords, np.double)
|
|
|
|
# determine triangle index for each coordinate
|
|
simplex = self._tesselation.find_simplex(coords)
|
|
|
|
# coordinates outside of mesh
|
|
out[simplex == -1, :] = -1
|
|
|
|
for index in range(len(self._tesselation.vertices)):
|
|
# affine transform for triangle
|
|
affine = self.affines[index]
|
|
# all coordinates within triangle
|
|
index_mask = simplex == index
|
|
|
|
out[index_mask, :] = affine(coords[index_mask, :])
|
|
|
|
return out
|
|
|
|
def inverse(self, coords):
|
|
"""Apply inverse transformation.
|
|
|
|
Coordinates outside of the mesh will be set to `- 1`.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
Source coordinates.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
|
|
out = np.empty_like(coords, np.double)
|
|
|
|
# determine triangle index for each coordinate
|
|
simplex = self._inverse_tesselation.find_simplex(coords)
|
|
|
|
# coordinates outside of mesh
|
|
out[simplex == -1, :] = -1
|
|
|
|
for index in range(len(self._inverse_tesselation.vertices)):
|
|
# affine transform for triangle
|
|
affine = self.inverse_affines[index]
|
|
# all coordinates within triangle
|
|
index_mask = simplex == index
|
|
|
|
out[index_mask, :] = affine(coords[index_mask, :])
|
|
|
|
return out
|
|
|
|
|
|
class EuclideanTransform(ProjectiveTransform):
|
|
"""2D Euclidean transformation of the form:
|
|
|
|
X = a0 * x - b0 * y + a1 =
|
|
= x * cos(rotation) - y * sin(rotation) + a1
|
|
|
|
Y = b0 * x + a0 * y + b1 =
|
|
= x * sin(rotation) + y * cos(rotation) + b1
|
|
|
|
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
|
|
Homogeneous transformation matrix.
|
|
scale : float, optional
|
|
Scale factor.
|
|
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, scale=None, rotation=None,
|
|
translation=None):
|
|
params = any(param is not None
|
|
for param in (scale, 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 scale is None:
|
|
scale = 1
|
|
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, 0:2] *= scale
|
|
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, True)
|
|
|
|
return True
|
|
|
|
@property
|
|
def scale(self):
|
|
if abs(math.cos(self.rotation)) < np.spacing(1):
|
|
# sin(self.rotation) == 1
|
|
scale = self.params[1, 0]
|
|
else:
|
|
scale = self.params[0, 0] / math.cos(self.rotation)
|
|
return scale
|
|
|
|
|
|
class PolynomialTransform(GeometricTransform):
|
|
"""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 ))
|
|
|
|
Parameters
|
|
----------
|
|
params : (2, N) array, optional
|
|
Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So,
|
|
a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`.
|
|
|
|
Attributes
|
|
----------
|
|
params : (2, N) array
|
|
Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So,
|
|
a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`.
|
|
|
|
"""
|
|
|
|
def __init__(self, params=None):
|
|
if params is None:
|
|
# default to transformation which preserves original coordinates
|
|
params = np.array([[0, 1, 0], [0, 0, 1]])
|
|
if params.shape[0] != 2:
|
|
raise ValueError("invalid shape of transformation parameters")
|
|
self.params = params
|
|
|
|
def estimate(self, src, dst, order=2):
|
|
"""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 = 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 ))
|
|
|
|
These equations can be transformed to the following form::
|
|
|
|
0 = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) - X
|
|
0 = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) - 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 = [[1 x y x**2 x*y y**2 ... 0 ... 0 -X]
|
|
[0 ... 0 1 x y x**2 x*y y**2 -Y]
|
|
...
|
|
...
|
|
]
|
|
x.T = [a00 a10 a11 a20 a21 a22 ... ann
|
|
b00 b10 b11 b20 b21 b22 ... bnn 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
|
|
Source coordinates.
|
|
dst : (N, 2) array
|
|
Destination coordinates.
|
|
order : int, optional
|
|
Polynomial order (number of coefficients is order + 1).
|
|
|
|
Returns
|
|
-------
|
|
success : bool
|
|
True, if model estimation succeeds.
|
|
|
|
"""
|
|
xs = src[:, 0]
|
|
ys = src[:, 1]
|
|
xd = dst[:, 0]
|
|
yd = dst[:, 1]
|
|
rows = src.shape[0]
|
|
|
|
# number of unknown polynomial coefficients
|
|
order = safe_as_int(order)
|
|
u = (order + 1) * (order + 2)
|
|
|
|
A = np.zeros((rows * 2, u + 1))
|
|
pidx = 0
|
|
for j in range(order + 1):
|
|
for i in range(j + 1):
|
|
A[:rows, pidx] = xs ** (j - i) * ys ** i
|
|
A[rows:, pidx + u // 2] = xs ** (j - i) * ys ** i
|
|
pidx += 1
|
|
|
|
A[:rows, -1] = xd
|
|
A[rows:, -1] = yd
|
|
|
|
_, _, V = np.linalg.svd(A)
|
|
|
|
# solution is right singular vector that corresponds to smallest
|
|
# singular value
|
|
params = - V[-1, :-1] / V[-1, -1]
|
|
|
|
self.params = params.reshape((2, u // 2))
|
|
|
|
return True
|
|
|
|
def __call__(self, coords):
|
|
"""Apply forward transformation.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
source coordinates
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
x = coords[:, 0]
|
|
y = coords[:, 1]
|
|
u = len(self.params.ravel())
|
|
# number of coefficients -> u = (order + 1) * (order + 2)
|
|
order = int((- 3 + math.sqrt(9 - 4 * (2 - u))) / 2)
|
|
dst = np.zeros(coords.shape)
|
|
|
|
pidx = 0
|
|
for j in range(order + 1):
|
|
for i in range(j + 1):
|
|
dst[:, 0] += self.params[0, pidx] * x ** (j - i) * y ** i
|
|
dst[:, 1] += self.params[1, pidx] * x ** (j - i) * y ** i
|
|
pidx += 1
|
|
|
|
return dst
|
|
|
|
def inverse(self, coords):
|
|
raise Exception(
|
|
'There is no explicit way to do the inverse polynomial '
|
|
'transformation. Instead, estimate the inverse transformation '
|
|
'parameters by exchanging source and destination coordinates,'
|
|
'then apply the forward transformation.')
|
|
|
|
|
|
TRANSFORMS = {
|
|
'euclidean': EuclideanTransform,
|
|
'similarity': SimilarityTransform,
|
|
'affine': AffineTransform,
|
|
'piecewise-affine': PiecewiseAffineTransform,
|
|
'projective': ProjectiveTransform,
|
|
'polynomial': PolynomialTransform,
|
|
}
|
|
|
|
HOMOGRAPHY_TRANSFORMS = (
|
|
SimilarityTransform,
|
|
AffineTransform,
|
|
ProjectiveTransform
|
|
)
|
|
|
|
|
|
def estimate_transform(ttype, src, dst, **kwargs):
|
|
"""Estimate 2D geometric transformation parameters.
|
|
|
|
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
|
|
----------
|
|
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`
|
|
'projective' `src, `dst`
|
|
'polynomial' `src, `dst`, `order` (polynomial order,
|
|
default order is 2)
|
|
|
|
Also see examples below.
|
|
|
|
Returns
|
|
-------
|
|
tform : :class:`GeometricTransform`
|
|
Transform object containing the transformation parameters and providing
|
|
access to forward and inverse transformation functions.
|
|
|
|
Examples
|
|
--------
|
|
>>> import numpy as np
|
|
>>> from skimage import transform as tf
|
|
|
|
>>> # estimate transformation parameters
|
|
>>> src = np.array([0, 0, 10, 10]).reshape((2, 2))
|
|
>>> dst = np.array([12, 14, 1, -20]).reshape((2, 2))
|
|
|
|
>>> tform = tf.estimate_transform('similarity', src, dst)
|
|
|
|
>>> np.allclose(tform.inverse(tform(src)), src)
|
|
True
|
|
|
|
>>> # warp image using the estimated transformation
|
|
>>> from skimage import data
|
|
>>> image = data.camera()
|
|
|
|
>>> warp(image, inverse_map=tform.inverse) # doctest: +SKIP
|
|
|
|
>>> # create transformation with explicit parameters
|
|
>>> tform2 = tf.SimilarityTransform(scale=1.1, rotation=1,
|
|
... translation=(10, 20))
|
|
|
|
>>> # unite transformations, applied in order from left to right
|
|
>>> tform3 = tform + tform2
|
|
>>> np.allclose(tform3(src), tform2(tform(src)))
|
|
True
|
|
|
|
"""
|
|
ttype = ttype.lower()
|
|
if ttype not in TRANSFORMS:
|
|
raise ValueError('the transformation type \'%s\' is not'
|
|
'implemented' % ttype)
|
|
|
|
tform = TRANSFORMS[ttype]()
|
|
tform.estimate(src, dst, **kwargs)
|
|
|
|
return tform
|
|
|
|
|
|
def matrix_transform(coords, matrix):
|
|
"""Apply 2D matrix transform.
|
|
|
|
Parameters
|
|
----------
|
|
coords : (N, 2) array
|
|
x, y coordinates to transform
|
|
matrix : (3, 3) array
|
|
Homogeneous transformation matrix.
|
|
|
|
Returns
|
|
-------
|
|
coords : (N, 2) array
|
|
Transformed coordinates.
|
|
|
|
"""
|
|
return ProjectiveTransform(matrix)(coords)
|
|
|
|
|
|
def _stackcopy(a, b):
|
|
"""Copy b into each color layer of a, such that::
|
|
|
|
a[:,:,0] = a[:,:,1] = ... = b
|
|
|
|
Parameters
|
|
----------
|
|
a : (M, N) or (M, N, P) ndarray
|
|
Target array.
|
|
b : (M, N)
|
|
Source array.
|
|
|
|
Notes
|
|
-----
|
|
Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
|
|
|
|
"""
|
|
if a.ndim == 3:
|
|
a[:] = b[:, :, np.newaxis]
|
|
else:
|
|
a[:] = b
|
|
|
|
|
|
def warp_coords(coord_map, shape, dtype=np.float64):
|
|
"""Build the source coordinates for the output of a 2-D image warp.
|
|
|
|
Parameters
|
|
----------
|
|
coord_map : callable like GeometricTransform.inverse
|
|
Return input coordinates for given output coordinates.
|
|
Coordinates are in the shape (P, 2), where P is the number
|
|
of coordinates and each element is a ``(row, col)`` pair.
|
|
shape : tuple
|
|
Shape of output image ``(rows, cols[, bands])``.
|
|
dtype : np.dtype or string
|
|
dtype for return value (sane choices: float32 or float64).
|
|
|
|
Returns
|
|
-------
|
|
coords : (ndim, rows, cols[, bands]) array of dtype `dtype`
|
|
Coordinates for `scipy.ndimage.map_coordinates`, that will yield
|
|
an image of shape (orows, ocols, bands) by drawing from source
|
|
points according to the `coord_transform_fn`.
|
|
|
|
Notes
|
|
-----
|
|
|
|
This is a lower-level routine that produces the source coordinates for 2-D
|
|
images used by `warp()`.
|
|
|
|
It is provided separately from `warp` to give additional flexibility to
|
|
users who would like, for example, to re-use a particular coordinate
|
|
mapping, to use specific dtypes at various points along the the
|
|
image-warping process, or to implement different post-processing logic
|
|
than `warp` performs after the call to `ndi.map_coordinates`.
|
|
|
|
|
|
Examples
|
|
--------
|
|
Produce a coordinate map that shifts an image up and to the right:
|
|
|
|
>>> from skimage import data
|
|
>>> from scipy.ndimage import map_coordinates
|
|
>>>
|
|
>>> def shift_up10_left20(xy):
|
|
... return xy - np.array([-20, 10])[None, :]
|
|
>>>
|
|
>>> image = data.astronaut().astype(np.float32)
|
|
>>> coords = warp_coords(shift_up10_left20, image.shape)
|
|
>>> warped_image = map_coordinates(image, coords)
|
|
|
|
"""
|
|
shape = safe_as_int(shape)
|
|
rows, cols = shape[0], shape[1]
|
|
coords_shape = [len(shape), rows, cols]
|
|
if len(shape) == 3:
|
|
coords_shape.append(shape[2])
|
|
coords = np.empty(coords_shape, dtype=dtype)
|
|
|
|
# Reshape grid coordinates into a (P, 2) array of (row, col) pairs
|
|
tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T
|
|
|
|
# Map each (row, col) pair to the source image according to
|
|
# the user-provided mapping
|
|
tf_coords = coord_map(tf_coords)
|
|
|
|
# Reshape back to a (2, M, N) coordinate grid
|
|
tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
|
|
|
|
# Place the y-coordinate mapping
|
|
_stackcopy(coords[1, ...], tf_coords[0, ...])
|
|
|
|
# Place the x-coordinate mapping
|
|
_stackcopy(coords[0, ...], tf_coords[1, ...])
|
|
|
|
if len(shape) == 3:
|
|
coords[2, ...] = range(shape[2])
|
|
|
|
return coords
|
|
|
|
|
|
def _convert_warp_input(image, preserve_range):
|
|
"""Convert input image to double image with the appropriate range."""
|
|
if preserve_range:
|
|
image = image.astype(np.double)
|
|
else:
|
|
image = img_as_float(image)
|
|
return image
|
|
|
|
|
|
def _clip_warp_output(input_image, output_image, order, mode, cval, clip):
|
|
"""Clip output image to range of values of input image.
|
|
|
|
Note that this function modifies the values of `output_image` in-place
|
|
and it is only modified if ``clip=True``.
|
|
|
|
Parameters
|
|
----------
|
|
input_image : ndarray
|
|
Input image.
|
|
output_image : ndarray
|
|
Output image, which is modified in-place.
|
|
|
|
Other parameters
|
|
----------------
|
|
order : int, optional
|
|
The order of the spline interpolation, default is 1. The order has to
|
|
be in the range 0-5. See `skimage.transform.warp` for detail.
|
|
mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
|
|
Points outside the boundaries of the input are filled according
|
|
to the given mode. Modes match the behaviour of `numpy.pad`.
|
|
cval : float, optional
|
|
Used in conjunction with mode 'constant', the value outside
|
|
the image boundaries.
|
|
clip : bool, optional
|
|
Whether to clip the output to the range of values of the input image.
|
|
This is enabled by default, since higher order interpolation may
|
|
produce values outside the given input range.
|
|
|
|
"""
|
|
if clip and order != 0:
|
|
min_val = input_image.min()
|
|
max_val = input_image.max()
|
|
|
|
preserve_cval = mode == 'constant' and not \
|
|
(min_val <= cval <= max_val)
|
|
|
|
if preserve_cval:
|
|
cval_mask = output_image == cval
|
|
|
|
np.clip(output_image, min_val, max_val, out=output_image)
|
|
|
|
if preserve_cval:
|
|
output_image[cval_mask] = cval
|
|
|
|
|
|
def warp(image, inverse_map, map_args={}, output_shape=None, order=1,
|
|
mode='constant', cval=0., clip=True, preserve_range=False):
|
|
"""Warp an image according to a given coordinate transformation.
|
|
|
|
Parameters
|
|
----------
|
|
image : ndarray
|
|
Input image.
|
|
inverse_map : transformation object, callable ``cr = f(cr, **kwargs)``, or ndarray
|
|
Inverse coordinate map, which transforms coordinates in the output
|
|
images into their corresponding coordinates in the input image.
|
|
|
|
There are a number of different options to define this map, depending
|
|
on the dimensionality of the input image. A 2-D image can have 2
|
|
dimensions for gray-scale images, or 3 dimensions with color
|
|
information.
|
|
|
|
- For 2-D images, you can directly pass a transformation object,
|
|
e.g. `skimage.transform.SimilarityTransform`, or its inverse.
|
|
- For 2-D images, you can pass a ``(3, 3)`` homogeneous
|
|
transformation matrix, e.g.
|
|
`skimage.transform.SimilarityTransform.params`.
|
|
- For 2-D images, a function that transforms a ``(M, 2)`` array of
|
|
``(col, row)`` coordinates in the output image to their
|
|
corresponding coordinates in the input image. Extra parameters to
|
|
the function can be specified through `map_args`.
|
|
- For N-D images, you can directly pass an array of coordinates.
|
|
The first dimension specifies the coordinates in the input image,
|
|
while the subsequent dimensions determine the position in the
|
|
output image. E.g. in case of 2-D images, you need to pass an array
|
|
of shape ``(2, rows, cols)``, where `rows` and `cols` determine the
|
|
shape of the output image, and the first dimension contains the
|
|
``(row, col)`` coordinate in the input image.
|
|
See `scipy.ndimage.map_coordinates` for further documentation.
|
|
|
|
Note, that a ``(3, 3)`` matrix is interpreted as a homogeneous
|
|
transformation matrix, so you cannot interpolate values from a 3-D
|
|
input, if the output is of shape ``(3,)``.
|
|
|
|
See example section for usage.
|
|
map_args : dict, optional
|
|
Keyword arguments passed to `inverse_map`.
|
|
output_shape : tuple (rows, cols), optional
|
|
Shape of the output image generated. By default the shape of the input
|
|
image is preserved. Note that, even for multi-band images, only rows
|
|
and columns need to be specified.
|
|
order : int, optional
|
|
The order of interpolation. The order has to be in the range 0-5:
|
|
- 0: Nearest-neighbor
|
|
- 1: Bi-linear (default)
|
|
- 2: Bi-quadratic
|
|
- 3: Bi-cubic
|
|
- 4: Bi-quartic
|
|
- 5: Bi-quintic
|
|
mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
|
|
Points outside the boundaries of the input are filled according
|
|
to the given mode. Modes match the behaviour of `numpy.pad`.
|
|
cval : float, optional
|
|
Used in conjunction with mode 'constant', the value outside
|
|
the image boundaries.
|
|
clip : bool, optional
|
|
Whether to clip the output to the range of values of the input image.
|
|
This is enabled by default, since higher order interpolation may
|
|
produce values outside the given input range.
|
|
preserve_range : bool, optional
|
|
Whether to keep the original range of values. Otherwise, the input
|
|
image is converted according to the conventions of `img_as_float`.
|
|
|
|
Returns
|
|
-------
|
|
warped : double ndarray
|
|
The warped input image.
|
|
|
|
Notes
|
|
-----
|
|
- The input image is converted to a `double` image.
|
|
- In case of a `SimilarityTransform`, `AffineTransform` and
|
|
`ProjectiveTransform` and `order` in [0, 3] this function uses the
|
|
underlying transformation matrix to warp the image with a much faster
|
|
routine.
|
|
|
|
Examples
|
|
--------
|
|
>>> from skimage.transform import warp
|
|
>>> from skimage import data
|
|
>>> image = data.camera()
|
|
|
|
The following image warps are all equal but differ substantially in
|
|
execution time. The image is shifted to the bottom.
|
|
|
|
Use a geometric transform to warp an image (fast):
|
|
|
|
>>> from skimage.transform import SimilarityTransform
|
|
>>> tform = SimilarityTransform(translation=(0, -10))
|
|
>>> warped = warp(image, tform)
|
|
|
|
Use a callable (slow):
|
|
|
|
>>> def shift_down(xy):
|
|
... xy[:, 1] -= 10
|
|
... return xy
|
|
>>> warped = warp(image, shift_down)
|
|
|
|
Use a transformation matrix to warp an image (fast):
|
|
|
|
>>> matrix = np.array([[1, 0, 0], [0, 1, -10], [0, 0, 1]])
|
|
>>> warped = warp(image, matrix)
|
|
>>> from skimage.transform import ProjectiveTransform
|
|
>>> warped = warp(image, ProjectiveTransform(matrix=matrix))
|
|
|
|
You can also use the inverse of a geometric transformation (fast):
|
|
|
|
>>> warped = warp(image, tform.inverse)
|
|
|
|
For N-D images you can pass a coordinate array, that specifies the
|
|
coordinates in the input image for every element in the output image. E.g.
|
|
if you want to rescale a 3-D cube, you can do:
|
|
|
|
>>> cube_shape = np.array([30, 30, 30])
|
|
>>> cube = np.random.rand(*cube_shape)
|
|
|
|
Setup the coordinate array, that defines the scaling:
|
|
|
|
>>> scale = 0.1
|
|
>>> output_shape = (scale * cube_shape).astype(int)
|
|
>>> coords0, coords1, coords2 = np.mgrid[:output_shape[0],
|
|
... :output_shape[1], :output_shape[2]]
|
|
>>> coords = np.array([coords0, coords1, coords2])
|
|
|
|
Assume that the cube contains spatial data, where the first array element
|
|
center is at coordinate (0.5, 0.5, 0.5) in real space, i.e. we have to
|
|
account for this extra offset when scaling the image:
|
|
|
|
>>> coords = (coords + 0.5) / scale - 0.5
|
|
>>> warped = warp(cube, coords)
|
|
|
|
"""
|
|
image = _convert_warp_input(image, preserve_range)
|
|
|
|
input_shape = np.array(image.shape)
|
|
|
|
if output_shape is None:
|
|
output_shape = input_shape
|
|
else:
|
|
output_shape = safe_as_int(output_shape)
|
|
|
|
warped = None
|
|
|
|
if order == 2:
|
|
# When fixing this issue, make sure to fix the branches further
|
|
# below in this function
|
|
warn("Bi-quadratic interpolation behavior has changed due "
|
|
"to a bug in the implementation of scikit-image. "
|
|
"The new version now serves as a wrapper "
|
|
"around SciPy's interpolation functions, which itself "
|
|
"is not verified to be a correct implementation. Until "
|
|
"skimage's implementation is fixed, we recommend "
|
|
"to use bi-linear or bi-cubic interpolation instead.")
|
|
|
|
if order in (0, 1, 3) and not map_args:
|
|
# use fast Cython version for specific interpolation orders and input
|
|
|
|
matrix = None
|
|
|
|
if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3):
|
|
# inverse_map is a transformation matrix as numpy array
|
|
matrix = inverse_map
|
|
|
|
elif isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
|
|
# inverse_map is a homography
|
|
matrix = inverse_map.params
|
|
|
|
elif (hasattr(inverse_map, '__name__')
|
|
and inverse_map.__name__ == 'inverse'
|
|
and get_bound_method_class(inverse_map) \
|
|
in HOMOGRAPHY_TRANSFORMS):
|
|
# inverse_map is the inverse of a homography
|
|
matrix = np.linalg.inv(six.get_method_self(inverse_map).params)
|
|
|
|
if matrix is not None:
|
|
matrix = matrix.astype(np.double)
|
|
if image.ndim == 2:
|
|
warped = _warp_fast(image, matrix,
|
|
output_shape=output_shape,
|
|
order=order, mode=mode, cval=cval)
|
|
elif image.ndim == 3:
|
|
dims = []
|
|
for dim in range(image.shape[2]):
|
|
dims.append(_warp_fast(image[..., dim], matrix,
|
|
output_shape=output_shape,
|
|
order=order, mode=mode, cval=cval))
|
|
warped = np.dstack(dims)
|
|
|
|
if warped is None:
|
|
# use ndi.map_coordinates
|
|
|
|
if (isinstance(inverse_map, np.ndarray)
|
|
and inverse_map.shape == (3, 3)):
|
|
# inverse_map is a transformation matrix as numpy array,
|
|
# this is only used for order >= 4.
|
|
inverse_map = ProjectiveTransform(matrix=inverse_map)
|
|
|
|
if isinstance(inverse_map, np.ndarray):
|
|
# inverse_map is directly given as coordinates
|
|
coords = inverse_map
|
|
else:
|
|
# inverse_map is given as function, that transforms (N, 2)
|
|
# destination coordinates to their corresponding source
|
|
# coordinates. This is only supported for 2(+1)-D images.
|
|
|
|
if image.ndim < 2 or image.ndim > 3:
|
|
raise ValueError("Only 2-D images (grayscale or color) are "
|
|
"supported, when providing a callable "
|
|
"`inverse_map`.")
|
|
|
|
def coord_map(*args):
|
|
return inverse_map(*args, **map_args)
|
|
|
|
if len(input_shape) == 3 and len(output_shape) == 2:
|
|
# Input image is 2D and has color channel, but output_shape is
|
|
# given for 2-D images. Automatically add the color channel
|
|
# dimensionality.
|
|
output_shape = (output_shape[0], output_shape[1],
|
|
input_shape[2])
|
|
|
|
coords = warp_coords(coord_map, output_shape)
|
|
|
|
# Pre-filtering not necessary for order 0, 1 interpolation
|
|
prefilter = order > 1
|
|
|
|
ndi_mode = _to_ndimage_mode(mode)
|
|
warped = ndi.map_coordinates(image, coords, prefilter=prefilter,
|
|
mode=ndi_mode, order=order, cval=cval)
|
|
|
|
|
|
_clip_warp_output(image, warped, order, mode, cval, clip)
|
|
|
|
return warped
|