octave_packages/m/geometry/delaunayn.m

   1 ## Copyright (C) 2007-2012 David Bateman
   2 ##
   3 ## This file is part of Octave.
   4 ##
   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.
   9 ##
  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.
  14 ##
  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/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{T} =} delaunayn (@var{pts})
  21 ## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{pts}, @var{options})
  22 ## Compute the Delaunay triangulation for an N-dimensional set of points.
  23 ## The Delaunay triangulation is a tessellation of the convex hull of a set
  24 ## of points such that no N-sphere defined by the N-triangles contains
  25 ## any other points from the set.
  26 ##
  27 ## The input matrix @var{pts} of size [n, dim] contains n points in a space of
  28 ## dimension dim.  The return matrix @var{T} has size [m, dim+1].  Each row
  29 ## of @var{T} contains a set of indices back into the original set of points
  30 ## @var{pts} which describes a simplex of dimension dim.  For example, a 2-D
  31 ## simplex is a triangle and 3-D simplex is a tetrahedron.
  32 ##
  33 ## An optional second argument, which must be a string or cell array of strings,
  34 ## contains options passed to the underlying qhull command.
  35 ## See the documentation for the Qhull library for details
  36 ## @url{http://www.qhull.org/html/qh-quick.htm#options}.
  37 ## The default options depend on the dimension of the input:
  38 ##
  39 ## @itemize
  40 ## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}}
  41 ##
  42 ## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}}
  43 ## @end itemize
  44 ##
  45 ## If @var{options} is not present or @code{[]} then the default arguments are
  46 ## used.  Otherwise, @var{options} replaces the default argument list.
  47 ## To append user options to the defaults it is necessary to repeat the
  48 ## default arguments in @var{options}.  Use a null string to pass no arguments.
  49 ##
  50 ## @seealso{delaunay, delaunay3, convhulln, voronoin}
  51 ## @end deftypefn
  52
  53 function T = delaunayn (pts, varargin)
  54
  55   if (nargin < 1)
  56     print_usage ();
  57   endif
  58
  59   T = __delaunayn__ (pts, varargin{:});
  60
  61   if (isa (pts, "single"))
  62     myeps = eps ("single");
  63   else
  64     myeps = eps;
  65   endif
  66
  67   ## Try to remove the zero volume simplices.  The volume of the i-th simplex is
  68   ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n)
  69   ## (reference http://en.wikipedia.org/wiki/Simplex).  Any simplex with a
  70   ## relative volume less than some arbitrary criteria is rejected.  The
  71   ## criteria we use is the volume of the simplex corresponding to an
  72   ## orthogonal simplex is equal edge length all equal to the edge length of
  73   ## the original simplex.  If the relative volume is 1e3*eps then the simplex
  74   ## is rejected.  Note division of the two volumes means that the factor
  75   ## prod(1:n) is dropped.
  76   idx = [];
  77   [nt, n] = size (T);
  78   ## FIXME: Vectorize this for loop or convert to delaunayn to .oct function
  79   for i = 1:nt
  80     X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
  81     if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
  82       idx(end+1) = i;
  83     endif
  84   endfor
  85   T(idx,:) = [];
  86
  87 endfunction
  88
  89
  90 %% FIXME: Need tests for delaunayn
  91
  92 %% FIXME: Need input validation tests
  93