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
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
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
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
Catmull-Rom spline and a non-spline calculation for
approximation invented by Chaikin and called the
CornerCutting algorithm. For a very readable article on the latter method see: Chaikin’s Algorithms for Curves by Kenneth I. Joy.
First of all it is convenient to have a
Point class with
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()
Public Sub New(nx As Double, ny As Double)
X = nx
Y = ny
Public Sub New(p As PointD)
X = p.X
Y = p.Y
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)
Public Shared Operator +(p As PointD, d As Double) As PointD
Return New PointD(p.X+d,p.Y+d)
Public Shared Operator +(d As Double, p As PointD) As PointD
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)
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
Dim dx As Double = spoints(1).X-spoints(0).X
Dim dy As Double = spoints(1).Y-spoints(0).Y
dx = spoints(spoints.Count-1).X-spoints(spoints.Count-2).X
dy = spoints(spoints.Count-1).Y-spoints(spoints.Count-2).Y
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)
Catch exc As Exception
The adaptions I have made are:
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+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
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.
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
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
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
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
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.
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
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
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
For i As Integer = 1 To nrOfIterations
nl = getSmootherChaikin(nl,cutdist)
Private Function getSmootherChaikin(points As List(Of PointD), cuttingDist As Double) As List(Of PointD)
Dim nl As New List(Of PointD)
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)
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.
Article written in april 2016.