home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / ICUBIC.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.7 KB  |  80 lines
  1. function F=icubic(x,y,u)
  2. %ICUBIC Cubic Interpolation of a 1-D function.
  3. %
  4. %    YI=ICUBIC(Y,XI) returns the value of the 1-D function Y at the
  5. %    points XI using cubic interpolation. length(YI)=length(XI). XI is
  6. %    an index into the vector Y. Y is the value of the function
  7. %    evaluated uniformly on a interval. If Y is a matrix, then
  8. %    the interpolation is performed for each column of Y.
  9. %
  10. %    If Y is of length N then XI must contain values between 1 and N.
  11. %    The value NaN is returned if this is not the case.
  12. %
  13. %    YI = ICUBIC(X,Y,XI) uses the vector X to specify the coordinates
  14. %    of the underlying interval. X must be equally spaced and
  15. %    monotonic. XI must lie within the coordinates in X.
  16.  
  17. %    Clay M. Thompson 7-4-91
  18. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  19.  
  20. %    Based on "Cubic Convolution Interpolation for Digital Image
  21. %    Processing", Robert G. Keys, IEEE Trans. on Acoustics, Speech, and
  22. %    Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160.
  23.  
  24. if nargin==2,    % No X specified.
  25.   u = y; y = x;
  26.  
  27. % Check for vector problem.  If so, make everything a column vector.
  28.   if min(size(y))==1, y = y(:); end
  29.   if min(size(u))==1, u = u(:); end
  30.   [nrows,ncols] = size(y);
  31.  
  32. elseif nargin==3, % X specified.
  33.   % Check for vector problem.  If so, make everything a column vector.
  34.   if min(size(y))==1, y = y(:); end
  35.   if min(size(x))==1, x = x(:); end
  36.   if min(size(u))==1, u = u(:); end
  37.   [nrows,ncols] = size(y);
  38.   % Scale and shift u to be indices into Y.
  39.   if (min(size(x))~=1), error('X must be a vector.'); end
  40.   x = x(:);
  41.   [m,n] = size(x);
  42.   if m ~= nrows, 
  43.     error('The length of X must match the number of rows of Y.');
  44.   end
  45.   u = 1 + (u-x(1))*((nrows-1)/(x(m)-x(1)));
  46.   
  47. else
  48.   error('Wrong number of input arguments.');
  49. end
  50.  
  51. if nrows<3, error('Y must have at least 3 rows.'); end
  52. [m,n] = size(u); 
  53. if n==1, u = u*ones(1,ncols); [m,n] = size(u); end    % Expand u 
  54. if n~=ncols, error('The number of columns in XI and Y must match.'); end
  55.  
  56. % Check for out of range values of u and set to 1
  57. uout = find((u<1)|(u>nrows));
  58. nuout = length(uout);
  59. if nuout>0, u(uout) = ones(nuout,1); end
  60.  
  61. % Interpolation parameters
  62. s = (u - floor(u));
  63. u = floor(u);
  64. d = find(u==nrows); if length(d)>0, u(d) = u(d)-1; s(d) = s(d)+1; end
  65.  
  66. % Expand y so interpolation is valid at the boundary.
  67. y = [3*y(1,:)-3*y(2,:)+y(3,:);y;3*y(nrows,:)-3*y(nrows-1,:)+y(nrows-2,:)];
  68. nrows = nrows + 2;
  69.  
  70. % Now interpolate using computationally efficient algorithm.
  71. s2 = s.*s; s3 = s.*s2;
  72. ndx = u+ones(m,1)*[0:n-1]*nrows;
  73. F     = y(ndx).*(-s3+2*s2-s) + y(ndx+1).*(3*s3-5*s2+2) + ...
  74.     y(ndx+2).*(-3*s3+4*s2+s) + y(ndx+3).*(s3-s2);
  75. F = F/2;
  76.  
  77. % Now set out of range values to NaN.
  78. if nuout>0, F(uout) = NaN*ones(nuout,1); end
  79.  
  80.