Ramer-Douglas-Peucker line simplification

Ramer–Douglas–Peucker algorithm is a line simplification algorithm for reducing the number of points used to define its shape.^[1]

Task

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) using the Ramer–Douglas–Peucker algorithm. The error threshold to use is 1.0. Display the remaining points.

C++

<lang cpp>#include <iostream>

include <cmath>
include <utility>
include <vector>
include <stdexcept>

using namespace std;

typedef std::pair<double, double> Point;

double PerpendicularDistance(const Point &pt, const Point &lineStart, const Point &lineEnd) { double dx = lineEnd.first - lineStart.first; double dy = lineEnd.second - lineStart.second;

//Normalise double mag = pow(pow(dx,2.0)+pow(dy,2.0),0.5); if(mag > 0.0) { dx /= mag; dy /= mag; }

double pvx = pt.first - lineStart.first; double pvy = pt.second - lineStart.second;

//Get dot product (project pv onto normalized direction) double pvdot = dx * pvx + dy * pvy;

//Scale line direction vector double dsx = pvdot * dx; double dsy = pvdot * dy;

//Subtract this from pv double ax = pvx - dsx; double ay = pvy - dsy;

return pow(pow(ax,2.0)+pow(ay,2.0),0.5); }

void RamerDouglasPeucker(const vector<Point> &pointList, double epsilon, vector<Point> &out) { if(pointList.size()<2) throw invalid_argument("Not enough points to simplify");

// Find the point with the maximum distance from line between start and end double dmax = 0.0; size_t index = 0; size_t end = pointList.size()-1; for(size_t i = 1; i < end; i++) { double d = PerpendicularDistance(pointList[i], pointList[0], pointList[end]); if (d > dmax) { index = i; dmax = d; } }

// If max distance is greater than epsilon, recursively simplify if(dmax > epsilon) { // Recursive call vector<Point> recResults1; vector<Point> recResults2; vector<Point> firstLine(pointList.begin(), pointList.begin()+index+1); vector<Point> lastLine(pointList.begin()+index, pointList.end()); RamerDouglasPeucker(firstLine, epsilon, recResults1); RamerDouglasPeucker(lastLine, epsilon, recResults2);

// Build the result list out.assign(recResults1.begin(), recResults1.end()-1); out.insert(out.end(), recResults2.begin(), recResults2.end()); if(out.size()<2) throw runtime_error("Problem assembling output"); } else { //Just return start and end points out.clear(); out.push_back(pointList[0]); out.push_back(pointList[end]); } }

int main() { vector<Point> pointList; vector<Point> pointListOut;

pointList.push_back(Point(0.0, 0.0)); pointList.push_back(Point(1.0, 0.1)); pointList.push_back(Point(2.0, -0.1)); pointList.push_back(Point(3.0, 5.0)); pointList.push_back(Point(4.0, 6.0)); pointList.push_back(Point(5.0, 7.0)); pointList.push_back(Point(6.0, 8.1)); pointList.push_back(Point(7.0, 9.0)); pointList.push_back(Point(8.0, 9.0)); pointList.push_back(Point(9.0, 9.0));

RamerDouglasPeucker(pointList, 1.0, pointListOut);

cout << "result" << endl; for(size_t i=0;i< pointListOut.size();i++) { cout << pointListOut[i].first << "," << pointListOut[i].second << endl; }

return 0; }</lang>

Output:

result
0,0
2,-0.1
3,5
7,9
9,9

Perl 6

Works with: Rakudo version 2016.11

Translation of: C++

<lang perl6>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) «/=» ($Δx² + $Δy²).sqrt;
   ((($Δpx, $Δpy) «-» ($Δx, $Δy) «*» ($Δx*$Δpx + $Δy*$Δpy)) «**» 2).sum.sqrt;

}

sub Ramer-Douglas-Peucker( @points where * > 1, \ε = 1 ) {

   return @points if @points == 2;
   my @d = (^@points).map:
      { perpendicular-distance(@points[0], @points[*-1], @points[$_]) };
   my ($index, $dmax) = @d.first: @d.max, :kv;

   if $dmax > ε {
       return Ramer-Douglas-Peucker( @points[0..$index]  , ε )[0..*-2],
              Ramer-Douglas-Peucker( @points[$index..*-1], ε )
   } else {
       return @points[0, *-1];
   }

}

TESTING

say flat Ramer-Douglas-Peucker(

  [(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]

);</lang>

Output:

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

Python

An approach using the shapely library:

<lang python>from __future__ import print_function from shapely.geometry import LineString

if __name__=="__main__": line = LineString([(0,0),(1,0.1),(2,-0.1),(3,5),(4,6),(5,7),(6,8.1),(7,9),(8,9),(9,9)]) print (line.simplify(1.0, preserve_topology=False))</lang>

Output:

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

Sidef

Translation of: Perl 6

<lang ruby>func perpendicular_distance(Arr start, Arr end, Arr point) {

   ((point == start) || (point == end)) && return 0
   var (Δx,  Δy ) = (  end »-« start)...
   var (Δpx, Δpy) = (point »-« start)...
   var h = hypot(Δx, Δy)
   [\Δx, \Δy].map { *_ /= h }
   (([Δpx, Δpy] »-« ([Δx, Δy] »*» (Δx*Δpx + Δy*Δpy))) »**» 2).sum.sqrt

}

func Ramer_Douglas_Peucker(Arr points { .all { .len > 1 } }, ε = 1) {

   points.len == 2 && return points

   var d = (^points -> map {
       perpendicular_distance(points[0], points[-1], points[_])
   })

   if (d.max > ε) {
       var i = d.index(d.max)
       return [Ramer_Douglas_Peucker(points.ft(0, i), ε).ft(0, -2)...,
               Ramer_Douglas_Peucker(points.ft(i),    ε)...]
   }

   return [points[0,-1]]

}

say Ramer_Douglas_Peucker(

   [[0,0],[1,0.1],[2,-0.1],[3,5],[4,6],[5,7],[6,8.1],[7,9],[8,9],[9,9]]

)</lang>

Output:

[[0, 0], [2, -1/10], [3, 5], [7, 9], [9, 9]]

zkl

Translation of: Perl 6

<lang zkl>fcn perpendicularDistance(start,end, point){ // all are tuples: (x,y) -->|d|

  dx,dy   := end  .zipWith('-,start);	// deltas
  dpx,dpy := point.zipWith('-,start);
  mag     := (dx*dx + dy*dy).sqrt();
  if(mag>0.0){ dx/=mag; dy/=mag; }
  p,dsx,dsy := dx*dpx + dy*dpy, p*dx, p*dy;
  ((dpx - dsx).pow(2) + (dpy - dsy).pow(2)).sqrt()

}

fcn RamerDouglasPeucker(points,epsilon=1.0){ // list of tuples --> same

  if(points.len()==2) return(points);  // but we'll do one point
  d:=points.pump(List,  // first result/element is always zero
     fcn(p, s,e){ perpendicularDistance(s,e,p) }.fp1(points[0],points[-1]));
  index,dmax := (0.0).minMaxNs(d)[1], d[index]; // minMaxNs-->index of min & max
  if(dmax>epsilon){
      return(RamerDouglasPeucker(points[0,index],epsilon)[0,-1].extend(
             RamerDouglasPeucker(points[index,*],epsilon)))
  } else return(points[0],points[-1]);

}</lang> <lang zkl>RamerDouglasPeucker(

  T( T(0.0, 0.0), T(1.0, 0.1), T(2.0, -0.1), T(3.0, 5.0), T(4.0, 6.0),
     T(5.0, 7.0), T(6.0, 8.1), T(7.0,  9.0), T(8.0, 9.0), T(9.0, 9.0) ))

.println();</lang>

Output:

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

References

↑ [1]

[1] [1]

[1]