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.

 First Prev Next
 How does it work ? YvesDaoust23-May-16 21:31 YvesDaoust 23-May-16 21:31
 I would have liked to see an explanation about how recursion is avoided, as this seems to be the purpose of this implementation. And doesn't that degrade the time complexity ?
 Re: How does it work ? veen_rp25-May-16 10:37 veen_rp 25-May-16 10:37
 Re: How does it work ? Member 1027179616-Mar-17 14:56 Member 10271796 16-Mar-17 14:56
 Re: How does it work ? veen_rp28-Mar-17 20:09 veen_rp 28-Mar-17 20:09
 Last Visit: 31-Dec-99 18:00     Last Update: 29-May-23 12:08 Refresh 1