From 84e18de02da8f01907e30bb0ffe1ded7e006a095 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 4 Jun 2012 15:22:30 +0200 Subject: [PATCH 01/41] improve and restructure geomtric transformations --- skimage/transform/__init__.py | 3 +- skimage/transform/_warp.py | 113 ----- skimage/transform/_warp_zoo.py | 75 --- skimage/transform/geometric.py | 558 ++++++++++++++++++++++ skimage/transform/project.py | 137 ------ skimage/transform/tests/test_geometric.py | 146 ++++++ skimage/transform/tests/test_project.py | 63 --- skimage/transform/tests/test_swirl.py | 17 - 8 files changed, 705 insertions(+), 407 deletions(-) delete mode 100644 skimage/transform/_warp.py delete mode 100644 skimage/transform/_warp_zoo.py create mode 100644 skimage/transform/geometric.py delete mode 100644 skimage/transform/project.py create mode 100644 skimage/transform/tests/test_geometric.py delete mode 100644 skimage/transform/tests/test_project.py delete mode 100644 skimage/transform/tests/test_swirl.py diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index fa4059ee..f5621444 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -4,5 +4,4 @@ from .finite_radon_transform import * from .project import * from ._project import homography as fast_homography from .integral import * -from ._warp import warp -from ._warp_zoo import swirl +from .geometric import warp, make_tform, swirl, homography diff --git a/skimage/transform/_warp.py b/skimage/transform/_warp.py deleted file mode 100644 index 6477c988..00000000 --- a/skimage/transform/_warp.py +++ /dev/null @@ -1,113 +0,0 @@ -__all__ = ['warp'] - -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 warp(image, reverse_map, 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 : 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*. 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 - 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) diff --git a/skimage/transform/_warp_zoo.py b/skimage/transform/_warp_zoo.py deleted file mode 100644 index 4df531d1..00000000 --- a/skimage/transform/_warp_zoo.py +++ /dev/null @@ -1,75 +0,0 @@ -from __future__ import division -import numpy as np - -from ._warp import warp - - -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) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py new file mode 100644 index 00000000..0306af59 --- /dev/null +++ b/skimage/transform/geometric.py @@ -0,0 +1,558 @@ +# coding: utf-8 +import math +import numpy as np +from scipy import ndimage +from skimage.util import img_as_float + + +EPS = np.spacing(1) + + +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 _make_similarity(src, dst): + """Determine parameters of the 2D similarity transformation: + X = a0*x - b0*y + a1 + Y = b0*x + a0*y + a2 + where the homogeneous transformation matrix is: + [[a1 -b1 a0] + [b1 a1 b0] + [0 0 1]] + """ + xs = src[:,0] + ys = src[:,1] + xd = dst[:,0] + yd = dst[:,1] + rows = src.shape[0] + + A = np.zeros((rows*2, 4)) + b = np.zeros((rows*2,)) + + 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[:rows] = xd + b[rows:] = yd + + a0, a1, b0, b1 = np.linalg.lstsq(A, b)[0] + matrix = np.eye(3) + matrix[0,0] = a0 + matrix[0,1] = - b0 + matrix[0,2] = a1 + matrix[1,0] = b0 + matrix[1,1] = a0 + matrix[1,2] = b1 + return matrix + +def _make_affine(src, dst): + """Determine parameters of the 2D affine transformation: + X = a0*x + a1*y + a3 + Y = b0*x + b1*y + b3 + where the homogeneous transformation matrix is: + [[a0 a1 a2] + [b0 b1 b2] + [0 0 1]] + """ + xs = src[:,0] + ys = src[:,1] + xd = dst[:,0] + yd = dst[:,1] + rows = src.shape[0] + + A = np.zeros((rows*2, 6)) + b = np.zeros((rows*2,)) + + 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[:rows] = xd + b[rows:] = yd + + params = np.linalg.lstsq(A, b)[0] + matrix = np.eye(3) + matrix[:2,:] = params.reshape((2, 3)) + return matrix + +def _make_projective(src, dst): + """Determine transformation matrix of the 2D projective transformation: + X = (a0 + a1*x + a2*y) / (c0*x + c1*y + c3) + Y = (b0 + b1*x + b2*y) / (c0*x + c1*y + c3) + where the homogeneous transformation matrix is: + [[a0 a1 a2] + [b0 b1 b2] + [c0 c1 c3]] + """ + xs = src[:,0] + ys = src[:,1] + xd = dst[:,0] + yd = dst[:,1] + rows = src.shape[0] + + A = np.zeros((rows*2, 8)) + b = np.zeros((rows*2,)) + + + 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[:rows] = dst[:,0] + b[rows:] = dst[:,1] + + matrix = np.eye(3).flatten() + matrix[:8] = np.linalg.lstsq(A, b)[0] + return matrix.reshape((3, 3)) + +def _make_polynomial(src, dst, order): + """Determine parameters of 2D polynomial transformation of order n: + 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 )) + """ + 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)) + b = np.zeros((rows*2,)) + + 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[:rows] = xd + b[rows:] = yd + + return np.linalg.lstsq(A, b)[0] + +def _make_rotation(angle): + """Determine homogeneous transformation matrix of 2D rotation: + [[cos(angle) -sin(angle) 0] + [sin(angle) cos(angle) 0] + [0 0 1]] + """ + R = [ + [math.cos(angle), -math.sin(angle), 0], + [math.sin(angle), math.cos(angle), 0], + [0, 0, 1], + ] + return np.array(R) + +def _transform(coords, matrix): + src = np.vstack((coords[:,0], coords[:,1], np.ones((coords.shape[0],)))) + dst = np.dot(src.transpose(), matrix.transpose()) + # rescale to homogeneous coordinates + dst[:,0] *= 1 / dst[:,2] + dst[:,1] *= 1 / dst[:,2] + dst[np.abs(dst) < EPS] = 0 + return dst[:,:2] + +def _transform_polynomial(coords, matrix): + x = coords[:,0] + y = coords[:,1] + u = len(matrix) + # 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] += matrix[pidx] * x ** (j - i) * y ** i + dst[:,1] += matrix[pidx+u/2] * x ** (j - i) * y ** i + pidx += 1 + + return dst + + +TRANSFORMATIONS = { + 'similarity': (_make_similarity, _transform), + 'affine': (_make_affine, _transform), + 'projective': (_make_projective, _transform), + 'polynomial': (_make_polynomial, _transform_polynomial), + 'rotation': (_make_rotation, _transform), +} + + +class Transformation(object): + + def __init__(self, ttype, matrix): + """Create transformation which contains the transformation parameters + and can perform forward and inverse transformations. + + Parameters + ---------- + ttype : str + one of similarity, affine, projective, polynomial, rotation + matrix : 3x3 array + homogeneous transformation matrix + """ + self.ttype = ttype + self.matrix = matrix + + def fwd(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : Nx2 array + source coordinates + + Returns + ------- + coords : Nx2 array + transformed coordinates + """ + return TRANSFORMATIONS[self.ttype][1](coords, self.matrix) + + def inv(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : Nx2 array + source coordinates + + Returns + ------- + coords : Nx2 array + transformed coordinates + """ + if self.ttype == 'polynomial': + raise Exception( + 'There is no explicit way to do the inverse polynomial ' + 'transformation. Instead determine the inverse transformation ' + 'parameters and use the forward transformation instead.') + matrix = np.linalg.inv(self.matrix) + return TRANSFORMATIONS[self.ttype][1](coords, matrix) + + +def make_tform(ttype, **kwargs): + """Create geometric transformation. + + 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, rotation + 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) + 'rotation' `angle` + + Alternatively you can explicitly define a 3x3 homogeneous transformation + matrix with the `matrix` parameter. + + See examples section below for usage. + + Returns + ------- + tform : :class:`Transformation` + tform object containing the transformation parameters + """ + + ttype = ttype.lower() + if ttype not in TRANSFORMATIONS: + raise NotImplemented('the transformation type \'%s\' is not' + 'implemented' % ttype) + if 'matrix' in kwargs: + matrix = kwargs['matrix'] + else: + matrix = TRANSFORMATIONS[ttype][0](**kwargs) + return Transformation(ttype, matrix) + +def warp(image, reverse_map=None, map_args={}, tform=None, + 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 : 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*. Also see examples below. + map_args : dict, optional + Keyword arguments passed to `reverse_map`. + tform : :class:`Transformation` object + The inverse transformation will be used to transform coordinates in the + *output image* into their corresponding coordinates in the + *source image*. + 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(reverse_map): + tf_coords = reverse_map(tf_coords, **map_args) + else: + tf_coords = tform.inv(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, ...]) + + # 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 = make_tform('projective', matrix=H) + return warp(image, tform=tform, output_shape=output_shape, order=order, + mode=mode, cval=cval) diff --git a/skimage/transform/project.py b/skimage/transform/project.py deleted file mode 100644 index a140e03e..00000000 --- a/skimage/transform/project.py +++ /dev/null @@ -1,137 +0,0 @@ -"""Image projection. - -""" - -import numpy as np -from scipy.ndimage import interpolation as ndii -from ._warp import _stackcopy - -__all__ = ['homography'] - -eps = np.finfo(float).eps - - -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) - - """ - if image.ndim < 2: - raise ValueError("Input must have more than 1 dimension.") - - image = np.atleast_3d(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) - - # TODO: Refactor this method to use transform.warp instead. - - # Construct transformed coordinates - rows, cols = output_shape[:2] - rows, cols = np.mgrid[:rows, :cols] - tf_coords = np.empty(shape=cols.shape, - dtype=[('cols', float), - ('rows', float), - ('z', float)]) - tf_coords['cols'], tf_coords['rows'] = cols, rows - tf_coords['z'] = 1 - tf_coords = tf_coords.view((float, 3)) - - tf_coords = np.dot(tf_coords, np.linalg.inv(H).transpose()) - tf_coords[np.absolute(tf_coords) < eps] = 0. - - # normalize coordinates - tf_coords[..., :2] /= tf_coords[..., 2, np.newaxis] - - # y-coordinate mapping - _stackcopy(coords[0, ...], tf_coords[..., 1]) - - # x-coordinate mapping - _stackcopy(coords[1, ...], tf_coords[..., 0]) - - # colour-coordinate mapping - coords[2, ...] = range(bands) - - # Prefilter not necessary for order 1 interpolation - prefilter = order > 1 - mapped = ndii.map_coordinates(image, coords, prefilter=prefilter, - mode=mode, order=order, cval=cval) - - return mapped.squeeze() diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py new file mode 100644 index 00000000..a227f597 --- /dev/null +++ b/skimage/transform/tests/test_geometric.py @@ -0,0 +1,146 @@ +import numpy as np +from numpy.testing import assert_array_almost_equal + +from skimage.transform.geometric import _stackcopy +from skimage.transform import make_tform +from skimage.transform import homography, fast_homography +from skimage import transform as tf, data, img_as_float +from skimage.color import rgb2gray + + +SRC = np.array([ + [-12.3705, -10.5075], + [-10.7865, 15.4305], + [8.6985, 10.8675], + [11.4975, -9.5715], + [7.8435, 7.4835], + [-5.3325, 6.5025], + [6.7905, -6.3765], + [-6.1695, -0.8235], +]) +DST = np.array([ + [0, 0], + [0, 5800], + [4900, 5800], + [4900, 0], + [4479, 4580], + [1176, 3660], + [3754, 790], + [1024, 1931], +]) + + +def test_stackcopy(): + layers = 4 + x = np.empty((3, 3, layers)) + y = np.eye(3, 3) + _stackcopy(x, y) + for i in range(layers): + assert_array_almost_equal(x[...,i], y) + +def test_similarity(): + #: exact solution + tform = make_tform('similarity', src=SRC[:2,:], dst=DST[:2,:]) + assert_array_almost_equal(tform.fwd(SRC[:2,:]), DST[:2,:]) + assert_array_almost_equal(tform.inv(tform.fwd(SRC)), SRC) + + #: over-determined + tform = make_tform('similarity', src=SRC, dst=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.inv(tform.fwd(SRC)), SRC) + +def test_affine(): + #: exact solution + tform = make_tform('affine', src=SRC[:3,:], dst=DST[:3,:]) + assert_array_almost_equal(tform.fwd(SRC[:3,:]), DST[:3,:]) + assert_array_almost_equal(tform.inv(tform.fwd(SRC)), SRC) + + #: over-determined + tform = make_tform('affine', src=SRC, dst=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.inv(tform.fwd(SRC)), SRC) + +def test_projective(): + #: exact solution + tform = make_tform('projective', src=SRC[:4,:], dst=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.inv(tform.fwd(SRC)), SRC) + + #: over-determined + tform = make_tform('projective', src=SRC[:4,:], dst=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.inv(tform.fwd(SRC)), SRC) + +def test_polynomial(): + tform = make_tform('polynomial', src=SRC, dst=DST, order=10) + assert_array_almost_equal(tform.fwd(SRC), DST, 6) + +def test_homography(): + x = img_as_float(np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1) + 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) + 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__": + from numpy.testing import run_module_suite + run_module_suite() diff --git a/skimage/transform/tests/test_project.py b/skimage/transform/tests/test_project.py deleted file mode 100644 index 25ebda20..00000000 --- a/skimage/transform/tests/test_project.py +++ /dev/null @@ -1,63 +0,0 @@ -import numpy as np -from numpy.testing import assert_array_almost_equal - -from skimage.transform._warp import _stackcopy -from skimage.transform import homography, fast_homography -from skimage import data -from skimage.color import rgb2gray - - -def test_stackcopy(): - layers = 4 - x = np.empty((3, 3, layers)) - y = np.eye(3, 3) - _stackcopy(x, y) - for i in range(layers): - assert_array_almost_equal(x[..., i], y) - - -def test_homography(): - x = np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1 - 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) - 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 - - -if __name__ == "__main__": - from numpy.testing import run_module_suite - run_module_suite() diff --git a/skimage/transform/tests/test_swirl.py b/skimage/transform/tests/test_swirl.py deleted file mode 100644 index d71f8231..00000000 --- a/skimage/transform/tests/test_swirl.py +++ /dev/null @@ -1,17 +0,0 @@ -import numpy as np -from numpy.testing import assert_array_almost_equal - -from skimage import transform as tf, data, img_as_float - - -def test_roundtrip(): - 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__": - np.testing.run_module_suite() From da3239d0ff879ad4823e3fcb468a34cf215551cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 4 Jun 2012 15:27:11 +0200 Subject: [PATCH 02/41] remove old import --- skimage/transform/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index f5621444..ad179cbb 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -1,7 +1,6 @@ from .hough_transform import * from .radon_transform import * from .finite_radon_transform import * -from .project import * from ._project import homography as fast_homography from .integral import * from .geometric import warp, make_tform, swirl, homography From bbd9280c2fc16cbef49afae436a174639a0770b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 4 Jun 2012 15:37:31 +0200 Subject: [PATCH 03/41] add example to make_tform doc string --- skimage/transform/geometric.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 0306af59..7578b6fb 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -298,6 +298,16 @@ def make_tform(ttype, **kwargs): ------- tform : :class:`Transformation` tform object containing the transformation parameters + + Examples + -------- + >>> import numpy as np + >>> from skimage.transform import make_tform + >>> src = np.array([0, 0, 10, 10]).reshape((2, 2)) + >>> dst = np.array([12, 14, 1, -20]).reshape((2, 2)) + >>> tform = make_tform('similarity', src=src, dst=dst) + >>> print tform.matrix + >>> print tform.inv(tform.fwd(src)) # == src """ ttype = ttype.lower() From 33a7ddec0e986b28bec4cc3ce282f94221e7afa5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 4 Jun 2012 18:14:38 +0200 Subject: [PATCH 04/41] update contributors --- CONTRIBUTORS.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index baaf064b..5961407d 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -107,3 +107,4 @@ Polygon, circle and ellipse drawing functions Adaptive thresholding Implementation of Matlab's `regionprops` + Estimation of geometric transformation parameters From fdf1b6dac1ceade98bf7df9fafff0c9ab493961f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 25 Jun 2012 19:57:02 +0200 Subject: [PATCH 05/41] fixe index errors and improve setup of some matrices --- skimage/transform/geometric.py | 89 +++++++++++++++------------------- 1 file changed, 38 insertions(+), 51 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 7578b6fb..89e8cd75 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -35,8 +35,8 @@ def _make_similarity(src, dst): X = a0*x - b0*y + a1 Y = b0*x + a0*y + a2 where the homogeneous transformation matrix is: - [[a1 -b1 a0] - [b1 a1 b0] + [[a0 -b0 a1] + [b0 a0 b1] [0 0 1]] """ xs = src[:,0] @@ -45,32 +45,27 @@ def _make_similarity(src, dst): yd = dst[:,1] rows = src.shape[0] + #: params: a0, a1, b0, b1 A = np.zeros((rows*2, 4)) - b = np.zeros((rows*2,)) - 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[:rows] = xd - b[rows:] = yd + + b = np.hstack([xd, yd]) a0, a1, b0, b1 = np.linalg.lstsq(A, b)[0] - matrix = np.eye(3) - matrix[0,0] = a0 - matrix[0,1] = - b0 - matrix[0,2] = a1 - matrix[1,0] = b0 - matrix[1,1] = a0 - matrix[1,2] = b1 + matrix = np.array([[a0, -b0, a1], + [b0, a0, b1], + [ 0, 0, 1]]) return matrix def _make_affine(src, dst): """Determine parameters of the 2D affine transformation: - X = a0*x + a1*y + a3 - Y = b0*x + b1*y + b3 + X = a0*x + a1*y + a2 + Y = b0*x + b1*y + b2 where the homogeneous transformation matrix is: [[a0 a1 a2] [b0 b1 b2] @@ -82,9 +77,8 @@ def _make_affine(src, dst): yd = dst[:,1] rows = src.shape[0] + #: params: a0, a1, a2, b0, b1, b2 A = np.zeros((rows*2, 6)) - b = np.zeros((rows*2,)) - A[:rows,0] = xs A[:rows,1] = ys A[:rows,2] = 1 @@ -92,22 +86,22 @@ def _make_affine(src, dst): A[rows:,4] = ys A[rows:,5] = 1 - b[:rows] = xd - b[rows:] = yd + b = np.hstack([xd, yd]) - params = np.linalg.lstsq(A, b)[0] - matrix = np.eye(3) - matrix[:2,:] = params.reshape((2, 3)) + a0, a1, a2, b0, b1, b2 = np.linalg.lstsq(A, b)[0] + matrix = np.array([[a0, a1, a2], + [b0, b1, b2], + [0, 0, 1]]) return matrix def _make_projective(src, dst): """Determine transformation matrix of the 2D projective transformation: - X = (a0 + a1*x + a2*y) / (c0*x + c1*y + c3) - Y = (b0 + b1*x + b2*y) / (c0*x + c1*y + c3) + 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 c3]] + [c0 c1 1]] """ xs = src[:,0] ys = src[:,1] @@ -115,10 +109,8 @@ def _make_projective(src, dst): yd = dst[:,1] rows = src.shape[0] + #: params: a0, a1, a2, b0, b1, b2, c0, c1 A = np.zeros((rows*2, 8)) - b = np.zeros((rows*2,)) - - A[:rows,0] = xs A[:rows,1] = ys A[:rows,2] = 1 @@ -129,12 +121,14 @@ def _make_projective(src, dst): A[rows:,5] = 1 A[rows:,6] = - yd * xs A[rows:,7] = - yd * ys - b[:rows] = dst[:,0] - b[rows:] = dst[:,1] - matrix = np.eye(3).flatten() - matrix[:8] = np.linalg.lstsq(A, b)[0] - return matrix.reshape((3, 3)) + b = np.hstack([xd, yd]) + + a0, a1, a2, b0, b1, b2, c0, c1 = np.linalg.lstsq(A, b)[0] + matrix = np.array([[a0, a1, a2], + [b0, b1, b2], + [c0, c1, 1]]) + return matrix def _make_polynomial(src, dst, order): """Determine parameters of 2D polynomial transformation of order n: @@ -149,17 +143,16 @@ def _make_polynomial(src, dst, order): # number of unknown polynomial coefficients u = (order + 1) * (order + 2) - A = np.zeros((rows*2, u)) - b = np.zeros((rows*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[:rows] = xd - b[rows:] = yd + + b = np.hstack([xd, yd]) return np.linalg.lstsq(A, b)[0] @@ -171,8 +164,8 @@ def _make_rotation(angle): """ R = [ [math.cos(angle), -math.sin(angle), 0], - [math.sin(angle), math.cos(angle), 0], - [0, 0, 1], + [math.sin(angle), math.cos(angle), 0], + [0, 0, 1], ] return np.array(R) @@ -180,8 +173,9 @@ def _transform(coords, matrix): src = np.vstack((coords[:,0], coords[:,1], np.ones((coords.shape[0],)))) dst = np.dot(src.transpose(), matrix.transpose()) # rescale to homogeneous coordinates - dst[:,0] *= 1 / dst[:,2] - dst[:,1] *= 1 / dst[:,2] + dst[:,0] /= dst[:,2] + dst[:,1] /= dst[:,2] + # values close to zero because of limited numerical precision dst[np.abs(dst) < EPS] = 0 return dst[:,:2] @@ -334,10 +328,6 @@ def warp(image, reverse_map=None, map_args={}, tform=None, coordinates in the *source image*. Also see examples below. map_args : dict, optional Keyword arguments passed to `reverse_map`. - tform : :class:`Transformation` object - The inverse transformation will be used to transform coordinates in the - *output image* into their corresponding coordinates in the - *source image*. output_shape : tuple (rows, cols) Shape of the output image generated. order : int @@ -385,10 +375,7 @@ def warp(image, reverse_map=None, map_args={}, tform=None, # Map each (x, y) pair to the source image according to # the user-provided mapping - if callable(reverse_map): - tf_coords = reverse_map(tf_coords, **map_args) - else: - tf_coords = tform.inv(tf_coords) + 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) @@ -564,5 +551,5 @@ def homography(image, H, output_shape=None, order=1, category=DeprecationWarning) tform = make_tform('projective', matrix=H) - return warp(image, tform=tform, output_shape=output_shape, order=order, - mode=mode, cval=cval) + return warp(image, reverse_map=tform.inv, output_shape=output_shape, + order=order, mode=mode, cval=cval) From acb1d71cd548728a0ae3e475361817001e7164fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 25 Jun 2012 20:03:59 +0200 Subject: [PATCH 06/41] apply PEP8 guideline --- skimage/transform/geometric.py | 136 ++++++++++++++++++--------------- 1 file changed, 74 insertions(+), 62 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 89e8cd75..7b59242e 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -30,6 +30,7 @@ def _stackcopy(a, b): else: a[:] = b + def _make_similarity(src, dst): """Determine parameters of the 2D similarity transformation: X = a0*x - b0*y + a1 @@ -39,20 +40,20 @@ def _make_similarity(src, dst): [b0 a0 b1] [0 0 1]] """ - xs = src[:,0] - ys = src[:,1] - xd = dst[:,0] - yd = dst[:,1] + 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 + 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]) @@ -62,6 +63,7 @@ def _make_similarity(src, dst): [ 0, 0, 1]]) return matrix + def _make_affine(src, dst): """Determine parameters of the 2D affine transformation: X = a0*x + a1*y + a2 @@ -71,20 +73,20 @@ def _make_affine(src, dst): [b0 b1 b2] [0 0 1]] """ - xs = src[:,0] - ys = src[:,1] - xd = dst[:,0] - yd = dst[:,1] + 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 + 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]) @@ -94,6 +96,7 @@ def _make_affine(src, dst): [0, 0, 1]]) return matrix + def _make_projective(src, dst): """Determine transformation matrix of the 2D projective transformation: X = (a0 + a1*x + a2*y) / (c0*x + c1*y + 1) @@ -103,24 +106,24 @@ def _make_projective(src, dst): [b0 b1 b2] [c0 c1 1]] """ - xs = src[:,0] - ys = src[:,1] - xd = dst[:,0] - yd = dst[:,1] + 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 + 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]) @@ -130,32 +133,34 @@ def _make_projective(src, dst): [c0, c1, 1]]) return matrix + def _make_polynomial(src, dst, order): """Determine parameters of 2D polynomial transformation of order n: 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 )) """ - xs = src[:,0] - ys = src[:,1] - xd = dst[:,0] - yd = dst[:,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)) + 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 + 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]) return np.linalg.lstsq(A, b)[0] + def _make_rotation(angle): """Determine homogeneous transformation matrix of 2D rotation: [[cos(angle) -sin(angle) 0] @@ -169,29 +174,32 @@ def _make_rotation(angle): ] return np.array(R) + def _transform(coords, matrix): - src = np.vstack((coords[:,0], coords[:,1], np.ones((coords.shape[0],)))) + 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] + dst[:, 0] /= dst[:, 2] + dst[:, 1] /= dst[:, 2] # values close to zero because of limited numerical precision dst[np.abs(dst) < EPS] = 0 - return dst[:,:2] + return dst[:, :2] + def _transform_polynomial(coords, matrix): - x = coords[:,0] - y = coords[:,1] + x = coords[:, 0] + y = coords[:, 1] u = len(matrix) # number of coefficients -> u = (order + 1) * (order + 2) - order = int((-3 + math.sqrt(9 - 4 * (2 - u))) / 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] += matrix[pidx] * x ** (j - i) * y ** i - dst[:,1] += matrix[pidx+u/2] * x ** (j - i) * y ** i + for j in xrange(order + 1): + for i in xrange(j + 1): + dst[:, 0] += matrix[pidx] * x ** (j - i) * y ** i + dst[:, 1] += matrix[pidx + u / 2] * x ** (j - i) * y ** i pidx += 1 return dst @@ -283,8 +291,8 @@ def make_tform(ttype, **kwargs): 'polynomial' `src, `dst`, `order` (polynomial order) 'rotation' `angle` - Alternatively you can explicitly define a 3x3 homogeneous transformation - matrix with the `matrix` parameter. + Alternatively you can explicitly define a 3x3 homogeneous + transformation matrix with the `matrix` parameter. See examples section below for usage. @@ -314,8 +322,9 @@ def make_tform(ttype, **kwargs): matrix = TRANSFORMATIONS[ttype][0](**kwargs) return Transformation(ttype, matrix) -def warp(image, reverse_map=None, map_args={}, tform=None, - output_shape=None, order=1, mode='constant', cval=0.): + +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 @@ -398,10 +407,11 @@ def warp(image, reverse_map=None, map_args={}, tform=None, # 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) + rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) # Ensure that the transformation decays to approximately 1/1000-th # within the specified radius. @@ -416,6 +426,7 @@ def _swirl_mapping(xy, center, rotation, strength, radius): 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. @@ -467,6 +478,7 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0, 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. From 8bde92b66cb2d296fd86b9af61a794fc8399b91b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 2 Jul 2012 09:27:35 +0200 Subject: [PATCH 07/41] remove inconsistent numeric correction and fix test case --- skimage/transform/geometric.py | 5 ----- skimage/transform/tests/test_geometric.py | 6 ++++-- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 7b59242e..71ec7752 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -5,9 +5,6 @@ from scipy import ndimage from skimage.util import img_as_float -EPS = np.spacing(1) - - def _stackcopy(a, b): """Copy b into each color layer of a, such that:: @@ -182,8 +179,6 @@ def _transform(coords, matrix): # rescale to homogeneous coordinates dst[:, 0] /= dst[:, 2] dst[:, 1] /= dst[:, 2] - # values close to zero because of limited numerical precision - dst[np.abs(dst) < EPS] = 0 return dst[:, :2] diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index a227f597..d0edace1 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -92,10 +92,12 @@ def test_polynomial(): assert_array_almost_equal(tform.fwd(SRC), DST, 6) def test_homography(): - x = img_as_float(np.arange(9, dtype=np.uint8).reshape((3, 3)) + 1) + 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),2], + [np.sin(theta), np.cos(theta),4], [0, 0, 1]]) x90 = homography(x, M, order=1) assert_array_almost_equal(x90, np.rot90(x)) From cb3c93a1107b0b7f0308072693e18d48cfef9c5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 9 Jul 2012 23:00:52 +0200 Subject: [PATCH 08/41] change API design and rename some functions --- skimage/transform/__init__.py | 3 +- skimage/transform/geometric.py | 145 ++++++++++++---------- skimage/transform/tests/test_geometric.py | 51 ++++---- 3 files changed, 108 insertions(+), 91 deletions(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index ad179cbb..296503d2 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -3,4 +3,5 @@ from .radon_transform import * from .finite_radon_transform import * from ._project import homography as fast_homography from .integral import * -from .geometric import warp, make_tform, swirl, homography +from .geometric import warp, estimate_transformation, geometric_transform, \ + swirl, homography diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 71ec7752..abe27d4c 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -28,7 +28,7 @@ def _stackcopy(a, b): a[:] = b -def _make_similarity(src, dst): +def _estimate_similarity(src, dst): """Determine parameters of the 2D similarity transformation: X = a0*x - b0*y + a1 Y = b0*x + a0*y + a2 @@ -36,6 +36,7 @@ def _make_similarity(src, dst): [[a0 -b0 a1] [b0 a0 b1] [0 0 1]] + """ xs = src[:, 0] ys = src[:, 1] @@ -61,7 +62,7 @@ def _make_similarity(src, dst): return matrix -def _make_affine(src, dst): +def _estimate_affine(src, dst): """Determine parameters of the 2D affine transformation: X = a0*x + a1*y + a2 Y = b0*x + b1*y + b2 @@ -69,6 +70,7 @@ def _make_affine(src, dst): [[a0 a1 a2] [b0 b1 b2] [0 0 1]] + """ xs = src[:, 0] ys = src[:, 1] @@ -94,7 +96,7 @@ def _make_affine(src, dst): return matrix -def _make_projective(src, dst): +def _estimate_projective(src, dst): """Determine transformation matrix of the 2D projective transformation: X = (a0 + a1*x + a2*y) / (c0*x + c1*y + 1) Y = (b0 + b1*x + b2*y) / (c0*x + c1*y + 1) @@ -102,6 +104,7 @@ def _make_projective(src, dst): [[a0 a1 a2] [b0 b1 b2] [c0 c1 1]] + """ xs = src[:, 0] ys = src[:, 1] @@ -131,10 +134,11 @@ def _make_projective(src, dst): return matrix -def _make_polynomial(src, dst, order): +def _estimate_polynomial(src, dst, order): """Determine parameters of 2D polynomial transformation of order n: 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 )) + """ xs = src[:, 0] ys = src[:, 1] @@ -158,21 +162,21 @@ def _make_polynomial(src, dst, order): return np.linalg.lstsq(A, b)[0] -def _make_rotation(angle): - """Determine homogeneous transformation matrix of 2D rotation: - [[cos(angle) -sin(angle) 0] - [sin(angle) cos(angle) 0] - [0 0 1]] +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 """ - R = [ - [math.cos(angle), -math.sin(angle), 0], - [math.sin(angle), math.cos(angle), 0], - [0, 0, 1], - ] - return np.array(R) - - -def _transform(coords, matrix): x, y = np.transpose(coords) src = np.vstack((x, y, np.ones_like(x))) dst = np.dot(src.transpose(), matrix.transpose()) @@ -200,32 +204,28 @@ def _transform_polynomial(coords, matrix): return dst -TRANSFORMATIONS = { - 'similarity': (_make_similarity, _transform), - 'affine': (_make_affine, _transform), - 'projective': (_make_projective, _transform), - 'polynomial': (_make_polynomial, _transform_polynomial), - 'rotation': (_make_rotation, _transform), -} +class GeometricTransformation(object): - -class Transformation(object): - - def __init__(self, ttype, matrix): - """Create transformation which contains the transformation parameters - and can perform forward and inverse transformations. + def __init__(self, ttype, params, transform_func): + """Create geometric transformation which contains the transformation + parameters and can perform forward and reverse transformations. Parameters ---------- ttype : str - one of similarity, affine, projective, polynomial, rotation - matrix : 3x3 array - homogeneous transformation matrix + transformation type - one of 'similarity', 'affine', 'projective', + 'polynomial' + params : array + transformation parameters + transform_func : callable = func(coords, matrix) + transformation function + """ self.ttype = ttype - self.matrix = matrix + self.params = params + self.transform_func = transform_func - def fwd(self, coords): + def forward(self, coords): """Apply forward transformation. Parameters @@ -237,11 +237,12 @@ class Transformation(object): ------- coords : Nx2 array transformed coordinates - """ - return TRANSFORMATIONS[self.ttype][1](coords, self.matrix) - def inv(self, coords): - """Apply inverse transformation. + """ + return self.transform_func(coords, self.params) + + def reverse(self, coords): + """Apply reverse transformation. Parameters ---------- @@ -252,30 +253,38 @@ class Transformation(object): ------- coords : Nx2 array transformed coordinates + """ if self.ttype == 'polynomial': raise Exception( - 'There is no explicit way to do the inverse polynomial ' - 'transformation. Instead determine the inverse transformation ' - 'parameters and use the forward transformation instead.') - matrix = np.linalg.inv(self.matrix) - return TRANSFORMATIONS[self.ttype][1](coords, matrix) + '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.') + inv_matrix = np.linalg.inv(self.params) + return self.transform_func(coords, inv_matrix) -def make_tform(ttype, **kwargs): - """Create geometric transformation. +ESTIMATED_TRANSFORMATIONS = { + 'similarity': (_estimate_similarity, geometric_transform), + 'affine': (_estimate_affine, geometric_transform), + 'projective': (_estimate_projective, geometric_transform), + 'polynomial': (_estimate_polynomial, _transform_polynomial), +} + + +def estimate_transformation(ttype, *args, **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 : str - one of similarity, affine, projective, polynomial, rotation + one of similarity, affine, projective, polynomial kwargs : array or int function parameters (src, dst, n, angle): @@ -284,17 +293,14 @@ def make_tform(ttype, **kwargs): 'affine' `src, `dst` 'projective' `src, `dst` 'polynomial' `src, `dst`, `order` (polynomial order) - 'rotation' `angle` - - Alternatively you can explicitly define a 3x3 homogeneous - transformation matrix with the `matrix` parameter. See examples section below for usage. Returns ------- - tform : :class:`Transformation` - tform object containing the transformation parameters + tform : :class:`GeometricTransformation` + tform object containing the transformation parameters and providing + access to forward and reverse transformation functions Examples -------- @@ -302,20 +308,23 @@ def make_tform(ttype, **kwargs): >>> from skimage.transform import make_tform >>> src = np.array([0, 0, 10, 10]).reshape((2, 2)) >>> dst = np.array([12, 14, 1, -20]).reshape((2, 2)) - >>> tform = make_tform('similarity', src=src, dst=dst) - >>> print tform.matrix - >>> print tform.inv(tform.fwd(src)) # == src - """ + >>> tform = estimate_transformation('similarity', src, dst) + >>> print tform.params + >>> print tform.reverse(tform.forward(src)) # == src + >>> # warp image using the transformation + >>> from skimage import data + >>> image = data.camera() + >>> warp(image, reverse_map=tform.forward) + >>> warp(image, reverse_map=tform.reverse) + """ ttype = ttype.lower() - if ttype not in TRANSFORMATIONS: + if ttype not in ESTIMATED_TRANSFORMATIONS: raise NotImplemented('the transformation type \'%s\' is not' 'implemented' % ttype) - if 'matrix' in kwargs: - matrix = kwargs['matrix'] - else: - matrix = TRANSFORMATIONS[ttype][0](**kwargs) - return Transformation(ttype, matrix) + matrix = ESTIMATED_TRANSFORMATIONS[ttype][0](*args, **kwargs) + transform_func = ESTIMATED_TRANSFORMATIONS[ttype][1] + return GeometricTransformation(ttype, matrix, transform_func) def warp(image, reverse_map=None, map_args={}, output_shape=None, order=1, @@ -557,6 +566,6 @@ def homography(image, H, output_shape=None, order=1, 'use the `warp` and `tform` function instead', category=DeprecationWarning) - tform = make_tform('projective', matrix=H) - return warp(image, reverse_map=tform.inv, output_shape=output_shape, + tform = GeometricTransformation('projective', H, geometric_transform) + 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 d0edace1..247637bb 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -2,7 +2,7 @@ import numpy as np from numpy.testing import assert_array_almost_equal from skimage.transform.geometric import _stackcopy -from skimage.transform import make_tform +from skimage.transform import estimate_transformation from skimage.transform import homography, fast_homography from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray @@ -36,60 +36,65 @@ def test_stackcopy(): y = np.eye(3, 3) _stackcopy(x, y) for i in range(layers): - assert_array_almost_equal(x[...,i], y) + assert_array_almost_equal(x[..., i], y) + def test_similarity(): #: exact solution - tform = make_tform('similarity', src=SRC[:2,:], dst=DST[:2,:]) - assert_array_almost_equal(tform.fwd(SRC[:2,:]), DST[:2,:]) - assert_array_almost_equal(tform.inv(tform.fwd(SRC)), SRC) + 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) #: over-determined - tform = make_tform('similarity', src=SRC, dst=DST) + 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.inv(tform.fwd(SRC)), SRC) + assert_array_almost_equal(tform.params, ref) + assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + def test_affine(): #: exact solution - tform = make_tform('affine', src=SRC[:3,:], dst=DST[:3,:]) - assert_array_almost_equal(tform.fwd(SRC[:3,:]), DST[:3,:]) - assert_array_almost_equal(tform.inv(tform.fwd(SRC)), SRC) + 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) #: over-determined - tform = make_tform('affine', src=SRC, dst=DST) + 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.inv(tform.fwd(SRC)), SRC) + assert_array_almost_equal(tform.params, ref) + assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + def test_projective(): #: exact solution - tform = make_tform('projective', src=SRC[:4,:], dst=DST[:4,:]) + 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.inv(tform.fwd(SRC)), SRC) + assert_array_almost_equal(tform.params, ref, 6) + assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) #: over-determined - tform = make_tform('projective', src=SRC[:4,:], dst=DST[:4,:]) + 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.inv(tform.fwd(SRC)), SRC) + assert_array_almost_equal(tform.params, ref, 6) + assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) + def test_polynomial(): - tform = make_tform('polynomial', src=SRC, dst=DST, order=10) - assert_array_almost_equal(tform.fwd(SRC), DST, 6) + tform = estimate_transformation('polynomial', SRC, DST, order=10) + assert_array_almost_equal(tform.forward(SRC), DST, 6) + def test_homography(): x = np.zeros((5,5), dtype=np.uint8) @@ -102,6 +107,7 @@ def test_homography(): 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] @@ -133,6 +139,7 @@ def test_fast_homography(): d = np.mean(np.abs(p0 - p1)) assert d < 0.2 + def test_swirl(): image = img_as_float(data.checkerboard()) From 2a75b78838112d214a3ac44ddde49239a099d05b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 19:36:32 +0200 Subject: [PATCH 09/41] change arguments of function estimate_transformation --- skimage/transform/geometric.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index abe27d4c..68f2dae8 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -273,7 +273,7 @@ ESTIMATED_TRANSFORMATIONS = { } -def estimate_transformation(ttype, *args, **kwargs): +def estimate_transformation(ttype, src, dst, order=None): """Estimate 2D geometric transformation parameters. You can determine the over-, well- and under-determined parameters @@ -322,7 +322,10 @@ def estimate_transformation(ttype, *args, **kwargs): if ttype not in ESTIMATED_TRANSFORMATIONS: raise NotImplemented('the transformation type \'%s\' is not' 'implemented' % ttype) - matrix = ESTIMATED_TRANSFORMATIONS[ttype][0](*args, **kwargs) + args = [src, dst] + if order is not None and ttype == 'polynomial': + args.append(order) + matrix = ESTIMATED_TRANSFORMATIONS[ttype][0](*args) transform_func = ESTIMATED_TRANSFORMATIONS[ttype][1] return GeometricTransformation(ttype, matrix, transform_func) From ec5c339c682fb2187d77d2c35ce12cf77a95bb6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 19:37:33 +0200 Subject: [PATCH 10/41] fix wrong exception type --- skimage/transform/geometric.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 68f2dae8..4519af88 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -320,8 +320,8 @@ def estimate_transformation(ttype, src, dst, order=None): """ ttype = ttype.lower() if ttype not in ESTIMATED_TRANSFORMATIONS: - raise NotImplemented('the transformation type \'%s\' is not' - 'implemented' % ttype) + raise ValueError('the transformation type \'%s\' is not' + 'implemented' % ttype) args = [src, dst] if order is not None and ttype == 'polynomial': args.append(order) From 640edc2a6261e54cf4a0108b3120ceffdf339d11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 19:38:19 +0200 Subject: [PATCH 11/41] fix doc string formatting of function estimate_transformation --- skimage/transform/geometric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 4519af88..ec631043 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -285,7 +285,7 @@ def estimate_transformation(ttype, src, dst, order=None): ---------- ttype : str one of similarity, affine, projective, polynomial - kwargs : array or int + kwargs :: array or int function parameters (src, dst, n, angle): NAME / TTYPE FUNCTION PARAMETERS From e2ce1d63de70243ad3f5da0959cbde90db40e780 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 22:58:10 +0200 Subject: [PATCH 12/41] redesign class interface --- skimage/transform/__init__.py | 3 +- skimage/transform/geometric.py | 523 ++++++++++++++-------- skimage/transform/tests/test_geometric.py | 40 +- 3 files changed, 374 insertions(+), 192 deletions(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 296503d2..b0790992 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -4,4 +4,5 @@ from .finite_radon_transform import * from ._project import homography as fast_homography from .integral import * from .geometric import warp, estimate_transformation, geometric_transform, \ - swirl, homography + SimilarityTransformation, AffineTransformation, ProjectiveTransformation, \ + PolynomialTransformation, swirl, homography diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index ec631043..d5f50b17 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -28,140 +28,6 @@ def _stackcopy(a, b): a[:] = b -def _estimate_similarity(src, dst): - """Determine parameters of the 2D similarity transformation: - X = a0*x - b0*y + a1 - Y = b0*x + a0*y + a2 - where the homogeneous transformation matrix is: - [[a0 -b0 a1] - [b0 a0 b1] - [0 0 1]] - - """ - 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] - matrix = np.array([[a0, -b0, a1], - [b0, a0, b1], - [ 0, 0, 1]]) - return matrix - - -def _estimate_affine(src, dst): - """Determine parameters of the 2D affine transformation: - X = a0*x + a1*y + a2 - Y = b0*x + b1*y + b2 - where the homogeneous transformation matrix is: - [[a0 a1 a2] - [b0 b1 b2] - [0 0 1]] - - """ - 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] - matrix = np.array([[a0, a1, a2], - [b0, b1, b2], - [0, 0, 1]]) - return matrix - - -def _estimate_projective(src, dst): - """Determine transformation matrix of the 2D projective transformation: - 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]] - - """ - 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] - matrix = np.array([[a0, a1, a2], - [b0, b1, b2], - [c0, c1, 1]]) - return matrix - - -def _estimate_polynomial(src, dst, order): - """Determine parameters of 2D polynomial transformation of order n: - 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 )) - - """ - 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]) - - return np.linalg.lstsq(A, b)[0] - - def geometric_transform(coords, matrix): """Apply 2D geometric transformation. @@ -186,44 +52,20 @@ def geometric_transform(coords, matrix): return dst[:, :2] -def _transform_polynomial(coords, matrix): - x = coords[:, 0] - y = coords[:, 1] - u = len(matrix) - # 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] += matrix[pidx] * x ** (j - i) * y ** i - dst[:, 1] += matrix[pidx + u / 2] * x ** (j - i) * y ** i - pidx += 1 - - return dst - - class GeometricTransformation(object): - def __init__(self, ttype, params, transform_func): + def __init__(self, matrix=None): """Create geometric transformation which contains the transformation parameters and can perform forward and reverse transformations. Parameters ---------- - ttype : str - transformation type - one of 'similarity', 'affine', 'projective', - 'polynomial' - params : array - transformation parameters - transform_func : callable = func(coords, matrix) - transformation function + matrix : 3x3 array, optional + homogeneous transformation matrix """ - self.ttype = ttype - self.params = params - self.transform_func = transform_func + self.matrix = matrix + self.inverse_matrix = None def forward(self, coords): """Apply forward transformation. @@ -239,7 +81,9 @@ class GeometricTransformation(object): transformed coordinates """ - return self.transform_func(coords, self.params) + 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. @@ -255,21 +99,332 @@ class GeometricTransformation(object): transformed coordinates """ - if self.ttype == 'polynomial': - 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.') - inv_matrix = np.linalg.inv(self.params) - return self.transform_func(coords, inv_matrix) + 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 union(self, other): + return GeometricTransformation(self.matrix.dot(other.matrix)) + + def __mul__(self, other): + return self.union(self, other) + + def __add__(self, other): + return self.union(self, other) -ESTIMATED_TRANSFORMATIONS = { - 'similarity': (_estimate_similarity, geometric_transform), - 'affine': (_estimate_affine, geometric_transform), - 'projective': (_estimate_projective, geometric_transform), - 'polynomial': (_estimate_polynomial, _transform_polynomial), +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): + + def estimate(self, src, dst): + """Estimate transformation matrix of the 2D projective transformation: + 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]] + + 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): + + 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): + """Estimate parameters of 2D polynomial transformation of order n: + 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 )) + + 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): + 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.') + + def union(self, other): + raise Exception('Cannot unite polynomial transformations.') + + def __mul__(self, other): + return self.union(self, other) + + def __add__(self, other): + return self.union(self, other) + + +TRANSFORMATIONS = { + 'similarity': SimilarityTransformation, + 'affine': AffineTransformation, + 'projective': ProjectiveTransformation, + 'polynomial': PolynomialTransformation, } @@ -298,7 +453,7 @@ def estimate_transformation(ttype, src, dst, order=None): Returns ------- - tform : :class:`GeometricTransformation` + tform : subclass of :class:`GeometricTransformation` tform object containing the transformation parameters and providing access to forward and reverse transformation functions @@ -309,7 +464,7 @@ def estimate_transformation(ttype, src, dst, order=None): >>> src = np.array([0, 0, 10, 10]).reshape((2, 2)) >>> dst = np.array([12, 14, 1, -20]).reshape((2, 2)) >>> tform = estimate_transformation('similarity', src, dst) - >>> print tform.params + >>> print tform.matrix >>> print tform.reverse(tform.forward(src)) # == src >>> # warp image using the transformation >>> from skimage import data @@ -319,15 +474,15 @@ def estimate_transformation(ttype, src, dst, order=None): """ ttype = ttype.lower() - if ttype not in ESTIMATED_TRANSFORMATIONS: + 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) - matrix = ESTIMATED_TRANSFORMATIONS[ttype][0](*args) - transform_func = ESTIMATED_TRANSFORMATIONS[ttype][1] - return GeometricTransformation(ttype, matrix, transform_func) + tform = TRANSFORMATIONS[ttype]() + tform.estimate(*args) + return tform def warp(image, reverse_map=None, map_args={}, output_shape=None, order=1, @@ -569,6 +724,6 @@ def homography(image, H, output_shape=None, order=1, 'use the `warp` and `tform` function instead', category=DeprecationWarning) - tform = GeometricTransformation('projective', H, geometric_transform) + 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 247637bb..0285a53d 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -2,7 +2,9 @@ import numpy as np from numpy.testing import assert_array_almost_equal from skimage.transform.geometric import _stackcopy -from skimage.transform import estimate_transformation +from skimage.transform import estimate_transformation, \ + SimilarityTransformation, AffineTransformation, ProjectiveTransformation, \ + PolynomialTransformation from skimage.transform import homography, fast_homography from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray @@ -39,7 +41,7 @@ def test_stackcopy(): assert_array_almost_equal(x[..., i], y) -def test_similarity(): +def test_similarity_estimation(): #: exact solution tform = estimate_transformation('similarity', SRC[:2, :], DST[:2, :]) assert_array_almost_equal(tform.forward(SRC[:2, :]), DST[:2, :]) @@ -51,11 +53,22 @@ def test_similarity(): [[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.params, ref) + assert_array_almost_equal(tform.matrix, ref) assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) -def test_affine(): +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, :]) @@ -67,10 +80,23 @@ def test_affine(): [[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.params, ref) + 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) + + def test_projective(): #: exact solution tform = estimate_transformation('projective', SRC[:4, :], DST[:4, :]) @@ -78,7 +104,7 @@ def test_projective(): [[ 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.params, ref, 6) + assert_array_almost_equal(tform.matrix, ref, 6) assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) #: over-determined @@ -87,7 +113,7 @@ def test_projective(): [[ 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.params, ref, 6) + assert_array_almost_equal(tform.matrix, ref, 6) assert_array_almost_equal(tform.reverse(tform.forward(SRC)), SRC) From 234810be1070f3f3012cb14a192dae16c45e5bee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 23:01:07 +0200 Subject: [PATCH 13/41] fix inconsistent doc strings --- skimage/transform/geometric.py | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index d5f50b17..03e088a2 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -292,14 +292,19 @@ class AffineTransformation(GeometricTransformation): 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): - """Estimate transformation matrix of the 2D projective transformation: - 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]] + """Set the transformation matrix with the explicit transformation + parameters. Parameters ---------- @@ -338,6 +343,12 @@ class ProjectiveTransformation(GeometricTransformation): 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. @@ -351,9 +362,8 @@ class PolynomialTransformation(GeometricTransformation): self.coeffs = coeffs def estimate(self, src, dst, order): - """Estimate parameters of 2D polynomial transformation of order n: - 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 )) + """Set the transformation matrix with the explicit transformation + parameters. Parameters ---------- From b2ca8335094fb5787cc3a56d00239376dc5db38b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 10 Jul 2012 23:20:32 +0200 Subject: [PATCH 14/41] fix transformation union and add test case --- skimage/transform/geometric.py | 10 +++++++--- skimage/transform/tests/test_geometric.py | 21 +++++++++++++++++++++ 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 03e088a2..06b0839e 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -106,13 +106,17 @@ class GeometricTransformation(object): return geometric_transform(coords, self.inverse_matrix) def union(self, other): - return GeometricTransformation(self.matrix.dot(other.matrix)) + if type(self) == type(other): + transformation = self.__class__ + else: + transformation = GeometricTransformation + return transformation(self.matrix.dot(other.matrix)) def __mul__(self, other): - return self.union(self, other) + return self.union(other) def __add__(self, other): - return self.union(self, other) + return self.union(other) class SimilarityTransformation(GeometricTransformation): diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 0285a53d..14e839f0 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -122,6 +122,27 @@ def test_polynomial(): assert_array_almost_equal(tform.forward(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) + + tform = tform1.union(tform2) + tform = tform1 + tform2 + tform = tform1 * tform2 + + assert_array_almost_equal(tform.scale, scale1 * scale2) + assert_array_almost_equal(tform.rotation, rotation1 + rotation2) + + def test_homography(): x = np.zeros((5,5), dtype=np.uint8) x[1, 1] = 255 From 4dcf8528bfc982f049d13d7074de75a6ef9c9e4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 14 Jul 2012 09:30:52 +0200 Subject: [PATCH 15/41] change interface of transformation merging --- skimage/transform/geometric.py | 16 ++-------------- skimage/transform/tests/test_geometric.py | 1 - 2 files changed, 2 insertions(+), 15 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 06b0839e..6dffeb8a 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -105,18 +105,15 @@ class GeometricTransformation(object): self.inverse_matrix = np.linalg.inv(self.matrix) return geometric_transform(coords, self.inverse_matrix) - def union(self, other): + def __mul__(self, other): if type(self) == type(other): transformation = self.__class__ else: transformation = GeometricTransformation return transformation(self.matrix.dot(other.matrix)) - def __mul__(self, other): - return self.union(other) - def __add__(self, other): - return self.union(other) + return self.__mul__(other) class SimilarityTransformation(GeometricTransformation): @@ -424,15 +421,6 @@ class PolynomialTransformation(GeometricTransformation): 'parameters by exchanging source and destination coordinates.' 'Then apply the forward transformation.') - def union(self, other): - raise Exception('Cannot unite polynomial transformations.') - - def __mul__(self, other): - return self.union(self, other) - - def __add__(self, other): - return self.union(self, other) - TRANSFORMATIONS = { 'similarity': SimilarityTransformation, diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 14e839f0..b797942e 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -135,7 +135,6 @@ def test_union(): translation2 = (0, 0) tform2.from_params(scale2, rotation2, translation2) - tform = tform1.union(tform2) tform = tform1 + tform2 tform = tform1 * tform2 From 9dbad0023caf6f523c5fd3471ab670653cd4ed62 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 14 Jul 2012 09:50:36 +0200 Subject: [PATCH 16/41] add support for using transformation objects in warp function --- skimage/transform/geometric.py | 10 +++++++--- skimage/transform/tests/test_geometric.py | 24 ++++++++++++++++++----- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 6dffeb8a..2e8e3ab8 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -495,10 +495,12 @@ def warp(image, reverse_map=None, map_args={}, output_shape=None, order=1, ---------- image : 2-D array Input image. - reverse_map : callable xy = f(xy, **kwargs) - Reverse coordinate map. A function that transforms a Px2 array of + 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*. Also see examples below. + 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) @@ -548,6 +550,8 @@ def warp(image, reverse_map=None, map_args={}, output_shape=None, order=1, # 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 diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index b797942e..dd911529 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -2,10 +2,9 @@ import numpy as np from numpy.testing import assert_array_almost_equal from skimage.transform.geometric import _stackcopy -from skimage.transform import estimate_transformation, \ - SimilarityTransformation, AffineTransformation, ProjectiveTransformation, \ - PolynomialTransformation -from skimage.transform import homography, fast_homography +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 @@ -142,8 +141,23 @@ def test_union(): 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 = np.zeros((5, 5), dtype=np.uint8) x[1, 1] = 255 x = img_as_float(x) theta = -np.pi/2 From 5feafee2203d743b68b6c3c0204629cfcad8314d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 14 Jul 2012 10:06:21 +0200 Subject: [PATCH 17/41] extend doc string example for geometric transformations --- skimage/transform/geometric.py | 31 ++++++++++++++--------- skimage/transform/tests/test_geometric.py | 1 - 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 2e8e3ab8..3e67cbeb 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -105,15 +105,12 @@ class GeometricTransformation(object): self.inverse_matrix = np.linalg.inv(self.matrix) return geometric_transform(coords, self.inverse_matrix) - def __mul__(self, other): + def __add__(self, other): if type(self) == type(other): transformation = self.__class__ else: transformation = GeometricTransformation - return transformation(self.matrix.dot(other.matrix)) - - def __add__(self, other): - return self.__mul__(other) + return transformation(other.matrix.dot(self.matrix)) class SimilarityTransformation(GeometricTransformation): @@ -462,17 +459,27 @@ def estimate_transformation(ttype, src, dst, order=None): Examples -------- >>> import numpy as np - >>> from skimage.transform import make_tform + >>> 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 = estimate_transformation('similarity', src, dst) - >>> print tform.matrix - >>> print tform.reverse(tform.forward(src)) # == src - >>> # warp image using the transformation + >>> 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() - >>> warp(image, reverse_map=tform.forward) - >>> warp(image, reverse_map=tform.reverse) + >>> 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() diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index dd911529..1648e6f1 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -135,7 +135,6 @@ def test_union(): tform2.from_params(scale2, rotation2, translation2) tform = tform1 + tform2 - tform = tform1 * tform2 assert_array_almost_equal(tform.scale, scale1 * scale2) assert_array_almost_equal(tform.rotation, rotation1 + rotation2) From d7b2c5b51b28195f75c4c767ddbb214af4522e9f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 14 Jul 2012 10:13:10 +0200 Subject: [PATCH 18/41] add missing doc string for polynomial forward transformation --- skimage/transform/geometric.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index 3e67cbeb..edb3d014 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -395,6 +395,19 @@ class PolynomialTransformation(GeometricTransformation): 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) From afb479d766ccb97291acae27afd484a63c28c7e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sun, 15 Jul 2012 17:51:34 +0200 Subject: [PATCH 19/41] geometric_transform can transform single coordinate tuple --- skimage/transform/geometric.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/skimage/transform/geometric.py b/skimage/transform/geometric.py index edb3d014..2051e241 100644 --- a/skimage/transform/geometric.py +++ b/skimage/transform/geometric.py @@ -43,13 +43,22 @@ def geometric_transform(coords, matrix): 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] - return dst[:, :2] + + if shape == (2,): + return dst[0, :2] + else: + return dst[:, :2] class GeometricTransformation(object): From 8e8e2b99a0fc0ce572b31451078628709a80e7d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sun, 15 Jul 2012 19:03:44 +0200 Subject: [PATCH 20/41] add short tutorial for geometric transformations --- doc/examples/applications/plot_geometric.py | 134 ++++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 doc/examples/applications/plot_geometric.py diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/applications/plot_geometric.py new file mode 100644 index 00000000..90fbce51 --- /dev/null +++ b/doc/examples/applications/plot_geometric.py @@ -0,0 +1,134 @@ +""" +=============================== +Using geometric transformations +=============================== + +In this example, we will see how to use geometric transformations in the context +of image processing. +""" + +import math +import numpy as np +import matplotlib.pyplot as plt + +from skimage import data +from skimage import transform as tf + +margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) + +""" +Basics +====== + +Several different geometric transformation types are supported: similarity, +affine, projective and polynomial. + +Geometric transformations can either be created using the explicit parameters +(e.g. scale, shear, rotation and translation) or the transformation matrix: +""" + +#: create using explicit parameters +tform = tf.SimilarityTransformation() +scale = 1 +rotation = math.pi/2 +translation = (0, 1) +tform.from_params(scale, rotation, translation) +print tform.matrix + +#: create using transformation matrix +matrix = tform.matrix.copy() +matrix[1, 2] = 2 +tform2 = tf.SimilarityTransformation(matrix) + +""" +These transformation objects can be used to forward and reverse transform +coordinates between the source and destination coordinate systems: +""" + +coord = [1, 0] +print tform2.forward(coord) +print tform2.reverse(tform.forward(coord)) + +""" +Image warping +============= + +Geometric transformations can also be used to warp images: +""" + +text = data.text() +tform.from_params(1, math.pi/4, (text.shape[0] / 2, -100)) + +# uses tform.reverse, alternatively use tf.warp(text, tform.reverse) +rotated = tf.warp(text, tform) +back_rotated = tf.warp(rotated, tform.forward) + +plt.figure(figsize=(8, 3)) +plt.subplot(131) +plt.imshow(text) +plt.axis('off') +plt.gray() +plt.subplot(132) +plt.imshow(rotated) +plt.axis('off') +plt.gray() +plt.subplot(133) +plt.imshow(back_rotated) +plt.axis('off') +plt.gray() +plt.subplots_adjust(**margins) + +""" +.. image:: PLOT2RST.current_figure + +Parameter estimation +==================== + +In addition to the basic functionality mentioned above you can also estimate the +parameters of a geometric transformation using the least-squares method. + +This can amongst other things be used for image registration or rectification, +where you have a set of control points or homologous points in two images. + +Let's assume we want to recognize letters on a photograph which was not taken +from the front but at a certain angle. In the simplest case of a plane paper +surface the letters are projectively distorted. Simple matching algorithms would +not be able to match such symbols. One solution to this problem would be to warp +the image so that the distortion is removed and then apply a matching algorithm: +""" + +text = data.text() + +src = np.array(( + (155, 15), + (65, 40), + (260, 130), + (360, 95) +)) +dst = np.array(( + (0, 0), + (0, 50), + (300, 50), + (300, 0) +)) + +tform3 = tf.estimate_transformation('projective', src, dst) +warped = tf.warp(text, tform3, output_shape=(50, 300)) + +plt.figure(figsize=(8, 3)) +plt.subplot(211) +plt.imshow(text) +plt.plot(src[:, 0], src[:, 1], '.r') +plt.axis('off') +plt.gray() +plt.subplot(212) +plt.imshow(warped) +plt.axis('off') +plt.gray() +plt.subplots_adjust(**margins) + +""" +.. image:: PLOT2RST.current_figure +""" + +plt.show() From a130b8d2d9dbe9558552a2e7e0f4fb03f49e8b5b Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Fri, 20 Jul 2012 17:57:15 -0400 Subject: [PATCH 21/41] Refactor geometric transforms. --- skimage/transform/__init__.py | 7 +- skimage/transform/_geometric.py | 563 ++++++++++++++++ skimage/transform/_warps.py | 160 +++++ skimage/transform/geometric.py | 764 ---------------------- skimage/transform/tests/test_geometric.py | 171 +---- skimage/transform/tests/test_warps.py | 79 +++ 6 files changed, 830 insertions(+), 914 deletions(-) create mode 100644 skimage/transform/_geometric.py create mode 100644 skimage/transform/_warps.py delete mode 100644 skimage/transform/geometric.py create mode 100644 skimage/transform/tests/test_warps.py 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() From f971c25a0b8a1e0a29276c6e899c6a5f9f00b44c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 08:21:31 +0200 Subject: [PATCH 22/41] remove some blank lines --- skimage/transform/_geometric.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 0af7ffb1..fe829109 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -52,7 +52,6 @@ class GeometricTransform(object): """ raise NotImplementedError() - def inverse(self, coords): """Apply inverse transformation. @@ -69,7 +68,6 @@ class GeometricTransform(object): """ raise NotImplementedError() - def __add__(self, other): """Combine this transformation with another. @@ -131,7 +129,6 @@ class ProjectiveTransform(GeometricTransform): 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. @@ -185,7 +182,6 @@ class ProjectiveTransform(GeometricTransform): "types.") - class AffineTransform(ProjectiveTransform): """2D affine transformation of the form:: @@ -560,4 +556,3 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, # The spline filters sometimes return results outside [0, 1], # so clip to ensure valid data return np.clip(mapped.squeeze(), 0, 1) - From 0f5d6153d24ca7cbe5dda1a1bac991cf862a566a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 08:58:37 +0200 Subject: [PATCH 23/41] fix bug in estimation of similarity transformation --- skimage/transform/_geometric.py | 34 +++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index fe829109..70584664 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -271,6 +271,40 @@ class SimilarityTransform(AffineTransform): shear=0, translation=translation) + 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, 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]]) + class PolynomialTransform(GeometricTransform): """2D transformation of the form:: From 87c57770ba51f9727f455ebe16fe0b2753182f1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 16:24:58 +0200 Subject: [PATCH 24/41] reimplement implicit parameter functionality of transformations --- skimage/transform/_geometric.py | 142 +++++++++++++++------- skimage/transform/tests/test_geometric.py | 39 +++++- 2 files changed, 134 insertions(+), 47 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 70584664..a47b51e4 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -30,11 +30,6 @@ def _stackcopy(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. @@ -99,6 +94,11 @@ class ProjectiveTransform(GeometricTransform): [0 1 20] [0 0 1 ]]. + Parameters + ---------- + matrix : 3x3 array, optional + Homogeneous transformation matrix. + """ _coefs = range(8) @@ -201,22 +201,30 @@ class AffineTransform(ProjectiveTransform): 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. + matrix : 3x3 array, optional + Homogeneous transformation matrix. """ _coefs = range(6) - def __init__(self, scale=None, rotation=None, shear=None, translation=None): - ProjectiveTransform.__init__(self) + def compose_implicit(self, scale=None, rotation=None, shear=None, + translation=None): + """Set the transformation matrix with the implicit 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 + + """ if scale is None: scale = (1, 1) if rotation is None: @@ -226,18 +234,35 @@ class AffineTransform(ProjectiveTransform): 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] + [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 SimilarityTransform(AffineTransform): +class SimilarityTransform(ProjectiveTransform): """2D similarity transformation of the form:: X = a0*x + b0*y + a1 = @@ -254,23 +279,11 @@ class SimilarityTransform(AffineTransform): 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 + matrix : 3x3 array, optional + Homogeneous transformation matrix. """ - 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) - def estimate(self, src, dst): """Set the transformation matrix with the explicit transformation parameters. @@ -305,6 +318,52 @@ class SimilarityTransform(AffineTransform): [b0, a0, b1], [ 0, 0, 1]]) + def compose_implicit(self, scale=None, rotation=None, translation=None): + """Set the transformation matrix with the implicit transformation + parameters. + + Parameters + ---------- + scale : float, optional + scale factor + rotation : float, optional + rotation angle in counter-clockwise direction + translation : (tx, ty) as array, list or tuple, optional + x, y translation parameters + + """ + if scale is None: + scale = (1, 1) + if rotation is None: + rotation = 0 + if translation is None: + translation = (0, 0) + + 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): + if math.cos(self.rotation) == 0: + # sin(self.rotation) == 1 + scale = self._matrix[0, 1] + else: + scale = self._matrix[0, 0] / math.cos(self.rotation) + return scale + + @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 PolynomialTransform(GeometricTransform): """2D transformation of the form:: @@ -449,7 +508,7 @@ def estimate_transform(ttype, src, dst, **kwargs): >>> tform = tf.estimate_transform('similarity', src, dst) - >>> tform.inverse(tform.forward(src)) # == src + >>> tform.inverse(tform(src)) # == src >>> # warp image using the estimated transformation >>> from skimage import data @@ -458,15 +517,12 @@ def estimate_transform(ttype, src, dst, **kwargs): >>> 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) + >>> tform2 = tf.SimilarityTransform() + >>> tform2.compose_implicit(scale=1.1, rotation=1, translation=(10, 20)) >>> # unite transformations, applied in order from left to right >>> tform3 = tform + tform2 - >>> tform3.forward(src) # == tform2.forward(tform.forward(src)) + >>> tform3(src) # == tform2(tform(src)) """ ttype = ttype.lower() diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 4430673c..f99bab7a 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -1,5 +1,5 @@ import numpy as np -from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_equal, assert_array_almost_equal from skimage.transform._geometric import _stackcopy from skimage.transform import (estimate_transform, SimilarityTransform, @@ -43,12 +43,26 @@ def test_similarity_estimation(): 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) + assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) + assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) #: over-determined tform = estimate_transform('similarity', SRC, DST) assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) + assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) + assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) +def test_similarity_implicit(): + tform = SimilarityTransform() + scale = 0.1 + rotation = 1 + translation = (1, 1) + tform.compose_implicit(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 @@ -61,6 +75,19 @@ def test_affine_estimation(): assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) +def test_affine_implicit(): + tform = AffineTransform() + scale = (0.1, 0.13) + rotation = 1 + shear = 0.1 + translation = (1, 1) + tform.compose_implicit(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) + + def test_projective(): #: exact solution tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) @@ -77,9 +104,13 @@ def test_polynomial(): def test_union(): - tform1 = SimilarityTransform(1, 0.3) - tform2 = SimilarityTransform(1, 0.6) - tform3 = SimilarityTransform(1, 0.9) + tform1 = SimilarityTransform() + tform1.compose_implicit(scale=0.1, rotation=0.3) + tform2 = SimilarityTransform() + tform2.compose_implicit(scale=0.1, rotation=0.9) + tform3 = SimilarityTransform() + tform3.compose_implicit(scale=0.1**2, rotation=0.3+0.9) + tform = tform1 + tform2 From be6bb0c809da448a846f6092cfc9085c74380f61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 16:38:25 +0200 Subject: [PATCH 25/41] combination of two transformations of the same type result in this type again --- skimage/transform/_geometric.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index a47b51e4..087e4ae6 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -176,7 +176,13 @@ class ProjectiveTransform(GeometricTransform): """ if isinstance(other, ProjectiveTransform): - return ProjectiveTransform(np.dot(other._matrix, self._matrix)) + # 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._matrix.dot(self._matrix)) else: raise TypeError("Cannot combine transformations of differing " "types.") From 2ae4dd4551801b92b740de19f8f383022a13bf28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 16:41:58 +0200 Subject: [PATCH 26/41] fix scale initialization in implicit composition of similarity --- skimage/transform/_geometric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 087e4ae6..1058315c 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -339,7 +339,7 @@ class SimilarityTransform(ProjectiveTransform): """ if scale is None: - scale = (1, 1) + scale = 1 if rotation is None: rotation = 0 if translation is None: From d9a88c95b5c5f19f374fcac2dd69da3ae8759a3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 22 Jul 2012 16:52:54 +0200 Subject: [PATCH 27/41] add doc for and restructure polynomial coefficients --- skimage/transform/_geometric.py | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 1058315c..c1c39b2b 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -374,24 +374,19 @@ class SimilarityTransform(ProjectiveTransform): 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 )) + 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 )) - TODO: Describe structure of coefficients. - Shall we store it as a (2, M) ndarray? + Parameters + ---------- + coeffs : 2xN array, optional + Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So, + a_ji is defined in `coeffs[0, :]` and b_ji in `coeffs[1, :]`. """ def __init__(self, coeffs=None): - """Create polynomial transformation. - - Parameters - ---------- - coeffs : array, optional - polynomial coefficients - - """ - self.coeffs = coeffs + self._coeffs = coeffs def estimate(self, src, dst, order): """Set the transformation matrix with the explicit transformation @@ -426,7 +421,7 @@ class PolynomialTransform(GeometricTransform): b = np.hstack([xd, yd]) - self.coeffs = np.linalg.lstsq(A, b)[0] + self._coeffs = np.linalg.lstsq(A, b)[0].reshape((2, u / 2)) def __call__(self, coords): """Apply forward transformation. @@ -444,7 +439,7 @@ class PolynomialTransform(GeometricTransform): """ x = coords[:, 0] y = coords[:, 1] - u = len(self.coeffs.ravel()) + 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) @@ -452,8 +447,8 @@ class PolynomialTransform(GeometricTransform): 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 + dst[:, 0] += self._coeffs[0, pidx] * x ** (j - i) * y ** i + dst[:, 1] += self._coeffs[1, pidx] * x ** (j - i) * y ** i pidx += 1 return dst From eb1e71114ca161615d7bce080ae9a2ccc0bf465a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 8 Aug 2012 19:13:45 +0200 Subject: [PATCH 28/41] fix and improve estimation of geometric transformation parameters Design matrix was not composed correctly as functional model was incorrect. Additionally estimation is now based on total least-squares method. --- skimage/transform/_geometric.py | 49 ++++++++++++++++++----- skimage/transform/tests/test_geometric.py | 6 +-- 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c1c39b2b..0c767cd7 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -133,6 +133,11 @@ class ProjectiveTransform(GeometricTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -148,7 +153,7 @@ class ProjectiveTransform(GeometricTransform): rows = src.shape[0] #: params: a0, a1, a2, b0, b1, b2, c0, c1 - A = np.zeros((rows * 2, 8)) + A = np.zeros((rows * 2, 9)) A[:rows, 0] = xs A[:rows, 1] = ys A[:rows, 2] = 1 @@ -159,14 +164,18 @@ class ProjectiveTransform(GeometricTransform): 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 coeffs - A = A[:, self._coefs] + A = A[:, self._coefs + [8]] - b = np.hstack([xd, yd]) + _, _, V = np.linalg.svd(A) H = np.zeros((3, 3)) - H.flat[self._coefs] = np.linalg.lstsq(A, b)[0] + # solution is eigen vector that corresponds to smallest eigen value + # and normed by c3 + H.flat[self._coefs + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 self._matrix = H @@ -294,6 +303,11 @@ class SimilarityTransform(ProjectiveTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -309,17 +323,20 @@ class SimilarityTransform(ProjectiveTransform): rows = src.shape[0] #: params: a0, a1, b0, b1 - A = np.zeros((rows * 2, 4)) + A = np.zeros((rows * 2, 5)) A[:rows, 0] = xs A[:rows, 2] = - ys A[:rows, 1] = 1 A[rows:, 2] = xs A[rows:, 0] = ys A[rows:, 3] = 1 + A[:rows, 4] = xd + A[rows:, 4] = yd - b = np.hstack([xd, yd]) + _, _, V = np.linalg.svd(A) + + a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1] - a0, a1, b0, b1 = np.linalg.lstsq(A, b)[0] self._matrix = np.array([[a0, -b0, a1], [b0, a0, b1], [ 0, 0, 1]]) @@ -392,6 +409,11 @@ class PolynomialTransform(GeometricTransform): """Set the transformation matrix with the explicit transformation parameters. + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source must match number of destination coordinates. + Parameters ---------- src : Nx2 array @@ -411,7 +433,7 @@ class PolynomialTransform(GeometricTransform): # number of unknown polynomial coefficients u = (order + 1) * (order + 2) - A = np.zeros((rows * 2, u)) + A = np.zeros((rows * 2, u + 1)) pidx = 0 for j in xrange(order + 1): for i in xrange(j + 1): @@ -419,9 +441,14 @@ class PolynomialTransform(GeometricTransform): A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i pidx += 1 - b = np.hstack([xd, yd]) + A[:rows, -1] = xd + A[rows:, -1] = yd - self._coeffs = np.linalg.lstsq(A, b)[0].reshape((2, u / 2)) + _, _, V = np.linalg.svd(A) + + coeffs = - V[-1, :-1] / V[-1, -1] + + self._coeffs = coeffs.reshape((2, u / 2)) def __call__(self, coords): """Apply forward transformation. @@ -473,7 +500,7 @@ 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. + with the total least-squares method. Number of source must match number of destination coordinates. diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index f99bab7a..c7305956 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -42,7 +42,6 @@ def test_similarity_estimation(): #: exact solution 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) assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) @@ -68,7 +67,6 @@ def test_affine_estimation(): #: exact solution 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_transform('affine', SRC, DST) @@ -91,10 +89,10 @@ def test_affine_implicit(): def test_projective(): #: exact solution tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) + assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :]) #: over-determined - tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) + tform = estimate_transform('projective', SRC, DST) assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) From 54452550f1d15a2a777fe7958020d6998d419fd1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 8 Aug 2012 19:19:18 +0200 Subject: [PATCH 29/41] fix incorrect comment --- skimage/transform/_geometric.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 0c767cd7..32d98d1c 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -173,8 +173,8 @@ class ProjectiveTransform(GeometricTransform): _, _, V = np.linalg.svd(A) H = np.zeros((3, 3)) - # solution is eigen vector that corresponds to smallest eigen value - # and normed by c3 + # solution is right singular vector that corresponds to smallest + # singular value and normed by c3 H.flat[self._coefs + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 @@ -335,6 +335,8 @@ class SimilarityTransform(ProjectiveTransform): _, _, V = np.linalg.svd(A) + # solution is right singular vector that corresponds to smallest + # singular value and normed by c3 a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1] self._matrix = np.array([[a0, -b0, a1], @@ -446,6 +448,8 @@ class PolynomialTransform(GeometricTransform): _, _, V = np.linalg.svd(A) + # solution is right singular vector that corresponds to smallest + # singular value and normed by c3 coeffs = - V[-1, :-1] / V[-1, -1] self._coeffs = coeffs.reshape((2, u / 2)) From e5ed1882d34a6234159f736bdbe5970fb3bbc9fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 07:51:58 +0200 Subject: [PATCH 30/41] fix and improve comments, doc strings, variable names for consistency reasons --- skimage/transform/_geometric.py | 99 ++++++++++++++++----------------- skimage/transform/_warps.py | 2 +- 2 files changed, 48 insertions(+), 53 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 32d98d1c..2db9c7d4 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -18,7 +18,7 @@ def _stackcopy(a, b): Notes ----- - Color images are stored as an ``MxNx3`` or ``MxNx4`` arrays. + Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays. """ if a.ndim == 3: @@ -36,12 +36,12 @@ class GeometricTransform(object): Parameters ---------- - coords : Nx2 array + coords : (N, 2) array source coordinates Returns ------- - coords : Nx2 array + coords : (N, 2) array transformed coordinates """ @@ -52,12 +52,12 @@ class GeometricTransform(object): Parameters ---------- - coords : Nx2 array + coords : (N, 2) array source coordinates Returns ------- - coords : Nx2 array + coords : (N, 2) array transformed coordinates """ @@ -78,17 +78,13 @@ class ProjectiveTransform(GeometricTransform): 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 - - :: + 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, - - :: + or, to translate x by 10 and y by 20:: [[1 0 10] [0 1 20] @@ -96,12 +92,12 @@ class ProjectiveTransform(GeometricTransform): Parameters ---------- - matrix : 3x3 array, optional + matrix : (3, 3) array, optional Homogeneous transformation matrix. """ - _coefs = range(8) + coeffs = range(8) def __init__(self, matrix=None): self._matrix = matrix @@ -136,13 +132,13 @@ class ProjectiveTransform(GeometricTransform): You can determine the over-, well- and under-determined parameters with the total least-squares method. - Number of source must match number of destination coordinates. + Number of source and destination coordinates must match. Parameters ---------- - src : Nx2 array + src : (N, 2) array source coordinates - dst : Nx2 array + dst : (N, 2) array destination coordinates """ @@ -152,7 +148,7 @@ class ProjectiveTransform(GeometricTransform): yd = dst[:, 1] rows = src.shape[0] - #: params: a0, a1, a2, b0, b1, b2, c0, c1 + # params: a0, a1, a2, b0, b1, b2, c0, c1 A = np.zeros((rows * 2, 9)) A[:rows, 0] = xs A[:rows, 1] = ys @@ -167,15 +163,15 @@ class ProjectiveTransform(GeometricTransform): A[:rows, 8] = xd A[rows:, 8] = yd - # Select relevant columns, depending on coeffs - A = A[:, self._coefs + [8]] + # Select relevant columns, depending on params + A = A[:, self.coeffs + [8]] _, _, V = np.linalg.svd(A) H = np.zeros((3, 3)) # solution is right singular vector that corresponds to smallest # singular value and normed by c3 - H.flat[self._coefs + [8]] = - V[-1, :-1] / V[-1, -1] + H.flat[self.coeffs + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 self._matrix = H @@ -216,12 +212,12 @@ class AffineTransform(ProjectiveTransform): Parameters ---------- - matrix : 3x3 array, optional + matrix : (3, 3) array, optional Homogeneous transformation matrix. """ - _coefs = range(6) + coeffs = range(6) def compose_implicit(self, scale=None, rotation=None, shear=None, translation=None): @@ -294,25 +290,24 @@ class SimilarityTransform(ProjectiveTransform): Parameters ---------- - matrix : 3x3 array, optional + matrix : (3, 3) array, optional Homogeneous transformation matrix. """ def estimate(self, src, dst): - """Set the transformation matrix with the explicit transformation - parameters. + """Set the transformation matrix with the explicit parameters. You can determine the over-, well- and under-determined parameters with the total least-squares method. - Number of source must match number of destination coordinates. + Number of source and destination coordinates must match. Parameters ---------- - src : Nx2 array + src : (N, 2) array source coordinates - dst : Nx2 array + dst : (N, 2) array destination coordinates """ @@ -322,7 +317,7 @@ class SimilarityTransform(ProjectiveTransform): yd = dst[:, 1] rows = src.shape[0] - #: params: a0, a1, b0, b1 + # params: a0, a1, b0, b1 A = np.zeros((rows * 2, 5)) A[:rows, 0] = xs A[:rows, 2] = - ys @@ -398,14 +393,14 @@ class PolynomialTransform(GeometricTransform): Parameters ---------- - coeffs : 2xN array, optional + params : (2, N) array, optional Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So, - a_ji is defined in `coeffs[0, :]` and b_ji in `coeffs[1, :]`. + a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`. """ - def __init__(self, coeffs=None): - self._coeffs = coeffs + def __init__(self, params=None): + self._params = params def estimate(self, src, dst, order): """Set the transformation matrix with the explicit transformation @@ -414,13 +409,13 @@ class PolynomialTransform(GeometricTransform): You can determine the over-, well- and under-determined parameters with the total least-squares method. - Number of source must match number of destination coordinates. + Number of source and destination coordinates must match. Parameters ---------- - src : Nx2 array + src : (N, 2) array source coordinates - dst : Nx2 array + dst : (N, 2) array destination coordinates order : int polynomial order (number of coefficients is order + 1) @@ -437,8 +432,8 @@ class PolynomialTransform(GeometricTransform): A = np.zeros((rows * 2, u + 1)) pidx = 0 - for j in xrange(order + 1): - for i in xrange(j + 1): + 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 @@ -450,36 +445,36 @@ class PolynomialTransform(GeometricTransform): # solution is right singular vector that corresponds to smallest # singular value and normed by c3 - coeffs = - V[-1, :-1] / V[-1, -1] + params = - V[-1, :-1] / V[-1, -1] - self._coeffs = coeffs.reshape((2, u / 2)) + self._params = params.reshape((2, u / 2)) def __call__(self, coords): """Apply forward transformation. Parameters ---------- - coords : Nx2 array + coords : (N, 2) array source coordinates Returns ------- - coords : Nx2 array + coords : (N, 2) array transformed coordinates """ x = coords[:, 0] y = coords[:, 1] - u = len(self._coeffs.ravel()) + 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 xrange(order + 1): - for i in xrange(j + 1): - dst[:, 0] += self._coeffs[0, pidx] * x ** (j - i) * y ** i - dst[:, 1] += self._coeffs[1, pidx] * x ** (j - i) * y ** i + 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 @@ -506,7 +501,7 @@ def estimate_transform(ttype, src, dst, **kwargs): You can determine the over-, well- and under-determined parameters with the total least-squares method. - Number of source must match number of destination coordinates. + Number of source and destination coordinates must match. Parameters ---------- @@ -573,14 +568,14 @@ def matrix_transform(coords, matrix): Parameters ---------- - coords : Nx2 array + coords : (N, 2) array x, y coordinates to transform - matrix : 3x3 array + matrix : (3, 3) array Homogeneous transformation matrix. Returns ------- - coords : Nx2 array + coords : (N, 2) array transformed coordinates """ @@ -596,7 +591,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, 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 + Inverse coordinate map. A function that transforms a (N, 2) 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 diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 347b8440..f09f7944 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -152,7 +152,7 @@ def homography(image, H, output_shape=None, order=1, """ import warnings warnings.warn('the homography function is deprecated; ' - 'use the `warp` and `tform` function instead', + 'use the `warp` and `ProjectiveTransform` class instead', category=DeprecationWarning) tform = ProjectiveTransform(H) From 5085ded9943cf6c3aa52b078eebd2856ed1c5ad9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 07:57:44 +0200 Subject: [PATCH 31/41] fix geometric transformation example after refactoring the module --- doc/examples/applications/plot_geometric.py | 23 +++++++++++---------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/applications/plot_geometric.py index 90fbce51..4e0a0cea 100644 --- a/doc/examples/applications/plot_geometric.py +++ b/doc/examples/applications/plot_geometric.py @@ -28,17 +28,17 @@ Geometric transformations can either be created using the explicit parameters """ #: create using explicit parameters -tform = tf.SimilarityTransformation() +tform = tf.SimilarityTransform() scale = 1 rotation = math.pi/2 translation = (0, 1) -tform.from_params(scale, rotation, translation) -print tform.matrix +tform.compose_implicit(scale, rotation, translation) +print tform._matrix #: create using transformation matrix -matrix = tform.matrix.copy() +matrix = tform._matrix.copy() matrix[1, 2] = 2 -tform2 = tf.SimilarityTransformation(matrix) +tform2 = tf.SimilarityTransform(matrix) """ These transformation objects can be used to forward and reverse transform @@ -46,8 +46,8 @@ coordinates between the source and destination coordinate systems: """ coord = [1, 0] -print tform2.forward(coord) -print tform2.reverse(tform.forward(coord)) +print tform2(coord) +print tform2.inverse(tform(coord)) """ Image warping @@ -57,11 +57,11 @@ Geometric transformations can also be used to warp images: """ text = data.text() -tform.from_params(1, math.pi/4, (text.shape[0] / 2, -100)) +tform.compose_implicit(1, math.pi/4, (text.shape[0] / 2, -100)) -# uses tform.reverse, alternatively use tf.warp(text, tform.reverse) +# uses tform.inverse, alternatively use tf.warp(text, tform.inverse) rotated = tf.warp(text, tform) -back_rotated = tf.warp(rotated, tform.forward) +back_rotated = tf.warp(rotated, tform) plt.figure(figsize=(8, 3)) plt.subplot(131) @@ -112,7 +112,8 @@ dst = np.array(( (300, 0) )) -tform3 = tf.estimate_transformation('projective', src, dst) +tform3 = tf.ProjectiveTransform() +tform3.estimate(src, dst) warped = tf.warp(text, tform3, output_shape=(50, 300)) plt.figure(figsize=(8, 3)) From 1bb896e48735b1e61541f99277d377128bf62ae9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 08:19:37 +0200 Subject: [PATCH 32/41] remove confusing comment --- skimage/transform/_geometric.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 2db9c7d4..f5ad71f4 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -170,7 +170,7 @@ class ProjectiveTransform(GeometricTransform): H = np.zeros((3, 3)) # solution is right singular vector that corresponds to smallest - # singular value and normed by c3 + # singular value H.flat[self.coeffs + [8]] = - V[-1, :-1] / V[-1, -1] H[2, 2] = 1 @@ -331,7 +331,7 @@ class SimilarityTransform(ProjectiveTransform): _, _, V = np.linalg.svd(A) # solution is right singular vector that corresponds to smallest - # singular value and normed by c3 + # singular value a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1] self._matrix = np.array([[a0, -b0, a1], @@ -444,7 +444,7 @@ class PolynomialTransform(GeometricTransform): _, _, V = np.linalg.svd(A) # solution is right singular vector that corresponds to smallest - # singular value and normed by c3 + # singular value params = - V[-1, :-1] / V[-1, -1] self._params = params.reshape((2, u / 2)) From 76931fa61be74581910d91f936564a5bc0a8335f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 09:21:49 +0200 Subject: [PATCH 33/41] handle composition of transformations with implicit parameters with __init__ --- skimage/transform/_geometric.py | 117 +++++++++++----------- skimage/transform/tests/test_geometric.py | 18 ++-- 2 files changed, 63 insertions(+), 72 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index f5ad71f4..d394e0cf 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -214,44 +214,41 @@ class AffineTransform(ProjectiveTransform): ---------- 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, optional + shear : float, optional + shear angle in counter-clockwise direction + translation : (tx, ty) as array, list or tuple, optional + translation parameters """ coeffs = range(6) - def compose_implicit(self, scale=None, rotation=None, shear=None, - translation=None): - """Set the transformation matrix with the implicit transformation - parameters. + def __init__(self, matrix=None, scale=None, rotation=None, shear=None, + translation=None): + params = (scale, rotation, shear, translation) + if matrix is not None: + self._matrix = matrix + elif any(param is not None for param in 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) - 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 - - """ - 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._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 + 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): @@ -292,9 +289,36 @@ class SimilarityTransform(ProjectiveTransform): ---------- matrix : (3, 3) array, optional Homogeneous transformation matrix. + scale : float, optional + scale factor + rotation : float, optional + rotation angle in counter-clockwise direction + translation : (tx, ty) as array, list or tuple, optional + x, y translation parameters """ + def __init__(self, matrix=None, scale=None, rotation=None, + translation=None): + params = (scale, rotation, translation) + if matrix is not None: + self._matrix = matrix + elif any(param is not None for param in params): + if scale is None: + scale = 1 + if rotation is None: + rotation = 0 + if translation is None: + translation = (0, 0) + + 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 + def estimate(self, src, dst): """Set the transformation matrix with the explicit parameters. @@ -338,35 +362,6 @@ class SimilarityTransform(ProjectiveTransform): [b0, a0, b1], [ 0, 0, 1]]) - def compose_implicit(self, scale=None, rotation=None, translation=None): - """Set the transformation matrix with the implicit transformation - parameters. - - Parameters - ---------- - scale : float, optional - scale factor - rotation : float, optional - rotation angle in counter-clockwise direction - translation : (tx, ty) as array, list or tuple, optional - x, y translation parameters - - """ - if scale is None: - scale = 1 - if rotation is None: - rotation = 0 - if translation is None: - translation = (0, 0) - - 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): if math.cos(self.rotation) == 0: diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index c7305956..b25b1aef 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -53,11 +53,11 @@ def test_similarity_estimation(): def test_similarity_implicit(): - tform = SimilarityTransform() scale = 0.1 rotation = 1 translation = (1, 1) - tform.compose_implicit(scale, rotation, translation) + tform = SimilarityTransform(scale=scale, rotation=rotation, + translation=translation) assert_array_almost_equal(tform.scale, scale) assert_array_almost_equal(tform.rotation, rotation) assert_array_almost_equal(tform.translation, translation) @@ -74,12 +74,12 @@ def test_affine_estimation(): def test_affine_implicit(): - tform = AffineTransform() scale = (0.1, 0.13) rotation = 1 shear = 0.1 translation = (1, 1) - tform.compose_implicit(scale, rotation, shear, translation) + tform = AffineTransform(scale=scale, rotation=rotation, shear=shear, + translation=translation) assert_array_almost_equal(tform.scale, scale) assert_array_almost_equal(tform.rotation, rotation) assert_array_almost_equal(tform.shear, shear) @@ -102,13 +102,9 @@ def test_polynomial(): def test_union(): - tform1 = SimilarityTransform() - tform1.compose_implicit(scale=0.1, rotation=0.3) - tform2 = SimilarityTransform() - tform2.compose_implicit(scale=0.1, rotation=0.9) - tform3 = SimilarityTransform() - tform3.compose_implicit(scale=0.1**2, rotation=0.3+0.9) - + tform1 = SimilarityTransform(scale=0.1, rotation=0.3) + tform2 = SimilarityTransform(scale=0.1, rotation=0.9) + tform3 = SimilarityTransform(scale=0.1 ** 2, rotation=0.3 + 0.9) tform = tform1 + tform2 From 8e242d0436a8e49038640a559b881bb3d1aad8e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 10:28:05 +0200 Subject: [PATCH 34/41] add mathematical description of estimation in doc strings --- skimage/transform/_geometric.py | 97 ++++++++++++++++++++++++++++++++- 1 file changed, 94 insertions(+), 3 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index d394e0cf..cb6392bc 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -77,8 +77,13 @@ class ProjectiveTransform(GeometricTransform): 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:: + :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] @@ -134,6 +139,41 @@ class ProjectiveTransform(GeometricTransform): 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 solutions 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 @@ -273,7 +313,7 @@ class AffineTransform(ProjectiveTransform): class SimilarityTransform(ProjectiveTransform): """2D similarity transformation of the form:: - X = a0*x + b0*y + a1 = + X = a0*x - b0*y + a1 = = m*x*cos(rotation) + m*y*sin(rotation) + a1 Y = b0*x + a0*y + b1 = @@ -327,6 +367,31 @@ class SimilarityTransform(ProjectiveTransform): Number of source and destination coordinates must match. + The transformation is defined as:: + + X = a0*x - b0*y + a1 + Y = b0*x + a0*y + b1 + + These equations can be transformed to the following form:: + + 0 = a0*x - b0*y + a1 - X + 0 = b0*x + a0*y + b1 - Y + + which exist for each set of corresponding points, so we have a set of + N * 2 equations. The coefficients appear linearly so we can write + A x = 0, where:: + + A = [[x 1 -y 0 -X] + [y 0 x 1 -Y] + ... + ... + ] + x.T = [a0 a1 b0 b1 c3] + + In case of total least-squares the solutions 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 @@ -406,6 +471,32 @@ class PolynomialTransform(GeometricTransform): 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 solutions 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 From ff1c8c4376c5a542bb3b44f0bbf707f40ef95a92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 9 Aug 2012 10:30:40 +0200 Subject: [PATCH 35/41] fix typo --- skimage/transform/_geometric.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index cb6392bc..efd11a29 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -160,7 +160,7 @@ class ProjectiveTransform(GeometricTransform): ] x.T = [a0 a1 a2 b0 b1 b2 c0 c1 c3] - In case of total least-squares the solutions of this homogeneous system + 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. @@ -388,7 +388,7 @@ class SimilarityTransform(ProjectiveTransform): ] x.T = [a0 a1 b0 b1 c3] - In case of total least-squares the solutions of this homogeneous system + 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. @@ -493,7 +493,7 @@ class PolynomialTransform(GeometricTransform): x.T = [a00 a10 a11 a20 a21 a22 ... ann b00 b10 b11 b20 b21 b22 ... bnn c3] - In case of total least-squares the solutions of this homogeneous system + 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. From a3ec8e04829dec1cdc969f7817fb30d073af348d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 10 Aug 2012 07:35:16 +0200 Subject: [PATCH 36/41] change inverse_map parameter handling of warp function --- skimage/transform/_geometric.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index efd11a29..c2bd7cb3 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -679,9 +679,8 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, inverse_map : transformation object, callable xy = f(xy, **kwargs) Inverse coordinate map. A function that transforms a (N, 2) 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. + coordinates in the *source image* (e.g. a transformation object or its + inverse). map_args : dict, optional Keyword arguments passed to `inverse_map`. output_shape : tuple (rows, cols) @@ -735,8 +734,6 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, # 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 From 724a931d4217175a85a7aea02c77a48ec7d222f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 10 Aug 2012 07:49:04 +0200 Subject: [PATCH 37/41] adapt geometric example script to new API and improve some expressions --- doc/examples/applications/plot_geometric.py | 81 ++++++++++----------- 1 file changed, 39 insertions(+), 42 deletions(-) diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/applications/plot_geometric.py index 4e0a0cea..337ecad7 100644 --- a/doc/examples/applications/plot_geometric.py +++ b/doc/examples/applications/plot_geometric.py @@ -25,24 +25,27 @@ affine, projective and polynomial. Geometric transformations can either be created using the explicit parameters (e.g. scale, shear, rotation and translation) or the transformation matrix: + +First we create a transformation using explicit parameters: """ -#: create using explicit parameters -tform = tf.SimilarityTransform() -scale = 1 -rotation = math.pi/2 -translation = (0, 1) -tform.compose_implicit(scale, rotation, translation) +tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2, + translation=(0, 1)) print tform._matrix -#: create using transformation matrix +""" +Alternatively you can define a transformation by the transformation matrix +itself: +""" + matrix = tform._matrix.copy() matrix[1, 2] = 2 tform2 = tf.SimilarityTransform(matrix) """ -These transformation objects can be used to forward and reverse transform -coordinates between the source and destination coordinate systems: +These transformation objects can then be used to apply forward and inverse +coordinate transformations between the source and destination coordinate +systems: """ coord = [1, 0] @@ -57,26 +60,22 @@ Geometric transformations can also be used to warp images: """ text = data.text() -tform.compose_implicit(1, math.pi/4, (text.shape[0] / 2, -100)) -# uses tform.inverse, alternatively use tf.warp(text, tform.inverse) +tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 4, + translation=(text.shape[0] / 2, -100)) + rotated = tf.warp(text, tform) -back_rotated = tf.warp(rotated, tform) +back_rotated = tf.warp(rotated, tform.inverse) -plt.figure(figsize=(8, 3)) -plt.subplot(131) -plt.imshow(text) -plt.axis('off') +fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) +fig.subplots_adjust(**margins) plt.gray() -plt.subplot(132) -plt.imshow(rotated) -plt.axis('off') -plt.gray() -plt.subplot(133) -plt.imshow(back_rotated) -plt.axis('off') -plt.gray() -plt.subplots_adjust(**margins) +ax1.imshow(text) +ax1.axis('off') +ax2.imshow(rotated) +ax2.axis('off') +ax3.imshow(back_rotated) +ax3.axis('off') """ .. image:: PLOT2RST.current_figure @@ -88,7 +87,8 @@ In addition to the basic functionality mentioned above you can also estimate the parameters of a geometric transformation using the least-squares method. This can amongst other things be used for image registration or rectification, -where you have a set of control points or homologous points in two images. +where you have a set of control points or homologous/corresponding points in two +images. Let's assume we want to recognize letters on a photograph which was not taken from the front but at a certain angle. In the simplest case of a plane paper @@ -100,33 +100,30 @@ the image so that the distortion is removed and then apply a matching algorithm: text = data.text() src = np.array(( - (155, 15), - (65, 40), - (260, 130), - (360, 95) -)) -dst = np.array(( (0, 0), (0, 50), (300, 50), (300, 0) )) +dst = np.array(( + (155, 15), + (65, 40), + (260, 130), + (360, 95) +)) tform3 = tf.ProjectiveTransform() tform3.estimate(src, dst) warped = tf.warp(text, tform3, output_shape=(50, 300)) -plt.figure(figsize=(8, 3)) -plt.subplot(211) -plt.imshow(text) -plt.plot(src[:, 0], src[:, 1], '.r') -plt.axis('off') +fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(8, 3)) +fig.subplots_adjust(**margins) plt.gray() -plt.subplot(212) -plt.imshow(warped) -plt.axis('off') -plt.gray() -plt.subplots_adjust(**margins) +ax1.imshow(text) +ax1.plot(dst[:, 0], dst[:, 1], '.r') +ax1.axis('off') +ax2.imshow(warped) +ax2.axis('off') """ .. image:: PLOT2RST.current_figure From 8ab0c9a32958e7a03946473c5dc6ea8a574765f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 10 Aug 2012 08:15:06 +0200 Subject: [PATCH 38/41] add test for shape of transformation matrix --- skimage/transform/_geometric.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c2bd7cb3..54865b02 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -105,6 +105,8 @@ class ProjectiveTransform(GeometricTransform): coeffs = range(8) def __init__(self, matrix=None): + if matrix is not None and matrix.shape != (3, 3): + raise ValueError("invalid shape of transformation matrix") self._matrix = matrix @property @@ -269,10 +271,12 @@ class AffineTransform(ProjectiveTransform): def __init__(self, matrix=None, scale=None, rotation=None, shear=None, translation=None): + if matrix is not None and matrix.shape != (3, 3): + raise ValueError("invalid shape of transformation matrix") + self._matrix = matrix + params = (scale, rotation, shear, translation) - if matrix is not None: - self._matrix = matrix - elif any(param is not None for param in params): + if any(param is not None for param in params): if scale is None: scale = (1, 1) if rotation is None: @@ -340,10 +344,12 @@ class SimilarityTransform(ProjectiveTransform): def __init__(self, matrix=None, scale=None, rotation=None, translation=None): + if matrix is not None and matrix.shape != (3, 3): + raise ValueError("invalid shape of transformation matrix") + self._matrix = matrix + params = (scale, rotation, translation) - if matrix is not None: - self._matrix = matrix - elif any(param is not None for param in params): + if any(param is not None for param in params): if scale is None: scale = 1 if rotation is None: @@ -460,6 +466,8 @@ class PolynomialTransform(GeometricTransform): """ def __init__(self, params=None): + if params is not None and params.shape[0] != 2: + raise ValueError("invalid shape of transformation parameters") self._params = params def estimate(self, src, dst, order): From c7cac1cb0ffef4d0627837a17c808ef33cff36c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 10 Aug 2012 08:16:05 +0200 Subject: [PATCH 39/41] fix, improve and extend test cases related to geometric transformations --- skimage/transform/tests/test_geometric.py | 77 ++++++++++++++++++----- skimage/transform/tests/test_warps.py | 4 +- 2 files changed, 62 insertions(+), 19 deletions(-) diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index b25b1aef..d3e7f20c 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -39,20 +39,26 @@ def test_stackcopy(): def test_similarity_estimation(): - #: exact solution + # exact solution tform = estimate_transform('similarity', SRC[:2, :], DST[:2, :]) assert_array_almost_equal(tform(SRC[:2, :]), DST[:2, :]) assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) - #: over-determined - tform = estimate_transform('similarity', SRC, DST) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) - assert_equal(tform._matrix[0, 0], tform._matrix[1, 1]) - assert_equal(tform._matrix[0, 1], - tform._matrix[1, 0]) + # over-determined + tform2 = estimate_transform('similarity', SRC, DST) + assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + assert_equal(tform2._matrix[0, 0], tform2._matrix[1, 1]) + assert_equal(tform2._matrix[0, 1], - tform2._matrix[1, 0]) + + # via estimate method + tform3 = SimilarityTransform() + tform3.estimate(SRC, DST) + assert_array_almost_equal(tform3._matrix, tform2._matrix) -def test_similarity_implicit(): +def test_similarity_init(): + # init with implicit parameters scale = 0.1 rotation = 1 translation = (1, 1) @@ -62,18 +68,30 @@ def test_similarity_implicit(): assert_array_almost_equal(tform.rotation, rotation) assert_array_almost_equal(tform.translation, translation) + # init with transformation matrix + tform2 = SimilarityTransform(tform._matrix) + assert_array_almost_equal(tform2.scale, scale) + assert_array_almost_equal(tform2.rotation, rotation) + assert_array_almost_equal(tform2.translation, translation) + def test_affine_estimation(): - #: exact solution + # exact solution tform = estimate_transform('affine', SRC[:3, :], DST[:3, :]) assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :]) - #: over-determined - tform = estimate_transform('affine', SRC, DST) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) + # over-determined + tform2 = estimate_transform('affine', SRC, DST) + assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + + # via estimate method + tform3 = AffineTransform() + tform3.estimate(SRC, DST) + assert_array_almost_equal(tform3._matrix, tform2._matrix) -def test_affine_implicit(): +def test_affine_init(): + # init with implicit parameters scale = (0.1, 0.13) rotation = 1 shear = 0.1 @@ -85,21 +103,46 @@ def test_affine_implicit(): assert_array_almost_equal(tform.shear, shear) assert_array_almost_equal(tform.translation, translation) + # init with transformation matrix + tform2 = AffineTransform(tform._matrix) + assert_array_almost_equal(tform2.scale, scale) + assert_array_almost_equal(tform2.rotation, rotation) + assert_array_almost_equal(tform2.shear, shear) + assert_array_almost_equal(tform2.translation, translation) -def test_projective(): - #: exact solution + +def test_projective_estimation(): + # exact solution tform = estimate_transform('projective', SRC[:4, :], DST[:4, :]) assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :]) - #: over-determined + # over-determined + tform2 = estimate_transform('projective', SRC, DST) + assert_array_almost_equal(tform2.inverse(tform2(SRC)), SRC) + + # via estimate method + tform3 = ProjectiveTransform() + tform3.estimate(SRC, DST) + assert_array_almost_equal(tform3._matrix, tform2._matrix) + + +def test_projective_init(): tform = estimate_transform('projective', SRC, DST) - assert_array_almost_equal(tform.inverse(tform(SRC)), SRC) + # init with transformation matrix + tform2 = ProjectiveTransform(tform._matrix) + assert_array_almost_equal(tform2._matrix, tform._matrix) -def test_polynomial(): +def test_polynomial_estimation(): + # over-determined tform = estimate_transform('polynomial', SRC, DST, order=10) assert_array_almost_equal(tform(SRC), DST, 6) + # via estimate method + tform2 = PolynomialTransform() + tform2.estimate(SRC, DST, order=10) + assert_array_almost_equal(tform2._params, tform._params) + def test_union(): tform1 = SimilarityTransform(scale=0.1, rotation=0.3) diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 047b477b..7c2e52f2 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -11,8 +11,8 @@ 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)) + theta = - np.pi / 2 + tform = SimilarityTransform(scale=1, rotation=theta, translation=(0, 4)) x90 = warp(x, tform, order=1) assert_array_almost_equal(x90, np.rot90(x)) From 08f4379e0e38fb3cd7bae927a3d5a32185636f1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 10 Aug 2012 08:19:48 +0200 Subject: [PATCH 40/41] add information about unit of angles to doc strings --- skimage/transform/_geometric.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 54865b02..19e3cb2d 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -259,9 +259,9 @@ class AffineTransform(ProjectiveTransform): scale : (sx, sy) as array, list or tuple, optional scale factors rotation : float, optional - rotation angle in counter-clockwise direction, optional + rotation angle in counter-clockwise direction as radians shear : float, optional - shear angle in counter-clockwise direction + shear angle in counter-clockwise direction as radians translation : (tx, ty) as array, list or tuple, optional translation parameters @@ -334,9 +334,9 @@ class SimilarityTransform(ProjectiveTransform): matrix : (3, 3) array, optional Homogeneous transformation matrix. scale : float, optional - scale factor + scale factor rotation : float, optional - rotation angle in counter-clockwise direction + rotation angle in counter-clockwise direction as radians translation : (tx, ty) as array, list or tuple, optional x, y translation parameters From 4541146561d43ce2ad4ce9c3236ff60241e2b5b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 12 Aug 2012 10:13:29 +0200 Subject: [PATCH 41/41] raise exception if matrix and implicit parameter arguments provided --- skimage/transform/_geometric.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 19e3cb2d..13bcc7f6 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -271,12 +271,16 @@ class AffineTransform(ProjectiveTransform): def __init__(self, matrix=None, scale=None, rotation=None, shear=None, translation=None): - if matrix is not None and matrix.shape != (3, 3): - raise ValueError("invalid shape of transformation matrix") self._matrix = matrix + params = any(param is not None + for param in (scale, rotation, shear, translation)) - params = (scale, rotation, shear, translation) - if any(param is not None for param in params): + 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 and matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + elif params: if scale is None: scale = (1, 1) if rotation is None: @@ -344,12 +348,16 @@ class SimilarityTransform(ProjectiveTransform): def __init__(self, matrix=None, scale=None, rotation=None, translation=None): - if matrix is not None and matrix.shape != (3, 3): - raise ValueError("invalid shape of transformation matrix") self._matrix = matrix + params = any(param is not None + for param in (scale, rotation, translation)) - params = (scale, rotation, translation) - if any(param is not None for param in params): + 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 and matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + elif params: if scale is None: scale = 1 if rotation is None: