function [fout, vout] = isosurface(varargin)
%ISOSURFACE  Isosurface extractor.
%   FV = ISOSURFACE(X,Y,Z,V,ISOVALUE) computes isosurface geometry for
%   data V at isosurface value ISOVALUE. Arrays (X,Y,Z) specify the points
%   at which the data V is given. The struct FV contains the faces and
%   vertices of the isosurface and can be passed directly to the PATCH
%   command.
%   
%   FV = ISOSURFACE(V,ISOVALUE) assumes [X Y Z] = meshgrid(1:N, 1:M, 1:P) 
%        where [M,N,P]=SIZE(V). 
%   
%   FV = ISOSURFACE(X,Y,Z,V) or FV = ISOSURFACE(V) selects an isosurface
%        value automatically using the histogram of the data.
%   
%   FV = ISOSURFACE(..., 'noshare') does not attempt to create shared
%        vertices. This is faster, but produces a larger set of vertices.
%   
%   FV = ISOSURFACE(..., 'verbose') prints progress messages to the command
%        window as the computation progresses. 
%   
%   [F, V] = ISOSURFACE(...) returns the faces and vertices in two arrays
%            instead of a struct.
%   
%   ISOSURFACE(...) With no output arguments, a patch is created with the
%   computed faces and vertices.
%   
%   
%   Example:
%      [x y z v] = flow;
%      p = patch(isosurface(x, y, z, v, -3));
%      isonormals(x,y,z,v, p)
%      set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
%      daspect([1 1 1])
%      view(3)
%      camlight 
%      lighting phong
%
%   See also ISONORMALS, ISOCAPS, SMOOTH3, SUBVOLUME, REDUCEVOLUME,
%            REDUCEPATCH, SHRINKFACES.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 1998/05/13 23:26:58 $

[x y z data value noshare verbose] = parseargs(nargin,varargin);

if length(value)>1
  error('Isovalue must be a scalar.'); 
end

% Take this out when other data types are handled
data = double(data);

[msg x y z] = xyzvcheck(x,y,z,data);  error(msg)

if isempty(value)
  value = isovalue(data);
end

[v f] = isosurf(data, value, noshare, verbose);
v = v';
f = f';
if isempty(v)
  v = [];
  f = [];
end

if ~isempty(x) & ~isempty(v)
  sz = size(x);
  if ~(isequal(x, 1:sz(2)) & isequal(y, 1:sz(1)) & isequal(z, 1:sz(3)))
    nv(:,1) = interp3( x,  v(:,1), v(:,2), v(:,3));
    nv(:,2) = interp3( y,  v(:,1), v(:,2), v(:,3));
    nv(:,3) = interp3( z,  v(:,1), v(:,2), v(:,3));
    v = nv;
  end
end

if nargout==0
  p=patch('faces', f, 'vertices', v, 'facevertexcdata', value, ...
          'facecolor', 'flat', 'edgecolor', 'none', 'userdata', value);
  if ~isempty(x)
    isonormals(x,y,z,data, p);
  else
    isonormals(data, p);
  end
elseif nargout==1
  fout.vertices = v;
  fout.faces = f;
else
  fout = f;
  vout = v;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y, z, data, value, noshare, verbose] = parseargs(nin, vargin)


x = [];
y = [];
z = [];
data = [];
value = [];
noshare = 0;
verbose = 0;

for j = 1:2
  if nin>0
    lastarg = vargin{nin};
    if isstr(lastarg)
      if ~isempty(lastarg)
	if lower(lastarg(1))=='n' % noshare
	  noshare = 1;
	end
	if lower(lastarg(1))=='v' % verbose
	  verbose = 1;
	end
      end
      nin = nin - 1;
    end
  end
end

if nin==1                % isosurface(v)
  data = vargin{1};
elseif nin==2            % isosurface(v, isoval)
  data = vargin{1};
  value = vargin{2};
elseif nin==4            % isosurface(x,y,z,v)
  x = vargin{1};
  y = vargin{2};
  z = vargin{3};
  data = vargin{4};
elseif nin==5            % isosurface(x,y,z,v, isovalue)
  x = vargin{1};
  y = vargin{2};
  z = vargin{3};
  data = vargin{4};
  value = vargin{5};
else
  error('Wrong number of input arguments.'); 
end

