########################################################################
##
## Copyright (C) 2007-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## 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
## .
##
########################################################################
## -*- 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 with (ndim-1) edge lengths equal to the edge lengths of
## the original simplex. If the relative volume is smaller than 1e3*eps, the
## simplex is rejected. Note division of the two volumes means that the
## factor factorial(ndim+1) is dropped from the volume calculations.
[nt, nd] = size (T); ## nt = simplex count, nd = # of simplex points
dim = nd - 1;
## calculate common origin edge vectors for each simplex (p2-p1, p3-p1, ...)
## stored in [nt x (nd-1), nd-1] array, grouped by simplex.
edg_vec = pts(T(:,2:nd).'(:),:) - kron (pts(T(:,1),:), ones (dim, 1));
## Calculate simplex volumes according to dimensionality of problem
if (any(dim == [2 3]))%[2, 3]
## faster vector-product code paths for 2-D & 3-D cases
## reshape edge vectors to 3D array such that
## rows = nt tris, cols = dim coordinates, pages = dim vectors
edg_vec = permute(reshape(edg_vec, [dim nt dim]), [2 3 1]);
if (dim == 2)
#2D Use simple component cross product to calculate 2D triangle area
vol = abs (edg_vec(:,:,1)(:,1) .* edg_vec(:,:,2)(:,2)...
- edg_vec(:,:,2)(:,1) .* edg_vec(:,:,1)(:,2));
else
#3D: Use scalar triple product to calculate 3D tetrahedron volume
vol = abs (dot (edg_vec(:,:,1), cross (edg_vec(:,:,2), edg_vec(:,:,3), 2), 2));
endif
else
## 4D and higher dimensions
## Place edge vectors into a block diagonal matrix
eqs = sparse (dim * nt, dim * nt);
eqs(logical (kron (speye (nt, nt), true (dim)))) = edg_vec.'(:);
## Extract diagonal of LU factorization of that block diagonal matrix
[~, u, p, ~] = lu (eqs, "vector");
reordered_tri_idx = kron (1:nt, ones (1, dim))(p);
R = abs (diag (u));
## extract simplex volumes as product of diagonal elements of u
## preserving order relative to simplexes from delaunay triangulation
##TODO: figure out why
## kron (R, ones(1,nt)) and ()(reordered_tri_idx.'==(1:nt)) parts
## this gets 100x slower for larger problems
## if can't redo/work around will need to work with R instead of vol
## and make a separate delaunayn ()