add polygon approximation algorithm

CONTRIBUTORS.txt

@@ -105,11 +105,8 @@
   Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks``
 
 - Johannes Schönberger
-  Polygon, circle and ellipse drawing functions
-  Adaptive thresholding
-  Implementation of Matlab's `regionprops`
-  Estimation of geometric transformation parameters
-  Local binary pattern texture classification
+  Drawing functions, adaptive thresholding, regionprops, geometric
+  transformations, LBPs, polygon approximations, and more.
 
 - Pavel Campr
   Fixes and tests for Histograms of Oriented Gradients.

skimage/measure/__init__.py

@@ -1,3 +1,4 @@
 from .find_contours import find_contours
 from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
+from ._polygon import approximate_polygon

skimage/measure/_polygon.py

@@ -0,0 +1,88 @@
+import numpy as np
+
+
+def approximate_polygon(coords, tolerance):
+    """Approximate a polygonal chain with the specified tolerance.
+
+    It is based on the Douglas-Peucker algorithm.
+
+    Parameters
+    ----------
+    coords : (N, 2) array
+        Coordinate array.
+    tolerance : float
+        Maximum distance from original points of polygon to approximated
+        polygonal chain. If tolerance is 0, the original coordinate array
+        is returned.
+
+    Returns
+    -------
+    coords : (M, 2) array
+        Approximated polygonal chain where M <= N.
+
+    References
+    ----------
+    .. [1] http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
+    """
+    if tolerance == 0:
+        return coords
+
+    chain = np.zeros(coords.shape[0], 'bool')
+    # pre-allocate distance array for all points
+    dists = np.zeros(coords.shape[0])
+    chain[0] = True
+    chain[-1] = True
+    pos_stack = [(0, chain.shape[0] - 1)]
+
+    while 1:
+        start, end = pos_stack.pop()
+        # determine properties of current line segment
+        r0, c0 = coords[start, :]
+        r1, c1 = coords[end, :]
+        dr = r1 - r0
+        dc = c1 - c0
+        segment_angle = - np.arctan2(dr, dc)
+        segment_dist = c0 * np.sin(segment_angle) + r0 * np.cos(segment_angle)
+
+        # select points in-between line segment
+        segment_coords = coords[start + 1:end, :]
+        segment_dists = dists[start + 1:end]
+
+
+        # check whether to take perpendicular or euclidean distance with
+        # inner product of vectors
+
+        # vectors from points -> start and end
+        dr0 = segment_coords[:, 0] - r0
+        dc0 = segment_coords[:, 1] - c0
+        dr1 = segment_coords[:, 0] - r1
+        dc1 = segment_coords[:, 1] - c1
+        # vectors points -> start and end projected on start -> end vector
+        projected_lengths0 = dr0 * dr + dc0 * dc
+        projected_lengths1 = - dr1 * dr - dc1 * dc
+        perp = np.logical_and(projected_lengths0 > 0,
+                              projected_lengths1 > 0)
+        eucl = np.logical_not(perp)
+        segment_dists[perp] = np.abs(
+            segment_coords[perp, 0] * np.cos(segment_angle)
+            + segment_coords[perp, 1] * np.sin(segment_angle)
+            - segment_dist
+        )
+        segment_dists[eucl] = np.minimum(
+            # distance to start point
+            np.sqrt(dc0[eucl] ** 2 + dr0[eucl] ** 2),
+            # distance to end point
+            np.sqrt(dc1[eucl] ** 2 + dr1[eucl] ** 2)
+        )
+
+        if np.any(segment_dists > tolerance):
+            # select point with maximum distance to line
+            new_end = start + np.argmax(segment_dists) + 1
+            pos_stack.append((new_end, end))
+            pos_stack.append((start, new_end))
+            chain[new_end] = True
+
+        if len(pos_stack) == 0:
+            break
+
+    return coords[chain, :]

skimage/measure/tests/test_polygon.py

@@ -0,0 +1,27 @@
+import numpy as np
+from skimage.measure import approximate_polygon
+
+
+def test_approximate_polygon():
+    square = np.array([
+        [0, 0], [0, 1], [0, 2], [0, 3],
+        [1, 3], [2, 3], [3, 3],
+        [3, 2], [3, 1], [3, 0],
+        [2, 0], [1, 0], [0, 0]
+    ])
+
+    out = approximate_polygon(square, 0.1)
+    np.testing.assert_array_equal(out, square[(0, 3, 6, 9, 12), :])
+
+    out = approximate_polygon(square, 2.2)
+    np.testing.assert_array_equal(out, square[(0, 6, 12), :])
+
+    out = approximate_polygon(square[(0, 1, 3, 4, 5, 6, 7, 9, 11, 12), :], 0.1)
+    np.testing.assert_array_equal(out, square[(0, 3, 6, 9, 12), :])
+
+    out = approximate_polygon(square, -1)
+    np.testing.assert_array_equal(out, square)
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()