function z = CubicInterp2D(xc,yc,a,b,c,d,fA)
% Bicubic interpolation on a rectangular grid
% of a function f(x,y)
%
% a <= xc <= b; c <= yc <= d
% [n,m] = size(fA)
% dx=(b-a)/(n-1); dy=(d-c)/(m-1)
% fA(i,j)=f(a+(i-1)dx,c+(j-1)dy)
%
% z=S(xc,yc) where S is the bicubic interpolant
% of f using the nearest 4X4 subgrid
 
[n,m,lc] = size(fA);
 
% Check for consistency
 
if (n<4) | (m<4)
   disp(' insufficent data ')
   return
end
 
dx=(b-a)/(n-1);
dy=(d-c)/(m-1);
 
% Now we find the indices of the lower left corner
% of the subgrid
 
i=ceil((xc-a)/dx)-1;
if i<1
   i=1;
end
if i+3 > n
   i=n-3;
end
j=ceil((yc-c)/dy)-1;
if j<1
   j=1;
end
if j+3 > m
   j=m-3;
end
 
fI=zeros(4,1);
 
% First compute interpolants in x for fixed y
xd=a*ones(4,1)+dx*linspace(i-1,i+2,4)';
for k=1:4
   fx=double(fA(i:i+3,j+k-1));
   cc=InterpNRecur(xd,fx);
   fI(k,1)=HornerN(cc,xd,xc);
end
 
%Now interpolate in y
yd=c*ones(4,1)+dy*linspace(j-1,j+2,4)';
cc=InterpNRecur(yd,fI);
z=uint8(HornerN(cc,yd,yc));