我剛寫下了返回係數的函數。假設16個值來自(x = -dx,0,dx,2dx)和(y = -dx,0,dx,2dx)。我還寫了函數'Deriv',它涉及我的數值函數關於x,y,xy的導數:
function [out]=Deriv(in,str,dx)
switch str
case 'x'
out=(circshift(in,[0 -1])-circshift(in,[0 1]))./(2*dx);
case 'y'
out=(circshift(in,[-1 0])-circshift(in,[1 0]))./(2*dx);
case 'xy'
out=(circshift(in,[-1 -1])+circshift(in,[1 1])-circshift(in,[-1 1])-circshift(in,[1 -1]))./(4*dx*dx);
end
function out=bl(in1,in2)
out=[in1^3, in1^2, in1, 1,in2*[in1^3, in1^2, in1, 1], in2^2*[in1^3, in1^2, in1, 1]...
in2^3*[in1^3, in1^2, in1, 1]];
end
function out=blx(in1,in2)
out=[3*in1^2, 2*in1, 1, 0,in2*[3*in1^2, 2*in1, 1, 0], in2^2*[3*in1^2, 2*in1, 1, 0]...
in2^3*[3*in1^2, 2*in1, 1, 0]];
end
function out=bly(in1,in2)
out=[0,0,0,0,[in1^3, in1^2, in1, 1], 2*in2*[in1^3, in1^2, in1, 1]...
3*in2^2*[in1^3, in1^2, in1, 1]];
end
function out=MyBicubic(in,dx)
A=zeros(16,16); C=zeros(16,1);
A(1,:)=bl(0,0); C(1)=in(2,2);
A(2,:)=bl(dx,0); C(2)=in(2,3);
A(3,:)=bl(0,dx); C(3)=in(3,2);
A(4,:)=bl(dx,dx); C(4)=in(3,3);
A(5,:)=blx(0,0); temp=Deriv(in,'x',dx); C(5)=temp(2,2);
A(6,:)=blx(dx,0); C(6)=temp(2,3);
A(7,:)=blx(0,dx); C(7)=temp(3,2);
A(8,:)=blx(dx,dx); C(8)=temp(3,3);
A(9,:)=bly(0,0); temp=Deriv(in,'y',dx); C(9)=temp(2,2);
A(10,:)=bly(dx,0); C(10)=temp(2,3);
A(11,:)=bly(0,dx); C(11)=temp(3,2);
A(12,:)=bly(dx,dx); C(12)=temp(3,3);
A(13,:)=blxy(0,0); temp=Deriv(in,'xy',dx); C(13)=temp(2,2);
A(14,:)=blxy(dx,0); C(14)=temp(2,3);
A(15,:)=blxy(0,dx); C(15)=temp(3,2);
A(16,:)=blxy(dx,dx); C(16)=temp(3,3);
out=A\C;
end