Алгоритм рамера дугласа пекера python

Python implementation of the Ramer-Douglas-Peucker Algorithm

I recently implemented the RDP polygon approximation algorithm in Python and I’m skeptical of whether or not I implemented it correctly of with the greatest efficiency. The algorithm runs in around 0.0003 seconds for a polygon with 30 points on my computer (an Intel Core i5 with 3.8 GHz of RAM), so I’m worried about how it may run on a slower computer. Also, there seems to be a cap as to the number of points that can be removed by the algorithm, or at least there’s a cap in my implementation. No matter how high I set the tolerance, the approximation always caps at about \$\frac<2n>\$ points where \$N\$ is the number of points in the input polygon. Could I be doing something wrong?

NegInf = float('-inf') def distance(v1, v2): """ Calculate the distance between two points. PARAMETERS ========== v1, v2 >> The first and second vertices respectively. """ dx = v2[0] - v1[0] dy = v2[1] - v1[1] return math.sqrt(dx*dx + dy*dy) def perpendicular_distance(point, line_start, line_end): """ Calculate the perpendicular distance from a point to a line. PARAMETERS ========== point >> The point of which to calculate the distance from the line (must be an (x, y) tuple). line_start, line_end >> Start and end points defining the line respectively (must each be an (x, y) tuple). """ x1, y1 = line_start x2, y2 = line_end vx, vy = point if x1 == x2: return abs(x1 - vx) m = (y2 - y1)/(x2 - x1) b = y1 - m*x1 return abs(m * vx - vy + b)/math.sqrt(m*m + 1) def _rdp_approx(points, tolerance, depth): """ Internal Function: Recursively perform the RDP algorithm. """ if not points: # In case the furthest point index discovered is equal to the length of the # list of points, leading to points[furthest:] sending in an empty list. return [] elif len(points) max_dist: max_dist = dist furthest = i # Recursively calculate the RDP approximation of the two polygonal chains formed by # slicing at the index of the furthest discovered point. prev_points = _rdp_approx(points[:furthest+1], tolerance, depth+1) next_points = _rdp_approx(points[furthest:], tolerance, depth+1) new_points = [] for point in prev_points + next_points: # Filter out the duplicate points whilst maintaining order. # TODO:: There's probably some fancy slicing trick I just haven't figured out # that can be applied to prev_points and next_points so that we don't have to # do this, but it's not a huge bottleneck so we needn't worry about it now. if point not in new_points: new_points.append(point) return new_points def rdp_polygon_approximate(coordinates, tolerance): """ Use the Ramer-Douglas-Peucker algorithm to approximate the points on a polygon. The RDP algorithm recursively cuts away parts of a polygon that stray from the average of the edges. It is a great algorithm for maintaining the overall form of the input polygon, however one should be careful when using this for larger polygons as the algorithm has an average complexity of T(n) = 2T(n/2) + O(n) and a worst case complexity of O(n^2). PARAMETERS ========== coordinates >> The coordinates of the polygon to approximate. tolerance >> The amount of tolerance the algorithm will use. The tolerance determines the minimum distance a point has to sway from the average before it gets deleted from the polygon. Thus, setting the tolerance to be higher should delete more points on the final polygon. That said, due to how the algorithm works there is a limit to the number of vertices that can be removed on a polygon. Setting the tolerance to float('inf') or sys.maxsize will not approximate the polygon to a single point. Usually the minimum points an approximated polygon can have if the original polygon had N points is between 2N/3 and N/3. FURTHER READING =============== For further reading on the Ramer-Douglas-Peucker algorithm, see http://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm """ return _rdp_approx(coordinates, tolerance, 1) if __name__ == '__main__': poly = [(3, 0), (4, 2), (5, 2), (5.5, 3), (5, 4), (4, 5), (5, 6), (7, 5), (7, 3), (8, 2.5), (8, 4), (9, 5), (8, 7), (7, 8), (6, 7), (4, 7.75), (3.5, 7.5), (3, 8), (3, 8.5), (2.5, 9), (1, 9), (0, 8), (2, 7), (1, 7), (0, 6), (1, 4), (2, 5), (2, 2), (3, 3), (2, 1)] print(rdp_polygon_approximate(poly, 3)) print(rdp_polygon_approximate(poly, float('inf'))) 

Источник

Python : Ramer-Douglas-Peucker (RDP) algorithm with number of points instead of epsilon

I would like to modify this following python script for RDP algorithm with the purpose of not using epsilon but to choose the number of points I want to keep at the final :

class DPAlgorithm(): def distance(self, a, b): return sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) def point_line_distance(self, point, start, end): if (start == end): return self.distance(point, start) else: n = abs( (end[0] - start[0]) * (start[1] - point[1]) - (start[0] - point[0]) * (end[1] - start[1]) ) d = sqrt( (end[0] - start[0]) ** 2 + (end[1] - start[1]) ** 2 ) return n / d def rdp(self, points, epsilon): """ Reduces a series of points to a simplified version that loses detail, but maintains the general shape of the series. """ dmax = 0.0 index = 0 i=1 for i in range(1, len(points) - 1): d = self.point_line_distance(points[i], points[0], points[-1]) if d > dmax : index = i dmax = d if dmax >= epsilon : results = self.rdp(points[:index+1], epsilon)[:-1] + self.rdp(points[index:], epsilon) else: results = [points[0], points[-1]] return results 

I found a Java script in that spirit : https://gist.github.com/msbarry/9152218 Does anyone know a version for Python 3.X ? Thanks Momow

1 Answer 1

Ported JS code from the link above to Python [2.7]:

# -*- coding: utf-8 -*- import math import time def timenow(): return int(time.time() * 1000) def sqr(x): return x*x def distSquared(p1, p2): return sqr(p1[0] - p2[0]) + sqr(p1[1] - p2[1]) class Line(object): def __init__(self, p1, p2): self.p1 = p1 self.p2 = p2 self.lengthSquared = distSquared(self.p1, self.p2) def getRatio(self, point): segmentLength = self.lengthSquared if segmentLength == 0: return distSquared(point, p1); return ((point[0] - self.p1[0]) * (self.p2[0] - self.p1[0]) + \ (point[1] - self.p1[1]) * (self.p2[1] - self.p1[1])) / segmentLength def distanceToSquared(self, point): t = self.getRatio(point) if t < 0: return distSquared(point, self.p1) if t >1: return distSquared(point, self.p2) return distSquared(point, [ self.p1[0] + t * (self.p2[0] - self.p1[0]), self.p1[1] + t * (self.p2[1] - self.p1[1]) ]) def distanceTo(self, point): return math.sqrt(self.distanceToSquared(point)) def simplifyDouglasPeucker(points, pointsToKeep): weights = [] length = len(points) def douglasPeucker(start, end): if (end > start + 1): line = Line(points[start], points[end]) maxDist = -1 maxDistIndex = 0 for i in range(start + 1, end): dist = line.distanceToSquared(points[i]) if dist > maxDist: maxDist = dist maxDistIndex = i weights.insert(maxDistIndex, maxDist) douglasPeucker(start, maxDistIndex) douglasPeucker(maxDistIndex, end) douglasPeucker(0, length - 1) weights.insert(0, float("inf")) weights.append(float("inf")) weightsDescending = weights weightsDescending = sorted(weightsDescending, reverse=True) maxTolerance = weightsDescending[pointsToKeep - 1] result = [ point for i, point in enumerate(points) if weights[i] >= maxTolerance ] return result 

Источник

Ramer-Douglas-Peucker Algorithm¶

The Ramer–Douglas–Peucker algorithm (RDP) is an algorithm for reducing the number of points in a curve that is approximated by a series of points.

An interactive version of this algorithm can be found in this blog post.

This implementation works on 2D and 3D data.

Installation¶

The rdp package is available via pip:

The code of this package is hosted at GitHub.

Usage¶

Simplifies a given array of points using the Ramer-Douglas-Peucker algorithm.

>>> from rdp import rdp >>> rdp([[1, 1], [2, 2], [3, 3], [4, 4]]) [[1, 1], [4, 4]] 

This is a convenience wrapper around both rdp.rdp_iter() and 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]]) 

• M (numpy array with shape (n,d) where n is the number of points and d their dimension) – a series of points
• epsilon (float) – epsilon in the rdp algorithm
• dist (function with signature f(point, start, end) – see rdp.pldist() ) – distance function
• algo (string) – either iter for an iterative algorithm or rec for a recursive algorithm
• return_mask (bool) – return mask instead of simplified array

Simplifies a given array of points.

Simplifies a given array of points.

Calculates the distance from point to the line given by the points start and end .

© Copyright 2016, Fabian Hirschmann. Revision 715a436f .

Versions latest stable Downloads pdf htmlzip epub On Read the Docs Project Home Builds Free document hosting provided by Read the Docs.

Источник

Реализации алгоритмов/Алгоритм Рамера — Дугласа — Пекера

Алгоритм Рамера — Дугласа — Пекера — это алгоритм, позволяющий уменьшить число точек кривой, аппроксимированной большей серией точек. Алгоритм был независимо открыт Урсом Рамером в 1972 и Давидом Дугласом и Томасом Пекером в 1973. Также алгоритм известен под следующими именами: алгоритм Дугласа — Пекера, алгоритм итеративной ближайшей точки и алгоритм разбиения и слияния.

Суть алгоритма состоит в том, чтобы по данной ломаной, аппроксимирующей кривую, построить ломаную с меньшим числом точек. Алгоритм определяет расхождение, которое вычисляется по максимальному расстоянию между исходной и упрощённой кривыми. Упрощенная кривая состоит из подмножества точек, которые определяются из исходной кривой.

Псевдокод [ править ]

function DouglasPeucker(PointList[], epsilon) // Находим точку с максимальным расстоянием от прямой между первой и последней точками набора dmax = 0 index = 0 for i = 2 to (length(PointList) - 1) d = PerpendicularDistance(PointList[i], Line(PointList[1], PointList[end])) if d > dmax index = i dmax = d end end // Если максимальная дистанция больше, чем epsilon, то рекурсивно вызываем её на участках if dmax >= epsilon //Recursive call recResults1[] = DouglasPeucker(PointList[1. index], epsilon) recResults2[] = DouglasPeucker(PointList[index. end], epsilon) // Строим итоговый набор точек ResultList[] = else ResultList[] = end // Возвращаем результат return ResultList[] end

Python [ править ]

import math def simplify_points (pts, tolerance): anchor = 0 floater = len(pts) - 1 stack = [] keep = set() stack.append((anchor, floater)) while stack: anchor, floater = stack.pop() # инициализация отрезка if pts[floater] != pts[anchor]: anchorX = float(pts[floater][0] - pts[anchor][0]) anchorY = float(pts[floater][1] - pts[anchor][1]) seg_len = math.sqrt(anchorX ** 2 + anchorY ** 2) # get the unit vector anchorX /= seg_len anchorY /= seg_len else: anchorX = anchorY = seg_len = 0.0 # внутренний цикл: max_dist = 0.0 farthest = anchor + 1 for i in range(anchor + 1, floater): dist_to_seg = 0.0 # compare to anchor vecX = float(pts[i][0] - pts[anchor][0]) vecY = float(pts[i][1] - pts[anchor][1]) seg_len = math.sqrt( vecX ** 2 + vecY ** 2 ) # dot product: proj = vecX * anchorX + vecY * anchorY if proj  0.0: dist_to_seg = seg_len else: # compare to floater vecX = float(pts[i][0] - pts[floater][0]) vecY = float(pts[i][1] - pts[floater][1]) seg_len = math.sqrt( vecX ** 2 + vecY ** 2 ) # dot product: proj = vecX * (-anchorX) + vecY * (-anchorY) if proj  0.0: dist_to_seg = seg_len else: # расстояние от до прямой по теореме Пифагора: dist_to_seg = math.sqrt(abs(seg_len ** 2 - proj ** 2)) if max_dist  dist_to_seg: max_dist = dist_to_seg farthest = i if max_dist  tolerance: # использование отрезка keep.add(anchor) keep.add(floater) else: stack.append((anchor, farthest)) stack.append((farthest, floater)) keep = list(keep) keep.sort() return [pts[i] for i in keep] 

>>> line = [(0,0),(1,0),(2,0),(2,1),(2,2),(1,2),(0,2),(0,1),(0,0)] >>> simplify_points(line, 1.0) [(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)] >>> line = [(0,0),(0.5,0.5),(1,0),(1.25,-0.25),(1.5,.5)] >>> simplify_points(line, 0.25) [(0, 0), (0.5, 0.5), (1.25, -0.25), (1.5, 0.5)] 

Примечания [ править ]

1. ↑Pure-Python Douglas-Peucker line simplification/generalization This code was written by Schuyler Erle and is made available in the public domain. The code was ported from a freely-licensed example at http://www.3dsoftware.com/Cartography/Programming/PolyLineReduction/ The original page is no longer available, but is mirrored at http://www.mappinghacks.com/code/PolyLineReduction/

Источник