Ramer-Douglas-Peucker line simplification: Difference between revisions

Content added Content deleted

Inline

Revision as of 11:49, 28 August 2022

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

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)]

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

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}]

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)]

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]]

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)

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

The computation for the perpendicular distance was taken from the GO example.

/*REXX program uses the  Ramer─Douglas─Peucker (RDP)  line simplification algorithm  for*/
/*───────────────────────────── reducing the number of points used to define its shape. */
parse arg epsilon pts                            /*obtain optional arguments from the CL*/
if epsilon='' | epsilon=","   then epsilon= 1    /*Not specified?  Then use the default.*/
if pts=''  then pts= '(0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)'
pts= space(pts)                                  /*elide all superfluous blanks.        */
            say '  error threshold: '   epsilon  /*echo the error threshold to the term.*/
            say ' points specified: '   pts      /*  "   "    shape points   "  "    "  */
$= RDP(pts)                                      /*invoke Ramer─Douglas─Peucker function*/
            say 'points simplified: '   rez($)   /*display points with () ───► terminal.*/
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
bld: parse arg _; #= words(_); dMax=-#; idx=1; do j=1  for #; @.j= word(_, j); end; return
px:  parse arg _;          return word( translate(_, , ','),  1)   /*obtain the X coörd.*/
py:  parse arg _;          return word( translate(_, , ','),  2)   /*   "    "  Y   "   */
reb: parse arg a,b,,_;                  do k=a  to b;  _= _ @.k;    end;   return strip(_)
rez: parse arg z,_;   do k=1  for words(z); _= _ '('word(z, k)") "; end;   return strip(_)
/*──────────────────────────────────────────────────────────────────────────────────────*/
RDP: procedure expose epsilon;    call bld  space( translate(arg(1), , ')(][}{') )
     L= px(@.#) - px(@.1)
     H= py(@.#) - py(@.1)                        /* [↓] find point IDX with max distance*/
                         do i=2  to #-1
                         d= abs(H*px(@.i) - L*py(@.i) + px(@.#)*py(@.1) - py(@.#)*px(@.1))
                         if d>dMax  then do;   idx= i;   dMax= d
                                         end
                         end   /*i*/             /* [↑]  D is the perpendicular distance*/

     if dMax>epsilon  then do;   r= RDP( reb(1, idx) )
                                 return subword(r, 1, words(r) - 1)     RDP( reb(idx, #) )
                           end
     return @.1  @.#

output when using the default inputs:

  error threshold:  1
 points specified:  (0,0) (1,0.1) (2,-0.1) (3,5) (4,6) (5,7) (6,8.1) (7,9) (8,9) (9,9)
points simplified:  (0,0)  (2,-0.1)  (3,5)  (7,9)  (9,9)

Rust

An implementation of the algorithm

#[derive(Copy, Clone)]
struct Point {
    x: f64,
    y: f64,
}

use std::fmt;

impl fmt::Display for Point {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "({}, {})", self.x, self.y)
    }
}

// Returns the distance from point p to the line between p1 and p2
fn perpendicular_distance(p: &Point, p1: &Point, p2: &Point) -> f64 {
    let dx = p2.x - p1.x;
    let dy = p2.y - p1.y;
    (p.x * dy - p.y * dx + p2.x * p1.y - p2.y * p1.x).abs() / dx.hypot(dy)
}

fn rdp(points: &[Point], epsilon: f64, result: &mut Vec<Point>) {
    let n = points.len();
    if n < 2 {
        return;
    }
    let mut max_dist = 0.0;
    let mut index = 0;
    for i in 1..n - 1 {
        let dist = perpendicular_distance(&points[i], &points[0], &points[n - 1]);
        if dist > max_dist {
            max_dist = dist;
            index = i;
        }
    }
    if max_dist > epsilon {
        rdp(&points[0..=index], epsilon, result);
        rdp(&points[index..n], epsilon, result);
    } else {
        result.push(points[n - 1]);
    }
}

fn ramer_douglas_peucker(points: &[Point], epsilon: f64) -> Vec<Point> {
    let mut result = Vec::new();
    if points.len() > 0 && epsilon >= 0.0 {
        result.push(points[0]);
        rdp(points, epsilon, &mut result);
    }
    result
}

fn main() {
    let points = vec![
        Point { x: 0.0, y: 0.0 },
        Point { x: 1.0, y: 0.1 },
        Point { x: 2.0, y: -0.1 },
        Point { x: 3.0, y: 5.0 },
        Point { x: 4.0, y: 6.0 },
        Point { x: 5.0, y: 7.0 },
        Point { x: 6.0, y: 8.1 },
        Point { x: 7.0, y: 9.0 },
        Point { x: 8.0, y: 9.0 },
        Point { x: 9.0, y: 9.0 },
    ];
    for p in ramer_douglas_peucker(&points, 1.0) {
        println!("{}", p);
    }
}

Output:

(0, 0)
(2, -0.1)
(3, 5)
(7, 9)
(9, 9)

Using the geo crate

The geo crate contains an implementation of the Ramer-Douglas-Peucker line simplification algorithm.

// [dependencies]
// geo = "0.14"

use geo::algorithm::simplify::Simplify;
use geo::line_string;

fn main() {
    let points = line_string![
        (x: 0.0, y: 0.0),
        (x: 1.0, y: 0.1),
        (x: 2.0, y: -0.1),
        (x: 3.0, y: 5.0),
        (x: 4.0, y: 6.0),
        (x: 5.0, y: 7.0),
        (x: 6.0, y: 8.1),
        (x: 7.0, y: 9.0),
        (x: 8.0, y: 9.0),
        (x: 9.0, y: 9.0),
    ];
    for p in points.simplify(&1.0) {
        println!("({}, {})", p.x, p.y);
    }
}

Output:

(0, 0)
(2, -0.1)
(3, 5)
(7, 9)
(9, 9)

Sidef

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
}

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]]
)

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)]

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