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
<lang perl6>sub perpendicular-distance (@start, @end where @end !eqv @start , @point) {
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 ($index, $dmax) = (0,0); for (0 ^..^ @points) { my $d = perpendicular-distance(@points[0], @points[*-1], @points[$_]); if $d > $dmax { $dmax = $d; $index = $_ } } 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, 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)] );</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=True))</lang>
- Output:
LINESTRING (0 0, 2 -0.1, 3 5, 7 9, 9 9)