15,663,557 members
Articles / General Programming / Computational Geometry
Tip/Trick
Posted 18 Apr 2016

16.4K views
2 bookmarked

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

Rate me:

## Introduction

This tip on the Douglas-Peucker algorithm is an add-on to two already published articles on Code Project: Polyline Simplification by and A C# Implementation of Douglas Peucker by . 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 .

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.

VB.NET
```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)

'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.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
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`:

VB.NET
```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
End Function

Public Function ToStringXY(fmt As String) As String
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

Written By
Engineer VeeTools
Netherlands
A (usually exploring) geologist who sometimes develops software for fun and colleagues... Check out my new website at www.veetools.xyz for online mapping and coordinate conversion.

 The recursion, which is used in the code of Craig Selbert, is avoided by iteratively (i.e. sequentially and not nested) calling function `iterateDPT`. With the iteration a list is build of indexes to points which cannot be 'simplified away' based on the tolerance. Once all the points have been analyzed and either kept in the index list `nli As List(Of Integer)`, or discarded (i.e. 'simplified away'), the iteration signals this by returning a `0`. At this moment there is a final `nli As List(Of Integer)` containing all the indexes to all the points which should be kept of the original polyline. Within my code there is a maximum of 3 loops nested: 1. the outer `While busyI` loop, 2. the middle `For i As Integer = 0 To nli.Count-2` loop to grow the new `nli` list and 3. the inner `For i As Integer = firstpoint To lastpoint` loop in the `iterateDPT` function With recursion every time a new valid `farpoint` is encountered a new search for a valid `farpoint` is started withing the parent search. If you have a long and complex polyline the level of nesting can become very deep. Moreover you cannot control the nesting depth from you own code; it manages itself (because it is recursive, obviously ). I'm not sure if I make my point clear with the above explanation; anyway it's a try. The best way to really grasp it is to put Craig Selbert's code next to mine and spot the differences... Hope this helps