## 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{reduced_fv} =} reducepatch (@var{fv})
## @deftypefnx {} {@var{reduced_fv} =} reducepatch (@var{faces}, @var{vertices})
## @deftypefnx {} {@var{reduced_fv} =} reducepatch (@var{patch_handle})
## @deftypefnx {} reducepatch (@var{patch_handle})
## @deftypefnx {} {@var{reduced_fv} =} reducepatch (@dots{}, @var{reduction_factor})
## @deftypefnx {} {@var{reduced_fv} =} reducepatch (@dots{}, "fast")
## @deftypefnx {} {@var{reduced_fv} =} reducepatch (@dots{}, "verbose")
## @deftypefnx {} {[@var{reduced_faces}, @var{reduces_vertices}] =} reducepatch (@dots{})
##
## Reduce number of patch faces and vertices
##
## @command{reducepatch} can be called with a structure @var{fv} with the fields
## @code{faces} and @code{vertices}, with two matrices @var{faces},
## @var{vertices} (see e.g. @command{isosurface}) or a valid handle to a patch
## object @var{patch_handle} (see: @command{patch}).
##
## The number of faces and vertices in the patch is reduced maintaining the
## overall shape of the patch. This is done by iteratively collapsing the
## shortest edge of the patch to its midpoint (as discussed e.g. here:
## @url{http://libigl.github.io/libigl/tutorial/tutorial.html#meshdecimation}).
##
## Currently, only patches consisting of triangles are supported. The resulting
## patch consists of only triangles, too.
##
## If @command{reducepatch} is called with a handle to a valid patch
## @var{patch_handle} and without output arguments, the given patch is updated
## immediately in its axes.
##
## If the @var{reduction_factor} is omitted, the resulting structure
## @var{reduced_fv} includes approximately 50% of the faces of the original
## patch.  If @var{reduction_factor} is a fraction between 0 (excluded) and 1
## (included), a patch with the corresponding fraction of faces is determined.  
## If @var{reduction_factor} is an integer larger than one, the resulting patch
## has approximately @var{reduction_factor} faces.  Depending on the geometry of
## the patch, the resulting number of faces can differ from the given value of
## @var{reduction_factor}. This is especially true when many shared vertices are
## detected.
##
## For the reduction, it is necessary that vertices of touching faces are
## shared. Shared vertices are detected automatically. This detection can be
## skipped by passing the optional string argument @qcode{"fast"}.  
## With the optional string arguments @qcode{"verbose"}, additional status
## messages are printed to the command window.
## The string input arguments must be passed after the other arguments.
##
## If called with two output arguments, the reduced faces and vertices are
## returned in two separate matrices @var{reduced_faces} and
## @var{reduces_vertices}.
##
## @seealso{isosurface, isocaps, isonormals, reducevolume, smooth3, patch}
## @end deftypefn

## Author: mmuetzel

## FIXME: Convert faces to only triangles if necessary

function varargout = reducepatch (varargin)

  if (nargin < 1 || nargin > 5 || nargout > 2)
    print_usage ();
  endif

  [faces, vertices, max_faces, patch_handle, fast, verbose] = __get_check_reducepatch_args__ (varargin{:});

  if verbose
    fprintf ("reducepatch: before reduction of faces:           %5d faces, %5d vertices\n", size (faces, 1), size (vertices, 1))
  endif

  if (!fast)
    [faces, vertices] = __unite_shared_vertices__ (faces, vertices);

    if verbose
      fprintf ("reducepatch: after detection of shared vertices:  %5d faces, %5d vertices\n", size (faces, 1), size (vertices, 1))
    endif
  endif

  [faces_reduced, vertices_reduced] = __reducepatch__ (vertices, faces, max_faces, verbose);

  if verbose
    fprintf ("reducepatch: after reduction of faces:            %5d faces, %5d vertices\n", size (faces_reduced, 1), size (vertices_reduced, 1))
  endif

  ## output
  if (!isempty (patch_handle) && nargout == 0)
    ## update patch object
    set (patch_handle, "Faces", faces_reduced, "Vertices", vertices_reduced);
  elseif (nargout == 2) # faces, vertices
    varargout{1} = faces_reduced;
    varargout{2} = vertices_reduced;
  else # fv structure
    varargout{1}.faces = faces_reduced;
    varargout{1}.vertices = vertices_reduced;
  endif

endfunction

function [faces, vertices, max_faces, patch_handle, fast, verbose] = __get_check_reducepatch_args__ (varargin)
## assign input parameters and check their validity

  ## default values
  faces = vertices = max_faces = patch_handle = [];
  reduction_factor = 0.5;
  fast = verbose = false;

  num_string_inputs = 0;
  for i_arg = (nargin:-1:nargin-1)
    ## check whether last maximum 2 input arguments are strings and assign parameters
    if (!ischar (varargin{i_arg}) || i_arg < 1)
      ## string arguments must be at the end
      break
    end
    switch (tolower (varargin{i_arg}))
      case {"f", "fast"}
        fast = true;
        num_string_inputs++;
      case {"v", "verbose"}
        verbose = true;
        num_string_inputs++;
      otherwise
        error ("reducepatch: parameter '%s' not supported", varargin{i_arg})
    endswitch
  endfor

  ## get faces and vertices
  i_fv = 1;
  if (isstruct (varargin{1}))
    if (all(isfield (varargin{1}, {"vertices", "faces"})))
      vertices = varargin{1}.vertices;
      faces = varargin{1}.faces;
    else
      error ("reducepatch: The struct FV must contain the fields 'vertices' and 'faces'.")
    endif
  elseif (isscalar (varargin{1}))
    patch_handle = varargin{1};
    if (ishandle (patch_handle) && strcmpi (get (patch_handle, "Type"), "patch"))
      vertices = get (patch_handle, "Vertices");
      faces = get (patch_handle, "Faces");
    else
      error ("reducepatch: PATCH_HANDLE must be a valid handle to a patch")
    endif
  elseif (ismatrix (varargin{1}))
    faces = varargin{1};
    if (nargin - num_string_inputs > 1 && ismatrix (varargin{2}))
      vertices = varargin{2};
      i_fv = 2;
    else
      error ("reducepatch: If first argument is a matrix containing FACES, second argument must be a matrix containing VERTICES")
    endif
  else
    print_usage ();
  endif

  ## get reduction_factor
  if (nargin - num_string_inputs > i_fv) 
    if (isscalar (varargin{i_fv+1}) && varargin{i_fv+1} > 0)
      reduction_factor = varargin{i_fv+1};
    else
      error ("reducepatch: REDUCTION_FACTOR must be a scalar > 0")
    endif
  endif

  ## check faces and vertices
  if (size (vertices, 2) != 3)
    error ("reducepatch: VERTICES must be a Mx3 matrix")
  endif
  if (size (faces, 2) != 3)
    ## FIXME: Convert faces to only triangles if necessary
    error ("reducepatch: Currently patches must consist of triangles only")
  endif
  if (any (mod (faces(:), 1)) || any (faces(:) < 1))
    error ("reducepatch: FACES must consist of positiv integer indices only")
  endif
  if (max (faces(:)) > size (vertices, 1))
    error ("reducepatch: not enough VERTICES for given FACES")
  endif

  ## get max_faces
  num_faces = size (faces, 1);
  if (reduction_factor <= 1)
    max_faces = floor (num_faces * reduction_factor);
  else
    max_faces = reduction_factor;
  endif

endfunction

function [faces_reduced, vertices_reduced] = __reducepatch__ (vertices, faces, max_faces, verbose)
## actual reduction of faces

  num_faces = size (faces, 1);
  faces = sort (faces, 2);

  ## get all unique edges
  all_edges = sort([faces(:,1) faces(:,2); faces(:,2) faces(:,3); faces(:,3) faces(:,1)], 2);
  edges = unique (all_edges, "rows");
  num_edges = size (edges, 1);

  ## calculate length of edges
  edge_length = norm (vertices(edges(:,1),:) - vertices(edges(:,2),:), 2, "rows");

  if (verbose)
    fprintf ("reducepatch: reducing number of faces")
    ## do not spam the command window with dots if many faces are collapsed
    verbose_stepwidth = (num_faces - max_faces) / 100; # max. 100 dots should be enough
    verbose_counter = 0;
  endif

  ## reduce number of faces
  clean_finish = true;
  while (max_faces < num_faces - sum (isinf (faces(:,1))))
    if (verbose)
      verbose_counter++;
      if (verbose_counter > verbose_stepwidth)
        fprintf (".")
        verbose_counter = 0;
      endif
    endif

    ## find shortest edge
    i_shortest_edge = find (edge_length == min (edge_length), 1, "first");
    if ( any ( isinf ( vertices(edges(i_shortest_edge,:),:) ) ) )
      warning ("reducepatch: shortest edge is marked as deleted. This should never happen.")
      clean_finish = false;
      break
    endif

    start_point_idx = edges(i_shortest_edge,1);
    end_point_idx = edges(i_shortest_edge,2);

    ## collapse the edge to its midpoint
    vertices(start_point_idx,:) = 0.5* ( vertices(start_point_idx,:) + vertices(end_point_idx,:) );
    # endpoint will be deleted later

    ## remove collapsed edge info
    edges(i_shortest_edge,:) = [];
    edge_length(i_shortest_edge) = [];

    ## unite collapsed edges of neighboring faces to one edge
    faces_with_collapsed_edge = any (faces == start_point_idx, 2) & any (faces == end_point_idx, 2);
    edges_of_collapsed_faces = faces(faces_with_collapsed_edge,:);
    ## get other vertex of collapsed faces
    third_vertices = edges_of_collapsed_faces(faces(faces_with_collapsed_edge,:) != start_point_idx & faces(faces_with_collapsed_edge,:) != end_point_idx);
    if (!isempty (third_vertices))
      ## delete edge of the collapsed faces which also has end_point (keep the one with start_point)
      edges_to_delete = any (edges == end_point_idx, 2) & any (any (bsxfun (@eq, edges, reshape(third_vertices, 1, 1, [])), 3), 2);
      edges(edges_to_delete,:) = [];
      edge_length(edges_to_delete) = [];
    endif

    ## mark the faces that collapsed for later removal
    faces(faces_with_collapsed_edge,:) = Inf;

    ## update the lengths of the moved edges
    edges_with_moved_point = any (edges == start_point_idx, 2);
    edge_length(edges_with_moved_point) = norm (vertices(edges(edges_with_moved_point,1),:) - vertices(edges(edges_with_moved_point,2),:), 2, "rows");

    ## replace all vertices that use end_point to use start_point
    if (start_point_idx != end_point_idx)
      edges(edges == end_point_idx) = start_point_idx;
      faces(faces == end_point_idx) = start_point_idx;
      vertices(end_point_idx,:) = Inf; # mark vertex for later removal
    endif

  endwhile

  if (verbose)
    fprintf("\n")
  endif

  if (!clean_finish)
    ## continue anyway?
    ## FIXME: What should be done in this case?
  endif

  ## finally, remove vertices and faces
  removed_vertices = isinf (vertices(:,1));
  vertices_reduced = vertices(!removed_vertices,:);

  faces_reduced = faces(!isinf (faces(:,1)),:);

  ## renumber vertices in faces with subsequent numbers
  vertex_renum = ones (size (vertices, 1), 1);
  vertex_renum(removed_vertices) = 0;
  vertex_renum = cumsum (vertex_renum);
  faces_reduced = vertex_renum(faces_reduced);

endfunction


%!demo
%! [x,y,z] = meshgrid (-2:0.5:2, -2:0.5:2, -2:0.5:2);
%! val = x.^2 + y.^2 + z.^2;
%! fv = isosurface (x, y, z, val, 1);
%! figure;
%! subplot (1, 2, 1)
%! patch (fv, "FaceColor", "g"), view (3), axis equal
%! title ("Sphere with all faces")
%! subplot(1, 2, 2)
%! patch (reducepatch (fv), "FaceColor", "g"), view (3), axis equal
%! title ("Sphere with reduced number of faces")

%!shared faces, vertices, fv, num_faces
%! [x,y,z] = meshgrid (-2:0.5:2, -2:0.5:2, -2:0.5:2);
%! val = x.^2 + y.^2 + z.^2;
%! [faces, vertices] = isosurface (x, y, z, val, 1);
%! fv.faces = faces; fv.vertices = vertices;
%! num_faces = size (faces, 1);

## one input (structure), one output
%!test
%! fv_reduced = reducepatch (fv);
%! assert (size (fv_reduced.faces, 1), num_faces * .5, 3)
%! assert (size (fv_reduced.faces, 2), 3)
%! assert (size (fv_reduced.vertices, 2), 3)

## two inputs (faces, vertices), one output
%!test
%! fv_reduced = reducepatch (faces, vertices);
%! assert (size (fv_reduced.faces, 1), num_faces * .5, 3)
%! assert (size (fv_reduced.faces, 2), 3)
%! assert (size (fv_reduced.vertices, 2), 3)

## two inputs (structure, reduction_factor > 1), two outputs
%!test
%! [faces_reduced, vertices_reduced] = reducepatch (fv, 20);
%! assert (size (faces_reduced, 1), 20, 3)
%! assert (size (faces_reduced, 2), 3)
%! assert (size (vertices_reduced, 2), 3)

## three inputs (faces, vertices, reduction_factor < 1), two outputs
%!test
%! [faces_reduced, vertices_reduced] = reducepatch (faces, vertices, .5);
%! assert (size (faces_reduced, 1), num_faces * .5, 3)
%! assert (size (faces_reduced, 2), 3)
%! assert (size (vertices_reduced, 2), 3)

## two inputs (structure, reduction_factor < 1), no outputs (update patch object in figure)
%!test
%! h_fig = figure ("Visible", "off");
%! p = patch (fv);
%! reducepatch (p, .35, "v");
%! assert (size (get (p, "Faces"), 1), num_faces * .35, 3)
%! assert (size (get (p, "Faces"), 2), 3)
%! assert (size (get (p, "Vertices"), 2), 3)
%! close (h_fig)

## three inputs (faces, vertices, reduction_factor > 1) + one string, one output
%!test
%! fv_reduced = reducepatch (faces, vertices, 52, "v");
%! assert (size (fv_reduced.faces, 1), 52, 3)
%! assert (size (fv_reduced.faces, 2), 3)
%! assert (size (fv_reduced.vertices, 2), 3)

## two inputs (structure, reduction_factor < 1) + one string, two outputs
%!test
%! [faces_reduced, vertices_reduced] = reducepatch (fv, .4, "fast");
%! assert (size (faces_reduced, 1), num_faces * .4, 3)
%! assert (size (faces_reduced, 2), 3)
%! assert (size (vertices_reduced, 2), 3)

## one input (structure) + two strings, two outputs
%!test
%! [faces_reduced, vertices_reduced] = reducepatch (fv, "f", "verbose");
%! assert (size (faces_reduced, 1), num_faces * .5, 3)
%! assert (size (faces_reduced, 2), 3)
%! assert (size (vertices_reduced, 2), 3)

## test for each error
%!test
%!error <Invalid call to reducepatch.> reducepatch ();
%!error <Invalid call to reducepatch.> reducepatch (fv, faces, vertices, .5, "f", "v");
%!error <Invalid call to reducepatch.> [f, v, c] = reducepatch (fv);
%!error <reducepatch: parameter 'test' not supported> fv_reduced = reducepatch (faces, vertices, .7, "test");
%!error <reducepatch: The struct FV must contain the fields 'vertices' and 'faces'.> fv_incomplete.faces = faces; reducepatch (fv_incomplete, .7);
%!error <reducepatch: PATCH_HANDLE must be a valid handle to a patch> reducepatch (15, .7);
%!error <reducepatch: If first argument is a matrix containing FACES, second argument must be a matrix containing VERTICES> reducepatch (faces);
%!error <reducepatch: REDUCTION_FACTOR must be a scalar > 0> reducepatch (fv, [.7 .5]);
%!error <reducepatch: REDUCTION_FACTOR must be a scalar > 0> reducepatch (fv, -5);
%!error <reducepatch: VERTICES must be a Mx3 matrix> reducepatch (faces, .7, "v");
%!error <reducepatch: Currently patches must consist of triangles only> faces_new = NaN (size (faces, 1), 4); faces_new(:,1:3) = faces; faces_new(1,4) = 5; reducepatch (faces_new, vertices);
%!error <reducepatch: not enough VERTICES for given FACES> faces_new = faces; faces_new(1,3) = size (vertices, 1) + 1; reducepatch (faces_new, vertices);
