Introduction
This tip on the Douglas-Peucker algorithm is an add-on to two already published articles on Code Project: Polyline Simplification by Elmar de Koning and A C# Implementation of Douglas Peucker by CraigSelbert. The below posted code is a more or less straight conversion from Craig Selbert's code which is short and sweet.
The only thing about the Douglas-Peucker algorithm I do not really like is the recursion because this can potentially cause your program to overload the stack beyond the control of your managed code. As Murphy stated: anything that can go wrong, will go wrong at the worst possible moment.
So in addition to Craig's code, I have written VB.NET code for the Douglas-Peucker algorithm with iterations instead of recursion. I post the code here without much further words. For an excellent explanation on the algorithm, please refer to the article by Elmar de Koning.
The class PointD I use in this code snippet (and a lot of other methods as well) is a helper class (see also my article: 2D Polyline Vertex Smoothing), the code is given at the bottom of this tip.
The Code
The public function GetSimplifyDouglasPeucker is the wrapper for the iteration. With respect to Craig's code, I have replaced his PerpendicularDistance method with three methods giving the area, base and height of a triangle.
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
Dim pointfirst As Integer = 0
Dim pointlast As Integer = points.Count-1
While points(pointfirst).IsCoincident(points(pointlast))
pointlast -= 1
End While
Dim nli As New List(Of Integer)
nli.Add(pointfirst)
nli.Add(pointlast)
Dim farPoint, newFarPoint As Integer
Dim busyI As Boolean = True
While busyI
newFarPoint = 0
For i As Integer = 0 To nli.Count-2
farPoint = iterateDPT(points,nli(i),nli(i+1),tolerance)
If farPoint > newFarPoint Then
newFarPoint = farPoint
Exit For
End If
Next
If newFarPoint > 0 Then
nli.Add(newFarPoint)
nli.Sort()
Else
busyI = False
End If
End While
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 Function iterateDPT(points As List(Of PointD), _
firstpoint As Integer, lastpoint As Integer, tolerance As Double) As Integer
Dim maxDist As Double = 0
Dim idxFarthest As Integer = 0
Dim distance As Double = 0
For i As Integer = firstpoint To lastpoint
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
Return idxFarthest
Else
Return 0
End If
End Function
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
And the code for the helper class PointD:
Public Class PointD
Public Sub New()
End Sub
Public Sub New(nx As Double, ny As Double)
X = nx
Y = ny
End Sub
Public Sub New(p As Veet.Geometry.PointD)
X = p.X
Y = p.Y
End Sub
Public X As Double = 0
Public Y As Double = 0
Public Function ToStringXY() As String
Return ToStringXY(""," ",False)
End Function
Public Function ToStringXY(fmt As String) As String
Return ToStringXY(fmt," ",True)
End Function
Public Function ToStringXY(fmt As String, seperator As String, withLabel As Boolean) As String
Dim xstr As String = X.ToString
Dim ystr As String = Y.ToString
If fmt <> "" Then
xstr = X.ToString(fmt)
ystr = Y.ToString(fmt)
End If
If withLabel Then
xstr = "X=" & xstr
ystr = "Y=" & ystr
End If
Return String.Format("{0}{1}{2}",xstr,seperator,ystr)
End Function
Public Shared Operator +(p1 As PointD, p2 As PointD) As PointD
Return New PointD(p1.X+p2.X,p1.Y+p2.Y)
End Operator
Public Shared Operator +(p As PointD, d As Double) As PointD
Return New PointD(p.X+d,p.Y+d)
End Operator
Public Shared Operator +(d As Double, p As PointD) As PointD
Return p+d
End Operator
Public Shared Operator -(p1 As PointD, p2 As PointD) As PointD
Return New PointD(p1.X-p2.X,p1.Y-p2.Y)
End Operator
Public Shared Operator -(p As PointD, d As Double) As PointD
Return New PointD(p.X-d,p.Y-d)
End Operator
Public Shared Operator -(d As Double, p As PointD) As PointD
Return p-d
End Operator
Public Shared Operator *(p1 As PointD, p2 As PointD) As PointD
Return New PointD(p1.X*p2.X,p1.Y*p2.Y)
End Operator
Public Shared Operator *(p As PointD, d As Double) As PointD
Return New PointD(p.X*d,p.Y*d)
End Operator
Public Shared Operator *(d As Double, p As PointD) As PointD
Return p*d
End Operator
Public Shared Operator /(p1 As PointD, p2 As PointD) As PointD
Return New PointD(p1.X/p2.X,p1.Y/p2.Y)
End Operator
Public Shared Operator /(p As PointD, d As Double) As PointD
Return New PointD(p.X/d,p.Y/d)
End Operator
Public Shared Operator /(d As Double, p As PointD) As PointD
Return New PointD(d/p.X,d/p.Y)
End Operator
End Class
And that is it...
History
