1## Copyright (C) 2007, 2008, 2009 David Bateman
3## This file is part of Octave.
5## Octave is free software; you can redistribute it and/or modify it
6## under the terms of the GNU General Public License as published by
7## the Free Software Foundation; either version 3 of the License, or (at
8## your option) any later version.
10## Octave is distributed in the hope that it will be useful, but
11## WITHOUT ANY WARRANTY; without even the implied warranty of
12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13## General Public License for more details.
15## You should have received a copy of the GNU General Public License
16## along with Octave; see the file COPYING. If not, see
17## <http://www.gnu.org/licenses/>.
20## @deftypefn {Function File} {@var{T} =} delaunayn (@var{P})
21## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{P}, @var{opt})
22## Form the Delaunay triangulation for a set of points.
23## The Delaunay triangulation is a tessellation of the convex hull of the
24## points such that no n-sphere defined by the n-triangles contains
25## any other points from the set.
26## The input matrix @var{P} of size @code{[n, dim]} contains @var{n}
27## points in a space of dimension dim. The return matrix @var{T} has the
28## size @code{[m, dim+1]}. It contains for each row a set of indices to
29## the points, which describes a simplex of dimension dim. For example,
30## a 2d simplex is a triangle and 3d simplex is a tetrahedron.
32## Extra options for the underlying Qhull command can be specified by the
33## second argument. This argument is a cell array of strings. The default
34## options depend on the dimension of the input:
37## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}}
38## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}}
41## If @var{opt} is [], then the default arguments are used. If @var{opt}
42## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull.
43## See the Qhull documentation for the available options.
45## All options can also be specified as single string, for example
46## @code{"Qt Qbb Qc Qz"}.
50function t = delaunayn (x, varargin)
55 t = __delaunayn__ (x, varargin{:});
57 if (isa (x, "single"))
58 myeps = eps ("single");
63 ## Try to remove the zero volume simplices. The volume of the i-th simplex is
64 ## given by abs(det(x(t(i,1:end-1),:)-x(t(i,2:end),:)))/prod(1:n)
65 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a
66 ## relative volume less than some arbitrary criteria is rejected. The
67 ## criteria we use is the volume of the simplex corresponding to an
68 ## orthogonal simplex is equal edge length all equal to the edge length of
69 ## the original simplex. If the relative volume is 1e3*eps then the simplex
70 ## is rejected. Note division of the two volumes means that the factor
71 ## prod(1:n) is dropped.
75 X = x(t(i,1:end-1),:) - x(t(i,2:end),:);
76 if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
