function plywrite(filename,points,varargin)
  fid = fopen(filename,'w');

  % pointsTriangulation = delaunayTriangulation(points);
  % points = pointsTriangulation.Points;
  % facets = getFacetsFromTetrahedrons(pointsTriangulation.ConnectivityList);
  % facets = getFacetsFromTetrahedrons(alphaTriangulation(alphaShape(points),1));
  facets = boundary(points,1) - [1 1 1];

  writeHeader(fid,size(points,1),size(facets,1));
  if nargin == 2
    writePoints(fid,points);
  else
    colormap = varargin{1};
    writePointsWithColormap(fid,points,colormap);
  end
  writeFacets(fid,facets);
  fclose(fid);
  disp('Ply successfully written.');
end


function output = getFacetsFromTetrahedrons(connectivity)
  connectivity = connectivity - [1 1 1 1];
  facets = [connectivity(:,[1 2 3]);connectivity(:,[1 2 4]);connectivity(:,[1 3 4]);connectivity(:,[2 3 4])];
  facets = sort(facets,2);
  facets = sortrows(facets);
  duplicates = find(all(diff(facets) == 0,2));
  facets([duplicates;duplicates + 1],:) = [];
  output = facets;
end


function writeHeader(fid,numberOfPoints,numberOfFacets)
  fprintf(fid,...
    ['ply\n'...
    'format ascii 1.0\n'...
    'element vertex %d\n'...
    'property float32 x\n'...
    'property float32 y\n'...
    'property float32 z\n'...
    'element face %d\n'...
    'property list uint8 int32 vertex_indices\n'...
    'end_header\n'],...
    numberOfPoints,numberOfFacets);
end


function writePoints(fid,points)
  for i = 1:size(points,1)
    fprintf(fid,...
      '%.6f %.6f %.6f\n',...
      points(i,1),points(i,2),points(i,3));
  end
end


function writePointsWithColormap(fid,points,colormap)
  for i = 1:size(points,1)
    fprintf(fid,...
      '%.6f %.6f %.6f %u %u %u\n',...
      points(i,1),points(i,2),points(i,3),...
      colormap(i,1),colormap(i,2),colormap(i,3));
  end
end


function writeFacets(fid,facets)
  for i = 1:size(facets,1)
    fprintf(fid,...
      '%d %d %d %d\n',...
      size(facets,2),facets(i,1),facets(i,2),facets(i,3));
  end
end