"""
rdp
~~~
Python implementation of the Ramer-Douglas-Peucker algorithm.
:copyright: 2014-2016 Fabian Hirschmann <fabian@hirschmann.email>
:license: MIT, see LICENSE.txt for more details.
"""
from math import sqrt
from functools import partial
import numpy as np
import sys
if sys.version_info[0] >= 3:
xrange = range
[docs]def pldist(point, start, end):
"""
Calculates the distance from ``point`` to the line given
by the points ``start`` and ``end``.
:param point: a point
:type point: numpy array
:param start: a point of the line
:type start: numpy array
:param end: another point of the line
:type end: numpy array
"""
if np.all(np.equal(start, end)):
return np.linalg.norm(point - start)
return np.divide(
np.abs(np.linalg.norm(np.cross(end - start, start - point))),
np.linalg.norm(end - start))
[docs]def rdp_rec(M, epsilon, dist=pldist):
"""
Simplifies a given array of points.
Recursive version.
:param M: an array
:type M: numpy array
:param epsilon: epsilon in the rdp algorithm
:type epsilon: float
:param dist: distance function
:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
"""
dmax = 0.0
index = -1
for i in xrange(1, M.shape[0]):
d = dist(M[i], M[0], M[-1])
if d > dmax:
index = i
dmax = d
if dmax > epsilon:
r1 = rdp_rec(M[:index + 1], epsilon, dist)
r2 = rdp_rec(M[index:], epsilon, dist)
return np.vstack((r1[:-1], r2))
else:
return np.vstack((M[0], M[-1]))
def _rdp_iter(M, start_index, last_index, epsilon, dist=pldist):
stk = []
stk.append([start_index, last_index])
global_start_index = start_index
indices = np.ones(last_index - start_index + 1, dtype=bool)
while stk:
start_index, last_index = stk.pop()
dmax = 0.0
index = start_index
for i in xrange(index + 1, last_index):
if indices[i - global_start_index]:
d = dist(M[i], M[start_index], M[last_index])
if d > dmax:
index = i
dmax = d
if dmax > epsilon:
stk.append([start_index, index])
stk.append([index, last_index])
else:
for i in xrange(start_index + 1, last_index):
indices[i - global_start_index] = False
return indices
[docs]def rdp_iter(M, epsilon, dist=pldist, return_mask=False):
"""
Simplifies a given array of points.
Iterative version.
:param M: an array
:type M: numpy array
:param epsilon: epsilon in the rdp algorithm
:type epsilon: float
:param dist: distance function
:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
:param return_mask: return the mask of points to keep instead
:type return_mask: bool
"""
mask = _rdp_iter(M, 0, len(M) - 1, epsilon, dist)
if return_mask:
return mask
return M[mask]
[docs]def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):
"""
Simplifies a given array of points using the Ramer-Douglas-Peucker
algorithm.
Example:
>>> from rdp import rdp
>>> rdp([[1, 1], [2, 2], [3, 3], [4, 4]])
[[1, 1], [4, 4]]
This is a convenience wrapper around both :func:`rdp.rdp_iter`
and :func:`rdp.rdp_rec` that detects if the input is a numpy array
in order to adapt the output accordingly. This means that
when it is called using a Python list as argument, a Python
list is returned, and in case of an invocation using a numpy
array, a NumPy array is returned.
The parameter ``return_mask=True`` can be used in conjunction
with ``algo="iter"`` to return only the mask of points to keep. Example:
>>> from rdp import rdp
>>> import numpy as np
>>> arr = np.array([1, 1, 2, 2, 3, 3, 4, 4]).reshape(4, 2)
>>> arr
array([[1, 1],
[2, 2],
[3, 3],
[4, 4]])
>>> mask = rdp(arr, algo="iter", return_mask=True)
>>> mask
array([ True, False, False, True], dtype=bool)
>>> arr[mask]
array([[1, 1],
[4, 4]])
:param M: a series of points
:type M: numpy array with shape ``(n,d)`` where ``n`` is the number of points and ``d`` their dimension
:param epsilon: epsilon in the rdp algorithm
:type epsilon: float
:param dist: distance function
:type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
:param algo: either ``iter`` for an iterative algorithm or ``rec`` for a recursive algorithm
:type algo: string
:param return_mask: return mask instead of simplified array
:type return_mask: bool
"""
if algo == "iter":
algo = partial(rdp_iter, return_mask=return_mask)
elif algo == "rec":
if return_mask:
raise NotImplementedError("return_mask=True not supported with algo=\"rec\"")
algo = rdp_rec
if "numpy" in str(type(M)):
return algo(M, epsilon, dist)
return algo(np.array(M), epsilon, dist).tolist()