diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index b0790992..4721e7d1 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -3,6 +3,7 @@ from .radon_transform import * from .finite_radon_transform import * from ._project import homography as fast_homography from .integral import * -from .geometric import warp, estimate_transformation, geometric_transform, \ - SimilarityTransformation, AffineTransformation, ProjectiveTransformation, \ - PolynomialTransformation, swirl, homography +from ._geometric import (warp, estimate_transform, + SimilarityTransform, AffineTransform, + ProjectiveTransform, PolynomialTransform) +from ._warps import swirl, homography diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py new file mode 100644 index 00000000..0af7ffb1 --- /dev/null +++ b/skimage/transform/_geometric.py @@ -0,0 +1,563 @@ +import math +import numpy as np +from scipy import ndimage +from skimage.util import img_as_float + + +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 ``MxNx3`` or ``MxNx4`` arrays. + + """ + if a.ndim == 3: + a[:] = b[:, :, np.newaxis] + else: + a[:] = b + + +class GeometricTransform(object): + """Perform geometric transformations on a set of coordinates. + + Parameters + ---------- + matrix : 3x3 array, optional + Homogeneous transformation matrix. + + """ + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : Nx2 array + source coordinates + + Returns + ------- + coords : Nx2 array + transformed coordinates + + """ + raise NotImplementedError() + + + def inverse(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : Nx2 array + source coordinates + + Returns + ------- + coords : Nx2 array + transformed coordinates + + """ + raise NotImplementedError() + + + 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}`. 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 ]]. + + """ + + _coefs = range(8) + + def __init__(self, matrix=None): + self._matrix = matrix + + @property + def _inv_matrix(self): + return np.linalg.inv(self._matrix) + + 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._matrix) + + def inverse(self, coords): + return self._apply_mat(coords, self._inv_matrix) + + + def estimate(self, src, dst): + """Set the transformation matrix with the explicit transformation + parameters. + + Parameters + ---------- + src : Nx2 array + source coordinates + dst : Nx2 array + destination coordinates + + """ + 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, 8)) + 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 + + # Select relevant columns, depending on coeffs + A = A[:, self._coefs] + + b = np.hstack([xd, yd]) + + H = np.zeros((3, 3)) + H.flat[self._coefs] = np.linalg.lstsq(A, b)[0] + H[2, 2] = 1 + + self._matrix = H + + def __add__(self, other): + """Combine this transformation with another. + + """ + if isinstance(other, ProjectiveTransform): + return ProjectiveTransform(np.dot(other._matrix, self._matrix)) + 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 zoom factors in the x and y directions, + and the homogeneous transformation matrix is:: + + [[a0 a1 a2] + [b0 b1 b2] + [0 0 1]] + + Parameters + ---------- + scale : (sx, sy), floats + Scale factors. + rotation : float + Rotation angle in radians, counter-clockwise direction. + shear : float + Shear angle in radians, counter-clockwise direction. + translation : (tx, ty), floats + Translation in x and y. + + """ + + _coefs = range(6) + + def __init__(self, scale=None, rotation=None, shear=None, translation=None): + ProjectiveTransform.__init__(self) + + 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) + + a = rotation + sx, sy = scale + tx, ty = translation + + self._matrix = np.array([ + [sx * math.cos(a), - sy * math.sin(a + shear), tx], + [sx * math.sin(a), sy * math.cos(a + shear), ty], + [0, 0, 1] + ]) + + +class SimilarityTransform(AffineTransform): + """2D similarity transformation of the form:: + + X = a0*x + b0*y + a1 = + = m*x*cos(rotation) + m*y*sin(rotation) + a1 + + Y = b0*x + a0*y + b1 = + = m*x*sin(rotation) + m*y*cos(rotation) + b1 + + where ``m`` is a zoom factor and the homogeneous transformation matrix is:: + + [[a0 b0 a1] + [b0 a0 b1] + [0 0 1]] + + Parameters + ---------- + scale : float, optional + Scale / zoom factor. + rotation : float, optional + Rotation angle, counter-clockwise, in radians. + translation : (tx, ty) of float + x, y translation parameters + + """ + + def __init__(self, scale=None, rotation=None, translation=None): + if scale is not None: + scale = (scale, scale) + AffineTransform.__init__(self, scale=scale, + rotation=rotation, + shear=0, + translation=translation) + + +class PolynomialTransform(GeometricTransform): + """2D transformation of the form:: + + X = sum[j=0:n]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) + Y = sum[j=0:n]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) + + TODO: Describe structure of coefficients. + Shall we store it as a (2, M) ndarray? + + """ + + def __init__(self, coeffs=None): + """Create polynomial transformation. + + Parameters + ---------- + coeffs : array, optional + polynomial coefficients + + """ + self.coeffs = coeffs + + def estimate(self, src, dst, order): + """Set the transformation matrix with the explicit transformation + parameters. + + Parameters + ---------- + src : Nx2 array + source coordinates + dst : Nx2 array + destination coordinates + order : int + polynomial order (number of coefficients is order + 1) + + """ + xs = src[:, 0] + ys = src[:, 1] + xd = dst[:, 0] + yd = dst[:, 1] + rows = src.shape[0] + + # number of unknown polynomial coefficients + u = (order + 1) * (order + 2) + + A = np.zeros((rows * 2, u)) + pidx = 0 + for j in xrange(order + 1): + for i in xrange(j + 1): + A[:rows, pidx] = xs ** (j - i) * ys ** i + A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i + pidx += 1 + + b = np.hstack([xd, yd]) + + self.coeffs = np.linalg.lstsq(A, b)[0] + + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : Nx2 array + source coordinates + + Returns + ------- + coords : Nx2 array + transformed coordinates + + """ + x = coords[:, 0] + y = coords[:, 1] + u = len(self.coeffs.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 xrange(order + 1): + for i in xrange(j + 1): + dst[:, 0] += self.coeffs[pidx] * x ** (j - i) * y ** i + dst[:, 1] += self.coeffs[pidx + u / 2] * 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.') + + +TRANSFORMATIONS = { + 'similarity': SimilarityTransform, + 'affine': AffineTransform, + 'projective': ProjectiveTransform, + 'polynomial': PolynomialTransform, +} + + +def estimate_transform(ttype, src, dst, **kwargs): + """Estimate 2D geometric transformation parameters. + + You can determine the over-, well- and under-determined parameters + with the least-squares method. + + Number of source must match number of destination coordinates. + + Parameters + ---------- + ttype : {'similarity', 'affine', 'projective', 'polynomial'} + Type of transform. + kwargs : array or int + Function parameters (src, dst, n, angle):: + + NAME / TTYPE FUNCTION PARAMETERS + 'similarity' `src, `dst` + 'affine' `src, `dst` + 'projective' `src, `dst` + 'polynomial' `src, `dst`, `order` (polynomial order) + + 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) + + >>> tform.inverse(tform.forward(src)) # == src + + >>> # warp image using the estimated transformation + >>> from skimage import data + >>> image = data.camera() + + >>> warp(image, inverse_map=tform.inverse) + + >>> # create transformation with explicit parameters + >>> scale = 1.1 + >>> rotation = 1 + >>> translation = (10, 20) + >>> + >>> tform2 = tf.SimilarityTransform(scale, rotation, translation) + + >>> # unite transformations, applied in order from left to right + >>> tform3 = tform + tform2 + >>> tform3.forward(src) # == tform2.forward(tform.forward(src)) + + """ + ttype = ttype.lower() + if ttype not in TRANSFORMATIONS: + raise ValueError('the transformation type \'%s\' is not' + 'implemented' % ttype) + + tform = TRANSFORMATIONS[ttype]() + tform.estimate(src, dst, **kwargs) + + return tform + + +def matrix_transform(coords, matrix): + """Apply 2D matrix transform. + + Parameters + ---------- + coords : Nx2 array + x, y coordinates to transform + matrix : 3x3 array + Homogeneous transformation matrix. + + Returns + ------- + coords : Nx2 array + transformed coordinates + + """ + return ProjectiveTransform(matrix)(coords) + + +def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, + mode='constant', cval=0., reverse_map=None): + """Warp an image according to a given coordinate transformation. + + Parameters + ---------- + image : 2-D array + Input image. + inverse_map : transformation object, callable xy = f(xy, **kwargs) + Inverse coordinate map. A function that transforms a Px2 array of + ``(x, y)`` coordinates in the *output image* into their corresponding + coordinates in the *source image*. In case of a transformation object + its `inverse` method will be used as transformation function. Also see + examples below. + map_args : dict, optional + Keyword arguments passed to `inverse_map`. + output_shape : tuple (rows, cols) + Shape of the output image generated. + order : int + Order of splines used in interpolation. See + `scipy.ndimage.map_coordinates` for detail. + mode : string + How to handle values outside the image borders. See + `scipy.ndimage.map_coordinates` for detail. + cval : string + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + Examples + -------- + Shift an image to the right: + + >>> from skimage import data + >>> image = data.camera() + >>> + >>> def shift_right(xy): + ... xy[:, 0] -= 10 + ... return xy + >>> + >>> warp(image, shift_right) + + """ + # Backward API compatibility + if reverse_map is not None: + inverse_map = reverse_map + + if image.ndim < 2: + raise ValueError("Input must have more than 1 dimension.") + + image = np.atleast_3d(img_as_float(image)) + ishape = np.array(image.shape) + bands = ishape[2] + + if output_shape is None: + output_shape = ishape + + coords = np.empty(np.r_[3, output_shape], dtype=float) + + ## Construct transformed coordinates + + rows, cols = output_shape[:2] + + # Reshape grid coordinates into a (P, 2) array of (x, y) pairs + tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T + + # Map each (x, y) pair to the source image according to + # the user-provided mapping + if callable(getattr(inverse_map, 'inverse', None)): + inverse_map = inverse_map.inverse + tf_coords = inverse_map(tf_coords, **map_args) + + # 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, ...]) + + # colour-coordinate mapping + coords[2, ...] = range(bands) + + # Prefilter not necessary for order 1 interpolation + prefilter = order > 1 + mapped = ndimage.map_coordinates(image, coords, prefilter=prefilter, + mode=mode, order=order, cval=cval) + + # The spline filters sometimes return results outside [0, 1], + # so clip to ensure valid data + return np.clip(mapped.squeeze(), 0, 1) + diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py new file mode 100644 index 00000000..347b8440 --- /dev/null +++ b/skimage/transform/_warps.py @@ -0,0 +1,160 @@ +from ._geometric import warp, ProjectiveTransform +import numpy as np + +def _swirl_mapping(xy, center, rotation, strength, radius): + x, y = xy.T + x0, y0 = center + rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) + + # Ensure that the transformation decays to approximately 1/1000-th + # within the specified radius. + radius = radius / 5 * np.log(2) + + theta = rotation + strength * \ + np.exp(-rho / radius) + \ + np.arctan2(y - y0, x - x0) + + xy[..., 0] = x0 + rho * np.cos(theta) + xy[..., 1] = y0 + rho * np.sin(theta) + + return xy + + +def swirl(image, center=None, strength=1, radius=100, rotation=0, + output_shape=None, order=1, mode='constant', cval=0): + """Perform a swirl transformation. + + Parameters + ---------- + image : ndarray + Input image. + center : (x,y) tuple or (2,) ndarray + Center coordinate of transformation. + strength : float + The amount of swirling applied. + radius : float + The extent of the swirl in pixels. The effect dies out + rapidly beyond `radius`. + rotation : float + Additional rotation applied to the image. + + Returns + ------- + swirled : ndarray + Swirled version of the input. + + Other parameters + ---------------- + output_shape : tuple or ndarray + Size of the generated output image. + order : int + Order of splines used in interpolation. See + `scipy.ndimage.map_coordinates` for detail. + mode : string + How to handle values outside the image borders. See + `scipy.ndimage.map_coordinates` for detail. + cval : string + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + """ + + if center is None: + center = np.array(image.shape)[:2] / 2 + + warp_args = {'center': center, + 'rotation': rotation, + 'strength': strength, + 'radius': radius} + + return warp(image, _swirl_mapping, map_args=warp_args, + output_shape=output_shape, + order=order, mode=mode, cval=cval) + + +def homography(image, H, output_shape=None, order=1, + mode='constant', cval=0.): + """Perform a projective transformation (homography) on an image. + + For each pixel, given its 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}`. + 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 + ---------- + image : 2-D array + Input image. + H : array of shape ``(3, 3)`` + Transformation matrix H that defines the homography. + output_shape : tuple (rows, cols) + Shape of the output image generated. + order : int + Order of splines used in interpolation. + mode : string + How to handle values outside the image borders. Passed as-is + to ndimage. + cval : string + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + Examples + -------- + >>> # rotate by 90 degrees around origin and shift down by 2 + >>> x = np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1 + >>> x + array([[1, 2, 3], + [4, 5, 6], + [7, 8, 9]], dtype=uint8) + >>> theta = -np.pi/2 + >>> M = np.array([[np.cos(theta),-np.sin(theta),0], + ... [np.sin(theta), np.cos(theta),2], + ... [0, 0, 1]]) + >>> x90 = homography(x, M, order=1) + >>> x90 + array([[3, 6, 9], + [2, 5, 8], + [1, 4, 7]], dtype=uint8) + >>> # translate right by 2 and down by 1 + >>> y = np.zeros((5,5), dtype=np.uint8) + >>> y[1, 1] = 255 + >>> y + array([[ 0, 0, 0, 0, 0], + [ 0, 255, 0, 0, 0], + [ 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0]], dtype=uint8) + >>> M = np.array([[ 1., 0., 2.], + ... [ 0., 1., 1.], + ... [ 0., 0., 1.]]) + >>> y21 = homography(y, M, order=1) + >>> y21 + array([[ 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0], + [ 0, 0, 0, 255, 0], + [ 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0]], dtype=uint8) + + """ + import warnings + warnings.warn('the homography function is deprecated; ' + 'use the `warp` and `tform` function instead', + category=DeprecationWarning) + + tform = ProjectiveTransform(H) + return warp(image, inverse_map=tform.inverse, output_shape=output_shape, + order=order, mode=mode, cval=cval) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py deleted file mode 100644 index 2051e241..00000000 --- a/skimage/transform/geometric.py +++ /dev/null @@ -1,764 +0,0 @@ -# coding: utf-8 -import math -import numpy as np -from scipy import ndimage -from skimage.util import img_as_float - - -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 ``MxNx3`` or ``MxNx4`` arrays. - - """ - if a.ndim == 3: - a[:] = b[:, :, np.newaxis] - else: - a[:] = b - - -def geometric_transform(coords, matrix): - """Apply 2D geometric transformation. - - Parameters - ---------- - ttype : Nx2 array - x, y coordinates to transform - matrix : 3x3 array - homogeneous transformation matrix - - Returns - ------- - coords : Nx2 array - transformed coordinates - """ - coords = np.asarray(coords) - shape = coords.shape - if shape == (2,): - coords = np.array([coords]) - - 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] - - if shape == (2,): - return dst[0, :2] - else: - return dst[:, :2] - - -class GeometricTransformation(object): - - def __init__(self, matrix=None): - """Create geometric transformation which contains the transformation - parameters and can perform forward and reverse transformations. - - Parameters - ---------- - matrix : 3x3 array, optional - homogeneous transformation matrix - - """ - self.matrix = matrix - self.inverse_matrix = None - - def forward(self, coords): - """Apply forward transformation. - - Parameters - ---------- - coords : Nx2 array - source coordinates - - Returns - ------- - coords : Nx2 array - transformed coordinates - - """ - if self.matrix is None: - raise Exception('Transformation matrix must be set or estimated.') - return geometric_transform(coords, self.matrix) - - def reverse(self, coords): - """Apply reverse transformation. - - Parameters - ---------- - coords : Nx2 array - source coordinates - - Returns - ------- - coords : Nx2 array - transformed coordinates - - """ - if self.matrix is None: - raise Exception('Transformation matrix must be set or estimated.') - if self.inverse_matrix is None: - self.inverse_matrix = np.linalg.inv(self.matrix) - return geometric_transform(coords, self.inverse_matrix) - - def __add__(self, other): - if type(self) == type(other): - transformation = self.__class__ - else: - transformation = GeometricTransformation - return transformation(other.matrix.dot(self.matrix)) - - -class SimilarityTransformation(GeometricTransformation): - - """2D similarity transformation of the following form: - X = a0*x - b0*y + a1 = - = m*x*cos(rotation) - m*y*sin(rotation) + a1 - Y = b0*x + a0*y + b1 = - = m*x*sin(rotation) + m*y*cos(rotation) + b1 - where the homogeneous transformation matrix is: - [[a0 -b0 a1] - [b0 a0 b1] - [0 0 1]] - - """ - - def estimate(self, src, dst): - """Set the transformation matrix with the estimated parameters of the - given control points. - - Parameters - ---------- - src : Nx2 array - source coordinates - dst : Nx2 array - destination coordinates - - """ - xs = src[:, 0] - ys = src[:, 1] - xd = dst[:, 0] - yd = dst[:, 1] - rows = src.shape[0] - - #: params: a0, a1, b0, b1 - A = np.zeros((rows * 2, 4)) - A[:rows, 0] = xs - A[:rows, 2] = - ys - A[:rows, 1] = 1 - A[rows:, 2] = xs - A[rows:, 0] = ys - A[rows:, 3] = 1 - - b = np.hstack([xd, yd]) - - a0, a1, b0, b1 = np.linalg.lstsq(A, b)[0] - self.matrix = np.array([[a0, -b0, a1], - [b0, a0, b1], - [ 0, 0, 1]]) - - def from_params(self, scale, rotation, translation): - """Set the transformation matrix with the explicit transformation - parameters. - - Parameters - ---------- - scale : float - scale factor - rotation : float - rotation angle in counter-clockwise direction - translation : (tx, ty) as array, list or tuple - x, y translation parameters - - """ - self.matrix = np.array([ - [math.cos(rotation), - math.sin(rotation), 0], - [math.sin(rotation), math.cos(rotation), 0], - [ 0, 0, 1] - ]) - self.matrix *= scale - self.matrix[0:2, 2] = translation - - @property - def scale(self): - return self.matrix[0, 0] / math.cos(self.rotation) - - @property - def rotation(self): - return math.atan2(self.matrix[1, 0], self.matrix[1, 1]) - - @property - def translation(self): - return self.matrix[0:2, 2] - - -class AffineTransformation(GeometricTransformation): - - """2D affine transformation of the following 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 the homogeneous transformation matrix is: - [[a0 a1 a2] - [b0 b1 b2] - [0 0 1]] - - """ - - def estimate(self, src, dst): - """Set the transformation matrix with the estimated parameters of the - given control points. - - Parameters - ---------- - src : Nx2 array - source coordinates - dst : Nx2 array - destination coordinates - - """ - xs = src[:, 0] - ys = src[:, 1] - xd = dst[:, 0] - yd = dst[:, 1] - rows = src.shape[0] - - #: params: a0, a1, a2, b0, b1, b2 - A = np.zeros((rows * 2, 6)) - A[:rows, 0] = xs - A[:rows, 1] = ys - A[:rows, 2] = 1 - A[rows:, 3] = xs - A[rows:, 4] = ys - A[rows:, 5] = 1 - - b = np.hstack([xd, yd]) - - a0, a1, a2, b0, b1, b2 = np.linalg.lstsq(A, b)[0] - self.matrix = np.array([[a0, a1, a2], - [b0, b1, b2], - [0, 0, 1]]) - - def from_params(self, scale, rotation, shear, translation): - """Set the transformation matrix with the explicit transformation - parameters. - - Parameters - ---------- - scale : (sx, sy) as array, list or tuple - scale factors - rotation : float - rotation angle in counter-clockwise direction - shear : float - shear angle in counter-clockwise direction - translation : (tx, ty) as array, list or tuple - translation parameters - - """ - sx, sy = scale - self.matrix = 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.matrix[0:2, 2] = translation - - @property - def scale(self): - sx = math.sqrt(self.matrix[0, 0] ** 2 + self.matrix[1, 0] ** 2) - sy = math.sqrt(self.matrix[0, 1] ** 2 + self.matrix[1, 1] ** 2) - return sx, sy - - @property - def rotation(self): - return math.atan2(self.matrix[1, 0], self.matrix[0, 0]) - - @property - def shear(self): - beta = math.atan2(- self.matrix[0, 1], self.matrix[1, 1]) - return beta - self.rotation - - @property - def translation(self): - return self.matrix[0:2, 2] - - -class ProjectiveTransformation(GeometricTransformation): - - """2D projective transformation of the following form - X = (a0 + a1*x + a2*y) / (c0*x + c1*y + 1) - Y = (b0 + b1*x + b2*y) / (c0*x + c1*y + 1) - where the homogeneous transformation matrix is: - [[a0 a1 a2] - [b0 b1 b2] - [c0 c1 1]] - - """ - - def estimate(self, src, dst): - """Set the transformation matrix with the explicit transformation - parameters. - - Parameters - ---------- - src : Nx2 array - source coordinates - dst : Nx2 array - destination coordinates - - """ - 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, 8)) - 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 - - b = np.hstack([xd, yd]) - - a0, a1, a2, b0, b1, b2, c0, c1 = np.linalg.lstsq(A, b)[0] - self.matrix = np.array([[a0, a1, a2], - [b0, b1, b2], - [c0, c1, 1]]) - - -class PolynomialTransformation(GeometricTransformation): - - """2D affine transformation of the following form - X = sum[j=0:n]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) - Y = sum[j=0:n]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) - - """ - - def __init__(self, coeffs=None): - """Create polynomial transformation which contains the transformation - parameters and can perform forward and reverse transformations. - - Parameters - ---------- - coeffs : array, optional - polynomial coefficients - - """ - self.coeffs = coeffs - - def estimate(self, src, dst, order): - """Set the transformation matrix with the explicit transformation - parameters. - - Parameters - ---------- - src : Nx2 array - source coordinates - dst : Nx2 array - destination coordinates - order : int - polynomial order (number of coefficients is order + 1) - - """ - xs = src[:, 0] - ys = src[:, 1] - xd = dst[:, 0] - yd = dst[:, 1] - rows = src.shape[0] - - # number of unknown polynomial coefficients - u = (order + 1) * (order + 2) - - A = np.zeros((rows * 2, u)) - pidx = 0 - for j in xrange(order + 1): - for i in xrange(j + 1): - A[:rows, pidx] = xs ** (j - i) * ys ** i - A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i - pidx += 1 - - b = np.hstack([xd, yd]) - - self.coeffs = np.linalg.lstsq(A, b)[0] - - def forward(self, coords): - """Apply forward transformation. - - Parameters - ---------- - coords : Nx2 array - source coordinates - - Returns - ------- - coords : Nx2 array - transformed coordinates - - """ - x = coords[:, 0] - y = coords[:, 1] - u = len(self.coeffs) - # 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 xrange(order + 1): - for i in xrange(j + 1): - dst[:, 0] += self.coeffs[pidx] * x ** (j - i) * y ** i - dst[:, 1] += self.coeffs[pidx + u / 2] * x ** (j - i) * y ** i - pidx += 1 - - return dst - - def reverse(self, coords): - raise Exception( - 'There is no explicit way to do the reverse polynomial ' - 'transformation. Instead determine the reverse transformation ' - 'parameters by exchanging source and destination coordinates.' - 'Then apply the forward transformation.') - - -TRANSFORMATIONS = { - 'similarity': SimilarityTransformation, - 'affine': AffineTransformation, - 'projective': ProjectiveTransformation, - 'polynomial': PolynomialTransformation, -} - - -def estimate_transformation(ttype, src, dst, order=None): - """Estimate 2D geometric transformation parameters. - - You can determine the over-, well- and under-determined parameters - with the least-squares method. - - Number of source must match number of destination coordinates. - - Parameters - ---------- - ttype : str - one of similarity, affine, projective, polynomial - kwargs :: array or int - function parameters (src, dst, n, angle): - - NAME / TTYPE FUNCTION PARAMETERS - 'similarity' `src, `dst` - 'affine' `src, `dst` - 'projective' `src, `dst` - 'polynomial' `src, `dst`, `order` (polynomial order) - - See examples section below for usage. - - Returns - ------- - tform : subclass of :class:`GeometricTransformation` - tform object containing the transformation parameters and providing - access to forward and reverse 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_transformation('similarity', src, dst) - >>> tform.matrix - >>> tform.reverse(tform.forward(src)) # == src - >>> # warp image using the estimated transformation - >>> from skimage import data - >>> image = data.camera() - >>> tf.warp(image, tform) # == warp(image, reverse_map=tform.reverse) - >>> tf.warp(image, reverse_map=tform.forward) - >>> # create transformation with explicit parameters - >>> tform2 = tf.SimilarityTransformation() - >>> scale = 1.1 - >>> rotation = 1 - >>> translation = (10, 20) - >>> tform2.from_params(scale, rotation, translation) - >>> # unite transformations, applied in order from left to right - >>> tform3 = tform + tform2 - >>> tform3.forward(src) # == tform2.forward(tform.forward(src)) - - """ - ttype = ttype.lower() - if ttype not in TRANSFORMATIONS: - raise ValueError('the transformation type \'%s\' is not' - 'implemented' % ttype) - args = [src, dst] - if order is not None and ttype == 'polynomial': - args.append(order) - tform = TRANSFORMATIONS[ttype]() - tform.estimate(*args) - return tform - - -def warp(image, reverse_map=None, map_args={}, output_shape=None, order=1, - mode='constant', cval=0.): - """Warp an image according to a given coordinate transformation. - - Parameters - ---------- - image : 2-D array - Input image. - reverse_map : transformation object, callable xy = f(xy, **kwargs) - Reverse coordinate map. A function that transforms a Px2 array of - ``(x, y)`` coordinates in the *output image* into their corresponding - coordinates in the *source image*. In case of a transformation object - its `reverse` method will be used as transformation function. Also see - examples below. - map_args : dict, optional - Keyword arguments passed to `reverse_map`. - output_shape : tuple (rows, cols) - Shape of the output image generated. - order : int - Order of splines used in interpolation. See - `scipy.ndimage.map_coordinates` for detail. - mode : string - How to handle values outside the image borders. See - `scipy.ndimage.map_coordinates` for detail. - cval : string - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - Examples - -------- - Shift an image to the right: - - >>> from skimage import data - >>> image = data.camera() - >>> - >>> def shift_right(xy): - ... xy[:, 0] -= 10 - ... return xy - >>> - >>> warp(image, shift_right) - - """ - if image.ndim < 2: - raise ValueError("Input must have more than 1 dimension.") - - image = np.atleast_3d(img_as_float(image)) - ishape = np.array(image.shape) - bands = ishape[2] - - if output_shape is None: - output_shape = ishape - - coords = np.empty(np.r_[3, output_shape], dtype=float) - - ## Construct transformed coordinates - - rows, cols = output_shape[:2] - - # Reshape grid coordinates into a (P, 2) array of (x, y) pairs - tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T - - # Map each (x, y) pair to the source image according to - # the user-provided mapping - if callable(getattr(reverse_map, 'reverse', None)): - reverse_map = reverse_map.reverse - tf_coords = reverse_map(tf_coords, **map_args) - - # 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, ...]) - - # colour-coordinate mapping - coords[2, ...] = range(bands) - - # Prefilter not necessary for order 1 interpolation - prefilter = order > 1 - mapped = ndimage.map_coordinates(image, coords, prefilter=prefilter, - mode=mode, order=order, cval=cval) - - # The spline filters sometimes return results outside [0, 1], - # so clip to ensure valid data - return np.clip(mapped.squeeze(), 0, 1) - - -def _swirl_mapping(xy, center, rotation, strength, radius): - x, y = xy.T - x0, y0 = center - rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) - - # Ensure that the transformation decays to approximately 1/1000-th - # within the specified radius. - radius = radius / 5 * np.log(2) - - theta = rotation + strength * \ - np.exp(-rho / radius) + \ - np.arctan2(y - y0, x - x0) - - xy[..., 0] = x0 + rho * np.cos(theta) - xy[..., 1] = y0 + rho * np.sin(theta) - - return xy - - -def swirl(image, center=None, strength=1, radius=100, rotation=0, - output_shape=None, order=1, mode='constant', cval=0): - """Perform a swirl transformation. - - Parameters - ---------- - image : ndarray - Input image. - center : (x,y) tuple or (2,) ndarray - Center coordinate of transformation. - strength : float - The amount of swirling applied. - radius : float - The extent of the swirl in pixels. The effect dies out - rapidly beyond `radius`. - rotation : float - Additional rotation applied to the image. - - Returns - ------- - swirled : ndarray - Swirled version of the input. - - Other parameters - ---------------- - output_shape : tuple or ndarray - Size of the generated output image. - order : int - Order of splines used in interpolation. See - `scipy.ndimage.map_coordinates` for detail. - mode : string - How to handle values outside the image borders. See - `scipy.ndimage.map_coordinates` for detail. - cval : string - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - """ - - if center is None: - center = np.array(image.shape)[:2] / 2 - - warp_args = {'center': center, - 'rotation': rotation, - 'strength': strength, - 'radius': radius} - - return warp(image, _swirl_mapping, map_args=warp_args, - output_shape=output_shape, - order=order, mode=mode, cval=cval) - - -def homography(image, H, output_shape=None, order=1, - mode='constant', cval=0.): - """Perform a projective transformation (homography) on an image. - - For each pixel, given its 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}`. - 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 - ---------- - image : 2-D array - Input image. - H : array of shape ``(3, 3)`` - Transformation matrix H that defines the homography. - output_shape : tuple (rows, cols) - Shape of the output image generated. - order : int - Order of splines used in interpolation. - mode : string - How to handle values outside the image borders. Passed as-is - to ndimage. - cval : string - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - Examples - -------- - >>> # rotate by 90 degrees around origin and shift down by 2 - >>> x = np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1 - >>> x - array([[1, 2, 3], - [4, 5, 6], - [7, 8, 9]], dtype=uint8) - >>> theta = -np.pi/2 - >>> M = np.array([[np.cos(theta),-np.sin(theta),0], - ... [np.sin(theta), np.cos(theta),2], - ... [0, 0, 1]]) - >>> x90 = homography(x, M, order=1) - >>> x90 - array([[3, 6, 9], - [2, 5, 8], - [1, 4, 7]], dtype=uint8) - >>> # translate right by 2 and down by 1 - >>> y = np.zeros((5,5), dtype=np.uint8) - >>> y[1, 1] = 255 - >>> y - array([[ 0, 0, 0, 0, 0], - [ 0, 255, 0, 0, 0], - [ 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0]], dtype=uint8) - >>> M = np.array([[ 1., 0., 2.], - ... [ 0., 1., 1.], - ... [ 0., 0., 1.]]) - >>> y21 = homography(y, M, order=1) - >>> y21 - array([[ 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0], - [ 0, 0, 0, 255, 0], - [ 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0]], dtype=uint8) - - """ - import warnings - warnings.warn('the homography function is deprecated; ' - 'use the `warp` and `tform` function instead', - category=DeprecationWarning) - - tform = ProjectiveTransformation(H) - return warp(image, reverse_map=tform.reverse, output_shape=output_shape, - order=order, mode=mode, cval=cval) diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 1648e6f1..4430673c 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -1,12 +1,10 @@ import numpy as np from numpy.testing import assert_array_almost_equal -from skimage.transform.geometric import _stackcopy -from skimage.transform import estimate_transformation, homography, warp, \ - fast_homography, SimilarityTransformation, AffineTransformation, \ - ProjectiveTransformation, PolynomialTransformation -from skimage import transform as tf, data, img_as_float -from skimage.color import rgb2gray +from skimage.transform._geometric import _stackcopy +from skimage.transform import (estimate_transform, SimilarityTransform, + AffineTransform, ProjectiveTransform, + PolynomialTransform) SRC = np.array([ @@ -42,171 +40,50 @@ def test_stackcopy(): def test_similarity_estimation(): #: exact solution - tform = estimate_transformation('similarity', SRC[:2, :], DST[:2, :]) - assert_array_almost_equal(tform.forward(SRC[:2, :]), DST[:2, :]) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + tform = estimate_transform('similarity', SRC[:2, :], DST[:2, :]) + assert_array_almost_equal(tform(SRC[:2, :]), DST[:2, :]) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) #: over-determined - tform = estimate_transformation('similarity', SRC, DST) - ref = np.array( - [[2.3632898110e+02, -5.5876792257e+00, 2.5331569391e+03], - [5.5876792257e+00, 2.3632898110e+02, 2.4358232635e+03], - [0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00]]) - assert_array_almost_equal(tform.matrix, ref) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + tform = estimate_transform('similarity', SRC, DST) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) -def test_similarity_explicit(): - tform = SimilarityTransformation() - scale = 0.1 - rotation = 1 - translation = (1, 1) - tform.from_params(scale, rotation, translation) - assert_array_almost_equal(tform.scale, scale) - assert_array_almost_equal(tform.rotation, rotation) - assert_array_almost_equal(tform.translation, translation) - def test_affine_estimation(): #: exact solution - tform = estimate_transformation('affine', SRC[:3, :], DST[:3, :]) - assert_array_almost_equal(tform.forward(SRC[:3, :]), DST[:3, :]) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + tform = estimate_transform('affine', SRC[:3, :], DST[:3, :]) + assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :]) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) #: over-determined - tform = estimate_transformation('affine', SRC, DST) - ref = np.array( - [[2.2573930047e+02, 7.1588596765e+00, 2.5126622012e+03], - [2.1234856855e+01, 2.4931019555e+02, 2.4143862183e+03], - [0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00]]) - assert_array_almost_equal(tform.matrix, ref) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) - - -def test_affine_explicit(): - tform = AffineTransformation() - scale = (0.1, 0.13) - rotation = 1 - shear = 0.1 - translation = (1, 1) - tform.from_params(scale, rotation, shear, 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) + tform = estimate_transform('affine', SRC, DST) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) def test_projective(): #: exact solution - tform = estimate_transformation('projective', SRC[:4, :], DST[:4, :]) - ref = np.array( - [[ 1.9466901291e+02, -1.1888183994e+01, 2.2832379309e+03], - [ -8.6910077540e+00, 2.2162069773e+02, 2.2211673699e+03], - [ -1.2695966735e-02, -9.6053624285e-03, 1.0000000000e+00]]) - assert_array_almost_equal(tform.matrix, ref, 6) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) #: over-determined - tform = estimate_transformation('projective', SRC[:4, :], DST[:4, :]) - ref = np.array( - [[ 1.9466901291e+02, -1.1888183994e+01, 2.2832379309e+03], - [ -8.6910077540e+00, 2.2162069773e+02, 2.2211673699e+03], - [ -1.2695966735e-02, -9.6053624285e-03, 1.0000000000e+00]]) - assert_array_almost_equal(tform.matrix, ref, 6) - assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) + assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) def test_polynomial(): - tform = estimate_transformation('polynomial', SRC, DST, order=10) - assert_array_almost_equal(tform.forward(SRC), DST, 6) + tform = estimate_transform('polynomial', SRC, DST, order=10) + assert_array_almost_equal(tform(SRC), DST, 6) def test_union(): - tform1 = SimilarityTransformation() - scale1 = 0.1 - rotation1 = 1 - translation1 = (0, 0) - tform1.from_params(scale1, rotation1, translation1) - - tform2 = SimilarityTransformation() - scale2 = 0.1 - rotation2 = 1 - translation2 = (0, 0) - tform2.from_params(scale2, rotation2, translation2) + tform1 = SimilarityTransform(1, 0.3) + tform2 = SimilarityTransform(1, 0.6) + tform3 = SimilarityTransform(1, 0.9) tform = tform1 + tform2 - assert_array_almost_equal(tform.scale, scale1 * scale2) - assert_array_almost_equal(tform.rotation, rotation1 + rotation2) - - -def test_warp(): - x = np.zeros((5, 5), dtype=np.uint8) - x[2, 2] = 255 - x = img_as_float(x) - theta = -np.pi/2 - tform = SimilarityTransformation() - tform.from_params(1, theta, (0, 4)) - - x90 = warp(x, tform, order=1) - assert_array_almost_equal(x90, np.rot90(x)) - - x90 = warp(x, tform.reverse, order=1) - assert_array_almost_equal(x90, np.rot90(x)) - - -def test_homography(): - x = np.zeros((5, 5), dtype=np.uint8) - x[1, 1] = 255 - x = img_as_float(x) - theta = -np.pi/2 - M = np.array([[np.cos(theta),-np.sin(theta),0], - [np.sin(theta), np.cos(theta),4], - [0, 0, 1]]) - x90 = homography(x, M, order=1) - assert_array_almost_equal(x90, np.rot90(x)) - - -def test_fast_homography(): - img = rgb2gray(data.lena()).astype(np.uint8) - img = img[:, :100] - - theta = np.deg2rad(30) - scale = 0.5 - tx, ty = 50, 50 - - H = np.eye(3) - S = scale * np.sin(theta) - C = scale * np.cos(theta) - - H[:2, :2] = [[C, -S], [S, C]] - H[:2, 2] = [tx, ty] - - for mode in ('constant', 'mirror', 'wrap'): - p0 = homography(img, H, mode=mode, order=1) - p1 = fast_homography(img, H, mode=mode) - p1 = np.round(p1) - - ## import matplotlib.pyplot as plt - ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4) - ## ax0.imshow(img) - ## ax1.imshow(p0, cmap=plt.cm.gray) - ## ax2.imshow(p1, cmap=plt.cm.gray) - ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) - ## plt.show() - - d = np.mean(np.abs(p0 - p1)) - assert d < 0.2 - - -def test_swirl(): - image = img_as_float(data.checkerboard()) - - swirl_params = {'radius': 80, 'rotation': 0, 'order': 2, 'mode': 'reflect'} - swirled = tf.swirl(image, strength=10, **swirl_params) - unswirled = tf.swirl(swirled, strength=-10, **swirl_params) - - assert np.mean(np.abs(image - unswirled)) < 0.01 + assert_array_almost_equal(tform._matrix, tform3._matrix) if __name__ == "__main__": diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py new file mode 100644 index 00000000..047b477b --- /dev/null +++ b/skimage/transform/tests/test_warps.py @@ -0,0 +1,79 @@ +from numpy.testing import assert_array_almost_equal, run_module_suite +import numpy as np + +from skimage.transform import (warp, homography, fast_homography, + SimilarityTransform) +from skimage import transform as tf, data, img_as_float +from skimage.color import rgb2gray + + +def test_warp(): + x = np.zeros((5, 5), dtype=np.uint8) + x[2, 2] = 255 + x = img_as_float(x) + theta = -np.pi/2 + tform = SimilarityTransform(1, theta, (0, 4)) + + x90 = warp(x, tform, order=1) + assert_array_almost_equal(x90, np.rot90(x)) + + x90 = warp(x, tform.inverse, order=1) + assert_array_almost_equal(x90, np.rot90(x)) + + +def test_homography(): + x = np.zeros((5, 5), dtype=np.uint8) + x[1, 1] = 255 + x = img_as_float(x) + theta = -np.pi/2 + M = np.array([[np.cos(theta),-np.sin(theta),0], + [np.sin(theta), np.cos(theta),4], + [0, 0, 1]]) + x90 = homography(x, M, order=1) + assert_array_almost_equal(x90, np.rot90(x)) + + +def test_fast_homography(): + img = rgb2gray(data.lena()).astype(np.uint8) + img = img[:, :100] + + theta = np.deg2rad(30) + scale = 0.5 + tx, ty = 50, 50 + + H = np.eye(3) + S = scale * np.sin(theta) + C = scale * np.cos(theta) + + H[:2, :2] = [[C, -S], [S, C]] + H[:2, 2] = [tx, ty] + + for mode in ('constant', 'mirror', 'wrap'): + p0 = homography(img, H, mode=mode, order=1) + p1 = fast_homography(img, H, mode=mode) + p1 = np.round(p1) + + ## import matplotlib.pyplot as plt + ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4) + ## ax0.imshow(img) + ## ax1.imshow(p0, cmap=plt.cm.gray) + ## ax2.imshow(p1, cmap=plt.cm.gray) + ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) + ## plt.show() + + d = np.mean(np.abs(p0 - p1)) + assert d < 0.2 + + +def test_swirl(): + image = img_as_float(data.checkerboard()) + + swirl_params = {'radius': 80, 'rotation': 0, 'order': 2, 'mode': 'reflect'} + swirled = tf.swirl(image, strength=10, **swirl_params) + unswirled = tf.swirl(swirled, strength=-10, **swirl_params) + + assert np.mean(np.abs(image - unswirled)) < 0.01 + + +if __name__ == "__main__": + run_module_suite()