Ramer-Douglas-Peucker line simplification: Difference between revisions
(Added C solution) |
|||
(31 intermediate revisions by 14 users not shown) | |||
Line 18: | Line 18: | ||
:* the Wikipedia article: [https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm Ramer-Douglas-Peucker algorithm]. |
:* the Wikipedia article: [https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm Ramer-Douglas-Peucker algorithm]. |
||
<br><br> |
<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|BASIC256}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="vbnet">arraybase 1 |
|||
global matriz |
|||
dim matriz(10, 3) |
|||
matriz[ 1, 1] = 0 : matriz[ 1, 2] = 0 : matriz[ 1, 3] = 0 |
|||
matriz[ 2, 1] = 1 : matriz[ 2, 2] = 0.1 : matriz[ 2, 3] = 0 |
|||
matriz[ 3, 1] = 2 : matriz[ 3, 2] = -0.1 : matriz[ 3, 3] = 0 |
|||
matriz[ 4, 1] = 3 : matriz[ 4, 2] = 5 : matriz[ 4, 3] = 0 |
|||
matriz[ 5, 1] = 4 : matriz[ 5, 2] = 6 : matriz[ 5, 3] = 0 |
|||
matriz[ 6, 1] = 5 : matriz[ 6, 2] = 7 : matriz[ 6, 3] = 0 |
|||
matriz[ 7, 1] = 6 : matriz[ 7, 2] = 8.1 : matriz[ 7, 3] = 0 |
|||
matriz[ 8, 1] = 7 : matriz[ 8, 2] = 9 : matriz[ 8, 3] = 0 |
|||
matriz[ 9, 1] = 8 : matriz[ 9, 2] = 9 : matriz[ 9, 3] = 0 |
|||
matriz[10, 1] = 9 : matriz[10, 2] = 9 : matriz[10, 3] = 0 |
|||
call DRDP(matriz, 1, 10, 1) |
|||
print "Remaining points after simplifying:" |
|||
matriz[1, 3] = true |
|||
matriz[10, 3] = true |
|||
for i = 1 to matriz[?][] |
|||
if matriz[i, 3] = true then print "(";matriz[i, 1];", "; matriz[i, 2];") "; |
|||
next i |
|||
end |
|||
function DistanciaPerpendicular(tabla, i, ini, fin) |
|||
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 then dx /= mag : 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 function |
|||
subroutine DRDP(matriz, ini, fin, epsilon) |
|||
dmax = 0 |
|||
# Find the point with the maximum distance |
|||
if matriz[?][] < 2 then print "Not enough points to simplify": end |
|||
for i = ini + 1 to fin |
|||
d = DistanciaPerpendicular(matriz, i, ini, fin) |
|||
if d > dmax then indice = i : dmax = d |
|||
next |
|||
# If max distance is greater than epsilon, recursively simplify |
|||
if dmax > epsilon then |
|||
matriz[indice, 3] = true |
|||
# Recursive call |
|||
call DRDP(matriz, ini, indice, epsilon) |
|||
call DRDP(matriz, indice, fin, epsilon) |
|||
end if |
|||
end subroutine</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Remaining points after simplifying: |
|||
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9) </pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <assert.h> |
||
#include <math.h> |
#include <math.h> |
||
#include <stdio.h> |
#include <stdio.h> |
||
Line 39: | Line 221: | ||
// Simplify an array of points using the Ramer–Douglas–Peucker algorithm. |
// Simplify an array of points using the Ramer–Douglas–Peucker algorithm. |
||
// Returns the number of output points. |
// Returns the number of output points. |
||
size_t |
size_t douglas_peucker(const point_t* points, size_t n, double epsilon, |
||
point_t* dest, size_t |
point_t* dest, size_t destlen) { |
||
assert(n >= 2); |
assert(n >= 2); |
||
assert(epsilon >= 0); |
assert(epsilon >= 0); |
||
Line 53: | Line 235: | ||
} |
} |
||
if (max_dist > epsilon) { |
if (max_dist > epsilon) { |
||
size_t n1 = |
size_t n1 = douglas_peucker(points, index + 1, epsilon, dest, destlen); |
||
if ( |
if (destlen >= n1 - 1) { |
||
destlen -= n1 - 1; |
|||
dest += n1 - 1; |
dest += n1 - 1; |
||
} else { |
} else { |
||
destlen = 0; |
|||
} |
} |
||
size_t n2 = |
size_t n2 = douglas_peucker(points + index, n - index, epsilon, dest, destlen); |
||
return n1 + n2 - 1; |
return n1 + n2 - 1; |
||
} |
} |
||
if ( |
if (destlen >= 2) { |
||
dest[0] = points[0]; |
dest[0] = points[0]; |
||
dest[1] = points[n - 1]; |
dest[1] = points[n - 1]; |
||
Line 86: | Line 268: | ||
const size_t len = sizeof(points)/sizeof(points[0]); |
const size_t len = sizeof(points)/sizeof(points[0]); |
||
point_t out[len]; |
point_t out[len]; |
||
size_t n = |
size_t n = douglas_peucker(points, len, 1.0, out, len); |
||
print_points(out, n); |
print_points(out, n); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 98: | Line 280: | ||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using System.Collections.Generic; |
using System.Collections.Generic; |
||
using System.Linq; |
using System.Linq; |
||
Line 187: | Line 369: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Points remaining after simplification: |
<pre>Points remaining after simplification: |
||
Line 198: | Line 380: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
< |
<syntaxhighlight lang="cpp">#include <iostream> |
||
#include <cmath> |
#include <cmath> |
||
#include <utility> |
#include <utility> |
||
Line 306: | Line 488: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 318: | Line 500: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang="d">import std.algorithm; |
||
import std.exception : enforce; |
import std.exception : enforce; |
||
import std.math; |
import std.math; |
||
Line 406: | Line 588: | ||
output ~= pointList[end]; |
output ~= pointList[end]; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 415: | Line 597: | ||
7,9 |
7,9 |
||
9,9</pre> |
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}}== |
=={{header|Go}}== |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 449: | Line 755: | ||
fmt.Println(RDP([]point{{0, 0}, {1, 0.1}, {2, -0.1}, |
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)) |
{3, 5}, {4, 6}, {5, 7}, {6, 8.1}, {7, 9}, {8, 9}, {9, 9}}, 1)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 457: | Line 763: | ||
=={{header|J}}== |
=={{header|J}}== |
||
'''Solution:''' |
'''Solution:''' |
||
< |
<syntaxhighlight lang="j">mp=: +/ .* NB. matrix product |
||
norm=: +/&.:*: NB. vector norm |
norm=: +/&.:*: NB. vector norm |
||
normalize=: (% norm)^:(0 < norm) |
normalize=: (% norm)^:(0 < norm) |
||
Line 480: | Line 786: | ||
({. ,: {:) points |
({. ,: {:) points |
||
end. |
end. |
||
)</ |
)</syntaxhighlight> |
||
'''Example Usage:''' |
'''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 |
1.0 rdp Points |
||
0 0 |
0 0 |
||
Line 488: | Line 794: | ||
3 5 |
3 5 |
||
7 9 |
7 9 |
||
9 9</ |
9 9</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
{{works with|Java|9}} |
{{works with|Java|9}} |
||
< |
<syntaxhighlight lang="java">import javafx.util.Pair; |
||
import java.util.ArrayList; |
import java.util.ArrayList; |
||
Line 587: | Line 893: | ||
pointListOut.forEach(System.out::println); |
pointListOut.forEach(System.out::println); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Points remaining after simplification: |
<pre>Points remaining after simplification: |
||
Line 598: | Line 904: | ||
=={{header|JavaScript}}== |
=={{header|JavaScript}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<syntaxhighlight lang="javascript">/** |
|||
<lang JavaScript>/** |
|||
* @typedef {{ |
* @typedef {{ |
||
* x: (!number), |
* x: (!number), |
||
Line 642: | Line 948: | ||
{x: 9, y: 9}]; |
{x: 9, y: 9}]; |
||
console.log(RDP(points, 1));</ |
console.log(RDP(points, 1));</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>[{x: 0, y: 0}, |
<pre>[{x: 0, y: 0}, |
||
Line 649: | Line 955: | ||
{x: 7, y: 9}, |
{x: 7, y: 9}, |
||
{x: 9, y: 9}]</pre> |
{x: 9, y: 9}]</pre> |
||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
'''Works with jq, the C implementation of jq''' |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''Works with jaq, the Rust implementation of jq''' |
|||
<syntaxhighlight lang="jq"> |
|||
def Point($x;$y): {x:$x, y:$y}; # jq and gojq allow: {$x,$y} |
|||
def ppArray: |
|||
def pp: "(\(.x), \(.y))"; |
|||
"[" + (map(pp) | join(", ")) + "]"; |
|||
def rdp($eps): |
|||
. as $l |
|||
| {x: 0, |
|||
dMax: -1} |
|||
| .p1 = $l[0] # first |
|||
| .p2 = $l[-1] # last |
|||
| (.p2.x - .p1.x) as $x21 |
|||
| (.p2.y - .p1.y) as $y21 |
|||
| reduce range(1; ($l|length)) as $i (.; |
|||
.p = $l[$i] |
|||
| ( ($y21*.p.x - $x21*.p.y + .p2.x*.p1.y - .p2.y*.p1.x) | length) as $d # abs ~ length |
|||
| if $d > .dMax then .x += 1 | .dMax = $d end ) |
|||
| if .dMax > $eps |
|||
then ($l[0: 1+.x]|rdp($eps)) + ($l[.x:]|rdp(eps))[1:] |
|||
else [$l[0], $l[-1]] |
|||
end; |
|||
def points: [ |
|||
Point(0; 0), Point(1; 0.1), Point(2; -0.1), Point(3; 5), Point(4; 6), |
|||
Point(5; 7), Point(6; 8.1), Point(7; 9), Point(8; 9), Point(9; 9) |
|||
]; |
|||
points | rdp(1) | ppArray |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)] |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 654: | Line 1,005: | ||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang="julia">const Point = Vector{Float64} |
||
function perpdist(pt::Point, lnstart::Point, lnend::Point) |
function perpdist(pt::Point, lnstart::Point, lnend::Point) |
||
Line 695: | Line 1,046: | ||
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]] |
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 plist |
||
@show rdp(plist)</ |
@show rdp(plist)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 703: | Line 1,054: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang="scala">// version 1.1.0 |
||
typealias Point = Pair<Double, Double> |
typealias Point = Pair<Double, Double> |
||
Line 778: | Line 1,129: | ||
println("Points remaining after simplification:") |
println("Points remaining after simplification:") |
||
for (p in pointListOut) println(p) |
for (p in pointListOut) println(p) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 791: | Line 1,142: | ||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
< |
<syntaxhighlight lang="nim">import math |
||
type |
type |
||
Line 828: | Line 1,179: | ||
var line: seq[Point] = @[(0.0, 0.0), (1.0, 0.1), (2.0, -0.1), (3.0, 5.0), (4.0, 6.0), |
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)] |
(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)</ |
echo rdp(line, line.low, line.high)</syntaxhighlight> |
||
{{out}} |
{{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> |
<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}}== |
=={{header|Openscad}}== |
||
Line 838: | Line 1,296: | ||
{{works with|Openscad|2019.05}} |
{{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 |
// Find the distance from the point to the line |
||
Line 913: | Line 1,371: | ||
$fn=16; |
$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]]; |
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);</ |
demo(points);</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 920: | Line 1,378: | ||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
use feature 'say'; |
use feature 'say'; |
||
Line 963: | Line 1,421: | ||
say '(' . join(' ', @$_) . ') ' |
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] )</ |
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}} |
{{out}} |
||
<pre>(0 0) |
<pre>(0 0) |
||
Line 973: | Line 1,431: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>function rdp(sequence l, atom e) |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
if length(l)<2 then crash("not enough points to simplify") end if |
|||
<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> |
|||
integer idx := 0 |
|||
<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> |
|||
atom dMax := -1, |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">idx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">0</span> |
|||
{p1x,p1y} := l[1], |
|||
<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> |
|||
{p2x,p2y} := l[$], |
|||
<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> |
|||
x21 := p2x - p1x, |
|||
<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> |
|||
y21 := p2y - p1y |
|||
<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> |
|||
for i=1 to length(l) do |
|||
<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> |
|||
atom {px,py} = l[i], |
|||
<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> |
|||
d = abs(y21*px-x21*py + p2x*p1y-p2y*p1x) |
|||
<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> |
|||
if d>dMax then |
|||
<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> |
|||
idx = i |
|||
<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> |
|||
dMax = d |
|||
<span style="color: #000000;">idx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span> |
|||
end if |
|||
<span style="color: #000000;">dMax</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
if dMax>e then |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
return rdp(l[1..idx], e) & rdp(l[idx..$], e)[2..$] |
|||
<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> |
|||
end if |
|||
<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> |
|||
return {l[1], l[$]} |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end function |
|||
<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> |
|||
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}} |
|||
<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> |
|||
?rdp(points, 1)</lang> |
|||
<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}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,005: | Line 1,466: | ||
=={{header|PHP}}== |
=={{header|PHP}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang="php">function perpendicular_distance(array $pt, array $line) { |
||
// Calculate the normalized delta x and y of the line. |
// Calculate the normalized delta x and y of the line. |
||
$dx = $line[1][0] - $line[0][0]; |
$dx = $line[1][0] - $line[0][0]; |
||
Line 1,077: | Line 1,538: | ||
foreach ($result as $point) { |
foreach ($result as $point) { |
||
print $point[0] . ',' . $point[1] . "\n"; |
print $point[0] . ',' . $point[1] . "\n"; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,092: | Line 1,553: | ||
{{libheader|Shapely}} |
{{libheader|Shapely}} |
||
An approach using the shapely library: |
An approach using the shapely library: |
||
< |
<syntaxhighlight lang="python">from __future__ import print_function |
||
from shapely.geometry import LineString |
from shapely.geometry import LineString |
||
if __name__=="__main__": |
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)]) |
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))</ |
print (line.simplify(1.0, preserve_topology=False))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)</pre> |
<pre>LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)</pre> |
||
=={{header| |
=={{header|QBasic}}== |
||
{{works with|QBasic|1.1}} |
|||
{{works with|QuickBasic|4.5}} |
|||
{{works with|QB64}} |
|||
<syntaxhighlight lang="qbasic">DECLARE SUB DRDP (ListaDePuntos() AS DOUBLE, ini AS INTEGER, fin AS INTEGER, epsilon AS DOUBLE) |
|||
DECLARE FUNCTION DistanciaPerpendicular! (tabla() AS DOUBLE, i AS INTEGER, ini AS INTEGER, fin AS INTEGER) |
|||
CONST True = -1 |
|||
<lang racket>#lang racket |
|||
DIM matriz(1 TO 10, 1 TO 3) AS DOUBLE |
|||
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 = LBOUND(matriz) TO UBOUND(matriz) |
|||
READ matriz(i, 1), matriz(i, 2) |
|||
NEXT i |
|||
CALL DRDP(matriz(), 1, 10, 1) |
|||
PRINT "Remaining points after simplifying:" |
|||
matriz(1, 3) = True: matriz(10, 3) = True |
|||
FOR i = LBOUND(matriz) TO UBOUND(matriz) |
|||
IF matriz(i, 3) THEN PRINT "("; matriz(i, 1); ", "; matriz(i, 2); ") "; |
|||
NEXT i |
|||
END |
|||
FUNCTION DistanciaPerpendicular (tabla() AS DOUBLE, i AS INTEGER, ini AS INTEGER, fin AS INTEGER) |
|||
DIM dx AS DOUBLE, dy AS DOUBLE, mag AS DOUBLE, pvx AS DOUBLE, pvy AS DOUBLE |
|||
DIM pvdot AS DOUBLE, dsx AS DOUBLE, dsy AS DOUBLE, ax AS DOUBLE, ay AS DOUBLE |
|||
dx = tabla(fin, 1) - tabla(ini, 1) |
|||
dy = tabla(fin, 2) - tabla(ini, 2) |
|||
'Normalise |
|||
mag = (dx ^ 2 + dy ^ 2) ^ .5 |
|||
IF mag > 0 THEN 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 |
|||
DistanciaPerpendicular = (ax ^ 2 + ay ^ 2) ^ .5 |
|||
END FUNCTION |
|||
SUB DRDP (ListaDePuntos() AS DOUBLE, ini AS INTEGER, fin AS INTEGER, epsilon AS DOUBLE) |
|||
DIM dmax AS DOUBLE, d AS DOUBLE |
|||
DIM indice AS INTEGER, i AS INTEGER |
|||
' Find the point with the maximum distance |
|||
IF UBOUND(ListaDePuntos) < 2 THEN PRINT "Not enough points to simplify": END |
|||
FOR i = ini + 1 TO fin |
|||
d = DistanciaPerpendicular(ListaDePuntos(), i, ini, fin) |
|||
IF d > dmax THEN indice = i: dmax = d |
|||
NEXT |
|||
' If max distance is greater than epsilon, recursively simplify |
|||
IF dmax > epsilon THEN |
|||
ListaDePuntos(indice, 3) = True |
|||
' Recursive call |
|||
CALL DRDP(ListaDePuntos(), ini, indice, epsilon) |
|||
CALL DRDP(ListaDePuntos(), indice, fin, epsilon) |
|||
END IF |
|||
END SUB</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Remaining points after simplifying: |
|||
( 0 , 0 ) ( 2 , -.1 ) ( 3 , 5 ) ( 7 , 9 ) ( 9 , 9 )</pre> |
|||
=={{header|QB64}}== |
|||
The [[#QBasic|QBasic]] solution works without any changes. |
|||
=={{header|Racket}}== |
|||
<syntaxhighlight lang="racket">#lang racket |
|||
(require math/flonum) |
(require math/flonum) |
||
;; points are lists of x y (maybe extensible to z) |
;; points are lists of x y (maybe extensible to z) |
||
Line 1,147: | Line 1,685: | ||
(module+ test |
(module+ test |
||
(require rackunit) |
(require rackunit) |
||
(check-= ((⊥-distance '(0 0) '(0 1)) '(1 0)) 1 epsilon.0))</ |
(check-= ((⊥-distance '(0 0) '(0 1)) '(1 0)) 1 epsilon.0))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,157: | Line 1,695: | ||
{{trans|C++}} |
{{trans|C++}} |
||
<lang |
<syntaxhighlight lang="raku" line>sub norm (*@list) { @list»².sum.sqrt } |
||
sub perpendicular-distance (@start, @end where @end !eqv @start, @point) { |
sub perpendicular-distance (@start, @end where @end !eqv @start, @point) { |
||
Line 1,181: | Line 1,719: | ||
say Ramer-Douglas-Peucker( |
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)] |
[(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}} |
{{out}} |
||
<pre>((0 0) (2 -0.1) (3 5) (7 9) (9 9))</pre> |
<pre>((0 0) (2 -0.1) (3 5) (7 9) (9 9))</pre> |
||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
===Version 1=== |
|||
The computation for the ''perpendicular distance'' was taken from the '''GO''' example. |
|||
The computation for the ''perpendicular distance'' was taken from |
|||
<lang rexx>/*REXX program uses the Ramer─Douglas─Peucker (RDP) line simplification algorithm for*/ |
|||
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. */ |
/*───────────────────────────── reducing the number of points used to define its shape. */ |
||
parse arg epsilon pts /*obtain optional arguments from the CL*/ |
parse arg epsilon pts /*obtain optional arguments from the CL*/ |
||
Line 1,193: | Line 1,733: | ||
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)' |
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. */ |
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*/ |
$= RDP(pts) /*invoke Ramer─Douglas─Peucker function*/ |
||
say 'points simplified: ' rez($) /*display points with () ───► terminal.*/ |
|||
exit |
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 |
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.*/ |
px: parse arg _; return word( translate(_, , ','), 1) /*obtain the X coörd.*/ |
||
py: parse arg _; return word( translate(_, , ','), 2) /* " " Y " */ |
py: parse arg _; return word( translate(_, , ','), 2) /* " " Y " */ |
||
reb: parse arg a,b,,_; do k=a to b; |
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(_) |
rez: parse arg z,_; do k=1 for words(z); _= _ '('word(z, k)") "; end; return strip(_) |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
RDP: procedure expose epsilon; |
RDP: procedure expose epsilon; call bld space( translate(arg(1), , ')(][}{') ) |
||
L= px(@.#)-px(@.1) |
L= px(@.#) - px(@.1) |
||
H= py(@.#)-py(@.1) |
H= py(@.#) - py(@.1) /* [↓] find point IDX with max distance*/ |
||
do i=2 to #-1 |
do i=2 to #-1 |
||
d= abs(H*px(@.i) - L*py(@.i) + px(@.#)*py(@.1) - py(@.#)*px(@.1) |
d= abs(H*px(@.i) - L*py(@.i) + px(@.#)*py(@.1) - py(@.#)*px(@.1)) |
||
if d>dMax then do; idx= i; dMax= d |
if d>dMax then do; idx= i; dMax= d |
||
end |
end |
||
end /*i*/ |
end /*i*/ /* [↑] D is the perpendicular distance*/ |
||
if dMax>epsilon then do; r= RDP( reb(1, idx) ) |
if dMax>epsilon then do; r= RDP( reb(1, idx) ) |
||
return subword(r, 1, words(r) - 1) RDP( reb(idx, #) ) |
return subword(r, 1, words(r) - 1) RDP( reb(idx, #) ) |
||
end |
end |
||
return @.1 @.#</ |
return @.1 @.#</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 1,224: | Line 1,764: | ||
points simplified: (0,0) (2,-0.1) (3,5) (7,9) (9,9) |
points simplified: (0,0) (2,-0.1) (3,5) (7,9) (9,9) |
||
</pre> |
</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}}== |
=={{header|Sidef}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="ruby">func perpendicular_distance(Arr start, Arr end, Arr point) { |
||
((point == start) || (point == end)) && return 0 |
((point == start) || (point == end)) && return 0 |
||
var (Δx, Δy ) = ( end |
var (Δx, Δy ) = ( end »-« start)... |
||
var (Δpx, Δpy) = (point |
var (Δpx, Δpy) = (point »-« start)... |
||
var h = hypot(Δx, Δy) |
var h = hypot(Δx, Δy) |
||
[\Δx, \Δy].map { *_ /= h } |
[\Δx, \Δy].map { *_ /= h } |
||
(([Δpx, Δpy] |
(([Δpx, Δpy] »-« ([Δx, Δy] »*» (Δx*Δpx + Δy*Δpy))) »**» 2).sum.sqrt.round(-20) |
||
} |
} |
||
Line 1,245: | Line 2,115: | ||
if (d.max > ε) { |
if (d.max > ε) { |
||
var i = d.index(d.max) |
var i = d.index(d.max) |
||
return [Ramer_Douglas_Peucker(points. |
return [Ramer_Douglas_Peucker(points.first(i), ε).first(-1)..., |
||
Ramer_Douglas_Peucker(points. |
Ramer_Douglas_Peucker(points.slice(i), ε)...] |
||
} |
} |
||
Line 1,254: | Line 2,124: | ||
say Ramer_Douglas_Peucker( |
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]] |
[[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}} |
{{out}} |
||
<pre> |
<pre> |
||
[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]] |
[[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> |
</pre> |
||
=={{header|Yabasic}}== |
=={{header|Yabasic}}== |
||
< |
<syntaxhighlight lang="yabasic">sub perpendicularDistance(tabla(), i, ini, fin) |
||
local dx, cy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay |
local dx, cy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay |
||
Line 1,320: | Line 2,347: | ||
for i = 1 to 10 |
for i = 1 to 10 |
||
if matriz(i, 3) print matriz(i, 1), matriz(i, 2) |
if matriz(i, 3) print matriz(i, 1), matriz(i, 2) |
||
next</ |
next</syntaxhighlight> |
||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="zkl">fcn perpendicularDistance(start,end, point){ // all are tuples: (x,y) -->|d| |
||
dx,dy := end .zipWith('-,start); // deltas |
dx,dy := end .zipWith('-,start); // deltas |
||
dpx,dpy := point.zipWith('-,start); |
dpx,dpy := point.zipWith('-,start); |
||
Line 1,342: | Line 2,369: | ||
RamerDouglasPeucker(points[index,*],epsilon))) |
RamerDouglasPeucker(points[index,*],epsilon))) |
||
} else return(points[0],points[-1]); |
} 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( 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) )) |
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();</ |
.println();</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 08:23, 12 June 2024
The Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Using the Ramer–Douglas–Peucker algorithm, simplify the 2D line defined by the points:
(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)
The error threshold to be used is: 1.0.
Display the remaining points here.
- Reference
-
- the Wikipedia article: Ramer-Douglas-Peucker algorithm.
11l
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)]
ALGOL 68
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
- Output:
[ ( 0, 0 ) ( 2, -0.1 ) ( 3, 5 ) ( 7, 9 ) ( 8, 9 ) ( 9, 9 ) ]
BASIC256
arraybase 1
global matriz
dim matriz(10, 3)
matriz[ 1, 1] = 0 : matriz[ 1, 2] = 0 : matriz[ 1, 3] = 0
matriz[ 2, 1] = 1 : matriz[ 2, 2] = 0.1 : matriz[ 2, 3] = 0
matriz[ 3, 1] = 2 : matriz[ 3, 2] = -0.1 : matriz[ 3, 3] = 0
matriz[ 4, 1] = 3 : matriz[ 4, 2] = 5 : matriz[ 4, 3] = 0
matriz[ 5, 1] = 4 : matriz[ 5, 2] = 6 : matriz[ 5, 3] = 0
matriz[ 6, 1] = 5 : matriz[ 6, 2] = 7 : matriz[ 6, 3] = 0
matriz[ 7, 1] = 6 : matriz[ 7, 2] = 8.1 : matriz[ 7, 3] = 0
matriz[ 8, 1] = 7 : matriz[ 8, 2] = 9 : matriz[ 8, 3] = 0
matriz[ 9, 1] = 8 : matriz[ 9, 2] = 9 : matriz[ 9, 3] = 0
matriz[10, 1] = 9 : matriz[10, 2] = 9 : matriz[10, 3] = 0
call DRDP(matriz, 1, 10, 1)
print "Remaining points after simplifying:"
matriz[1, 3] = true
matriz[10, 3] = true
for i = 1 to matriz[?][]
if matriz[i, 3] = true then print "(";matriz[i, 1];", "; matriz[i, 2];") ";
next i
end
function DistanciaPerpendicular(tabla, i, ini, fin)
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 then dx /= mag : 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 function
subroutine DRDP(matriz, ini, fin, epsilon)
dmax = 0
# Find the point with the maximum distance
if matriz[?][] < 2 then print "Not enough points to simplify": end
for i = ini + 1 to fin
d = DistanciaPerpendicular(matriz, i, ini, fin)
if d > dmax then indice = i : dmax = d
next
# If max distance is greater than epsilon, recursively simplify
if dmax > epsilon then
matriz[indice, 3] = true
# Recursive call
call DRDP(matriz, ini, indice, epsilon)
call DRDP(matriz, indice, fin, epsilon)
end if
end subroutine
- Output:
Remaining points after simplifying: (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
EasyLang
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[]
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}]
jq
Adapted from Wren
Works with jq, the C implementation of jq
Works with gojq, the Go implementation of jq
Works with jaq, the Rust implementation of jq
def Point($x;$y): {x:$x, y:$y}; # jq and gojq allow: {$x,$y}
def ppArray:
def pp: "(\(.x), \(.y))";
"[" + (map(pp) | join(", ")) + "]";
def rdp($eps):
. as $l
| {x: 0,
dMax: -1}
| .p1 = $l[0] # first
| .p2 = $l[-1] # last
| (.p2.x - .p1.x) as $x21
| (.p2.y - .p1.y) as $y21
| reduce range(1; ($l|length)) as $i (.;
.p = $l[$i]
| ( ($y21*.p.x - $x21*.p.y + .p2.x*.p1.y - .p2.y*.p1.x) | length) as $d # abs ~ length
| if $d > .dMax then .x += 1 | .dMax = $d end )
| if .dMax > $eps
then ($l[0: 1+.x]|rdp($eps)) + ($l[.x:]|rdp(eps))[1:]
else [$l[0], $l[-1]]
end;
def points: [
Point(0; 0), Point(1; 0.1), Point(2; -0.1), Point(3; 5), Point(4; 6),
Point(5; 7), Point(6; 8.1), Point(7; 9), Point(8; 9), Point(9; 9)
];
points | rdp(1) | ppArray
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 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)]
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
- Output:
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)
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)
QBasic
DECLARE SUB DRDP (ListaDePuntos() AS DOUBLE, ini AS INTEGER, fin AS INTEGER, epsilon AS DOUBLE)
DECLARE FUNCTION DistanciaPerpendicular! (tabla() AS DOUBLE, i AS INTEGER, ini AS INTEGER, fin AS INTEGER)
CONST True = -1
DIM matriz(1 TO 10, 1 TO 3) AS DOUBLE
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 = LBOUND(matriz) TO UBOUND(matriz)
READ matriz(i, 1), matriz(i, 2)
NEXT i
CALL DRDP(matriz(), 1, 10, 1)
PRINT "Remaining points after simplifying:"
matriz(1, 3) = True: matriz(10, 3) = True
FOR i = LBOUND(matriz) TO UBOUND(matriz)
IF matriz(i, 3) THEN PRINT "("; matriz(i, 1); ", "; matriz(i, 2); ") ";
NEXT i
END
FUNCTION DistanciaPerpendicular (tabla() AS DOUBLE, i AS INTEGER, ini AS INTEGER, fin AS INTEGER)
DIM dx AS DOUBLE, dy AS DOUBLE, mag AS DOUBLE, pvx AS DOUBLE, pvy AS DOUBLE
DIM pvdot AS DOUBLE, dsx AS DOUBLE, dsy AS DOUBLE, ax AS DOUBLE, ay AS DOUBLE
dx = tabla(fin, 1) - tabla(ini, 1)
dy = tabla(fin, 2) - tabla(ini, 2)
'Normalise
mag = (dx ^ 2 + dy ^ 2) ^ .5
IF mag > 0 THEN 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
DistanciaPerpendicular = (ax ^ 2 + ay ^ 2) ^ .5
END FUNCTION
SUB DRDP (ListaDePuntos() AS DOUBLE, ini AS INTEGER, fin AS INTEGER, epsilon AS DOUBLE)
DIM dmax AS DOUBLE, d AS DOUBLE
DIM indice AS INTEGER, i AS INTEGER
' Find the point with the maximum distance
IF UBOUND(ListaDePuntos) < 2 THEN PRINT "Not enough points to simplify": END
FOR i = ini + 1 TO fin
d = DistanciaPerpendicular(ListaDePuntos(), i, ini, fin)
IF d > dmax THEN indice = i: dmax = d
NEXT
' If max distance is greater than epsilon, recursively simplify
IF dmax > epsilon THEN
ListaDePuntos(indice, 3) = True
' Recursive call
CALL DRDP(ListaDePuntos(), ini, indice, epsilon)
CALL DRDP(ListaDePuntos(), indice, fin, epsilon)
END IF
END SUB
- Output:
Remaining points after simplifying: ( 0 , 0 ) ( 2 , -.1 ) ( 3 , 5 ) ( 7 , 9 ) ( 9 , 9 )
QB64
The QBasic solution works without any changes.
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
Version 1
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)
Version 2
/*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
- 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)
RPL
≪ 3 PICK - CONJ SWAP ROT - (0,1) * DUP ABS / * RE ABS ≫ ≫ 'DM→AB' 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 DM→AB 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 RDP points ROT nb SUB epsilon RDP 2 nb SUB + END END ≫ ≫ 'RDP' 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 RDP
- Output:
1: { (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)
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)
}
- 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)
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)]
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);
]
- 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