Click here to Skip to main content
12,301,374 members (64,150 online)
Click here to Skip to main content
Add your own
alternative version

Stats

24.8K views
1.4K downloads
55 bookmarked
Posted

Spline Interpolation - history, theory and implementation

, 2 Dec 2014 CPOL
Rate this:
Please Sign up or sign in to vote.
Implementation of Bezier curve, Derivative Bezier curve, Cathull-Rom spline, Bessel-Overhauser spline, Lagrange interpolation and convex hull

Introduction

There are many implementations of interpolation schema based on the Bezier curve around the web, but they all seem to be either specifically oriented about one curve, or have functions that are not general enough for a wide variety of curves.

The build in Bezier curve is just the quadratic and cubic curve, and should be enough for the simple use, but in many cases this simply won't do. The goal of this article is to give a program with reusable functions that you could easily import into your program.

Background

Let's say that you wanted to create a smooth curve that went through some points, and you are stuck in the 1950's; What you'd have to do is to place ducts at the points you wanted to stay still, get a spline with the desired stiffness, and climb up on the table like this guy from Boeing:

As Microsoft writes in their MSDN article "A physical spline is a thin piece of wood or other flexible material" and voila the term Spline later got coined by Isaac Jacob Schoenberg. Further history behind the development could be read here, and the basic is that it was used many places, but the mathematics behind it didn't take off until the 50's and 60's, and I quote:

The use of splines for modeling automobile bodies seems to have several independent beginnings. Credit is claimed on behalf of de Casteljau at Citroen, Bezier at Renault, and Birkhoff, Garabedian, and de Boor at General Motors, all for work occurring in the very early 1960's or late 1950's. At least one of de Casteljau's papers was published, but not widely, in 1959. De Boor's work at GM resulted in a number of papers being published in the early 60's, including some of the fundamental work on B-splines.

Bezier curves

The most commonly know way of constructing a Bezier curve is to use De Casteljau's algorithm, which is a linear interpolation to generate the Bezier curve. However, the Bezier curve can also be defined in a different way, through the use of Bernstein polynomials.

Bernstein polynomial

The Bernstein polynomial is linked with a theorem published by Karl Weierstrass in 1885 (He was 70 years old at the time, so the theory that mathematical innovations are only done by young people (20 or younger) should perhaps be laid to rest?). He basically found out that "any continuous function on the interval [0,1] may be uniformly approximated, arbitrarily closely, by polynomials". The Bernstein polynomial is simply a different proof of Weierstrass theorem, and it can be shown trough the standard way; with the use of the binomial function:

\begin{equation} B_{i,n}(t) = {n \choose i} t^i (1-t)^{n-i} \end{equation}
However If you have taken statistics you would (or should) notice that the formula is the exact same as the Binomial mass distribution. At first, the formula might seem a bit strange and exotic, but in reality, it is really simple. The first thing that one needs to understand is the phrase ({N \choose k}), and written out this formula is simply:
\begin{equation} {N \choose k} = \frac{N!}{k! \cdot (N - k) !} \end{equation}

Explain in layman's terms it just explains that if you have N elements and want to pick out k elements from it, there is exactly ({N \choose k} ) ways of doing it. The aforementioned formula really consists of two distinct parts, the first is which explains the number of ways you can pick out k elements from N: (\frac{N!}{(N-k)!}), and the second (\frac{1}{k!}) simply eliminates permutations, number of different ways you can combine (k) elements in.

Why the factorial appears in selecting elements is actually really neat. If you have a card deck of 52 cards and you want to choose 4 cards from it, I can break the statistics of it down of each draw from the deck as follows:

  • The first time I choose a card, there is exactly 52 cards to pick from, so I have 52 equally likely choices.
  • The next time I want to take a card from the deck, there is now just 51 cards left to choose from. The number of choices I had to get exactly the two first cards is then (52 \cdot 51).
  • The third draw, 50 cards to choose from, number of combinations is (52 \cdot 51 \cdot 50).
  • The fourth time the odds are (52 \cdot 51 \cdot 50 \cdot 49).

This is quite tedious to write out, but luckily there is a simple general way of creating this series on by simply writing out:

\begin{equation} \frac{N!}{(N-k)!} = \frac{52 \cdot 51 \cdots 2 \cdot 1}{48 \cdot 47 \cdots 2 \cdot 1} \end{equation}

 

This will give the number of ways to combine 4 cards in an explicit sequence, but normally you don't care about the sequence of events you draw them in, just the resulting 4 cards. It can be shown by this simple analogy; if you have to choose 2 cards, they can either come in the sequence {1,2} or {2,1} and the two ways are called permutations in math. The number of ways to combine the sequence and still get the cards you want is 2! or in our case with 4 cards 4!. We want to eliminate number of permutations in the way of getting k so we modify the formula:

\begin{equation}{N \choose k} = \frac{N!}{k!(N-k)!}\end{equation}

The number of ways to combine them is just half the story when calculating the probability of getting the desired outcome. You also need to determine the probability of getting the desired outcome vs not getting it. The easiest way to do this with the Bernstein polynomial is to choose 2 points and combine them. There is a demand that the probability of the desired outcome plus the probability of undesired outcome must be 1 for all t. It is just one way of achieving this and thus the formula becomes:

$B_{i,n}(t) = {n \choose i} t^i (1-t)^{n-i}$

We generally don't like to take factorial with computers, as the numbers can get frighteningly large very quickly and is thus highly subjectable to roundoff errors. The formula is instead implemented by the use of the recursive relation of the Bernstein polynomial:

$B_{i,n} = (1-t) \cdot B_{i,n-1}+ t \cdot B_{i-1,n-1}$

This can be translated into code as given below:

''' <summary>
''' The code uses the recursive relation B_[i,n](u)=(1-u)*B_[i,n-1](u) + u*B_[i-1,n-1](u) to compute all nth-degree Bernstein polynomials
''' </summary>
''' <param name="n">The sum of the start point, the end point and all the knot points between</param>
''' <param name="u">Ranges from 0 to 1, and represents the current position of the curve</param>
''' <returns></returns>
''' <remarks>This code is translated to VB from the original C++  code given on page 21 in "The NURBS Book" by Les Piegl and Wayne Tiller </remarks>
Private Function AllBernstein(ByVal n As Integer, ByVal u As Double) As Double()
    Dim B(n - 1) As Double
    B(0) = 1
    Dim u1 As Double
    u1 = 1 - u
    Dim saved, temp As Double
    For j As Integer = 1 To n - 1

        saved = 0
        For k As Integer = 0 To j - 1
            temp = B(k)
            B(k) = saved + u1 * temp
            saved = u * temp
        Next
        B(j) = saved
    Next

    Return B
End Function

The code will produce all the Bernstein polynomials that are required given the number of points in the input. The curves are shown below, ranging from 1 point to 6 points given to the function (image from MathWorld):

The implementation, given that you have gotten all your Bernstein polynomials, is pretty straight forward. It is just the sum of the Bernstein polynomial and the point that belongs to it:

''' <summary>
''' Create a Bezier curve from two points, together with knot points in between the startpoint and endpoint given in the pointcollection
''' </summary>
''' <param name="p">The two points and the knots in between</param>
''' <param name="StepSize">How close should the step in the curve be</param>
''' <returns></returns>
''' <remarks></remarks>
Private Function BezierFunction(ByVal p As PointCollection, ByVal StepSize As Double) As PointCollection
    Dim result As New PointCollection
    Dim B As Double()

    For k As Double = 0 To 1 Step StepSize

        B = AllBernstein(p.Count, k)

        Dim CX, CY As Double
        CX = 0
        CY = 0
        For j As Integer = 0 To p.Count - 1
            CX = CX + B(j) * p(j).X
            CY = CY + B(j) * p(j).Y
        Next
        result.Add(New Point(CX, CY))
    Next

    Return result
End Function

Derivative of the Beziers Curve

The derivative of the Bernstein polynomial could be used to generate the tangent of any Bezier curve, as is done here in an article on Silverlight. I'm going to do the same thing only with some slight changes.

Remember that it is only the Bernstein polynomial that is dependent on changes in u. This means that the knots as well as the start and endpoint is to be considered a constant. There are several articles about how to get the derivative of the Bezier curve, however few mentions that one can use the Bernstein polynomial to get the derivative. The relation that one is interested in is:

B'i,n(u) )= n ( Bi-1,n-1(u) - Bi,n-1(u))

One must of course remember that the Bernstein polynomial is defined to be zero if i is below zero, or i is higher than n. In mathematical (Read: In short) terms:

B-1,n-1(u) = Bn,n-1(u) = 0

The code for the derivative of the Bernstein polynomial could then easily be constructed:

Private Function AllDerivateBernstein(ByVal n As Integer, ByVal u As Double) As Double()

      Dim B As Double()
      B = AllBernstein(n - 1, u)

      Dim result(n - 1) As Double

      For i As Integer = 0 To n - 1
          If i = 0 Then
              result(i) = n * (-B(i))
          ElseIf i = n - 1 Then
              result(i) = n * (B(i - 1))
          Else
              result(i) = n * (B(i - 1) - B(i))
          End If
      Next

      Return result
  End Function

Since it is only the Bernstein polynomial that provides the line in the curve, together with the control points, one can use the series expansion (or power laws) to integrate the polynomial part by part. All we have to do now is to use the derivative to calculate the slope:

Private Function DerivateBezierFunction(ByVal p As PointCollection, ByVal StepSize As Double) As PointCollection
    Dim result As New PointCollection
    Dim B As Double()

    For k As Double = 0 To 1 Step StepSize

        B = AllDerivateBernstein(p.Count, k)

        Dim test() As Double = AllDerivateBernstein(p.Count, k)
        Dim CX, CY As Double
        CX = 0
        CY = 0
        For j As Integer = 0 To p.Count - 1
            CX = CX + B(j) * p(j).X
            CY = CY + B(j) * p(j).Y
        Next
        result.Add(New Point(CX, CY))
    Next

    result.Add(p(p.Count - 1))

    Return result
End Function

All you need to do now, is to divide y'/x' and you'll get the slope at a given point. Since this is an analytical method it is valid in every step point from 0 to 1 in the Bezier curve.

So, I decided to nick the arrow picture from the article "Bezier curve angle calculation in Silverlight", and use it in my implementation. There is however, a problem, as the rotation should not revolve around the center, or rather it can, but it is easier to revolve it around its given point at 0,0 and move that point by half the Length/Width with a 90-degree turn. This ensures that the arrow would always be in the right position, and the code makes easy use of trigonometry:

''' <summary>
''' Given a startpoint, a length and the degrees in which it should travel
''' </summary>
''' <param name="StartPoint">The line would start in this point</param>
''' <param name="Length">How long should the line be?</param>
''' <param name="theta">The angle it should travel relative to a straight line along the X axis</param>
''' <returns></returns>
''' <remarks></remarks>
Private Function MovePoint(ByVal StartPoint As Point, ByVal Length As Double, Optional ByVal theta As Double = 90) As Point
    Dim EndPoint As New Point

    theta *= Math.PI / 180

    EndPoint.X = Math.Cos(theta) * (Length) - Math.Sin(theta) * (Length) + StartPoint.X
    EndPoint.Y = Math.Sin(theta) * (Length) + Math.Cos(theta) * (Length) + StartPoint.Y
    Return EndPoint
End Function

If you see in the comments section there is one that points out that one could use the Tan2 method to calculate the angle. It takes in two points, and calculates the angle in relation to a straight line on the X-axis:

Private Function AngleInDegree(ByVal start As Point, ByVal endd As Point) As Double
    Return Math.Atan2(start.Y - endd.Y, endd.X - start.X) * 180 / Math.PI
End Function

The code for placing the image and rotating it could also be simplified:

RemoveImage()
Dim image As New Image()
Dim uri As New Uri("arrow.jpg", UriKind.Relative)
Dim img As New System.Windows.Media.Imaging.BitmapImage(uri)


Dim pp As Point = MovePoint(New Point(ll.X1, ll.Y1), 10, angle - 90)
image.Margin = New Thickness(pp.X, pp.Y, 0, 0)

image.Source = img
image.Width = 20
image.Height = 20
image.Stretch = Stretch.Fill

Dim RotTransform As New RotateTransform
RotTransform.Angle = angle
image.RenderTransform = RotTransform

cnvMain.Children.Add(image)

WPF's Bezier methods

Now, there is also the possibility of using the predefined method in WPF to draw a cubic Bezier segment. MSDN gives the following example in XAML:

<Path Stroke="Black" StrokeThickness="1">
  <Path.Data>
    <PathGeometry>
      <PathGeometry.Figures>
        <PathFigureCollection>
          <PathFigure StartPoint="10,100">
            <PathFigure.Segments>
              <PathSegmentCollection>
                <BezierSegment Point1="100,0" Point2="200,200" Point3="300,100" />
              </PathSegmentCollection>
            </PathFigure.Segments>
          </PathFigure>
        </PathFigureCollection>
      </PathGeometry.Figures>
    </PathGeometry>
  </Path.Data>
</Path>

And there is a different code for the quadratic Bezier segment too (also from MSDN):

<Path Stroke="Black" StrokeThickness="1">
  <Path.Data>
    <PathGeometry>
      <PathGeometry.Figures>
        <PathFigureCollection>
          <PathFigure StartPoint="10,100">
            <PathFigure.Segments>
              <PathSegmentCollection>
                <QuadraticBezierSegment Point1="200,200" Point2="300,100" />
              </PathSegmentCollection>
            </PathFigure.Segments>
          </PathFigure>
        </PathFigureCollection>
      </PathGeometry.Figures>
    </PathGeometry>
  </Path.Data>
</Path>

Both of these versions is included if you decide to use the Bernstein polynomial instead, however one could also set up animations in XAML by using the code above.

Rational Bezier curves

If one uses the Bernstein polynomial alone, all points are equally important. If one, on the other hand, wants to assign different weight on each of the points (same as different mass in the solar system), one gets more control over the shape of the Bezier curve. With proper usage one can even replicate a circle with just three points.

So, instead of using the Bernstein polynomial directly, one multiply it by a weight, called wi . However, to ensure that the curve stays within boundaries it is important that the sum is equal to 1 for all u. This is quite simply called the rational basis function Ri,n(u) and is defined:

Private Function RationalBasisFunction(ByVal n As Integer, ByVal u As Double, ByVal weight() As Double) As Double()
    If weight.Length <> n Then Return Nothing

    Dim B() As Double
    B = AllBernstein(n, u)

    Dim result(n - 1) As Double

    Dim test As Double
    test = 0
    For j As Integer = 0 To n - 1
        test += B(j) * weight(j)
    Next

    For i As Integer = 0 To n - 1
        result(i) = B(i) * weight(i) / test

    Next

    Return result
End Function

Catmull-Rom interpolation

I know programmers can be kind of lazy, so they could just see that they need a Catmull-Rom interpolation, do a search, and find one. Like this one and grab the code:

''' <summary>
''' Calculates interpolated point between two points using Catmull-Rom Spline </summary>
''' <remarks>
''' Points calculated exist on the spline between points two and three. </remarks>
''' <param name="p0">First Point</param>
''' <param name="p1">Second Point</param>
''' <param name="p2">Third Point</param>
''' <param name="p3">Fourth Point</param>
''' <param name="t">
''' Normalised distance between second and third point where the spline point will be calculated </param>
''' <returns>Calculated Spline Point </returns>
Public Function PointOnCurve(ByVal p0 As System.Windows.Point, ByVal p1 As System.Windows.Point, ByVal p2 As System.Windows.Point, ByVal p3 As System.Windows.Point, ByVal t As Double) As System.Windows.Point
    Dim result As New System.Windows.Point()

    Dim t2 As Single = t * t
    Dim t3 As Single = t2 * t

    result.X = 0.5F * ((2.0F * p1.X) + (-p0.X + p2.X) * t + (2.0F * p0.X - 5.0F * p1.X + 4 * p2.X - p3.X) * t2 + (-p0.X + 3.0F * p1.X - 3.0F * p2.X + p3.X) * t3)
    result.Y = 0.5F * ((2.0F * p1.Y) + (-p0.Y + p2.Y) * t + (2.0F * p0.Y - 5.0F * p1.Y + 4 * p2.Y - p3.Y) * t2 + (-p0.Y + 3.0F * p1.Y - 3.0F * p2.Y + p3.Y) * t3)

    Return result
End Function

There is also another way, and that is to try and understand what the Catmull-Rom spline really is. The polynomial that is used in the code above is actually generated from the derivatives in the figure below. But it only uses the derivative to define two control points between point 1 and point 2.

The derivative is of course defined as usual:

And the points Pi+ and pi+1- is the control points designed by the derivative:

We can now calculate two new points between the points marked Pi and Pi+1, the clue now is to actually use a 3rd degree Bezier curve by the use of the start point Pi and the end point Pi+1, and two-knot points Pi+ and Pi+1-. We now understand that the Catmull-Rom spline is just a special form of 3rd degree Bezier curve. Useful links that explains this and where the figures are from is here and here.

''' <summary>
''' Calculates interpolated point between two points using Catmull-Rom Spline </summary>
''' <remarks>
''' Points calculated exist on the spline between points two and three. </remarks>
''' <param name="p0">First Point</param>
''' <param name="p1">Second Point</param>
''' <param name="p2">Third Point</param>
''' <param name="p3">Fourth Point</param>
''' <param name="t">
''' Normalised distance between second and third point where the spline point will be calculated </param>
''' <returns>Calculated Spline Point </returns>
Public Function PointOnCurve(ByVal p0 As System.Windows.Point, ByVal p1 As System.Windows.Point, ByVal p2 As System.Windows.Point, ByVal p3 As System.Windows.Point, ByVal t As Double) As System.Windows.Point
    Dim result As New System.Windows.Point()

    'Derivative calcualtions
    Dim lix1, liy1, lix2, liy2 As Double
    lix1 = 0.5 * (p2.X - p0.X)
    lix2 = 0.5 * (p3.X - p1.X)

    liy1 = 0.5 * (p2.Y - p0.Y)
    liy2 = 0.5 * (p3.Y - p1.Y)

    ' Location of Controlpoints
    Dim PointList As New PointCollection
    PointList.Add(p1)
    PointList.Add(New Point(p1.X + (1 / 3) * lix1, p1.Y + (1 / 3) * liy1))
    PointList.Add(New Point(p2.X - (1 / 3) * lix2, p2.Y - (1 / 3) * liy2))
    PointList.Add(p2)

    ' Get the calcualted value from the 3rd degree Bezier curve
    Return PointBezierFunction(PointList, t)
End Function

The problem I faced was with a single line, which I often had to deal with, so how should I interpolate those? I choose to use the logic that the weather tomorrow will be the same as the weather today, or rather half a day. I added one point prior and one point past the original line, and I said that there was a point at the same slope as between the two first points, and the two last points, only half as long away. It seem to generate a good agreement as long as you did not go completely mental in the beginning or end.

All that is left to construct is a wrapper that can take the standard PointCollection as input for the Catmull-Rom spline:

Public Function CatmullRomSpline(ByVal Points As PointCollection, Optional ByVal InterpolationStep As Double = 0.1, Optional ByVal IsPolygon As Boolean = False) As PointCollection
    Dim result As New PointCollection

    If Points.Count <= 2 Then
        Return Points
    End If

    If IsPolygon Then
        For i As Integer = 0 To Points.Count - 1
            If i = 0 Then
                For k As Double = 0 To (1 - InterpolationStep) Step InterpolationStep
                    result.Add(PointOnCurve(Points(Points.Count - 1), Points(i), Points(i + 1), Points(i + 2), k))
                Next
            ElseIf i = Points.Count - 1 Then
                For k As Double = 0 To (1 - InterpolationStep) Step InterpolationStep
                    result.Add(PointOnCurve(Points(i - 1), Points(i), Points(0), Points(1), k))
                Next
            ElseIf i = Points.Count - 2 Then
                For k As Double = 0 To (1 - InterpolationStep) Step InterpolationStep
                    result.Add(PointOnCurve(Points(i - 1), Points(i), Points(i + 1), Points(0), k))
                Next
            Else
                For k As Double = 0 To (1 - InterpolationStep) Step InterpolationStep
                    result.Add(PointOnCurve(Points(i - 1), Points(i), Points(i + 1), Points(i + 2), k))
                Next
            End If
        Next
    Else
        Dim yarray, xarray As New List(Of Double)
        xarray.Add(Points(0).X - (Points(1).X - Points(0).X) / 2)
        yarray.Add(Points(0).Y - (Points(1).Y - Points(0).Y) / 2)

        For Each ps As System.Windows.Point In Points
            xarray.Add(ps.X)
            yarray.Add(ps.Y)
        Next

        xarray.Add((Points(Points.Count - 1).X - (Points(Points.Count - 2).X) / 2 + Points(Points.Count - 1).X))
        yarray.Add((Points(Points.Count - 1).Y - (Points(Points.Count - 2).Y) / 2 + Points(Points.Count - 1).Y))

        Dim r As New PointCollection
        For i As Integer = 0 To yarray.Count - 1
            r.Add(New System.Windows.Point(xarray(i), yarray(i)))
        Next

        For i As Integer = 3 To r.Count - 1
            For k As Double = 0 To (1 - InterpolationStep) Step InterpolationStep
                result.Add(PointOnCurve(r(i - 3), r(i - 2), r(i - 1), r(i), k))
            Next
        Next
        result.Add(Points(Points.Count - 1))
    End If

    Return result
End Function

The Catmull-Rom spline is prone to rather nasty curves if the control points are very unevenly distributed. If you take a look at what happens to that curve if you take 4 points and place the two center points quite close together and have the start and endpoint quite a distance from the two center points. You'll get a sharp bend between the two center points.

It uses the Bessel tangent method or Overhausers method to generate the spline, by the use of a non-uniform weight ti. Incidentally, if you have a uniform weight it will produce the Catmull-Rom Spline instead. The values of ti is suppose to reduce the "velocity" so that a more smooth curve can appear.

I didn't find any implementation of the curve so it's created from an algorithm in the book "Mathematical tools in computer graphics with C# implementations", although you can find the implementation described on this site or in this pdf master thesis.

''' <summary>
''' The Bessel-Overhauser Spline interpolation
''' </summary>
''' <param name="p0">First point</param>
''' <param name="p1">Second point</param>
''' <param name="p2">Third point</param>
''' <param name="p3">Forth point</param>
''' <param name="t"></param>
''' <param name="u">The value from 0 - 1 where 0 is the start of the curve (globally) and 1 is the end of the curve (globally)</param>
''' <returns></returns>
''' <remarks></remarks>
Public Function PointOnBesselOverhauserCurve(ByVal p0 As System.Windows.Point, ByVal p1 As System.Windows.Point, ByVal p2 As System.Windows.Point, ByVal p3 As System.Windows.Point, ByVal u As Double, Optional ByVal t() As Double = Nothing) As System.Windows.Point
    Dim result As New System.Windows.Point()

    If t Is Nothing Then
        'Using these values will produce the Catmull-Rom Spline
        t = {1, 2, 3, 4}
    End If

    Dim ViXPlusHalf, VixMinusHalf, ViYPlusHalf, ViYMinusHalf, ViX, ViY As Double
    ViXPlusHalf = (p2.X - p1.X) / (t(2) - t(1))
    VixMinusHalf = (p1.X - p0.X) / (t(1) - t(0))
    ViYPlusHalf = (p2.Y - p1.Y) / (t(2) - t(1))
    ViYMinusHalf = (p1.Y - p0.Y) / (t(1) - t(0))

    ViX = ((t(2) - t(1)) * VixMinusHalf + (t(1) - t(0)) * ViXPlusHalf) / (t(2) - t(0))
    ViY = ((t(2) - t(1)) * ViYMinusHalf + (t(1) - t(0)) * ViYPlusHalf) / (t(2) - t(0))


    ' Location of Controlpoints
    Dim PointList As New PointCollection
    PointList.Add(p1)
    PointList.Add(New Point(p1.X + (1 / 3) * (t(2) - t(1)) * ViX, p1.Y + (1 / 3) * (t(2) - t(1)) * ViY))

    ViXPlusHalf = (p3.X - p2.X) / (t(3) - t(2))
    VixMinusHalf = (p2.X - p1.X) / (t(2) - t(1))
    ViYPlusHalf = (p3.Y - p2.Y) / (t(3) - t(2))
    ViYMinusHalf = (p2.Y - p1.Y) / (t(2) - t(1))

    ViX = ((t(3) - t(2)) * VixMinusHalf + (t(2) - t(1)) * ViXPlusHalf) / (t(3) - t(1))
    ViY = ((t(3) - t(2)) * ViYMinusHalf + (t(2) - t(1)) * ViYPlusHalf) / (t(3) - t(1))

    PointList.Add(New Point(p2.X - (1 / 3) * (t(3) - t(2)) * ViX, p2.Y - (1 / 3) * (t(3) - t(2)) * ViY))
    PointList.Add(p2)

    ' Get the calcualted value from the 3rd degree Bezier curve
    Return PointBezierFunction(PointList, u)
End Function

Lagrange Interpolation

This is one of the first properly known algorithms for drawing continuous curves , and this implementation uses code from "Mathematical tools in Computer graphics with C# implementations. It uses Neville's algorithm to calculate the basis function for a given t:

Private Function LagrangeBasisFunction(ByVal n As Integer, ByVal k As Integer, ByVal t As Double) As Double
    Dim l As Double = 1
    Dim tj, tk As Double
    tk = k / (n - 1)

    For j As Integer = 0 To n - 1
        If j <> k Then
            tj = j / (n - 1)
            l *= (t - tj) / (tk - tj)
        End If
    Next
    Return l
End Function

We must then use the same principle as the Bernstein polynomial to calculate the value at a given point on the curve:

Private Function GetLagrangianAtPointT(ByVal p As PointCollection, ByVal t As Double) As Point
    Dim n As Integer = p.Count

    Dim rx, ry As Double
    rx = 0
    ry = 0
    For i As Integer = 0 To n - 1
        rx += LagrangeBasisFunction(n, i, t) * p(i).X
        ry += LagrangeBasisFunction(n, i, t) * p(i).Y
    Next
    Return New Point(rx, ry)
End Function

Then the wrapper to translate the mathematics into a simple call that returns the curve:

Private Function LargrangianInterpolation(ByVal p As PointCollection, Optional ByVal StepSize As Double = 0.01) As PointCollection
    Dim result As New PointCollection

    For i As Double = 0 To 1 Step StepSize
        result.Add(GetLagrangianAtPointT(p, i))
    Next

    result.Add(p(p.Count - 1))
    Return result
End Function

Now that you have this, I must warn you that by using it you are quite confident that all you points are quite evenly spaced in both X and Y direction. If not, it might start to oscillate quite badly, so consider yourself warned.

Convex Hull

A convex hull could be explained quite easily by thinking of the two-dimensional points as strikes that firmly planted in the ground. The convex hull would be the rubber band that encloses all the points. The most commonly used computer algorithm for convex hull uses the cross product to evaluate which side a point is relative to a line.

But to use that algorithm most efficiently (not checking the same point twice) one need's to sort the 2D points lexicographically, that means first sort the points in the X direction, and if the X values are equal, the point will be further sorted on the Y value. To do this with an array of point's you first need a class that implements IComparer:

Private Class PointSort
      Implements IComparer
      Public Enum Mode
          X
          Y
      End Enum

      Private currentMode As Mode = Mode.X

      Public Sub New(ByVal mode As Mode)
          currentMode = mode
      End Sub

      'Comparing function
      'Returns one of three values - 0 (equal), 1 (greater than), -1 (less than)
      Private Function IComparer_Compare(ByVal a As Object, ByVal b As Object) As Integer Implements IComparer.Compare
          Dim point1 As Point = CType(a, Point)
          Dim point2 As Point = CType(b, Point)

          If currentMode = Mode.X Then
              If point1.X > point2.X Then
                  Return 1
              ElseIf point1.X < point2.X Then
                  Return -1
              Else
                  If point1.Y > point2.Y Then
                      Return 1
                  ElseIf point1.Y < point2.Y Then
                      Return -1
                  Else
                      Return 0
                  End If
              End If
          Else
              If point1.Y > point2.Y Then
                  Return 1
              ElseIf point1.Y < point2.Y Then
                  Return -1
              Else
                  If point1.X > point2.X Then
                      Return 1
                  ElseIf point1.X < point2.X Then
                      Return -1
                  Else
                      Return 0
                  End If
              End If
          End If
      End Function
  End Class

The IComparer is used with the Array.Sort method, and placed inside a wrapper function:

''' <summary>
''' Returns a sorted pointcollection based on Lexografically sorting
''' </summary>
''' <param name="samplepoints"></param>
''' <param name="SortByXdirection"></param>
''' <returns></returns>
''' <remarks></remarks>
Public Function SortPoints(ByVal samplepoints As PointCollection, ByVal SortByXdirection As Boolean) As PointCollection
    'Create another array so we can keep the original array out of order
    Dim copySamplePoints As Point() = New Point(samplepoints.Count - 1) {}
    samplepoints.CopyTo(copySamplePoints, 0)

    If SortByXdirection Then
        Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.X))
    Else
        Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.Y))
    End If

    Dim result As New PointCollection

    For Each p As Point In copySamplePoints
        result.Add(p)
    Next

    Return result
End Function

Checking which side a point lies relative to a line is implemented the following way:

''' <summary>
''' Finds which side of a line the point is
''' </summary>
''' <param name="PointToBeEvaluated">Evaluation point</param>
''' <param name="StartPointOnLine">Startpoint of line</param>
''' <param name="EndPointOnLine">Endpoint on line</param>
''' <returns>-1 for a point to the left, 0 for a point on the line, +1 for a point to the right</returns>
''' <remarks></remarks>
Private Function WhichSide(ByVal PointToBeEvaluated As Point, ByVal StartPointOnLine As Point, ByVal EndPointOnLine As Point) As Integer
    Dim ReturnvalueEquation As Double
    ReturnvalueEquation = ((PointToBeEvaluated.Y - StartPointOnLine.Y) _
                           * (EndPointOnLine.X - StartPointOnLine.X)) - ((EndPointOnLine.Y - StartPointOnLine.Y) _
                           * (PointToBeEvaluated.X - StartPointOnLine.X))

    If ReturnvalueEquation > 0 Then
        Return -1
    ElseIf ReturnvalueEquation = 0 Then
        Return 0
    Else
        Return 1
    End If

End Function

After creating all the basic functions needed, all that is left is to put it together. The code below should run on n*log n time, but I suspect that several improvements are possible, so that it could be even faster.

The code below does the following thing:
  1. Removes duplicate points
  2. Sorts the points
  3. Creates two arrays one sorts the points in descending order and the collection is called upper, and the other in ascending order and is called lower.
  4. Removes point in Upper if there is a point above it, and the same is done for lower
  5. Connect the two arrays again

The code is given below:

''' <summary>
''' Returns the polygon that uses the least amount of points to make the polygon that contains all the points
''' </summary>
''' <param name="points"></param>
''' <returns>PointCollection</returns>
''' <remarks>It removes duplicate points</remarks>
Public Function ConvexHull2D(ByVal points As PointCollection) As PointCollection
    Dim SortedList As New PointCollection
    Dim Upper As New PointCollection
    Dim Lower As New PointCollection
    Dim tempPoints As New PointCollection

    ' Add points only if they are unique
    For Each p As Point In points
        If Not tempPoints.Contains(p) Then tempPoints.Add(p)
    Next

    Dim res As New PointCollection

    'Sort the points lexografically
    SortedList = SortPoints(tempPoints, True)

    'Save an upper list for the top curve
    Upper = SortedList
    For i As Integer = SortedList.Count - 1 To 0 Step -1
        'Reverse order sorted list for the down side curve
        Lower.Add(SortedList(i))
    Next

    Dim j_set As Boolean = False
    Dim j As Integer = 0
    Dim r As Integer
    Do While j < Upper.Count - 2
        r = WhichSide(Upper(j + 1), Upper(j), Upper(j + 2))
        If r = -1 Or r = 0 Then
            Upper.RemoveAt(j + 1)
            j = 0
            j_set = True
        End If

        If Not j_set Then
            j += 1
        End If

        j_set = False
    Loop

    j = 0
    j_set = False
    Do While j < Lower.Count - 2
        r = WhichSide(Lower(j + 1), Lower(j), Lower(j + 2))
        If r = -1 Or r = 0 Then
            Lower.RemoveAt(j + 1)
            j = 0
            j_set = True
        End If

        If Not j_set Then
            j += 1
        End If

        j_set = False
    Loop

    'To connect the upper and lower curves stor them in a new variable
    For Each p As Point In Upper
        res.Add(p)
    Next

    For Each p As Point In Lower
        If Not res.Contains(p) Then
            res.Add(p)
        End If
    Next

    'Just to create a full circle, connect the first and last point
    res.Add(res(0))

    Return res
End Function

Epilogue

There is a lot of information that got lost in the implementation. It is not directly involved but it needs to be addressed.

First off it the term of continuity, as a curve which is piecewise interpolated have, but the question remains how many of the derivative are continuous. A definition of a cure is C0, as it just says that the curve is continuous . But a curve would be c1 continuous if it also satisfied that the first derivative is also conscious. C2 of course, means that the first and second derivative are both continuous in each connecting point.

There are also something called G1 continuity, which demands that the changes are constant, i.a. the derivative is strictly increasing. A typical variation of this requirement would be the clothoid spline which is used to design turns in the road, with an increasing inward acceleration and deceleration.

One function is missing, and that is the general implementation of the B-Spline curve, and it is planned to write that code and include it here as well later.

References

 

The bulk of information that is used in this article came from the first chapter in the book:

"The NURBS Book" 2nd Edition"

The code with Catmull-Rom and Bessel-Overhauser was described in the book below, and the Lagrange interpolation is just rewritten from it:

"Mathematical tools in computer graphics with C# implementation"

There are other sources as well, but they are mentioned in the main article.

"Harmonic interpolation for smooth curves and surfaces" by Alexandre Hardy (Ph.d thesis). Large parts of the previous book are taken from this thesis.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Kenneth Haugland
Engineer
Norway Norway
I hope that you like the stuff I have created and if you do wish to say thank you then a donation is always appreciated.
You can donate here[^].

You may also be interested in...

Comments and Discussions

 
GeneralMy vote of 5 Pin
John Schroedl5-Dec-14 2:28
memberJohn Schroedl5-Dec-14 2:28 
GeneralRe: My vote of 5 Pin
Kenneth Haugland5-Dec-14 2:58
professionalKenneth Haugland5-Dec-14 2:58 
GeneralExtraordinary! Pin
Member 103697863-Dec-14 22:39
memberMember 103697863-Dec-14 22:39 
GeneralRe: Extraordinary! Pin
Kenneth Haugland4-Dec-14 23:46
professionalKenneth Haugland4-Dec-14 23:46 
GeneralGood to see something of this quality Pin
loyal ginger3-Dec-14 2:56
memberloyal ginger3-Dec-14 2:56 
GeneralRe: Good to see something of this quality Pin
Kenneth Haugland4-Dec-14 23:41
professionalKenneth Haugland4-Dec-14 23:41 
QuestionSPline what ??? Pin
Teruteru31426-May-14 3:48
memberTeruteru31426-May-14 3:48 
AnswerRe: SPline what ??? Pin
Kenneth Haugland26-May-14 4:22
professionalKenneth Haugland26-May-14 4:22 
GeneralMy vote of 5 Pin
Renju Vinod14-May-14 21:21
professionalRenju Vinod14-May-14 21:21 
AnswerRe: My vote of 5 Pin
Kenneth Haugland14-May-14 21:39
professionalKenneth Haugland14-May-14 21:39 
QuestionGood Stuff! Pin
S Houghtelin12-May-14 1:36
professionalS Houghtelin12-May-14 1:36 
AnswerRe: Good Stuff! Pin
Kenneth Haugland12-May-14 6:19
professionalKenneth Haugland12-May-14 6:19 
Questionvery nice Pin
BillW339-May-14 9:28
professionalBillW339-May-14 9:28 
AnswerRe: very nice Pin
Kenneth Haugland9-May-14 20:20
professionalKenneth Haugland9-May-14 20:20 
QuestionThank you Pin
phil.o7-Apr-14 1:29
professionalphil.o7-Apr-14 1:29 
AnswerRe: Thank you Pin
Kenneth Haugland7-Apr-14 3:27
memberKenneth Haugland7-Apr-14 3:27 
QuestionVery cleanly done Pin
YvesDaoust6-Apr-14 23:32
memberYvesDaoust6-Apr-14 23:32 
AnswerRe: Very cleanly done Pin
Kenneth Haugland7-Apr-14 0:09
memberKenneth Haugland7-Apr-14 0:09 
QuestionWell Done Pin
rspercy656-Apr-14 1:51
memberrspercy656-Apr-14 1:51 
AnswerRe: Well Done Pin
Kenneth Haugland6-Apr-14 2:26
memberKenneth Haugland6-Apr-14 2:26 
QuestionNice Work Pin
sam.hill5-Apr-14 17:04
membersam.hill5-Apr-14 17:04 
AnswerRe: Nice Work Pin
Kenneth Haugland5-Apr-14 19:57
memberKenneth Haugland5-Apr-14 19:57 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.160525.2 | Last Updated 2 Dec 2014
Article Copyright 2014 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2016
Layout: fixed | fluid