Ramer-Douglas-Peucker line simplification: Difference between revisions
Thundergnat (talk | contribs) m (→{{header|Perl 6}}: minor simplification) |
|||
Line 235: | Line 235: | ||
if $dmax > ε { |
if $dmax > ε { |
||
return Ramer-Douglas-Peucker( @points[0..$index] , ε )[ |
return Ramer-Douglas-Peucker( @points[0..$index] , ε )[^(*-1)], |
||
Ramer-Douglas-Peucker( @points[$index..*-1], ε ) |
Ramer-Douglas-Peucker( @points[$index..*-1], ε ) |
||
} else { |
|||
⚫ | |||
} |
} |
||
⚫ | |||
} |
} |
||
# TESTING |
# TESTING |
||
say flat Ramer-Douglas-Peucker( |
say flat Ramer-Douglas-Peucker( |
Revision as of 17:47, 16 June 2017
Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.[1]
- Task
Simplify the 2D line defined by the points (0,0),(1,0.1)(2,-0.1)(3,5)(4,6)(5,7)(6,8.1)(7,9)(8,9)(9,9) using the Ramer–Douglas–Peucker algorithm. The error threshold to use is 1.0. Display the remaining points.
C++
<lang cpp>#include <iostream>
- include <cmath>
- include <utility>
- include <vector>
- include <stdexcept>
using namespace std;
typedef std::pair<double, double> Point;
double PerpendicularDistance(const Point &pt, const Point &lineStart, const Point &lineEnd) { double dx = lineEnd.first - lineStart.first; double dy = lineEnd.second - lineStart.second;
//Normalise double mag = pow(pow(dx,2.0)+pow(dy,2.0),0.5); if(mag > 0.0) { dx /= mag; dy /= mag; }
double pvx = pt.first - lineStart.first; double pvy = pt.second - lineStart.second;
//Get dot product (project pv onto normalized direction) double pvdot = dx * pvx + dy * pvy;
//Scale line direction vector double dsx = pvdot * dx; double dsy = pvdot * dy;
//Subtract this from pv double ax = pvx - dsx; double ay = pvy - dsy;
return pow(pow(ax,2.0)+pow(ay,2.0),0.5); }
void RamerDouglasPeucker(const vector<Point> &pointList, double epsilon, vector<Point> &out) { if(pointList.size()<2) throw invalid_argument("Not enough points to simplify");
// Find the point with the maximum distance from line between start and end double dmax = 0.0; size_t index = 0; size_t end = pointList.size()-1; for(size_t i = 1; i < end; i++) { double d = PerpendicularDistance(pointList[i], pointList[0], pointList[end]); if (d > dmax) { index = i; dmax = d; } }
// If max distance is greater than epsilon, recursively simplify if(dmax > epsilon) { // Recursive call vector<Point> recResults1; vector<Point> recResults2; vector<Point> firstLine(pointList.begin(), pointList.begin()+index+1); vector<Point> lastLine(pointList.begin()+index, pointList.end()); RamerDouglasPeucker(firstLine, epsilon, recResults1); RamerDouglasPeucker(lastLine, epsilon, recResults2);
// Build the result list out.assign(recResults1.begin(), recResults1.end()-1); out.insert(out.end(), recResults2.begin(), recResults2.end()); if(out.size()<2) throw runtime_error("Problem assembling output"); } else { //Just return start and end points out.clear(); out.push_back(pointList[0]); out.push_back(pointList[end]); } }
int main() { vector<Point> pointList; vector<Point> pointListOut;
pointList.push_back(Point(0.0, 0.0)); pointList.push_back(Point(1.0, 0.1)); pointList.push_back(Point(2.0, -0.1)); pointList.push_back(Point(3.0, 5.0)); pointList.push_back(Point(4.0, 6.0)); pointList.push_back(Point(5.0, 7.0)); pointList.push_back(Point(6.0, 8.1)); pointList.push_back(Point(7.0, 9.0)); pointList.push_back(Point(8.0, 9.0)); pointList.push_back(Point(9.0, 9.0));
RamerDouglasPeucker(pointList, 1.0, pointListOut);
cout << "result" << endl; for(size_t i=0;i< pointListOut.size();i++) { cout << pointListOut[i].first << "," << pointListOut[i].second << endl; }
return 0; }</lang>
- Output:
result 0,0 2,-0.1 3,5 7,9 9,9
Kotlin
<lang scala>// version 1.1.0
typealias Point = Pair<Double, Double>
fun perpendicularDistance(pt: Point, lineStart: Point, lineEnd: Point): Double {
var dx = lineEnd.first - lineStart.first var dy = lineEnd.second - lineStart.second
// Normalize val mag = Math.hypot(dx, dy) if (mag > 0.0) { dx /= mag; dy /= mag } val pvx = pt.first - lineStart.first val pvy = pt.second - lineStart.second
// Get dot product (project pv onto normalized direction) val pvdot = dx * pvx + dy * pvy // Scale line direction vector and substract it from pv val ax = pvx - pvdot * dx val ay = pvy - pvdot * dy return Math.hypot(ax, ay)
}
fun RamerDouglasPeucker(pointList: List<Point>, epsilon: Double, out: MutableList<Point>) {
if (pointList.size < 2) throw IllegalArgumentException("Not enough points to simplify") // Find the point with the maximum distance from line between start and end var dmax = 0.0 var index = 0 val end = pointList.size - 1 for (i in 1 until end) { val d = perpendicularDistance(pointList[i], pointList[0], pointList[end]) if (d > dmax) { index = i; dmax = d } }
// If max distance is greater than epsilon, recursively simplify if (dmax > epsilon) { val recResults1 = mutableListOf<Point>() val recResults2 = mutableListOf<Point>() val firstLine = pointList.take(index + 1) val lastLine = pointList.drop(index) RamerDouglasPeucker(firstLine, epsilon, recResults1) RamerDouglasPeucker(lastLine, epsilon, recResults2)
// build the result list out.addAll(recResults1.take(recResults1.size - 1)) out.addAll(recResults2) if (out.size < 2) throw RuntimeException("Problem assembling output") } else { // Just return start and end points out.clear() out.add(pointList.first()) out.add(pointList.last()) }
}
fun main(args: Array<String>) {
val pointList = listOf( Point(0.0, 0.0), Point(1.0, 0.1), Point(2.0, -0.1), Point(3.0, 5.0), Point(4.0, 6.0), Point(5.0, 7.0), Point(6.0, 8.1),
Point(7.0, 9.0), Point(8.0, 9.0),
Point(9.0, 9.0) ) val pointListOut = mutableListOf<Point>() RamerDouglasPeucker(pointList, 1.0, pointListOut) println("Points remaining after simplification:") for (p in pointListOut) println(p)
}</lang>
- Output:
Points remaining after simplification: (0.0, 0.0) (2.0, -0.1) (3.0, 5.0) (7.0, 9.0) (9.0, 9.0)
Perl 6
<lang perl6>sub perpendicular-distance (@start, @end where @end !eqv @start , @point) {
return 0 if @point eqv any(@start, @end); my ($Δx, $Δy ) = @end «-» @start; my ($Δpx, $Δpy) = @point «-» @start; ($Δx, $Δy) «/=» ($Δx² + $Δy²).sqrt; ((($Δpx, $Δpy) «-» ($Δx, $Δy) «*» ($Δx*$Δpx + $Δy*$Δpy)) «**» 2).sum.sqrt;
}
sub Ramer-Douglas-Peucker( @points where * > 1, \ε = 1 ) {
return @points if @points == 2; my @d = (^@points).map: { perpendicular-distance(@points[0], @points[*-1], @points[$_]) }; my ($index, $dmax) = @d.first: @d.max, :kv;
if $dmax > ε { return Ramer-Douglas-Peucker( @points[0..$index] , ε )[^(*-1)], Ramer-Douglas-Peucker( @points[$index..*-1], ε ) } @points[0, *-1];
}
- TESTING
say flat Ramer-Douglas-Peucker(
[(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]
);</lang>
- Output:
((0 0) (2 -0.1) (3 5) (7 9) (9 9))
Python
An approach using the shapely library:
<lang python>from __future__ import print_function from shapely.geometry import LineString
if __name__=="__main__": line = LineString([(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]) print (line.simplify(1.0, preserve_topology=False))</lang>
- Output:
LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)
Racket
<lang racket>#lang racket (require math/flonum)
- points are lists of x y (maybe extensible to z)
- x+y gets both parts as values
(define (x+y p) (values (first p) (second p)))
(define (⊥-distance P1 P2)
(let*-values ([(x1 y1) (x+y P1)] [(x2 y2) (x+y P2)] [(dx dy) (values (- x2 x1) (- y2 y1))] [(h) (sqrt (+ (sqr dy) (sqr dx)))]) (λ (P0) (let-values (((x0 y0) (x+y P0))) (/ (abs (+ (* dy x0) (* -1 dx y0) (* x2 y1) (* -1 y2 x1))) h)))))
(define (douglas-peucker points-in ϵ)
(let recur ((ps points-in)) ;; curried distance function which will be applicable to all points (let*-values ([(p0) (first ps)] [(pz) (last ps)] [(p-d) (⊥-distance p0 pz)] ;; Find the point with the maximum distance [(dmax index) (for/fold ((dmax 0) (index 0)) ((i (in-range 1 (sub1 (length ps))))) ; skips the first, stops before the last (define d (p-d (list-ref ps i))) (if (> d dmax) (values d i) (values dmax index)))]) ;; If max distance is greater than epsilon, recursively simplify (if (> dmax ϵ) ;; recursive call (let-values ([(l r) (split-at ps index)]) (append (drop-right (recur l) 1) (recur r))) (list p0 pz))))) ;; else we can return this simplification
(module+ main
(douglas-peucker '((0 0) (1 0.1) (2 -0.1) (3 5) (4 6) (5 7) (6 8.1) (7 9) (8 9) (9 9)) 1.0))
(module+ test
(require rackunit) (check-= ((⊥-distance '(0 0) '(0 1)) '(1 0)) 1 epsilon.0))</lang>
- Output:
'((0 0) (2 -0.1) (3 5) (7 9) (9 9))
Sidef
<lang ruby>func perpendicular_distance(Arr start, Arr end, Arr point) {
((point == start) || (point == end)) && return 0 var (Δx, Δy ) = ( end »-« start)... var (Δpx, Δpy) = (point »-« start)... var h = hypot(Δx, Δy) [\Δx, \Δy].map { *_ /= h } (([Δpx, Δpy] »-« ([Δx, Δy] »*» (Δx*Δpx + Δy*Δpy))) »**» 2).sum.sqrt
}
func Ramer_Douglas_Peucker(Arr points { .all { .len > 1 } }, ε = 1) {
points.len == 2 && return points
var d = (^points -> map { perpendicular_distance(points[0], points[-1], points[_]) })
if (d.max > ε) { var i = d.index(d.max) return [Ramer_Douglas_Peucker(points.ft(0, i), ε).ft(0, -2)..., Ramer_Douglas_Peucker(points.ft(i), ε)...] }
return [points[0,-1]]
}
say Ramer_Douglas_Peucker(
[[0,0],[1,0.1],[2,-0.1],[3,5],[4,6],[5,7],[6,8.1],[7,9],[8,9],[9,9]]
)</lang>
- Output:
[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]]
zkl
<lang zkl>fcn perpendicularDistance(start,end, point){ // all are tuples: (x,y) -->|d|
dx,dy := end .zipWith('-,start); // deltas dpx,dpy := point.zipWith('-,start); mag := (dx*dx + dy*dy).sqrt(); if(mag>0.0){ dx/=mag; dy/=mag; } p,dsx,dsy := dx*dpx + dy*dpy, p*dx, p*dy; ((dpx - dsx).pow(2) + (dpy - dsy).pow(2)).sqrt()
}
fcn RamerDouglasPeucker(points,epsilon=1.0){ // list of tuples --> same
if(points.len()==2) return(points); // but we'll do one point d:=points.pump(List, // first result/element is always zero fcn(p, s,e){ perpendicularDistance(s,e,p) }.fp1(points[0],points[-1])); index,dmax := (0.0).minMaxNs(d)[1], d[index]; // minMaxNs-->index of min & max if(dmax>epsilon){ return(RamerDouglasPeucker(points[0,index],epsilon)[0,-1].extend( RamerDouglasPeucker(points[index,*],epsilon))) } else return(points[0],points[-1]);
}</lang> <lang zkl>RamerDouglasPeucker(
T( T(0.0, 0.0), T(1.0, 0.1), T(2.0, -0.1), T(3.0, 5.0), T(4.0, 6.0), T(5.0, 7.0), T(6.0, 8.1), T(7.0, 9.0), T(8.0, 9.0), T(9.0, 9.0) ))
.println();</lang>
- Output:
L(L(0,0),L(2,-0.1),L(3,5),L(7,9),L(9,9))