正确实现三次样条插值

Question

我正在使用其中一种提议的算法,但结果非常糟糕.

我实现了wiki算法

在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);
        }
    }

}

我得到的结果如下:

但它应该类似于:

---

我还试图以另一种方式实现该算法:http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

首先,他们展示了如何做线性样条,这很容易.我创建了计算A和B系数的函数.然后,他们通过添加二阶导数来扩展线性样条.C和D系数是容易计算了.

但是当我试图计算二阶导数时,问题就出现了.我不明白他们是如何计算的.

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

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

Answer 1

最近,Unser,Thévenaz 等人在一系列论文中描述了立方b样条.,见其他人

M. Unser,A.Aldroubi,M.Eden,"用于连续图像表示和插值的快速B样条变换",IEEE Trans.模式肛门.机器智能.,第一卷 13,n.3,pp.277-285,1991年3月.

M. Unser,"样条,完美适合信号和图像处理",IEEE Signal Proc.MAG.,pp.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在结,即,多项式片的接合点.

如何构造样条线？

零阶样条曲线如下

所有其他样条曲线都可以构造为

卷积的n-1时间.

立方样条

最流行的样条是三次样条,其表达式为

样条插值问题

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

其中ck的是插值系数和s(k) = f(k).

样条预过滤

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

所以ck's不是插值系数.它们可以通过求解通过强迫获得的线性方程组来确定s(k) = f(k),即

这样的方程可以以卷积形式重铸,并在变换z空间中求解

哪里

因此,

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

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

哪里

第一部分代表因果过滤器,而第二部分代表反义过滤器.它们都在下图中说明.

因果过滤器

Anticausal过滤器

在上图中,

滤波器的输出可以用以下递归方程表示

上述等式可以通过首先确定"初始条件"来解决c-和c+.上假定周期性的,镜像的输入序列fk,使得

然后可以证明初始条件c+可以表示为

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

Answer 2

对不起,但你的源代码对我来说真的是一个难以理解的混乱,所以我坚持理论.以下是一些提示:

1. SPLINE立方体

   SPLINE不是插值,而是使用它们的近似值,您不需要任何推导.如果你有十个点:p0,p1,p2,p3,p4,p5,p6,p7,p8,p9则三次样条曲线以三重点开始/结束.如果你创建了'绘制' 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);
   

   不要忘记SPLINE曲线p0,p1,p2,p3只绘制曲线'之间' 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).如果你派生它,它将创建一个仅取决于输入点坐标的简单方程:

       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. 变量:

   • a0,... p0,... 是向量(维度的数量必须与输入点匹配)
   • 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(p.x,p.y);
        }