wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-02 15:40:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #197 from ahojnnes/projection

Browse Source 
ENH: Add geometric transformations.
This commit is contained in:
Stefan van der Walt
2012-08-12 18:21:44 -07:00
parent a75b7ca7f5 4541146561
commit 4e50fd2798
11 changed files with 1225 additions and 324 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+1
View File
@@ -107,3 +107,4 @@
Polygon, circle and ellipse drawing functions
Adaptive thresholding
Implementation of Matlab's `regionprops`
Estimation of geometric transformation parameters
doc/examples/applications/plot_geometric.py
+132
View File
@@ -0,0 +1,132 @@
"""
===============================
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:
First we create a transformation using explicit parameters:
"""
tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2,
translation=(0, 1))
print tform._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 then be used to apply forward and inverse
coordinate transformations between the source and destination coordinate
systems:
"""
coord = [1, 0]
print tform2(coord)
print tform2.inverse(tform(coord))
"""
Image warping
=============
Geometric transformations can also be used to warp images:
"""
text = data.text()
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.inverse)
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3))
fig.subplots_adjust(**margins)
plt.gray()
ax1.imshow(text)
ax1.axis('off')
ax2.imshow(rotated)
ax2.axis('off')
ax3.imshow(back_rotated)
ax3.axis('off')
"""
.. 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/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
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((
(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))
fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(8, 3))
fig.subplots_adjust(**margins)
plt.gray()
ax1.imshow(text)
ax1.plot(dst[:, 0], dst[:, 1], '.r')
ax1.axis('off')
ax2.imshow(warped)
ax2.axis('off')
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
skimage/transform/__init__.py
+4 -3
View File
@@ -1,8 +1,9 @@
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 ._warp import warp
from ._warp_zoo import swirl
from ._geometric import (warp, estimate_transform,
SimilarityTransform, AffineTransform,
ProjectiveTransform, PolynomialTransform)
from ._warps import swirl, homography
skimage/transform/_geometric.py
+774
View File
@@ -0,0 +1,774 @@
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 ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
"""
if a.ndim == 3:
a[:] = b[:, :, np.newaxis]
else:
a[:] = b
class GeometricTransform(object):
"""Perform geometric transformations on a set of coordinates.
"""
def __call__(self, coords):
"""Apply forward transformation.
Parameters
----------
coords : (N, 2) array
source coordinates
Returns
-------
coords : (N, 2) array
transformed coordinates
"""
raise NotImplementedError()
def inverse(self, coords):
"""Apply inverse transformation.
Parameters
----------
coords : (N, 2) array
source coordinates
Returns
-------
coords : (N, 2) array
transformed coordinates
"""
raise NotImplementedError()
def __add__(self, other):
"""Combine this transformation with another.
"""
raise NotImplementedError()
class ProjectiveTransform(GeometricTransform):
"""Matrix transformation.
Apply a projective transformation (homography) on coordinates.
For each homogeneous coordinate :math:`\mathbf{x} = [x, y, 1]^T`, its
target position is calculated by multiplying with the given matrix,
:math:`H`, to give :math:`H \mathbf{x}`::
[[a0 a1 a2]
[b0 b1 b2]
[c0 c1 1 ]].
E.g., to rotate by theta degrees clockwise, the matrix should be::
[[cos(theta) -sin(theta) 0]
[sin(theta) cos(theta) 0]
[0 0 1]]
or, to translate x by 10 and y by 20::
[[1 0 10]
[0 1 20]
[0 0 1 ]].
Parameters
----------
matrix : (3, 3) array, optional
Homogeneous transformation matrix.
"""
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
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.
You can determine the over-, well- and under-determined parameters
with the total least-squares method.
Number of source and destination coordinates must match.
The transformation is defined as::
X = (a0*x + a1*y + a2) / (c0*x + c1*y + 1)
Y = (b0*x + b1*y + b2) / (c0*x + c1*y + 1)
These equations can be transformed to the following form::
0 = a0*x + a1*y + a2 - c0*x*X - c1*y*X - X
0 = b0*x + b1*y + b2 - c0*x*Y - c1*y*Y - Y
which exist for each set of corresponding points, so we have a set of
N * 2 equations. The coefficients appear linearly so we can write
A x = 0, where::
A = [[x y 1 0 0 0 -x*X -y*X -X]
[0 0 0 x y 1 -x*Y -y*Y -Y]
...
...
]
x.T = [a0 a1 a2 b0 b1 b2 c0 c1 c3]
In case of total least-squares the solution of this homogeneous system
of equations is the right singular vector of A which corresponds to the
smallest singular value normed by the coefficient c3.
In case of the affine transformation the coefficients c0 and c1 are 0.
Thus the system of equations is::
A = [[x y 1 0 0 0 -X]
[0 0 0 x y 1 -Y]
...
...
]
x.T = [a0 a1 a2 b0 b1 b2 c3]
Parameters
----------
src : (N, 2) array
source coordinates
dst : (N, 2) array
destination coordinates
"""
xs = src[:, 0]
ys = src[:, 1]
xd = dst[:, 0]
yd = dst[:, 1]
rows = src.shape[0]
# params: a0, a1, a2, b0, b1, b2, c0, c1
A = np.zeros((rows * 2, 9))
A[:rows, 0] = xs
A[:rows, 1] = ys
A[:rows, 2] = 1
A[:rows, 6] = - xd * xs
A[:rows, 7] = - xd * ys
A[rows:, 3] = xs
A[rows:, 4] = ys
A[rows:, 5] = 1
A[rows:, 6] = - yd * xs
A[rows:, 7] = - yd * ys
A[:rows, 8] = xd
A[rows:, 8] = yd
# Select relevant columns, depending on params
A = A[:, self.coeffs + [8]]
_, _, V = np.linalg.svd(A)
H = np.zeros((3, 3))
# solution is right singular vector that corresponds to smallest
# singular value
H.flat[self.coeffs + [8]] = - V[-1, :-1] / V[-1, -1]
H[2, 2] = 1
self._matrix = H
def __add__(self, other):
"""Combine this transformation with another.
"""
if isinstance(other, ProjectiveTransform):
# combination of the same types result in a transformation of this
# type again, otherwise use general projective transformation
if type(self) == type(other):
tform = self.__class__
else:
tform = ProjectiveTransform
return tform(other._matrix.dot(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
----------
matrix : (3, 3) array, optional
Homogeneous transformation matrix.
scale : (sx, sy) as array, list or tuple, optional
scale factors
rotation : float, optional
rotation angle in counter-clockwise direction as radians
shear : float, optional
shear angle in counter-clockwise direction as radians
translation : (tx, ty) as array, list or tuple, optional
translation parameters
"""
coeffs = range(6)
def __init__(self, matrix=None, scale=None, rotation=None, shear=None,
translation=None):
self._matrix = matrix
params = any(param is not None
for param in (scale, rotation, shear, translation))
if params and matrix is not None:
raise ValueError("You cannot specify the transformation matrix and "
"the implicit parameters at the same time.")
elif matrix is not None 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:
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
@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(ProjectiveTransform):
"""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
----------
matrix : (3, 3) array, optional
Homogeneous transformation matrix.
scale : float, optional
scale factor
rotation : float, optional
rotation angle in counter-clockwise direction as radians
translation : (tx, ty) as array, list or tuple, optional
x, y translation parameters
"""
def __init__(self, matrix=None, scale=None, rotation=None,
translation=None):
self._matrix = matrix
params = any(param is not None
for param in (scale, rotation, translation))
if params and matrix is not None:
raise ValueError("You cannot specify the transformation matrix and "
"the implicit parameters at the same time.")
elif matrix is not None and matrix.shape != (3, 3):
raise ValueError("Invalid shape of transformation matrix.")
elif 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.
You can determine the over-, well- and under-determined parameters
with the total least-squares method.
Number of source and destination coordinates must match.
The transformation is defined as::
X = a0*x - 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 solution of this homogeneous system
of equations is the right singular vector of A which corresponds to the
smallest singular value normed by the coefficient c3.
Parameters
----------
src : (N, 2) array
source coordinates
dst : (N, 2) array
destination coordinates
"""
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, 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
_, _, V = np.linalg.svd(A)
# solution is right singular vector that corresponds to smallest
# singular value
a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1]
self._matrix = np.array([[a0, -b0, a1],
[b0, a0, b1],
[ 0, 0, 1]])
@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::
X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i ))
Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i ))
Parameters
----------
params : (2, N) array, optional
Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So,
a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`.
"""
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):
"""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 and destination coordinates must match.
The transformation is defined as::
X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i ))
Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i ))
These equations can be transformed to the following form::
0 = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) - X
0 = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) - Y
which exist for each set of corresponding points, so we have a set of
N * 2 equations. The coefficients appear linearly so we can write
A x = 0, where::
A = [[1 x y x**2 x*y y**2 ... 0 ... 0 -X]
[0 ... 0 1 x y x**2 x*y y**2 -Y]
...
...
]
x.T = [a00 a10 a11 a20 a21 a22 ... ann
b00 b10 b11 b20 b21 b22 ... bnn c3]
In case of total least-squares the solution of this homogeneous system
of equations is the right singular vector of A which corresponds to the
smallest singular value normed by the coefficient c3.
Parameters
----------
src : (N, 2) array
source coordinates
dst : (N, 2) array
destination coordinates
order : int
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 + 1))
pidx = 0
for j in range(order + 1):
for i in range(j + 1):
A[:rows, pidx] = xs ** (j - i) * ys ** i
A[rows:, pidx + u / 2] = xs ** (j - i) * ys ** i
pidx += 1
A[:rows, -1] = xd
A[rows:, -1] = yd
_, _, V = np.linalg.svd(A)
# solution is right singular vector that corresponds to smallest
# singular value
params = - V[-1, :-1] / V[-1, -1]
self._params = params.reshape((2, u / 2))
def __call__(self, coords):
"""Apply forward transformation.
Parameters
----------
coords : (N, 2) array
source coordinates
Returns
-------
coords : (N, 2) array
transformed coordinates
"""
x = coords[:, 0]
y = coords[:, 1]
u = len(self._params.ravel())
# number of coefficients -> u = (order + 1) * (order + 2)
order = int((- 3 + math.sqrt(9 - 4 * (2 - u))) / 2)
dst = np.zeros(coords.shape)
pidx = 0
for j in range(order + 1):
for i in range(j + 1):
dst[:, 0] += self._params[0, pidx] * x ** (j - i) * y ** i
dst[:, 1] += self._params[1, pidx] * x ** (j - i) * y ** i
pidx += 1
return dst
def inverse(self, coords):
raise Exception(
'There is no explicit way to do the inverse polynomial '
'transformation. Instead, estimate the inverse transformation '
'parameters by exchanging source and destination coordinates,'
'then apply the forward transformation.')
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 total least-squares method.
Number of source and destination coordinates must match.
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(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
>>> 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(src) # == tform2(tform(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 : (N, 2) array
x, y coordinates to transform
matrix : (3, 3) array
Homogeneous transformation matrix.
Returns
-------
coords : (N, 2) array
transformed coordinates
"""
return ProjectiveTransform(matrix)(coords)
def 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 (N, 2) array of
``(x, y)`` coordinates in the *output image* into their corresponding
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)
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
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/_warp.py
-113
View File
@@ -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)
skimage/transform/_warp_zoo.py
-75
View File
@@ -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)
skimage/transform/project.py → skimage/transform/_warps.py
+76 -53
View File
@@ -1,14 +1,75 @@
"""Image projection.
"""
from ._geometric import warp, ProjectiveTransform
import numpy as np
from scipy.ndimage import interpolation as ndii
from ._warp import _stackcopy
__all__ = ['homography']
def _swirl_mapping(xy, center, rotation, strength, radius):
x, y = xy.T
x0, y0 = center
rho = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
eps = np.finfo(float).eps
# 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,
@@ -89,49 +150,11 @@ def homography(image, H, output_shape=None, order=1,
[ 0, 0, 0, 0, 0]], dtype=uint8)
"""
if image.ndim < 2:
raise ValueError("Input must have more than 1 dimension.")
import warnings
warnings.warn('the homography function is deprecated; '
'use the `warp` and `ProjectiveTransform` class instead',
category=DeprecationWarning)
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()
tform = ProjectiveTransform(H)
return warp(image, inverse_map=tform.inverse, output_shape=output_shape,
order=order, mode=mode, cval=cval)
skimage/transform/tests/test_geometric.py
+159
View File
@@ -0,0 +1,159 @@
import numpy as np
from numpy.testing import assert_equal, assert_array_almost_equal
from skimage.transform._geometric import _stackcopy
from skimage.transform import (estimate_transform, SimilarityTransform,
AffineTransform, ProjectiveTransform,
PolynomialTransform)
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_estimation():
# 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
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_init():
# init with implicit parameters
scale = 0.1
rotation = 1
translation = (1, 1)
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)
# 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
tform = estimate_transform('affine', SRC[:3, :], DST[:3, :])
assert_array_almost_equal(tform(SRC[:3, :]), DST[:3, :])
# 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_init():
# init with implicit parameters
scale = (0.1, 0.13)
rotation = 1
shear = 0.1
translation = (1, 1)
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)
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_estimation():
# exact solution
tform = estimate_transform('projective', SRC[:4, :], DST[:4, :])
assert_array_almost_equal(tform(SRC[:4, :]), DST[:4, :])
# 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)
# init with transformation matrix
tform2 = ProjectiveTransform(tform._matrix)
assert_array_almost_equal(tform2._matrix, tform._matrix)
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)
tform2 = SimilarityTransform(scale=0.1, rotation=0.9)
tform3 = SimilarityTransform(scale=0.1 ** 2, rotation=0.3 + 0.9)
tform = tform1 + tform2
assert_array_almost_equal(tform._matrix, tform3._matrix)
if __name__ == "__main__":
from numpy.testing import run_module_suite
run_module_suite()
skimage/transform/tests/test_project.py
-63
View File
@@ -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()
skimage/transform/tests/test_swirl.py
-17
View File
@@ -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()
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(scale=1, rotation=theta, translation=(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()