## Copyright (C) 2016 Markus Muetzel
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn {} {[@var{faces}, @var{vertices}, @var{J}] =} __unite_shared_vertices__ (@var{faces}, @var{vertices})
##
## Detect and unite shared vertices in patches
##
## Vertices of neighboring faces are detected and united to shared vertices. For
## this, the mutual squared distances between all vertices are calculated. If
## these are smaller than @command{3*eps ()}, the vertices are united to one.
##
## @var{J} holds the indices of the deleted vertices.
##
## @seealso{isosurface, reducepatch}
## @end deftypefn

## Author: mmuetzel

function [faces, vertices, J] = __unite_shared_vertices__ (faces, vertices)
  ### unite shared vertices
  #[vertices, vertices_idx, faces_idx] = unique (vertices, "rows");
  #faces = faces_idx(faces);
  ## unique does not collapse points that are not exactly equal but differ less than eps

  ## Calculate the mutual squared distances between all vertices
  ## This can propably be done more efficiently:
  dx = vertices(:,1) .- vertices(:,1)';
  dy = vertices(:,2) .- vertices(:,2)';
  dz = vertices(:,3) .- vertices(:,3)';
  mutual_dist = dx.^2 + dy.^2 + dz.^2;
  is_very_close = mutual_dist < 6*eps (max (abs (vertices(:))));
  [I, J] = find (is_very_close & triu (true (size (is_very_close)), 1));
  ## always replace to the first unique vertice
  J = flipud(J);
  I = flipud(I);

  num_vertices = size (vertices, 1);
  vertices(J,:) = []; # delete multiple shared vertices
  ## renumber deleted vertices in faces to the one it is replaced by
  vertex_renum = 1:num_vertices;
  vertex_renum(J) = I;
  faces = vertex_renum(faces);
  ## renumber vertices in faces with subsequent numbers
  vertex_renum2 = ones (num_vertices, 1);
  vertex_renum2(J) = 0;
  vertex_renum2 = cumsum (vertex_renum2);
  faces = vertex_renum2(faces);

  ## eliminate identical faces
  faces = sort (faces, 2);
  faces = unique (faces, "rows");

  ## eliminate faces with zero area
  is_zero_area = (faces(:,1) == faces(:,2)) | (faces(:,2) == faces(:,3)); # vertices in faces are sorted
  faces = faces(!is_zero_area,:);

endfunction
