From 7983f354b80bc7cdc8c0ae39cbebd8ca3ac07060 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:17:46 +0200 Subject: [PATCH 01/12] Apply PEP8 and improve docs --- skimage/transform/_geometric.py | 164 ++++++++++++++++---------------- 1 file changed, 84 insertions(+), 80 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 030197da..54764742 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -303,6 +303,85 @@ 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. + + Parameters + ---------- + TODO + + """ + + def __init__(self): + self.tesselation = None + self.affines = [] + self._matrix = None + + 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. + + """ + + # 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) + + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : (N, 2) array + source coordinates + + Returns + ------- + coords : (N, 2) array + Transformed coordinates. + + """ + + out = - 1 * np.ones((coords.shape[0], 2)) + + for index, pt in enumerate(coords): + # determine which triangle contains the point + simplex_index = self.tesselation.find_simplex(pt) + + if simplex_index == - 1: + # this point is outside the hull of the control points + out[index, 0] = - 1 + out[index, 1] = - 1 + continue + + # calculate affine transformed position + affine = self.affines[simplex_index] + dst_pos = affine(pt) + out[index, 0] = dst_pos[0][0] + out[index, 1] = dst_pos[0][1] + + return out + + class SimilarityTransform(ProjectiveTransform): """2D similarity transformation of the form:: @@ -580,90 +659,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 +673,7 @@ HOMOGRAPHY_TRANSFORMS = ( ProjectiveTransform ) + def estimate_transform(ttype, src, dst, **kwargs): """Estimate 2D geometric transformation parameters. @@ -681,7 +684,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 +693,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) From c1f3515e0065e5aa2af51def5d5300a4aafef8f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:29:19 +0200 Subject: [PATCH 02/12] Speed up transformation of piecewise-affine --- skimage/transform/_geometric.py | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 54764742..6d64ec5c 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -361,23 +361,19 @@ class PiecewiseAffineTransform(ProjectiveTransform): """ - out = - 1 * np.ones((coords.shape[0], 2)) + out = np.empty_like(coords) - for index, pt in enumerate(coords): - # determine which triangle contains the point - simplex_index = self.tesselation.find_simplex(pt) + simplex = self.tesselation.find_simplex(coords) - if simplex_index == - 1: - # this point is outside the hull of the control points - out[index, 0] = - 1 - out[index, 1] = - 1 - continue + out[simplex == -1, :] = -1 - # calculate affine transformed position - affine = self.affines[simplex_index] - dst_pos = affine(pt) - out[index, 0] = dst_pos[0][0] - out[index, 1] = dst_pos[0][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 From 33151f4349bdc899b817942b255a7b15c49605f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:40:45 +0200 Subject: [PATCH 03/12] Add missing doc string for inverse transform --- skimage/transform/_geometric.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 6d64ec5c..ab530852 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -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): From d541ebac17d764f4ccaac440322c8d371a20edf7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:41:16 +0200 Subject: [PATCH 04/12] Implement inverse transformation for piecewise-affine --- skimage/transform/_geometric.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index ab530852..e4913d4b 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -390,6 +390,37 @@ class PiecewiseAffineTransform(ProjectiveTransform): return out + def inverse(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : (N, 2) array + Source coordinates. + + Returns + ------- + coords : (N, 2) array + Transformed coordinates. + + """ + + out = np.empty_like(coords) + + simplex = self.tesselation.find_simplex(coords) + + 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.inverse(coords[index_mask, :]) + + return out + class SimilarityTransform(ProjectiveTransform): """2D similarity transformation of the form:: From 3f7d962206b175265054c284dbe59d3af5a7e783 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:41:53 +0200 Subject: [PATCH 05/12] Remove unnecessary params section in doc string --- skimage/transform/_geometric.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index e4913d4b..6c0afbdb 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -324,16 +324,11 @@ class PiecewiseAffineTransform(ProjectiveTransform): 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.tesselation = None self.affines = [] - self._matrix = None def estimate(self, src, dst): """Set the control points with which to perform the piecewise mapping. From 65879a2cdef723a2ec870dd7d62d8e1fde01c17d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:46:14 +0200 Subject: [PATCH 06/12] Fix fast warping of images --- skimage/transform/_geometric.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 6c0afbdb..cfa12033 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -954,9 +954,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: From f79012fe2264642085aee80a176b27bf12602808 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 21:49:14 +0200 Subject: [PATCH 07/12] Add note about coordinates outside of mesh --- skimage/transform/_geometric.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index cfa12033..e2975def 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -357,10 +357,12 @@ class PiecewiseAffineTransform(ProjectiveTransform): def __call__(self, coords): """Apply forward transformation. + Coordinates outside of the mesh will be set to `- 1`. + Parameters ---------- coords : (N, 2) array - source coordinates + Source coordinates. Returns ------- @@ -388,6 +390,8 @@ class PiecewiseAffineTransform(ProjectiveTransform): def inverse(self, coords): """Apply inverse transformation. + Coordinates outside of the mesh will be set to `- 1`. + Parameters ---------- coords : (N, 2) array From 677478873c2d2e6a439b0bf443e3524f6d3cba41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 22:04:57 +0200 Subject: [PATCH 08/12] Fix inverse piecewise affine --- skimage/transform/_geometric.py | 41 +++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index e2975def..f81e6e22 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -327,8 +327,10 @@ class PiecewiseAffineTransform(ProjectiveTransform): """ def __init__(self): - self.tesselation = None + 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. @@ -344,16 +346,27 @@ class PiecewiseAffineTransform(ProjectiveTransform): """ + # forward piecewise affine # triangulate input positions into mesh - self.tesselation = spatial.Delaunay(src) - + self._tesselation = spatial.Delaunay(src) # find affine mapping from source positions to destination self.affines = [] - for tri in self.tesselation.vertices: + 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. @@ -371,13 +384,15 @@ class PiecewiseAffineTransform(ProjectiveTransform): """ - out = np.empty_like(coords) + out = np.empty_like(coords, np.double) - simplex = self.tesselation.find_simplex(coords) + # 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)): + for index in range(len(self._tesselation.vertices)): # affine transform for triangle affine = self.affines[index] # all coordinates within triangle @@ -404,19 +419,21 @@ class PiecewiseAffineTransform(ProjectiveTransform): """ - out = np.empty_like(coords) + out = np.empty_like(coords, np.double) - simplex = self.tesselation.find_simplex(coords) + # 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.tesselation.vertices)): + for index in range(len(self._inverse_tesselation.vertices)): # affine transform for triangle - affine = self.affines[index] + affine = self.inverse_affines[index] # all coordinates within triangle index_mask = simplex == index - out[index_mask, :] = affine.inverse(coords[index_mask, :]) + out[index_mask, :] = affine(coords[index_mask, :]) return out From 05ae17e0d5cc999e88a3e73993c495b72e2e3383 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 22:05:12 +0200 Subject: [PATCH 09/12] Add test cases for piecewise affine --- skimage/transform/tests/test_geometric.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index 57e95d46..d3a75f7b 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -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, :]) From e16a6d32dc4919db711bc5963601179aa4e3019c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 22:07:10 +0200 Subject: [PATCH 10/12] Fix typo in polygon example script doc --- doc/examples/plot_polygon.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_polygon.py b/doc/examples/plot_polygon.py index 5b745fed..05ca5359 100644 --- a/doc/examples/plot_polygon.py +++ b/doc/examples/plot_polygon.py @@ -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. From 53edb67dd3a7a88d9133a85cd2fe1a9d22a3a436 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 22:49:09 +0200 Subject: [PATCH 11/12] Add example script for piecewise affine transform --- doc/examples/plot_piecewise_affine.py | 38 +++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 doc/examples/plot_piecewise_affine.py diff --git a/doc/examples/plot_piecewise_affine.py b/doc/examples/plot_piecewise_affine.py new file mode 100644 index 00000000..9db0e402 --- /dev/null +++ b/doc/examples/plot_piecewise_affine.py @@ -0,0 +1,38 @@ +""" +=============================== +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, 20) +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) + +output_shape = (image.shape[0] - 1.5 * 50, image.shape[1]) +out = warp(image, tform, output_shape=output_shape) + +plt.imshow(out) +plt.show() From aff606998eccb328a48323f79d26d6c96ad4900a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 3 Sep 2012 23:06:59 +0200 Subject: [PATCH 12/12] Add mesh points to plot --- doc/examples/plot_piecewise_affine.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_piecewise_affine.py b/doc/examples/plot_piecewise_affine.py index 9db0e402..2dcbd9f1 100644 --- a/doc/examples/plot_piecewise_affine.py +++ b/doc/examples/plot_piecewise_affine.py @@ -16,7 +16,7 @@ image = data.lena() rows, cols = image.shape[0], image.shape[1] src_cols = np.linspace(0, cols, 20) -src_rows = np.linspace(0, rows, 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] @@ -31,8 +31,11 @@ dst = np.vstack([dst_cols, dst_rows]).T tform = PiecewiseAffineTransform() tform.estimate(src, dst) -output_shape = (image.shape[0] - 1.5 * 50, image.shape[1]) -out = warp(image, tform, output_shape=output_shape) +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()