thiele插值算法

1点插值算法

function [C,c]=thiele(X,Y,Z)%X为插值点横坐标，Y为插值点纵坐标,Z为待插值点横坐标.输出为逆商系数和待插点纵坐标. C=nichashang(X,Y); c=display(Z,C,X) end function c=nichashang(a,b)%计算逆插商的值 [x,y]=size(a); d=zeros(y,y); c=ones(1,y) for j=1:y if(j==1) for i=1:y d(i,1)=b(1,i); end else for i=j:y d(i,j)=(a(1,i)-a(1,j-1))/(d(i,j-1)-d(j-1,j-1)); end end end for k=1:y c(1,k)=d(k,k); end end function c=display(X,Y,Z)%x为待求点值，y为thiele多项式系数，z为原插值点值 [x,y]=size(Y); [x1,y1]=size(X); c=ones(1,y1); for i=1:y1 A=ones(y+1,1); A(1,1)=1;A(2,1)=Y(1,1); B=ones(y+1,1); B(1,1)=0;B(2,1)=1; for k=3:y+1 A(k,1)=A(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*A(k-2,1); B(k,1)=B(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*B(k-2,1); end c(1,i)=A(y+1,1)/B(y+1,1); end end
示例
插值函数为sin(x)的四个点：pi/6 pi/4 pi/3 pi/2
待求点为Z=pi/5 pi/5.5 pi/6.5
x =

0.5236    0.7854    1.0472    1.5708

y =

0.5000    0.7071    0.8660    1.0000

输出结果为：
逆商： 0.5000 1.2641 1.5731 -0.8348
插值点：0.5881 0.5409 0.4643

牛顿插值：

function [H,Z]=newtoncz(x,y,z)%x,y为插值坐标矩阵，z为待插值点 x=x';y=y';z=z';%转致 hangx=size(x,1);%获取x的行数; H=ones(hangx,hangx);%生成行数等于x的矩阵，用来存取系数 for j=1:hangx%这个循环是算法核心，即那个倒三角的表 if(j==1) H(:,1)=y; else for i=j:hangx H(i,j)=(H(i,j-1)-H(i-1,j-1))/(x(i,1)-x(i-j+1,1)); end end end %输出系数 fprintf('系数为：'); for i=1:hangx fprintf('%f',H(i,i)); end fprintf('\n'); %计算待插点的值; hangz=size(z,1); Z=ones(hangz,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+H(k,k)*S; S=S.*(z-x(k,1)); end end end
用法和示例同上：
下次分析这两种插值的特点和优势。
基于一维牛顿插值的图像放缩：
function A=tuxiangchazhi(Y,apha)%进行图像的一维插值运算，基于thiele和netwon两种算法。X为灰色图像矩阵图像，apha为缩放系数，暂时假设aoha为2; [x,y]=size(Y); X=im2double(Y); x1=x*apha;y1=y*apha; A=zeros(x1,y1);%生成插值后矩阵模板 for i=apha:apha:x1 for j=apha:apha:y1 A(i,j)=X(i/apha,j/apha); end end%将原图像的值对应赋予新图像的位置 a=floor(y1/6);a1=mod(y1,6); b=round(x1/6);b1=mod(x1,6); for i=apha:apha:x1%先对行进行变化，使用newton方法,采取每八个点为一段，若段内变化率大于30%，直接赋值相邻两点中值 k=0; for j=1:6:y1 k=k+1; if (k==a+1)%检验图像的每段中是否为8个点，每推进一段，重合两点 if(a1~=0) for l=j:j+a1-1 A(i,l)=A(i,j); end end elseif(k==1) if(bianhualv(A(i,apha),A(i,2*apha))&&bianhualv(A(i,2*apha),A(i,3*apha))&&bianhualv(A(i,3*apha),A(i,4*apha)))%如果变化率小于30%，进行插值 s=newtoncz([0.1 0.2 0.3 0.4],[A(i,apha),A(i,2*apha),A(i,3*apha),A(i,4*apha)],[0 0.15 0.25 0.35]); A(i,1)=s(1,1);A(i,3)=s(2,1);A(i,5)=s(3,1);A(i,7)=s(4,1); else A(i,1)=A(i,2);A(i,3)=(A(i,2)+A(i,4))/2;A(i,5)=(A(i,4)+A(i,6))/2;A(i,7)=(A(i,6)+A(i,8))/2; end elseif(k==a) break; else if(bianhualv(A(i,j),A(i,j+1))&&bianhualv(A(i,j+1),A(i,j+3))&&bianhualv(A(i,j+3),A(i,j+5))&&bianhualv(A(i,j+5),A(i,j+7)))%如果变化率小于30%，进行插值 s=newtoncz([0.15 0.2 0.3 0.4 0.5],[A(i,j),A(i,j+1),A(i,j+3),A(i,j+5),A(i,j+7)],[0.25 0.35 0.45]); A(i,j+2)=s(1,1);A(i,j+4)=s(2,1);A(i,j+6)=s(3,1); else A(i,j+2)=(A(i,j+1)+A(i,j+3))/2;A(i,j+4)=(A(i,j+3)+A(i,j+5))/2;A(i,j+6)=(A(i,j+5)+A(i,j+7))/2; end end end end for j=apha:apha:y1%对列进行变化，使用newton方法,采取每八个点为一段，若段内变化率大于30%，直接赋值相邻两点中值 k=0; for i=1:6:x1 k=k+1; if (k==b+1)%检验图像的每段中是否为8个点，每推进一段，重合两点 if(b1~=0) for l=j:j+b1-1 A(i,j)=A(i+l,j); end end elseif(k==1) if(bianhualv(A(apha,i),A(2*apha,i))&&bianhualv(A(2*apha,i),A(3*apha,i))&&bianhualv(A(3*apha,i),A(4*apha,i)))%如果变化率小于30%，进行插值 s=newtoncz([0.1 0.2 0.3 0.4],[A(apha,i),A(2*apha,i),A(3*apha,i),A(4*apha,i)],[0 0.15 0.25 0.35]); A(1,j)=s(1,1);A(3,j)=s(2,1);A(5,j)=s(3,1);A(7,j)=s(4,1); else A(1,j)=A(2,j);A(3,j)=(A(2,j)+A(4,j))/2;A(5,j)=(A(4,j)+A(6,j))/2;A(7,j)=(A(6,j)+A(8,j))/2; end elseif(k==a) break; else if(bianhualv(A(i,j),A(i+1,j))&&bianhualv(A(i+1,j),A(i+3,j))&&bianhualv(A(i+3,j),A(i+5,j))&&bianhualv(A(i+5,j),A(i+7,j)))%如果变化率小于30%，进行插值 s=newtoncz([0.15 0.2 0.3 0.4 0.5],[A(i,j),A(i+1,j),A(i+3,j),A(i+5,j),A(i+7,j)],[0.25,0.35,0.45]); A(i+2,j)=s(1,1);A(i+4,j)=s(2,1);A(i+6,j)=s(3,1); else A(i+2,j)=(A(i+1,j)+A(i+3,j))/2;A(j+4,i)=(A(i+3,j)+A(i+5,j))/2;A(i+6,j)=(A(i+5,j)+A(i+7,j))/2; end end end end A=im2uint8(A); end function d=bianhualv(a,b)%计算相邻像素变化率 c=abs(a-b); if(c/a<=0.3) d=1; else d=0; end end
处理示例：

004.jpg 005.jpg

效果太差了，不知道为什么越到下面越模糊，大家看看这程序，找找问题。

Netwon-Thiele算法

function B2=netwon2(Z,T)%Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 M=Z'; [n,m]=size(M);%取得插值矩阵的大小 %首先对x方向newton递推 M1=zeros(n,m); y1=1:m; for i1=1:n L=M(i1,:); S=netwon(y1,L)%对每一行进行牛顿递推 M1(i1,:)=S; end %然后对y方向进行thiele递推 M2=zeros(n,m); x1=1:n; for j1=1:m L1=M1(:,j1); L1=L1'; S1=nichashang(x1,L1); M2(:,j1)=S1; end %利用M2表中的值来计算y待插值点值 [n1,m1]=size(T); A1=T(:,2); A1=A1';%获得y的坐标 A2=zeros(m,n1); for j1=1:m L1=M2(:,j1); L1=L1'; L2=thieledis(A1,L1,x1); A2(j1,:)=L2; end %计算（x,y）待插值点; B1=T(:,1); B1=B1'; B2=zeros(n1,n1); for i1=n1 l1=A2(:,n1); l1=l1'; Z=netwondis(B1,y1,l1); B2(n1,:)=Z; end end function S=netwon(x,y) x=x';y=y';%转致 hangx=size(x,1);%获取x的行数; H=ones(hangx,hangx);%生成行数等于x的矩阵，用来存取系数 for j=1:hangx%这个循环是算法核心，即那个倒三角的表 if(j==1) H(:,1)=y; else for i=j:hangx H(i,j)=(H(i,j-1)-H(i-1,j-1))/(x(i,1)-x(i-j+1,1)); end end end S=zeros(1,hangx); for i=1:hangx S(1,i)=H(i,i); end end function c=nichashang(a,b)%计算逆插商的值 [x,y]=size(a); d=zeros(y,y); c=ones(1,y); for j=1:y if(j==1) for i=1:y d(i,1)=b(1,i); end else for i=j:y d(i,j)=(a(1,i)-a(1,j-1))/(d(i,j-1)-d(j-1,j-1)); end end end for k=1:y c(1,k)=d(k,k); end c=c'; end function c=thieledis(X,Y,Z)%x为待求点值，y为thiele多项式系数，z为原插值点值 [x,y]=size(Y); [x1,y1]=size(X); c=ones(1,y1); for i=1:y1 A=ones(y+1,1); A(1,1)=1;A(2,1)=Y(1,1); B=ones(y+1,1); B(1,1)=0;B(2,1)=1; for k=3:y+1 A(k,1)=A(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*A(k-2,1); B(k,1)=B(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*B(k-2,1); end c(1,i)=A(y+1,1)/B(y+1,1); end end function Z=netwondis(z,x,y)%输出插商的值，z为待插值点，x为原插值点，y为插商值 hangz=size(z,1); Z=ones(hangz,1); x=x'; hangx=size(x,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+y(1,k)*S; S=S.*(z-x(k,1)); end end Z=Z'; end
基于自适应的超分辨率图像重建
function B1=tuxiangcj(X,apha,p)%X为待处理图片，apha为放大倍数，先假设为2,区分平坦块系数 %预处理，先将X规范化，能进行4*4的块运算 X=im2double(X); [x,y]=size(X); A=X; if(mod(x,2)~=0) A=[A;A(x,:)]; end if(mod(y,2)~=0) A=[A,A(:,y)]; end %对图像进行放大 [x1,y1]=size(A); x2=x1*apha;y2=y1*apha; B=zeros(x2,y2);%生成插值后矩阵模板 for i=apha:apha:x2 for j=apha:apha:y2 B(i,j)=A(i/apha,j/apha); end end%将原图像的值对应赋予新图像的位置 B1=B; %对图像进行纹理块或平坦块判断; for i=1:2:x1-2 for j=1:2:y1-2 L=[A(i,j),A(i,j+1),A(i,j+2),A(i,j+3); A(i+1,j),A(i+1,j+1),A(i+1,j+2),A(i+1,j+3); A(i+2,j),A(i+2,j+1),A(i+2,j+2),A(i+2,j+3); A(i+3,j),A(i+3,j+1),A(i+3,j+2),A(i+3,j+3)]; x=panduan(L,p); if(x==1)%对平坦块使用牛顿牛顿法 L1=netnetwon([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 else L1=netwon2([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 end if(i==1&&j==1) B1(2*i-1:2*i+6,2*j-1:2*j+6)=L1; end if(i==1&&j~=1) L2=(L1(1:8,1:4)+B1(2*i-1:2*i+6,2*j-1:2*j+2))./2; L1=[L2,L1(1:8,5:8)]; B1(2*i-1:2*i+6,2*j-1:2*j+6)=L1; end if(i~=1&&j==1) L2=(L1(1:4,1:8)+B1(2*i-1:2*i+2,2*j-1:2*j+6))./2; L1=[L2;L1(5:8,1:8)]; B1(2*i-1:2*i+6,2*j-1:2*j+6)=L1; end if(i~=1&&j~=1) L2=(L1(1:4,1:8)+B1(2*i-1:2*i+2,2*j-1:2*j+6))./2; L3=(L1(1:8,1:4)+B1(2*i-1:2*i+6,2*j-1:2*j+2))./2; L4=(L2(1:4,1:4)+L3(1:4,1:4))./2; L1=[L4,L2(1:4,5:8);L3(5:8,1:4),L1(5:8,5:8)]; B1(2*i-1:2*i+6,2*j-1:2*j+6)=L1; end end end B1=im2uint8(B1); imshow(B1); end function X=panduan(A,p) B=sum(sum(A))/16; s=0; for i1=1:4 for j1=1:4 s=s+(B-A(i1,j1))^2; end end if(s<p) X=1; else X=0; end end
function B2=netwon2(X,Y,Z,T)%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 M=Z'; [n,m]=size(M);%取得插值矩阵的大小 %首先对x方向newton递推 M1=zeros(n,m); for i1=1:n L=M(i1,:); S=netwon(X,L);%对每一行进行牛顿递推 M1(i1,:)=S; end %然后对y方向进行thiele递推 M2=zeros(n,m); k2=0; for j1=1:m L1=M1(:,j1); L1=L1'; [S1,k1]=nichashang(Y,L1); M2(:,j1)=S1; if(k1==1) k2=1; break; end end if(k2==1) M2=zeros(n,m); for j1=1:m L1=M1(:,j1); L1=L1'; S1=netwon(Y,L1); M2(:,j1)=S1; end end %利用M2表中的值来计算y待插值点值 [n1,m1]=size(T); A1=T(:,2); A1=A1';%获得y的坐标 A2=zeros(m,n1); if(k2~=1) for j1=1:m L1=M2(:,j1); L1=L1'; L2=thieledis(A1,L1,Y); A2(j1,:)=L2; end else for j1=1:m L1=M2(:,j1); L1=L1'; L2=netwondis(A1,Y,L1); A2(j1,:)=L2; end end %计算（x,y）待插值点; B1=T(:,1); B1=B1'; B2=zeros(n1,n1); for i1=1:n1 l1=A2(:,i1); l1=l1'; Z=netwondis(B1,X,l1); Z=Z'; B2(i1,:)=Z; end end function S=netwon(x,y) x=x';y=y';%转致 hangx=size(x,1);%获取x的行数; H=ones(hangx,hangx);%生成行数等于x的矩阵，用来存取系数 for j=1:hangx%这个循环是算法核心，即那个倒三角的表 if(j==1) H(:,1)=y; else for i=j:hangx H(i,j)=(H(i,j-1)-H(i-1,j-1))/(x(i,1)-x(i-j+1,1)); end end end S=zeros(1,hangx); for i=1:hangx S(1,i)=H(i,i); end end function [c,k1]=nichashang(a,b)%计算逆插商的值 [x,y]=size(a); d=zeros(y,y); c=ones(1,y); k1=0; for j=1:y if(j==1) for i=1:y d(i,1)=b(1,i); end else for i=j:y if(abs(d(i,j-1)-d(j-1,j-1))<0.0001) k1=1; break; end d(i,j)=(a(1,i)-a(1,j-1))/(d(i,j-1)-d(j-1,j-1)); end if(k1==1) break; end end if(k1==1) break; end end for k=1:y c(1,k)=d(k,k); end c=c'; end function c=thieledis(X,Y,Z)%x为待求点值，y为thiele多项式系数，z为原插值点值 y=size(Y,2); y1=size(X,2); c=ones(1,y1); for i=1:y1 A=ones(y+1,1); A(1,1)=1;A(2,1)=Y(1,1); B=ones(y+1,1); B(1,1)=0;B(2,1)=1; for k=3:y+1 A(k,1)=A(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*A(k-2,1); B(k,1)=B(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*B(k-2,1); end c(1,i)=A(y+1,1)/B(y+1,1); end end function Z=netwondis(z,x,y)%输出插商的值，z为待插值点，x为原插值点，y为插商值 hangz=size(z,1); Z=ones(hangz,1); x=x'; hangx=size(x,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+y(1,k)*S; S=S.*(z-x(k,1)); end end Z=Z'; end

function B2=netnetwon(X,Y,Z,T) M=Z'; [n,m]=size(M);%取得插值矩阵的大小 %首先对x方向newton递推 M1=zeros(n,m); for i1=1:n L=M(i1,:); S=netwon(X,L);%对每一行进行牛顿递推 M1(i1,:)=S; end %然后对y方向进行netwon递推 M2=zeros(n,m); for j1=1:m L1=M1(:,j1); L1=L1'; S1=netwon(Y,L1); M2(:,j1)=S1; end %利用M2表中的值来计算y待插值点值 [n1,m1]=size(T); A1=T(:,2); A1=A1';%获得y的坐标 A2=zeros(m,n1); for j1=1:m L1=M2(:,j1); L1=L1'; L2=netwondis(A1,Y,L1); A2(j1,:)=L2; end %计算（x,y）待插值点; B1=T(:,1); B1=B1'; B2=zeros(n1,n1); for i1=1:n1 l1=A2(:,i1); l1=l1'; Z=netwondis(B1,X,l1); B2(i1,:)=Z; end end function S=netwon(x,y) x=x';y=y';%转致 hangx=size(x,1);%获取x的行数; H=ones(hangx,hangx);%生成行数等于x的矩阵，用来存取系数 for j=1:hangx%这个循环是算法核心，即那个倒三角的表 if(j==1) H(:,1)=y; else for i=j:hangx H(i,j)=(H(i,j-1)-H(i-1,j-1))/(x(i,1)-x(i-j+1,1)); end end end S=zeros(1,hangx); for i=1:hangx S(1,i)=H(i,i); end end function Z=netwondis(z,x,y)%输出插商的值，z为待插值点，x为原插值点，y为插商值 hangz=size(z,1); Z=ones(hangz,1); x=x'; hangx=size(x,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+y(1,k)*S; S=S.*(z-x(k,1)); end end end
由于这程序比较简陋，很多地方没有重构和优化，故运行有点慢，效果一般。下一步优化算法，改成对彩色图像处理。
示例：
修复前：

test2.jpg

修复后：

test4.jpg test3.jpg

由于算法细节没有弄好，细节上与matlab有些差别，大家一起来优化算法。
2016.08.05
再更新下彩色图像版，只有主程序变动,但是效果很差
function H=tuxiangcj(X,apha,p)%X为待处理图片，apha为放大倍数，先假设为2,区分平坦块系数 %预处理，先将X规范化，能进行4*4的块运算 R=X(:,:,1); G=X(:,:,2); B=X(:,:,3); R=im2double(R); G=im2double(G); B=im2double(B); [x,y]=size(R); R1=R; if(mod(x,2)~=0) R1=[R1;R1(x,:)]; end if(mod(y,2)~=0) R1=[R1,R1(:,y)]; end G1=G; if(mod(x,2)~=0) G1=[G1;G1(x,:)]; end if(mod(y,2)~=0) G1=[G1,G1(:,y)]; end B1=B; if(mod(x,2)~=0) B1=[B1;B1(x,:)]; end if(mod(y,2)~=0) B1=[B1,B1(:,y)]; end %对图像进行放大 [x1,y1]=size(R1); x2=x1*apha;y2=y1*apha; R2=zeros(x2,y2);%生成插值后矩阵模板 for i=apha:apha:x2 for j=apha:apha:y2 R2(i,j)=R1(i/apha,j/apha); end end%将原图像的值对应赋予新图像的位置 B2=zeros(x2,y2);%生成插值后矩阵模板 for i=apha:apha:x2 for j=apha:apha:y2 B2(i,j)=B1(i/apha,j/apha); end end%将原图像的值对应赋予新图像的位置 G2=zeros(x2,y2);%生成插值后矩阵模板 for i=apha:apha:x2 for j=apha:apha:y2 G2(i,j)=G1(i/apha,j/apha); end end%将原图像的值对应赋予新图像的位置 %对图像进行纹理块或平坦块判断; for i=1:2:x1-2 for j=1:2:y1-2 L1=[R1(i,j),R1(i,j+1),R1(i,j+2),R1(i,j+3); R1(i+1,j),R1(i+1,j+1),R1(i+1,j+2),R1(i+1,j+3); R1(i+2,j),R1(i+2,j+1),R1(i+2,j+2),R1(i+2,j+3); R1(i+3,j),R1(i+3,j+1),R1(i+3,j+2),R1(i+3,j+3)]; L2=[G1(i,j),G1(i,j+1),G1(i,j+2),G1(i,j+3); G1(i+1,j),G1(i+1,j+1),G1(i+1,j+2),G1(i+1,j+3); G1(i+2,j),G1(i+2,j+1),G1(i+2,j+2),G1(i+2,j+3); G1(i+3,j),G1(i+3,j+1),G1(i+3,j+2),G1(i+3,j+3)]; L3=[B1(i,j),B1(i,j+1),B1(i,j+2),B1(i,j+3); B1(i+1,j),B1(i+1,j+1),B1(i+1,j+2),B1(i+1,j+3); B1(i+2,j),B1(i+2,j+1),B1(i+2,j+2),B1(i+2,j+3); B1(i+3,j),B1(i+3,j+1),B1(i+3,j+2),B1(i+3,j+3)]; x=panduan(L1,L2,L3,p); if(x==1)%对平坦块使用牛顿牛顿法 RL=netnetwon([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L1,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 GL=netnetwon([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L2,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 BL=netnetwon([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L3,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 else RL=netwon2([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L1,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 GL=netwon2([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L2,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 BL=netwon2([0.1 0.2 0.3 0.4],[0.1 0.2 0.3 0.4],L3,[0.05 0.05;0.1 0.1;0.15 0.15;0.2 0.2;0.25 0.25;0.3 0.3;0.35 0.35;0.4 0.4]);%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 end if(i==1&&j==1) R2(2*i-1:2*i+6,2*j-1:2*j+6)=RL; G2(2*i-1:2*i+6,2*j-1:2*j+6)=GL; B2(2*i-1:2*i+6,2*j-1:2*j+6)=BL; end if(i==1&&j~=1) RL2=(RL(1:8,1:4)+R2(2*i-1:2*i+6,2*j-1:2*j+2))./2; RL1=[RL2,RL(1:8,5:8)]; R2(2*i-1:2*i+6,2*j-1:2*j+6)=RL1; GL2=(GL(1:8,1:4)+G2(2*i-1:2*i+6,2*j-1:2*j+2))./2; GL1=[GL2,GL(1:8,5:8)]; G2(2*i-1:2*i+6,2*j-1:2*j+6)=GL1; BL2=(BL(1:8,1:4)+B2(2*i-1:2*i+6,2*j-1:2*j+2))./2; BL1=[BL2,BL(1:8,5:8)]; B2(2*i-1:2*i+6,2*j-1:2*j+6)=BL1; end if(i~=1&&j==1) RL2=(RL(1:4,1:8)+R2(2*i-1:2*i+2,2*j-1:2*j+6))./2; RL1=[RL2;RL(5:8,1:8)]; R2(2*i-1:2*i+6,2*j-1:2*j+6)=RL1; GL2=(GL(1:4,1:8)+G2(2*i-1:2*i+2,2*j-1:2*j+6))./2; GL1=[GL2;GL(5:8,1:8)]; G2(2*i-1:2*i+6,2*j-1:2*j+6)=GL1; BL2=(BL(1:4,1:8)+B2(2*i-1:2*i+2,2*j-1:2*j+6))./2; BL1=[BL2;BL(5:8,1:8)]; B2(2*i-1:2*i+6,2*j-1:2*j+6)=BL1; end if(i~=1&&j~=1) RL2=(RL(1:4,1:8)+R2(2*i-1:2*i+2,2*j-1:2*j+6))./2; RL3=(RL(1:8,1:4)+R2(2*i-1:2*i+6,2*j-1:2*j+2))./2; RL4=(RL2(1:4,1:4)+RL3(1:4,1:4))./2; RL1=[RL4,RL2(1:4,5:8);RL3(5:8,1:4),RL1(5:8,5:8)]; R2(2*i-1:2*i+6,2*j-1:2*j+6)=RL1; GL2=(GL(1:4,1:8)+G2(2*i-1:2*i+2,2*j-1:2*j+6))./2; GL3=(GL(1:8,1:4)+G2(2*i-1:2*i+6,2*j-1:2*j+2))./2; GL4=(GL2(1:4,1:4)+GL3(1:4,1:4))./2; GL1=[GL4,GL2(1:4,5:8);GL3(5:8,1:4),GL1(5:8,5:8)]; G2(2*i-1:2*i+6,2*j-1:2*j+6)=GL1; BL2=(BL(1:4,1:8)+B2(2*i-1:2*i+2,2*j-1:2*j+6))./2; BL3=(BL(1:8,1:4)+B2(2*i-1:2*i+6,2*j-1:2*j+2))./2; BL4=(BL2(1:4,1:4)+BL3(1:4,1:4))./2; BL1=[BL4,BL2(1:4,5:8);BL3(5:8,1:4),BL1(5:8,5:8)]; B2(2*i-1:2*i+6,2*j-1:2*j+6)=BL1; end end end R2=im2uint8(R2); G2=im2uint8(G2); B2=im2uint8(B2); H(:,:,1)=R2; H(:,:,2)=G2; H(:,:,3)=B2; imshow(H); end function X=panduan(R,G,B,p) S1=sum(sum(R))/16; S2=sum(sum(G))/16; S3=sum(sum(B))/16; s=0; for i1=1:4 for j1=1:4 s=s+(S1-R(i1,j1))^2+(S2-G(i1,j1))^2+(S3-B(i1,j1))^2; end end if(s<p) X=1; else X=0; end end

N-T算法优化----考虑奇异点的情况
function Z1=NT(X,Y,Z,X1,Y1)%X,Y,Z为插值点坐标，T为待插值点坐标横纵坐标，P为T的插值结果 M=Z'; [n,m]=size(M);%取得插值矩阵的大小 %首先对x方向newton递推 x=X(1,:); M1=zeros(n,m); for i1=1:n L=M(i1,:); S=netwon(x,L);%对每一行进行牛顿递推 M1(i1,:)=S; end %然后对y方向进行thiele递推 y=Y(:,1)'; M2=zeros(n,m); P=zeros(1,m); for j1=1:m L1=M1(:,j1); L1=L1'; [S1,p]=thiele1(y,L1); P(1,j1)=p; M2(:,j1)=S1; end %利用M2表中的值来计算y待插值点值 y1=Y1(:,1)'; n1=size(y1,2); %获得y的坐标 A2=zeros(m,n1); for j1=1:m L1=M2(:,j1); L1=L1'; L2=thieledis1(y1,L1,y,P(1,j1)); A2(j1,:)=L2; end %计算（x,y）待插值点; x1=X1(1,:); m1=size(x1,2); Z1=zeros(n1,m1); %获取x坐标 for i1=1:n1 l1=A2(:,i1); l1=l1'; Z=netwondis(x1,l1,x); Z=Z'; Z1(i1,:)=Z; end end function S=netwon(x,y) x=x';y=y';%转致 hangx=size(x,1);%获取x的行数; H=ones(hangx,hangx);%生成行数等于x的矩阵，用来存取系数 for j=1:hangx%这个循环是算法核心，即那个倒三角的表 if(j==1) H(:,1)=y; else for i=j:hangx H(i,j)=(H(i,j-1)-H(i-1,j-1))/(x(i,1)-x(i-j+1,1)); end end end S=zeros(1,hangx); for i=1:hangx S(1,i)=H(i,i); end end
function C=thieledis1(X,Y,Z,p)%X为待插值点，Y为逆差商，Z为原插值节点横坐标，p为奇异点坐标 XH=size(X,2); YH=size(Y,2); if(p~=0) c=netwondis(X,Y(1,p+1:YH),Z(1,p+1:YH)); c=c'; c1=(X-Z(1,p)).*c+Y(1,p); for i=1:XH if(c1(1,i)>0.0001) C=thieledis(X,[Y(1,1:p-1),c1(1,i)],Z(1,1:p)); else C=thieledis(X,Y(1,1:p-1),Z(1,1:(p-1))); end end else C=thieledis(X,Y,Z); end end function Z=netwondis(z,y,x)%输出插商的值，z为待插值点，x为原插值点，y为插商值 hangz=size(z,1); Z=ones(hangz,1); x=x'; hangx=size(x,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+y(1,k)*S; S=S.*(z-x(k,1)); end end Z=Z'; end function c=thieledis(X,Y,Z)%x为待求点值，y为thiele多项式系数，z为原插值点值 y=size(Y,2); y1=size(X,2); c=ones(1,y1); for i=1:y1 A=ones(y+1,1); A(1,1)=1;A(2,1)=Y(1,1); B=ones(y+1,1); B(1,1)=0;B(2,1)=1; for k=3:y+1 A(k,1)=A(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*A(k-2,1); B(k,1)=B(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*B(k-2,1); end c(1,i)=A(y+1,1)/B(y+1,1); end end

function c=thieledis(X,Y,Z)%x为待求点值，y为thiele多项式系数，z为原插值点值 y=size(Y,2); y1=size(X,2); c=ones(1,y1); for i=1:y1 A=ones(y+1,1); A(1,1)=1;A(2,1)=Y(1,1); B=ones(y+1,1); B(1,1)=0;B(2,1)=1; for k=3:y+1 A(k,1)=A(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*A(k-2,1); B(k,1)=B(k-1,1)*Y(1,k-1)+(X(1,i)-Z(k-2))*B(k-2,1); end c(1,i)=A(y+1,1)/B(y+1,1); end end

function Z=netwondis(z,y,x)%输出插商的值，z为待插值点，x为原插值点，y为插商值 hangz=size(z,1); Z=ones(hangz,1); x=x'; hangx=size(x,1); for k=1:hangx if(k==1) Z(:,1)=y(1,1); S=z-x(1,1); else Z=Z+y(1,k)*S; S=S.*(z-x(k,1)); end end Z=Z'; end
function C=thieledis1(X,Y,Z,p)%X为待插值点，Y为逆差商，Z为原插值节点横坐标，p为奇异点坐标 XH=size(X,2); YH=size(Y,2); if(p~=0) c=netwondis(X,Y(1,p+1:YH),Z(1,p+1:YH)); c=c'; c1=(X-Z(1,p)).*c+Y(1,p); for i=1:XH if(c1(1,i)>0.0001) C=thieledis(X,[Y(1,1:p-1),c1(1,i)],Z(1,1:p)); else C=thieledis(X,Y(1,1:p-1),Z(1,1:(p-1))); end end else C=thieledis(X,Y,Z); end end
function [c,count]=thiele1(a,b)%a,b为待插值点的坐标 y=size(a,2); d=zeros(y,y); c=ones(1,y); count=0; p=0; for j=1:y if(j==1) for i=1:y d(i,1)=b(1,i); end else %若逆差商为0，转换为差商 for i=j:y if(abs(d(i,j-1)-d(j-1,j-1))>0.0000001&&p==0) d(i,j)=(a(1,i)-a(1,j-1))/(d(i,j-1)-d(j-1,j-1)); else d(i,j)=(d(i,j-1)-d(i-1,j-1))/(a(1,i)-a(1,i-j+1)); if(p==0) count=j; p=1; end end end end end for k=1:y c(1,k)=d(k,k); end end