function [faces,newVertices] = taubin(srf,rot1,rot2,lambda,mu,levelSet,iterations,reduce)

if(reduce==1)
 fv = reducepatch(isosurface(srf,levelSet))
else
 fv = isosurface(srf,levelSet)
end

szf = size(fv.faces);
szv = size(fv.vertices);

vertices = fv.vertices;
faces = fv.faces;
newVertices = vertices;
for i=1:iterations
for step=1:2

  if(step==1)
    scaleFactor=lambda
  else
    scaleFactor=mu
  end

  for v=1:szv(1)
  vertex = vertices(v,:);
    %find indices of all faces containing the vth vertex
    [row col] = find(faces==v);
    occurences = size(row);
    numNeighbors=1;
    for n = 1:occurences(1)
        if(col(n)==1)
         neighbors(numNeighbors)=faces(row(n),2);
         neighbors(numNeighbors+1)=faces(row(n),3);
        elseif(col(n)==2)
         neighbors(numNeighbors)=faces(row(n),1);
         neighbors(numNeighbors+1)=faces(row(n),3);
        else
         neighbors(numNeighbors)=faces(row(n),1);
         neighbors(numNeighbors+1)=faces(row(n),1);
        end
        numNeighbors = numNeighbors+2;
    end
    neighbors = unique(neighbors);

   numNeighbors = size(neighbors);
   points = [];
   for nextV = 1:numNeighbors(2)
     points(nextV,:) = vertices(neighbors(nextV),:)-vertex;
   end
   
   displacement = (sum(points))./(numNeighbors(2));
   newVertices(v,:) = vertices(v,:) + displacement.*scaleFactor;
  end
  vertices = newVertices;
end
end

fv.vertices = newVertices;
figure(2);
hiso=patch(fv,'FaceColor',[1,.75,.65],'EdgeColor','none'); 
daspect([1,1,1]);
view(rot1,rot2); lightangle(rot1,rot2); set(gcf,'Renderer','zbuffer'); lighting phong