function F = interp5(x,y,z,u,v)
%INTERP5 2-D bicubic data interpolation.
%	INTERP5(...) is the same as INTERP4(....) except that it uses
%	bicubic interpolation.
%	
%	See also INTERP4.

%	Copyright (c) 1984-93 by The MathWorks, Inc.
%	Clay M. Thompson 4-26-91, revised 7-3-91 by CMT.

%	Based on "Cubic Convolution Interpolation for Digital Image
%	Processing", Robert G. Keys, IEEE Trans. on Acoustics, Speech, and
%	Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160.

if nargin==1, % Expand Z
  z = x;
  [nrows,ncols] = size(z);
  u = ones(2*nrows-1,1)*[1:.5:ncols];
  v = [1:.5:nrows]'*ones(1,2*ncols-1);

elseif nargin==2, % Expand Z ntimes
  z = x;
  ntimes = floor(y);
  [nrows,ncols] = size(z);
  u = ones(2*ntimes*(nrows-1)+1,1)*[1:1/(2*ntimes):ncols];
  v = [1:1/(2*ntimes):nrows]'*ones(1,2*ntimes*(ncols-1)+1);

elseif nargin==3, % No X or Y specified.
  v = z; u = y; z = x;
  [nrows,ncols] = size(z);

elseif nargin==5, % X and Y specified.
  [nrows,ncols] = size(z);
  mx = prod(size(x)); my = prod(size(y));
  if any([mx my] ~= [ncols nrows]) & (size(x)~=size(z) | size(y)~=size(z)), 
    error('The lengths of the X and Y vectors must match Z.');
  end
  u = 1 + (u-x(1))*((ncols-1)/(x(mx)-x(1)));
  v = 1 + (v-y(1))*((nrows-1)/(y(my)-y(1)));
  
else
  error('Wrong number of input arguments.');
end

if any(size(z)<[3 3]), error('Z must be at least 3-by-3.'); end
if size(u)~=size(v), error('X1 and Y1 must be the same size.'); end

% Check for out of range values of u and set to 1
uout = (u<1)|(u>ncols);
nuout = sum(uout(:));
if any(uout(:)), u(uout) = ones(nuout,1); end

% Check for out of range values of v and set to 1
vout = (v<1)|(v>nrows);
nvout = sum(vout(:));
if any(vout(:)), v(vout) = ones(nvout,1); end

% Interpolation parameters
s = (u - floor(u));  t = (v - floor(v));
u = floor(u); v = floor(v);
d = (u==ncols); if any(d(:)), u(d) = u(d)-1; s(d) = s(d)+1; end
d = (v==nrows); if any(d(:)), v(d) = v(d)-1; t(d) = t(d)+1; end

% Expand z so interpolation is valid at the boundaries.
z = [3*z(1,:)-3*z(2,:)+z(3,:);z;3*z(nrows,:)-3*z(nrows-1,:)+z(nrows-2,:)];
z = [3*z(:,1)-3*z(:,2)+z(:,3),z,3*z(:,ncols)-3*z(:,ncols-1)+z(:,ncols-2)];
nrows = nrows+2; ncols = ncols+2;

% Now interpolate using computationally efficient algorithm.
s2 = s.*s; s3 = s.*s2;
t2 = t.*t; t3 = t.*t2;
ndx = v+(u-1)*nrows;
F     = ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
    z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
    (-s3 + 2*s2 - s);
ndx = ndx + nrows;
F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
    z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
    (3*s3 - 5*s2 + 2);
ndx = ndx + nrows;
F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
    z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
    (-3*s3 + 4*s2 + 1*s);
ndx = ndx + nrows;
F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
    z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
    (s3 - s2);
F = F/4;

% Now set out of range values to NaN.
if any(uout(:)), F(uout) = NaN*ones(nuout,1); end
if any(vout(:)), F(vout) = NaN*ones(nvout,1); end


