home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / INTERP5.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.1 KB  |  95 lines
  1. function F = interp5(x,y,z,u,v)
  2. %INTERP5 2-D bicubic data interpolation.
  3. %    INTERP5(...) is the same as INTERP4(....) except that it uses
  4. %    bicubic interpolation.
  5. %    
  6. %    See also INTERP4.
  7.  
  8. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  9. %    Clay M. Thompson 4-26-91, revised 7-3-91 by CMT.
  10.  
  11. %    Based on "Cubic Convolution Interpolation for Digital Image
  12. %    Processing", Robert G. Keys, IEEE Trans. on Acoustics, Speech, and
  13. %    Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160.
  14.  
  15. if nargin==1, % Expand Z
  16.   z = x;
  17.   [nrows,ncols] = size(z);
  18.   u = ones(2*nrows-1,1)*[1:.5:ncols];
  19.   v = [1:.5:nrows]'*ones(1,2*ncols-1);
  20.  
  21. elseif nargin==2, % Expand Z ntimes
  22.   z = x;
  23.   ntimes = floor(y);
  24.   [nrows,ncols] = size(z);
  25.   u = ones(2*ntimes*(nrows-1)+1,1)*[1:1/(2*ntimes):ncols];
  26.   v = [1:1/(2*ntimes):nrows]'*ones(1,2*ntimes*(ncols-1)+1);
  27.  
  28. elseif nargin==3, % No X or Y specified.
  29.   v = z; u = y; z = x;
  30.   [nrows,ncols] = size(z);
  31.  
  32. elseif nargin==5, % X and Y specified.
  33.   [nrows,ncols] = size(z);
  34.   mx = prod(size(x)); my = prod(size(y));
  35.   if any([mx my] ~= [ncols nrows]) & (size(x)~=size(z) | size(y)~=size(z)), 
  36.     error('The lengths of the X and Y vectors must match Z.');
  37.   end
  38.   u = 1 + (u-x(1))*((ncols-1)/(x(mx)-x(1)));
  39.   v = 1 + (v-y(1))*((nrows-1)/(y(my)-y(1)));
  40.   
  41. else
  42.   error('Wrong number of input arguments.');
  43. end
  44.  
  45. if any(size(z)<[3 3]), error('Z must be at least 3-by-3.'); end
  46. if size(u)~=size(v), error('X1 and Y1 must be the same size.'); end
  47.  
  48. % Check for out of range values of u and set to 1
  49. uout = (u<1)|(u>ncols);
  50. nuout = sum(uout(:));
  51. if any(uout(:)), u(uout) = ones(nuout,1); end
  52.  
  53. % Check for out of range values of v and set to 1
  54. vout = (v<1)|(v>nrows);
  55. nvout = sum(vout(:));
  56. if any(vout(:)), v(vout) = ones(nvout,1); end
  57.  
  58. % Interpolation parameters
  59. s = (u - floor(u));  t = (v - floor(v));
  60. u = floor(u); v = floor(v);
  61. d = (u==ncols); if any(d(:)), u(d) = u(d)-1; s(d) = s(d)+1; end
  62. d = (v==nrows); if any(d(:)), v(d) = v(d)-1; t(d) = t(d)+1; end
  63.  
  64. % Expand z so interpolation is valid at the boundaries.
  65. z = [3*z(1,:)-3*z(2,:)+z(3,:);z;3*z(nrows,:)-3*z(nrows-1,:)+z(nrows-2,:)];
  66. z = [3*z(:,1)-3*z(:,2)+z(:,3),z,3*z(:,ncols)-3*z(:,ncols-1)+z(:,ncols-2)];
  67. nrows = nrows+2; ncols = ncols+2;
  68.  
  69. % Now interpolate using computationally efficient algorithm.
  70. s2 = s.*s; s3 = s.*s2;
  71. t2 = t.*t; t3 = t.*t2;
  72. ndx = v+(u-1)*nrows;
  73. F     = ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
  74.     z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
  75.     (-s3 + 2*s2 - s);
  76. ndx = ndx + nrows;
  77. F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
  78.     z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
  79.     (3*s3 - 5*s2 + 2);
  80. ndx = ndx + nrows;
  81. F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
  82.     z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
  83.     (-3*s3 + 4*s2 + 1*s);
  84. ndx = ndx + nrows;
  85. F(:)  = F + ( z(ndx).*(-t3+2*t2-t) + z(ndx+1).*(3*t3-5*t2+2) + ...
  86.     z(ndx+2).*(-3*t3+4*t2+t) + z(ndx+3).*(t3-t2) ) .* ...
  87.     (s3 - s2);
  88. F = F/4;
  89.  
  90. % Now set out of range values to NaN.
  91. if any(uout(:)), F(uout) = NaN*ones(nuout,1); end
  92. if any(vout(:)), F(vout) = NaN*ones(nvout,1); end
  93.  
  94.  
  95.