数值分析15样条插值

样条插值多元插值《数值分析》14

http://www.tspline/先是雷诺和雪铁龙工作的

PauldeCasteljau和PierreBézier,随后美国通用汽车的其它人一起推动了现在称为三次样条和Bézier

样条的建立。样条是通过很少的控制点就能够生成复杂平滑曲线的方法。Ref:/design/movies-before-after-green-screen-cgi/只给定离散数据,不知道函数,当然也不知道其导数(或导数值)。定义

1

给定区间[a,b]上的一个分划:a=x0<x1<…<xn=b分段三次多项式满足如下条件:

(1)Sk(xk-1)=yk-1,Sk(xk)=yk

(k=1,···,n)

(2)

S'k(xk)=S'k+1(xk)

(k=1,···,n-1)(3)S''k(xk)=S''k+1(xk)

(k=1,···,n-1)则称

S(x)为三次样条插值函数。当xk-1

≤x≤xk时,

Sk(x)=ak+bkx+ckx2+dkx3(k=1,…,n)条件:(1)Sk(xk-1)=yk-1,Sk(xk)=yk

(k=1,···,n)(2)S'k(xk)=S'k+1(xk)

(k=1,···,n-1)(3)S''k(xk)=S''k+1(xk)

(k=1,···,n-1)由样条插值函数满足条件,可以建立(4n-2)个方程！

n个三次多项式,共4n个待定系数。自然边界条件:S''(x0)=0,S''(xn)=0例1已知f(–1)=1,f(0)=0,f(1)=1。构造分段三次多项式是满足自然边界的样条函数。解:如何构造满足条件的样条插值函数？分段Hermite插值插值函数已经满足插值条件和一阶导数连续。怎么保证二阶导数连续

？

当xk-1

≤x≤xk时,

Sk(x)=ak+bkx+ckx2+dkx3(k=1,…,n)当

xk-1

≤x≤xk时,

Sk(x)由如下基函数组成当

xk≤x≤xk+1时,

Sk+1

(x)由如下基函数组成由(k=1,2,···,n-1)

待定未知数n+1个,建立n-1个方程。自然样条的导数值满足:设自然边界条件成立即(k=1,2,······,n-1)自然样条的导数值满足:(k=1,2,······,n-1)回顾:严格主对角占优矩阵一定是非奇异的。Demo

x=-5:5;y=1./(x.^2+1);u=-5:.01:5;v1=polyinterp(x,y,u);plot(x,y,'o',u,v1,'-')holdon,v2=piecelin(x,y,u);plot(u,v2,'r-')holdon,v3=spline(x,y,u);plot(u,v3,'b-')interpgui

样条插值是插值函数的光滑性与局部单调性之间的折衷方案。rowcol二维插值a bc daa bbaa bbcc ddcc ddzero-orderfirst-ordera (a+b)/2 b(a+c)/2

(a+b+c+d)/4(b+d)/2c (c+d)/2 d双线性插值

一维数据的插值interp1x=0:2:24;y=[22,21,19,18,20,24,27,32,31,28,26,23,22];xi=0:0.1:24;yi=interp1(x,y,xi,’spline’);plot(xi,yi,’-’,x,y,’o’);二维网格数据的插值interp2[x,y]=meshgrid(1:5,1:3);z=[8281808284;7963616581;8484828586];mesh(x,y,z)[xi,yi]=meshgrid(1:0.2:5,1:0.2:3);zi=interp2(x,y,z,xi,yi,'spline');mesh(xi,yi,zi);x=double(imread('fl_orig.pgm'));figure,subplot(221),imshow(x,[]);%helpinterp2(orimresize)%zero-orderinterpolation(replication)y0=interp2(x,'nearest');subplot(222),imshow(y0,[]);%first-orderinterpolation(bilinearinterpolation)y1=interp2(x,'linear');subplot(223),imshow(y1,[]);%bicubicinterpolationy2=interp2(x,'cubic');subplot(224),imshow(y2,[]);%%colorimagex=double(imread('fl_orig.ppm'));figure,subplot(221),imshow(x/255,[]);fori=1:3%zero-orderinterpolation(replication)y0(:,:,i)=interp2(x(:,:,i),'nearest');%first-orderinterpolation(bilinearinterpolation)y1(:,:,i)=interp2(x(:,:,i),'linear');%bicubicinterpolationy2(:,:,i)=interp2(x(:,:,i),'cubic');endsubplot(222),imshow(y0/255,[]);subplot(223),imshow(y1/255,[]);subplot(224),imshow(y2/255,[]);%Howtoimplementimageinterpolationbyyourself?figure,subplot(221),imshow(x,[]);%zero-orderinterpolation(replication)[M,N]=size(x);z0=zeros(2*M,2*N);z0(1:2:2*M,1:2:2*N)=x;h=[11;11];z0=filter2(h,z0);subplot(222),imshow(z0,[]);%first-orderinterpolation(bilinearinterpolation)z1=zeros(2*M,2*N);z1(1:2:2*M,1:2:2*N)=x;h=[121;242;121]/4;z1=filter2(h,z1);subplot(223),imshow(z1,[]);%cubicinterpolationz2=zeros(2*M,2*N);z2(1:2:2*M,1:2:2*N)=x;h1=[-1091690-1]/16;h=h1'*h1;z2=filter2(h,z2);subplot(224),imshow(z2,[]);Ref：Cubicconvolutioninterpolationfordigitalimageprocessing最近邻插值low-resolutionimage