## 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<sub>n</sub>`

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 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
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
If points.Count < 3 Then Throw New Exception("Catmull-Rom Spline requires at least 3 points")
If nrOfInterpolatedPoints < 1 Then nrOfInterpolatedPoints = 1
Dim spoints As New List(Of PointD)
For Each p As PointD In points
spoints.Add(New PointD(p))
Next
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))
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
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
spline.Add(spoints(spoints.Count-2))
Return spline
Catch exc As Exception
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<sub>n</sub>`

to `P<sub>n+1</sub>`

and insert it at the end of `spoints`

. The reason I do this is because the Catmull-Rom spline calculates point `P<sub>t</sub>`

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<sub>t</sub>`

'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 2 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)
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
Dim cutdist As Double = 0.05 + (tension*0.4)
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)
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
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. 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.

Catmull-Rom on a closed polyline

Chaikin on a closed polyline

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

Happy coding!

## History

Article written in april 2016.