Files
scikit-image/skimage/measure/_polygon.py
T

89 lines
2.8 KiB
Python

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, :]