2D Polyline Vertex Smoothing

• Download source - 10.7 KB

Note: The VB.NET code of this article has been ported to C# by xtos (Chris C) on GitHub, see PolylineSmoothCSharp, thanks Chris.

Introduction & Background

When working with map (GIS) or chart data, you will have objects in the shape of 2D points, lines, polylines and polygons. These objects are known under a lot of different names: shapes, paths, areas, regions, etc. In this article, I will consider a point as a single x,y pair or vertex or simply P; a line as a pair of vertices with startpoint P1 and endpoint P2); a polyline as a collection of multiple vertices with P1 to P_n where n > 2. In a polyline, the vertices are connected in the order they appear in the collection. A polygon is simply a polyline where the startpoint is connected to the endpoint or a closed polyline. See figure below for a graphic explanation.

When drawing these objects, there is often a need to smooth the vertices of a polyline. I realise with GDI there are the DrawCurve and DrawBezier methods, but this article will delve a bit deeper into the generation of the smoothed polylines and it will show how to deal with more than just a fixed number of vertices (which is 4 in the case of DrawBezier) and how to tweak the 'smoothness'.

Smoothing a polyline can be done in two ways: (a) by interpolation, meaning that the original polyline points will left unchanged and in the new smoothed polyline, and (b) by approximation, meaning that the new smoothed polyline will approximate the old polyline but the original points will not be preserved. The first is also known as a cardinal or canonical spline, the second is often solved by quadratic or cubic Bezier curves. Again, see the next figure for a graphic explanation.

There is lots of information to be found on the web on spline, Beziers, etc. Obvious starting point: Wikipedia on Splines and a blog a specifically liked: Splines and curves part I – Bézier Curves by Martin Doms. Here on CodeProject, there is an excellent in depth article: Spline Interpolation - history, theory and implementation by Kenneth Haugland (by the way, I am a typical specimen of what he refers to as a lazy programmer). Splines and Bezier curves are very powerful because they are parametric: they behave as a mathematical function and can be used to estimate, approximate or interpolate data points in between the original data points. In fact, they can serve as an alternative to for example Linear or Polynomial Regression of datapoints.

However, when smoothing a polyline for plotting purposes, you are not necessarily interested in the underlying function but more in performance and simplicity of use. There is also a drawback when using splines in the sense that you use 3 or 4 control points or vertices at a time in a piece-wise approach. If you want to apply it on a whole series of points, you somehow have to chain the individual curves together. It can and has been done, but it gets complex easily.

This article will show the use of a simple spline for interpolation called Catmull-Rom spline and a non-spline calculation for approximation invented by Chaikin and called the Chaikin or CornerCutting algorithm. For a very readable article on the latter method, see Chaikin’s Algorithms for Curves by Kenneth I. Joy.

The Code

First of all, it is convenient to have a Point class with X and Y As Double and also define standard mathematical operators for this class to be able to code PointC = PointA * PointB. This is my 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 PointD)
        X = p.X
        Y = p.Y
    End Sub
    
    Public X As Double = 0
    Public Y As Double = 0
    
    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
    
    'and similar operators for "-", "*" and "/"
    '...
    
End Class

The full class is provided in the source file. Secondly, when working with polyline drawing, I prefer using the Chart Windows forms control in the System.Windows.Forms.DataVisualization namespace which since .NET 4 is built-in. This has all been set up in the source code of the project provided for download.

1. Catmull-Rom Interpolation

There are numerous examples of coding the Catmull-Rom spline on the web, see for example my reference to the article of Kenneth Haugland on CodeProject above. I have taken the code from Nice Curves! and adapted it to:

Private Function getSplineInterpolationCatmullRom_
(points As List(Of PointD), nrOfInterpolatedPoints As Integer) As List(Of PointD)
    Try
        'The Catmull-Rom Spline, requires at least 4 points so it is possible
        'to extrapolate from 3 points, but not from 2. you would get a straight line anyway
        If points.Count < 3 Then Throw New Exception_
                            ("Catmull-Rom Spline requires at least 3 points")
        
        'could throw an error on the following, but it is easily fixed implicitly
        If nrOfInterpolatedPoints < 1 Then nrOfInterpolatedPoints = 1
        
        'part1 
        'create a new pointlist to do splining on
        'if you don't do this, the original pointlist gets extended 
        'with the exptrapolated points
        Dim spoints As New List(Of PointD)
        For Each p As PointD In points
            spoints.Add(New PointD(p))
        Next
        
        'always extrapolate the first and last point out
        Dim dx As Double = spoints(1).X-spoints(0).X
        Dim dy As Double = spoints(1).Y-spoints(0).Y
        spoints.Insert(0,New PointD(spoints(0).X-dx,spoints(0).Y-dy))
        dx = spoints(spoints.Count-1).X-spoints(spoints.Count-2).X
        dy = spoints(spoints.Count-1).Y-spoints(spoints.Count-2).Y
        spoints.Insert(spoints.Count, _
            New PointD(spoints(spoints.Count-1).X+dx,spoints(spoints.Count-1).Y+dy))
        
        'part2
        'Note the nrOfInterpolatedPoints acts as a kind of tension factor 
        'between 0 and 1 because
        'it is normalised to 1/nrOfInterpolatedPoints.
        Dim t As Double = 0
        Dim spoint As PointD
        Dim spline As New List(Of PointD)
        
        For i As Integer = 0 To spoints.Count-4
            spoint = New PointD()
            For intp As Integer = 0 To nrOfInterpolatedPoints-1
                'define t as a fraction between point 2 and 3 based on the nr of
                'interpolated points desired
                t = 1/nrOfInterpolatedPoints*intp
                
                spoint = 0.5*( _
                    2 * spoints(i+1) + (-1 * spoints(i) + spoints(i+2)) * t + _
                    (2 * spoints(i) - 5 * spoints(i+1) + _
                                      4 * spoints(i+2) - spoints(i+3)) * t^2 + _
                    (-1 * spoints(i) + 3 * spoints(i+1) - _
                                       3 * spoints(i+2) + spoints(i+3)) * t^3)
                
                spline.Add(New PointD(spoint))
            Next
            
        Next
        
        'part3
        'add the last point, but skip the interpolated last point, so second last...
        spline.Add(spoints(spoints.Count-2))
        Return spline
    Catch exc As Exception
        'Debug.Print(exc.ToString)
        Return Nothing
    End Try
End Function

The adaptions I have made are:

1. See part1 in the code section above: Before iterating through the polyline vertices to smooth by interpolation, I create an extrapolation from P1 of the polyline to P0 and insert it at the beginning of a copy of the original pointlist (i.e., spoints As New List(Of PointD)). And a similar extrapolation from P_n to P_n+1 and insert it at the end of spoints. The reason I do this is because the Catmull-Rom spline calculates point P_t on the spline between P2 and P3 taking piece-wise sections of 4 points from the original polyline. And by extrapolating the first and last point of the original polyline, I preserve these points in the smoothed interpolated output.
2. See part2 in the code section above: The actual number of P_t's calculated depends on the number of interpolated points desired. t is a normalised distance between P2 and P3 of which the coordinates are calculated by the spline function. So the outer For Next loop iterates through spoints taking 4 points at a time and the inner For Next loop calculates interpolated points at t = 1/nrOfInterpolatedPoints*intp distance between P2 and P3 of the current 4 points taken out of the polyline defined by spoints. Example: if I would need to calculate 3 interpolated points between P2 and P3 of the polyline, I would generate t = 1/3*0 = 0, 1/3*1 = 0.3333 and 1/3*2 = 0.6667 distance between P2 and P3, calculate the new point spoint and add it to the end of the new pointlist spline (I admit that nrOfInterpolatedPoints should rather have been called nrOfInterpolatedSegments). Note the input parameter serves as a sort of smoothing degree: the more points you interpolate between 2 original vertices, the smoother the line in between will be. However, this reaches an optimum and the smoothness of the polyline is not affected anymore beyond that optimum.

The last step at part3 in the code section above adds the final point to the new spline pointlist. Note that this is the second last (or spoints.Count-2) point of the spoints pointlist because we do not want to add our extrapolated last point of part1. Examples of the Catmull-Rom spline with varying number of interpolation points is shown in the two figures below:

2. Chaikin

Sometimes you do not want an interpolation of the vertices in a polyline because it is noisy, jagged, jittery or has a couple of strange outliers. In this case, an approximation is called for. Of course, you could use spline approximation but there is a more simple and, yes, less mathematical approach called Corner Cutting. From the referenced article I have given above by Kenneth I. Joy is the following figure:

From each segment of the polyline, you define a new segment at 1/4 from the startpoint and at 3/4 (or at (1-1/4) from the endpoint). If you iterate this a few times, your polyline will be approximated smoothly between the original points. Simple and elegant. Couple of notes though: (a) the method is not parametric in the sense that you can continously calculate each point on the new smooth polyline and (b) start- and endpoint P0 and P4 are cut from the new polyline. The code for this is relatively simple and again, I have made a few adaptions:

Private Function getCurveSmoothingChaikin(points As List(Of PointD), _
        tension As Double, nrOfIterations As Integer) As List(Of PointD)
    'checks
    If points Is Nothing OrElse points.Count < 3 Then Return Nothing
    
    If nrOfIterations < 1 Then nrOfIterations = 1
    
    If tension < 0 Then
        tension = 0
    ElseIf tension > 1
        tension = 1
    End If
    
    'the tension factor defines a scale between corner cutting distance in segment half length,
    'i.e. between 0.05 and 0.45. The opposite corner will be cut by the inverse
    '(i.e. 1-cutting distance) to keep symmetry.
    'with a tension value of 0.5 this amounts to 0.25 = 1/4 and 0.75 = 3/4,
    'the original Chaikin values

    Dim cutdist As Double = 0.05 + (tension*0.4)
    
    'make a copy of the pointlist and iterate it
    Dim nl As New List(Of PointD)
    For i As Integer = 0 To points.Count-1
        nl.Add(New PointD(points(i)))
    Next
    
    For i As Integer = 1 To nrOfIterations
        nl = getSmootherChaikin(nl,cutdist)
    Next
    
    Return nl
End Function

Private Function getSmootherChaikin(points As List(Of PointD), cuttingDist As Double) _
        As List(Of PointD)
    Dim nl As New List(Of PointD)
    
    'always add the first point
    nl.Add(New PointD(points(0)))
    
    Dim q, r As PointD
    
    For i As Integer = 0 To points.Count-2
        q = (1-cuttingDist)*points(i) + cuttingDist*points(i+1)
        r = cuttingDist*points(i) + (1-cuttingDist)*points(i+1)
        nl.Add(q)
        nl.Add(r)
    Next
    
    'always add the last point
    nl.Add(New PointD(points(points.Count-1)))
    
    Return nl
End Function

The adaptions I made are:

1. The code consists of the wrapper getCurveSmoothingChaikin() which defines the cutting distance or cutdist from the corner. In Chaikin's algorithm, this is at 0.25 and 0.75 relative distance. I have transformed this into a tension factor between 0 and 1 which defines a cutdist between 0.05 and 0.45 relative distance by calculating 0.05+(tension*0.4), the opposite or complementary corner is then 1-cutdist. A value of t=0.5 gives the original Chaikin cutting distance of 0.25.
2. 2. The second adaption is the iterated method getSmootherChaikin() where I do add the original start and endpoint of the polyline to the new smoothed polyline.

In the above way, the smoothness of the resulting polyline and the tension from the corners of the original polyline can be controlled. Examples of the Chaikin approximation varying iteration and tension are shown in the next 3 figures.

Finally, the next 2 figures show that both methods work for a closed polyline (or polygon) in a similar way. Note the difference here between interpolation and approximation. The sharp reader will notice that the start- and endpoint in the closed polygon (i.e., the same point) is not smoothed. You could solve this by checking if the first and last point fall on top of each other and, if they do, intrapolate a line segment around this point when copying the pointlist.

The above little Windows application is found in the supplied source files.

Happy coding!

History

• 17^th April 2016: Article written
• 2^nd April, 2019: Added a link to a C# port of this article on GitHub