C# 立方样条插值

• 下载源代码

引言

这是基于维基百科文章 样条插值 和 三对角矩阵算法 的三次样条插值实现。我创建此实现的目的是提供一个简单、清晰的实现，它与维基百科文章中的公式紧密匹配，而不是一个优化的实现。

背景

样条算法是一种用一组连接曲线（每条曲线称为一个样条）拟合数据点的方法，这样可以计算（插值）数据点之间的值。可以使用各种类型/阶数的方程来指定样条，包括线性、二次、三次等。

N 个点需要 N-1 个样条来连接它们。每个样条都由一个方程（例如，一个多项式）描述。这些多项式中的系数最初是未知的，样条算法会计算它们。一旦系数已知，要计算任何特定 X 的 Y，您需要找到与该 X 相关的样条，并使用这些系数来评估该 X 处的多项式。

为了计算样条的系数，该算法使用创建方程组然后求解该方程组的标准线性代数过程。这些方程体现了对样条的约束，即它们与要拟合的点相交，并且在样条之间的每个连接处，前一个样条的一阶和二阶导数等于后一个样条的一阶和二阶导数。这还不足以完全确定系统，因此又添加了两个约束，规定第一个和最后一个点的二阶导数为零。这些最终约束被称为“自然”样条样式。或者，用户 Quience（如下）提供了代码，以便调用者可以明确指定第一个和最后一个点的斜率。有了任何一组附加约束，系统就变得可解。

幸运的是，方程组被限制为三对角线，因此使用托马斯算法求解系统是 O(N)，而不是例如高斯消元的 O(N^3)。这是与上图对应的三对角矩阵，它用于解决维基百科文章的示例问题

2.0000, 1.0000, 0.0000
1.0000, 2.6667, 0.3333
0.0000, 0.3333, 0.6667

这个例子有点小，不能很好地说明三对角矩阵，但三对角线意味着只有主对角线（`row == col`）元素、主对角线下方（`row == col + 1`）的对角线和主对角线上方（`row == col - 1`）的对角线可以具有非零值。

链接的维基百科文章很好地解释了样条，因此我不会在这里提供任何冗余的解释。

Using the Code

我严格遵循了维基百科文章的实现，尽可能使用了变量名。有些注释引用了文章中的方程编号（不幸的是，当文章更改时，这些编号有时会失效）。

主要的类是 `CubicSpline` 和 `TriDiagonalMatrix`。`CubicSpline` 类使用 `TriDiagonalMatrix` 来求解方程组。通常，如果您只是使用此代码来拟合曲线，您将只使用 `CubicSpline`。

在源代码中，还有一个 *Program.cs*，其中包含一个 `Main()` 方法，该方法运行几个函数来练习这些类。

这是 *Program.cs* 中 `TestSplineOnWikipediaExample()` 方法的 `CubicSpline` 的一个示例用法。这部分代码设置了要拟合的测试数据

    // Create the test points.
    float[] x = new float[] { -1.0f, 0.0f, 3.0f };
    float[] y = new float[] { 0.5f, 0.0f, 3.0f };

    Console.WriteLine("x: {0}", ArrayUtil.ToString(x));
    Console.WriteLine("y: {0}", ArrayUtil.ToString(y));

下一步是创建一个 `float[]`（这里称为 `xs`），其中包含您想要确定 Y 值的 X 值。换句话说，这些是您想要在图表上绘制点的 X 值

    // Create the upsampled X values to interpolate
    int n = 20;
    float[] xs = new float[n];
    float stepSize = (x[x.Length - 1] - x[0]) / (n - 1);

    for (int i = 0; i < n; i++)
    {
        xs[i] = x[0] + i * stepSize;
    }

接下来是使用 `CubicSpline` 拟合数据，然后使用上面生成的 xs 值评估拟合的样条

    // Solve
    CubicSpline spline = new CubicSpline();
    float[] ys = spline.FitAndEval(x, y, xs);

    Console.WriteLine("xs: {0}", ArrayUtil.ToString(xs));
    Console.WriteLine("ys: {0}", ArrayUtil.ToString(ys));

    // Plot
    string path = @"..\..\spline-wikipedia.png";
    PlotSplineSolution("Cubic Spline Interpolation - Wikipedia Example", x, y, xs, ys, path);

这将显示以下控制台输出

x: -1, 0, 3
y: 0.5, 0, 3
xs: -1, -0.7894737, -0.5789474, -0.368421, -0.1578947, 0.05263158, ...
ys: 0.5, 0.3570127, 0.2245225, 0.1130267, 0.0330223, -0.005029888, ...

`PlotSplineSolution()` 方法使用 .NET Framework 的 `System.Windows.Forms.DataVisualization.Charting` 类绘制输入点和拟合插值样条曲线的图表。由此创建的图像是本文中包含的图像（上图）。

实现

样条拟合函数的核心是设置三对角矩阵，然后使用它来求解方程组。三对角矩阵不表示为矩阵，而是三个一维数组 A、B 和 C。数组 A 是下对角线，B 是主对角线，C 是上对角线，以匹配维基百科文章的名称。一旦设置好三对角矩阵，样条拟合函数就会调用 `TriDiagonalMatrix.Solve()`。

这是 *TriDiagonalMatrix.cs* 中的三对角矩阵求解函数

// Solve the system of equations this*x=d given the specified d.
public float[] Solve(float[] d)
{
    int n = this.N;

    if (d.Length != n)
    {
        throw new ArgumentException("The input d is not the same size as this matrix.");
    }

    // cPrime
    float[] cPrime = new float[n];
    cPrime[0] = C[0] / B[0];

    for (int i = 1; i < n; i++)
    {
        cPrime[i] = C[i] / (B[i] - cPrime[i-1] * A[i]);
    }

    // dPrime
    float[] dPrime = new float[n];
    dPrime[0] = d[0] / B[0];

    for (int i = 1; i < n; i++)
    {
        dPrime[i] = (d[i] - dPrime[i-1]*A[i]) / (B[i] - cPrime[i - 1] * A[i]);
    }

    // Back substitution
    float[] x = new float[n];
    x[n - 1] = dPrime[n - 1];

    for (int i = n-2; i >= 0; i--)
    {
        x[i] = dPrime[i] - cPrime[i] * x[i + 1];
    }

    return x;
}

以下是 `CubicSpline.Eval()` 方法的核心。此方法接受一个 `float[]` x，其中包含您想要使用拟合的样条计算 `y` 的 x 值。循环正在遍历 `x` 的每个值，计算相应的 `y` 值。

for (int i = 0; i < n; i++)
{
    // Find which spline can be used to compute this x
    int j = GetNextXIndex(x[i]);

    // Evaluate using j'th spline
    float t = (x[i] - xOrig[j]) / (xOrig[j + 1] - xOrig[j]);
    // equation 9 in the Wiki
    y[i] = (1 - t) * yOrig[j] + t * yOrig[j + 1] + t * (1 - t) * (a[j] * (1 - t) + b[j] * t); 
}

另一个例子

这是另一个图表。这个图表是由 *Program.cs* 中的 `TestSpline()` 方法创建的

计算斜率

除了计算每个 `X` 的 `Y` 之外，还可以计算样条的斜率。它在维基百科文章中表示为 q'，公式为方程 5。计算斜率的方法是 CubicSpline.EvalSlope()，它返回一个 `float` 数组（为您提供的每个 X 返回一个斜率值）。您必须在调用此方法之前调用以下方法之一Fit()（或 FitEndEval()）。这是一个示例用法

    // Try slope (spline is already computed at this point, see above code example)
    float[] slope = spline.EvalSlope(xs); // Same xs as first example above
    string slopePath = @"..\..\spline-wikipedia-slope.png";
    PlotSplineSolution("Cubic Spline Interpolation - Wikipedia Example - Slope", 
                       x, y, xs, ys, slopePath, slope);

这是结果图表

输入约束和参数拟合

正常的三次样条算法适用于二维点，其中 `y` 是 `x` 的函数，即 `y=f(x)`，并且每个 `x` 都有一个单独的 `y` 值。然而，下面的评论中用户 LutzL 指出了一种巧妙的方法，可以使用样条来拟合不符合此定义的点序列

[你]发明了一个单调递增的第三个时间坐标，T=0,1,2,3，或者从 T=0 开始，并按当前点到下一个点的距离递增。然后使用 (T,X) 和 (T,Y) 对来计算两个三次样条 x(t) 和 y(t)，并将曲线 (x(t),y(t)) 作为结果绘制出来。

我在 CubicSpline 类中将其实现为一个新的 `static` 方法 FitParametric()。这是一个示例用法

    float[] x = { 0.5f, 2.0f, 3.0f, 4.5f, 3.0f, 2.0f };
    float[] y = { 4.0f, 2.0f, 6.0f, 4.0f, 3.0f, 5.0f };
    float[] xs, ys;
    CubicSpline.FitParametric(x, y, 100, out xs, out ys);

这是结果图表

感谢用户 YvesDaoust 为此方法提供了正确的术语。

关注点

令人惊讶的是，在我基准测试的场景中，这个简单的实现比 C++ 中非常简洁的 Numerical Recipes（Press 等人，剑桥大学出版社，2002）版本的 C# 移植版本快两倍多。我只基准测试了一个场景，并没有对结果进行太多调查，所以不应该得出太多结论，但无论如何这很有趣。我使用 Visual Studio 2012 的 Release 构建，在 Win 7 Core i5 2500k 上运行，对实现进行了基准测试。我使用了 10,000 个随机点并插值了 100,000 个点，重复了 100 次。请参阅 *Program.cs* 方法 `TestPerf()`。由于许可和版权问题，我没有在代码中包含 C++ 版本的 Numerical Recipes。[更新：这种性能差异很可能是由于 NR in C 实现对每个评估都进行适当样条的搜索，而我要求评估点已排序，因此我可以进行同步遍历。]

`System.Windows.Forms.DataVisualization.Charting` 命名空间中的 `Chart` 类非常适合我的简单图表。其中的类使渲染图表并保存到文件变得非常容易。

不幸的是，C#/.NET 泛型（仍然）无法处理创建允许您使用相同代码同时处理 `float` 和 `double` 的泛型类型约束。我使用 `float` 实现了这一点，因为我通常从 `float` 开始，直到证明我需要 `double`。

历史

• 2013年3月10日 - 第一个版本
• 2013年9月20日 - 添加了第一个和最后一个点约束的描述，包括显式斜率参数，以及其他一些小改动
• 2014年7月23日 - 更新了代码以添加斜率计算，并在文章中添加了“计算斜率”部分
• 2014年7月25日 - 添加了“输入约束和几何拟合”部分，并更新了代码
• 2016年4月1日 - 将“几何”重命名为“参数”，并为参数拟合添加了起点和终点斜率规范