myshrinkfaces(varargin)

function [fout, vout] = myshrinkfaces(varargin)
%MYSHRINKFACES  Reduce size of patch faces.
%   MYSHRINKFACES(P, SF) shrinks the area of the faces in patch P to shrink
%   factor SF. If SF is 0.6, each face is shrunk to 60% of its original
%   area. If the patch contains shared vertices, non-shared vertices are
%   created before reduction.
%   
%   NFV = MYSHRINKFACES(P, SF) returns the faces and vertices but does not
%   set the Faces and Vertices properties of patch P.  The struct NFV
%   contains the new faces and vertices.
%   
%   NFV = MYSHRINKFACES(FV, SF) uses faces and vertices from struct FV.
%
%   MYSHRINKFACES(P) or MYSHRINKFACES(FV) assumes a shrink factor of .3.
%   
%   NFV = MYSHRINKFACES(F, V, SF) uses faces and vertices in arrays F and V.
%   
%   [NF, NV] = MYSHRINKFACES(...) returns the faces and vertices in two
%              arrays instead of a struct.
%   Example:
%      [x y z v] = flow;
%      [x y z v] = reducevolume(x,y,z,v, 2);
%      fv = isosurface(x, y, z, v, -3);
%      subplot(1,2,1)
%      p = patch(fv);
%      set(p, 'facecolor', [.5 .5 .5], 'EdgeColor', 'black');
%      daspect([1 1 1]); view(3); axis tight
%      title('Original')
%      subplot(1,2,2)
%      p = patch(myshrinkfaces(fv, .2)); % shrink faces to 20% of original 
%      set(p, 'facecolor', [.5 .5 .5], 'EdgeColor', 'black');
%      daspect([1 1 1]); view(3); axis tight
%      title('After Shrinking')
%   
%   See also ISOSURFACE, ISONORMALS, ISOCAPS, SMOOTH3, SUBVOLUME, 
%            REDUCEVOLUME, REDUCEPATCH.

%   Copyright 1984-2005 The MathWorks, Inc. 
%   $Revision: 1.7.4.1 $  $Date: 2005/04/28 19:56:57 $

[p faces verts sf] = parseargs(nargin,varargin);

oldNV = size(verts,1);

if length(sf)>1
  error(id('NonScalarFactor'),'Shrink factor must be a scalar.'); 
end
if sf<0
  error(id('NonPositiveFactor'),'Shrink factor must be positive.'); 
end


if isempty(sf)
  sf = sqrt(.3);
end

nanindex = isnan(faces);
[faces verts newVColorIdx] = facesvertsnoshare(faces, verts);

fcols = size(faces,2);
coords = verts(faces.',:);
facexyz = reshape(coords,fcols,numel(coords)/fcols);

av = nanmean(facexyz);
av = repmat(av,[fcols 1]);

facexyz = facexyz*sf - av*(sf-1);

verts(faces',:) = reshape(facexyz, size(coords));
faces(nanindex) = nan; 

if nargout==0
  if ~isempty(p)
	% if the original patch was not formed with faces/vertices properties,
	% the data in FaceVertexCdata has oldNV rows, and will be useless here
	% as we have now 3*new_num_faces vertices
	oldVCdat = get(p,'FaceVertexCData');
	if (oldNV~= size(verts,1))  % vert num has increased indeed
	    newVCdat = oldVCdat(newVColorIdx,:);
	    set(p, 'faces', faces, 'vertices', verts,'FaceVertexCData',newVCdat);
	else
	    set(p, 'faces', faces, 'vertices', verts);
	end 
    
  else
    fout.faces = faces;
    fout.vertices = verts;
  end
elseif nargout==1
  fout.faces = faces;
  fout.vertices = verts;
else
  fout = faces;
  vout = verts;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p, faces, verts, sf] = parseargs(nin, vargin)

p=[];
sf = [];

if nin==1 || nin==2           % myshrinkfaces(p), myshrinkfaces(fv), myshrinkfaces(arg, sf) 
  firstarg = vargin{1};
  if strcmp(class(firstarg), 'struct')
    faces = firstarg.faces;
    verts = firstarg.vertices;
  elseif all(ishandle(firstarg)) && all(strcmp(get(firstarg, 'type'), 'patch'))
    p = firstarg;
    faces = get(p, 'faces');
    verts = get(p, 'vertices');
  else
    error(id('InvalidFirstArgument'),'The first argument must be a patch handle or a struct containing faces and vertices.');
  end
  if nin==2
    sf = vargin{2};
  end
elseif nin==3            %myshrinkfaces(f, v, sf)
  faces = vargin{1};
  verts = vargin{2};
  sf = vargin{3};
else
  error(id('WrongNumberOfInputs'),'Wrong number of input arguments.'); 
end

if ~isempty(sf)
  if sf>=0
    sf = sqrt(sf);
  else
    sf = -sqrt(-sf);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [newf, newv, newVColorIdx]=facesvertsnoshare(f, v)

fcols = size(f,2);
fmax = 1+max(f(:));
nanindex = isnan(f);
f(nanindex)=fmax;
findex = f';
v(fmax,:) = nan*zeros(1,size(v,2));

newv = v(findex,:);
vrows = size(newv,1);
newf = reshape(1:vrows, fcols, vrows/fcols).';

oldc = 1:size(v,1);
newVColorIdx = oldc(findex);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nanmean(x)
%NANMEAN Average or mean ignoring NaNs.
%   NANMEAN(X) returns the average treating NaNs as missing values.  
%   For vectors, NANMEAN(X) is the mean value of the non-NaN
%   elements in X.  For matrices, NANMEAN(X) is a row vector
%   containing the mean value of each column, ignoring NaNs.

if isempty(x) % Check for empty input.
    y = NaN;
    return
end

% Replace NaNs with zeros.
nans = isnan(x);
i = find(nans);
x(i) = zeros(size(i));

if min(size(x))==1,
  count = length(x)-sum(nans);
else
  count = size(x,1)-sum(nans);
end

% Protect against a column of all NaNs
i = find(count==0);
count(i) = ones(size(i));
y = sum(x)./count;
y(i) = i + NaN;

function str=id(str)
str = ['MATLAB:myshrinkfaces:' str];