######################################################################## ## ## 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])) ## 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 = 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 ()