diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 531ab842..4abe9256 100644 --- a/CONTRIBUTORS.txt +++ b/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. diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index f5109698..d7397b2f 100755 --- a/skimage/measure/__init__.py +++ b/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 \ No newline at end of file diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py new file mode 100644 index 00000000..0b80599a --- /dev/null +++ b/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, :] diff --git a/skimage/measure/tests/test_polygon.py b/skimage/measure/tests/test_polygon.py new file mode 100644 index 00000000..e539f69c --- /dev/null +++ b/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()