########################################################################
##
## Copyright (C) 2007-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## 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
## <https://www.gnu.org/licenses/>.
##
########################################################################

## -*- texinfo -*-
## @deftypefn  {} {@var{T} =} delaunayn (@var{pts})
## @deftypefnx {} {@var{T} =} delaunayn (@var{pts}, @var{options})
## Compute the Delaunay triangulation for an N-dimensional set of points.
##
## The Delaunay triangulation is a tessellation of the convex hull of a set of
## points such that no N-sphere defined by the N-triangles contains any other
## points from the set.
##
## The input matrix @var{pts} of size [n, dim] contains n points in a space of
## dimension dim.  The return matrix @var{T} has size [m, dim+1].  Each row of
## @var{T} contains a set of indices back into the original set of points
## @var{pts} which describes a simplex of dimension dim.  For example, a 2-D
## simplex is a triangle and 3-D simplex is a tetrahedron.
##
## An optional second argument, which must be a string or cell array of
## strings, contains options passed to the underlying qhull command.  See the
## documentation for the Qhull library for details
## @url{http://www.qhull.org/html/qh-quick.htm#options}.
## The default options depend on the dimension of the input:
##
## @itemize
## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc"@}}
##
## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}}
## @end itemize
##
## If Qhull fails for 2-D input the triangulation is attempted again with
## the options @code{@{"Qt", "Qbb", "Qc", "Qz"@}} which may result in
## reduced accuracy.
##
## If @var{options} is not present or @code{[]} then the default arguments are
## used.  Otherwise, @var{options} replaces the default argument list.
## To append user options to the defaults it is necessary to repeat the
## default arguments in @var{options}.  Use a null string to pass no arguments.
##
## @seealso{delaunay, convhulln, voronoin, trimesh, tetramesh}
## @end deftypefn

function T = delaunayn (pts, varargin)

  if (nargin < 1)
    print_usage ();
  endif

  if (isempty (varargin) || isempty (varargin{1}))
    try
      T = __delaunayn__ (pts);
    catch
      if (columns (pts) <= 2)
        T = __delaunayn__ (pts, "Qt Qbb Qc Qz");
      endif
    end_try_catch
  else
    T = __delaunayn__ (pts, varargin{:});
  endif

  if (isa (pts, "single"))
    tol = 1e3 * eps ("single");
  else
    tol = 1e3 * eps;
  endif

  ## Try to remove the zero volume simplices.  The volume of the i-th simplex is
  ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/factorial(ndim+1)
  ## (reference http://en.wikipedia.org/wiki/Simplex).  Any simplex with a
  ## relative volume less than some arbitrary criteria is rejected.  The
  ## criteria we use is the volume of the simplex corresponding to an
  ## orthogonal simplex is equal edge length all equal to the edge length of
  ## the original simplex.  If the relative volume is 1e3*eps then the simplex
  ## is rejected.  Note division of the two volumes means that the factor
  ## factorial(ndim+1) is dropped.
  [nt, nd] = size (T);
  if (nd == 3)
    ## 2-D case
    np = rows (pts);
    ptsz = [pts, zeros(np, 1)];
    p1 = ptsz(T(:,1), :);
    p2 = ptsz(T(:,2), :);
    p3 = ptsz(T(:,3), :);
    p12 = p1 - p2;
    p23 = p2 - p3;
    det = cross (p12, p23, 2);
    idx = abs (det (:,3) ./ sqrt (sumsq (p12, 2))) < tol & ...
          abs (det (:,3) ./ sqrt (sumsq (p23, 2))) < tol;
  else
    edgvec = pts(T(:, 2:nd).'(:), :) - kron (pts(T(:, 1), :),ones (nd - 1, 1));
    eqs = sparse ((nd - 1) * nt,(nd - 1) * nt);
    eqs(logical (kron (speye (nt, nt), true (nd-1))))=edgvec.'(:);
    [l u p q] = lu (eqs, "vector");
    R = abs (diag (u));
    reorderdtriidx = kron (1:nt, ones (1,nd-1))(p)(q);
    idx = unique (reorderdtriidx(R < 100 * (nd - 1) * nt * eps (max (R))));
  endif
  T(idx,:) = [];

endfunction


%!testif HAVE_QHULL
%! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0];
%! assert (sortrows (sort (delaunayn (x), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]);

## Test 3-D input
%!testif HAVE_QHULL
%! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
%! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)), [1,2,3,4;1,2,4,5]);

## FIXME: Need tests for delaunayn

## Input validation tests
%!error <Invalid call> delaunayn ()