Jump to content

Ramer-Douglas-Peucker line simplification

From Rosetta Code


Task
Ramer-Douglas-Peucker line simplification
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



11l

Translation of: Go
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

Translation of: Go
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

Translation of: FreeBASIC
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#

Translation of: Java
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

Translation of: C++
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

Translation of: C
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

Translation of: Yabasic
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

Translation of: Kotlin
Works with: Java version 9
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

Translation of: Go
/**
 * @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

Works with: Julia version 0.6
Translation of: C++
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

Translation of: C++
// 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

Translation of: C++
Works with: Openscad version 2015.03
Works with: Openscad version 2019.05
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]]

PascalABC.NET

Translation of: C#
type
  Point = (real, real);

function PerpendicularDistance(pt, lineStart, lineEnd: Point): real;
begin
  var dx := lineEnd[0] - lineStart[0];
  var dy := lineEnd[1] - lineStart[1];
  var mag := sqrt(dx * dx + dy * dy);
  if mag > 0 then begin
    dx /= mag;
    dy /= mag;
  end;
  var pvx := pt[0] - linestart[0];
  var pvy := pt[1] - linestart[1];
  var pvdot := dx * pvx + dy * pvy;
  var ax := pvx - pvdot * dx;
  var ay := pvy - pvdot * dy;
  result := sqrt(ax * ax + ay * ay);
end;

function rdp(Pointlist: list<Point>; ε: real): list<Point>;
begin
  if Pointlist.Count < 2 then exit;
  var dmax := 0.0;
  var index := 0;
  var endindex := Pointlist.Count - 1;
  for var i := 1 to endindex do
  begin
    var d := PerpendicularDistance(Pointlist[i], Pointlist[0], Pointlist[endindex]);
    if d > dmax then begin
      index := i;
      dmax := d;
    end;
  end;
  var output := new list<point>; 
  if dmax > ε then begin
    var firstline := Pointlist.Take(index + 1).ToList;
    var lastline := Pointlist.Skip(index).ToList;
    var recresults1 := rdp(firstline, ε);
    var recresults2 := rdp(lastline, ε);
    output.AddRange(recresults1.Take(recresults1.Count - 1));
    output.AddRange(recresults2);
  end
  else begin
    output.Clear();
    output.Add(Pointlist[0]);
    output.Add(Pointlist[pointlist.Count - 1]);
  end;
  result := output;
end;

begin
  var pointlist: List<Point> := 
    Lst((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));
  Println(rdp(pointlist, 1.0));
end.
Output:
[(0,0),(2,-0.1),(3,5),(7,9),(9,9)]

Perl

Translation of: Raku
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

Translation of: Go
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

Translation of: C++
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

Library: Shapely

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

Works with: QBasic version 1.1
Works with: QuickBasic version 4.5
Works with: QB64
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)

Works with: Rakudo version 2017.05
Translation of: C++
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

Translation of: 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=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

Works with: Halcyon Calc version 4.2.9
≪ 
   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

Translation of: Kotlin
// 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

Translation of: Raku
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

Translation of: Go
Library: Wren-dynamic
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

Translation of: C
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

Translation of: Raku
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

Cookies help us deliver our services. By using our services, you agree to our use of cookies.