Ramer-Douglas-Peucker line simplification

From Rosetta Code
Ramer-Douglas-Peucker line simplification is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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++[edit]

#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;
}
Output:
result
0,0
2,-0.1
3,5
7,9
9,9

Kotlin[edit]

Translation of: C++
// 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)
}
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[edit]

Works with: Rakudo version 2017.05
Translation of: C++
sub norm (*@list) { @list»².sum.sqrt }
 
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) «/=» norm $Δx, $Δy;
norm ($Δpx, $Δpy) «-» (x, $Δy) «*» (x*$Δpx + $Δy*$Δpy);
}
 
sub Ramer-Douglas-Peucker(@points where * > 1,= 1) {
return @points if @points == 2;
my @d = (^@points).map: { perpendicular-distance |@points[0, *-1, $_] };
my ($index, $dmax) = @d.first: @d.max, :kv;
return flat
Ramer-Douglas-Peucker( @points[0..$index], ε )[^(*-1)],
Ramer-Douglas-Peucker( @points[$index..*], ε )
if $dmax > ε;
@points[0, *-1];
}
 
# TESTING
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)]
);
Output:
((0 0) (2 -0.1) (3 5) (7 9) (9 9))

Python[edit]

An approach using the shapely library:

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))
Output:
LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)

Racket[edit]

#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)))
 
;; https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
(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))
Output:
'((0 0) (2 -0.1) (3 5) (7 9) (9 9))

Sidef[edit]

Translation of: Perl 6
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]]
)
Output:
[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]]

zkl[edit]

Translation of: Perl 6
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]);
}
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();
Output:
L(L(0,0),L(2,-0.1),L(3,5),L(7,9),L(9,9))

References[edit]

  1. [1]