正确实施三次样条插值

我正在使用其中一种algorithm，但结果非常糟糕。

我实现了wikialgorithm

在Java中（代码如下）。 x(0)是points.get(0) ， y(0)是values[points.get(0)] ， α是alfa ， μ是mi 。 其余部分与wiki伪代码相同。

  public void createSpline(double[] values, ArrayList<Integer> points){ a = new double[points.size()+1]; for (int i=0; i <points.size();i++) { a[i] = values[points.get(i)]; } b = new double[points.size()]; d = new double[points.size()]; h = new double[points.size()]; for (int i=0; i<points.size()-1; i++){ h[i] = points.get(i+1) - points.get(i); } alfa = new double[points.size()]; for (int i=1; i <points.size()-1; i++){ alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]); } c = new double[points.size()+1]; l = new double[points.size()+1]; mi = new double[points.size()+1]; z = new double[points.size()+1]; l[0] = 1; mi[0] = z[0] = 0; for (int i =1; i<points.size()-1;i++){ l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1]; mi[i] = h[i]/l[i]; z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i]; } l[points.size()] = 1; z[points.size()] = c[points.size()] = 0; for (int j=points.size()-1; j >0; j--) { c[j] = z[j] - mi[j]*c[j-1]; b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ; d[j] = (c[j+1]-c[j])/((double)3*h[j]); } for (int i=0; i<points.size()-1;i++){ for (int j = points.get(i); j<points.get(i+1);j++){ // fk[j] = values[points.get(i)]; functionResult[j] = a[i] + b[i] * (j - points.get(i)) + c[i] * Math.pow((j - points.get(i)),2) + d[i] * Math.pow((j - points.get(i)),3); } } } 

我得到的结果如下：

但应该是这样的：

---

我也试图按照另一种方式实现algorithm： http : //www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

起初，他们展示了如何做线性样条，这很容易。 我创build了计算A和B系数的函数。 然后通过添加二阶导数来扩展线性样条。 C和D系数也很容易计算。

但是，当我试图计算二阶导数时，问题就开始了。 我不明白他们如何计算他们。

所以我只实现了线性插值。 结果是：

有谁知道如何解决第一个algorithm或解释我如何计算第二个algorithm中的二阶导数？

Answer 1

最近，Unser，Thévenaz 等人在一系列论文中描述了三次b样条。 ，除其他之外

M.User，A.Aldroubi，M.Eden，“用于连续图像表示和插值的快速B样条变换”， IEEE Trans。 模式肛门。 机器智能。 ，vol。 13，n。 3卷，277-285页，1991年3月。

M.Lunser，“Splines，一种完美适合信号和image processing”， IEEE Signal Proc。 MAG。 ，第22-38页，1999年11月。

和

P.Thévenaz，T. Blu，M. Unser，“Interpolation Revisited”， IEEE Trans。 在医学影像 ，第一卷。 19，没有。 7，pp.739-758，2000年7月。

这里有一些指导方针。

什么样条？

样条曲线是平滑连接在一起的分段多项式。 对于n样条，每个分段是n多项式。 这些片段被连接起来，使得样条在连接点上是连续的，直到n-1阶导数，即多项式片段的连接点。

如何构build样条？

零阶样条如下

所有其他样条可以被构造为

其中卷积是n-1次。

立方样条

最stream行的样条是三次样条，其expression式是

样条插值问题

给定在离散整数点k处采样的函数f(x) ，样条插值问题是确定以下列方式表示的近似s(x)到f(x)

其中ck是插值系数， s(k) = f(k) 。

样条预过滤

不幸的是，从n=3开始，

所以ck不是插值系数。 它们可以通过求解通过强制s(k) = f(k)得到的线性方程组来确定，

这样的方程可以以卷积forms进行重构，并在变换的z空间中求解

哪里

因此，

以这种方式进行总是优于通过例如LU分解提供线性方程组的解。

上述等式的解可以通过注意到来确定

哪里

第一部分代表一个因果filter ，而第二部分代表一个抗因filter 。 它们都在下面的图中显示。

因果filter

防止filter

在最后一幅图中，

滤波器的输出可以用下面的recursion方程来表示

上述方程可以通过首先确定c-和c+ “初始条件”来解决。 假定一个周期性的镜像input序列fk ，

那么可以certificatec+的初始条件可以表示为

而c+的初始条件可以表示为

Answer 2

对不起，你的源代码对我来说真的是一团糟，所以我坚持理论。 这里有一些提示：

1. SPLINE立方体

   SPLINE不是插值，而是近似使用它们，不需要任何派生。 如果你有十个点： p0,p1,p2,p3,p4,p5,p6,p7,p8,p9则三次样条曲线以三重点开始/结束。 如果你创build了“绘制” SPLINE三次曲线补丁的函数，那么为了保证连续性，调用顺序将如下所示：

     spline(p0,p0,p0,p1); spline(p0,p0,p1,p2); spline(p0,p1,p2,p3); spline(p1,p2,p3,p4); spline(p2,p3,p4,p5); spline(p3,p4,p5,p6); spline(p4,p5,p6,p7); spline(p5,p6,p7,p8); spline(p6,p7,p8,p9); spline(p7,p8,p9,p9); spline(p8,p9,p9,p9); 

   不要忘记， p0,p1,p2,p3 SPLINE曲线只绘制' p1,p2之间的曲线！

2. BEZIER立方体

   4点BEZIER系数可以这样计算：

     a0= ( p0); a1= (3.0*p1)-(3.0*p0); a2= (3.0*p2)-(6.0*p1)+(3.0*p0); a3=( p3)-(3.0*p2)+(3.0*p1)-( p0); 

3. 插值

   要使用插值，您必须使用插值多项式。 在那里有很多，但我更喜欢使用立方体…类似于BEZIER / SPLINE / NURBS …所以

   • p(t) = a0+a1*t+a2*t*t+a3*t*t*t其中t = <0,1>

   剩下唯一要做的就是计算a0,a1,a2,a3 。 你有2个方程（ p(t)及其由t导出）和4个点的数据集。 你也必须保证连续性…因此，对于相邻曲线（ t=0,t=1 ），连接点的一阶导数必须相同。 这导致了4个线性方程组（使用t=0和t=1 ）。 如果你得到它，它将创build一个简单的方程只依赖于input点坐标：

     double d1,d2; d1=0.5*(p2-p0); d2=0.5*(p3-p1); a0=p1; a1=d1; a2=(3.0*(p2-p1))-(2.0*d1)-d2; a3=d1+d2+(2.0*(-p2+p1)); 

   呼叫顺序与SPLINE相同

[笔记]

1. 插值和逼近的区别在于：

   插值曲线始终通过控制点，但高阶多项式倾向于振荡，近似刚好接近控制点（在某些情况下可以跨越它们，但通常不会）。

2. variables：

   • a0,... p0,...是vector（维数必须与input点匹配）
   • t是标量

3. 从系数a0,..a3画出立方

   只是做这样的事情：

     MoveTo(a0.x,a0.y); for (t=0.0;t<=1.0;t+=0.01) { tt=t*t; ttt=tt*t; p=a0+(a1*t)+(a2*tt)+(a3*ttt); LineTo(px,py); } 

Answer 3

尽pipe它们只给出了一个可用的3×3示例，但请参见样条插值 。 对于更多的采样点，说N + 1枚举x[0..N]的值为y[0..N]必须解决以下系统的未知k[0..N]

  2* c[0] * k[0] + c[0] * k[1] == 3* d[0]; c[0] * k[0] + 2*(c[0]+c[1]) * k[1] + c[1] * k[2] == 3*(d[0]+d[1]); c[1] * k[1] + 2*(c[1]+c[2]) * k[2] + c[2] * k[3] == 3*(d[1]+d[2]); c[2] * k[2] + 2*(c[2]+c[3]) * k[3] + c[3] * k[4] == 3*(d[2]+d[3]); ... c[N-2]*k[N-2]+2*(c[N-2]+c[N-1])*k[N-1]+c[N-1]*k[N] == 3*(d[N-2]+d[N-1]); c[N-2]*k[N-1] + 2*c[N-1] * k[N] == 3* d[N-1]; 

哪里

 c[k]=1/(x[k+1]-x[k]); d[k]=(y[k+1]-y[k])*c[k]*c[k]; 

这可以使用Gauß-Seidel迭代来解决（这不是为这个系统准确发明的），还是你最喜欢的Krylov空间解算器。