Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/__voronoi__.cc @ 9003:0631d397fbe0
replace lo_ieee_isnan by xisnan, add missing includes
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Sat, 21 Mar 2009 20:46:21 +0100 |
parents | eb63fbe60fab |
children | 2c279308f6ab |
rev | line source |
---|---|
6823 | 1 /* |
2 | |
8920 | 3 Copyright (C) 2000, 2007, 2008 Kai Habel |
6823 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
6823 | 11 |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
6823 | 20 |
21 */ | |
22 | |
23 /* | |
24 20. Augiust 2000 - Kai Habel: first release | |
25 */ | |
26 | |
27 /* | |
28 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net> | |
29 Added optional second argument to pass options to the underlying | |
30 qhull command | |
31 */ | |
32 | |
33 #include <iostream> | |
34 #ifdef HAVE_CONFIG_H | |
35 #include <config.h> | |
36 #endif | |
9003
0631d397fbe0
replace lo_ieee_isnan by xisnan, add missing includes
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
37 #include "lo-ieee.h" |
6823 | 38 #include "oct.h" |
39 #include "Cell.h" | |
40 | |
41 #ifdef HAVE_QHULL | |
42 extern "C" { | |
43 #include <qhull/qhull_a.h> | |
44 } | |
45 | |
46 #ifdef NEED_QHULL_VERSION | |
47 char qh_version[] = "__voronoi__.oct 2007-07-24"; | |
48 #endif | |
49 #endif | |
50 | |
51 DEFUN_DLD (__voronoi__, args, , | |
52 "-*- texinfo -*-\n\ | |
53 @deftypefn {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\ | |
54 @deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\ | |
6945 | 55 Undocumented internal function.\n\ |
6823 | 56 @end deftypefn") |
57 { | |
58 octave_value_list retval; | |
6945 | 59 |
6823 | 60 #ifdef HAVE_QHULL |
6945 | 61 |
6823 | 62 retval(0) = 0.0; |
63 | |
6945 | 64 int nargin = args.length (); |
6823 | 65 if (nargin < 1 || nargin > 2) |
66 { | |
67 print_usage (); | |
68 return retval; | |
69 } | |
70 | |
6880 | 71 const char *options; |
72 | |
6823 | 73 if (nargin == 2) |
74 { | |
6945 | 75 if (! args (1).is_string ()) |
6823 | 76 { |
77 error ("__voronoi__: second argument must be a string"); | |
78 return retval; | |
79 } | |
6945 | 80 |
81 options = args (1).string_value().c_str (); | |
6823 | 82 } |
83 else | |
84 options = ""; | |
85 | |
6945 | 86 Matrix p (args(0).matrix_value ()); |
6823 | 87 |
6945 | 88 const octave_idx_type dim = p.columns (); |
89 const octave_idx_type np = p.rows (); | |
90 p = p.transpose (); | |
6823 | 91 |
6945 | 92 double *pt_array = p.fortran_vec (); |
93 | |
6823 | 94 //double pt_array[dim * np]; |
95 //for (int i = 0; i < np; i++) | |
96 // { | |
97 // for (int j = 0; j < dim; j++) | |
98 // { | |
99 // pt_array[j+i*dim] = p(i,j); | |
100 // } | |
101 // } | |
102 | |
103 boolT ismalloc = false; | |
104 | |
6945 | 105 OCTAVE_LOCAL_BUFFER (char, flags, 250); |
6880 | 106 |
6823 | 107 // hmm lot's of options for qhull here |
6945 | 108 sprintf (flags, "qhull v Fv T0 %s", options); |
6880 | 109 |
6945 | 110 // If you want some debugging information replace the 0 pointer |
111 // with stdout or some other file open for writing. | |
112 | |
113 FILE *outfile = 0; | |
6880 | 114 FILE *errfile = stderr; |
115 | |
6945 | 116 if (! qh_new_qhull (dim, np, pt_array, ismalloc, flags, outfile, errfile)) |
6823 | 117 { |
118 facetT *facet; | |
119 vertexT *vertex; | |
120 | |
121 octave_idx_type i = 0, n = 1, k = 0, m = 0, fidx = 0, j = 0, r = 0; | |
6945 | 122 OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, np); |
6823 | 123 |
124 for (i = 0; i < np; i++) | |
125 ni[i] = 0; | |
6945 | 126 qh_setvoronoi_all (); |
6823 | 127 bool infinity_seen = false; |
6945 | 128 facetT *neighbor, **neighborp; |
6823 | 129 coordT *voronoi_vertex; |
6945 | 130 |
6823 | 131 FORALLfacets |
132 { | |
133 facet->seen = false; | |
134 } | |
6945 | 135 |
6823 | 136 FORALLvertices |
137 { | |
138 if (qh hull_dim == 3) | |
6945 | 139 qh_order_vertexneighbors (vertex); |
6823 | 140 infinity_seen = false; |
6945 | 141 |
142 FOREACHneighbor_ (vertex) | |
6823 | 143 { |
6945 | 144 if (! neighbor->upperdelaunay) |
6823 | 145 { |
6945 | 146 if (! neighbor->seen) |
6823 | 147 { |
148 n++; | |
6945 | 149 neighbor->seen = true; |
6823 | 150 } |
151 ni[k]++; | |
152 } | |
6945 | 153 else if (! infinity_seen) |
6823 | 154 { |
155 infinity_seen = true; | |
156 ni[k]++; | |
157 } | |
158 } | |
159 k++; | |
160 } | |
161 | |
6945 | 162 Matrix v (n, dim); |
6823 | 163 for (octave_idx_type d = 0; d < dim; d++) |
164 v(0,d) = octave_Inf; | |
165 | |
6945 | 166 boolMatrix AtInf (np, 1); |
6823 | 167 for (i = 0; i < np; i++) |
168 AtInf(i) = false; | |
8415
fa12c67683d3
Reszie third arg of __voronoi__ before retuning it
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
169 octave_value_list F (np, octave_value ()); |
6823 | 170 k = 0; |
171 i = 0; | |
6945 | 172 |
6823 | 173 FORALLfacets |
174 { | |
175 facet->seen = false; | |
176 } | |
6945 | 177 |
6823 | 178 FORALLvertices |
179 { | |
180 if (qh hull_dim == 3) | |
181 qh_order_vertexneighbors(vertex); | |
182 infinity_seen = false; | |
6945 | 183 RowVector facet_list (ni[k++]); |
6823 | 184 m = 0; |
6945 | 185 |
6823 | 186 FOREACHneighbor_(vertex) |
187 { | |
188 if (neighbor->upperdelaunay) | |
189 { | |
6945 | 190 if (! infinity_seen) |
6823 | 191 { |
192 infinity_seen = true; | |
193 facet_list(m++) = 1; | |
194 AtInf(j) = true; | |
195 } | |
196 } | |
197 else | |
198 { | |
6945 | 199 if (! neighbor->seen) |
6823 | 200 { |
201 voronoi_vertex = neighbor->center; | |
202 fidx = neighbor->id; | |
203 i++; | |
204 for (octave_idx_type d = 0; d < dim; d++) | |
205 { | |
206 v(i,d) = *(voronoi_vertex++); | |
207 } | |
208 neighbor->seen = true; | |
209 neighbor->visitid = i; | |
210 } | |
6945 | 211 facet_list(m++) = neighbor->visitid + 1; |
6823 | 212 } |
213 } | |
214 F(r++) = facet_list; | |
215 j++; | |
216 } | |
217 | |
218 Cell C(r, 1); | |
219 for (i = 0; i < r; i++) | |
220 C.elem (i) = F(i); | |
221 | |
222 retval(0) = v; | |
223 retval(1) = C; | |
8415
fa12c67683d3
Reszie third arg of __voronoi__ before retuning it
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
224 AtInf.resize (r, 1); |
6823 | 225 retval(2) = AtInf; |
226 | |
6945 | 227 // free long memory |
228 qh_freeqhull (! qh_ALL); | |
6823 | 229 |
6945 | 230 // free short memory and memory allocator |
6823 | 231 int curlong, totlong; |
232 qh_memfreeshort (&curlong, &totlong); | |
233 | |
6945 | 234 if (curlong || totlong) |
235 warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong); | |
6823 | 236 } |
6945 | 237 else |
238 error ("__voronoi__: qhull failed"); | |
239 | |
6823 | 240 #else |
241 error ("__voronoi__: not available in this version of Octave"); | |
242 #endif | |
6945 | 243 |
6823 | 244 return retval; |
245 } |