Java 三次样条插值的正确实现
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/20303172/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Proper implementation of cubic spline interpolation
提问by Mr Jedi
I was using one of the proposed algorithms out there but the results are very bad.
我正在使用其中一种建议的算法，但结果非常糟糕。
I implemented the wiki algorithm
我实现了维基算法
in Java (code below). x(0)is points.get(0), y(0)is values[points.get(0)], αis alfaand μis mi. The rest is the same as in the wiki pseudocode.
在 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);
}
}
}
The result that I get is the following:
我得到的结果如下：
but it should be similar to this:
但它应该类似于：
I'm also trying to implement the algorithm in another way according to: http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf
我还尝试根据以下内容以另一种方式实现该算法：http: //www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf
At first they show how to do linear spline and it's pretty easy. I create functions that calculate Aand Bcoefficients. Then they extend linear spline by adding second derivative. Cand Dcoefficients are easy to calculate too.
起初他们展示了如何做线性样条，这很容易。我创建了计算A和B系数的函数。然后他们通过添加二阶导数来扩展线性样条。C并且D系数也很容易计算。
But the problems starts when I attempt to calculate the second derivative. I do not understand how they calculate them.
但是当我尝试计算二阶导数时问题就开始了。我不明白他们是如何计算的。
So I implemented only linear interpolation. The result is:
所以我只实现了线性插值。结果是：
Does anyone know how to fix the first algoritm or explain me how to calculate the second derivative in the second algorithm?
有谁知道如何修复第一个算法或向我解释如何计算第二个算法中的二阶导数？
回答by Spektre
Sorry but Your source code is really a unreadable mess to me so I stick to theory. Here are some hints:
抱歉，你的源代码对我来说真的是一团糟，所以我坚持理论。这里有一些提示：
SPLINE cubics
SPLINE is not interpolation but approximationto use them you do not need any derivation. If you have ten points:
p0,p1,p2,p3,p4,p5,p6,p7,p8,p9then cubic spline starts/ends with triple point. If you create function to 'draw' SPLINEcubic 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,p3draw only curve 'between'p1,p2!!!BEZIER cubics
4-point BEZIERcoefficients 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);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*twheret = <0,1>
The only thing left to do is compute
a0,a1,a2,a3. You have 2 equations (p(t)and its derivation byt) 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 (uset=0andt=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
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！！！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);插值
要使用插值，您必须使用插值多项式。那里有很多，但我更喜欢使用立方体......类似于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相同
[Notes]
[笔记]
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).
variables:
a0,... p0,...are vectors (number of dimensions must match the input points)tis scalar
to draw cubic from coefficients
a0,..a3just 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); }
插值和近似的区别在于：
插值曲线总是通过控制点，但高阶多项式往往会振荡，逼近只是接近控制点（在某些情况下可以穿过它们但通常不会）。
变量：
a0,... p0,...是向量（维数必须与输入点匹配）t是标量
从系数中绘制三次方
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); }
回答by Lutz Lehmann
See spline interpolation, although they give only a usable 3x3 example. For more sample points, say N+1 enumerated x[0..N]with values y[0..N]one has to solve the following system for the unknown k[0..N]
请参阅样条插值，尽管它们仅给出了一个可用的 3x3 示例。对于更多的样本点，假设 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];
where
在哪里
c[k]=1/(x[k+1]-x[k]); d[k]=(y[k+1]-y[k])*c[k]*c[k];
This can be solved using the Gau?-Seidel iteration (was this not invented exactly for this system?) or your favorite Krylov space solver.
这可以使用 Gau?-Seidel 迭代（这不是专门为这个系统发明的吗？）或您最喜欢的 Krylov 空间求解器来解决。
回答by HymanOLantern
Cubic b-spline has been recently descried in a series of papers by Unser, Thévenaz et al., see among others
最近，Unser、Thévenaz等人在一系列论文中描述了三次 b 样条。，见其他
M. Unser, A. Aldroubi, M. Eden, "Fast B-Spline Transforms for Continuous Image Representation and Interpolation", IEEE Trans. Pattern Anal. Machine Intell., vol. 13, n. 3, pp. 277-285, Mar. 1991.
M. Unser, "Splines, a Perfect Fit for Signal and Image Processing", IEEE Signal Proc. Mag., pp. 22- 38, Nov. 1999.
M. Unser、A. Aldroubi、M. Eden，“用于连续图像表示和插值的快速 B 样条变换”，IEEE Trans。模式肛门。机器智能。，卷。13, n. 3，第 277-285 页，1991 年 3 月。
M. Unser，“样条，非常适合信号和图像处理”，IEEE Signal Proc。玛格。，第 22-38 页，1999 年 11 月。
and
和
P. Thévenaz, T. Blu, M. Unser, "Interpolation Revisited," IEEE Trans. on Medical Imaging, vol. 19, no. 7, pp. 739-758, July 2000.
P. Thévenaz、T. Blu、M. Unser，“重新审视插值法”，IEEE Trans。关于医学影像，第一卷。19，没有。7，第 739-758 页，2000 年 7 月。
Here are some guidelines.
这里有一些指导方针。
What are splines?
什么是样条？
Splines are piecewise polynomials that are smoothly connected together. For a spline of degree n, each segment is a polynomial of degree n. The pieces are connected so that the spline is continuous up to its derivative of degree n-1at the knots, namely, the joining points of the polynomial pieces.
样条是平滑连接在一起的分段多项式。对于度的样条n，每段都是度的多项式n。这些部分被连接起来，使得样条n-1在节点处连续直到它的阶导数，即多项式部分的连接点。
How can splines be constructed?
如何构造样条？
The zero-th order spline is the following
零阶样条如下
All the other splines can be constructed as
所有其他样条都可以构造为
where the convolution is taken n-1times.
其中卷积是采取n-1时间。
Cubic splines
三次样条
The most popular splines are cubic splines, whose expression is
最流行的样条是三次样条，其表达式为
Spline interpolation problem
样条插值问题
Given a function f(x)sampled at the discrete integer points k, the spline interpolation problem is to determine an approximation s(x)to f(x)expressed in the following way
给定的函数f(x)在所述离散整数点采样k，样条内插的问题是确定的近似s(x)，以f(x)下列方式表示
where the ck's are interpolation coefficients and s(k) = f(k).
其中ck是插值系数和s(k) = f(k)。
Spline prefiltering
样条预过滤
Unfortunately, starting from n=3on,
不幸的是，从此n=3开始，
so that the ck's are not the interpolation coefficients. They can be determined by solving the linear system of equations obtained by forcing s(k) = f(k), namely,
所以ck's 不是插值系数。它们可以通过求解强制获得的线性方程组来确定s(k) = f(k)，即，
Such an equation can be recast in a convolution form and solved in the transformed z-space as
这样的方程可以以卷积形式重铸并在变换z空间中求解为
where
在哪里
Accordingly,
因此，
Proceeding this way is always preferable than affording the solution of a linear system of equations by, for example, LUdecomposition.
以这种方式进行总是比提供线性方程组的解（例如，LU分解）更可取。
The solution to the above equation can be determined by noticing that
上述方程的解可以通过注意到
where
在哪里
The first fraction is representative of a causal filter, while the second one is representative of an anticausal filter. Both of them are illustrated in the figures below.
第一个分数代表因果过滤器，而第二个分数代表反因果过滤器。两者都如下图所示。
Causal filter
因果过滤器
Anticausal filter
反因果过滤器
In the last figure,
在最后一张图中，
The output of the filters can be expressed by the following recursive equations
滤波器的输出可以用以下递归方程表示
The above equations can be solved by first determining "initial conditions" for c-and c+. On assuming a periodic, mirroredinput sequence fksuch that
上述等式可以通过首先确定“初始条件”来解决c-和c+。假设一个周期性的、镜像的输入序列fk使得
then it can be shown that the initial condition for c+can be expressed as
那么可以证明 的初始条件c+可以表示为
while the initial condition for c+can be expressed as
而 的初始条件c+可以表示为
回答by warazdr
// File: cbspline-test.cpp
// 文件：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'
*/
// ========================================================
