THESE/codes/snake/gvf_dist_v4.2c/interp4.m

   1 function F = interp4(arg1,arg2,arg3,arg4,arg5)
   2 %INTERP4 2-D bilinear data interpolation.
   3 %       ZI = INTERP4(X,Y,Z,XI,YI) uses bilinear interpolation to
   4 %       find ZI, the values of the underlying 2-D function in Z at the points
   5 %       in matrices XI and YI.  Matrices X and Y specify the points at which
   6 %       the data Z is given.  X and Y can also be vectors specifying the
   7 %       abscissae for the matrix Z as for MESHGRID. In both cases, X
   8 %       and Y must be equally spaced and monotonic.
   9 %
  10 %       Values of NaN are returned in ZI for values of XI and YI that are
  11 %       outside of the range of X and Y.
  12 %
  13 %       If XI and YI are vectors, INTERP4 returns vector ZI containing
  14 %       the interpolated values at the corresponding points (XI,YI).
  15 %
  16 %       ZI = INTERP4(Z,XI,YI) assumes X = 1:N and Y = 1:M, where
  17 %       [M,N] = SIZE(Z).
  18 %
  19 %       ZI = INTERP4(Z,NTIMES) returns the matrix Z expanded by interleaving
  20 %       bilinear interpolates between every element, working recursively
  21 %       for NTIMES.  INTERP4(Z) is the same as INTERP4(Z,1).
  22 %
  23 %       This function needs about 4 times SIZE(XI) memory to be available.
  24 %
  25 %       See also INTERP2, INTERP5.
  26
  27 %       Copyright (c) 1984-94 by The MathWorks, Inc.
  28 %       Clay M. Thompson 4-26-91, revised 7-3-91, 3-22-93 by CMT.
  29
  30 if nargin==1, % interp4(z), Expand Z
  31   [nrows,ncols] = size(arg1);
  32   s = 1:.5:ncols; sizs = size(s);
  33   t = (1:.5:nrows)'; sizt = size(t);
  34   s = s(ones(sizt),:);
  35   t = t(:,ones(sizs));
  36
  37 elseif nargin==2, % interp4(z,n), Expand Z n times
  38   [nrows,ncols] = size(arg1);
  39   ntimes = floor(arg2);
  40   s = 1:1/(2*ntimes):ncols; sizs = size(s);
  41   t = (1:1/(2*ntimes):nrows)'; sizt = size(t);
  42   s = s(ones(sizt),:);
  43   t = t(:,ones(sizs));
  44
  45 elseif nargin==3, % interp4(z,s,t), No X or Y specified.
  46   [nrows,ncols] = size(arg1);
  47   s = arg2; t = arg3;
  48
  49 elseif nargin==4,
  50   error('Wrong number of input arguments.');
  51
  52 elseif nargin==5, % interp4(x,y,z,s,t), X and Y specified.
  53   [nrows,ncols] = size(arg3);
  54   mx = prod(size(arg1)); my = prod(size(arg2));
  55   if any([mx my] ~= [ncols nrows]) & any(size(arg1)~=size(arg3) | ...
  56      size(arg2)~=size(arg3)),
  57     error('The lengths of the X and Y vectors must match Z.');
  58   end
  59   s = 1 + (arg4-arg1(1))*((ncols-1)/(arg1(mx)-arg1(1)));
  60   t = 1 + (arg5-arg2(1))*((nrows-1)/(arg2(my)-arg2(1)));
  61
  62 end
  63
  64 if any([nrows ncols]<[2 2]), error('Z must be at least 2-by-2.'); end
  65 if any(size(s)~=size(t)),
  66   error('XI and YI must be the same size.');
  67 end
  68
  69 % Check for out of range values of s and set to 1
  70 sout = find((s<1)|(s>ncols));
  71 if length(sout)>0, s(sout) = ones(size(sout)); end
  72
  73 % Check for out of range values of t and set to 1
  74 tout = find((t<1)|(t>nrows));
  75 if length(tout)>0, t(tout) = ones(size(tout)); end
  76
  77 % Matrix element indexing
  78 ndx = floor(t)+floor(s-1)*nrows;
  79
  80 % Compute intepolation parameters, check for boundary value.
  81 d = find(s==ncols);
  82 s(:) = (s - floor(s));
  83 if length(d)>0, s(d) = s(d)+1; ndx(d) = ndx(d)-nrows; end
  84
  85 % Compute intepolation parameters, check for boundary value.
  86 d = find(t==nrows);
  87 t(:) = (t - floor(t));
  88 if length(d)>0, t(d) = t(d)+1; ndx(d) = ndx(d)-1; end
  89 d = [];
  90
  91 % Now interpolate, reuse u and v to save memory.
  92 if nargin==5,
  93   F =  ( arg3(ndx).*(1-t) + arg3(ndx+1).*t ).*(1-s) + ...
  94        ( arg3(ndx+nrows).*(1-t) + arg3(ndx+(nrows+1)).*t ).*s;
  95 else
  96   F =  ( arg1(ndx).*(1-t) + arg1(ndx+1).*t ).*(1-s) + ...
  97        ( arg1(ndx+nrows).*(1-t) + arg1(ndx+(nrows+1)).*t ).*s;
  98 end
  99
 100 % Now set out of range values to NaN.
 101 if length(sout)>0, F(sout) = NaN*ones(size(sout)); end
 102 if length(tout)>0, F(tout) = NaN*ones(size(tout)); end