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