Ramer-Douglas-Peucker line simplification
The Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.
You are encouraged to solve this task according to the task description, using any language you may know.
- 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.
11l
<lang 11l>F rdp(l, ε) -> [(Float, Float)]
V x = 0 V dMax = -1.0 V p1 = l[0] V p2 = l.last V p21 = p2 - p1
L(p) l[1.<(len)-1]
V d = abs(cross(p, p21) + cross(p2, p1))
I d > dMax
x = L.index + 1
dMax = d
I dMax > ε
R rdp(l[0..x], ε) [+] rdp(l[x..], ε)[1..]
R [l[0], l.last]
print(rdp([(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)], 1.0))</lang>
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
C
<lang c>#include <assert.h>
- include <math.h>
- include <stdio.h>
typedef struct point_tag {
double x; double y;
} point_t;
// Returns the distance from point p to the line between p1 and p2 double perpendicular_distance(point_t p, point_t p1, point_t p2) {
double dx = p2.x - p1.x; double dy = p2.y - p1.y; double d = sqrt(dx * dx + dy * dy); return fabs(p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x)/d;
}
// Simplify an array of points using the Ramer–Douglas–Peucker algorithm. // Returns the number of output points. size_t douglas_peucker(const point_t* points, size_t n, double epsilon,
point_t* dest, size_t destlen) {
assert(n >= 2);
assert(epsilon >= 0);
double max_dist = 0;
size_t index = 0;
for (size_t i = 1; i + 1 < n; ++i) {
double dist = perpendicular_distance(points[i], points[0], points[n - 1]);
if (dist > max_dist) {
max_dist = dist;
index = i;
}
}
if (max_dist > epsilon) {
size_t n1 = douglas_peucker(points, index + 1, epsilon, dest, destlen);
if (destlen >= n1 - 1) {
destlen -= n1 - 1;
dest += n1 - 1;
} else {
destlen = 0;
}
size_t n2 = douglas_peucker(points + index, n - index, epsilon, dest, destlen);
return n1 + n2 - 1;
}
if (destlen >= 2) {
dest[0] = points[0];
dest[1] = points[n - 1];
}
return 2;
}
void print_points(const point_t* points, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (i > 0)
printf(" ");
printf("(%g, %g)", points[i].x, points[i].y);
}
printf("\n");
}
int main() {
point_t points[] = {
{0,0}, {1,0.1}, {2,-0.1}, {3,5}, {4,6},
{5,7}, {6,8.1}, {7,9}, {8,9}, {9,9}
};
const size_t len = sizeof(points)/sizeof(points[0]);
point_t out[len];
size_t n = douglas_peucker(points, len, 1.0, out, len);
print_points(out, n);
return 0;
}</lang>
- Output:
(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)
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
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
FreeBASIC
<lang freebasic> Function DistanciaPerpendicular(tabla() As Double, i As Integer, ini As Integer, fin As Integer) As Double
Dim As Double dx, dy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay dx = tabla(fin, 1) - tabla(ini, 1) dy = tabla(fin, 2) - tabla(ini, 2) 'Normaliza mag = (dx^2 + dy^2)^0.5 If mag > 0 Then dx /= mag : dy /= mag pvx = tabla(i, 1) - tabla(ini, 1) pvy = tabla(i, 2) - tabla(ini, 2) 'Producto escalado (proyecto pv en dirección normalizada) pvdot = dx * pvx + dy * pvy 'Vector de dirección de línea de escala dsx = pvdot * dx dsy = pvdot * dy 'Reste esto de pv ax = pvx - dsx ay = pvy - dsy Return (ax^2 + ay^2)^0.5
End Function
Sub DRDP(ListaDePuntos() As Double, ini As Integer, fin As Integer, epsilon As Double)
Dim As Double dmax, d
Dim As Integer indice, i
' Encuentra el punto con la distancia máxima
If Ubound(ListaDePuntos) < 2 Then Print "No hay suficientes puntos para simplificar " : Sleep : End
For i = ini + 1 To fin
d = DistanciaPerpendicular(ListaDePuntos(), i, ini, fin)
If d > dmax Then indice = i : dmax = d
Next
' Si la distancia máxima es mayor que épsilon, simplifique de forma recursiva
If dmax > epsilon Then
ListaDePuntos(indice, 3) = True
DRDP(ListaDePuntos(), ini, indice, epsilon)
DRDP(ListaDePuntos(), indice, fin, epsilon)
End If
End Sub
Dim As Double matriz(1 To 10, 1 To 3) Data 0, 0, 1, 0.1, 2, -0.1, 3, 5, 4, 6, 5, 7, 6, 8.1, 7, 9, 8, 9, 9, 9 For i As Integer = Lbound(matriz) To Ubound(matriz)
Read matriz(i, 1), matriz(i, 2)
Next i
DRDP(matriz(), 1, 10, 1)
Print "Puntos restantes tras de simplificar:" matriz(1, 3) = True : matriz(10, 3) = True For i As Integer = Lbound(matriz) To Ubound(matriz)
If matriz(i, 3) Then Print "(";matriz(i, 1);", "; matriz(i, 2);") ";
Next i Sleep </lang>
- Output:
Puntos restantes tras de simplificar: ( 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)
JavaScript
<lang JavaScript>/**
* @typedef {{
* x: (!number),
* y: (!number)
* }}
*/
let pointType;
/**
* @param {!Array<pointType>} l
* @param {number} eps
*/
const RDP = (l, eps) => {
const last = l.length - 1; const p1 = l[0]; const p2 = l[last]; const x21 = p2.x - p1.x; const y21 = p2.y - p1.y;
const [dMax, x] = l.slice(1, last)
.map(p => Math.abs(y21 * p.x - x21 * p.y + p2.x * p1.y - p2.y * p1.x))
.reduce((p, c, i) => {
const v = Math.max(p[0], c);
return [v, v === p[0] ? p[1] : i + 1];
}, [-1, 0]);
if (dMax > eps) {
return [...RDP(l.slice(0, x + 1), eps), ...RDP(l.slice(x), eps).slice(1)];
}
return [l[0], l[last]]
};
const points = [
{x: 0, y: 0},
{x: 1, y: 0.1},
{x: 2, y: -0.1},
{x: 3, y: 5},
{x: 4, y: 6},
{x: 5, y: 7},
{x: 6, y: 8.1},
{x: 7, y: 9},
{x: 8, y: 9},
{x: 9, y: 9}];
console.log(RDP(points, 1));</lang>
- Output:
[{x: 0, y: 0},
{x: 2, y: -0.1},
{x: 3, y: 5},
{x: 7, y: 9},
{x: 9, y: 9}]
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)]
Openscad
<lang openscad>function slice(a, v) = [ for (i = v) a[i] ];
// Find the distance from the point to the line function perpendicular_distance(pt, start, end) =
let (
d = end - start,
dn = d / norm(d),
dp = pt - start,
// Dot-product of two vectors
ddot = dn * dp,
ds = dn * ddot
)
norm(dp - ds);
// Find the point with the maximum distance from line between start and end. // Returns a vector of [dmax, point_index] function max_distance(pts, cindex=1, iresult=[0,0]) =
let (
d = perpendicular_distance(pts[cindex], pts[0],pts[len(pts)-1]),
result = (d > iresult[0] ? [d, cindex] : iresult)
)
(cindex == len(pts) - 1 ?
iresult :
max_distance(pts, cindex+1, result));
// Tail recursion optimization is not possible with tree recursion. // // It's probably possible to optimize this with tail recursion by converting to // iterative approach, see // https://namekdev.net/2014/06/iterative-version-of-ramer-douglas-peucker-line-simplification-algorithm/ . // It's not possible to use iterative approach without recursion at all because // of static nature of OpenSCAD arrays. function ramer_douglas_peucker(points, epsilon=1) =
let (dmax = max_distance(points))
// Simplify the line recursively if the max distance is greater than epsilon
(dmax[0] > epsilon ?
let (
lo = ramer_douglas_peucker(
slice(points, [0:dmax[1]]), epsilon),
hi = ramer_douglas_peucker(
slice(points, [dmax[1]:len(points)-1]), epsilon)
)
concat(slice(lo, [0:len(lo)-2]), hi) :
[points[0],points[len(points)-1]]);
/* -- Demo -- */
module line(points, thickness=1, dot=2) {
for (idx = [0:len(points)-2]) {
translate(points[idx])
sphere(d=dot);
hull() {
translate(points[idx])
sphere(d=thickness);
translate(points[idx+1])
sphere(d=thickness);
}
}
translate(points[len(points)-1])
sphere(d=dot);
}
function addz(pts, z=0) = [ for (p = pts) [p.x, p.y, z] ];
module demo(pts) {
rdp = ramer_douglas_peucker(points=pts, epsilon=1);
echo(rdp);
color("gray")
line(points=addz(pts, 0), thickness=0.1, dot=0.3);
color("red")
line(points=addz(rdp, 1), thickness=0.1, dot=0.3);
}
$fn=16; points = [[0,0], [1,0.1], [2,-0.1], [3,5], [4,6], [5,7], [6,8.1], [7,9], [8,9], [9,9]]; demo(points);</lang>
- Output:
ECHO: [[0, 0], [2, -0.1], [3, 5], [7, 9], [9, 9]]
Perl
<lang perl>use strict; use warnings; use feature 'say'; use List::MoreUtils qw(firstidx minmax);
my $epsilon = 1;
sub norm {
my(@list) = @_; my $sum; $sum += $_**2 for @list; sqrt($sum)
}
sub perpendicular_distance {
our(@start,@end,@point);
local(*start,*end,*point) = (shift, shift, shift);
return 0 if $start[0]==$point[0] && $start[1]==$point[1]
or $end[0]==$point[0] && $end[1]==$point[1];
my ( $dx, $dy) = ( $end[0]-$start[0], $end[1]-$start[1]);
my ($dpx, $dpy) = ($point[0]-$start[0],$point[1]-$start[1]);
my $t = norm($dx, $dy);
$dx /= $t;
$dy /= $t;
norm($dpx - $dx*($dx*$dpx + $dy*$dpy), $dpy - $dy*($dx*$dpx + $dy*$dpy));
}
sub Ramer_Douglas_Peucker {
my(@points) = @_;
return @points if @points == 2;
my @d;
push @d, perpendicular_distance(@points[0, -1, $_]) for 0..@points-1;
my(undef,$dmax) = minmax @d;
my $index = firstidx { $_ == $dmax } @d;
if ($dmax > $epsilon) {
my @lo = Ramer_Douglas_Peucker( @points[0..$index]);
my @hi = Ramer_Douglas_Peucker( @points[$index..$#points]);
return @lo[0..@lo-2], @hi;
}
@points[0, -1];
}
say '(' . join(' ', @$_) . ') '
for 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)
Phix
<lang Phix>function rdp(sequence l, atom e)
if length(l)<2 then crash("not enough points to simplify") end if
integer idx := 0
atom dMax := -1,
{p1x,p1y} := l[1],
{p2x,p2y} := l[$],
x21 := p2x - p1x,
y21 := p2y - p1y
for i=1 to length(l) do
atom {px,py} = l[i],
d = abs(y21*px-x21*py + p2x*p1y-p2y*p1x)
if d>dMax then
idx = i
dMax = d
end if
end for
if dMax>e then
return rdp(l[1..idx], e) & rdp(l[idx..$], e)[2..$]
end if
return {l[1], l[$]}
end function
sequence points = {{0, 0}, {1, 0.1}, {2, -0.1}, {3, 5}, {4, 6},
{5, 7}, {6, 8.1}, {7, 9}, {8, 9}, {9, 9}}
?rdp(points, 1)</lang>
- Output:
{{0,0},{2,-0.1},{3,5},{7,9},{9,9}}
PHP
<lang php>function perpendicular_distance(array $pt, array $line) {
// Calculate the normalized delta x and y of the line.
$dx = $line[1][0] - $line[0][0];
$dy = $line[1][1] - $line[0][1];
$mag = sqrt($dx * $dx + $dy * $dy);
if ($mag > 0) {
$dx /= $mag;
$dy /= $mag;
}
// Calculate dot product, projecting onto normalized direction. $pvx = $pt[0] - $line[0][0]; $pvy = $pt[1] - $line[0][1]; $pvdot = $dx * $pvx + $dy * $pvy;
// Scale line direction vector and subtract from pv. $dsx = $pvdot * $dx; $dsy = $pvdot * $dy; $ax = $pvx - $dsx; $ay = $pvy - $dsy;
return sqrt($ax * $ax + $ay * $ay);
}
function ramer_douglas_peucker(array $points, $epsilon) {
if (count($points) < 2) {
throw new InvalidArgumentException('Not enough points to simplify');
}
// Find the point with the maximum distance from the line between start/end.
$dmax = 0;
$index = 0;
$end = count($points) - 1;
$start_end_line = [$points[0], $points[$end]];
for ($i = 1; $i < $end; $i++) {
$dist = perpendicular_distance($points[$i], $start_end_line);
if ($dist > $dmax) {
$index = $i;
$dmax = $dist;
}
}
// If max distance is larger than epsilon, recursively simplify.
if ($dmax > $epsilon) {
$new_start = ramer_douglas_peucker(array_slice($points, 0, $index + 1), $epsilon);
$new_end = ramer_douglas_peucker(array_slice($points, $index), $epsilon);
array_pop($new_start);
return array_merge($new_start, $new_end);
}
// Max distance is below epsilon, so return a line from with just the // start and end points. return [ $points[0], $points[$end]];
}
$polyline = [
[0,0], [1,0.1], [2,-0.1], [3,5], [4,6], [5,7], [6,8.1], [7,9], [8,9], [9,9],
];
$result = ramer_douglas_peucker($polyline, 1.0); print "Result:\n"; foreach ($result as $point) {
print $point[0] . ',' . $point[1] . "\n";
}</lang>
- Output:
Result: 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))
Raku
(formerly 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))
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 0 /*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; call bld space( translate(arg(1), , ')(][}{') )
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)
Rust
An implementation of the algorithm
<lang rust>#[derive(Copy, Clone)] struct Point {
x: f64, y: f64,
}
use std::fmt;
impl fmt::Display for Point {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "({}, {})", self.x, self.y)
}
}
// Returns the distance from point p to the line between p1 and p2 fn perpendicular_distance(p: &Point, p1: &Point, p2: &Point) -> f64 {
let dx = p2.x - p1.x; let dy = p2.y - p1.y; (p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x).abs() / dx.hypot(dy)
}
fn rdp(points: &[Point], epsilon: f64, result: &mut Vec<Point>) {
let n = points.len();
if n < 2 {
return;
}
let mut max_dist = 0.0;
let mut index = 0;
for i in 1..n - 1 {
let dist = perpendicular_distance(&points[i], &points[0], &points[n - 1]);
if dist > max_dist {
max_dist = dist;
index = i;
}
}
if max_dist > epsilon {
rdp(&points[0..=index], epsilon, result);
rdp(&points[index..n], epsilon, result);
} else {
result.push(points[n - 1]);
}
}
fn ramer_douglas_peucker(points: &[Point], epsilon: f64) -> Vec<Point> {
let mut result = Vec::new();
if points.len() > 0 && epsilon >= 0.0 {
result.push(points[0]);
rdp(points, epsilon, &mut result);
}
result
}
fn main() {
let points = vec![
Point { x: 0.0, y: 0.0 },
Point { x: 1.0, y: 0.1 },
Point { x: 2.0, y: -0.1 },
Point { x: 3.0, y: 5.0 },
Point { x: 4.0, y: 6.0 },
Point { x: 5.0, y: 7.0 },
Point { x: 6.0, y: 8.1 },
Point { x: 7.0, y: 9.0 },
Point { x: 8.0, y: 9.0 },
Point { x: 9.0, y: 9.0 },
];
for p in ramer_douglas_peucker(&points, 1.0) {
println!("{}", p);
}
}</lang>
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
Using the geo crate
The geo crate contains an implementation of the Ramer-Douglas-Peucker line simplification algorithm. <lang rust>// [dependencies] // geo = "0.14"
use geo::algorithm::simplify::Simplify; use geo::line_string;
fn main() {
let points = line_string![
(x: 0.0, y: 0.0),
(x: 1.0, y: 0.1),
(x: 2.0, y: -0.1),
(x: 3.0, y: 5.0),
(x: 4.0, y: 6.0),
(x: 5.0, y: 7.0),
(x: 6.0, y: 8.1),
(x: 7.0, y: 9.0),
(x: 8.0, y: 9.0),
(x: 9.0, y: 9.0),
];
for p in points.simplify(&1.0) {
println!("({}, {})", p.x, p.y);
}
}</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]]
Swift
<lang swift>struct Point: CustomStringConvertible {
let x: Double, y: Double
var description: String {
return "(\(x), \(y))"
}
}
// Returns the distance from point p to the line between p1 and p2 func perpendicularDistance(p: Point, p1: Point, p2: Point) -> Double {
let dx = p2.x - p1.x let dy = p2.y - p1.y let d = (p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x) return abs(d)/(dx * dx + dy * dy).squareRoot()
}
func ramerDouglasPeucker(points: [Point], epsilon: Double) -> [Point] {
var result : [Point] = []
func rdp(begin: Int, end: Int) {
guard end > begin else {
return
}
var maxDist = 0.0
var index = 0
for i in begin+1..<end {
let dist = perpendicularDistance(p: points[i], p1: points[begin],
p2: points[end])
if dist > maxDist {
maxDist = dist
index = i
}
}
if maxDist > epsilon {
rdp(begin: begin, end: index)
rdp(begin: index, end: end)
} else {
result.append(points[end])
}
}
if points.count > 0 && epsilon >= 0.0 {
result.append(points[0])
rdp(begin: 0, end: points.count - 1)
}
return result
}
let points = [
Point(x: 0.0, y: 0.0), Point(x: 1.0, y: 0.1), Point(x: 2.0, y: -0.1), Point(x: 3.0, y: 5.0), Point(x: 4.0, y: 6.0), Point(x: 5.0, y: 7.0), Point(x: 6.0, y: 8.1), Point(x: 7.0, y: 9.0), Point(x: 8.0, y: 9.0), Point(x: 9.0, y: 9.0)
] print("\(ramerDouglasPeucker(points: points, epsilon: 1.0))")</lang>
- Output:
[(0.0, 0.0), (2.0, -0.1), (3.0, 5.0), (7.0, 9.0), (9.0, 9.0)]
Wren
<lang ecmascript>import "/dynamic" for Tuple
var Point = Tuple.create("Point", ["x", "y"])
var rdp // recursive rdp = Fn.new { |l, eps|
var x = 0
var dMax = -1
var p1 = l[0]
var p2 = l[-1]
var x21 = p2.x - p1.x
var y21 = p2.y - p1.y
var i = 0
for (p in l[1..-1]) {
var d = (y21*p.x - x21*p.y + p2.x*p1.y - p2.y*p1.x).abs
if (d > dMax) {
x = i + 1
dMax = d
}
i = i + 1
}
if (dMax > eps) {
return rdp.call(l[0..x], eps) + rdp.call(l[x..-1], eps)[1..-1]
}
return [l[0], l[-1]]
}
var points = [
Point.new(0, 0), Point.new(1, 0.1), Point.new(2, -0.1), Point.new(3, 5), Point.new(4, 6), Point.new(5, 7), Point.new(6, 8.1), Point.new(7, 9), Point.new(8, 9), Point.new(9, 9)
] System.print(rdp.call(points, 1))</lang>
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
Yabasic
<lang Yabasic>sub perpendicularDistance(tabla(), i, ini, fin)
local dx, cy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay
dx = tabla(fin, 1) - tabla(ini, 1) dy = tabla(fin, 2) - tabla(ini, 2) //Normalise mag = (dx^2 + dy^2)^0.5 if mag > 0 dx = dx / mag : dy = dy / mag pvx = tabla(i, 1) - tabla(ini, 1) pvy = tabla(i, 2) - tabla(ini, 2) //Get dot product (project pv onto normalized direction) pvdot = dx * pvx + dy * pvy //Scale line direction vector dsx = pvdot * dx dsy = pvdot * dy //Subtract this from pv ax = pvx - dsx ay = pvy - dsy return (ax^2 + ay^2)^0.5
end sub
sub DouglasPeucker(PointList(), ini, fin, epsilon)
local dmax, index, i, d // Find the point with the maximum distance
for i = ini + 1 to fin
d = perpendicularDistance(PointList(), i, ini, fin)
if d > dmax index = i : dmax = d
next
// If max distance is greater than epsilon, recursively simplify
if dmax > epsilon then
PointList(index, 3) = true
// Recursive call
DouglasPeucker(PointList(), ini, index, epsilon)
DouglasPeucker(PointList(), index, fin, epsilon)
end if
end sub
data 0,0, 1,0.1, 2,-0.1, 3,5, 4,6, 5,7, 6,8.1, 7,9, 8,9, 9,9
dim matriz(10, 3)
for i = 1 to 10
read matriz(i, 1), matriz(i, 2)
next
DouglasPeucker(matriz(), 1, 10, 1)
matriz(1, 3) = true : matriz(10, 3) = true for i = 1 to 10
if matriz(i, 3) print matriz(i, 1), matriz(i, 2)
next</lang>
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