wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-29 21:09:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1 from stefanv/projection

Browse Source 
Refactor geometric transforms.
This commit is contained in:
Johannes Schönberger
2012-07-21 15:20:17 -07:00
parent 8e8e2b99a0 a130b8d2d9
commit c07a76bda7
6 changed files with 830 additions and 914 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/transform/__init__.py
+4 -3
View File
@@ -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
skimage/transform/_geometric.py
+563
View File
@@ -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)
skimage/transform/_warps.py
+160
View File
@@ -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)
skimage/transform/geometric.py
-764
View File
@@ -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)
skimage/transform/tests/test_geometric.py
+24 -147
View File
@@ -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__":
skimage/transform/tests/test_warps.py
+79
View File
@@ -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()