wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-23 13:10:18 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #297 from ahojnnes/piecewise-affine

Browse Source 
ENH: Improve piecewise affine transforms.
This commit is contained in:
Stefan van der Walt
2012-09-03 14:32:38 -07:00
parent 87a3025196 aff606998e
commit ef5e31d03d
4 changed files with 197 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_piecewise_affine.py
+41
View File
@@ -0,0 +1,41 @@
"""
===============================
Piecewise Affine Transformation
===============================
This example shows how to use the Piecewise Affine Transformation.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import PiecewiseAffineTransform, warp
from skimage import data
image = data.lena()
rows, cols = image.shape[0], image.shape[1]
src_cols = np.linspace(0, cols, 20)
src_rows = np.linspace(0, rows, 10)
src_rows, src_cols = np.meshgrid(src_rows, src_cols)
src = np.dstack([src_cols.flat, src_rows.flat])[0]
# add sinusoidal oscillation to row coordinates
dst_rows = src[:, 1] - np.sin(np.linspace(0, 3 * np.pi, src.shape[0])) * 50
dst_cols = src[:, 0]
dst_rows *= 1.5
dst_rows -= 1.5 * 50
dst = np.vstack([dst_cols, dst_rows]).T
tform = PiecewiseAffineTransform()
tform.estimate(src, dst)
out_rows = image.shape[0] - 1.5 * 50
out_cols = cols
out = warp(image, tform, output_shape=(out_rows, out_cols))
plt.imshow(out)
plt.plot(tform.inverse(src)[:, 0], tform.inverse(src)[:, 1], '.b')
plt.axis((0, out_cols, out_rows, 0))
plt.show()
doc/examples/plot_polygon.py
+3 -3
View File
@@ -1,7 +1,7 @@
"""
===================================
Approximatea and subdivide polygons
===================================
==================================
Approximate and subdivide polygons
==================================
This example shows how to approximate (Douglas-Peucker algorithm) and subdivide
(B-Splines) polygonal chains.
skimage/transform/_geometric.py
+143 -82
View File
@@ -111,6 +111,19 @@ class ProjectiveTransform(GeometricTransform):
return self._apply_mat(coords, self._matrix)
def inverse(self, coords):
"""Apply inverse transformation.
Parameters
----------
coords : (N, 2) array
Source coordinates.
Returns
-------
coords : (N, 2) array
Transformed coordinates.
"""
return self._apply_mat(coords, self._inv_matrix)
def estimate(self, src, dst):
@@ -303,6 +316,128 @@ class AffineTransform(ProjectiveTransform):
return self._matrix[0:2, 2]
class PiecewiseAffineTransform(ProjectiveTransform):
"""2D piecewise affine transformation.
Control points are used to define the mapping. The transform is based on
a Delaunay triangulation of the points to form a mesh. Each triangle is
used to find a local affine transform.
"""
def __init__(self):
self._tesselation = None
self._inverse_tesselation = None
self.affines = []
self.inverse_affines = []
def estimate(self, src, dst):
"""Set the control points with which to perform the piecewise mapping.
Number of source and destination coordinates must match.
Parameters
----------
src : (N, 2) array
Source coordinates.
dst : (N, 2) array
Destination coordinates.
"""
# forward piecewise affine
# triangulate input positions into mesh
self._tesselation = spatial.Delaunay(src)
# find affine mapping from source positions to destination
self.affines = []
for tri in self._tesselation.vertices:
affine = AffineTransform()
affine.estimate(src[tri, :], dst[tri, :])
self.affines.append(affine)
# inverse piecewise affine
# triangulate input positions into mesh
self._inverse_tesselation = spatial.Delaunay(dst)
# find affine mapping from source positions to destination
self.inverse_affines = []
for tri in self._inverse_tesselation.vertices:
affine = AffineTransform()
affine.estimate(dst[tri, :], src[tri, :])
self.inverse_affines.append(affine)
def __call__(self, coords):
"""Apply forward transformation.
Coordinates outside of the mesh will be set to `- 1`.
Parameters
----------
coords : (N, 2) array
Source coordinates.
Returns
-------
coords : (N, 2) array
Transformed coordinates.
"""
out = np.empty_like(coords, np.double)
# determine triangle index for each coordinate
simplex = self._tesselation.find_simplex(coords)
# coordinates outside of mesh
out[simplex == -1, :] = -1
for index in range(len(self._tesselation.vertices)):
# affine transform for triangle
affine = self.affines[index]
# all coordinates within triangle
index_mask = simplex == index
out[index_mask, :] = affine(coords[index_mask, :])
return out
def inverse(self, coords):
"""Apply inverse transformation.
Coordinates outside of the mesh will be set to `- 1`.
Parameters
----------
coords : (N, 2) array
Source coordinates.
Returns
-------
coords : (N, 2) array
Transformed coordinates.
"""
out = np.empty_like(coords, np.double)
# determine triangle index for each coordinate
simplex = self._inverse_tesselation.find_simplex(coords)
# coordinates outside of mesh
out[simplex == -1, :] = -1
for index in range(len(self._inverse_tesselation.vertices)):
# affine transform for triangle
affine = self.inverse_affines[index]
# all coordinates within triangle
index_mask = simplex == index
out[index_mask, :] = affine(coords[index_mask, :])
return out
class SimilarityTransform(ProjectiveTransform):
"""2D similarity transformation of the form::
@@ -580,90 +715,13 @@ class PolynomialTransform(GeometricTransform):
'parameters by exchanging source and destination coordinates,'
'then apply the forward transformation.')
class PiecewiseAffineTransform(ProjectiveTransform):
"""2D piecewise affine transformation.
Control points are used to define the mapping. The transform is based on
a Delaunay triangulation of the points to form a mesh. Each triangle is
used to find a local affine transform.
Parameters
----------
TODO
"""
def __init__(self):
self.tess = None
self.triAffines = []
self._matrix = None
def estimate(self, src, dst):
"""Set the control points with which to perform the piecewise affine mapping.
Number of source and destination coordinates must match.
Parameters
----------
src : (N, 2) array
Source coordinates.
dst : (N, 2) array
Destination coordinates.
"""
#Triangulate input positions into mesh
self.tess = spatial.Delaunay(src)
#Find affine mapping from source positions to destination
self.triAffines = []
for tri in self.tess.vertices:
affine = AffineTransform()
affine.estimate(src[tri,:], dst[tri,:])
self.triAffines.append(affine)
def __call__(self, coords):
"""Apply forward transformation.
Parameters
----------
coords : (N, 2) array
source coordinates
Returns
-------
coords : (N, 2) array
Transformed coordinates.
"""
out = np.ones((coords.shape[0], 2)) * -1
for ptNum, pt in enumerate(coords):
#Determine which triangle contains the point
simplexIndex = self.tess.find_simplex(pt)
if simplexIndex == -1:
#This point is outside the hull of the control points
out[ptNum,0] = -1
out[ptNum,1] = -1
continue
#Calculate affine transformed position
affine = self.triAffines[simplexIndex]
destPos = affine(pt)
out[ptNum,0] = destPos[0][0]
out[ptNum,1] = destPos[0][1]
return out
TRANSFORMS = {
'similarity': SimilarityTransform,
'affine': AffineTransform,
'piecewise-affine': PiecewiseAffineTransform,
'projective': ProjectiveTransform,
'polynomial': PolynomialTransform,
'piecewiseaffine': PiecewiseAffineTransform,
}
HOMOGRAPHY_TRANSFORMS = (
SimilarityTransform,
@@ -671,6 +729,7 @@ HOMOGRAPHY_TRANSFORMS = (
ProjectiveTransform
)
def estimate_transform(ttype, src, dst, **kwargs):
"""Estimate 2D geometric transformation parameters.
@@ -681,7 +740,8 @@ def estimate_transform(ttype, src, dst, **kwargs):
Parameters
----------
ttype : {'similarity', 'affine', 'piecewiseaffine', 'projective', 'polynomial'}
ttype : {'similarity', 'affine', 'piecewise-affine', 'projective', \
'polynomial'}
Type of transform.
kwargs : array or int
Function parameters (src, dst, n, angle)::
@@ -689,7 +749,7 @@ def estimate_transform(ttype, src, dst, **kwargs):
NAME / TTYPE FUNCTION PARAMETERS
'similarity' `src, `dst`
'affine' `src, `dst`
'piecewiseaffine' `src, `dst`
'piecewise-affine' `src, `dst`
'projective' `src, `dst`
'polynomial' `src, `dst`, `order` (polynomial order)
@@ -915,9 +975,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
# use fast Cython version for specific interpolation orders
if order in range(4) and not map_args:
matrix = None
if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
if inverse_map in HOMOGRAPHY_TRANSFORMS:
matrix = inverse_map._matrix
elif inverse_map.__name__ == 'inverse' \
elif hasattr(inverse_map, '__name__') \
and inverse_map.__name__ == 'inverse' \
and inverse_map.im_class in HOMOGRAPHY_TRANSFORMS:
matrix = np.linalg.inv(inverse_map.im_self._matrix)
if matrix is not None:
skimage/transform/tests/test_geometric.py
+10 -1
View File
@@ -3,7 +3,8 @@ 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)
ProjectiveTransform, PolynomialTransform,
PiecewiseAffineTransform)
SRC = np.array([
@@ -110,6 +111,14 @@ def test_affine_init():
assert_array_almost_equal(tform2.translation, translation)
def test_piecewise_affine():
tform = PiecewiseAffineTransform()
tform.estimate(SRC, DST)
# make sure each single affine transform is exactly estimated
assert_array_almost_equal(tform(SRC), DST)
assert_array_almost_equal(tform.inverse(DST), SRC)
def test_projective_estimation():
# exact solution
tform = estimate_transform('projective', SRC[:4, :], DST[:4, :])