A VB.NET Implementation of the Douglas-Peucker Algorithm without Recursion

• Download source - 1.5 KB

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 'nothing to do
    
    Dim pointfirst As Integer = 0
    Dim pointlast As Integer = points.Count-1
    
    'first and last point cannot be the same
    While points(pointfirst).IsCoincident(points(pointlast))
        pointlast -= 1
    End While

    Dim nli As New List(Of Integer)
    nli.Add(pointfirst)
    nli.Add(pointlast)
    
    'loop DPT
    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
                'we can skip the rest of the for next loop
                Exit For
            End If
        Next
        
        If newFarPoint > 0 Then
            nli.Add(newFarPoint)
            nli.Sort()
        Else
            busyI = False
        End If
        
    End While
    
    'process and output the results
    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
    
    'the iteration should stop if either the tolerance is not exceeded or
    'when there are no more far (outlier) points
    'this is indicated by the output being either 0 or a valid farpoint
    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

• Written in April 2016