Ramer-Douglas-Peucker line simplification: Difference between revisions
Ramer-Douglas-Peucker line simplification (view source)
Revision as of 15:28, 29 January 2024
, 3 months ago→{{header|Wren}}: Minor tidy
Thundergnat (talk | contribs) (→{{header|Perl 6}}: Add Perl 6 example) |
m (→{{header|Wren}}: Minor tidy) |
||
| (74 intermediate revisions by 28 users not shown) | |||
Line 1:
[[Category:Geometry]]
[[Category:Recursion]]
{{task}}
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: [https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm Ramer-Douglas-Peucker algorithm].
<br><br>
=={{header|11l}}==
{{trans|Go}}
<syntaxhighlight 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))</syntaxhighlight>
{{out}}
<pre>
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
</pre>
=={{header|ALGOL 68}}==
{{trans|Go}}
<syntaxhighlight lang="algol68">
BEGIN # Ramer Douglas Peucker algotithm - translated from the Go sample #
MODE POINT = STRUCT( REAL x, y );
PRIO APPEND = 1;
OP APPEND = ( []POINT a, b )[]POINT: # append two POINT arrays #
BEGIN
INT len a = ( UPB a - LWB a ) + 1;
INT len b = ( UPB b - LWB b ) + 1;
[ 1 : lena + len b ]POINT result;
result[ : len a ] := a;
result[ len a + 1 : ] := b;
result
END # APPEND # ;
PROC rdp = ( []POINT l, REAL epsilon )[]POINT: # Ramer Douglas Peucker #
BEGIN # algorithm #
INT x := 0;
REAL d max := -1;
POINT p1 = l[ LWB l ];
POINT p2 = l[ UPB l ];
REAL x21 = x OF p2 - x OF p1;
REAL y21 = y OF p2 - y OF p1;
REAL p2xp1y = x OF p2 * y OF p1;
REAL p2yp1x = y OF p2 * x OF p1;
FOR i FROM LWB l TO UPB l DO
POINT p = l[ i ];
IF REAL d := ABS ( y21 * x OF p - x21 * y OF p + p2xp1y - p2yp1x); d > d max
THEN
x := i;
d max := d
FI
OD;
IF d max > epsilon THEN
rdp( l[ : x ], epsilon ) APPEND rdp( l[ x + 1 : ], epsilon )
ELSE
[]POINT( l[ LWB l ], l[ UPB l ] )
FI
END # rdp # ;
OP FMT = ( REAL v )STRING: # formsts v with up to 2 decimals #
BEGIN
STRING result := fixed( ABS v, 0, 2 );
IF result[ LWB result ] = "." THEN "0" +=: result FI;
WHILE result[ UPB result ] = "0" DO result := result[ : UPB result - 1 ] OD;
IF result[ UPB result ] = "." THEN result := result[ : UPB result - 1 ] FI;
IF v < 0 THEN "-" + result ELSE result FI
END # FMT # ;
OP SHOW = ( []POINT a )VOID: # prints an array of points #
BEGIN
print( ( "[" ) );
FOR i FROM LWB a TO UPB a DO
print( ( " ( ", FMT x OF a[ i ], ", ", FMT y OF a[ i ], " )" ) )
OD;
print( ( " ]" ) )
END # SHOW # ;
SHOW rdp( ( ( 0, 0 ), ( 1, 0.1 ), ( 2, -0.1 ), ( 3, 5 ), ( 4, 6 )
, ( 5, 7 ), ( 6, 8.1 ), ( 7, 9 ), ( 8, 9 ), ( 9, 9 )
)
, 1
)
END
</syntaxhighlight>
{{out}}
<pre>
[ ( 0, 0 ) ( 2, -0.1 ) ( 3, 5 ) ( 7, 9 ) ( 8, 9 ) ( 9, 9 ) ]
</pre>
=={{header|C}}==
<syntaxhighlight 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;
}</syntaxhighlight>
{{out}}
<pre>
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
</pre>
=={{header|C sharp|C#}}==
{{trans|Java}}
<syntaxhighlight 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));
}
}
}</syntaxhighlight>
{{out}}
<pre>Points remaining after simplification:
(0, 0)
(2, -0.1)
(3, 5)
(7, 9)
(9, 9)</pre>
=={{header|C++}}==
<
#include <cmath>
#include <utility>
Line 117 ⟶ 415:
return 0;
}</
{{out}}
Line 127 ⟶ 425:
9,9</pre>
=={{header|
{{trans|C++}}
<syntaxhighlight lang="d">import std.algorithm;
import std.exception : enforce;
import std.math;
import std.stdio;
void main() {
creal[] pointList = [
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) {
//Normalise
real mag = hypot(d.re, d.im);
d
}
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
output.length = 0;
output ~= pointList[0];
output ~= pointList[end];
}
}</syntaxhighlight>
{{out}}
<pre>result
0,0
2,-0.1
3,5
7,9
9,9</pre>
=={{header|EasyLang}}==
{{trans|C}}
<syntaxhighlight>
func perpdist px py p1x p1y p2x p2y .
dx = p2x - p1x
dy = p2y - p1y
d = sqrt (dx * dx + dy * dy)
return abs (px * dy - py * dx + p2x * p1y - p2y * p1x) / d
.
proc peucker eps . ptx[] pty[] outx[] outy[] .
n = len ptx[]
idx = 1
for i = 2 to n - 1
dist = perpdist ptx[i] pty[i] ptx[1] pty[1] ptx[n] pty[n]
if dist > dmax
dmax = dist
idx = i
.
.
if dmax > eps
for i to idx
px[] &= ptx[i]
py[] &= pty[i]
.
peucker eps px[] py[] outx[] outy[]
#
for i = idx to n
p2x[] &= ptx[i]
p2y[] &= pty[i]
.
peucker eps p2x[] p2y[] ox[] oy[]
for i = 2 to len ox[]
outx[] &= ox[i]
outy[] &= oy[i]
.
else
outx[] = [ ptx[1] ptx[n] ]
outy[] = [ pty[1] pty[n] ]
.
.
proc prpts px[] py[] . .
for i to len px[]
write "(" & px[i] & " " & py[i] & ") "
.
print ""
.
px[] = [ 0 1 2 3 4 5 6 7 8 9 ]
py[] = [ 0 0.1 -0.1 5 6 7 8.1 9 9 9 ]
peucker 1 px[] py[] ox[] oy[]
prpts ox[] oy[]
</syntaxhighlight>
=={{header|FreeBASIC}}==
{{trans|Yabasic}}
<syntaxhighlight 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
</syntaxhighlight>
{{out}}
<pre>
Puntos restantes tras de simplificar:
( 0, 0) ( 2, -0.1) ( 3, 5) ( 7, 9) ( 9, 9)
</pre>
=={{header|Go}}==
<syntaxhighlight 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,
{3, 5}, {4, 6}, {5, 7}, {6, 8.1}, {7, 9}, {8, 9}, {9, 9}}, 1))
}</syntaxhighlight>
{{out}}
<pre>
[{0 0} {2 -0.1} {3 5} {7 9} {9 9}]
</pre>
=={{header|J}}==
'''Solution:'''
<syntaxhighlight 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.
)</syntaxhighlight>
'''Example Usage:'''
<syntaxhighlight 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</syntaxhighlight>
=={{header|Java}}==
{{trans|Kotlin}}
{{works with|Java|9}}
<syntaxhighlight 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);
}
}</syntaxhighlight>
{{out}}
<pre>Points remaining after simplification:
(0.000000, 0.000000)
(2.000000, -0.100000)
(3.000000, 5.000000)
(7.000000, 9.000000)
(9.000000, 9.000000)</pre>
=={{header|JavaScript}}==
{{trans|Go}}
<syntaxhighlight 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));</syntaxhighlight>
{{out}}
<pre>[{x: 0, y: 0},
{x: 2, y: -0.1},
{x: 3, y: 5},
{x: 7, y: 9},
{x: 9, y: 9}]</pre>
=={{header|Julia}}==
{{works with|Julia|0.6}}
{{trans|C++}}
<syntaxhighlight 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)</syntaxhighlight>
{{out}}
<pre>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]]</pre>
=={{header|Kotlin}}==
{{trans|C++}}
<syntaxhighlight 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)
}</syntaxhighlight>
{{out}}
<pre>
Points remaining after simplification:
(0.0, 0.0)
(2.0, -0.1)
(3.0, 5.0)
(7.0, 9.0)
(9.0, 9.0)
</pre>
=={{header|Nim}}==
<syntaxhighlight 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)</syntaxhighlight>
{{out}}
<pre>@[(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)]</pre>
=={{header|ooRexx}}==
<syntaxhighlight lang="oorexx">
/*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 pl /*obtain optional arguments from the CL*/
If epsilon='' | epsilon=',' then epsilon= 1 /*Not specified? Then use the default.*/
If pl='' Then pl= '(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)'
Say ' error threshold:' epsilon
Say ' points specified:' pl
Say 'points simplified:' dlp(pl)
Exit
dlp: Procedure Expose epsilon
Parse Arg pl
plc=pl
Do i=1 By 1 While plc<>''
Parse Var plc '(' X ',' y ')' plc
p.i=.point~new(x,y)
End
end=i-1
dmax=0
index=0
Do i=2 To end-1
d=distpg(p.i,.line~new(p.1,p.end))
If d>dmax Then Do
index=i
dmax=d
End
End
If dmax>epsilon Then Do
rla=dlp(subword(pl,1,index))
rlb=dlp(subword(pl,index,end))
rl=subword(rla,1,words(rla)-1) rlb
End
Else
rl=word(pl,1) word(pl,end)
Return rl
::CLASS point public
::ATTRIBUTE X
::ATTRIBUTE Y
::METHOD init PUBLIC
Expose X Y
Use Arg X,Y
::CLASS line public
::ATTRIBUTE A
::ATTRIBUTE B
::METHOD init PUBLIC
Expose A B
Use Arg A,B
If A~x=B~x & A~y=B~y Then Do
Say 'not a line'
Return .nil
End
::METHOD k PUBLIC
Expose A B
ax=A~x; ay=A~y; bx=B~x; by=B~y
If ax=bx Then
res='infinite'
Else
res=(by-ay)/(bx-ax)
Return res
::METHOD d PUBLIC
Expose A B
ax=A~x; ay=A~y
If self~k='infinite' Then
res='indeterminate'
Else
res=round(ay-ax*self~k)
Return res
::ROUTINE distpg PUBLIC --Compute the distance from a point to a line
/***********************************************************************
* Compute the distance from a point to a line
***********************************************************************/
Use Arg A,g
ax=A~x; ay=A~y
k=g~k
If k='infinite' Then Do
Parse Value g~kxd With 'x='gx
res=gx-ax
End
Else
res=(ay-k*ax-g~d)/rxCalcsqrt(1+k**2)
Return abs(res)
::ROUTINE round PUBLIC --Round a number to 3 decimal digits
/***********************************************************************
* Round a nmber to 3 decimal digits
***********************************************************************/
Use Arg z,d
Numeric Digits 30
res=z+0
If d>'' Then
res=format(res,9,6)
Return strip(res)
::requires rxMath library
</syntaxhighlight>
{{out}}
<pre> 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)</pre>
=={{header|Openscad}}==
{{trans|C++}}
{{works with|Openscad|2015.03}}
{{works with|Openscad|2019.05}}
<syntaxhighlight 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);</syntaxhighlight>
{{out}}
<pre>ECHO: [[0, 0], [2, -0.1], [3, 5], [7, 9], [9, 9]]</pre>
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight 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] )</syntaxhighlight>
{{out}}
<pre>(0 0)
(2 -0.1)
(3 5)
(7 9)
(9 9)</pre>
=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rdp</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">2</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"not enough points to simplify"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">idx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">0</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">dMax</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">p1x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p1y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">p2x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p2y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">[$],</span>
<span style="color: #000000;">x21</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">p2x</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">p1x</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">y21</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">p2y</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">p1y</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">px</span><span style="color: #0000FF;">,</span><span style="color: #000000;">py</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y21</span><span style="color: #0000FF;">*</span><span style="color: #000000;">px</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x21</span><span style="color: #0000FF;">*</span><span style="color: #000000;">py</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">p2x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">p1y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p2y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">p1x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">></span><span style="color: #000000;">dMax</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">idx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span>
<span style="color: #000000;">dMax</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">dMax</span><span style="color: #0000FF;">></span><span style="color: #000000;">e</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rdp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">&</span> <span style="color: #000000;">rdp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">..$],</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">[$]}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">points</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">0.1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8.1</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">}}</span>
<span style="color: #0000FF;">?</span><span style="color: #000000;">rdp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">points</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
{{0,0},{2,-0.1},{3,5},{7,9},{9,9}}
</pre>
=={{header|PHP}}==
{{trans|C++}}
<syntaxhighlight 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";
}</syntaxhighlight>
{{out}}
<pre>
Result:
0,0
2,-0.1
3,5
7,9
9,9
</pre>
=={{header|Python}}==
{{libheader|Shapely}}
An approach using the shapely library:
<syntaxhighlight 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=
{{out}}
<pre>LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)</pre>
=={{header|Racket}}==
<syntaxhighlight 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)))
;; 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))</syntaxhighlight>
{{out}}
<pre>'((0 0) (2 -0.1) (3 5) (7 9) (9 9))</pre>
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2017.05}}
{{trans|C++}}
<syntaxhighlight lang="raku" line>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)]
);</syntaxhighlight>
{{out}}
<pre>((0 0) (2 -0.1) (3 5) (7 9) (9 9))</pre>
=={{header|REXX}}==
===Version 1===
The computation for the ''perpendicular distance'' was taken from
the '''GO''' example.
<syntaxhighlight 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 @.#</syntaxhighlight>
{{out|output|text= when using the default inputs:}}
<pre>
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)
</pre>
===Version 2===
{{trans|ooRexx}}
<syntaxhighlight 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 pl /*obtain optional arguments from the CL*/
If epsilon='' | epsilon=',' then epsilon= 1 /*Not specified? Then use the default.*/
If pl='' Then pl= '(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)'
Say ' error threshold:' epsilon
Say ' points specified:' pl
Say 'points simplified:' dlp(pl)
Exit
dlp: Procedure Expose epsilon
Parse Arg pl
plc=pl
Do i=1 By 1 While plc<>''
Parse Var plc '(' x ',' y ')' plc
p.i=x y
End
end=i-1
dmax=0
index=0
Do i=2 To end-1
d=distpg(p.i,p.1,p.end)
If d>dmax Then Do
index=i
dmax=d
End
End
If dmax>epsilon Then Do
rla=dlp(subword(pl,1,index))
rlb=dlp(subword(pl,index,end))
rl=subword(rla,1,words(rla)-1) rlb
End
Else
rl=word(pl,1) word(pl,end)
Return rl
distpg: Procedure
/**********************************************************************
* compute the distance of point P from the line giveb by A and B
**********************************************************************/
Parse Arg P,A,B
Parse Var P px py
Parse Var A ax ay
Parse Var B bx by
If ax=bx Then
res=px-ax
Else Do
k=(by-ay)/(bx-ax)
d=(ay-ax*k)
res=(py-k*px-d)/sqrt(1+k**2)
End
Return abs(res)
sqrt: Procedure
/* REXX ***************************************************************
* EXEC to calculate the square root of a = 2 with high precision
**********************************************************************/
Parse Arg x,prec
If prec<9 Then prec=9
prec1=2*prec
eps=10**(-prec1)
k = 1
Numeric Digits 3
r0= x
r = 1
Do i=1 By 1 Until r=r0 | (abs(r*r-x)<eps)
r0 = r
r = (r + x/r) / 2
k = min(prec1,2*k)
Numeric Digits (k + 5)
End
Numeric Digits prec
r=r+0
Return r
</syntaxhighlight>
{{out|output|text= when using the default inputs:}}
<pre>
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)
</pre>
=={{header|RPL}}==
{{works with|Halcyon Calc|4.2.9}}
≪
3 PICK - CONJ SWAP ROT
- (0,1) * DUP ABS /
* RE ABS
≫ ≫ '<span style="color:blue">DM→AB</span>' STO
≪ OVER SIZE 0 → points epsilon nb dmax
≪ '''IF''' nb 2 ≤ '''THEN''' points '''ELSE'''
0 points 1 GET points nb GET
2 nb 1 - '''FOR''' j
DUP2 points j GET <span style="color:blue">DM→AB</span>
'''IF''' DUP dmax < '''THEN''' DROP '''ELSE'''
'dmax' STO ROT DROP j ROT ROT '''END'''
'''NEXT'''
'''IF''' dmax epsilon < '''THEN'''
{ } + + SWAP DROP
'''ELSE'''
DROP2 points 1 3 PICK SUB epsilon <span style="color:blue">RDP</span>
points ROT nb SUB epsilon <span style="color:blue">RDP</span> 2 nb SUB +
'''END END'''
≫ ≫ '<span style="color:blue">RDP</span>' STO
{ (0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9) } 1 <span style="color:blue">RDP</span>
{{out}}
<pre>
1: { (0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9) }
</pre>
=={{header|Rust}}==
===An implementation of the algorithm===
<syntaxhighlight 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);
}
}</syntaxhighlight>
{{out}}
<pre>
(0, 0)
(2, -0.1)
(3, 5)
(7, 9)
(9, 9)
</pre>
===Using the geo crate===
The geo crate contains an implementation of the Ramer-Douglas-Peucker line
simplification algorithm.
<syntaxhighlight 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);
}
}</syntaxhighlight>
{{out}}
<pre>
(0, 0)
(2, -0.1)
(3, 5)
(7, 9)
(9, 9)
</pre>
=={{header|Scala}}==
{{trans|Kotlin}}
<syntaxhighlight lang="scala">// version 1.1.0
object SimplifyPolyline extends App {
type Point = (Double, Double)
def perpendicularDistance(
pt: Point,
lineStart: Point,
lineEnd: Point
): Double = {
var dx = lineEnd._1 - lineStart._1
var dy = lineEnd._2 - lineStart._2
// Normalize
val mag = math.hypot(dx, dy)
if (mag > 0.0) { dx /= mag; dy /= mag }
val pvx = pt._1 - lineStart._1
val pvy = pt._2 - lineStart._2
// Get dot product (project pv onto normalized direction)
val pvdot = dx * pvx + dy * pvy
// Scale line direction vector and subtract it from pv
val ax = pvx - pvdot * dx
val ay = pvy - pvdot * dy
math.hypot(ax, ay)
}
def RamerDouglasPeucker(
pointList: List[Point],
epsilon: Double
): List[Point] = {
if (pointList.length < 2)
throw new IllegalArgumentException("Not enough points to simplify")
// Check if there are points to process
if (pointList.length > 2) {
// Find the point with the maximum distance from the line between the first and last
val (dmax, index) = pointList.zipWithIndex
.slice(1, pointList.length - 1)
.map { case (point, i) =>
(perpendicularDistance(point, pointList.head, pointList.last), i)
}
.maxBy(_._1)
// If max distance is greater than epsilon, recursively simplify
if (dmax > epsilon) {
val firstLine = pointList.take(index + 1)
val lastLine = pointList.drop(index)
val recResults1 = RamerDouglasPeucker(firstLine, epsilon)
val recResults2 = RamerDouglasPeucker(lastLine, epsilon)
// Combine the results
(recResults1.init ++ recResults2).distinct
} else {
// If no point is further than epsilon, return the endpoints
List(pointList.head, pointList.last)
}
} else {
// If there are only two points, just return them
pointList
}
}
val pointList = List(
(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)
)
val simplifiedPointList = RamerDouglasPeucker(pointList, 1.0)
println("Points remaining after simplification:")
simplifiedPointList.foreach(println)
}
</syntaxhighlight>
{{out}}
<pre>
Points remaining after simplification:
(0.0,0.0)
(2.0,-0.1)
(3.0,5.0)
(7.0,9.0)
(9.0,9.0)
</pre>
=={{header|Sidef}}==
{{trans|Raku}}
<syntaxhighlight 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.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]]
)</syntaxhighlight>
{{out}}
<pre>
[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]]
</pre>
=={{header|Swift}}==
<syntaxhighlight 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))")</syntaxhighlight>
{{out}}
<pre>
[(0.0, 0.0), (2.0, -0.1), (3.0, 5.0), (7.0, 9.0), (9.0, 9.0)]
</pre>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-dynamic}}
<syntaxhighlight lang="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))</syntaxhighlight>
{{out}}
<pre>
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
</pre>
=={{header|XPL0}}==
{{trans|C}}
<syntaxhighlight lang "XPL0">include xpllib; \for PerpDist and Print
func RDP(Points, Size, Epsilon, RemPoints);
\Reduce array Points using Ramer-Douglas-Peucker algorithm
\Return array of remaining points in RemPoints and its size
real Points; int Size; real Epsilon, RemPoints;
real Dist, DistMax;
int Index, I, Size1, Size2;
[\Find Index of point farthest from line between first and last points
DistMax:= 0.; Index:= 0;
for I:= 1 to Size-2 do
[Dist:= PerpDist(Points(I), Points(0), Points(Size-1));
if Dist > DistMax then
[DistMax:= Dist;
Index:= I;
];
];
if DistMax <= Epsilon then \replace in-between points with first and last points
[RemPoints(0):= Points(0);
RemPoints(1):= Points(Size-1);
return 2;
]
else \recursively check sub-line-segments
[Size1:= RDP(Points, Index+1, Epsilon, RemPoints);
RemPoints:= @RemPoints(Size1-1); \add retained points to array
Size2:= RDP(@Points(Index), Size-Index, Epsilon, RemPoints);
return Size1 + Size2 - 1; \for point at Index
];
];
real Points, RemPoints;
int Size, I;
[Points:= [ [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] ];
Size:= 10;
RemPoints:= RlRes(Size);
Size:= RDP(Points, Size, 1.0, RemPoints);
for I:= 0 to Size-1 do
[if I > 0 then Print(", ");
Print("[%1.1g, %1.1g]", RemPoints(I,0), RemPoints(I,1));
];
CrLf(0);
]</syntaxhighlight>
{{out}}
<pre>
[0, 0], [2, -0.1], [3, 5], [7, 9], [9, 9]
</pre>
=={{header|Yabasic}}==
<syntaxhighlight 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</syntaxhighlight>
=={{header|zkl}}==
{{trans|Raku}}
<syntaxhighlight 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]);
}</syntaxhighlight>
<syntaxhighlight 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();</syntaxhighlight>
{{out}}
<pre>
L(L(0,0),L(2,-0.1),L(3,5),L(7,9),L(9,9))
</pre>
'''References'''
<references/>
| |||