正确实施三次样条插值

我正在使用其中一种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了计算AB系数的函数。 然后通过添加二阶导数来扩展线性样条。 CD系数也很容易计算。

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

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

在这里输入图像描述

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

最近,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+的初始条件可以表示为

在这里输入图像描述

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

  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=0t=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); } 

尽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空间解算器。