立体样条插值的正确实现

4
我使用了其中一个提出的算法,但结果非常糟糕。
我在Java中实现了维基算法(代码如下)。x(0)points.get(0)y(0)values[points.get(0)]αalfaμmi。其余部分与维基伪代码相同。
    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);
        }
    }

}

我得到的结果如下: enter image description here 但应该与这个类似: enter image description here 同时,我也试图按照http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf中的方式以另一种方法实现算法。
首先他们展示了如何进行线性样条插值,这很容易。我创建了计算A和B系数的函数。然后他们通过添加二阶导数来扩展线性样条插值。C和D系数也很容易计算。
但是当我尝试计算二阶导数时,问题就开始了。我不知道他们是如何计算的。
因此,我只实现了线性插值。结果如下: enter image description here 有人知道如何修复第一个算法或者解释一下第二个算法中如何计算二阶导数吗?

在这里尝试一下:http://codereview.stackexchange.com/ - ceving
@ceving 他们告诉我,他们正在审查代码质量而不是代码结果。 - Mr Jedi
4
由于这里的结果不是期望的结果,所以这个问题不适合进行代码审查。 - Simon Forsberg
很高兴能在您的图表中看到控制点,以便了解真正的问题和非问题。 - Spektre
也许不应该使用被标记为混乱的过时维基百科文章。最近关于样条插值的文章令人惊讶地在spline interpolation上。 - Lutz Lehmann
4个回答

7

很抱歉,您的源代码对我来说真的是一团糟,我只能坚持理论。以下是一些提示:

  1. SPLINE cubics

    SPLINE is not interpolation but approximation to use them you do not need any derivation. If you have ten points: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9 then cubic spline starts/ends with triple point. If you create function to 'draw' SPLINE cubic curve patch then to assure continuity the call sequence will be like this:

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

    do not forget that SPLINE curve for p0,p1,p2,p3 draw only curve 'between' p1,p2 !!!

  2. BEZIER cubics

    4-point BEZIER coefficients can be computed like this:

        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. Interpolation

    to use interpolation you must use interpolation polynomials. There are many out there but I prefer to use cubics ... similar to BEZIER/SPLINE/NURBS... so

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

    The only thing left to do is compute a0,a1,a2,a3. You have 2 equations (p(t) and its derivation by t) and 4 points from the data set. You also must ensure the continuity ... So first derivation for join points must be the same for neighboring curves (t=0,t=1). This leads to system of 4 linear equations (use t=0 and t=1). If you derive it it will create an simple equation depended only on input point coordinates:

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

    the call sequence is the same as for SPLINE

[注意事项]

  1. the difference between interpolation and approximation is that:

    interpolation curve goes always through the control points but high order polynomials tend to oscillate and approximation just approaches to control points (in some cases can cross them but usually not).

  2. variables:

    • a0,... p0,... are vectors (number of dimensions must match the input points)
    • t is scalar
  3. to draw cubic from coefficients a0,..a3

    just do something like this:

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

2
可以使用插值三次样条来实现,在插值点处函数值、导数和二阶导数匹配。只需要解决一个大型但带状的线性方程组即可。 - Lutz Lehmann
为什么要这样反其道而行之呢?立方曲线具有c1连接限制(值和一阶导数匹配),如果需要更高的连通性,则使用更多点,5个点-> c2,6个点-> c3等等...使用4个点来实现c2连通性,那么它就不再是SPLINE了...甚至不是4点立方体,因为你也在解决附近的点...当然,我可能会错过一些最近发现的新方程,但我强烈怀疑。 - Spektre
2
我相信你错过了一些非常古老的公式。是的,如果给定值和导数,您可以构建一个分段三次函数。如果没有给出导数,则在选择它们时有很大的自由度。通过要求最小化二阶导数来减少这种自由度。实际上,这会导致具有线性端点的两次可微分的分段三次函数,该函数由三对角线性系统唯一确定。请参见spline interpolation - Lutz Lehmann
2
通常,“样条”一词是指那些在给定约束条件下具有最小曲率(或二阶导数)的分段三次函数。使用三次样条,C2是最好的选择,第三个导数几乎总是分段常数跳跃。如果没有跳跃,则第三个导数是常数,函数是一个三次多项式。 - Lutz Lehmann

7
最近由Unser、Thévenaz等人在一系列论文中描述了立方B样条,详见以下内容:
M. Unser, A. Aldroubi, M. Eden,“用于连续图像表示和插值的快速B样条变换”,IEEE Trans。Pattern Anal.Machine Intell.,vol。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,“插值再探讨”,IEEE Trans.on Medical Imaging,vol.19,no.7,pp.739-758,2000年7月。
以下是一些指南。
什么是样条?
样条是平滑连接在一起的分段多项式。对于度数为n的样条,每个片段都是一个n次多项式。这些部分相互连接,使得样条在节点上连续到它的导数为n-1,即多项式部分的连接点。
如何构建样条?
零阶样条如下:
所有其他样条都可以构造为:
其中卷积取n-1次。
三次样条
最流行的样条是三次样条,其表达式为:
样条插值问题
给定在离散整数点k处采样的函数f(x),样条插值问题是确定一个近似f(x)的逼近s(x),表示如下:
其中ck是插值系数,s(k)=f(k)。
样条预处理
不幸的是,从n=3开始:
因此,ck 不是插值系数。它们可以通过解决由强制 s(k) = f(k) 得到的线性方程组来确定,即,
这样的等式可以重构成卷积形式,并在变换的 z 空间中解决,如下所示:
其中
因此,
按这种方式进行处理总是比通过例如 LU 分解来解决线性方程组的解决方案更可取。
可以通过注意到
其中
第一个分数代表一个 因果滤波器,而第二个分数代表一个 反向因果滤波器。它们都在下面的图中说明。 因果滤波器 反向因果滤波器 在最后一个图中,
滤波器的输出可以用以下递归方程式表示
可以通过首先确定 c-c+ 的 "初始条件" 来解决上述方程式。假设周期为 镜像 输入序列 fk
则可以证明 c+ 的初始条件可以表示为
c+ 的初始条件可以表示为

1
请查看样条插值,尽管它们只提供一个可用的3x3示例。对于更多样本点,例如枚举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];

使用高斯-塞德尔迭代法可以解决这个问题(难道不是专门为这个系统发明的吗?)或者您喜欢的Krylov空间求解器。

0

// 文件:cbspline-test.cpp

// C++ Cubic Spline Interpolation using the Boost library.
// Interpolation is an estimation of a value within two known values in a sequence of values. 

// From a selected test function y(x) = (5.0/(x))*sin(5*x)+(x-6), 21 numbers of (x,y) equal-spaced sequential points were generated. 
// Using the known 21 (x,y) points, 1001 numbers of (x,y) equal-spaced cubic spline interpolated points were calculated.
// Since the test function y(x) is known exactly, the exact y-value for any x-value can be calculated.
// For each point, the interpolation error is calculated as the difference between the exact (x,y) value and the interpolated (x,y) value. 
// This program writes outputs as tab delimited results to both screen and named text file

// COMPILATION: $ g++ -lgmp -lm -std=c++11 -o cbspline-test.xx cbspline-test.cpp
// EXECUTION:   $ ./cbspline-test.xx 

// GNUPLOT 1:   gnuplot> load "cbspline-versus-exact-function-plots.gpl"
// GNUPLOT 2:   gnuplot> load "cbspline-interpolated-point-errors.gpl"

#include <iostream>
#include <fstream>
#include <boost/math/interpolators/cubic_b_spline.hpp>
using namespace std;

// ========================================================
int main(int argc, char* argv[]) {
// ========================================================

    // Vector size 21 = 20 space intervals, size 1001 = 1000 space intervals  
    std::vector<double> x_knot(21), x_exact(1001), x_cbspline(1001);
    std::vector<double> y_knot(21), y_exact(1001), y_cbspline(1001); 
    std::vector<double> y_diff(1001);

    double x;   double x0 = 0.0;    double t0 = 0.0;    
    double xstep1 = 0.5;    double xstep2 = 0.01;       

    ofstream outfile;                                   // Output data file tab delimited values
    outfile.open("cbspline-errors-1000-points.dat");    // y_diff = (y_exact - y_cbspline)

    // Handling zero-point infinity (1/x) when x = 0
    x0 = 1e-18; 
    x_knot[0] = x0; 
    y_knot[0] = (5.0/(x0))*sin(5*x0)+(x0-6);            // Selected test function

    for (int i = 1; i <= 20; i++) { // Note: Loop i begins with 1 not 0, 20 points 
        x = (i*xstep1); 
        x_knot[i] = x;
        y_knot[i] = (5.0/(x))*sin(5*x)+(x-6); 
    }

    // Execution of the cubic spline interpolation from Boost library
    // NOTE: Using xstep1 = 0.5 based on knot points stepsize (for 21 known points) 
    boost::math::cubic_b_spline<double> spline(y_knot.begin(), y_knot.end(), t0, xstep1);

    // Again handling zero-point infinity (1/x) when x = 0
    x_cbspline[0] = x_knot[0];  
    y_cbspline[0] = y_knot[0];

    for (int i = 1; i <= 1000; ++i) {       // 1000 points.
        x = (i*xstep2);  
        x_cbspline[i] = x;
        // NOTE: y = spline(x) is valid for any value of x in xrange [0:10] 
        // meaning, any x within range [x_knot.begin() and x_knot.end()]
        y_cbspline[i] = spline(x);   
    }

    // Write headers for screen display and output file
    std::cout << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  
    outfile   << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  

    // Again handling zero-point infinity (1/x) when x = 0
    x_exact[0] = x_knot[0]; 
    y_exact[0] = y_knot[0];
     y_diff[0] = (y_exact[0] - y_cbspline[0]);
    std::cout << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to screen
    outfile   << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to text file

    for (int i = 1; i <= 1000; ++i) {       // 1000 points
        x = (i*xstep2);  
        x_exact[i] = x;
        y_exact[i] = (5.0/(x))*sin(5*x)+(x-6); 
         y_diff[i] = (y_exact[i] - y_cbspline[i]);
        std::cout << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to screen
        outfile   << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to text file
    }
    outfile.close();

return(0);
}
// ========================================================
/*
# GNUPLOT script 1
# File: cbspline-versus-exact-function-plots.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-12.0:20.0]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 2.0
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y(x)"
set title "Function points y(x) = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:2 with lines linetype 3 linewidth 1.0 linecolor 'blue' title 'exact function curve', 'cbspline-errors-1000-points.dat' using 1:3 with lines linecolor 'red' linetype 3 linewidth 1.0 title 'cubic spline interpolated curve'

*/
// ========================================================
/*
# GNUPLOT script 2
# File: cbspline-interpolated-point-errors.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-2.50:2.5]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 0.5
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y"
set title "Function points y = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:4 with lines linetype 3 linewidth 1.0 linecolor 'red' title 'cubic spline interpolated errors'

*/
// ========================================================

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接