样条插值 - 历史、理论与实现

Kenneth Haugland

4.99/5 (37投票s)

2014年4月5日

CPOL

14分钟阅读

71947

2891

贝塞尔曲线、贝塞尔曲线导数、Catmull-Rom样条、Bessel-Overhauser样条、拉格朗日插值和凸包的实现

引言

网上有很多基于贝塞尔曲线的插值方案的实现，但它们似乎要么专门针对一种曲线，要么其函数不够通用，无法满足多种曲线的需求。

内置的贝塞尔曲线只有二次和三次曲线，对于简单应用来说应该足够了，但在许多情况下这根本不够。本文的目的是提供一个包含可重用函数的程序，您可以轻松地将其导入您的程序中。

背景

假设您想创建一个通过一些点的平滑曲线，而您还 stuck 在20世纪50年代；您需要做的是在想要保持静止的点上放置样条，获取所需硬度的样条，然后像波音公司的这位先生一样爬到桌子上。

正如微软在其MSDN文章中所写：“物理样条是一根细木条或其他柔性材料”，然后就有了样条这个词，后来由Isaac Jacob Schoenberg命名。有关发展的更多历史可以在这里阅读，基本情况是它在很多地方都有使用，但直到50年代和60年代数学理论才开始流行起来，我引用一下：

用于汽车车身建模的样条似乎有几个独立的起源。雪铁龙公司的de Casteljau、雷诺公司的Bezier、以及通用汽车公司的Birkhoff、Garabedian和de Boor都声称拥有各自的功劳，这些工作都发生在20世纪60年代初或50年代末。de Casteljau至少有一篇论文发表于1959年，但未广泛流传。de Boor在通用汽车公司的工作导致了20世纪60年代初发表了多篇论文，包括一些关于B样条的基础性工作。

贝塞尔曲线

构建贝塞尔曲线最常见的方法是使用De Casteljau算法，它是一种生成贝塞尔曲线的线性插值。然而，贝塞尔曲线也可以通过另一种方式定义，即使用伯恩斯坦多项式。

伯恩斯坦多项式

伯恩斯坦多项式与卡尔·魏尔斯特拉斯（Karl Weierstrass）于1885年发表的一项定理有关（他当时70岁，所以“数学创新只由年轻人（20岁或以下）完成”的理论可能应该被抛弃了？）。他基本上发现“区间[0,1]上的任何连续函数都可以被多项式任意精确地均匀逼近。”。伯恩斯坦多项式只是魏尔斯特拉斯定理的一种不同证明，可以通过标准方法证明；使用二项式函数

\begin{equation} B_{i,n}(t) = {n \choose i} t^i (1-t)^{n-i} \end{equation}

然而，如果您学过统计学，您会（或应该）注意到这个公式与二项分布完全相同。起初，这个公式可能看起来有些奇怪和陌生，但实际上它非常简单。一个人首先需要理解的短语是${N \choose k}$，展开来说就是

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

用通俗的话来说，它只是解释了如果您有N个元素并想从中选出k个元素，那么恰好有${N \choose k} $种方法。上述公式实际上包含两个不同的部分，第一部分解释了从N个元素中选出k个元素的方法数：$\frac{N!}{(N-k)!}$，第二部分$\frac{1}{k!}$只是消除了排列数，即组合k个元素的不同方法数。

为什么在选择元素时会用到阶乘实际上非常巧妙。如果您有一副52张牌，您想从中选出4张牌，我可以将每次从牌堆中抽取牌的统计分解如下：

第一次选择牌时，有52张牌可供选择，所以我有52个等可能的机会。
下一次我想从牌堆中取牌时，现在只剩下51张牌可供选择。因此，恰好拿到前两张牌的机会数是$52 \cdot 51$。
第三次抽取，有50张牌可供选择，组合数是$52 \cdot 51 \cdot 50$。
第四次，概率是$52 \cdot 51 \cdot 50 \cdot 49$。

写出来会很繁琐，但幸运的是，有一种简单的通用方法可以生成这个序列，只需写出

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

这将给出按明确顺序组合4张牌的方法数，但通常您不关心抽取的事件顺序，只关心最终的4张牌。可以通过这个简单的类比来说明；如果您需要选择2张牌，它们可以以序列{1,2}或{2,1}的形式出现，这两种方式在数学上称为排列。在获得您想要的牌的情况下，组合序列的方法数是2!，在我们有4张牌的情况下是4!。我们想消除获得k的方式中的排列数，所以我们修改公式

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

组合它们的方法数只是计算获得所需结果概率的一半。您还需要确定获得所需结果与未获得结果的概率。使用伯恩斯坦多项式最简单的方法是选择2个点并将它们组合起来。有一个要求，即对于所有t，期望结果的概率加上非期望结果的概率必须等于1。这只是一种实现方式，因此公式变为

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

我们通常不喜欢用计算机处理阶乘，因为数字很快就会变得非常大，并且非常容易出现舍入误差。公式而是通过使用伯恩斯坦多项式的递推关系来实现的

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

这可以翻译成如下代码

    ''' <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

该代码将生成给定输入点数所需的所有伯恩斯坦多项式。曲线如下所示，范围从输入函数的1个点到6个点（图片来自MathWorld）

实现起来相当简单，假设您已经获得了所有伯恩斯坦多项式。它只是伯恩斯坦多项式与其对应的点的总和

    ''' <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

贝塞尔曲线的导数

伯恩斯坦多项式的导数可用于生成任何贝塞尔曲线的切线，正如在Silverlight上的文章此处所做的那样。我将做同样的事情，只做一些小的改动。

请记住，只有伯恩斯坦多项式才依赖于u的变化。这意味着节点以及起点和终点都被视为常数。有几篇关于如何获得贝塞尔曲线导数的文章，但很少提到可以使用伯恩斯坦多项式来获得导数。我们感兴趣的关系是

B'_i,n(u) )= n ( B_i-1,n-1(u) - B_i,n-1(u))

当然，必须记住，伯恩斯坦多项式被定义为：如果i小于零，或者i大于n，则为零。用数学（简洁地说）术语来说

B_-1,n-1(u) = B_n,n-1(u) = 0

那么伯恩斯坦多项式导数的代码就可以轻松构建

  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

由于只有伯恩斯坦多项式与控制点一起提供了曲线的线条，因此可以使用级数展开（或幂律）逐段积分该多项式。现在我们只需使用导数来计算斜率

    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

您现在需要做的就是用y'/x'进行除法，您将得到给定点的斜率。由于这是一种解析方法，因此它在贝塞尔曲线从0到1的每个步点都有效。

所以，我决定从文章“Silverlight中的贝塞尔曲线角度计算”中借用箭头图片，并在我的实现中使用它。但是，有一个问题，旋转不应该围绕中心进行，或者说可以，但更容易围绕0,0处的给定点旋转，然后将该点移动一半长度/宽度并进行90度转弯。这确保了箭头总是在正确的位置，并且代码可以轻松使用三角学：

    ''' <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

如果您看到评论部分，其中一条指出可以使用Tan2方法来计算角度。它接受两个点，并计算与X轴上的直线相关的角度。

    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

放置图像和旋转图像的代码也可以简化

        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的贝塞尔方法

现在，也有可能使用WPF中预定义的 Cubic Bezier 段绘制方法。MSDN在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>

对于 Quadratic Bezier 段也有不同的代码（同样来自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>

如果您选择使用伯恩斯坦多项式而不是这些，那么这两种版本都包含在内，但您也可以使用上面的代码在XAML中设置动画。

有理贝塞尔曲线

如果单独使用伯恩斯坦多项式，所有点都具有同等的重要性。另一方面，如果您想为每个点分配不同的权重（类似于太阳系的质量不同），您将获得对贝塞尔曲线形状的更多控制。正确使用甚至可以用三个点复制一个圆。

因此，不是直接使用伯恩斯坦多项式，而是将其乘以一个称为w_i的权重。但是，为了确保曲线保持在边界内，所有u的总和等于1非常重要。这简称为有理基函数R_i,n(u)并定义为

    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插值

我知道程序员可能有点懒惰，所以他们可能会看到需要Catmull-Rom插值，搜索一下，然后找到一个。就像这个一样，然后获取代码。

    ''' <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

还有另一种方法，那就是尝试理解Catmull-Rom样条到底是什么。上面代码中使用的多项式实际上是从下图的导数生成的。但它只使用导数来定义点1和点2之间的两个控制点。

导数当然按照惯例定义

而点P_i⁺和p_i+1^-是由导数定义的控制点

我们现在可以在标记为Pi和Pi+1的点之间计算两个新点，关键是实际使用3次贝塞尔曲线，使用起点Pi和终点Pi+1，以及两个节点Pi+和Pi+1-。现在我们明白Catmull-Rom样条只是3次贝塞尔曲线的一种特殊形式。解释这一点以及图形来源的有用链接是这里和这里。

    ''' <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

我遇到的问题是一条单直线，我经常需要处理它，那么我应该如何插值呢？我选择使用这样的逻辑：明天的天气和今天的天气一样，或者说半天。我在原始线的前面和后面各加了一个点，然后我说有一个点与前两个点之间的斜率相同，与最后两个点之间的斜率相同，只是距离减半。只要你一开始或结尾没有完全疯掉，它似乎都能很好地拟合。

唯一需要构建的是一个包装器，它可以接受标准的PointCollection作为Catmull-Rom样条的输入。

    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

如果控制点分布非常不均匀，Catmull-Rom样条很容易产生相当糟糕的曲线。如果您看一下当您取4个点，将两个中心点放在非常近的位置，而起点和终点与两个中心点有很大的距离时，曲线会变成什么样。您会在两个中心点之间得到一个急剧的弯曲。

它使用Bessel切线方法或Overhauser方法通过非均匀权重t_i来生成样条。顺便说一句，如果您使用均匀权重，它将产生Catmull-Rom样条。t_i的值旨在减小“速度”，以便出现更平滑的曲线。

我没有找到该曲线的任何实现，因此它是根据《计算机图形学中的数学工具及C#实现》一书中的算法创建的，尽管您可以在此网站或这篇PDF硕士论文中找到实现。

    ''' <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

拉格朗日插值

这是最早为人所知的绘制连续曲线的算法之一，此实现使用了《计算机图形学中的数学工具及C#实现》中的代码。它使用Neville算法来计算给定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

然后，我们必须使用与伯恩斯坦多项式相同的原理来计算曲线上给定点的值。

    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

然后是包装器，将数学转换为一个简单的调用，返回曲线。

    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

现在您有了这个，我必须警告您，使用它意味着您很确定您所有的点在X和Y方向上都相当均匀地分布。否则，它可能会开始剧烈振荡，所以请做好心理准备。

凸包

可以很容易地解释凸包，想象一下二维点就像钉在地上的桩子。凸包就是包围所有点的橡皮筋。最常用的凸包计算机算法使用叉积来评估一个点相对于一条线的哪一侧。

但要最高效地使用该算法（避免重复检查同一点），您需要按字典序对二维点进行排序，这意味着首先按X方向排序点，如果X值相等，则按Y值进一步排序。要做到这一点，您需要一个实现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

IComparer与Array.Sort方法一起使用，并放在一个包装函数中。

    ''' <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

检查一个点相对于一条线的哪一侧的实现方式如下：

   ''' <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

创建所有基本函数后，剩下要做的就是将它们组合起来。下面的代码应该以n*log n的时间运行，但我怀疑可以进行多项改进，使其速度更快。

下面的代码执行以下操作：

移除重复的点

对点进行排序

创建两个数组，一个按降序排序点，集合称为upper，另一个按升序排序，称为lower。

如果有一个点在它上面，则从Upper中移除该点，对Lower也执行相同的操作。

重新连接两个数组

代码如下

    ''' <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

结语

实现中丢失了很多信息。它不直接涉及，但需要解决。

首先是连续性的概念，因为分段插值曲线具有连续性，但问题仍然是导数的连续性有多少。曲线的定义是C0，因为它只是说曲线是连续的。但如果一条曲线也满足其一阶导数也是连续的，那么它就是c1连续的。C2当然意味着一阶和二阶导数在每个连接点都是连续的。

还有一个叫做G1连续性，它要求变化是恒定的，即导数严格递增。这种要求的典型变化是Clothoid样条，它用于设计道路的转弯，具有递增的向内加速和减速。

缺少一个函数，那就是B样条曲线的通用实现，并且计划稍后编写该代码并包含在此处。

参考文献

本文使用的大部分信息来自该书的第一章。

"NURBS Book" 第二版"

Catmull-Rom和Bessel-Overhauser的代码在下面的书中有所描述，而拉格朗日插值只是重写了它。

"计算机图形学中的数学工具及C#实现"

还有其他来源，但它们在主文章中已经提到。

"谐波插值用于光滑曲线和曲面" ，作者：Alexandre Hardy（博士论文）。前一本书的大部分内容摘自这篇论文。