In this article, you will see the use of a simple spline for interpolation called Catmull-Rom spline and a non-spline calculation for approximation called the Chaikin or CornerCutting algorithm.
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 Pn 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
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:
- 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 Pn to Pn+1 and insert it at the end of spoints
. The reason I do this is because the Catmull-Rom spline calculates point Pt 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. - See
part2
in the code section above: The actual number of Pt'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)
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:
- 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
- 17th April 2016: Article written
- 2nd April, 2019: Added a link to a C# port of this article on GitHub