Ramer-Douglas-Peucker line simplification
You are encouraged to solve this task according to the task description, using any language you may know.
The Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.
- Task
Using the Ramer–Douglas–Peucker algorithm, 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)
The error threshold to be used is: 1.0.
Display the remaining points here.
- Reference
-
- the Wikipedia article: Ramer-Douglas-Peucker algorithm.
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
C#
<lang csharp>using System; using System.Collections.Generic; using System.Linq;
namespace LineSimplification {
using Point = Tuple<double, double>;
class Program {
static double PerpendicularDistance(Point pt, Point lineStart, Point lineEnd) {
double dx = lineEnd.Item1 - lineStart.Item1;
double dy = lineEnd.Item2 - lineStart.Item2;
// Normalize
double mag = Math.Sqrt(dx * dx + dy * dy);
if (mag > 0.0) {
dx /= mag;
dy /= mag;
}
double pvx = pt.Item1 - lineStart.Item1;
double pvy = pt.Item2 - lineStart.Item2;
// Get dot product (project pv onto normalized direction)
double pvdot = dx * pvx + dy * pvy;
// Scale line direction vector and subtract it from pv
double ax = pvx - pvdot * dx;
double ay = pvy - pvdot * dy;
return Math.Sqrt(ax * ax + ay * ay);
}
static void RamerDouglasPeucker(List<Point> pointList, double epsilon, List<Point> output) {
if (pointList.Count < 2) {
throw new ArgumentOutOfRangeException("Not enough points to simplify");
}
// Find the point with the maximum distance from line between the start and end
double dmax = 0.0;
int index = 0;
int end = pointList.Count - 1;
for (int 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) {
List<Point> recResults1 = new List<Point>();
List<Point> recResults2 = new List<Point>();
List<Point> firstLine = pointList.Take(index + 1).ToList();
List<Point> lastLine = pointList.Skip(index).ToList();
RamerDouglasPeucker(firstLine, epsilon, recResults1);
RamerDouglasPeucker(lastLine, epsilon, recResults2);
// build the result list
output.AddRange(recResults1.Take(recResults1.Count - 1));
output.AddRange(recResults2);
if (output.Count < 2) throw new Exception("Problem assembling output");
}
else {
// Just return start and end points
output.Clear();
output.Add(pointList[0]);
output.Add(pointList[pointList.Count - 1]);
}
}
static void Main(string[] args) {
List<Point> pointList = new List<Point>() {
new Point(0.0,0.0),
new Point(1.0,0.1),
new Point(2.0,-0.1),
new Point(3.0,5.0),
new Point(4.0,6.0),
new Point(5.0,7.0),
new Point(6.0,8.1),
new Point(7.0,9.0),
new Point(8.0,9.0),
new Point(9.0,9.0),
};
List<Point> pointListOut = new List<Point>();
RamerDouglasPeucker(pointList, 1.0, pointListOut);
Console.WriteLine("Points remaining after simplification:");
pointListOut.ForEach(p => Console.WriteLine(p));
}
}
}</lang>
- Output:
Points remaining after simplification: (0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
D
<lang D>import std.algorithm; import std.exception : enforce; import std.math; import std.stdio;
void main() {
creal[] pointList = [
0.0 + 0.0i,
1.0 + 0.1i,
2.0 + -0.1i,
3.0 + 5.0i,
4.0 + 6.0i,
5.0 + 7.0i,
6.0 + 8.1i,
7.0 + 9.0i,
8.0 + 9.0i,
9.0 + 9.0i
];
creal[] pointListOut;
ramerDouglasPeucker(pointList, 1.0, pointListOut);
writeln("result");
for (size_t i=0; i< pointListOut.length; i++) {
writeln(pointListOut[i].re, ",", pointListOut[i].im);
}
}
real perpendicularDistance(const creal pt, const creal lineStart, const creal lineEnd) {
creal d = lineEnd - lineStart;
//Normalise
real mag = hypot(d.re, d.im);
if (mag > 0.0) {
d /= mag;
}
creal pv = pt - lineStart;
//Get dot product (project pv onto normalized direction) real pvdot = d.re * pv.re + d.im * pv.im;
//Scale line direction vector creal ds = pvdot * d;
//Subtract this from pv creal a = pv - ds;
return hypot(a.re, a.im);
}
void ramerDouglasPeucker(const creal[] pointList, real epsilon, ref creal[] output) {
enforce(pointList.length >= 2, "Not enough points to simplify");
// Find the point with the maximum distance from line between start and end
real dmax = 0.0;
size_t index = 0;
size_t end = pointList.length-1;
for (size_t i=1; i<end; i++) {
real 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
creal[] firstLine = pointList[0..index+1].dup;
creal[] lastLine = pointList[index+1..$].dup;
creal[] recResults1;
ramerDouglasPeucker(firstLine, epsilon, recResults1);
creal[] recResults2;
ramerDouglasPeucker(lastLine, epsilon, recResults2);
// Build the result list
output = recResults1 ~ recResults2;
enforce(output.length>=2, "Problem assembling output");
} else {
//Just return start and end points
output.length = 0;
output ~= pointList[0];
output ~= pointList[end];
}
}</lang>
- Output:
result 0,0 2,-0.1 3,5 7,9 9,9
Go
<lang go>package main
import (
"fmt" "math"
)
type point struct{ x, y float64 }
func RDP(l []point, ε float64) []point {
x := 0
dMax := -1.
last := len(l) - 1
p1 := l[0]
p2 := l[last]
x21 := p2.x - p1.x
y21 := p2.y - p1.y
for i, p := range l[1:last] {
if d := math.Abs(y21*p.x - x21*p.y + p2.x*p1.y - p2.y*p1.x); d > dMax {
x = i + 1
dMax = d
}
}
if dMax > ε {
return append(RDP(l[:x+1], ε), RDP(l[x:], ε)[1:]...)
}
return []point{l[0], l[len(l)-1]}
}
func main() {
fmt.Println(RDP([]point{{0, 0}, {1, 0.1}, {2, -0.1},
{3, 5}, {4, 6}, {5, 7}, {6, 8.1}, {7, 9}, {8, 9}, {9, 9}}, 1))
}</lang>
- Output:
[{0 0} {2 -0.1} {3 5} {7 9} {9 9}]
J
Solution: <lang j>mp=: +/ .* NB. matrix product norm=: +/&.:*: NB. vector norm normalize=: (% norm)^:(0 < norm)
dxy=. normalize@({: - {.) pv=. -"1 {. NB.*perpDist v Calculate perpendicular distance of points from a line perpDist=: norm"1@(pv ([ -"1 mp"1~ */ ]) dxy) f.
rdp=: verb define
1 rdp y : points=. ,:^:(2 > #@$) y epsilon=. x if. 2 > # points do. points return. end.
NB. point with the maximum distance from line between start and end
'imax dmax'=. ((i. , ]) >./) perpDist points
if. dmax > epsilon do.
epsilon ((}:@rdp (1+imax)&{.) , (rdp imax&}.)) points
else.
({. ,: {:) points
end.
)</lang> Example Usage: <lang j> Points=: 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 rdp Points
0 0 2 _0.1 3 5 7 9 9 9</lang>
Java
<lang Java>import javafx.util.Pair;
import java.util.ArrayList; import java.util.List;
public class LineSimplification {
private static class Point extends Pair<Double, Double> {
Point(Double key, Double value) {
super(key, value);
}
@Override
public String toString() {
return String.format("(%f, %f)", getKey(), getValue());
}
}
private static double perpendicularDistance(Point pt, Point lineStart, Point lineEnd) {
double dx = lineEnd.getKey() - lineStart.getKey();
double dy = lineEnd.getValue() - lineStart.getValue();
// Normalize
double mag = Math.hypot(dx, dy);
if (mag > 0.0) {
dx /= mag;
dy /= mag;
}
double pvx = pt.getKey() - lineStart.getKey();
double pvy = pt.getValue() - lineStart.getValue();
// Get dot product (project pv onto normalized direction)
double pvdot = dx * pvx + dy * pvy;
// Scale line direction vector and subtract it from pv
double ax = pvx - pvdot * dx;
double ay = pvy - pvdot * dy;
return Math.hypot(ax, ay); }
private static void ramerDouglasPeucker(List<Point> pointList, double epsilon, List<Point> out) {
if (pointList.size() < 2) throw new IllegalArgumentException("Not enough points to simplify");
// Find the point with the maximum distance from line between the start and end
double dmax = 0.0;
int index = 0;
int end = pointList.size() - 1;
for (int i = 1; i < end; ++i) {
double d = perpendicularDistance(pointList.get(i), pointList.get(0), pointList.get(end));
if (d > dmax) {
index = i;
dmax = d;
}
}
// If max distance is greater than epsilon, recursively simplify
if (dmax > epsilon) {
List<Point> recResults1 = new ArrayList<>();
List<Point> recResults2 = new ArrayList<>();
List<Point> firstLine = pointList.subList(0, index + 1);
List<Point> lastLine = pointList.subList(index, pointList.size());
ramerDouglasPeucker(firstLine, epsilon, recResults1);
ramerDouglasPeucker(lastLine, epsilon, recResults2);
// build the result list
out.addAll(recResults1.subList(0, recResults1.size() - 1));
out.addAll(recResults2);
if (out.size() < 2) throw new RuntimeException("Problem assembling output");
} else {
// Just return start and end points
out.clear();
out.add(pointList.get(0));
out.add(pointList.get(pointList.size() - 1));
}
}
public static void main(String[] args) {
List<Point> pointList = List.of(
new Point(0.0, 0.0),
new Point(1.0, 0.1),
new Point(2.0, -0.1),
new Point(3.0, 5.0),
new Point(4.0, 6.0),
new Point(5.0, 7.0),
new Point(6.0, 8.1),
new Point(7.0, 9.0),
new Point(8.0, 9.0),
new Point(9.0, 9.0)
);
List<Point> pointListOut = new ArrayList<>();
ramerDouglasPeucker(pointList, 1.0, pointListOut);
System.out.println("Points remaining after simplification:");
pointListOut.forEach(System.out::println);
}
}</lang>
- Output:
Points remaining after simplification: (0.000000, 0.000000) (2.000000, -0.100000) (3.000000, 5.000000) (7.000000, 9.000000) (9.000000, 9.000000)
Julia
<lang julia>const Point = Vector{Float64}
function perpdist(pt::Point, lnstart::Point, lnend::Point)
d = normalize!(lnend .- lnstart)
pv = pt .- lnstart # Get dot product (project pv onto normalized direction) pvdot = dot(d, pv) # Scale line direction vector ds = pvdot .* d # Subtract this from pv return norm(pv .- ds)
end
function rdp(plist::Vector{Point}, ϵ::Float64 = 1.0)
if length(plist) < 2
throw(ArgumentError("not enough points to simplify"))
end
# Find the point with the maximum distance from line between start and end distances = collect(perpdist(pt, plist[1], plist[end]) for pt in plist) dmax, imax = findmax(distances)
# If max distance is greater than epsilon, recursively simplify
if dmax > ϵ
fstline = plist[1:imax]
lstline = plist[imax:end]
recrst1 = rdp(fstline, ϵ)
recrst2 = rdp(lstline, ϵ)
out = vcat(recrst1, recrst2)
else
out = [plist[1], plist[end]]
end
return out
end
plist = Point[[0.0, 0.0], [1.0, 0.1], [2.0, -0.1], [3.0, 5.0], [4.0, 6.0], [5.0, 7.0], [6.0, 8.1], [7.0, 9.0], [8.0, 9.0], [9.0, 9.0]] @show plist @show rdp(plist)</lang>
- Output:
plist = Array{Float64,1}[[0.0, 0.0], [1.0, 0.1], [2.0, -0.1], [3.0, 5.0], [4.0, 6.0], [5.0, 7.0], [6.0, 8.1], [7.0, 9.0], [8.0, 9.0], [9.0, 9.0]]
rdp(plist) = Array{Float64,1}[[0.0, 0.0], [2.0, -0.1], [2.0, -0.1], [3.0, 5.0], [3.0, 5.0], [7.0, 9.0], [7.0, 9.0], [9.0, 9.0]]
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)
Nim
<lang nim>import math
type
Point = tuple[x, y: float64]
proc pointLineDistance(pt, lineStart, lineEnd: Point): float64 =
var n, d, dx, dy: float64 dx = lineEnd.x - lineStart.x dy = lineEnd.y - lineStart.y n = abs(dx * (lineStart.y - pt.y) - (lineStart.x - pt.x) * dy) d = sqrt(dx * dx + dy * dy) n / d
proc rdp(points: seq[Point], startIndex, lastIndex: int, ε: float64 = 1.0): seq[Point] =
var dmax = 0.0 var index = startIndex
for i in index+1..<lastIndex:
var d = pointLineDistance(points[i], points[startIndex], points[lastIndex])
if d > dmax:
index = i
dmax = d
if dmax > ε:
var res1 = rdp(points, startIndex, index, ε)
var res2 = rdp(points, index, lastIndex, ε)
var finalRes: seq[Point] = @[] finalRes.add(res1[0..^2]) finalRes.add(res2[0..^1]) result = finalRes else: result = @[points[startIndex], points[lastIndex]]
var line: seq[Point] = @[(0.0, 0.0), (1.0, 0.1), (2.0, -0.1), (3.0, 5.0), (4.0, 6.0),
(5.0, 7.0), (6.0, 8.1), (7.0, 9.0), (8.0, 9.0), (9.0, 9.0)]
echo rdp(line, line.low, line.high)</lang>
- Output:
@[(x: 0.0, y: 0.0), (x: 2.0, y: -0.1), (x: 3.0, y: 5.0), (x: 7.0, y: 9.0), (x: 9.0, y: 9.0)]
Perl 6
<lang perl6>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)]
);</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))
REXX
The computation for the perpendicular distance was taken from the GO example. <lang rexx>/*REXX program uses the Ramer─Douglas─Peucker (RDP) line simplification algorithm for*/ /*───────────────────────────── reducing the number of points used to define its shape. */ parse arg epsilon pts /*obtain optional arguments from the CL*/ if epsilon= | epsilon="," then epsilon= 1 /*Not specified? Then use the default.*/ if pts= then pts= '(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)' pts= space(pts) /*elide all superfluous blanks. */ say ' error threshold: ' epsilon /*echo the error threshold to the term.*/ say ' points specified: ' pts /* " " shape points " " " */ $= RDP(pts) /*invoke Ramer─Douglas─Peucker function*/ say 'points simplified: ' rez($) /*display points with () ───► terminal.*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ bld: parse arg _; #= words(_); dMax=-#; idx=1; do j=1 for #; @.j= word(_, j); end; return px: parse arg _; return word( translate(_, , ','), 1) /*obtain the X coörd.*/ py: parse arg _; return word( translate(_, , ','), 2) /* " " Y " */ reb: parse arg a,b,,_; do k=a to b; _=_ @.k; end; return strip(_) rez: parse arg z,_; do k=1 for words(z); _= _ '('word(z, k)") "; end; return strip(_) /*──────────────────────────────────────────────────────────────────────────────────────*/ RDP: procedure expose epsilon; parse arg PT; call bld space( translate(PT, , ')(][}{') )
L= px(@.#)-px(@.1)
H= py(@.#)-py(@.1) /* [↓] find point IDX with max distance*/
do i=2 to #-1
d= abs(H*px(@.i) - L*py(@.i) + px(@.#)*py(@.1) - py(@.#)*px(@.1) )
if d>dMax then do; idx= i; dMax= d
end
end /*i*/ /* [↑] D is the perpendicular distance*/
if dMax>epsilon then do; r= RDP( reb(1, idx) )
return subword(r, 1, words(r) - 1) RDP( reb(idx, #) )
end
return @.1 @.#</lang>
- output when using the default inputs:
error threshold: 1 points specified: (0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9) points simplified: (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))
References