A C# Implementation of Douglas-Peucker Line Approximation Algorithm

CraigSelbert

Rate me:

4.73/5 (24 votes)

6 Jun 2007MIT

156.8K

3.5K

DP Line approximation algorithm is a well-known method to approximate 2D lines. It is quite fast, O(nlog_2(n)) for a n-points line and can drastically compress a data curve. Here, a fully OOP implementation is given.

Download source - 22.49 KB

Introduction

I have found multiple implementations of the Douglas-Peucker algorithm but not in any .NET language, so I decided to port it over. Jonathan de Halleux has a wonderful explanation here.

Background

I needed to reduce polygons size to display on a map based on zoom levels.

Using the Code

The code included is complete and should run out of the box in Visual Studio 2005. If it does not, please let me know.

/// <summary>
/// Uses the Douglas Peucker algorithm to reduce the number of points.
/// </summary>
/// <param name="Points">The points.</param>
/// <param name="Tolerance">The tolerance.</param>
/// <returns></returns>
public static List<Point> DouglasPeuckerReduction
    (List<Point> Points, Double Tolerance)
{
    if (Points == null || Points.Count < 3)
    return Points;

    Int32 firstPoint = 0;
    Int32 lastPoint = Points.Count - 1;
    List<Int32> pointIndexsToKeep = new List<Int32>();

    //Add the first and last index to the keepers
    pointIndexsToKeep.Add(firstPoint);
    pointIndexsToKeep.Add(lastPoint);

    //The first and the last point cannot be the same
    while (Points[firstPoint].Equals(Points[lastPoint]))
    {
        lastPoint--;
    }

    DouglasPeuckerReduction(Points, firstPoint, lastPoint, 
    Tolerance, ref pointIndexsToKeep);

    List<Point> returnPoints = new List<Point>();
    pointIndexsToKeep.Sort();
    foreach (Int32 index in pointIndexsToKeep)
    {
        returnPoints.Add(Points[index]);
    }

    return returnPoints;
}
    
/// <summary>
/// Douglases the peucker reduction.
/// </summary>
/// <param name="points">The points.</param>
/// <param name="firstPoint">The first point.</param>
/// <param name="lastPoint">The last point.</param>
/// <param name="tolerance">The tolerance.</param>
/// <param name="pointIndexsToKeep">The point index to keep.</param>
private static void DouglasPeuckerReduction(List<Point> 
    points, Int32 firstPoint, Int32 lastPoint, Double tolerance, 
    ref List<Int32> pointIndexsToKeep)
{
    Double maxDistance = 0;
    Int32 indexFarthest = 0;
    
    for (Int32 index = firstPoint; index < lastPoint; index++)
    {
        Double distance = PerpendicularDistance
            (points[firstPoint], points[lastPoint], points[index]);
        if (distance > maxDistance)
        {
            maxDistance = distance;
            indexFarthest = index;
        }
    }

    if (maxDistance > tolerance && indexFarthest != 0)
    {
        //Add the largest point that exceeds the tolerance
        pointIndexsToKeep.Add(indexFarthest);
    
        DouglasPeuckerReduction(points, firstPoint, 
        indexFarthest, tolerance, ref pointIndexsToKeep);
        DouglasPeuckerReduction(points, indexFarthest, 
        lastPoint, tolerance, ref pointIndexsToKeep);
    }
}

/// <summary>
/// The distance of a point from a line made from point1 and point2.
/// </summary>
/// <param name="pt1">The PT1.</param>
/// <param name="pt2">The PT2.</param>
/// <param name="p">The p.</param>
/// <returns></returns>
public static Double PerpendicularDistance
    (Point Point1, Point Point2, Point Point)
{
    //Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)|   *Area of triangle
    //Base = v((x1-x2)²+(x1-x2)²)                               *Base of Triangle*
    //Area = .5*Base*H                                          *Solve for height
    //Height = Area/.5/Base

    Double area = Math.Abs(.5 * (Point1.X * Point2.Y + Point2.X * 
    Point.Y + Point.X * Point1.Y - Point2.X * Point1.Y - Point.X * 
    Point2.Y - Point1.X * Point.Y));
    Double bottom = Math.Sqrt(Math.Pow(Point1.X - Point2.X, 2) + 
    Math.Pow(Point1.Y - Point2.Y, 2));
    Double height = area / bottom * 2;

    return height;
    
    //Another option
    //Double A = Point.X - Point1.X;
    //Double B = Point.Y - Point1.Y;
    //Double C = Point2.X - Point1.X;
    //Double D = Point2.Y - Point1.Y;
    
    //Double dot = A * C + B * D;
    //Double len_sq = C * C + D * D;
    //Double param = dot / len_sq;
    
    //Double xx, yy;
    
    //if (param < 0)
    //{
    //    xx = Point1.X;
    //    yy = Point1.Y;
    //}
    //else if (param > 1)
    //{
    //    xx = Point2.X;
    //    yy = Point2.Y;
    //}
    //else
    //{
    //    xx = Point1.X + param * C;
    //    yy = Point1.Y + param * D;
    //}
    
    //Double d = DistanceBetweenOn2DPlane(Point, new Point(xx, yy));
}

Points of Interest

The code is not overly complicated. It was fun to port this algorithm, since all I feel I do nowadays is hold business' hands to help them try and solve business problems that there is no consensus on.

History

Beta 1

License

This article, along with any associated source code and files, is licensed under The MIT License

Written By

CraigSelbert

Web Developer

United States

This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

More efficient version in Loyc.Math + Bug fix

Qwertie22-Nov-23 9:03

Qwertie

22-Nov-23 9:03

The above code seems to have a bug in its PerpendicularDistance method. For example, consider

PerpendicularDistance(new Point(0, 0), new Point(0, 100), new Point(20, 199))

In order for the line-simplification algorithm to work correctly, this should return 101 (because the closest point to (20, 199) on the line segment is (0, 100), and the distance between those is 101). Instead, it returns 20, which is the what the distance to the line would be if the line were infinitely long.

The version in this article also has a couple of inefficiencies:

The algorithm itself is done in three stages: (1) figure out which indexes to keep, but out-of-order, (2) sort the list of indexes, (3) built the output list of points. But in fact, the Douglas-Peucker algorithm can be done in a single stage that directly outputs the list of points in the correct order
It's not necessary to compute the distance from each point to the reference line. It's more efficient to compute the quadrance instead (quadrance is the square of distance) to avoid square root calculations.
It uses Math.Pow(x, 2) ... I'm pretty sure it's faster to multiply x * x.

The Loyc.Math package on NuGet avoids the bug and these inefficiencies in its LineMath.SimplifyPolyline static methods (although to be honest, the modified distance-to-line calculation, while avoiding the bug I mentioned, does use a division, which is not ideal - let me know if you figure out how to avoid that).

The Loyc version is also generic (designed to be compatible with any kind of point and any collection type) because .NET has several different Point types (e.g. WinForms points, WPF points). However, Loyc.Math doesn't take any of those GUI libraries as a dependency and therefore doesn't support them directly.

Here's the core algorithm (which supports any point type):

public static int SimplifyPolyline<List, Point, T>(
    List points, int iStart, int iStop,
    ICollection<Point> output, T tolerance,
    Func<Point, Point, Point, T> distanceToLine)
    where List : IReadOnlyList<Point> where T : IComparable<T>
{
    int iLast = iStop - 1;
    if (iStart >= iLast) {
        if (iStart == iLast)
            output.Add(points[iStart]);

        return iStop - iStart;
    }

    // Run Douglas-Peucker algorithm
    int count = SimplifyRecursively(points, iStart, iLast, output, tolerance, distanceToLine);

    // The last point is not included by <code>SimplifyRecursively</code>, so add it afterward.
    output.Add(points[iLast]);
    return count + 1;

    static int SimplifyRecursively(
        List points, int iFirst, int iLast, 
        ICollection<Point> output, T tolerance,
        Func<Point, Point, Point, T> distanceToLine)
    {
        Debug.Assert(iFirst < iLast);

        Point first = points[iFirst];
        int i = iFirst + 1;
        if (i < iLast) {
            Point last = points[iLast];
            T maxDist = tolerance;
            int iFarthest = -1;
            do {
                var dist = distanceToLine(points[i], first, last);
                if (maxDist.CompareTo(dist) < 0) {
                    maxDist = dist;
                    iFarthest = i;
                }
            } while (++i < iLast);

            if (iFarthest != -1) {
                int count = SimplifyRecursively(points, iFirst, iFarthest, output, tolerance, distanceToLine);
                count += SimplifyRecursively(points, iFarthest, iLast, output, tolerance, distanceToLine);
                return count;
            }
        }

        output.Add(first);
        return 1;
    }
}

(That's weird, CodeProject replaced `SimplifyRecursively` with <code>SimplifyRecursively</code> in a comment)

That method then needs to be specialized with a distance-to-line or quadrance-to-line calculation. For example, one specialization looks like this:

/// <summary>Simplifies a polyline using the Douglas-Peucker line 
///   simplification algorithm. This algorithm removes points that are 
///   deemed unimportant, so the output is a subset of the input.</summary>
/// <typeparam name="List">Original unsimplified polyline</typeparam>
/// <param name="output">The output polyline is added in order to this collection</param>
/// <param name="tolerance">The distance between the input polyline and the 
///   output polyline will never exceed this distance. Increase this value to 
///   simplify more aggressively.</param>
/// <returns>The number of output points.</returns>
/// <remarks>
///   The average time complexity of this algorithm is O(N log N). 
///   The worst-case time complexity is O(N^2).
/// </remarks>
public static int SimplifyPolyline<List>(List points, ICollection<Point<double>> output, double tolerance)
    where List : IReadOnlyList<Point<double>>
{
    return SimplifyPolyline(points, 0, points.Count, output, tolerance * tolerance, _quadranceToLineD);
}

public static List<Point<double>> SimplifyPolyline<List>(List points, double tolerance)
    where List : IReadOnlyList<Point<double>>
{
    var output = new List<Point<double>>();
    SimplifyPolyline(points, 0, points.Count, output, tolerance * tolerance, _quadranceToLineD);
    return output;
}

static readonly Func<Point<double>, Point<double>, Point<double>, double> _quadranceToLineD
    = (p, p0, p1) => (p - ProjectOnto(p, new LineSegment<double>(p0, p1), LineType.Segment, out _)).Quadrance();

But this is specialized for Point<double> which is not among the various standard .NET point types. If you want to use this code but make it compatible with some other kind of Point type (for example, System.Windows.Point), you'll need to copy the methods that the above code depends on (such as ProjectOnto() and Quadrance()), and change those other methods to meet your needs. For reference, here is the other code that the above code depends on:

public enum LineType { 
    Segment, // Finite line
    Ray,     // Line extended to infinity beyond end point
    Infinite // Line extended to infinity in both directions
}

/// <summary>Holds a 2D line segment.</summary>
/// <typeparam name="T">Coordinate type</typeparam>
public struct LineSegment<T> where T : IConvertible, IEquatable<T>
{
    public LineSegment(Point<T> a, Point<T> b) { A = a; B = b; }

    public Point<T> A, B;
    ...
}

// Meanwhile, in a different C# source file...
using T = System.Double;
using LineSegment = LineSegment<double>;
using Point = Point<double>;
using Vector = Vector<double>;
using System;

public static partial class LineMath
{
    /// <summary>Performs projection, which finds the point on a line segment 
    ///   or infinite line that is nearest to a specified point.</summary>
    /// <param name="seg">The line segment</param>
    /// <param name="p">The test point to be projected</param>
    /// <param name="type">Whether to treat the line segment as extended to infinite length.</param>
    /// <param name="end">Set to 0 if the point is on the line segment (including
    ///   one of the endpoints), -1 if the point is before seg.A, 1 if the point is 
    ///   after seg.B, and null if the line segment is degenerate (seg.A==seg.B)</param>
    /// <returns>The projected point.</returns>
    /// <remarks>
    ///   This algorithm is fast and accurate, and can be easily adapted to 3D.
    ///   A special advantage of this approach is that it runs fastest when the 
    ///   point is projected onto one of the endpoints (when infiniteLine is 
    /// false).
    /// <para/>
    ///   Algorithm comes from: <a href="http://geomalgorithms.com/a02-_lines.html">http://geomalgorithms.com/a02-_lines.html</a>
    ///   See section "Distance of a Point to a Ray or Segment"
    /// </remarks>
    public static Point ProjectOnto(this Point p, in LineSegment seg, LineType type, out int? end)
    {
        end = 0;
        Vector v = seg.Vector();
        Vector w = p.Sub(seg.A);
        T c1 = w.Dot(v); // c1 == |w|*|v|*cos(angle between them)
        if (c1 <= 0) { // angle between line segment and (p-seg.A) is negative (-180..0)?
            if (v.X == 0 && v.Y == 0) {
                // seg.A == seg.B
                end = null;
                return seg.A;
            } else if (c1 < 0)
                end = -1;
            if (type != LineType.Infinite)
                return seg.A;
        }
        T c2 = v.Quadrance(); // == |v|*|v|
        if (c1 >= c2) { // quadrance from seg.A to projected point >= quadrance of seg
            if (c1 > c2)
                end = 1;
            if (type == LineType.Segment)
                return seg.B;
        }
        if (c2 == 0) {
            // seg.A and seg.B are infitessimally close together; c2 was truncated to zero
            end = null;
            return seg.A;
        }

        T frac = c1 / c2;                         // == |w|/|v|*cos(angle)
        Point projected = seg.A.Add(v.Mul(frac)); // == p0 + v/|v|*|w|*cos(angle)
        return projected;
    }

    /// <summary>Returns seg.B - seg.A.</summary>
    public static Vector Vector(this in LineSegment seg) { return seg.B.Sub(seg.A); }
    ...
}

public static partial class PointMath
{
    public static T Quadrance(this Vector v) { return v.X*v.X + v.Y*v.Y; }

    public static Vector Sub(this Point a, Point b)   { return new Vector(a.X-b.X, a.Y-b.Y); }

    public static Point  Add(this Point a, Vector b)  { return new Point(a.X+b.X, a.Y+b.Y); }

    public static Vector Mul(this Vector p, T factor) { return new Vector(p.X*factor, p.Y*factor); }
}

(Yes, eventually I'll update these methods to support the new .NET Generic Numerics support in .NET 7, but for the time being it's still backward compatible all the way to .NET 4.5, even though the library is configured to use C# 10 which is much newer than .NET 4.5)

modified 22-Nov-23 17:14pm.

Addition to get rid of possible StackOverflow Exception

Alexander Panfilenok12-Oct-23 4:08

Alexander Panfilenok

12-Oct-23 4:08

add this If statement to get rid of StackOverflow exception generated because of float format (for those who uses float)

private static float PerpendicularDistance(Point lineStart, Point lineEnd, Point point) {
     if (point.Equals(lineStart) || point.Equals(lineEnd)) return 0;

Thanks!

Keith Rule21-May-22 8:19

Keith Rule

21-May-22 8:19

I'm trying to turn rotoscoped images into ila (laser vector files). I've found I need two open-source libraries to do this. CSPtrace and this one. The results are pretty good and this article help decimate the polygons into a number that a laser projector is likely able to project.

So thanks for writing this. It saved me from having to code the algorithm from the Wikipedia article.

Keith

Keith Rule

rejecting recursion + distance calc acceleration

adevder18-Feb-17 23:12

adevder

18-Feb-17 23:12

public static List<Point> DouglasPeuckerReduction(List<Point> Points, Double Tolerance)
        {
            double Tolerancesqrd = Tolerance * Tolerance;
            if (Points == null || Points.Count < 3)
                return Points;

            Int32 firstPoint = 0;
            Int32 lastPoint = Points.Count - 1;
            List<Int32> pointIndexsToKeep = new List<Int32>();

            //Add the first and last index to the keepers
            pointIndexsToKeep.Add(firstPoint);
            pointIndexsToKeep.Add(lastPoint);


            //The first and the last point can not be the same
            while (lastPoint >= 0 && Points[firstPoint].Equals(Points[lastPoint]))
            {
                lastPoint--;
            }
            if (lastPoint == 0) { return Points; }

            SortedDictionary<Int32, Int32> PairsIndexesToCheck = new SortedDictionary<Int32, Int32>();
            PairsIndexesToCheck.Add(firstPoint, lastPoint);
            Int32 checkcnt = 0;
            while (PairsIndexesToCheck.Count > 0&& checkcnt< lastPoint)
            {
                Double maxDistancesqrd = 0;
                Int32 indexFarthest = 0, currentFirstPoint= PairsIndexesToCheck.First().Key, currentLastPoint= PairsIndexesToCheck.First().Value;
                double deltax = Points[currentFirstPoint].X - Points[currentLastPoint].X,
                deltay = Points[currentFirstPoint].Y - Points[currentLastPoint].Y;
                Double oneoverbottomsqrd = 1 / (deltax * deltax + deltay * deltay),
                 x1y2 = Points[currentFirstPoint].X * Points[currentLastPoint].Y,
                 x2y1 = Points[currentFirstPoint].Y * Points[currentLastPoint].X,
                 x1y2_diff_x2y1 = x1y2 - x1y2;
                 ;
                for (Int32 index = currentFirstPoint+1; index < currentLastPoint; index++)
                {
                    Double distancesqrd = PerpendicularDistance(Points[currentFirstPoint], Points[currentLastPoint], Points[index], oneoverbottomsqrd, x1y2_diff_x2y1, deltax, deltay);
                    if (distancesqrd > maxDistancesqrd)
                    {
                        maxDistancesqrd = distancesqrd;
                        indexFarthest = index;
                    }
                }

                if (maxDistancesqrd > Tolerancesqrd )
                {
                    //Add the largest point that exceeds the tolerance
                    pointIndexsToKeep.Add(indexFarthest);
                    //split current pair
                    PairsIndexesToCheck[currentFirstPoint] = indexFarthest;
                    PairsIndexesToCheck.Add(indexFarthest, currentLastPoint);

                }
                else
                {
                    //pair is checked
                    PairsIndexesToCheck.Remove(currentFirstPoint);
                }
                checkcnt++;
            }
            if (checkcnt == lastPoint) { return Points; }


            List<Point> returnPoints = new List<Point>();
            pointIndexsToKeep.Sort();

            foreach (Int32 index in pointIndexsToKeep)
            {
                returnPoints.Add(Points[index]);
            }

            return returnPoints;
        }


<pre>public static Double PerpendicularDistance(Point Point1, Point Point2, Point Point, Double oneoverbottomsqrd, Double x1y2_diff_x2y1, Double deltax, Double deltay)
        {
            //Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)|   *Area of triangle вект произв 
            //Base = √((x1-x2)²+(x1-x2)²)                               *Base of Triangle*
            //Area = .5*Base*H                                          *Solve for height
            //Height = Area/.5/Base

            Double area = Math.Abs(x1y2_diff_x2y1 -  Point.Y* deltax + Point.X * deltay );
            Double areasqrd = area * area;
          
            Double heightsqrd = areasqrd * oneoverbottomsqrd ;

            return heightsqrd;

            

        }

Re: rejecting recursion + distance calc acceleration

TerrySpiderman21-Mar-17 2:05

TerrySpiderman

21-Mar-17 2:05

Thanks for making a non-recursive version. You have a typo though

x1y2_diff_x2y1 = x1y2 - x1y2;

should be

x1y2_diff_x2y1 = x1y2 - x2y1;

to save stack memory with great amount of points and tolerance near zero i prefer to reject recursion

adevder18-Feb-17 0:22

adevder

18-Feb-17 0:22

public static List<Point> DouglasPeuckerReduction(List<Point> Points, Double Tolerance)
        {

            if (Points == null || Points.Count < 3)
                return Points;

            Int32 firstPoint = 0;
            Int32 lastPoint = Points.Count - 1;
            List<Int32> pointIndexsToKeep = new List<Int32>();

            //Add the first and last index to the keepers
            pointIndexsToKeep.Add(firstPoint);
            pointIndexsToKeep.Add(lastPoint);


            //The first and the last point can not be the same
            while (lastPoint >= 0 && Points[firstPoint].Equals(Points[lastPoint]))
            {
                lastPoint--;
            }
            if (lastPoint == 0) { return Points; }

            SortedDictionary<Int32, Int32> PairsIndexesToCheck = new SortedDictionary<Int32, Int32>();
            PairsIndexesToCheck.Add(firstPoint, lastPoint);
            while (PairsIndexesToCheck.Count > 0)
            {
                Double maxDistance = 0;
                Int32 indexFarthest = 0, currentFirstPoint= PairsIndexesToCheck.First().Key, currentLastPoint= PairsIndexesToCheck.First().Value;

                for (Int32 index = currentFirstPoint+1; index < currentLastPoint; index++)
                {
                    Double distance = PerpendicularDistance(Points[currentFirstPoint], Points[currentLastPoint], Points[index]);
                    if (distance > maxDistance)
                    {
                        maxDistance = distance;
                        indexFarthest = index;
                    }
                }

                if (maxDistance > Tolerance )
                {
                    //Add the largest point that exceeds the tolerance
                    pointIndexsToKeep.Add(indexFarthest);
                    //split current pair
                    PairsIndexesToCheck[currentFirstPoint] = indexFarthest;
                    PairsIndexesToCheck.Add(indexFarthest, currentLastPoint);

                }
                else
                {
                    //pair is ok
                    PairsIndexesToCheck.Remove(currentFirstPoint);
                }
            }

        

            List<Point> returnPoints = new List<Point>();
            pointIndexsToKeep.Sort();
            foreach (Int32 index in pointIndexsToKeep)
            {
                returnPoints.Add(Points[index]);
            }

            return returnPoints;
        }

Nice

veen_rp18-Apr-16 18:41

veen_rp

18-Apr-16 18:41

Nice concise code. Here's my VB.NET conversion to whom it may concern:

Public Function GetSimplifyDouglasPeucker(points As List(Of PointD), tolerance As Double) As List(Of PointD)
	If points Is Nothing OrElse points.Count < 3 Then Return points 'nothing to do
	
	Dim pi1 As Integer = 0
	Dim pi2 As Integer = points.Count-1
	Dim nli As New List(Of Integer)
	
	nli.Add(pi1)
	nli.Add(pi2)
	
	'first and last point cannot be the same
	While points(pi1).IsCoincident(points(pi2))
		pi2 -= 1
	End While
	
	recurseDouglasPeuckerTrimming(points,pi1,pi2,tolerance,nli)
	
	Dim nl As New List(Of PointD)
	nli.Sort()
	For Each idx As Integer In nli
		nl.Add(points(idx))
	Next
	
	Return nl
	
End Function

Private Sub recurseDouglasPeuckerTrimming(points As List(Of PointD), firstpoint As Integer, lastpoint As Integer, tolerance As Double, ByRef pIndexKeep As List(Of Integer))
	Dim maxDist As Double = 0
	Dim idxFarthest As Integer = 0
	Dim distance As Double = 0
	
	For i As Integer = firstpoint To lastpoint-1
		distance = GetTriangleHeight(points(firstPoint),points(lastPoint),points(i))
		If distance > maxDist Then
			maxDist = distance
			idxFarthest = i
		End If
	Next
	
	If maxDist > tolerance AndAlso idxFarthest <> 0 Then
		'Add the largest point that exceeds the tolerance
		pIndexKeep.Add(idxFarthest)
		
		'recurse
		recurseDouglasPeuckerTrimming(points,firstpoint,idxFarthest,tolerance,pIndexKeep)
		recurseDouglasPeuckerTrimming(points,idxFarthest,lastpoint,tolerance,pIndexKeep)
	End If
End Sub

Public Function GetTriangleArea(base1 As PointD, base2 As PointD, apex As PointD) As Double
	Try
		Return Math.Abs(.5*(base1.X*base2.Y + base2.X*apex.Y + apex.X*base1.Y - base2.X*base1.Y - apex.X*base2.Y - base1.X*apex.Y))
	Catch
		Return Double.NaN
	End Try
End Function

Public Function GetTriangleBase(base1 As PointD, base2 As PointD, apex As PointD) As Double
	Try
		Return ((base1.X-base2.X)^2 + (base1.Y-base2.Y)^2)^0.5
	Catch
		Return Double.NaN
	End Try
End Function

Public Function GetTriangleHeight(base1 As PointD, base2 As PointD, apex As PointD) As Double
	Return GetTriangleArea(base1,base2,apex)/GetTriangleBase(base1,base2,apex)*2
End Function

If (first = last) Bug!

Giovani Luigi14-Dec-15 11:44

Giovani Luigi

14-Dec-15 11:44

There is a problem in your code.

The code:

//The first and the last point cannot be the same
while (Points[firstPoint].Equals(Points[lastPoint]))
{
    lastPoint--;
}

Should be executed before:

//Add the first and last index to the keepers
pointIndexsToKeep.Add(firstPoint);
pointIndexsToKeep.Add(lastPoint);

Try with a polygon with start and end points equal. (closed polygon)

Do we agree?

Re: If (first = last) Bug!

Keith Rule21-May-22 12:40

Keith Rule

21-May-22 12:40

Thanks for this comment. It saved me from having to troubleshoot that weird last point being off in some random direction issue. This fix solved the problem. Thanks for saving me a bunch of troubleshooting.

Keith Rule

DouglasPeuckerN() conversion?

cbordeman11-Dec-15 16:42

cbordeman

11-Dec-15 16:42

This is pretty cool, but could you also convert the DouglasPeuckerN() variant method? Instead of taking a threshold, it takes the exact number of points you need in the output. Most of the code is identical to what you've done already.

Of course, one could TRY to do a search algorithm using the algorithm you converted, but it's very easy for this to fail. I've tried converting the Doublas Peucker variant, but I just can't, my C++ knowledge is too limited.

PLEASE!

Re: DouglasPeuckerN() conversion?

Mr Lorenz27-Apr-17 6:44

Mr Lorenz

27-Apr-17 6:44

Hi,
I have the same problem and I develop a first version using a max numbers of points to show.

public static List<Point> DouglasPeuckerReduction(List<Point> points, double nPoints)
        {
            if (points == null || points.Count < 3)
                return points;
            int firstIndex = 0;
            int lastIndex = points.Count - 1;
            List<int> pointIndexsToKeep = new List<int>();
            //Add the first and last index to the keepers
            pointIndexsToKeep.Add(firstIndex);
            pointIndexsToKeep.Add(lastIndex);
            //The first and the last point can not be the same
            while (points[firstIndex].Equals(points[lastIndex]))
                lastIndex--;
            //start check
            List<SubPoly> polyList = new List<SubPoly>();
            PolyOrderBy listOrder = new PolyOrderBy();
            SubPoly subPoly = new SubPoly(firstIndex, lastIndex);
            subPoly.KeyInfo = DouglasPeuckerReduction(points, firstIndex, lastIndex);
            polyList.Add(subPoly);
            while (polyList.Count != 0)
            {
                subPoly = polyList[0];
                polyList.RemoveAt(0);
                pointIndexsToKeep.Add(subPoly.KeyInfo.Index);
                if (pointIndexsToKeep.Count == nPoints)
                    break;
                //left
                SubPoly left = new SubPoly(subPoly.FirstIndex, subPoly.KeyInfo.Index);
                left.KeyInfo = DouglasPeuckerReduction(points, left.FirstIndex, left.LastIndex);
                if (left.KeyInfo.Index != 0)
                    polyList.Add(left);
                //right
                SubPoly right = new SubPoly(subPoly.KeyInfo.Index, subPoly.LastIndex);
                right.KeyInfo = DouglasPeuckerReduction(points, right.FirstIndex, right.LastIndex);
                if (right.KeyInfo.Index != 0)
                    polyList.Add(right);
                //order queue by distance
                polyList.Sort(listOrder);
            }
            List<Point> returnPoints = new List<Point>();
            pointIndexsToKeep.Sort();
            foreach (int index in pointIndexsToKeep)
                returnPoints.Add(points[index]);
            //
            return returnPoints;
        }

        private static KeyInfo DouglasPeuckerReduction(List<Point> points, int firstIndex, int lastIndex)
        {
            double maxDistance = 0;
            int indexFarthest = 0;
            for (int index = firstIndex; index < lastIndex; index++)
            {
                double distance = PerpendicularDistance(points[firstIndex], points[lastIndex], points[index]);
                if (distance > maxDistance)
                {
                    maxDistance = distance;
                    indexFarthest = index;
                }
            }
            KeyInfo kInfo = new KeyInfo();
            kInfo.Distance = maxDistance;
            kInfo.Index = indexFarthest;
            return kInfo;
        }

        private class PolyOrderBy : IComparer<SubPoly>
        {
            public int Compare(SubPoly x, SubPoly y)
            {
                return y.KeyInfo.Distance.CompareTo(x.KeyInfo.Distance);
            }
        }

        public class KeyInfo
        {
            public int Index { get; set; }
            public double Distance { get; set; }
        }

        public class SubPoly
        {
            public int FirstIndex { get; set; }
            public int LastIndex { get; set; }
            public KeyInfo KeyInfo { get; set; }

            public SubPoly(int firstIndex, int lastIndex)
            {
                this.FirstIndex = firstIndex;
                this.LastIndex = lastIndex;
                this.KeyInfo = new KeyInfo();
            }
        }

Lorenzo

Useful

Dineshkumar_MCA22-Apr-15 1:18

Dineshkumar_MCA

22-Apr-15 1:18

BUG tolerance

Member 1017584518-Feb-15 1:08

Member 10175845

18-Feb-15 1:08

i have an issue regarding tolerance that,i get tolerance from tileX,tileY and zoom level,it works well up to zoom level 15 but after that
"An unhandled exception of type 'System.StackOverflowException' occurred in RDP_ GUI.exe".
this error appear.

Thanks + pruned code

Stoyan Haralampiev23-Oct-14 2:30

Stoyan Haralampiev

23-Oct-14 2:30

Thanks for the useful work!

In order to use it in my project I made the following changes:
* Improved readability by renaming some variables and updating comments
* Added the two checks suggested in matthew minke's comment below
* Removed Math.Pow() which seems to hurt performance. For details see http://codereview.stackexchange.com/questions/29002/ramer-douglas-peucker-algorithm[^]).

Using this algorithm on data from database

Saket Surya2-Jun-14 2:38

Saket Surya

2-Jun-14 2:38

I want to use this algorithm to reduce data points on a plot where data is fetched from database to plot. Can someone give an efficient way to do this?
My table has records in millions, and I would like to fetch only a few thousand to plot the chart.

Point projected on the line but behind the segment? Bug?

Sergiy Tkachuk25-Mar-14 13:31

Sergiy Tkachuk

25-Mar-14 13:31

I have not found the check if point is projected to line which contains segment, but actual projection lies behind segment.

Do you perform such check?
In such case tolerance should be calculated as a distance between point and one of the ends of segment, and should _not_ be considered as a length of the projection to the line.

Without such check wrong behaviour could appear if point is close to the line but far behind segment at the same time.
In such case point could be removed by algorithms, but in fact it should be preserved.

About Tolerance

kamlesh Prajapati8-Jan-14 22:55

kamlesh Prajapati

8-Jan-14 22:55

can you just tell me value about tolerance as default for google maps data and jsut explain about x and y of point? is it lattitude and longitude value?

Why do you create your own Point?

leiyangge3-Jan-14 21:37

leiyangge

3-Jan-14 21:37

In .net there's already System.Drawing.PointF, which I think is precise enough for you demo.

License

Diego Catalano14-Dec-13 1:42

Diego Catalano

14-Dec-13 1:42

Great work you have here!

I was deeply interested on it, so as much as that I would like to know if, perhaps, you would be willing to share your work under different licenses. For example, would you like to share your work under the LGPL? In this way I would be able to use it under other applications I have been working lately!

Thanks!
Diego

1 bug, + 1 potential bug

gyuri1023-Apr-12 22:45

gyuri10

23-Apr-12 22:45

First of all: GREAT WORK!

1 bug:

In the public DouglasPeuckerReduction method should have the loop like this:

//The first and the last point cannot be the same
while (Points[firstPoint].Equals(Points[lastPoint]))
{
    lastPoint--;
    pointIndexsToKeep.Add(lastPoint);  // <- this line has been inserted
}

If you do not add the very last point again (and it is a significant point), you loose it.

+1 potential bug:
In the PerpendicaularDistance calculation, System.Drawing.Point is used, where X and Y are Int32. Should X and Y be a not that big number (e.g. 1 million), during calculation one might get easily a bit overflow. It is nicer if X and Y are treated as Int64.

modified 24-Apr-12 5:52am.

Re: 1 bug, + 1 potential bug

tplokas26-Mar-13 13:46

tplokas

26-Mar-13 13:46

I think you are mistaken on both counts.

1) The loop was eliminating identical points. The fact that your modification adds them into the list of indexes to keep is working against the very purpose of reducing the number of points. Nothing is lost because the original value of lastIndex was already added to the list of indexes to keep, before this loop.

2) I don't see where System.Drawing.Point is being used at all. This code is using a custom class named Point whose X and Y are of type Double.

Useful,Thanks!

masy19112-Mar-12 22:43

masy191

12-Mar-12 22:43

I really need it,thank you very much!

Great Job!

Christopher Drake19-Dec-11 7:32

Christopher Drake

19-Dec-11 7:32

I implemented this in a pet silverlight project. It works well.

how calculate an epsilon value

unravel-the-sky23-May-11 5:52

unravel-the-sky

23-May-11 5:52

thank you very much for the douglas-peucker algorithm, i have been working on a project that uses that. i have a question, is there way to calculate an optimum epsilon value within the algorithm? or does it have to be user defined?

Some Issues with your implementation

MattMi15-Jul-10 10:57

MattMi

15-Jul-10 10:57

I have found some issues with your implementation of this code. here are the changes that i have made.

Double maxDistance = 0;
Int32 indexFarthest = 0;

if (lastPointIndex - firstPointIndex > 1)//ADDITION: need to have more than two points in the set we are looking through
{
    for (Int32 index = firstPointIndex; index < lastPointIndex; index++)
    {
        Double distance = PerpendicularDistance(points[firstPointIndex], points[lastPointIndex], points[index]);
        if (distance > maxDistance)
        {
            maxDistance = distance;
            indexFarthest = index;
        }
    }

    if (maxDistance > Tolerance && indexFarthest != firstPointIndex)//CHANGE: condition was wrong.
    {
        //Add the largest point that exceeds the tolerance
        pointIndexsToKeep.Add(indexFarthest);

        Reduction(points, firstPointIndex, indexFarthest, ref pointIndexsToKeep);
        Reduction(points, indexFarthest, lastPointIndex, ref pointIndexsToKeep);
    }
}

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.