Mercurial > octave
annotate scripts/geometry/delaunayn.m @ 7795:df9519e9990c
Handle single precision eps values
author | David Bateman <dbateman@free.fr> |
---|---|
date | Mon, 12 May 2008 22:57:11 +0200 |
parents | a730e47fda4d |
children | eb63fbe60fab |
rev | line source |
---|---|
7017 | 1 ## Copyright (C) 2007 David Bateman |
6823 | 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 | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
6823 | 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 | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
6823 | 18 |
19 ## -*- texinfo -*- | |
6846 | 20 ## @deftypefn {Function File} {@var{T} =} delaunayn (@var{P}) |
21 ## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{P}, @var{opt}) | |
6823 | 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 | |
6826 | 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. | |
6823 | 31 ## |
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: | |
35 ## | |
36 ## @itemize | |
6826 | 37 ## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}} |
38 ## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} | |
39 ## @end itemize | |
6823 | 40 ## |
41 ## If @var{opt} is [], then the default arguments are used. If @var{opt} | |
6826 | 42 ## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull. |
6823 | 43 ## See the Qhull documentation for the available options. |
44 ## | |
45 ## All options can also be specified as single string, for example | |
6826 | 46 ## @code{"Qt Qbb Qc Qz"}. |
6823 | 47 ## |
48 ## @end deftypefn | |
49 | |
50 function t = delaunayn (x, varargin) | |
51 if (nargin < 1) | |
6826 | 52 print_usage (); |
6823 | 53 endif |
54 | |
7208 | 55 t = __delaunayn__ (x, varargin{:}); |
6823 | 56 |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
57 if (isa (x, "single")) |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
58 myeps = eps ("single"); |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
59 else |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
60 myeps = eps; |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
61 endif |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
62 |
6823 | 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. | |
72 idx = []; | |
6826 | 73 [nt, n] = size (t); |
74 for i = 1:nt | |
6823 | 75 X = x(t(i,1:end-1),:) - x(t(i,2:end),:); |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7208
diff
changeset
|
76 if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) |
6823 | 77 idx = [idx, i]; |
78 endif | |
79 endfor | |
80 t(idx,:) = []; | |
81 endfunction |