數(shù)值分析（擬合、插值和逼近）之?dāng)?shù)據(jù)插值方法（線性插值、二次插值、Cubic插值、埃米爾特、拉格朗日多項(xiàng)式插值、牛頓插值、樣條插值）（含opengl程序）

戰(zhàn)神之家 2019-02-13

展開全文

插值、擬合和逼近的區(qū)別

據(jù)維基百科，科學(xué)和工程問題可以通過諸如采樣、實(shí)驗(yàn)等方法獲得若干離散的數(shù)據(jù)，根據(jù)這些數(shù)據(jù)，我們往往希望得到一個(gè)連續(xù)的函數(shù)（也就是曲線）或者更加密集的離散方程與已知數(shù)據(jù)相吻合，這過程就叫做擬合。通過擬合得到的函數(shù)獲得未知點(diǎn)的數(shù)據(jù)的方法，叫做插值。其中，擬合函數(shù)經(jīng)過所有已知點(diǎn)的插值方法，叫做內(nèi)插。
擬合是已知點(diǎn)列，從整體上靠近它們；插值是已知點(diǎn)列并且完全經(jīng)過點(diǎn)列；逼近是已知曲線，或者點(diǎn)列，通過逼近使得構(gòu)造的函數(shù)無限靠近它們。

最小二乘意義下的擬合，是要求擬合函數(shù)與原始數(shù)據(jù)的均方誤差達(dá)到極小，是一種整體意義的逼近，對局部性質(zhì)沒有要求。而所謂“插值”，就是要在原有離散數(shù)據(jù)之間“插入”一些值，這就要求插值函數(shù)必須通過所有的離散點(diǎn)，插值函數(shù)在離散點(diǎn)之外的那些點(diǎn)都相當(dāng)于“插入”的值。插值有全局插值，也有局部插值（比如分段線性插值），插值誤差通?？紤]的是逐點(diǎn)誤差或最大模誤差，插值的好壞往往通過某些局部的性質(zhì)來體現(xiàn)，比如龍格現(xiàn)象或吉布斯振蕩。

[知乎擬合與插值的區(qū)別？]

插值方法

多項(xiàng)式插值

對于大部分多項(xiàng)式插值函數(shù)，插值點(diǎn)的高度值可以視為所有（或某些）節(jié)點(diǎn)高度值的線性組合，而線性組合的系數(shù)一般是x坐標(biāo)的多項(xiàng)式函數(shù)，稱作基函數(shù)。對于一個(gè)節(jié)點(diǎn)的基函數(shù)，它在x等于該節(jié)點(diǎn)的x時(shí)等于1，在x等于其他節(jié)點(diǎn)的x時(shí)等于0。這就保證曲線必定經(jīng)過所有節(jié)點(diǎn)，所以屬于內(nèi)插方法。

在本小節(jié)，均以一組隨機(jī)數(shù)作為已知的高度值，使它們對應(yīng)于間隔固定的x坐標(biāo)，使用不同的插值函數(shù)獲得各已知點(diǎn)（稱為插值函數(shù)的節(jié)點(diǎn)）之外其它x坐標(biāo)所對應(yīng)的高度值，畫出這些點(diǎn)所對應(yīng)的曲線。再把所有高度值轉(zhuǎn)換成灰度值，以顏色的變化比較各插值函數(shù)。

原點(diǎn)列如圖：（假定橫向?yàn)閤，縱向?yàn)閥。各點(diǎn)x坐標(biāo)的間隔是固定的，但y坐標(biāo)是隨機(jī)的）

線性插值

線性插值是用一系列首尾相連的線段依次連接相鄰各點(diǎn)，每條線段內(nèi)的點(diǎn)的高度作為插值獲得的高度值。

以(xi,yi)表示某條線段的前一個(gè)端點(diǎn)，(x(i+1),y(i+1))表示該線段的后一個(gè)端點(diǎn)，則對于在[xi,x(i+1)]范圍內(nèi)的橫坐標(biāo)為x的點(diǎn)，其高度y為：

$p(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0). \,\!$

為便于與后面各函數(shù)比較，寫成比較對稱的形式：

其中，yi和y(i+1)的兩個(gè)參數(shù)稱為基函數(shù)，二者之和為1，分別代表yi和y(i+1)對插值點(diǎn)高度的權(quán)值。

[wikipedia線性插值]

插值圖像如下：

將高度轉(zhuǎn)化為灰度，得到如下條帶：

線性插值的特點(diǎn)是計(jì)算簡便，但光滑性很差。如果用線性插值擬合一條光滑曲線，對每一段線段，原曲線在該段內(nèi)二階導(dǎo)數(shù)絕對值的最大值越大，擬合的誤差越大。

二次插值

如果按照線性插值的形式，以每3個(gè)相鄰點(diǎn)做插值，就得到了二次插值：

OpenGL實(shí)現(xiàn)代碼如下：

void quadratic(float p[20][2])
{
float x,y;
int i;
float x01,x02,x12;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(i=0;i<20;i+=2)
{
x01=p[i][0]-p[i+1][0];
x02=p[i][0]-p[i+2][0];
x12=p[i+1][0]-p[i+2][0];
for(x=p[i][0];x<=p[i+2][0];x+=1.0)
{
y=(x-p[i+1][0])*(x-p[i+2][0])/x01/x02*p[i][1]-(x-p[i][0])*(x-p[i+2][0])/x01/x12*p[i+1][1]+(x-p[i][0])*(x-p[i+1][0])/x02/x12*p[i+2][1];
glVertex2f(x,y);
}
}
glEnd();
}

二次（分段）插值圖像如下：

轉(zhuǎn)換成灰度值如圖：

二次插值在每段二次曲線內(nèi)是光滑的，但在每條曲線的連接處其光滑性可能甚至比線性插值還差。二次插值只適合3個(gè)節(jié)點(diǎn)的情形，當(dāng)節(jié)點(diǎn)數(shù)超過3個(gè)時(shí)，就需要分段插值了。

Cubic插值

如果想要保證各段曲線連接處光滑（一階導(dǎo)數(shù)相同），并且不想使用除法運(yùn)算，可以考慮Cubic插值函數(shù)：

其中，v代表插值點(diǎn)，v0、v1、v2、v3代表4個(gè)連續(xù)的節(jié)點(diǎn)。t取值為[0,1]，將會產(chǎn)生一段連接v1和v2的曲線。也就是說，如果有n個(gè)節(jié)點(diǎn)，Cubic插值函數(shù)將會產(chǎn)生(n-2)段曲線，位于首尾兩端的節(jié)點(diǎn)不會納入曲線。

實(shí)現(xiàn)代碼如下：

float cubic(float v0,float v1,float v2,float v3,float x)
{
float v32=v3-v2;
float v01=v0-v1;
float v20=v2-v0;
float temp=(v32-v01)*x;
temp=(temp+v01+v01-v32)*x;
temp=(temp+v20)*x+v1;
return temp;
}
void drawCubic(float p[20])
{
float x,y;
int step;
float delta;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(step=0;step<17;step++)
{
for(delta=0.0;delta<=1.0;delta+=0.05)
{
x=delta*(p[step+1][0]-p[step][0])+p[step][0];
y=cubic(p[step],p[step+1],p[step+2],p[step+3],delta);
glVertex2f(x,y);
}
}
glEnd();
}

Cuibc插值圖像如下：

轉(zhuǎn)化成灰度如圖：

Cubic插值可以產(chǎn)生整體上光滑的曲線，但容易產(chǎn)生較劇烈的波動，使得曲線的最高點(diǎn)比最高的節(jié)點(diǎn)還高、曲線的最低點(diǎn)比最低的節(jié)點(diǎn)還低。所以對顏色等取值有嚴(yán)格的上下界的數(shù)據(jù)進(jìn)行插值時(shí)，會造成曲線的截取，破壞其光滑性。比如顏色（RGB三個(gè)分量取值范圍都是[0,255]），如果最高的節(jié)點(diǎn)是255，最低的節(jié)點(diǎn)是0，那么插值后負(fù)數(shù)會被截取成0，大于255的數(shù)會被截取成255。

拉格朗日多項(xiàng)式插值

解決插值函數(shù)的直接構(gòu)造問題以及插值誤差的估計(jì)問題，但是同時(shí)也使得插值函數(shù)值的計(jì)算變得復(fù)雜。

依照線性插值和二次插值的思路，可以增加基函數(shù)分子和分母的階數(shù)，構(gòu)造拉格朗日插值多項(xiàng)式：

一個(gè)n次的拉格朗日插值函數(shù)可以繪制經(jīng)過(n+1)個(gè)節(jié)點(diǎn)的曲線，但運(yùn)算量非常大。而且在次數(shù)比較高時(shí)，容易產(chǎn)生劇烈的震蕩（龍格現(xiàn)象）。所以要選擇位置特殊的節(jié)點(diǎn)（比如切比雪夫多項(xiàng)式的零點(diǎn)）進(jìn)行插值，或使用多個(gè)次數(shù)較低的拉格朗日函數(shù)分段插值。（關(guān)于拉格朗日多項(xiàng)式和龍格現(xiàn)象，詳見維基百科鏈接）

分段插值實(shí)現(xiàn)代碼如下：

//n為次數(shù)，nmax為節(jié)點(diǎn)的總數(shù)。n不大于nmax
void lagrange(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=0.0;
for(j=0;j<n;j++)
{
temp=1.0;
for(t=0;t<n;t++)
{
if(t==j)
continue;
temp*=(x-p[i+t][0])/(p[i+j][0]-p[i+t][0]);
}
y+=temp*p[i+j][1];
}
glVertex2f(x,y);
}
glEnd();
}
}

使用4次拉格朗日多項(xiàng)式分段插值：

轉(zhuǎn)化為灰度：

拉格朗日多項(xiàng)式插值也存在連接處不光滑的問題。

如果直接使用20次的拉格朗日插值，得到的圖像如下：

這樣的插值曲線顯然是不能容忍的~

牛頓插值

從算法角度考慮，提出便于計(jì)算的插值方法，也就是牛頓插值公式。拉格朗日插值每增加一個(gè)節(jié)點(diǎn)，整個(gè)函數(shù)要重新計(jì)算，計(jì)算量巨大。而牛頓插值每增加一個(gè)點(diǎn)只需要在多項(xiàng)式的最后增加一項(xiàng)，而且各基函數(shù)的系數(shù)可以遞歸計(jì)算，減少了很多計(jì)算量。

均差（Divided differences）

是遞歸除法過程。在數(shù)值分析中，也稱差商（Difference quotient），可用于計(jì)算牛頓多項(xiàng)式形式的多項(xiàng)式插值的系數(shù)。

定義

[數(shù)值分析鐘爾杰電子科大]

差商表（高階差商是兩個(gè)低一階差商的差商）
	$0$ 階差商	階差商	階差商	階差商	$\ldots$	階差商
$x_{0}$	$f[x_{{0}}]$
$x_{{1}}$	$f[x_{{1}}]$	$f[x_{{0}},x_{{1}}]$
$x_{{2}}$	$f[x_{{2}}]$	$f[x_{{1}},x_{{2}}]$	$f[x_{{0}},x_{{1}},x_{{2}}]$
$x_{{3}}$	$f[x_{{3}}]$	$f[x_{{2}},x_{{3}}]$	$f[x_{{1}},x_{{2}},x_{{3}}]$	$f[x_{{0}},x_{{1}},x_{{2}},x_{{3}}]$
$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
$x_{{k}}$	$f[x_{{k}}]$	$f[x_{{k-1}},x_{{k}}]$	$f[x_{{k-2}},x_{{k-1}},x_{{k}}]$	$f[x_{{k-3}},x_{{k-2}},x_{{k-1}},x_{{k}}]$	$\ldots$	$f[x_{{0}},\ldots ,x_{{k}}]$

由函數(shù)表出發(fā)，從上往下從左往右依次計(jì)算出一階，二階，。。。，各階均差。

[wikipedia均差]

牛頓多項(xiàng)式

將n次插值多項(xiàng)式寫作如下形式：

$N(x)=[y_{0}]+[y_{0},y_{1}](x-x_{0})+\cdots +[y_{0},\ldots ,y_{k}](x-x_{0})(x-x_{1})\cdots (x-x_{{k-1}})$

也就是牛頓插值公式中的各項(xiàng)系數(shù)就是均差表中計(jì)算出的各階均差的第一個(gè)數(shù)據(jù)！

即牛頓插值公式為：

[wikipedia牛頓多項(xiàng)式]

[wikipedia均差小節(jié)：牛頓插值法]

實(shí)現(xiàn)代碼如下：

void newton(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
float f[20][20];
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
for(t=0;t<n;t++)
{
f[t][0]=(p[i+t][1]-p[i+t+1][1])/(p[i+t][0]-p[i+t+1][0]);
for(j=1;j<=t;j++)
f[t][j]=(f[t-1][j-1]-f[t][j-1])/(p[i+t-j][0]-p[i+t+1][0]);
}
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=p[i][1];
temp=1.0;
for(j=0;j<n;j++)
{
temp*=x-p[i+j][0];
y+=temp*f[j][j];
}
glVertex2f(x,y);
}
glEnd();
}
}

可能由于計(jì)算精度的原因，牛頓插值繪制的曲線與拉格朗日插值的曲線略有不同。但次數(shù)較高時(shí)，牛頓插值也會產(chǎn)生劇烈的震蕩。分段4次牛頓插值圖像如下：

轉(zhuǎn)化成灰度如下：

牛頓插值也存在連接處不光滑的缺陷。

埃爾米特插值

以上各多項(xiàng)式插值方法的插值條件都是各節(jié)點(diǎn)的坐標(biāo)，在以低階函數(shù)分段插值時(shí)雖然可以保持曲線的穩(wěn)定（比較平緩），但在各分段曲線的連接處無法保持光滑（一階導(dǎo)數(shù)相等）。埃爾米特插值方法不但規(guī)定了各節(jié)點(diǎn)的坐標(biāo)值，還規(guī)定了曲線在每個(gè)節(jié)點(diǎn)的各階導(dǎo)數(shù)，這樣就可以既保持曲線的穩(wěn)定，又保證在連接處足夠光滑。

以3次二重（"m重"就是規(guī)定坐標(biāo)和曲線在所有節(jié)點(diǎn)處1到m-1階導(dǎo)數(shù)的值）埃爾米特插值為例：

4個(gè)基函數(shù)滿足分別只在y0，y1，y0的導(dǎo)數(shù)，y1的導(dǎo)數(shù)處等于1，而在其他3個(gè)條件下等于0?？梢园寻柮滋夭逯悼醋鲗ψ鴺?biāo)和導(dǎo)數(shù)的插值的組合。

曲線在每個(gè)節(jié)點(diǎn)的導(dǎo)數(shù)可以根據(jù)相鄰節(jié)點(diǎn)和它的相對位置確定，也可以完全隨機(jī)生成。只要位置比較高和比較低的節(jié)點(diǎn)的導(dǎo)數(shù)絕對值不是很大，就可以使整條曲線落在最高與最低節(jié)點(diǎn)定義的帶狀區(qū)域內(nèi)，這樣就可以避免對插值的截取。

對于本節(jié)的示例節(jié)點(diǎn)，一種可能的導(dǎo)數(shù)值如下：grad[20]={-4.0,4.0,0.5,-2.0,1.0,-2.0,-2.0,2.0,1.0,1.0,0.0,-1.0,-4.0,3.0,0.0,-3.0,3.0,0.0,-4.0,-4.0};

實(shí)現(xiàn)代碼如下：

//p[i][0],p[i][1]為點(diǎn)i的坐標(biāo)，p[i][2]為曲線在點(diǎn)i的導(dǎo)數(shù)
//nmax為節(jié)點(diǎn)的總數(shù)
void hermite(float p[][3],int nmax)
{
float x,y;
float a1,a0,b1,b0;
float x00,x11;
int i;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<nmax;i++)
{
glBegin(GL_LINE_STRIP);
x00=p[i][0];
x11=p[i+1][0];
for(x=p[i][0];x<=p[i+1][0];x+=1.0)
{
temp=(x11-x00)*(x11-x00);
a0=(x11-3*x00+2*x)*(x11-x)*(x11-x)/temp/(x11-x00);
a1=(-x00+3*x11-2*x)*(x-x00)*(x-x00)/temp/(x11-x00);
b0=(x-x00)*(x-x11)*(x-x11)/temp;
b1=(x-x00)*(x-x00)*(x-x11)/temp;
y=a0*p[i][1]+a1*p[i+1][1]+b0*p[i][2]+b1*p[i+1][2];
glVertex2f(x,y);
}
glEnd();
}
}

分段3次埃爾米特插值圖像如下：

轉(zhuǎn)化為灰度如下：

可見雖然n節(jié)點(diǎn)的埃爾米特插值是由(n-1)段曲線構(gòu)成，但在每一個(gè)連接處都是光滑的。

樣條插值

B-樣條

B-樣條是Bezier樣條的一般化，也可推廣為NURBS樣條。先來介紹一下B-樣條：

該式定義了一個(gè)n次的B-樣條，它有(m+1)個(gè)控制點(diǎn)（樣條曲線的“節(jié)點(diǎn)”稱作控制點(diǎn)，因?yàn)檫@些點(diǎn)控制曲線的彎曲方向和程度，但曲線不一定經(jīng)過這些點(diǎn)），也就有(m+1)個(gè)基函數(shù)。一般繪制的是完整的曲線，u(min)取0，u(max)取1。當(dāng)u取值均勻時(shí)，該樣條稱作均勻B-樣條。當(dāng)m=n時(shí)，B-樣條退化為Bezier樣條。

函數(shù)繪制的曲線始終落在其控制點(diǎn)的凸包（包含一個(gè)點(diǎn)集所有點(diǎn)的凸多邊形，而且該凸多邊形所有頂點(diǎn)都來自這個(gè)點(diǎn)集）中。對于一個(gè)n次的B-樣條，改變一個(gè)控制點(diǎn)的位置，只改變它所在的n段曲線（由n+1個(gè)控制點(diǎn)定義，且以該點(diǎn)起始）的形狀，而不對其余的(m-n)段曲線產(chǎn)生影響。

Bezier樣條

Bezier（貝塞爾）樣條是法國工程師皮埃爾·貝塞爾（Pierre Bézier）在1962年為了設(shè)計(jì)汽車主體外形曲線而提出的。其一般式表達(dá)式為：

其中，u取值為[0,1]，pk(k=0,...,n)為(n+1)個(gè)節(jié)點(diǎn)，n稱為階數(shù)。

Bezier樣條還可以遞歸定義為：

含義是n階Bezier樣條是兩條(n-1)階Bezier樣條的插值。

當(dāng)階數(shù)降為1時(shí)，Bezier插值退化成線性插值。改變?nèi)我庖粋€(gè)控制點(diǎn)的位置，整條曲線的形狀都會發(fā)生變化。

比較常用的Bezier樣條是3次Bezier：

Beizer樣條在首尾端的切線是前兩個(gè)點(diǎn)和最后兩個(gè)點(diǎn)的連線。除了第一個(gè)點(diǎn)和最后一個(gè)點(diǎn)，其他控制點(diǎn)一般都不在曲線上。

3階（4控制點(diǎn)）的Bezier函數(shù)圖像如下：（黑色曲線為Bezier曲線，藍(lán)色折線為控制點(diǎn)的連線）

Bezier樣條可以實(shí)現(xiàn)非常快速的運(yùn)算。為了方便說明，將3次Bezier樣條表示成如下形式：

對于任意給定的u，可以通過合并的方式將原來的7次乘法、4次加法減少為3次乘法、3次加法：

一般情況下，應(yīng)用Bezier樣條繪制曲線時(shí)，都是先給定一個(gè)很小的步長t（步長足夠小才能保證Bezier曲線的精確），讓t從0取到1，從頭到尾繪制整條曲線。在t不變的條件下，可以使用更快速的差分方法，利用前一個(gè)點(diǎn)計(jì)算出下一個(gè)點(diǎn)的值，將每步的計(jì)算量減小到只有3次加法：

只需要在繪制曲線之前計(jì)算4個(gè)常數(shù)，就可以很快地計(jì)算出曲線上的所有點(diǎn)：

采用這種方式，Bezier樣條的運(yùn)算量只隨階數(shù)線性增長。

實(shí)現(xiàn)代碼如下：

void bezier3(float xx0,float yy0,float xx1,float yy1,float xx2,float yy2,float xx3,float yy3)
{
GLint i;
GLfloat x,y;
GLdouble ax,bx,cx,ay,by,cy;
GLdouble t;
GLdouble dx,d2x,dy,d2y,delta;
GLdouble d3x,d3y;
ax=-xx0+3.0*xx1-3.0*xx2+xx3;
bx=3.0*xx0-6.0*xx1+3.0*xx2;
cx=-3.0*xx0+3.0*xx1;
ay=-yy0+3.0*yy1-3.0*yy2+yy3;
by=3.0*yy0-6.0*yy1+3.0*yy2;
cy=-3.0*yy0+3.0*yy1;
delta=0.0005;
d3x=6.0*ax*delta*delta*delta;
d3y=6.0*ay*delta*delta*delta;
x=xx0;
y=yy0;
glBegin(GL_LINE_STRIP);
glVertex2f(x,y+200);
d2x=6.0*ax*delta*delta*delta+2.0*bx*delta*delta;
dx=ax*delta*delta*delta+bx*delta*delta+cx*delta;
d2y=6.0*ay*delta*delta*delta+2.0*by*delta*delta;
dy=ay*delta*delta*delta+by*delta*delta+cy*delta;
for(t=0.0;t<=1.0;t+=delta)
{
d2x+=d3x;
dx+=d2x;
x+=(float)dx;
d2y+=d3y;
dy+=d2y;
y+=(float)dy;
glVertex2f(x,y);
}
glEnd();
}

NURBS樣條

NURBS（Non-Uniform Rational B-Splines 非均勻有理B-樣條），是貝塞爾樣條的推廣?！胺蔷鶆颉钡囊馑际强刂泣c(diǎn)的間隔可以是不均勻的，“有理”的意思是各控制點(diǎn)帶有不同的權(quán)值。和Bezier樣條相比，它對曲線形狀的控制更自由：

其基函數(shù)B(k,d)與B-樣條的基函數(shù)相同，w(k)為各點(diǎn)的權(quán)因子。和B-樣條一樣，改變一個(gè)控制點(diǎn)的位置，只改變它所在的n段曲線的形狀，而不對其余的(m-n)段曲線產(chǎn)生影響。

插值的一些應(yīng)用

噪聲

噪聲圖像一般需要二維或三維插值函數(shù)，不過了解了各一維插值函數(shù)，就很容易擴(kuò)展到二維和三維了。關(guān)于二維噪聲圖像的產(chǎn)生，請參見《二維噪聲圖像》。

水波

下面這個(gè)圖像就是用插值函數(shù)繪制的模擬水波的三維網(wǎng)格（靜態(tài)圖）：

雖然這張網(wǎng)格看起來似乎需要二維插值，但這張網(wǎng)格是由6個(gè)不同波長和頻率的Gestner合成的，每個(gè)獨(dú)立的波形只需要一維插值生成。

本站是提供個(gè)人知識管理的網(wǎng)絡(luò)存儲空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn)。請注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購買等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請點(diǎn)擊一鍵舉報(bào)。

轉(zhuǎn)藏 分享

QQ空間 QQ好友新浪微博微信

獻(xiàn)花（0） +1

來自：戰(zhàn)神之家 > 《待分類》

舉報(bào)/認(rèn)領(lǐng)