Ramer-Douglas-Peucker line simplification
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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.
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))
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
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;
}
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
C#
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));
}
}
}
- Output:
Points remaining after simplification: (0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
C++
#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
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];
}
}
- Output:
result 0,0 2,-0.1 3,5 7,9 9,9
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
- Output:
Puntos restantes tras de simplificar: ( 0, 0) ( 2, -0.1) ( 3, 5) ( 7, 9) ( 9, 9)
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))
}
- Output:
[{0 0} {2 -0.1} {3 5} {7 9} {9 9}]
J
Solution:
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.
)
Example Usage:
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
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);
}
}
- 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
/**
* @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));
- Output:
[{x: 0, y: 0}, {x: 2, y: -0.1}, {x: 3, y: 5}, {x: 7, y: 9}, {x: 9, y: 9}]
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)
- 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
// 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)
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)
- 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
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);
- Output:
ECHO: [[0, 0], [2, -0.1], [3, 5], [7, 9], [9, 9]]
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] )
- Output:
(0 0) (2 -0.1) (3 5) (7 9) (9 9)
Phix
with javascript_semantics 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)
- Output:
{{0,0},{2,-0.1},{3,5},{7,9},{9,9}}
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";
}
- Output:
Result: 0,0 2,-0.1 3,5 7,9 9,9
Python
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
#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))
Raku
(formerly Perl 6)
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))
REXX
The computation for the perpendicular distance was taken from the GO example.
/*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 @.#
- 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
#[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);
}
}
- 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.
// [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);
}
}
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
Sidef
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.round(-20)
}
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.first(i), ε).first(-1)...,
Ramer_Douglas_Peucker(points.slice(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]]
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))")
- Output:
[(0.0, 0.0), (2.0, -0.1), (3.0, 5.0), (7.0, 9.0), (9.0, 9.0)]
Wren
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))
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
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
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]);
}
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