A C# Implementation of Douglas-Peucker Line Approximation Algorithm

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

• 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.

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

/// <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();

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); }
}

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

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;

            

        }

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

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

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

//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--;
    pointIndexsToKeep.Add(lastPoint);  // <- this line has been inserted
}

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);
    }
}