Mercurial > octave

comparison scripts/geometry/delaunayn.m @ 13746:7ff0bdc3dc4c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Revamp geometry functions dependent on Qhull (Bug #34604, Bug #33346) * NEWS : Document new options being passed to Qhull * convhull.m, delaunay.m, delaunay3.m, delaunayn.m, voronoi.m, voronoin.m: Update docstrings. Put input validation first. Use same variable names as Matlab. Restore random state altered in demos. * __delaunayn__.cc: Use common syntax for parsing OPTIONS input. Add 'Qz' option to qhull command for 2D,3D data. Correctly free all Qhull memory and avoid segfault with non-simplicial facets. * __voronoi__.cc: Use common syntax for parsing OPTIONS input. Correctly free all Qhull memory. * convhulln.cc: Use common syntax for parsing OPTIONS input. Use Matlab-compatible options for qhull command. Correctly free all Qhull memory. Allow return of non-simplicial facets without causing a segfault.
author Rik <octave@nomad.inbox5.com>
date Tue, 25 Oct 2011 10:17:23 -0700
parents c792872f8942
children 72c96de7a403
comparison
equal deleted inserted replaced
13745:3c3b74677fa0 13746:7ff0bdc3dc4c
15 ## You should have received a copy of the GNU General Public License 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {@var{t} =} delaunayn (@var{p}) 20 ## @deftypefn {Function File} {@var{T} =} delaunayn (@var{pts})
21 ## @deftypefnx {Function File} {@var{t} =} delaunayn (@var{p}, @var{opt}) 21 ## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{pts}, @var{options})
22 ## Form the Delaunay triangulation for a set of points. 22 ## Compute the Delaunay triangulation for an N-dimensional set of points.
23 ## The Delaunay triangulation is a tessellation of the convex hull of the 23 ## The Delaunay triangulation is a tessellation of the convex hull of a set
24 ## points such that no n-sphere defined by the n-triangles contains 24 ## of points such that no N-sphere defined by the N-triangles contains
25 ## any other points from the set. 25 ## any other points from the set.
26 ## The input matrix @var{p} of size @code{[n, dim]} contains @math{n} 26 ##
27 ## points in a space of dimension dim. The return matrix @var{t} has the 27 ## The input matrix @var{pts} of size [n, dim] contains n points in a space of
28 ## size @code{[m, dim+1]}. It contains for each row a set of indices to 28 ## dimension dim. The return matrix @var{T} has size [m, dim+1]. Each row
29 ## the points, which describes a simplex of dimension dim. For example, 29 ## of @var{T} contains a set of indices back into the original set of points
30 ## a 2-D simplex is a triangle and 3-D simplex is a tetrahedron. 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.
31 ## 32 ##
32 ## Extra options for the underlying Qhull command can be specified by the 33 ## An optional second argument, which must be a string or cell array of strings,
33 ## second argument. This argument is a cell array of strings. The default 34 ## contains options passed to the underlying qhull command.
34 ## options depend on the dimension of the input: 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:
35 ## 38 ##
36 ## @itemize 39 ## @itemize
37 ## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}} 40 ## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}}
38 ## 41 ##
39 ## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} 42 ## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}}
40 ## @end itemize 43 ## @end itemize
41 ## 44 ##
42 ## If @var{opt} is [], then the default arguments are used. If @var{opt} 45 ## If @var{options} is not present or @code{[]} then the default arguments are
43 ## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull. 46 ## used. Otherwise, @var{options} replaces the default argument list.
44 ## See the Qhull documentation for the available options. 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.
45 ## 49 ##
46 ## All options can also be specified as single string, for example 50 ## @seealso{delaunay, delaunay3, convhulln, voronoin}
47 ## @code{"Qt Qbb Qc Qz"}.
48 ##
49 ## @end deftypefn 51 ## @end deftypefn
50 52
51 function t = delaunayn (p, varargin) 53 function T = delaunayn (pts, varargin)
54
52 if (nargin < 1) 55 if (nargin < 1)
53 print_usage (); 56 print_usage ();
54 endif 57 endif
55 58
56 t = __delaunayn__ (p, varargin{:}); 59 T = __delaunayn__ (pts, varargin{:});
57 60
58 if (isa (p, "single")) 61 if (isa (pts, "single"))
59 myeps = eps ("single"); 62 myeps = eps ("single");
60 else 63 else
61 myeps = eps; 64 myeps = eps;
62 endif 65 endif
63 66
64 ## Try to remove the zero volume simplices. The volume of the i-th simplex is 67 ## Try to remove the zero volume simplices. The volume of the i-th simplex is
65 ## given by abs(det(p(t(i,1:end-1),:)-p(t(i,2:end),:)))/prod(1:n) 68 ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n)
66 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a 69 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a
67 ## relative volume less than some arbitrary criteria is rejected. The 70 ## relative volume less than some arbitrary criteria is rejected. The
68 ## criteria we use is the volume of the simplex corresponding to an 71 ## criteria we use is the volume of the simplex corresponding to an
69 ## orthogonal simplex is equal edge length all equal to the edge length of 72 ## orthogonal simplex is equal edge length all equal to the edge length of
70 ## the original simplex. If the relative volume is 1e3*eps then the simplex 73 ## the original simplex. If the relative volume is 1e3*eps then the simplex
71 ## is rejected. Note division of the two volumes means that the factor 74 ## is rejected. Note division of the two volumes means that the factor
72 ## prod(1:n) is dropped. 75 ## prod(1:n) is dropped.
73 idx = []; 76 idx = [];
74 [nt, n] = size (t); 77 [nt, n] = size (T);
78 ## FIXME: Vectorize this for loop or convert to delaunayn to .oct function
75 for i = 1:nt 79 for i = 1:nt
76 X = p(t(i,1:end-1),:) - p(t(i,2:end),:); 80 X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
77 if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) 81 if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
78 idx = [idx, i]; 82 idx(end+1) = i;
79 endif 83 endif
80 endfor 84 endfor
81 t(idx,:) = []; 85 T(idx,:) = [];
86
82 endfunction 87 endfunction
88
89
90 %% FIXME: Need tests for delaunayn
91
92 %% FIXME: Need input validation tests
93