Ramer-Douglas-Peucker line simplification: Difference between revisions
Ramer-Douglas-Peucker line simplification en FreeBASIC |
m →{{header|Phix}}: syntax coloured |
||
Line 1,083: | Line 1,083: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<lang Phix> |
<!--<lang Phix>(phixonline)--> |
||
<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> |
|||
<!--</lang>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Revision as of 03:21, 28 March 2022
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
The Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.
- Task
Using the Ramer–Douglas–Peucker algorithm, simplify the 2D line defined by the points:
(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)
The error threshold to be used is: 1.0.
Display the remaining points here.
- Reference
-
- the Wikipedia article: Ramer-Douglas-Peucker algorithm.
11l
<lang 11l>F rdp(l, ε) -> [(Float, Float)]
V x = 0 V dMax = -1.0 V p1 = l[0] V p2 = l.last V p21 = p2 - p1
L(p) l[1.<(len)-1] V d = abs(cross(p, p21) + cross(p2, p1)) I d > dMax x = L.index + 1 dMax = d I dMax > ε R rdp(l[0..x], ε) [+] rdp(l[x..], ε)[1..]
R [l[0], l.last]
print(rdp([(0.0, 0.0),
(1.0, 0.1), (2.0,-0.1), (3.0, 5.0), (4.0, 6.0), (5.0, 7.0), (6.0, 8.1), (7.0, 9.0), (8.0, 9.0), (9.0, 9.0)], 1.0))</lang>
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
C
<lang c>#include <assert.h>
- include <math.h>
- include <stdio.h>
typedef struct point_tag {
double x; double y;
} point_t;
// Returns the distance from point p to the line between p1 and p2 double perpendicular_distance(point_t p, point_t p1, point_t p2) {
double dx = p2.x - p1.x; double dy = p2.y - p1.y; double d = sqrt(dx * dx + dy * dy); return fabs(p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x)/d;
}
// Simplify an array of points using the Ramer–Douglas–Peucker algorithm. // Returns the number of output points. size_t douglas_peucker(const point_t* points, size_t n, double epsilon,
point_t* dest, size_t destlen) { assert(n >= 2); assert(epsilon >= 0); double max_dist = 0; size_t index = 0; for (size_t i = 1; i + 1 < n; ++i) { double dist = perpendicular_distance(points[i], points[0], points[n - 1]); if (dist > max_dist) { max_dist = dist; index = i; } } if (max_dist > epsilon) { size_t n1 = douglas_peucker(points, index + 1, epsilon, dest, destlen); if (destlen >= n1 - 1) { destlen -= n1 - 1; dest += n1 - 1; } else { destlen = 0; } size_t n2 = douglas_peucker(points + index, n - index, epsilon, dest, destlen); return n1 + n2 - 1; } if (destlen >= 2) { dest[0] = points[0]; dest[1] = points[n - 1]; } return 2;
}
void print_points(const point_t* points, size_t n) {
for (size_t i = 0; i < n; ++i) { if (i > 0) printf(" "); printf("(%g, %g)", points[i].x, points[i].y); } printf("\n");
}
int main() {
point_t points[] = { {0,0}, {1,0.1}, {2,-0.1}, {3,5}, {4,6}, {5,7}, {6,8.1}, {7,9}, {8,9}, {9,9} }; const size_t len = sizeof(points)/sizeof(points[0]); point_t out[len]; size_t n = douglas_peucker(points, len, 1.0, out, len); print_points(out, n); return 0;
}</lang>
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
C#
<lang csharp>using System; using System.Collections.Generic; using System.Linq;
namespace LineSimplification {
using Point = Tuple<double, double>;
class Program { static double PerpendicularDistance(Point pt, Point lineStart, Point lineEnd) { double dx = lineEnd.Item1 - lineStart.Item1; double dy = lineEnd.Item2 - lineStart.Item2;
// Normalize double mag = Math.Sqrt(dx * dx + dy * dy); if (mag > 0.0) { dx /= mag; dy /= mag; } double pvx = pt.Item1 - lineStart.Item1; double pvy = pt.Item2 - lineStart.Item2;
// Get dot product (project pv onto normalized direction) double pvdot = dx * pvx + dy * pvy;
// Scale line direction vector and subtract it from pv double ax = pvx - pvdot * dx; double ay = pvy - pvdot * dy;
return Math.Sqrt(ax * ax + ay * ay); }
static void RamerDouglasPeucker(List<Point> pointList, double epsilon, List<Point> output) { if (pointList.Count < 2) { throw new ArgumentOutOfRangeException("Not enough points to simplify"); }
// Find the point with the maximum distance from line between the start and end double dmax = 0.0; int index = 0; int end = pointList.Count - 1; for (int i = 1; i < end; ++i) { double d = PerpendicularDistance(pointList[i], pointList[0], pointList[end]); if (d > dmax) { index = i; dmax = d; } }
// If max distance is greater than epsilon, recursively simplify if (dmax > epsilon) { List<Point> recResults1 = new List<Point>(); List<Point> recResults2 = new List<Point>(); List<Point> firstLine = pointList.Take(index + 1).ToList(); List<Point> lastLine = pointList.Skip(index).ToList(); RamerDouglasPeucker(firstLine, epsilon, recResults1); RamerDouglasPeucker(lastLine, epsilon, recResults2);
// build the result list output.AddRange(recResults1.Take(recResults1.Count - 1)); output.AddRange(recResults2); if (output.Count < 2) throw new Exception("Problem assembling output"); } else { // Just return start and end points output.Clear(); output.Add(pointList[0]); output.Add(pointList[pointList.Count - 1]); } }
static void Main(string[] args) { List<Point> pointList = new List<Point>() { new Point(0.0,0.0), new Point(1.0,0.1), new Point(2.0,-0.1), new Point(3.0,5.0), new Point(4.0,6.0), new Point(5.0,7.0), new Point(6.0,8.1), new Point(7.0,9.0), new Point(8.0,9.0), new Point(9.0,9.0), }; List<Point> pointListOut = new List<Point>(); RamerDouglasPeucker(pointList, 1.0, pointListOut); Console.WriteLine("Points remaining after simplification:"); pointListOut.ForEach(p => Console.WriteLine(p)); } }
}</lang>
- Output:
Points remaining after simplification: (0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
C++
<lang cpp>#include <iostream>
- include <cmath>
- include <utility>
- include <vector>
- include <stdexcept>
using namespace std;
typedef std::pair<double, double> Point;
double PerpendicularDistance(const Point &pt, const Point &lineStart, const Point &lineEnd) { double dx = lineEnd.first - lineStart.first; double dy = lineEnd.second - lineStart.second;
//Normalise double mag = pow(pow(dx,2.0)+pow(dy,2.0),0.5); if(mag > 0.0) { dx /= mag; dy /= mag; }
double pvx = pt.first - lineStart.first; double pvy = pt.second - lineStart.second;
//Get dot product (project pv onto normalized direction) double pvdot = dx * pvx + dy * pvy;
//Scale line direction vector double dsx = pvdot * dx; double dsy = pvdot * dy;
//Subtract this from pv double ax = pvx - dsx; double ay = pvy - dsy;
return pow(pow(ax,2.0)+pow(ay,2.0),0.5); }
void RamerDouglasPeucker(const vector<Point> &pointList, double epsilon, vector<Point> &out) { if(pointList.size()<2) throw invalid_argument("Not enough points to simplify");
// Find the point with the maximum distance from line between start and end double dmax = 0.0; size_t index = 0; size_t end = pointList.size()-1; for(size_t i = 1; i < end; i++) { double d = PerpendicularDistance(pointList[i], pointList[0], pointList[end]); if (d > dmax) { index = i; dmax = d; } }
// If max distance is greater than epsilon, recursively simplify if(dmax > epsilon) { // Recursive call vector<Point> recResults1; vector<Point> recResults2; vector<Point> firstLine(pointList.begin(), pointList.begin()+index+1); vector<Point> lastLine(pointList.begin()+index, pointList.end()); RamerDouglasPeucker(firstLine, epsilon, recResults1); RamerDouglasPeucker(lastLine, epsilon, recResults2);
// Build the result list out.assign(recResults1.begin(), recResults1.end()-1); out.insert(out.end(), recResults2.begin(), recResults2.end()); if(out.size()<2) throw runtime_error("Problem assembling output"); } else { //Just return start and end points out.clear(); out.push_back(pointList[0]); out.push_back(pointList[end]); } }
int main() { vector<Point> pointList; vector<Point> pointListOut;
pointList.push_back(Point(0.0, 0.0)); pointList.push_back(Point(1.0, 0.1)); pointList.push_back(Point(2.0, -0.1)); pointList.push_back(Point(3.0, 5.0)); pointList.push_back(Point(4.0, 6.0)); pointList.push_back(Point(5.0, 7.0)); pointList.push_back(Point(6.0, 8.1)); pointList.push_back(Point(7.0, 9.0)); pointList.push_back(Point(8.0, 9.0)); pointList.push_back(Point(9.0, 9.0));
RamerDouglasPeucker(pointList, 1.0, pointListOut);
cout << "result" << endl; for(size_t i=0;i< pointListOut.size();i++) { cout << pointListOut[i].first << "," << pointListOut[i].second << endl; }
return 0; }</lang>
- Output:
result 0,0 2,-0.1 3,5 7,9 9,9
D
<lang D>import std.algorithm; import std.exception : enforce; import std.math; import std.stdio;
void main() {
creal[] pointList = [ 0.0 + 0.0i, 1.0 + 0.1i, 2.0 + -0.1i, 3.0 + 5.0i, 4.0 + 6.0i, 5.0 + 7.0i, 6.0 + 8.1i, 7.0 + 9.0i, 8.0 + 9.0i, 9.0 + 9.0i ]; creal[] pointListOut;
ramerDouglasPeucker(pointList, 1.0, pointListOut);
writeln("result"); for (size_t i=0; i< pointListOut.length; i++) { writeln(pointListOut[i].re, ",", pointListOut[i].im); }
}
real perpendicularDistance(const creal pt, const creal lineStart, const creal lineEnd) {
creal d = lineEnd - lineStart;
//Normalise real mag = hypot(d.re, d.im); if (mag > 0.0) { d /= mag; }
creal pv = pt - lineStart;
//Get dot product (project pv onto normalized direction) real pvdot = d.re * pv.re + d.im * pv.im;
//Scale line direction vector creal ds = pvdot * d;
//Subtract this from pv creal a = pv - ds;
return hypot(a.re, a.im);
}
void ramerDouglasPeucker(const creal[] pointList, real epsilon, ref creal[] output) {
enforce(pointList.length >= 2, "Not enough points to simplify");
// Find the point with the maximum distance from line between start and end real dmax = 0.0; size_t index = 0; size_t end = pointList.length-1; for (size_t i=1; i<end; i++) { real d = perpendicularDistance(pointList[i], pointList[0], pointList[end]); if (d > dmax) { index = i; dmax = d; } }
// If max distance is greater than epsilon, recursively simplify if (dmax > epsilon) { // Recursive call creal[] firstLine = pointList[0..index+1].dup; creal[] lastLine = pointList[index+1..$].dup;
creal[] recResults1; ramerDouglasPeucker(firstLine, epsilon, recResults1);
creal[] recResults2; ramerDouglasPeucker(lastLine, epsilon, recResults2);
// Build the result list output = recResults1 ~ recResults2;
enforce(output.length>=2, "Problem assembling output"); } else { //Just return start and end points output.length = 0; output ~= pointList[0]; output ~= pointList[end]; }
}</lang>
- Output:
result 0,0 2,-0.1 3,5 7,9 9,9
FreeBASIC
<lang freebasic> Function DistanciaPerpendicular(tabla() As Double, i As Integer, ini As Integer, fin As Integer) As Double
Dim As Double dx, dy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay dx = tabla(fin, 1) - tabla(ini, 1) dy = tabla(fin, 2) - tabla(ini, 2) 'Normaliza mag = (dx^2 + dy^2)^0.5 If mag > 0 Then dx /= mag : dy /= mag pvx = tabla(i, 1) - tabla(ini, 1) pvy = tabla(i, 2) - tabla(ini, 2) 'Producto escalado (proyecto pv en dirección normalizada) pvdot = dx * pvx + dy * pvy 'Vector de dirección de línea de escala dsx = pvdot * dx dsy = pvdot * dy 'Reste esto de pv ax = pvx - dsx ay = pvy - dsy Return (ax^2 + ay^2)^0.5
End Function
Sub DRDP(ListaDePuntos() As Double, ini As Integer, fin As Integer, epsilon As Double)
Dim As Double dmax, d Dim As Integer indice, i ' Encuentra el punto con la distancia máxima If Ubound(ListaDePuntos) < 2 Then Print "No hay suficientes puntos para simplificar " : Sleep : End For i = ini + 1 To fin d = DistanciaPerpendicular(ListaDePuntos(), i, ini, fin) If d > dmax Then indice = i : dmax = d Next ' Si la distancia máxima es mayor que épsilon, simplifique de forma recursiva If dmax > epsilon Then ListaDePuntos(indice, 3) = True DRDP(ListaDePuntos(), ini, indice, epsilon) DRDP(ListaDePuntos(), indice, fin, epsilon) End If
End Sub
Dim As Double matriz(1 To 10, 1 To 3) Data 0, 0, 1, 0.1, 2, -0.1, 3, 5, 4, 6, 5, 7, 6, 8.1, 7, 9, 8, 9, 9, 9 For i As Integer = Lbound(matriz) To Ubound(matriz)
Read matriz(i, 1), matriz(i, 2)
Next i
DRDP(matriz(), 1, 10, 1)
Print "Puntos restantes tras de simplificar:" matriz(1, 3) = True : matriz(10, 3) = True For i As Integer = Lbound(matriz) To Ubound(matriz)
If matriz(i, 3) Then Print "(";matriz(i, 1);", "; matriz(i, 2);") ";
Next i Sleep </lang>
- Output:
Puntos restantes tras de simplificar: ( 0, 0) ( 2, -0.1) ( 3, 5) ( 7, 9) ( 9, 9)
Go
<lang go>package main
import (
"fmt" "math"
)
type point struct{ x, y float64 }
func RDP(l []point, ε float64) []point {
x := 0 dMax := -1. last := len(l) - 1 p1 := l[0] p2 := l[last] x21 := p2.x - p1.x y21 := p2.y - p1.y for i, p := range l[1:last] { if d := math.Abs(y21*p.x - x21*p.y + p2.x*p1.y - p2.y*p1.x); d > dMax { x = i + 1 dMax = d } } if dMax > ε { return append(RDP(l[:x+1], ε), RDP(l[x:], ε)[1:]...) } return []point{l[0], l[len(l)-1]}
}
func main() {
fmt.Println(RDP([]point{{0, 0}, {1, 0.1}, {2, -0.1}, {3, 5}, {4, 6}, {5, 7}, {6, 8.1}, {7, 9}, {8, 9}, {9, 9}}, 1))
}</lang>
- Output:
[{0 0} {2 -0.1} {3 5} {7 9} {9 9}]
J
Solution: <lang j>mp=: +/ .* NB. matrix product norm=: +/&.:*: NB. vector norm normalize=: (% norm)^:(0 < norm)
dxy=. normalize@({: - {.) pv=. -"1 {. NB.*perpDist v Calculate perpendicular distance of points from a line perpDist=: norm"1@(pv ([ -"1 mp"1~ */ ]) dxy) f.
rdp=: verb define
1 rdp y : points=. ,:^:(2 > #@$) y epsilon=. x if. 2 > # points do. points return. end.
NB. point with the maximum distance from line between start and end 'imax dmax'=. ((i. , ]) >./) perpDist points if. dmax > epsilon do. epsilon ((}:@rdp (1+imax)&{.) , (rdp imax&}.)) points else. ({. ,: {:) points end.
)</lang> Example Usage: <lang j> Points=: 0 0,1 0.1,2 _0.1,3 5,4 6,5 7,6 8.1,7 9,8 9,:9 9
1.0 rdp Points
0 0 2 _0.1 3 5 7 9 9 9</lang>
Java
<lang Java>import javafx.util.Pair;
import java.util.ArrayList; import java.util.List;
public class LineSimplification {
private static class Point extends Pair<Double, Double> { Point(Double key, Double value) { super(key, value); }
@Override public String toString() { return String.format("(%f, %f)", getKey(), getValue()); } }
private static double perpendicularDistance(Point pt, Point lineStart, Point lineEnd) { double dx = lineEnd.getKey() - lineStart.getKey(); double dy = lineEnd.getValue() - lineStart.getValue();
// Normalize double mag = Math.hypot(dx, dy); if (mag > 0.0) { dx /= mag; dy /= mag; } double pvx = pt.getKey() - lineStart.getKey(); double pvy = pt.getValue() - lineStart.getValue();
// Get dot product (project pv onto normalized direction) double pvdot = dx * pvx + dy * pvy;
// Scale line direction vector and subtract it from pv double ax = pvx - pvdot * dx; double ay = pvy - pvdot * dy;
return Math.hypot(ax, ay); }
private static void ramerDouglasPeucker(List<Point> pointList, double epsilon, List<Point> out) { if (pointList.size() < 2) throw new IllegalArgumentException("Not enough points to simplify");
// Find the point with the maximum distance from line between the start and end double dmax = 0.0; int index = 0; int end = pointList.size() - 1; for (int i = 1; i < end; ++i) { double d = perpendicularDistance(pointList.get(i), pointList.get(0), pointList.get(end)); if (d > dmax) { index = i; dmax = d; } }
// If max distance is greater than epsilon, recursively simplify if (dmax > epsilon) { List<Point> recResults1 = new ArrayList<>(); List<Point> recResults2 = new ArrayList<>(); List<Point> firstLine = pointList.subList(0, index + 1); List<Point> lastLine = pointList.subList(index, pointList.size()); ramerDouglasPeucker(firstLine, epsilon, recResults1); ramerDouglasPeucker(lastLine, epsilon, recResults2);
// build the result list out.addAll(recResults1.subList(0, recResults1.size() - 1)); out.addAll(recResults2); if (out.size() < 2) throw new RuntimeException("Problem assembling output"); } else { // Just return start and end points out.clear(); out.add(pointList.get(0)); out.add(pointList.get(pointList.size() - 1)); } }
public static void main(String[] args) { List<Point> pointList = List.of( new Point(0.0, 0.0), new Point(1.0, 0.1), new Point(2.0, -0.1), new Point(3.0, 5.0), new Point(4.0, 6.0), new Point(5.0, 7.0), new Point(6.0, 8.1), new Point(7.0, 9.0), new Point(8.0, 9.0), new Point(9.0, 9.0) ); List<Point> pointListOut = new ArrayList<>(); ramerDouglasPeucker(pointList, 1.0, pointListOut); System.out.println("Points remaining after simplification:"); pointListOut.forEach(System.out::println); }
}</lang>
- Output:
Points remaining after simplification: (0.000000, 0.000000) (2.000000, -0.100000) (3.000000, 5.000000) (7.000000, 9.000000) (9.000000, 9.000000)
JavaScript
<lang JavaScript>/**
* @typedef {{ * x: (!number), * y: (!number) * }} */
let pointType;
/**
* @param {!Array<pointType>} l * @param {number} eps */
const RDP = (l, eps) => {
const last = l.length - 1; const p1 = l[0]; const p2 = l[last]; const x21 = p2.x - p1.x; const y21 = p2.y - p1.y;
const [dMax, x] = l.slice(1, last) .map(p => Math.abs(y21 * p.x - x21 * p.y + p2.x * p1.y - p2.y * p1.x)) .reduce((p, c, i) => { const v = Math.max(p[0], c); return [v, v === p[0] ? p[1] : i + 1]; }, [-1, 0]);
if (dMax > eps) { return [...RDP(l.slice(0, x + 1), eps), ...RDP(l.slice(x), eps).slice(1)]; } return [l[0], l[last]]
};
const points = [
{x: 0, y: 0}, {x: 1, y: 0.1}, {x: 2, y: -0.1}, {x: 3, y: 5}, {x: 4, y: 6}, {x: 5, y: 7}, {x: 6, y: 8.1}, {x: 7, y: 9}, {x: 8, y: 9}, {x: 9, y: 9}];
console.log(RDP(points, 1));</lang>
- Output:
[{x: 0, y: 0}, {x: 2, y: -0.1}, {x: 3, y: 5}, {x: 7, y: 9}, {x: 9, y: 9}]
Julia
<lang julia>const Point = Vector{Float64}
function perpdist(pt::Point, lnstart::Point, lnend::Point)
d = normalize!(lnend .- lnstart)
pv = pt .- lnstart # Get dot product (project pv onto normalized direction) pvdot = dot(d, pv) # Scale line direction vector ds = pvdot .* d # Subtract this from pv return norm(pv .- ds)
end
function rdp(plist::Vector{Point}, ϵ::Float64 = 1.0)
if length(plist) < 2 throw(ArgumentError("not enough points to simplify")) end
# Find the point with the maximum distance from line between start and end distances = collect(perpdist(pt, plist[1], plist[end]) for pt in plist) dmax, imax = findmax(distances)
# If max distance is greater than epsilon, recursively simplify if dmax > ϵ fstline = plist[1:imax] lstline = plist[imax:end]
recrst1 = rdp(fstline, ϵ) recrst2 = rdp(lstline, ϵ)
out = vcat(recrst1, recrst2) else out = [plist[1], plist[end]] end
return out
end
plist = Point[[0.0, 0.0], [1.0, 0.1], [2.0, -0.1], [3.0, 5.0], [4.0, 6.0], [5.0, 7.0], [6.0, 8.1], [7.0, 9.0], [8.0, 9.0], [9.0, 9.0]] @show plist @show rdp(plist)</lang>
- Output:
plist = Array{Float64,1}[[0.0, 0.0], [1.0, 0.1], [2.0, -0.1], [3.0, 5.0], [4.0, 6.0], [5.0, 7.0], [6.0, 8.1], [7.0, 9.0], [8.0, 9.0], [9.0, 9.0]] rdp(plist) = Array{Float64,1}[[0.0, 0.0], [2.0, -0.1], [2.0, -0.1], [3.0, 5.0], [3.0, 5.0], [7.0, 9.0], [7.0, 9.0], [9.0, 9.0]]
Kotlin
<lang scala>// version 1.1.0
typealias Point = Pair<Double, Double>
fun perpendicularDistance(pt: Point, lineStart: Point, lineEnd: Point): Double {
var dx = lineEnd.first - lineStart.first var dy = lineEnd.second - lineStart.second
// Normalize val mag = Math.hypot(dx, dy) if (mag > 0.0) { dx /= mag; dy /= mag } val pvx = pt.first - lineStart.first val pvy = pt.second - lineStart.second
// Get dot product (project pv onto normalized direction) val pvdot = dx * pvx + dy * pvy // Scale line direction vector and substract it from pv val ax = pvx - pvdot * dx val ay = pvy - pvdot * dy return Math.hypot(ax, ay)
}
fun RamerDouglasPeucker(pointList: List<Point>, epsilon: Double, out: MutableList<Point>) {
if (pointList.size < 2) throw IllegalArgumentException("Not enough points to simplify") // Find the point with the maximum distance from line between start and end var dmax = 0.0 var index = 0 val end = pointList.size - 1 for (i in 1 until end) { val d = perpendicularDistance(pointList[i], pointList[0], pointList[end]) if (d > dmax) { index = i; dmax = d } }
// If max distance is greater than epsilon, recursively simplify if (dmax > epsilon) { val recResults1 = mutableListOf<Point>() val recResults2 = mutableListOf<Point>() val firstLine = pointList.take(index + 1) val lastLine = pointList.drop(index) RamerDouglasPeucker(firstLine, epsilon, recResults1) RamerDouglasPeucker(lastLine, epsilon, recResults2)
// build the result list out.addAll(recResults1.take(recResults1.size - 1)) out.addAll(recResults2) if (out.size < 2) throw RuntimeException("Problem assembling output") } else { // Just return start and end points out.clear() out.add(pointList.first()) out.add(pointList.last()) }
}
fun main(args: Array<String>) {
val pointList = listOf( Point(0.0, 0.0), Point(1.0, 0.1), Point(2.0, -0.1), Point(3.0, 5.0), Point(4.0, 6.0), Point(5.0, 7.0), Point(6.0, 8.1),
Point(7.0, 9.0), Point(8.0, 9.0),
Point(9.0, 9.0) ) val pointListOut = mutableListOf<Point>() RamerDouglasPeucker(pointList, 1.0, pointListOut) println("Points remaining after simplification:") for (p in pointListOut) println(p)
}</lang>
- Output:
Points remaining after simplification: (0.0, 0.0) (2.0, -0.1) (3.0, 5.0) (7.0, 9.0) (9.0, 9.0)
Nim
<lang nim>import math
type
Point = tuple[x, y: float64]
proc pointLineDistance(pt, lineStart, lineEnd: Point): float64 =
var n, d, dx, dy: float64 dx = lineEnd.x - lineStart.x dy = lineEnd.y - lineStart.y n = abs(dx * (lineStart.y - pt.y) - (lineStart.x - pt.x) * dy) d = sqrt(dx * dx + dy * dy) n / d
proc rdp(points: seq[Point], startIndex, lastIndex: int, ε: float64 = 1.0): seq[Point] =
var dmax = 0.0 var index = startIndex
for i in index+1..<lastIndex: var d = pointLineDistance(points[i], points[startIndex], points[lastIndex]) if d > dmax: index = i dmax = d if dmax > ε: var res1 = rdp(points, startIndex, index, ε) var res2 = rdp(points, index, lastIndex, ε)
var finalRes: seq[Point] = @[] finalRes.add(res1[0..^2]) finalRes.add(res2[0..^1]) result = finalRes else: result = @[points[startIndex], points[lastIndex]]
var line: seq[Point] = @[(0.0, 0.0), (1.0, 0.1), (2.0, -0.1), (3.0, 5.0), (4.0, 6.0),
(5.0, 7.0), (6.0, 8.1), (7.0, 9.0), (8.0, 9.0), (9.0, 9.0)]
echo rdp(line, line.low, line.high)</lang>
- Output:
@[(x: 0.0, y: 0.0), (x: 2.0, y: -0.1), (x: 3.0, y: 5.0), (x: 7.0, y: 9.0), (x: 9.0, y: 9.0)]
Openscad
<lang openscad>function slice(a, v) = [ for (i = v) a[i] ];
// Find the distance from the point to the line function perpendicular_distance(pt, start, end) =
let ( d = end - start, dn = d / norm(d), dp = pt - start, // Dot-product of two vectors ddot = dn * dp, ds = dn * ddot ) norm(dp - ds);
// Find the point with the maximum distance from line between start and end. // Returns a vector of [dmax, point_index] function max_distance(pts, cindex=1, iresult=[0,0]) =
let ( d = perpendicular_distance(pts[cindex], pts[0],pts[len(pts)-1]), result = (d > iresult[0] ? [d, cindex] : iresult) ) (cindex == len(pts) - 1 ? iresult : max_distance(pts, cindex+1, result));
// Tail recursion optimization is not possible with tree recursion. // // It's probably possible to optimize this with tail recursion by converting to // iterative approach, see // https://namekdev.net/2014/06/iterative-version-of-ramer-douglas-peucker-line-simplification-algorithm/ . // It's not possible to use iterative approach without recursion at all because // of static nature of OpenSCAD arrays. function ramer_douglas_peucker(points, epsilon=1) =
let (dmax = max_distance(points)) // Simplify the line recursively if the max distance is greater than epsilon (dmax[0] > epsilon ? let ( lo = ramer_douglas_peucker( slice(points, [0:dmax[1]]), epsilon), hi = ramer_douglas_peucker( slice(points, [dmax[1]:len(points)-1]), epsilon) ) concat(slice(lo, [0:len(lo)-2]), hi) : [points[0],points[len(points)-1]]);
/* -- Demo -- */
module line(points, thickness=1, dot=2) {
for (idx = [0:len(points)-2]) { translate(points[idx]) sphere(d=dot); hull() { translate(points[idx]) sphere(d=thickness); translate(points[idx+1]) sphere(d=thickness); } } translate(points[len(points)-1]) sphere(d=dot);
}
function addz(pts, z=0) = [ for (p = pts) [p.x, p.y, z] ];
module demo(pts) {
rdp = ramer_douglas_peucker(points=pts, epsilon=1); echo(rdp); color("gray") line(points=addz(pts, 0), thickness=0.1, dot=0.3); color("red") line(points=addz(rdp, 1), thickness=0.1, dot=0.3);
}
$fn=16; points = [[0,0], [1,0.1], [2,-0.1], [3,5], [4,6], [5,7], [6,8.1], [7,9], [8,9], [9,9]]; demo(points);</lang>
- Output:
ECHO: [[0, 0], [2, -0.1], [3, 5], [7, 9], [9, 9]]
Perl
<lang perl>use strict; use warnings; use feature 'say'; use List::MoreUtils qw(firstidx minmax);
my $epsilon = 1;
sub norm {
my(@list) = @_; my $sum; $sum += $_**2 for @list; sqrt($sum)
}
sub perpendicular_distance {
our(@start,@end,@point); local(*start,*end,*point) = (shift, shift, shift); return 0 if $start[0]==$point[0] && $start[1]==$point[1] or $end[0]==$point[0] && $end[1]==$point[1]; my ( $dx, $dy) = ( $end[0]-$start[0], $end[1]-$start[1]); my ($dpx, $dpy) = ($point[0]-$start[0],$point[1]-$start[1]); my $t = norm($dx, $dy); $dx /= $t; $dy /= $t; norm($dpx - $dx*($dx*$dpx + $dy*$dpy), $dpy - $dy*($dx*$dpx + $dy*$dpy));
}
sub Ramer_Douglas_Peucker {
my(@points) = @_; return @points if @points == 2; my @d; push @d, perpendicular_distance(@points[0, -1, $_]) for 0..@points-1; my(undef,$dmax) = minmax @d; my $index = firstidx { $_ == $dmax } @d; if ($dmax > $epsilon) { my @lo = Ramer_Douglas_Peucker( @points[0..$index]); my @hi = Ramer_Douglas_Peucker( @points[$index..$#points]); return @lo[0..@lo-2], @hi; } @points[0, -1];
}
say '(' . join(' ', @$_) . ') '
for Ramer_Douglas_Peucker( [0,0],[1,0.1],[2,-0.1],[3,5],[4,6],[5,7],[6,8.1],[7,9],[8,9],[9,9] )</lang>
- Output:
(0 0) (2 -0.1) (3 5) (7 9) (9 9)
Phix
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
<lang php>function perpendicular_distance(array $pt, array $line) {
// Calculate the normalized delta x and y of the line. $dx = $line[1][0] - $line[0][0]; $dy = $line[1][1] - $line[0][1]; $mag = sqrt($dx * $dx + $dy * $dy); if ($mag > 0) { $dx /= $mag; $dy /= $mag; }
// Calculate dot product, projecting onto normalized direction. $pvx = $pt[0] - $line[0][0]; $pvy = $pt[1] - $line[0][1]; $pvdot = $dx * $pvx + $dy * $pvy;
// Scale line direction vector and subtract from pv. $dsx = $pvdot * $dx; $dsy = $pvdot * $dy; $ax = $pvx - $dsx; $ay = $pvy - $dsy;
return sqrt($ax * $ax + $ay * $ay);
}
function ramer_douglas_peucker(array $points, $epsilon) {
if (count($points) < 2) { throw new InvalidArgumentException('Not enough points to simplify'); }
// Find the point with the maximum distance from the line between start/end. $dmax = 0; $index = 0; $end = count($points) - 1; $start_end_line = [$points[0], $points[$end]]; for ($i = 1; $i < $end; $i++) { $dist = perpendicular_distance($points[$i], $start_end_line); if ($dist > $dmax) { $index = $i; $dmax = $dist; } }
// If max distance is larger than epsilon, recursively simplify. if ($dmax > $epsilon) { $new_start = ramer_douglas_peucker(array_slice($points, 0, $index + 1), $epsilon); $new_end = ramer_douglas_peucker(array_slice($points, $index), $epsilon); array_pop($new_start); return array_merge($new_start, $new_end); }
// Max distance is below epsilon, so return a line from with just the // start and end points. return [ $points[0], $points[$end]];
}
$polyline = [
[0,0], [1,0.1], [2,-0.1], [3,5], [4,6], [5,7], [6,8.1], [7,9], [8,9], [9,9],
];
$result = ramer_douglas_peucker($polyline, 1.0); print "Result:\n"; foreach ($result as $point) {
print $point[0] . ',' . $point[1] . "\n";
}</lang>
- Output:
Result: 0,0 2,-0.1 3,5 7,9 9,9
Python
An approach using the shapely library: <lang python>from __future__ import print_function from shapely.geometry import LineString
if __name__=="__main__": line = LineString([(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]) print (line.simplify(1.0, preserve_topology=False))</lang>
- Output:
LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)
Racket
<lang racket>#lang racket (require math/flonum)
- points are lists of x y (maybe extensible to z)
- x+y gets both parts as values
(define (x+y p) (values (first p) (second p)))
(define (⊥-distance P1 P2)
(let*-values ([(x1 y1) (x+y P1)] [(x2 y2) (x+y P2)] [(dx dy) (values (- x2 x1) (- y2 y1))] [(h) (sqrt (+ (sqr dy) (sqr dx)))]) (λ (P0) (let-values (((x0 y0) (x+y P0))) (/ (abs (+ (* dy x0) (* -1 dx y0) (* x2 y1) (* -1 y2 x1))) h)))))
(define (douglas-peucker points-in ϵ)
(let recur ((ps points-in)) ;; curried distance function which will be applicable to all points (let*-values ([(p0) (first ps)] [(pz) (last ps)] [(p-d) (⊥-distance p0 pz)] ;; Find the point with the maximum distance [(dmax index) (for/fold ((dmax 0) (index 0)) ((i (in-range 1 (sub1 (length ps))))) ; skips the first, stops before the last (define d (p-d (list-ref ps i))) (if (> d dmax) (values d i) (values dmax index)))]) ;; If max distance is greater than epsilon, recursively simplify (if (> dmax ϵ) ;; recursive call (let-values ([(l r) (split-at ps index)]) (append (drop-right (recur l) 1) (recur r))) (list p0 pz))))) ;; else we can return this simplification
(module+ main
(douglas-peucker '((0 0) (1 0.1) (2 -0.1) (3 5) (4 6) (5 7) (6 8.1) (7 9) (8 9) (9 9)) 1.0))
(module+ test
(require rackunit) (check-= ((⊥-distance '(0 0) '(0 1)) '(1 0)) 1 epsilon.0))</lang>
- Output:
'((0 0) (2 -0.1) (3 5) (7 9) (9 9))
Raku
(formerly Perl 6)
<lang perl6>sub norm (*@list) { @list»².sum.sqrt }
sub perpendicular-distance (@start, @end where @end !eqv @start, @point) {
return 0 if @point eqv any(@start, @end); my ( $Δx, $Δy ) = @end «-» @start; my ($Δpx, $Δpy) = @point «-» @start; ($Δx, $Δy) «/=» norm $Δx, $Δy; norm ($Δpx, $Δpy) «-» ($Δx, $Δy) «*» ($Δx*$Δpx + $Δy*$Δpy);
}
sub Ramer-Douglas-Peucker(@points where * > 1, \ε = 1) {
return @points if @points == 2; my @d = (^@points).map: { perpendicular-distance |@points[0, *-1, $_] }; my ($index, $dmax) = @d.first: @d.max, :kv; return flat Ramer-Douglas-Peucker( @points[0..$index], ε )[^(*-1)], Ramer-Douglas-Peucker( @points[$index..*], ε ) if $dmax > ε; @points[0, *-1];
}
- TESTING
say Ramer-Douglas-Peucker(
[(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]
);</lang>
- Output:
((0 0) (2 -0.1) (3 5) (7 9) (9 9))
REXX
The computation for the perpendicular distance was taken from the GO example. <lang rexx>/*REXX program uses the Ramer─Douglas─Peucker (RDP) line simplification algorithm for*/ /*───────────────────────────── reducing the number of points used to define its shape. */ parse arg epsilon pts /*obtain optional arguments from the CL*/ if epsilon= | epsilon="," then epsilon= 1 /*Not specified? Then use the default.*/ if pts= then pts= '(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)' pts= space(pts) /*elide all superfluous blanks. */
say ' error threshold: ' epsilon /*echo the error threshold to the term.*/ say ' points specified: ' pts /* " " shape points " " " */
$= RDP(pts) /*invoke Ramer─Douglas─Peucker function*/
say 'points simplified: ' rez($) /*display points with () ───► terminal.*/
exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ bld: parse arg _; #= words(_); dMax=-#; idx=1; do j=1 for #; @.j= word(_, j); end; return px: parse arg _; return word( translate(_, , ','), 1) /*obtain the X coörd.*/ py: parse arg _; return word( translate(_, , ','), 2) /* " " Y " */ reb: parse arg a,b,,_; do k=a to b; _= _ @.k; end; return strip(_) rez: parse arg z,_; do k=1 for words(z); _= _ '('word(z, k)") "; end; return strip(_) /*──────────────────────────────────────────────────────────────────────────────────────*/ RDP: procedure expose epsilon; call bld space( translate(arg(1), , ')(][}{') )
L= px(@.#) - px(@.1) H= py(@.#) - py(@.1) /* [↓] find point IDX with max distance*/ do i=2 to #-1 d= abs(H*px(@.i) - L*py(@.i) + px(@.#)*py(@.1) - py(@.#)*px(@.1)) if d>dMax then do; idx= i; dMax= d end end /*i*/ /* [↑] D is the perpendicular distance*/
if dMax>epsilon then do; r= RDP( reb(1, idx) ) return subword(r, 1, words(r) - 1) RDP( reb(idx, #) ) end return @.1 @.#</lang>
- output when using the default inputs:
error threshold: 1 points specified: (0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9) points simplified: (0,0) (2,-0.1) (3,5) (7,9) (9,9)
Rust
An implementation of the algorithm
<lang rust>#[derive(Copy, Clone)] struct Point {
x: f64, y: f64,
}
use std::fmt;
impl fmt::Display for Point {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { write!(f, "({}, {})", self.x, self.y) }
}
// Returns the distance from point p to the line between p1 and p2 fn perpendicular_distance(p: &Point, p1: &Point, p2: &Point) -> f64 {
let dx = p2.x - p1.x; let dy = p2.y - p1.y; (p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x).abs() / dx.hypot(dy)
}
fn rdp(points: &[Point], epsilon: f64, result: &mut Vec<Point>) {
let n = points.len(); if n < 2 { return; } let mut max_dist = 0.0; let mut index = 0; for i in 1..n - 1 { let dist = perpendicular_distance(&points[i], &points[0], &points[n - 1]); if dist > max_dist { max_dist = dist; index = i; } } if max_dist > epsilon { rdp(&points[0..=index], epsilon, result); rdp(&points[index..n], epsilon, result); } else { result.push(points[n - 1]); }
}
fn ramer_douglas_peucker(points: &[Point], epsilon: f64) -> Vec<Point> {
let mut result = Vec::new(); if points.len() > 0 && epsilon >= 0.0 { result.push(points[0]); rdp(points, epsilon, &mut result); } result
}
fn main() {
let points = vec![ Point { x: 0.0, y: 0.0 }, Point { x: 1.0, y: 0.1 }, Point { x: 2.0, y: -0.1 }, Point { x: 3.0, y: 5.0 }, Point { x: 4.0, y: 6.0 }, Point { x: 5.0, y: 7.0 }, Point { x: 6.0, y: 8.1 }, Point { x: 7.0, y: 9.0 }, Point { x: 8.0, y: 9.0 }, Point { x: 9.0, y: 9.0 }, ]; for p in ramer_douglas_peucker(&points, 1.0) { println!("{}", p); }
}</lang>
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
Using the geo crate
The geo crate contains an implementation of the Ramer-Douglas-Peucker line simplification algorithm. <lang rust>// [dependencies] // geo = "0.14"
use geo::algorithm::simplify::Simplify; use geo::line_string;
fn main() {
let points = line_string![ (x: 0.0, y: 0.0), (x: 1.0, y: 0.1), (x: 2.0, y: -0.1), (x: 3.0, y: 5.0), (x: 4.0, y: 6.0), (x: 5.0, y: 7.0), (x: 6.0, y: 8.1), (x: 7.0, y: 9.0), (x: 8.0, y: 9.0), (x: 9.0, y: 9.0), ]; for p in points.simplify(&1.0) { println!("({}, {})", p.x, p.y); }
}</lang>
- Output:
(0, 0) (2, -0.1) (3, 5) (7, 9) (9, 9)
Sidef
<lang ruby>func perpendicular_distance(Arr start, Arr end, Arr point) {
((point == start) || (point == end)) && return 0 var (Δx, Δy ) = ( end »-« start)... var (Δpx, Δpy) = (point »-« start)... var h = hypot(Δx, Δy) [\Δx, \Δy].map { *_ /= h } (([Δpx, Δpy] »-« ([Δx, Δy] »*» (Δx*Δpx + Δy*Δpy))) »**» 2).sum.sqrt
}
func Ramer_Douglas_Peucker(Arr points { .all { .len > 1 } }, ε = 1) {
points.len == 2 && return points
var d = (^points -> map { perpendicular_distance(points[0], points[-1], points[_]) })
if (d.max > ε) { var i = d.index(d.max) return [Ramer_Douglas_Peucker(points.ft(0, i), ε).ft(0, -2)..., Ramer_Douglas_Peucker(points.ft(i), ε)...] }
return [points[0,-1]]
}
say Ramer_Douglas_Peucker(
[[0,0],[1,0.1],[2,-0.1],[3,5],[4,6],[5,7],[6,8.1],[7,9],[8,9],[9,9]]
)</lang>
- Output:
[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]]
Swift
<lang swift>struct Point: CustomStringConvertible {
let x: Double, y: Double
var description: String { return "(\(x), \(y))" }
}
// Returns the distance from point p to the line between p1 and p2 func perpendicularDistance(p: Point, p1: Point, p2: Point) -> Double {
let dx = p2.x - p1.x let dy = p2.y - p1.y let d = (p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x) return abs(d)/(dx * dx + dy * dy).squareRoot()
}
func ramerDouglasPeucker(points: [Point], epsilon: Double) -> [Point] {
var result : [Point] = [] func rdp(begin: Int, end: Int) { guard end > begin else { return } var maxDist = 0.0 var index = 0 for i in begin+1..<end { let dist = perpendicularDistance(p: points[i], p1: points[begin], p2: points[end]) if dist > maxDist { maxDist = dist index = i } } if maxDist > epsilon { rdp(begin: begin, end: index) rdp(begin: index, end: end) } else { result.append(points[end]) } } if points.count > 0 && epsilon >= 0.0 { result.append(points[0]) rdp(begin: 0, end: points.count - 1) } return result
}
let points = [
Point(x: 0.0, y: 0.0), Point(x: 1.0, y: 0.1), Point(x: 2.0, y: -0.1), Point(x: 3.0, y: 5.0), Point(x: 4.0, y: 6.0), Point(x: 5.0, y: 7.0), Point(x: 6.0, y: 8.1), Point(x: 7.0, y: 9.0), Point(x: 8.0, y: 9.0), Point(x: 9.0, y: 9.0)
] print("\(ramerDouglasPeucker(points: points, epsilon: 1.0))")</lang>
- Output:
[(0.0, 0.0), (2.0, -0.1), (3.0, 5.0), (7.0, 9.0), (9.0, 9.0)]
Wren
<lang ecmascript>import "/dynamic" for Tuple
var Point = Tuple.create("Point", ["x", "y"])
var rdp // recursive rdp = Fn.new { |l, eps|
var x = 0 var dMax = -1 var p1 = l[0] var p2 = l[-1] var x21 = p2.x - p1.x var y21 = p2.y - p1.y var i = 0 for (p in l[1..-1]) { var d = (y21*p.x - x21*p.y + p2.x*p1.y - p2.y*p1.x).abs if (d > dMax) { x = i + 1 dMax = d } i = i + 1 } if (dMax > eps) { return rdp.call(l[0..x], eps) + rdp.call(l[x..-1], eps)[1..-1] } return [l[0], l[-1]]
}
var points = [
Point.new(0, 0), Point.new(1, 0.1), Point.new(2, -0.1), Point.new(3, 5), Point.new(4, 6), Point.new(5, 7), Point.new(6, 8.1), Point.new(7, 9), Point.new(8, 9), Point.new(9, 9)
] System.print(rdp.call(points, 1))</lang>
- Output:
[(0, 0), (2, -0.1), (3, 5), (7, 9), (9, 9)]
Yabasic
<lang Yabasic>sub perpendicularDistance(tabla(), i, ini, fin)
local dx, cy, mag, pvx, pvy, pvdot, dsx, dsy, ax, ay
dx = tabla(fin, 1) - tabla(ini, 1) dy = tabla(fin, 2) - tabla(ini, 2) //Normalise mag = (dx^2 + dy^2)^0.5 if mag > 0 dx = dx / mag : dy = dy / mag pvx = tabla(i, 1) - tabla(ini, 1) pvy = tabla(i, 2) - tabla(ini, 2) //Get dot product (project pv onto normalized direction) pvdot = dx * pvx + dy * pvy //Scale line direction vector dsx = pvdot * dx dsy = pvdot * dy //Subtract this from pv ax = pvx - dsx ay = pvy - dsy return (ax^2 + ay^2)^0.5
end sub
sub DouglasPeucker(PointList(), ini, fin, epsilon)
local dmax, index, i, d // Find the point with the maximum distance
for i = ini + 1 to fin d = perpendicularDistance(PointList(), i, ini, fin) if d > dmax index = i : dmax = d next
// If max distance is greater than epsilon, recursively simplify if dmax > epsilon then PointList(index, 3) = true // Recursive call DouglasPeucker(PointList(), ini, index, epsilon) DouglasPeucker(PointList(), index, fin, epsilon) end if
end sub
data 0,0, 1,0.1, 2,-0.1, 3,5, 4,6, 5,7, 6,8.1, 7,9, 8,9, 9,9
dim matriz(10, 3)
for i = 1 to 10
read matriz(i, 1), matriz(i, 2)
next
DouglasPeucker(matriz(), 1, 10, 1)
matriz(1, 3) = true : matriz(10, 3) = true for i = 1 to 10
if matriz(i, 3) print matriz(i, 1), matriz(i, 2)
next</lang>
zkl
<lang zkl>fcn perpendicularDistance(start,end, point){ // all are tuples: (x,y) -->|d|
dx,dy := end .zipWith('-,start); // deltas dpx,dpy := point.zipWith('-,start); mag := (dx*dx + dy*dy).sqrt(); if(mag>0.0){ dx/=mag; dy/=mag; } p,dsx,dsy := dx*dpx + dy*dpy, p*dx, p*dy; ((dpx - dsx).pow(2) + (dpy - dsy).pow(2)).sqrt()
}
fcn RamerDouglasPeucker(points,epsilon=1.0){ // list of tuples --> same
if(points.len()==2) return(points); // but we'll do one point d:=points.pump(List, // first result/element is always zero fcn(p, s,e){ perpendicularDistance(s,e,p) }.fp1(points[0],points[-1])); index,dmax := (0.0).minMaxNs(d)[1], d[index]; // minMaxNs-->index of min & max if(dmax>epsilon){ return(RamerDouglasPeucker(points[0,index],epsilon)[0,-1].extend( RamerDouglasPeucker(points[index,*],epsilon))) } else return(points[0],points[-1]);
}</lang> <lang zkl>RamerDouglasPeucker(
T( T(0.0, 0.0), T(1.0, 0.1), T(2.0, -0.1), T(3.0, 5.0), T(4.0, 6.0), T(5.0, 7.0), T(6.0, 8.1), T(7.0, 9.0), T(8.0, 9.0), T(9.0, 9.0) ))
.println();</lang>
- Output:
L(L(0,0),L(2,-0.1),L(3,5),L(7,9),L(9,9))
References