Lagrange、分段线性和三次样条三种插值方法计算
数值分析上机实验报告(插值)
数值分析第一次上机练习实验报告——Lagrange 插值与三次样条插值一、 问题的描述设()2119f x x =+, []1,1x ∈-,取15iix =-+,0,1,2,...,10i =.试求出10次Lagrange 插值多项式()10L x 和三次样条插值函数()S x (采用自然边界条件),并用图画出()f x ,()10L x ,()S x .二、 方法描述——Lagrange 插值与三次样条插值我们取15i ix =-+,0,1,2,...,10i =，通过在i x 点的函数值()2119i i f x x =+来对原函数进行插值，我们记插值函数为()g x ，要求它满足如下条件：()()21,0,1,2,...,1019i i i g x f x i x ===+ (1)我们在此处要分别通过Lagrange 插值(即多项式插值)与三次样条插值的方法对原函数()2119f x x=+进行插值，看两种方法的插值结果，并进行结果的比较。
10次的Lagrange 插值多项式为：()()10100i i i L x y l x ==∑ (2)其中：()21,0,1,2,...,1019i i iy f x i x ===+ 以及()()()()()()()()()011011......,0,1,2,...,10......i i n i i i i i i i n x x x x x x x x l x i x x x x x x x x -+-+----==----我们根据(2)进行程序的编写，我们可以通过几个循环很容易实现函数的Lagrange 插值。
理论上我们根据区间[]1,1-上给出的节点做出的插值多项式()n L x 近似于()f x ，而多项式()n L x 的次数n 越高逼近()f x 的精度就越好。
但实际上并非如此，而是对任意的插值节点，当n →+∞的时候()n L x 不一定收敛到()f x ；而是有时会在插值区间的两端点附近会出现严重的()n L x 偏离()f x 的现象，即所谓的Runge 现象。
拉格朗日(Lagrange)插值
l0 (x)
再由另一条件 l0 (x0
c(x
) 1
x1 )( x x2
确定系数
)
c
(x0
1 x1)( x0
x2
)
从而导出
l0 (x)
(x (x0
x1)( x x2 ) x1 )( x0 x2 )
类似地可以构造出满足条件： l1(x1) 1, l1(x0 ) 0,
的插值多项式
l1 ( x)
li (x的) 插值
lk (x0 ) 0,,lk (xk1) 0,lk (xk ) 1,lk (xk1 ) 0,,lk (xn ) 0
即
lk
(xi )
ki
1 0
(i k) (i k)
由条件 lk (xi ) 0 ( i k)知, x0 , x1,, xk1, xk1,, xn
都是n次 lk (x) 的零点,故可设
即
a0
,
a1
,
a
满足下面
2
的代数方程组：
a0 a0
a1 x0 a1 x1
a2 x02 a2 x12
y0 y2
y0
y1
y1
O
x0
x1
x2
该三元一 次方程组 的系数矩阵
1 1
x0 x1
1 x2
y=L2(x)
y=f(x) x
x02 x12
x
2 2
的行列式是范德蒙行列式，当 x0 x1 x2 时，
(x ( x1
x0 )( x x2 ) x0 )( x1 x2 )
l1(x2 ) 0
及满足条件：l2 (x2 ) 1, l2 (x0 ) 0, l2 (x1) 0 的插值多项式
Matlab插值方法学习Hermite,lagrange,Newton
大家都知道插值在数学建模中很重要，现在介绍几种常用插值下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。
1.拉格朗日多项式插值拉格朗日插值就是给定n个数，让你用不超过n-1次的多项式你逼近它，当然这n个点要能满足多项式。
这是一种最基本的思想，计算很简单，先计算n个基函数，基函数可以自己上网搜一下，因为这里打出公式有点麻烦。
然后就是把每个点的y值乘以他的基函数，把这n个式子相加，最后化简就ok了。
下面我把代码写出来，我这些代码全是自己写的，注释比较详细，这里只以lagrange为例，其余都放在附件里了。
%定义myLagrange函数，参数为向量x，y，由用户调用该函数时输入function L=myLagrange (x,y)%n 插值结点的个数n=length(x);%L myLagrange函数计算的多项式系数行列式L=zeros(1,n);%%使用双重for循环，第一个for循环是fori=1:n%aa=1;%ww=1;%for循环for j=1:n%如果i不等于jif j~=i%累加法计算aa=a*(x(i)-x(j));%用向量乘法函数conv计算ww=conv(w,[1,-x(j)]);%if语句结束符end%第二个for循环结束符end%递归法计算L，其中y(i)/a*w表示第i个元素L=y(i)/a*w+L;%第一个for结束符end没错，就这么几句代码，所以很简单的。
2.牛顿插值牛顿插值其实是为了解决拉格朗日插值不能增加新的点来说的。
拉格朗日插值只能接受给定的那么多点，了然后插值。
如果你想再加一个点，它会重新开始计算，这个很费时间和内存。
因此牛顿插值就诞生了。
了解牛顿插值前要学习下差商和差分两个简单的概念。
Newton 插值的优点是：每增加一个节点，插值多项式只增加一项，即因而便于递推运算。
而且Newton 插值的计算量小于Lagrange 插值。
常见的插值方法及其原理
常见的插值方法及其原理1. 拉格朗日插值法（Lagrange Interpolation）拉格朗日插值法是一种基于多项式的插值方法，通过n+1个已知点的函数值来构造一个n次多项式。
具体的计算公式如下：L(x) = Σ[yk * lk(x)], k=0 to n其中yk为已知点（xi, yi）的函数值，lk(x)为拉格朗日基函数，定义为：lk(x) = Π[(x - xj)/(xi - xj)], j=0 to n, j≠k拉格朗日插值法的原理是通过构造一个通过已知点的n次多项式，来代替未知函数的近似值。
利用拉格朗日基函数的性质，可以保证插值多项式通过已知点。
2. 牛顿插值法（Newton Interpolation）牛顿插值法是一种递推的插值方法，通过已知点的函数值和差商来逐步构造插值多项式。
差商的定义如下：f[x0]=y0f[x1]=(f[x1]-f[x0])/(x1-x0)f[x2]=(f[x2]-f[x1])/(x2-x1)...f[xn] = (f[xn] - f[xn-1]) / (xn - xn-1)利用差商的定义，可以得到牛顿插值多项式的表达式：N(x) = f[x0] + f[x0, x1](x-x0) + f[x0, x1, x2](x-x0)(x-x1) + ... + f[x0, x1, ..., xn](x-x0)(x-x1)...(x-xn)牛顿插值法的原理是通过递推计算差商来得到插值多项式。
通过使用差商来处理已知点的函数值差异，可以得到更高次的插值多项式。
3. 样条插值法（Spline Interpolation）样条插值法是一种基于分段低次插值函数的插值方法，常用的是三次样条插值。
样条插值法通过寻找一组分段函数，使得满足原函数的插值条件，并要求函数在每个插值点处的函数值、一阶导数和二阶导数连续。
这样可以保证插值函数在每个插值点处的平滑性。
三次样条插值法的原理是将整个插值区间划分为多个小区间，在每个小区间内使用三次多项式进行插值。
第二章 插值法_课堂演示实验范文
第二章 插值法_课堂演示实验问题: 分别以函数()sin f x x =和21()1f x x =+为例,在区间[5,5]x ∈-上,取105ix i n=-+ (0,...,)i n =为节点,用本章所学的Lagrange 、Newton 、Hermite 插值法以及分段线性、分段三次Hermite 插值法、Spline 插值法作插值计算.具体要求如下： 1．取2,4,8n =，通过在同一坐标系中作出被插函数与插值函数的图形的方法，来观察插值函数的图象与被插函数的位置关系，并给出观察到的的结论.2．通过计算 4.5x =时的插值误差,结合前面得到的直观结论,试对各种插值方法的应用作出述评.解答：一、各种插值方法的算法公式及MATLAB 通用程序设已知)(x f y =函数表为(,()),0,1,,,i i x f x i n =L 且当i j ≠时i j x x ≠. 1．Lagrange 插值法 算法公式 :∑==nk k k n x l x f x L 0)()()(, 其中),,1,0(,)(0n k x x x x x l jk ikj j nk =--∏=≠=余项为 ∏=-++=-=ni i n n x x n n f x L x f x R 0)()!1())(1()()()(ξ, (,)a b ξ∈与x 有关.通用程序1(此程序用得较多) :function yy=lagr1(x,y,xx)%用途：拉格朗日插值法求插值点xx （可以是多个）处的插值yy %格式：yy=lagr(x,y,xx), x 是节点向量，y 是节点对应的函数值向量， % xx 是插值点（可以是多个），yy 返回插值结果 m=length(x);n=length(y);if m~=n, error('向量x 与y 的长度必须一致');end s=0; for i=1:nt=ones(1,length(xx)); for j=1:n if j~=it=t.*(xx-x(j))/(x(i)-x(j)); end ends=s+t*y(i);endyy=s;通用程序2 :function [L ,C, l ,L1]= lagr2(X,Y)%输入的量:n+1个节点(xi,yi)的横坐标向量X，纵坐标向量Y；%输出的量:n次拉格朗日插值多项式L及其系数向量C，基函数l及其系数矩阵L1m=length(X); L=ones(m,m);for k=1: mV=1;for i=1:mif k~=iV=conv(V,poly(X(i)))/(X(k)-X(i));endendL1(k,:)=V; l(k,:)=poly2sym (V);endC=Y*L1;L=Y*l;l=vpa(l,4);L=vpa(L,4);通用程序3 :function [y,R]=lagr3(X,Y,x,M)%输入的量：X 是n+1个节点的横坐标向量，Y是纵坐标向量, x是以向量形式输入的m个插值点，M是被插函数在[a,b]区间上的n+1阶导数的最大值.%输出的量：y为m个插值构成的向量，R是误差限.n=length(X); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0; q1=1.0; c1=1.0;for j=1:nif j~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;ends=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;湖北民族学院理学院《数值分析》课程教案教案编写人：陈以平2．Newton 插值法 算法公式 :001001011()()[,]()[,,,]()()()n n n N x f x f x x x x f x x x x x x x x x -=+-++---L L L余项为()()()n n R x f x N x =-)(],,,,[010i ni n x x x x x x f -∏== 0(1)()()(1)!ni i f n x x n ξ=+=∏-+其中 (,)a b ξ∈与x 有关. 通用程序1(此程序用得较多) :function s=newton1(x,y,x0,nn)%Newton 插值,x 与y 为已知的插值点及其函数值 %x0为需要求的插值点向是,s 返回插对应插值%nn 为newton 插值多项式的次数,即nn 次newton 插值多项式 nx=length(x); ny=length(y); if nx~=nywarning('向量x 与y 的长度应该相同') return; endm=length(x0);%按照公式,对需要求的插值点x0的每个元素进行计算 for i=1:m t=0.0; %j=1; yy=y; kk=1;%求各级均差 while(kk<=nn) kk=kk+1; for k=kk:nxyy(k)=(yy(k)-yy(kk-1))/(x(k)-x(kk-1)); end end%求插值结果 t=yy(1); for k=2:nx u=1.0; jj=1;while(jj<k)u=u*(x0(i)-x(jj)); jj=jj+1; endt=t+yy(k)*u;ends(i)=t;endreturn通用程序2 :function [y,R]=newton2(X,Y,x,M)%输入的量：X 是n+1个节点的横坐标向量，Y是纵坐标向量, x是以向量形式输入的m个插值点，M是被插函数在[a,b]区间上的n+1阶导数的最大值.%输出的量：y为m个插值构成的向量，R是误差限.n=length(X); m=length(x);for t=1:mz=x(t); A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; q1=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C); C(d)=C(d)+A(k,k);endy(t)= polyval(C, z);endR=M*q1/c1;通用程序3 :function [A,C,L,wcgs,Cw]=newton3(X,Y)%输入的量:n+1个节点的横坐标向量X，纵坐标向量Y, M是被插函数在[a,b]区间上的n+1阶导数的最大值.%输出的量:n阶牛顿插值多项式L及其系数向量C，差商的矩阵A，插值余项公式wcgs及多项式(x-x0)(x-x1)...(x-xn)/(n+1)!的系数向量Cwn=length(X); A=zeros(n,n); A(:,1)=Y';s=0.0; p=1.0; q=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));endb=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1;endC=A(n,n); b=poly(X(n)); q1=conv(q1,b);湖北民族学院理学院《数值分析》课程教案教案编写人：陈以平for k=(n-1):-1:1C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k); endL(k,:)=poly2sym(C); Q=poly2sym(q1); syms Mwcgs=M*Q/c1; Cw=q1/c1;L=vpa(L,4); wcgs =vpa(wcgs,4);3．Hermite 插值法已知],[)(1b a c x f y ∈=函数及导数表nn n y y y x f y y y x f x x x x ''''1)()(01010 其中i x 互异（),,,1,0n i =寻求12+n 次多项式)(12x H n +使满足插值条件⎩⎨⎧='='=++),,1,0()()(1212n i y x H y x H i n ii n算法公式 :∑=++=nj j j j j n y x y x x H 012)')()(()(βα,其中),,1,0(),(),(n j x x j j =βα为Hermite 插值基函数⎪⎪⎪⎩⎪⎪⎪⎨⎧--∏=-=---=≠==∑≠i ji n j i i j j j j j ni ij j j x x xx x l x l x x x x l x x x x x j i 0220)()()()()()1)(21()(βα 余项为)()()(1212x H x f x R n n ++-=(22)22201()()()()(22)!n n f x x x x x x n ξ+=---+L其中(,)a b ξ∈与x 有关. 通用程序1 :function f = hermite1(x,y,y_1,x0) syms t; f = 0.0;if(length(x) == length(y))if(length(y) == length(y_1))n = length(x);elsedisp('y和y的导数的维数不相等！');return;endelsedisp('x和y的维数不相等！');return;endfor i=1:nh = 1.0;a = 0.0;for j=1:nif( j ~= i)h = h*(t-x(j))^2/((x(i)-x(j))^2);a = a + 1/(x(i)-x(j));endendf = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));if(i==n)if(nargin == 4)f = subs(f,'t',x0);elsef = vpa(f,6);endendend通用程序2 :function [Hc, Hk,wcgs,Cw]= hermite2(X,Y,Y1)%输入的量:n+1个节点(xi,yi)的横坐标向量X，纵坐标向量Y;节点处的一阶导数向量Y1%输出的量:2n+1阶埃尔米特插值多项式Hk及其系数向量Hc,误差公式wcgs及其系数向量Cw m=length(X); n=m-1;s=0; H=0;q=1;c1=1; L=ones(m,m); G=ones(1,2);for k=1:n+1V=1;for i=1:n+1if k~=is=s+(1/(X(k)-X(i)));V=conv(V,poly(X(i)))/(X(k)-X(i));endh=poly(X(k)); g=(1-2*h*s); G=g*Y(k)+h*Y1(k);湖北民族学院理学院《数值分析》课程教案教案编写人：陈以平endH=H+conv(G,conv(V,V)); b=poly(X(k));b2=conv(b,b); q=conv(q,b2); t=2*n+2;Hc=H;Hk=poly2sym (H); Q=poly2sym(q); end for i=1:tc1=c1*i; endsyms M,wcgs=M*Q/c1; Cw=q/c1;注 : 1n =时称为二点三次Hermite 插值, 可叙述为解决如下如下问题给定],[)(1b a C x f ∈函数及导数表111)(')(+++k kk k k k m m x f y y x f x x x 寻求3次多项式)(3x H 使满足插值条件:⎩⎨⎧+==)1,()(')(33k k i m x H y x H i i ii)(3x H 表达式为)()()()()()(11113x x m x m x y x y x H k k k k k k k k +++++++=ββαα, 其中⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧---=---=----+=----+=++++++++++++2))(()())(()())(21()())(21()(111211211112111k k k k k k k k k k kk k k k k k k k k k k k kx x x x x x x x x x x x x x x x x x x x x x x a x x x x x x x x x a ββ ],[1+∈k k x x x余项公式为 (4)22331()()()()()()4!k k f R x f x H x x x x x ξ+=-=--4．分段线性插值法所谓分段线性插值就是通过插值点用折线段连接起来逼近)(x f 。
计算方法分段线性_三次样条插值
计算方法分段线性_三次样条插值分段线性和三次样条插值是两种常用的插值方法，在数值分析和插值问题中广泛使用。
1.分段线性插值分段线性插值是一种简单直观的插值方法，将插值区间划分为若干个子区间，在每个子区间上用线性函数进行插值。
假设给定的插值节点有n+1 个，节点为 (x0, y0), (x1, y1), ..., (xn, yn)，并且满足 x0 <x1 < ... < xn。
则对于任意 xx 使得 x 在 [xi, xi+1] 之间，可以通过线性插值得到其函数值 yy，即：yy = yi + (xx - xi) * (yi+1 - yi) / (xi+1 - xi)分段线性插值方法简单易懂，适用于一些较简单的插值问题。
但是由于插值函数在节点之间是线性的，可能不能准确地反映出数据的特征，因此不适用于一些需要高精度的插值问题。
三次样条插值是一种更复杂、更精确的插值方法，将插值区间划分为若干个子区间，在每个子区间上用三次多项式进行插值。
三次样条插值方法的基本思想是找到一组三次多项式，满足在每个子区间内插值点的函数值和一阶导数值相等，并且两个相邻多项式在节点处的二阶导数值也相等。
具体的求解步骤如下：(1) 假设有 n+1 个插值节点 (x0, y0), (x1, y1), ..., (xn, yn)，构造 n 个三次多项式，即每个多项式在 [xi, xi+1] 之间插值。
(2) 对每个子区间内的多项式进行插值，设第 i 个子区间的多项式为 Si(x) = ai + bi(x-xi) + ci(x-xi)^2 + di(x-xi)^3、将插值节点的函数值和一阶导数值代入多项式中，可以得到 n 个线性方程，利用这 n 个线性方程可以求解出 n 个子区间的系数。
(3)由于n个子区间的多项式必须在节点处一阶导数值相等，因此再设立n-1个方程，利用这些方程可以求解出n-1个子区间的二阶导数值。
(4)将求解得到的系数和二阶导数值代入每个子区间的多项式中，得到完整的三次样条插值函数。
拉格朗日（Lagrange）插值算法
拉格朗⽇（Lagrange）插值算法拉格朗⽇插值（Lagrange interpolation）是⼀种多项式插值⽅法，指插值条件中不出现被插函数导数值，过n+1个样点，满⾜如下图的插值条件的多项式。
也叫做拉格朗⽇公式。
这⾥以拉格朗⽇3次插值为例，利⽤C++进⾏实现：1//利⽤lagrange插值公式2 #include<iostream>3using namespace std;45double Lx(int i,double x,double* Arr)6 {7double fenzi=1,fenmu=1;8for (int k=0;k<4;k++)9 {10if (k==i)11continue;12 fenzi*=x-Arr[k];13 fenmu*=Arr[i]-Arr[k];14 }15return fenzi/fenmu;16 }1718int main()19 {20double xArr[4]={};21double yArr[4]={};22//输⼊4个节点坐标23 cout<<"请依次输⼊4个节点的坐标:"<<endl;24for (int i=0;i<4;i++)25 cin>>xArr[i]>>yArr[i];2627//输⼊要求解的节点的横坐标28 cout<<"请输⼊要求解的节点的横坐标:";29double x;30 cin>>x;31double y=0;32for (int i=0;i<4;i++)33 y+=Lx(i,x,xArr)*yArr[i];34 printf("x=%lf时，y=%lf\n",x,y);3536//分界，下⾯为已知y求x37 cout<<"请输⼊要求解的节点的纵坐标:";38 cin>>y;39 x=0;40for (int i=0;i<4;i++)41 x+=Lx(i,y,yArr)*xArr[i];42 printf("y=%lf时，x=%lf\n",y,x);4344 system("pause");45return0;46 }作者：耑新新，发布于转载请注明出处，欢迎邮件交流：zhuanxinxin@。
拉格朗日分段线性三次样条正弦函数
plot(x,y)
hold on
xi=0:pi/2:2*pi;
yi=[0 1 0 -1 0];
y1=spline(xi,yi,x);
y2=y-y1;
plot(x,y1,'r',x,y2,'k')
%n=9
>> x=0:pi/50:2*pi;
y=sin(x);
plot(x,y)
hold on
xi=0:pi/4:2*pi;
yi=[0 0.7071 1 0.7071 0 -0.7071 -1 -0.7071 0];
y1=spline(xi,yi,x);
y2=y-y1;
plot(x,y1,'r',x,y2,'k')
三种实验插值比较：n=9时
x=0:pi/50:2*pi;
y=sin(x);
plot(x,y)
=
数学应用软件大型实验实验报告
实验序号：B57日期：2011年8月19日
班级
信计0901
姓名
学号
实验
名称
选择一些函数，在n个节点上用拉格朗日，分段线性，三次样条三种插值方法，计算m个插值点的函数值。通过数值和图形，将三种插值结果与精确值进行比较，适当增加n，再做比较，下列函数可选择参考;(1) y=sinx ;
plot(x,y)
hold onx1=0:pi/:2*pi;y1=[0 0.7071 1 0.7071 0 -0.7071 -1 -0.7071 0];
y2=interp1(x1,y1,x);
y3=y-y2;
plot(x,y2,'r',x,y3,'k')
函数的数值逼近-插值
课程名称计算方法实验项目名称函数的数值逼近－插值实验成绩指导老师（签名）日期2011-9-16一. 实验目的和要求1．掌握用Matlab计算Lagrange、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析。
2．通过实例学习如何用插值方法解决实际问题。
二. 实验内容和原理1）编程题2-1要求写出Matlab源程序(m文件)，并对每一行语句加上适当的注释语句；2）分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab源程序、运行结果和结果的解释、算法的分析等写在实验报告上。
2-1分析应用题用12y x=在0,1,4,9,16x=产生5个节点15,,P P。
用以下五种不同的节点构造Lagrange插值公式来计算5x=处的插值，与精确值比较并进行分析。
function y=lagr(x0,y0,x)n=length(x0);m=length(x);L=zeros(1,n);y=zeros(1,m);for k=1:ms=0;for i=1:nL(i)=1;for j=1:nif j~=iL(i)=L(i)*(x(k)-x0(j))/(x0(i)-x0(j));endends=s+y0(i)*L(i);endy(k)=s;end1） 用34,P P 构造；>> x0=[4,9]; >> y0=[2,3]; >> lagr(x0,y0,5) ans =2.20002） 用234,,P P P 构造；>> x0=[1,4,9]; >> y0=[1,2,3]; >> lagr(x0,y0,5) ans =2.26673） 用2345,,,P P P P 构造；>> x0=[1,4,9,16]; >> y0=[1,2,3,4]; >> lagr(x0,y0,5) ans =2.25404） 用1245,,,P P P P 构造；>> x0=[0,1,9,16]; >> y0=[0,1,3,4]; >> lagr(x0,y0,5) ans =2.95245） 用全部插值节点12345,,,,P P P P P 构造。
数值分析上机报告插值法
0.0522 0.0000 0.0875 0.0000 0.1375 0.0000 -0.0625 0.0000 -0.0625 0.0000 0.1375 0.0000 0.0875 0.0000 0.0522 0.0000 0.0114 0.0206 0.0274 0.0318 0.0337 0.0328 0.0291 0.0225 0.0128 0.0000
0.0794 0.1000 0.1500 0.2000 0.3500 0.5000 0.7500 1.0000 0.7500 0.5000 0.3500 0.2000 0.1500 0.1000 0.0794 0.0588 0.0568 0.0548 0.0527 0.0507 0.0486 0.0466 0.0446 0.0425 0.0405 0.0385
二、方法描述
1. Lagrange 插值 Lagrange 插值是基于基函数的插值方法，其插值多项式可以表示为：
Ln ( x) = ∑ yi li ( x)
i =0
n
其中，
( xi ) = yi f=
i 次基函数
1 ， i = 0,1, 2, , n 1 + 25 xi 2
li ( x) =
同时满足
-3.9971 0.0000 0.8397 0.0000 -0.2351 0.0000 0.0543 0.0000 0.0543 0.0000 -0.2351 0.0000 0.8397 0.0000 -3.9971 0.0000 3.9907 9.4467 16.3087 24.2507 32.5478 39.9039 44.2322 42.3774 29.7702 0.0000
L(x) 0.0385 1.2303 1.8044 1.9590 1.8458 1.5787 1.2402 0.8881 0.5604 0.2802 0.0588
常见插值算法--拉格朗日插值、三次卷积插值、三次样条插值、兰克索斯插值
常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的python代码，本⽂公式使⽤latex编写，如有错误欢迎评论指出，如果谁知道如何修改latex字号也欢迎留⾔关于⼀维、⼆维和多维插值三次卷积插值、拉格朗⽇两点插值（线性插值）、兰克索斯插值在⼆维插值时改变x和y⽅向的计算顺序不影响最终结果，这三个也是图像缩放插值时常⽤的插值算法，⽽其他插值在改变计算顺序时会产⽣明显差异，多维的情况笔者没有尝试，读者可以⾃⾏尝试或推导最近邻插值法（Nearest Neighbour Interpolation）在待求像素的四邻像素中，将距离待求像素最近的像素值赋给待求像素p_{11}p_{12}pp_{21}p_{22}python代码1def NN_interpolation(srcImg, dstH, dstW):2 scrH, scrW, _ = srcImg.shape3 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)4for i in range(dstH - 1):5for j in range(dstW - 1):6 scrX = round(i * (scrH / dstH))7 scrY = round(j * (scrW / dstW))8 dstImg[i, j] = srcImg[scrX, scrY]9return dstImg拉格朗⽇插值（Lagrange Interpolation）拉格朗⽇插值法需要找到k个p_i(x)函数，使得每个函数分别在在x_i处取值为1，其余点取值为0，则y_ip_i(x)可以保证在x_i处取值为y_i，在其余点取值为0，因此L_k(x)能恰好经过所有点，这样的多项式被称为拉格朗⽇插值多项式，记为L_k(x)=\sum_{i=1}^ky_ip_i(x)p_i(x)=\prod_{j \neq i}^{1 \leq j \leq k}\frac{x-x_j}{x_i-x_j}以四点即三次图像插值为例，因为横坐标间隔为1，则设四个点横坐标为-1、0、1和2，可得p_1(x)、p_2(x)、p_3(x)和p_4(x)假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得拉格朗⽇函数如下图所⽰，待插值点横坐标范围为[0,1]在K=2时在k=2时，也被称为线性插值通⽤公式p_1=\frac{x-x_2}{x_1-x_2}p_2=\frac{x-x_1}{x_2-x_1}\begin{align} L_2x &= p_1y_1+p_2y_2 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}pp_{21}p_{22}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_2x &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \\ &= (x_2-x)y_1+(x-x_1)y_2 \nonumber \\ &= (1-dx)y_1+dxy_2 \nonumber \\ &= (y_2-y_1)dx+y_1 \nonumber \end{align}L_2'x=y_2-y_1在K=3时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\begin{align} L_3x &= p_1y_1+p_2y_2+p_3y_3 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1+\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2+\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}p_{13}p_{21}p_{22}p_{23}pp_{31}p_{32}p_{33}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_3x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \\ &= \frac{-dx(1-dx)}{(-1)\cdot(-2)}y_1 + \frac{-(1+dx)(1-dx)}{1\cdot(-1)}y_2 + \frac{(1+dx)dx}{2\cdot 1}y_3 \nonumber \\ &= (\frac{1}{2}d^2x-\frac{1}{2}dx)y_1 - (d^2x-1)y_2 + (\frac{1}{2}d^2x+\frac{1}{2}dx)y_3 \nonumber \\ &= d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{2}y_1+\frac{1}{2}y_3)+y_2 \nonumber \end{align}L_3'x=dx(y_1-2y_2+y_3)+(\frac{1}{2}y_3-\frac{1}{2}y_1)在K=4时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}p_4=\frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}\begin{align} L_4x &= p_1y_1+p_2y_2+p_3y_3+p_4y_4 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3} {x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4\nonumber \end{align}图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_4x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4 \nonumber \\ &= \frac{dx[-(1-dx)][-(2-dx)]}{(-1)\cdot(-2)\cdot(-3)}y_1 + \frac{(1+dx)[-(1-dx)][-(2-dx)]}{1\cdot(-1)\cdot(-2)}y_2 + \frac{(1+dx)dx[-(2-dx)]}{2\cdot 1\cdot(-1)}y_3 + \frac{(1+dx)dx[-(1-dx)]}{3\cdot 2\cdot 1}y_4 \nonumber \\ &= \frac{d^3x-3d^2x+2dx}{-6}y1 + \frac{d^3x-2d^2x-dx+2}{2}y_2 + \frac{d^3x-d^2x-2dx}{-2}y_3 + \frac{d^3x-dx}{6}y_4 \nonumber \\ &= d^3x(-\frac{1}{6}y_1+\frac{1}{2}y_2-\frac{1} {2}y_3+\frac{1}{6}y_4)+d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4)+y_2 \nonumber \end{align}\begin{align} L_4'x &= d^2x(-\frac{1}{2}y_1+\frac{3}{2}y_2-\frac{3}{2}y_3+\frac{1}{2}y_4)+dx(y_1-2y_2+y_3)+(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4) \nonumber \\ &= -[\frac{1}{2}d^2x(y_1-3y_2+3y_3-y_4)-dx(y_1-2y_2+y_3)+\frac{1}{6}(2y_1+3y_2-6y_3+y_4)] \nonumber \end{align}python代码插值核计算的时候乘法和加减法计算的顺序不同可能会导致结果存在细微的差异，读者可以⾃⾏研究⼀下1class BiLagrangeInterpolation:2 @staticmethod3def LagrangeInterpolation2(x, y1, y2):4 f1 = 1 - x5 f2 = x6 result = y1 * f1 + y2 * f27return result89 @staticmethod10def LagrangeInterpolation3(x, y1, y2, y3):11 f1 = (x ** 2 - x) / 2.012 f2 = 1 - x ** 213 f3 = (x ** 2 + x) / 2.014 result = y1 * f1 + y2 * f2 + y3 * f315return result1617 @staticmethod18def LagrangeInterpolation4(x, y1, y2, y3, y4):19 f1 = - (x ** 3 - 3 * x ** 2 + 2 * x) / 6.020 f2 = (x ** 3 - 2 * x ** 2 - x + 2) / 2.021 f3 = - (x ** 3 - x ** 2 - 2 * x) / 2.022 f4 = (x ** 3 - x) / 6.023 result = y1 * f1 + y2 * f2 + y3 * f3 + y4 * f424return result2526def biLag2_2(self, srcImg, dstH, dstW):27 dstH, dstW = int(dstH), int(dstW)28 srcH, srcW, _ = srcImg.shape29 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')30 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)31for dstY in range(dstH):32for dstX in range(dstW):33for channel in [0, 1, 2]:34# p11 p1235# p36# p21 p2237# 储存为 p(y, x)38 p = [dstY * srcH / dstH, dstX * srcW / dstW]39 p11 = [math.floor(p[0]), math.floor(p[1])]40 p12 = [p11[0], p11[1] + 1]4142 p21 = [p11[0] + 1, p11[1]]43 p22 = [p21[0], p12[1]]4445 diff_y, diff_x = p[0] - p11[0], p[1] - p11[1]46 r1 = grangeInterpolation2(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel])47 r2 = grangeInterpolation2(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel])4849 c = grangeInterpolation2(diff_y, r1, r2)5051 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)52return dstImg5354def biLag3_3(self, srcImg, dstH, dstW):55 dstH, dstW = int(dstH), int(dstW)56 srcH, srcW, _ = srcImg.shape57 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')58 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)59for dstY in range(dstH):60for dstX in range(dstW):61for channel in [0, 1, 2]:62# p11 p12 p1363#64# p21 p22 p2365# p66# p31 p32 p3367# 储存为 p(y, x)68 p = [dstY * srcH / dstH, dstX * srcW / dstW]69 p22 = [math.floor(p[0]), math.floor(p[1])]70 p21 = [p22[0], p22[1] - 1]71 p23 = [p22[0], p22[1] + 1]7273 p11 = [p21[0] - 1, p21[1]]74 p12 = [p11[0], p22[1]]75 p13 = [p11[0], p23[1]]7677 p31 = [p21[0] + 1, p21[1]]78 p32 = [p31[0], p22[1]]79 p33 = [p31[0], p23[1]]8081 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]82 r1 = grangeInterpolation3(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel])83 r2 = grangeInterpolation3(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel])84 r3 = grangeInterpolation3(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel]) 8586 c = grangeInterpolation3(diff_y, r1, r2, r3)8788 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)89return dstImg9091def biLag4_4(self, srcImg, dstH, dstW):92 dstH, dstW = int(dstH), int(dstW)93 srcH, srcW, _ = srcImg.shape94 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')95 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)96for dstY in range(dstH):97for dstX in range(dstW):98for channel in [0, 1, 2]:99# p11 p12 p13 p14100#101# p21 p22 p23 p24102# p103# p31 p32 p33 p34104#105# p41 p42 p43 p44106# 储存为 p(y, x)107 p = [dstY * srcH / dstH, dstX * srcW / dstW]108 p22 = [math.floor(p[0]), math.floor(p[1])]109 p21 = [p22[0], p22[1] - 1]110 p23 = [p22[0], p22[1] + 1]111 p24 = [p22[0], p22[1] + 2]112113 p11 = [p21[0] - 1, p21[1]]114 p12 = [p11[0], p22[1]]115 p13 = [p11[0], p23[1]]116 p14 = [p11[0], p24[1]]117118 p31 = [p21[0] + 1, p21[1]]119 p32 = [p31[0], p22[1]]120 p33 = [p31[0], p23[1]]121 p34 = [p31[0], p24[1]]122123 p41 = [p21[0] + 2, p21[1]]124 p42 = [p41[0], p22[1]]125 p43 = [p41[0], p23[1]]126 p44 = [p41[0], p24[1]]127128 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]129 r1 = grangeInterpolation4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 130 r2 = grangeInterpolation4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 131 r3 = grangeInterpolation4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 132 r4 = grangeInterpolation4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 133134 c = grangeInterpolation4(diff_y, r1, r2, r3, r4)135136 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)137return dstImg三次卷积插值法（Cubic Convolution Interpolation）使⽤上图中的卷积核进⾏加权平均计算，卷积核为u(s)，四个等距（距离为1）的采样点记为x_0、x_1、x_2和x_3，采样数值记为y_0、y_1、y_2和y_3，且保证四个点均在[-2,2]区间上，计算得到g(x)，假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得三次卷积插值函数如下图所⽰，待插值点横坐标范围为[0,1]公式推导设u(s)=\begin{cases} A_1|s|^3+B_1|s|^2+C_1|s|+D_1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\because函数在s=0,1,2处连续\therefore\begin{cases} 1=u(0^+)=D_1 \\ 0=u(1^-)=A_1+B_1+C_1+D_1 \\ 0=u(1^+)=A_2+B_2+C_2+D_2 \\ 0=u(2^-)=8A_2+4B_2+2C_2+D_2 \end{cases} (1)\because函数在s=0,1,2处导函数连续\therefore\begin{cases} u'(0^-)=u'(0+) \\ u'(1^-)=u'(1+) \\ u'(2^-)=u'(2+)\end{cases} \Rightarrow \begin{cases} -C_1=C_1 \\ 3A_1+2B_1+C_1=3A_2+2B_2+C_2\\ 12A_2+4B_2+C+2=0 \end{cases} ~~~~ (2)联⽴⽅程组(1)(2)，设A_2=a，解得\begin{cases} A_1=a+2 \\ B_1=-(a+3) \\ C_1=0 \\ D_1=1 \\ A_2=a \\ B_2=-5a \\ C_2=8a \\ D_2=-4a \end{cases}\Rightarrow u(s)=\begin{cases} (a+2)|s|^3-(a+3)|s|^2+1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2\\ 1, &s=0 \\ 0, &otherwise \end{cases}\because g(x)=\sum_kC_ku(s+j-k), ~~~~k=j-1,j, j+1,j+2且0<s<1⼜\because \begin{cases}\begin{align} u(s+1)&=as^3-2as^2+as \nonumber \\ u(s)&=(a+2)s^3-(a+3)s^2+1 \nonumber \\ u(s-1)&=-(a+2)s^3+(2a+3)s^2-as \nonumber \\ u(s-2)&=-as^3+as^2 \nonumber \end{align}\end{cases}\begin{align} \therefore g(x) &= C_{j-1}u(s+1)+C_{j}u(s)+C_{j+1}u(s-1)+C_{j+2}u(s-2) \nonumber \\ &= C_{j-1}(as^3-2as^2+as)+C_j[(a+2)s^3-(a+3)s^2+1]+C_{j+1}[-(a+2)s^3+ (2a+3)s^2-as]+C_{j+2}[-a^3+as^2] \nonumber \\ &= s^3[aC_{j-1}+(a+2)C_j-(a+2)C_{j+1}-aC_{j+2}]+s^2[-2aC_{j-1}-(a+3)C_j+(2a+3)C_{j+1}+aC_{j+2}]+s[aC_{j-1}-aC_{j+1}]+C_j \nonumber \end{align} ~~(3)f在x_j处泰勒展开得到f(x)=f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+\cdots\therefore \begin{cases} f(x_{j+1})=f(x_j)+f'(x_j)(x_{j+1}-x_j)+\frac{1}{2}f''(x_j)(x_{j+1}-x_j)^2+\cdots \\ f(x_{j+2})=f(x_j)+f'(x_j)(x_{j+2}-x_j)+\frac{1}{2}f''(x_j)(x_{j+2}-x_j)^2+\cdots \\ f(x_{j-1})=f(x_j)+f'(x_j)(x_{j-1}-x_j)+\frac{1}{2}f''(x_j)(x_{j-1}-x_j)^2+\cdots \end{cases}令x_{j+1}-x_j=h\because x_{i+1}=x_i+1\therefore x_{j+2}-x_j=2h，x_{j-1}-x_j=-h\therefore \begin{cases} f(x_{j+2})=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+\cdots \\ f(x_{j+1})=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \\ f(x_{j-1})=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \end{cases}\therefore \begin{cases} c_{j-1}=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3) \\ c_j=f(x_j) \\ c_{j+1}=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3)\\ c_{j+2}=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+o(h^3) \end{cases} ~~ (4)将(4)代⼊(3)，得g(x)=-(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3+[(6a+3)hf'(x_j)+\frac{4a+3}{2}h^2f''(x_j)]s^2-2ahf'(x_j)s+f(x_j)+o(h^3)\because h=1,s=x-x_J\therefore sh=x-x_j\begin{align}\therefore f(x)&= f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+o(h^3) \nonumber \\ &= f(x_j)+f'(x_j)sh+\frac{1}{2}f''(x_j)s^2h^2+o(h^3) \nonumber \end{align}\therefore f(x)-g(x)=(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3-(2a+1)[3hf'(x_j)+h^2f''(x_j)]s^2+[(2a+1)hf'(x_j)]s+o(h^3)\because 期望f(x)-g(x)趋于0\therefore 2a+1=0 \Rightarrow a=-\frac{1}{2}\therefore u(s)=\begin{cases} \frac{3}{2}|s|^3-\frac{5}{2}|s|^2+1, &0<|s|<1 \\ -\frac{1}{2}|s|^3+\frac{5}{2}|s|^2-4|s|+2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\therefore g(s)=s^3[-\frac{1}{2}c_{j-1}+\frac{3}{2}c_j-\frac{3}{2}c_{j+1}+\frac{1}{2}c_{j+2}]+s^2[c_{j-1}-\frac{5}{2}c_j+2c_{j+1}-\frac{1}{2}c_{j+2}]+s[-\frac{1}{2}c_{j-1}+\frac{1} {2}c_{j+1}]+c_j图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}python代码1class BiCubicConvInterpolation:2 @staticmethod3def CubicConvInterpolation1(p0, p1, p2, p3, s):4# ⽤g（s）公式计算，已经将四个u（s）计算完毕并整理5# as^3 + bs^2 + cs + d6 a = 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3)7 b = 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3)8 c = 0.5 * (-p0 + p2)9 d = p110return d + s * (c + s * (b + s * a))1112 @staticmethod13def CubicConvInterpolation2(s):14# ⽤u（s）公式计算15 s = abs(s)16if s <= 1:17return 1.5 * s ** 3 - 2.5 * s ** 2 + 118elif s <= 2:19return -0.5 * s ** 3 + 2.5 * s ** 2 - 4 * s + 220else:21return 02223def biCubic1(self, srcImg, dstH, dstW):24# p11 p12 p13 p1425#26# p21 p22 p23 p2427# p28# p31 p32 p33 p3429#30# p41 p42 p43 p4431 dstH, dstW = int(dstH), int(dstW)32 scrH, scrW, _ = srcImg.shape33 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')34 dstImg = np.zeros((dstH, dstW, 1), dtype=np.uint8)35for dstY in range(dstH):36for dstX in range(dstW):37for channel in [0]:38 y = dstY * scrH / dstH39 x = dstX * scrW / dstW40 y1 = math.floor(y)41 x1 = math.floor(x)4243 array = []44for i in [-1, 0, 1, 2]:45 temp = self.CubicConvInterpolation1(srcImg[y1 + i, x1 - 1, channel],46 srcImg[y1 + i, x1, channel],47 srcImg[y1 + i, x1 + 1, channel],48 srcImg[y1 + i, x1 + 2, channel],49 x - x1)50 array.append(temp)5152 temp = self.CubicConvInterpolation1(array[0], array[1], array[2], array[3], y - y1)53 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)5455return dstImg5657def biCubic2(self, srcImg, dstH, dstW):58# p11 p12 p13 p1459#60# p21 p22 p23 p2461# p62# p31 p32 p33 p3463#64# p41 p42 p43 p4465 dstH, dstW = int(dstH), int(dstW)66 scrH, scrW, _ = srcImg.shape67 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')68 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)69for dstY in range(dstH):70for dstX in range(dstW):71for channel in [0, 1, 2]:72 y = dstY * scrH / dstH73 x = dstX * scrW / dstW74 y1 = math.floor(y)75 x1 = math.floor(x)7677 array = []78for i in [-1, 0, 1, 2]:79 temp = 080for j in [-1, 0, 1, 2]:81 temp += srcImg[y1 + i, x1 + j, channel] * self.CubicConvInterpolation2(x - (x1 + j))82 array.append(temp)8384 temp = 085for i in [-1, 0, 1, 2]:86 temp += array[i + 1] * self.CubicConvInterpolation2(y - (y1 + i))87 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)8889return dstImg三次样条插值在n-1个区间上寻找n-1个三次曲线，使其满⾜相邻曲线在端点处值相等、⼀阶导数相等，⼆阶导数相等，在加以边界条件后可得每个曲线的⽅程，然后沿x轴依次偏移对应的距离即可得到插值结果，如仅需要特定范围内的结果，则可以⼤幅减少计算量公式推导设S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3, ~~~~i=0,1,...,n-1则 \begin{cases} S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\\ S_i''(x)=2c_i+6d_i(x-x_i)\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1设h_i(x)=x_{i+1}-x_i，可得\begin{cases} S_i(x)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3\\ S_i'(x)=b_i+2c_ih_i+3d_ih_i^2\\ S_i''(x)=2c_i+6d_ih_i\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1\because S_i(x)过点(x_i,y_i)\therefore S_i(x)=a_i=y+i, ~~~~i=0,1,...,n-1 ~~~~~~(1)\because S_i(x)与S_{i+1}(x)在X_{i+1}处相等\therefore S_i(x_{i+1})=S_{i+1}(x_{i+1})\Rightarrow a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}, ~~~~i=0,1,...,n-2~~~~~~(2)\because S_i'(x)与S_{i+1}'(x)在X_{i+1}处相等\therefore S_i'(x)-S_{i+1}'(x)=0\Rightarrow b_i+2c_ih_i+3d_ih_i^2-b_{i+1}=0~~~~~~(3)\because S_i''(x)与S_{i+1}''(x)在X_{i+1}处相等\therefore S_i''(x)-S_{i+1}''(x)=0\Rightarrow 2c_i+6d_ih_i-2c_{i+1}=0, ~~~~i=0,1,...,n-2~~~~~~(4)设m_i=S_i(x_i)=2c_i，即c_i=\frac{1}{2}m_i, ~~~~i=0,1,...,n-1~~~~~~(5)将(5)代⼊(4)，得2c_i+6d_ih_i-2c_{i+1}=0\Rightarrow m_i+6h_id_i-m_{i+1}=0\Rightarrow d_i=\frac{m_{i+1}-m_i}{6h_i}, ~~~~i=0,1,...,n-2~~~~~~(6)将(1)(5)(6)代⼊(2)，得\begin{align} &a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1} \nonumber \\ \Rightarrow&y_i+b_ih_i+\frac{1}{2}m_ih_i^2+\frac{m_{i+1}-m_i}{6h_i}h_i^3=y_{i+1} \nonumber \\\Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{2}m_ih_i-\frac{1}{6}(m_{i+1}-m_i)h_i \nonumber \\ \Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i, ~~~~i=0,1,...,n-2~~~~~~(7) \nonumber \end{align}将(5)(6)(7)代⼊(3)，得\begin{align} &\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+2\cdot\frac{1}{2}m_ih_i+3\frac{m_{i+1}-m_i}{6h_i}h_i^2-(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{1}{3}m_{i+1}h_{i+1}-\frac{1}{6}m_{i+2}h_{i+1})=0 \nonumber \\ \Rightarrow&\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+m_ih_i+\frac{1}{2}(m_{i+1}-m_i)h_i-\frac{y_{i+2}-y_{i+1}}{h_{i+1}}+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=0 \nonumber \\ \Rightarrow&m_ih_i(-\frac{1}{3}+1-\frac{1}{2})+m_{i+1}h_i(-\frac{1}{6}+\frac{1} {2})+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&\frac{1}{6}(m_ih_i+2m_{i+1}h_i+2m_{i+1}h_{i+1}+m_{i+2}h_{i+1})=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&m_ih_i+2m_{i+1}(h_i+h_{i+1})+m_{i+2}h_{i+1}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}), ~~~~i=0,1,...,n-2~~~~~~(8) \nonumber \end{align}由(8)可知i=0,1,...,n-2，则有m_0,m_1,...,m_n，需要两个额外条件⽅程组才有解⾃然边界(Natural)m_0=0，m_n=0\begin{bmatrix} \tiny 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}固定边界(Clamped)\begin{align} &\begin{cases} S_0'(x_0)=A\\ S_{n-1}'(x_n)=B \end{cases} \nonumber \\ \Rightarrow&\begin{cases} b_0=A\\ b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=B\end{cases} \nonumber \\ \Rightarrow&\begin{cases} A=\frac{y_1-y_0}{h_0}-\frac{h_0}{2}m_0-\frac{h_0}{6}(m_1-m_0)\\ B=\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{1}{3}m_{n-1}h_{n-1}+m_{n-1}h_{n-1}+\frac{1}{2}m_nh_{n-1}-\frac{1}{2}m_{n-1}h_{n-1} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 2h_0m_0+h_0m_1=6(\frac{y_1-y_0}{h_0}-A)\\ h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6(B-\frac{y_n-y_{n-1}}{h_{n-1}}) \end{cases} \nonumber \\ \end{align}\begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 1 & 2 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ B-\frac{y_n-y_{n-1}}{h_{n-1}} \end{bmatrix}⾮节点边界(Not-A-Knot)\begin{align} &\begin{cases} S_0'''(x_1)=S_1'''(x_1)\\ S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 6\cdot\frac{m_1-m_0}{6h_0}=6\cdot\frac{m_2-m_1}{6h_1}\\ 6\cdot\frac{m_{n-1}-m_{n-2}}{6h_{n-2}}=6\cdot\frac{m_n-m_{n-1}}{6h_{n-1}} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} -h_1m_0+(h_1+h_0)m_1-h_0m_2=0\\ -h_{n-1}m_{n-2}+(h_{n-1}+h_{n-2})m_{n-1}-h_{n-2}m_n=0 \end{cases} \nonumber \\ \end{align}\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 &2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & -h_{n-1} & h_{n-1}+h_{n-2} & -h_{n-2} \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}在n=4时通⽤公式⾃然边界\begin{bmatrix} 1 & 0 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}固定边界\begin{bmatrix} 2 & 1 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ B-\frac{y_3-y_2}{h_2} \end{bmatrix}⾮节点边界\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & -h_2 & h_1+h_2 & -h_1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}图像插值x_{i+1}-x_i=1 \Rightarrow h_i(x)=1在n=4时，即四个点时如下所⽰p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}⾃然边界（可⽤TDMA或化简计算）\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}固定边界（只能⽤TDMA计算）\begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} y_1-y_0-A\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ y_2-y_3+B \end{bmatrix}⾮节点边界（只能化简计算）\begin{bmatrix} -1 & 2 & -1 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & -1 & 2 & -1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}python代码1class BiSplineInterpolation:2 @staticmethod3 def TDMA(a, b, c, d):4 n = len(d)56 c[0] = c[0] / b[0]7 d[0] = d[0] / b[0]89for i in range(1, n):10 coef = 1.0 / (b[i] - a[i] * c[i - 1])11 c[i] = coef * c[i]12 d[i] = coef * (d[i] - a[i] * d[i - 1])1314for i in range(n - 2, -1, -1):15 d[i] = d[i] - c[i] * d[i + 1]1617return d1819 @staticmethod20 def Simplified_Natural4(y1, y2, y3, y4):21 # 四点⾃然边界化简公式22 d1 = y1 + y3 - 2 * y223 d2 = y2 + y4 - 2 * y32425 k0 = 026 k1 = (4 * d1 - d2) * 0.427 k2 = (4 * d2 - d1) * 0.428 k3 = 02930return [k0, k1, k2, k3]3132 @staticmethod33 def Simplified_Not_A_Knot4(y1, y2, y3, y4):34 # 四点⾮节点边界化简公式35 d1 = y1 + y3 - 2 * y236 d2 = y2 + y4 - 2 * y33738 k0 = 2 * d1 - d239 k1 = d140 k2 = d241 k3 = 2 * d2 - d14243return [k0, k1, k2, k3]4445 # TDMA矩阵说明46 # a0 和 c3 没有实际意义，占位⽤47 # a0 [b0 c0 00 ] [x0] [d0]48 # [a1 b1 c1 0 ] [x1] = [d1]49 # [0 a2 b2 c2] [x2] [d2]50 # [00 a3 b3] c3 [x3] [d3]5152 def SplineInterpolationNatural4(self, x, y1, y2, y3, y4):53 # ⽤TDMA计算54 # matrix_a = [0, 1, 1, 0]55 # matrix_b = [1, 4, 4, 1]56 # matrix_c = [0, 1, 1, 0]57 # matrix_d = [0, 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 0]58 # matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)5960 # 化简计算61 matrix_x = self.Simplified_Natural4(y1, y2, y3, y4)6263 a = y264 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.065 c = matrix_x[1] / 2.066 d = (matrix_x[2] - matrix_x[1]) / 6.06768 s = a + b * x + c * x * x + d * x * x * x69return s7071 def SplineInterpolationClamped4(self, x, y1, y2, y3, y4):72 # 仅有TDMA计算，⽆法化简73 A, B = 1, 17475 matrix_a = [0, 1, 1, 1]76 matrix_b = [2, 4, 4, 2]77 matrix_c = [1, 1, 1, 0]78 matrix_d = [6 * (y2 - y1 - A), 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 6 * (B - y4 + y3)]79 matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)8081 a = y282 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.083 c = matrix_x[1] / 2.084 d = (matrix_x[2] - matrix_x[1]) / 6.08586 s = a + b * x + c * x * x + d * x * x * x87return s8889 def SplineInterpolationNotAKnot4(self, x, y1, y2, y3, y4):90 # ⽆法使⽤TDMA计算91 matrix_x = self.Simplified_Not_A_Knot4(y1, y2, y3, y4)9293 a = y294 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.095 c = matrix_x[1] / 2.096 d = (matrix_x[2] - matrix_x[1]) / 6.09798 s = a + b * x + c * x * x + d * x * x * x99return s100101 def biSpline4(self, srcImg, dstH, dstW):102 dstH, dstW = int(dstH), int(dstW)103 srcH, srcW, _ = srcImg.shape104 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')105 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)106for dstY in range(dstH):107for dstX in range(dstW):108for channel in [0, 1, 2]:109 # p11 p12 p13 p14110 #111 # p21 p22 p23 p24112 # p113 # p31 p32 p33 p34114 #115 # p41 p42 p43 p44116 # 储存为 p(y, x)117 p = [dstY * srcH / dstH, dstX * srcW / dstW]118 p22 = [math.floor(p[0]), math.floor(p[1])]119 p21 = [p22[0], p22[1] - 1]120 p23 = [p22[0], p22[1] + 1]121 p24 = [p22[0], p22[1] + 2]122123 p11 = [p21[0] - 1, p21[1]]124 p12 = [p11[0], p22[1]]125 p13 = [p11[0], p23[1]]126 p14 = [p11[0], p24[1]]127128 p31 = [p21[0] + 1, p21[1]]129 p32 = [p31[0], p22[1]]130 p33 = [p31[0], p23[1]]131 p34 = [p31[0], p24[1]]132133 p41 = [p21[0] + 2, p21[1]]134 p42 = [p41[0], p22[1]]135 p43 = [p41[0], p23[1]]136 p44 = [p41[0], p24[1]]137138 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]139 r1 = self.SplineInterpolationNatural4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 140 r2 = self.SplineInterpolationNatural4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 141 r3 = self.SplineInterpolationNatural4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 142 r4 = self.SplineInterpolationNatural4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 143144 c = self.SplineInterpolationNatural4(diff_y, r1, r2, r3, r4)145146 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)。
插值法公式
插值法公式1. 什么是插值法？插值法是一种通过已知数据点之间的曲线进行估算或推测的数值方法。
它可以用来估计缺失点的数值，或者通过已知数据点之间的曲线来做出预测。
插值法在数学、统计学、计算机科学和工程等领域都有广泛的应用。
2. 常用的插值法在插值法中，有多种算法可供选择，下面介绍几种常用的插值法。
2.1 线性插值法线性插值法是一种简单但常用的插值法。
它假设两点之间的曲线是一条直线，根据已知的两个点（x₁, y₁）和（x₂, y₂）之间的线性关系，可以推断出任意两点之间的数值。
线性插值法的公式如下：y = y₁ + (y₂ - y₁) / (x₂ - x₁) * (x - x₁)其中，y是待估算的数值，x是已知的数据点。
2.2 拉格朗日插值法拉格朗日插值法是一种常用的多项式插值法。
它利用已知的数据点构造一个多项式，并通过该多项式来估算任意点的数值。
拉格朗日插值法的公式如下：L(x) = ∑[i=0~n] yᵢ * Lᵢ(x)其中，L(x)表示估算值，yᵢ表示已知数据点的y值，Lᵢ(x)表示拉格朗日基函数，定义如下：Lᵢ(x) = ∏[j=0~n, j≠i] (x - xₓ₊₀₋₀ⱼ) / (xₓ₊₀₋₀ᵢ - xₓ₊₀₋₀ⱼ)在这里，n是已知数据点的数量，xₓ₊₀₋₀ⱼ是第j个已知数据点的x值。
2.3 三次样条插值法三次样条插值法是一种更复杂的插值方法，它利用三次多项式来逼近已知数据点之间的曲线。
三次样条插值法的公式如下：S(x) = aⱼ(x - xₓ₊₂₋₂)³ + bⱼ(x - xₓ₊₂₋₂)² + cⱼ(x - xₓ₊₂₋₂) + dⱼ其中，S(x)表示估算值，aⱼ、bⱼ、cⱼ和dⱼ是通过已知数据点计算得到的系数。
3. 插值法的应用插值法在很多领域都有广泛的应用。
下面列举几个常见的应用场景：•图像处理：在图像处理中，插值法可以用来放大或缩小图像，通过已有像素点之间的颜色值来估算新的像素点的颜色值。
曲面插值算法
曲面插值算法
曲面插值算法（surface interpolation algorithm）是指根据给定的一组离散点数据，通过某种数学方法来拟合出一个连续的曲面模型的算法。
常用的曲面插值算法有以下几种：1. 三次样条插值算法（Cubic spline interpolation）：该算法基于三次多项式形式的曲线来实现曲面的插值。
它通过满足一些额外的条件（如节点间的光滑要求）来获得平滑的插值曲线。
2. Lagrange插值算法：该算法使用Lagrange插值多项式来拟合曲面。
Lagrange插值多项式是通过使用给定数据点上的拉格朗日插值基函数的线性组合来定义的。
3. 三角网格插值算法（Triangulated surface interpolation）：该算法使用一组三角形来构建曲面模型。
它通常通过将给定的离散点数据连接起来形成一个三角网格，并在每个三角形中使用线性插值来计算曲面上的其他点。
4. 回归分析算法（Regression analysis）：该算法使用回归分析方法来建立一个曲面模型。
它通过拟合某种确定性回归方程，使曲面与给定的离散点数据最符合。
这些算法在实际应用中都有各自的优缺点，并且适用于不同类型的曲面插值问题。
在选择曲面插值算法时需要根据实际问题的特点和要求，综合考虑各种因素来进行选择。
三次样条插值
一、问题提出
为给定的节点， 设 x0 , x1 xn 为给定的节点，yi = f ( xi ) ，i = 0,1, n 为相应的函数值， 为相应的函数值，求一个次数不超过 n 的多项式 Pn (x)， 使其满足
Pn ( xi ) = yi，
i = 0,1, n .
这类问题称为插值问题。 称为被插值函数 P 被插值函数， 这类问题称为插值问题。 f ( x) 称为被插值函数， n ( x) 称 插值问题 插值函数， 称为插值节点 为插值函数， x0 , x1 xn 称为插值节点
六、 分段插值
所谓分段插值，就是将被插值函数逐段多项式化。 所谓分段插值，就是将被插值函数逐段多项式化。在每 个 [ xi , xi +1 ] 子段上构造插值多项式，然后把它们装配在一， 子段上构造插值多项式，然后把它们装配在一， 作为整个区间 [ a, b ] 上的插值函数，即称为分段多项式。如果 上的插值函数，即称为分段多项式。 次式， 函数 Sk ( x ) 在每个子段上都是 k 次式，则称为 一般（低次： 一般（低次：k=1,2,3） ） 次式。 k 次式。
f [ x0 , x1 ] = 5, f [ x0 , x1, x2 , x3 ] = 1,
N n ( x) = 0 5( x 1) + 2( x 1)( x 2)
+ ( x 1)( x 2)( x 3)
= x3 4 x + 3
五、 Hermite插值多项式 插值多项式
给定的是节点上的函数值和导数值 问题： 问题：已知
∑
i=0
y i li ( x )
( x x0 ) ( x xi 1 )( x xi +1 ) ( x xn ) , i = 0,1, n ( xi x0 ) ( xi xi 1 )( xi xi +1 ) ( xi xn )
常见的插值方法及其原理
常见的插值方法及其原理插值是指在已知数据点的情况下，根据其中一种规则或算法，在这些数据点之间进行预测或估计。
常见的插值方法有：拉格朗日插值、牛顿插值、分段线性插值、样条插值和Kriging插值等。
1.拉格朗日插值方法：拉格朗日插值是一种基于多项式的插值方法。
它假设已知数据点的函数曲线可以由一个多项式来表示。
拉格朗日插值的原理是，通过确定多项式的系数，使多项式在已知数据点上满足给定的函数值。
具体地说，对于给定的一组已知数据点和对应的函数值，拉格朗日插值方法通过构造一个多项式，使得该多项式在每个数据点上的函数值等于给定的函数值。
然后，通过该多项式在插值点上的函数值来估计未知数据点的函数值。
2.牛顿插值方法：牛顿插值也是一种基于多项式的插值方法，其原理类似于拉格朗日插值。
它也是通过确定多项式的系数，使多项式在已知数据点上满足给定的函数值。
不同的是，牛顿插值使用了差商的概念，将插值多项式表示为一个累次求和的形式。
具体地说，对于给定的一组已知数据点和对应的函数值，牛顿插值方法通过差商的计算，得到一个多项式表达式。
然后，通过该多项式在插值点上的函数值来估计未知数据点的函数值。
3.分段线性插值方法：分段线性插值是一种简单而常用的插值方法。
它假设在两个相邻已知数据点之间的曲线是一条直线。
分段线性插值的原理是，通过连接相邻数据点之间的线段，构造一个连续的曲线。
具体地说，对于给定的一组已知数据点和对应的函数值，分段线性插值方法将曲线划分为若干小段，每一小段都是一条直线。
然后，在每个数据点之间的区域上，通过线性插值来估计未知数据点的函数值。
4.样条插值方法：样条插值是一种基于插值条件和光滑条件的插值方法。
它假设在两个相邻已知数据点之间的曲线是一个低次数的多项式。
样条插值的原理是，通过确定各个数据点之间的插值多项式系数，使得整个曲线在插值点上的各阶导数连续。
具体地说，对于给定的一组已知数据点和对应的函数值，样条插值方法将曲线划分为若干小段，每一小段都是一个低次数的多项式。
实验一拉格朗日插值、分段线性插值...
实验一拉格朗日插值、分段线性插值、三次样条插值的比较一、问题提出选择函数y=exp(-x2) (-2≤x≤2),在n个节点上（n不要太大，如5~11）用拉格朗日、分段线性、三次样条三种插值方法，计算m个插值点的函数值（m要适中，如50~100）。
二、要求通过数值和图形输出，将三种插值结果与精确值进行比较。
适当增加n，在作比较，由此作初步分析。
三、问题的解答为统一起见，认为题目中的节点指的是中间节点，插值时还要考虑两端点的函数值，故原题目中有 n+2 个插值节点。
以下给出题目所需的 MATLAB 函数，其中参数 count_knot 表示题目中的n ，count_dot 表示题目中的m 。
function result=campare3inter(count_knot,count_dot)%count_knot 中间节点的个数%count_dot 拟合的函数值的个数clfknot=linspace(-2,2,count_knot+2);x=-2:0.01:2;y=exp(-x.^2);y0=exp(-knot.^2);plot(x,y,knot,y0,'ro');%,hold on;x_new=linspace(-2,2,count_dot);y_real=exp(-x_new.^2);%Lagrange 插值y_lagrange=lagrange(knot,y0,x_new);plot(x_new,y_lagrange,'k');hold on;%分段线性插值y_line=zeros(1,length(x_new))count_1=1;for j=1:count_dotfor i=1:count_knot+1if((x_new(j)>=knot(i))&((x_new(j)<=knot(i+1))))%直线的点斜式方程y_line(count_1)=((y0(i)-y0(i+1))/(knot(i)-knot(i+1)))*(x_new(j)-knot(i))+y0(i);count_1=count_1+1;break;end endendplot(x_new,y_line,'b');hold on;%三次样条插值 y_temp=[0 y0 0];pp=csape(knot,y_temp,'second');[breaks,coefs,npolys,ncoefs,dim]=unmkpp(pp);count=1;for j=1:count_dotfor i=1:count_knot+1if((x_new(j)>=knot(i))&((x_new(j)<=knot(i+1))))y_spline(count)=polyval(coefs(i,:),x_new(j)-knot(i));count=count+1;break;end endendplot(x_new,y_spline,'g');%输出原函数值和三种插值函数值的比较结果result=[y_real' y_lagrange' y_line' y_spline']图形输出（n=5,m=50）10.90.80.70.60.50.40.30.20.1-2 -1.5 -1 -0.5 0 0.5 1 1.5 2绿色：节点和exp(-x2)。
数值分析常用的插值方法
数值分析常用的插值方法数值分析中常用的插值方法有线性插值、拉格朗日插值、分段线性插值、Newton插值、Hermite插值、样条插值等。
下面将对这些插值方法进行详细介绍。
一、线性插值（linear interpolation）线性插值是最简单的插值方法之一、假设已知函数在两个点上的函数值，通过这两个点之间的直线来估计中间点的函数值。
线性插值公式为：f(x)=f(x0)+(x-x0)*(f(x1)-f(x0))/(x1-x0)其中，f(x)表示要求的插值点的函数值，f(x0)和f(x1)是已知的两个点上的函数值，x0和x1是已知的两个点的横坐标。
二、拉格朗日插值（Lagrange interpolation）拉格朗日插值是一种基于多项式的插值方法。
它通过多个已知点的函数值构造一个多项式，并利用这个多项式来估计其他点的函数值。
拉格朗日插值多项式的一般形式为：f(x) = Σ[f(xi) * Li(x)] (i=0,1,2,...,n)其中，f(x)表示要求的插值点的函数值，f(xi)是已知的多个点的函数值，Li(x)是拉格朗日基函数。
拉格朗日基函数的表达式为：Li(x) = Π[(x-xj)/(xi-xj)] (i≠j，i,j=0,1,2,...,n)三、分段线性插值（piecewise linear interpolation）分段线性插值是一种逐段线性近似函数的方法。
通过将整个插值区间分成多个小段，在每个小段上使用线性插值来估计函数的值。
分段线性插值的过程分为两步：首先确定要插值的点所在的小段，在小段上进行线性插值来估计函数值。
四、Newton插值（Newton interpolation）Newton插值也是一种基于多项式的插值方法。
利用差商的概念来构造插值多项式。
Newton插值多项式的一般形式为：f(x)=f(x0)+(x-x0)*f[x0,x1]+(x-x0)*(x-x1)*f[x0,x1,x2]+...其中，f(x)表示要求的插值点的函数值，f(x0)是已知的一个点的函数值，f[xi,xi+1,...,xi+k]是k阶差商。
三种插值方法比较
17世界后牛顿,拉格朗日分别讨论了等距和非等距的一般插值公式.在近代,插值法仍然是数据处理和编制函数表的常用工具，又是数值积分、数值微分、非线性方程求根和微分方程数值解法的重要基础,许多求解计算公式都是以插值为基础导出的。
三种插值方法的比较：拉格朗日插值、分段线性插值与三次样条插值三种插值法在处理问题时的比较。
插值问题的提法是：已知f(x)(可能未知或非常复杂函数)在彼此不同的n+1个实点 0x ,1x ,…n x 处的函数值是f(0x )，f(1x )，…，f(n x )，这时我们简单的说f(x)有n+1个离散数据对｛（i x ，i y ）｝i n ＝0．要估算f(x)在其它点ｘ处的函数值，最常见的一种办法就是插值，即寻找一个相对简单的函数y(x)，使其满足下列插值条件：y (i x )＝f (i x )，i=0，1，…，n ．并以y (x)作为f (x)的近似值．其中y (x)称为插值函数，f (x)称为被插函数．[1,2,3] 选用不同类型的插值函数，逼近的效果不同，下面给出拉格朗日多项式插值、 分段线性插值及三次样条插值在处理问题时的应用比较分析．多项式插值是最常见的一种函数插值．在一般插值问题中，由插值条件可以唯一确定一个次数不超过ｎ的插值多项式满足上述条件．从几何上看可以理解为：已知平面上n+1个不同点，要寻找一条次数不超过ｎ的多项式曲线通过这些点．插值多项式一般有两种常见的表达形式，一个是拉格朗日(Lagrange)插值多项式，另一个是牛顿(Newton)插值多项式．且 Ｌａｇｒａｎｇｅ插值公式恒等于Ｎｅｗｔｏｎ插值公式．分段线性插值与三次样条插值可以避免高次插值可能出现的大幅度波动现象(龙格现象)，在实际应用中通常采用分段低次插值来提高近似程度，比如可用分线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的总体光滑性较差．为了克服这一缺点，一种全局化的分段插值方法———三次样条插值成为比较理想的工具。
计算方法分段线性_三次样条插值
0x
1
设 | f(4)(x) | 1, x (0,1)
则3次插值多项式的误差
|
R(x)
||
f(4)(ξ) (x 4!
x0 )(x
x1)(x
x2 )(x
x3 )
|
41！| (x x0 )(x x1)(x x2 )(x x3 ) |
41！
1 3
1 3
2 3
1
S(x),S' (x),S' ' (x)
在a, b 上连续)
三次样条函数插值
则称S(x)为三次样条插值函数。
比线性插值要求严 苛得多！
y
y=f(x) (黑色)
y=S(x) （蓝色）
x0
xi
xi+1
xn
x
问题：这个S(x)是三次样条插值函数吗？
回答：不是。 原因： 1）在每个区间内都是3次多项式； 2）在小区区间端点xi处连续并且：f(xi)=S(xi) ； 3）S”(x)在xi 点不连续，即不光滑。
x0 )(x
x1)
(x
xn )，
ξ (a, b)
例：设函数f(x)定义在区间[0, 1]上 ，并且满足 |f(4) (x)|<1，xϵ[0, 1]。在4个插值节点 x0=0, x1=1/3, x2=2/3, x3=1, 对f(x)进行插值得多项式P3(x)，估计误差。
下面讨论误差的情况：
x 0
1 f(0) 1
x 1
0 0
f(1)
1 x
x 2