function hout = slice(varargin)
%SLICE  Volumetric slice plot.
%   SLICE(X,Y,Z,V,Sx,Sy,Sz) draws slices along the x,y,z directions at
%   the points in the vectors Sx,Sy,Sz. The arrays X,Y,Z define the
%   coordinates for V and must be monotonic and 3-D plaid (as if
%   produced by MESHGRID).  The color at each point will be determined
%   by 3-D interpolation into the volume V.  V must be an M-by-N-by-P
%   volume array. 
%
%   SLICE(X,Y,Z,V,XI,YI,ZI) draws slices through the volume V along the
%   surface defined by the arrays XI,YI,ZI.
%
%   SLICE(V,Sx,Sy,Sz) or SLICE(V,XI,YI,ZI) assumes X=1:N, Y=1:M, Z=1:P. 
%
%   SLICE(...,'method') specifies the interpolation method to use.
%   'method' can be 'linear', 'cubic', or 'nearest'.  'linear' is the
%   default (see INTERP3).
%
%   H = SLICE(...) returns a vector of handles to SURFACE objects.
%
%   Example: To visualize the function x*exp(-x^2-y^2-z^2) over the
%   range -2 < x < 2, -2 < y < 2, -2 < z < 2, 
%
%      [x,y,z] = meshgrid(-2:.2:2, -2:.25:2, -2:.16:2);
%      v = x .* exp(-x.^2 - y.^2 - z.^2);
%      slice(x,y,z,v,[-1.2 .8 2],2,[-2 -.2])
%
%   See also MESHGRID, INTERP3.

%   J.N. Little 1-23-92
%   Revised 4-27-93, 2-10-94
%   Revised 6-17-94 by Clay M. Thompson
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.16 $  $Date: 1997/11/21 23:46:47 $

error(nargchk(4,8,nargin));

cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
axis_slice = 0; % Slice only along axes?
method = 'linear';

nin = nargin;
if isstr(varargin{nin}), % slice(...,'method')
  method = varargin{nin};
  nin = nin - 1;
end

if nin == 4 | nin == 5 % slice(v,xi,yi,zi,<nx>)
    v = varargin{1}; varargin{1} = [];

    if nin == 5,
      warning(sprintf([ ...
        'Specifying the number of columns is no longer necessary.\n',...
        '         Use slice(v,xi,yi,zi) instead.']))
      if ndims(v)==2, % Reshape the volume.
        v = reshape(v,size(v,1)/varargin{5},varargin{5},size(v,2));
      end
    end

    if ndims(v)~=3, error('V must be a 3-D array.'); end
    [ny,nx,nz] = size(v);
    if min(size(varargin{2}))<=1 & min(size(varargin{3}))<=1 & ...
       min(size(varargin{4}))<=1,
        axis_slice = 1;
    end
    sx = varargin{2};
    sy = varargin{3};
    sz = varargin{4};
    [x,y,z] = meshgrid(1:nx,1:ny,1:nz);
elseif nin == 6,
    error('Wrong number of input arguments.');
elseif nin == 7 | nin == 8, % slice(x,y,z,v,xi,yi,zi,<nx>)
    v = varargin{4}; varargin{4} = [];

    if nin == 8,
      warning(sprintf([ ...
        'Specifying the number of columns is no longer necessary.\n',...
        '         Use slice(x,y,z,v,xi,yi,zi) instead.']))
      if ndims(v)==2, % Reshape the volume
        v = reshape(v,size(v,1)/varargin{8},varargin{8},size(v,2));
      end
    end

    if ndims(v)~=3, error('V must be a 3-D array.'); end
    if min(size(varargin{1}))==1 | min(size(varargin{2}))==1 | ...
       min(size(varargin{3}))==1,
        [x,y,z] = meshgrid(varargin{1:3});
    else
        x = varargin{1};
        y = varargin{2};
        z = varargin{3};
    end
    [ny,nx,nz] = size(v);

    if min(size(varargin{5}))<=1 & min(size(varargin{6}))<=1 & ...
       min(size(varargin{7}))<=1,
        axis_slice = 1;
    end
    sx = varargin{5};
    sy = varargin{6};
    sz = varargin{7};
end

if axis_slice,
    h = [];
    [xi,yi,zi] = meshgrid(sx,y(:,1,1),z(1,1,:));
    vi = interp3(x,y,z,v,xi,yi,zi,method);
    for i = 1:length(sx)
        h = [h; surface( ...
            reshape(xi(:,i,:),[ny nz]),reshape(yi(:,i,:),[ny nz]), ...
            reshape(zi(:,i,:),[ny nz]),reshape(vi(:,i,:),[ny nz]))];
    end

    [xi,yi,zi] = meshgrid(x(1,:,1),sy,z(1,1,:));
    vi = interp3(x,y,z,v,xi,yi,zi,method);
    for i = 1:length(sy)
        h = [h; surface( ...
            reshape(xi(i,:,:),[nx nz]),reshape(yi(i,:,:),[nx nz]),...
            reshape(zi(i,:,:),[nx nz]),reshape(vi(i,:,:),[nx nz]))];
    end

    [xi,yi,zi] = meshgrid(x(1,:,1),y(:,1,1),sz);
    vi = interp3(x,y,z,v,xi,yi,zi,method);
    for i = 1:length(sz)
        h = [h; surface(xi(:,:,i),yi(:,:,i),zi(:,:,i),vi(:,:,i))];
    end
else
    vi = interp3(x,y,z,v,sx,sy,sz,method);
    h = surf(sx,sy,sz,vi);
end

if nargout > 0
    hout = h;
end
if ~hold_state
    view(3), grid on
end

% Use ISFINITE to make sure no NaNs or Infs get passed to CAXIS
u=v(isfinite(v)); u = u(:);
caxis([min(u) max(u)])
