function F = interp4(x,y,z,u,v)
%INTERP4 2-D bilinear data interpolation.
%	ZI = INTERP4(X,Y,Z,XI,YI) uses bilinear interpolation to
%	find ZI, the values of the underlying 2-D function in Z at the points
%	in matrices XI and YI.  Matrices X and Y specify the points at which 
%	the data Z is given.  X and Y can also be vectors specifying the 
%	abscissae for the matrix Z as for MESHGRID. In both cases, X
%	and Y must be equally spaced and monotonic.
%
%	Values of NaN are returned in ZI for values of XI and YI that are 
%	outside of the range of X and Y.
%
%	If XI and YI are vectors, INTERP4 returns vector ZI containing
%	the interpolated values at the corresponding points (XI,YI).
%
%	ZI = INTERP4(Z,XI,YI) assumes X = 1:N and Y = 1:M, where
%	[M,N] = SIZE(Z).
%
%	F = INTERP4(Z,NTIMES) returns the matrix Z expanded by interleaving
%	bilinear interpolates between every element, working recursively
%	for NTIMES.  INTERP4(Z) is the same as INTERP4(Z,1).
%
%	See also INTERP5, INTERP3.

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

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('XI and YI 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

% Now interpolate
ndx = v+(u-1)*nrows;
F =  ( z(ndx).*(1-t) + z(ndx+1).*t ).*(1-s) + ...
     ( z(ndx+nrows).*(1-t) + z(ndx+(nrows+1)).*t ).*s;

% 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
