bicubic interpolation of 3D surface

Question

I am writing an algorithm in matlab to for bicubic interpolation of a surface Psi(x,y). I have a bug in the code and cannot seem to track it down. I am trying a test case with Psi=X^2-0.25 so that its easier to track down the bug. It seems as if my interpolation has an offset. My comments are included in my code. Any help would be appreciated.

Plot of Psi=X^2 in blue and interpolation in red

Countour lines of Psi are plotted and the red dot is the point I am computing the interpolation about. The thick red line is the interpolation which is offset quite a bit from the red dot.

function main()

    epsilon=0.000001;
    xMin=-1+epsilon;
    xMax= 1+epsilon;
    yMin=-1+epsilon;
    yMax= 1+epsilon;

    dx=0.1;  Nx=ceil((xMax-xMin)/dx)+1;
    dy=0.1;  Ny=ceil((yMax-yMin)/dy)+1;

    x=xMin:dx:xMax;  x=x(1:Nx);
    y=yMin:dy:yMax;  y=y(1:Ny);

    [XPolInX,XPolInY]=GetGhostMatricies(Nx,Ny); %Linear extrapolation matrix
    [D0x,D0y]=GetDiffMatricies(Nx,Ny,dx,dy); %derivative matricies: D0x is central differencing in x

    [X,Y]=meshgrid(x,y);
    Psi=X.^2-0.25;    %Note that my algorithm is being written for a Psi that may not have a analytic representation. This Psi is only a test case. 

    psi=zeros(Nx+2,Ny+2);   %linearly extrapolate psi (for solving differential equation not shown here)
    psi(2:(Nx+1),2:(Ny+1))=Psi';
    psi=(XPolInY*(XPolInX*psi)')';

    %compute derivatives of psi
    psi_x =D0x*psi;             psi_x =(XPolInY*(XPolInX*psi_x)')';
    psi_y =(D0y*psi')';         psi_y =(XPolInY*(XPolInX*psi_y)')';
    psi_xy=D0x*psi_y;           psi_xy=(XPolInY*(XPolInX*psi_xy)')'; 
    % i have verified that my derivatives are computed correctly    

    biCubInv=GetBiCubicInverse(dx,dy);  

    i=5;  %lets compute the bicubic interpolation at this x(i), y(j)
    j=1;

    psiVoxel=[psi(   i,j),psi(   i+1,j),psi(   i,j+1),psi(   i+1,j+1),...
              psi_x( i,j),psi_x( i+1,j),psi_x( i,j+1),psi_x( i+1,j+1),...
              psi_y( i,j),psi_y( i+1,j),psi_y( i,j+1),psi_y( i+1,j+1),...
              psi_xy(i,j),psi_xy(i+1,j),psi_xy(i,j+1),psi_xy(i+1,j+1)]';

    a=biCubInv*psiVoxel; %a=[a00 a01 ... a33]; polynomial coefficients; 1st index is power of (x-xi), 2nd index is power of (y-yj)

    xi=x(5); yj=y(1);
    clear x y
    x=(xi-.2):.01:(xi+.2);   %this is a local region about the point we are interpolating
    y=(yj-.2):.01:(yj+.2);
    [dX,dY]=meshgrid(x,y);
    Psi=dX.^2-0.25;   

    figure(2)            %just plotting the 0 level contour of Psi here 
    plot(xi,yj,'.r','MarkerSize',20)
    hold on
    contour(x,y,Psi,[0 0],'r','LineWidth',2)
    set(gca,'FontSize',14)
    axis([x(1) x(end) y(1) y(end)])
    grid on
    set(gca,'xtick',(xi-.2):.1:(xi+.2));
    set(gca,'ytick',(yj-.2):.1:(yj+.2));
    xlabel('x')
    ylabel('y')

    [dX dY]=meshgrid(x-xi,y-yj);
    %P is my interpolating polynomial
    P =   a(1)       + a(5)       *dY + a(9)        *dY.^2 + a(13)       *dY.^3 ...
        + a(2)*dX    + a(6)*dX   .*dY + a(10)*dX   .*dY.^2 + a(14)*dX   .*dY.^3 ...
        + a(3)*dX.^2 + a(7)*dX.^2.*dY + a(11)*dX.^2.*dY.^2 + a(15)*dX.^2.*dY.^3 ...
        + a(4)*dX.^3 + a(8)*dX.^3.*dY + a(12)*dX.^3.*dY.^2 + a(16)*dX.^3.*dY.^3 ;
    [c h]=contour(x,y,P)  
    clabel(c,h)

    figure(3)
    plot(x,x.^2-.25)  %this is the exact function
    hold on
    plot(x,P(1,:),'-r*')

    %See there is some offset here 


end
%-------------------------------------------------------------------------
function [XPolInX,XPolInY]=GetGhostMatricies(Nx,Ny)

    XPolInX=diag(ones(1,Nx+2),0);
    XPolInY=diag(ones(1,Ny+2),0);

    XPolInX(1,1)      =0; XPolInX(1,2)      =2; XPolInX(1,3)    =-1;
    XPolInY(1,1)      =0; XPolInY(1,2)      =2; XPolInY(1,3)    =-1;
    XPolInX(Nx+2,Nx+2)=0; XPolInX(Nx+2,Nx+1)=2; XPolInX(Nx+2,Nx)=-1;
    XPolInY(Ny+2,Ny+2)=0; XPolInY(Ny+2,Ny+1)=2; XPolInY(Ny+2,Ny)=-1;

    fprintf('Done GetGhostMatricies\n')
end
%-------------------------------------------------------------------------
function [D0x,D0y]=GetDiffMatricies(Nx,Ny,dx,dy)

    D0x=diag(ones(1,Nx-1),1)-diag(ones(1,Nx-1),-1); 
    D0y=diag(ones(1,Ny-1),1)-diag(ones(1,Ny-1),-1);
    D0x(1,1)=-3; D0x(1,2)=4; D0x(1,3)=-1;
    D0y(1,1)=-3; D0y(1,2)=4; D0y(1,3)=-1;
    D0x(Nx,Nx)=3; D0x(Nx,Nx-1)=-4; D0x(Nx,Nx-2)=1;
    D0y(Ny,Ny)=3; D0y(Ny,Ny-1)=-4; D0y(Ny,Ny-2)=1;

    %pad with ghost cells which are simply zeros
    tmp=D0x;     D0x=zeros(Nx+2,Nx+2);   D0x(2:(Nx+1),2:(Nx+1))=tmp; tmp=0;
    tmp=D0y;     D0y=zeros(Ny+2,Ny+2);   D0y(2:(Ny+1),2:(Ny+1))=tmp; tmp=0;

    %scale appropriatley by dx & dy
    D0x=D0x/(2*dx);
    D0y=D0y/(2*dy);

end
%-------------------------------------------------------------------------
function biCubInv=GetBiCubicInverse(dx,dy)
    %p(x,y)=a00+a01(x-xi)+a02(x-xi)^2+...+a33(x-xi)^3(y-yj)^3
    %biCubic*a=[psi(i,j) psi(i+1,j) psi(i,j+1) psi(i+1,j+1) psi_x(i,j) ... psi_y(i,j) ... psi_xy(i,j) ... psi_xy(i+1,j+1)] 
    %here, psi_x is the x derivative of psi
    %I verified that this matrix is correct by setting dx=dy=1 and comparing to the inverse here http://en.wikipedia.org/wiki/Bicubic_interpolation

    biCubic=[ 
       %00      10      20      30      01      11      21      31        02      12      22        32          03        13        23          33
        1       0       0       0       0       0       0       0         0       0       0         0           0         0         0           0;
        1       dx      dx^2    dx^3    0       0       0       0         0       0       0         0           0         0         0           0;
        1       0       0       0       dy      0       0       0         dy^2    0       0         0           dy^3      0         0           0;
        1       dx      dx^2    dx^3    dy      dx*dy   dx^2*dy dx^3*dy   dy^2    dx*dy^2 dx^2*dy^2 dx^3*dy^2   dy^3      dx*dy^3   dx^2*dy^3   dx^3*dy^3;

        0       1       0       0       0       0       0       0         0       0       0         0           0         0         0           0;
        0       1       2*dx    3*dx^2  0       0       0       0         0       0       0         0           0         0         0           0;
        0       1       0       0       0       dy      0       0         0       dy^2    0         0           0         dy^3      0           0;
        0       1       2*dx    3*dx^2  0       dy      2*dx*dy 3*dx^2*dy 0       dy^2    2*dx*dy^2 3*dx^2*dy^2 0         dy^3      2*dx*dy^3   3*dx^2*dy^3;

        0       0       0       0       1       0       0       0         0       0       0         0           0         0         0           0;
        0       0       0       0       1       dx      dx^2    dx^3      0       0       0         0           0         0         0           0;
        0       0       0       0       1       0       0       0         2*dy    0       0         0           3*dy^2    0         0           0;
        0       0       0       0       1       dx      dx^2    dx^3      2*dy    2*dx*dy 2*dx^2*dy 2*dx^3*dy   3*dy^2    3*dx*dy^2 3*dx^2*dy^2 3*dx^3*dy^2;

        0       0       0       0       0       1       0       0         0       0       0         0           0         0         0           0;
        0       0       0       0       0       1       2*dx    3*dx^2    0       0       0         0           0         0         0           0;
        0       0       0       0       0       1       0       0         0       2*dy    0         0           0         3*dy^2    0           0;
        0       0       0       0       0       1       2*dx    3*dx^2    0       2*dy    4*dx*dy   6*dx^2*dy   0         3*dy^2    6*dx*dy^2   9*dx^2*dy^2];   

    biCubInv=inv(biCubic);

end

%-------------------------------------------------------------------------

Answer 1

I found my error. I pad my matricies with ghost cells, however I forget that now the i in Psi without ghost cells is i+1 in psi with ghost cells. Hence, I should be evaluating my interpolating polynomial P at xi=x(6); yj=y(2), not xi=x(5); yj=y(1).