|
6823
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2000 Kai Habel |
|
|
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 |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
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 |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 /* |
|
|
25 20. Augiust 2000 - Kai Habel: first release |
|
|
26 */ |
|
|
27 |
|
|
28 /* |
|
|
29 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net> |
|
|
30 Added optional second argument to pass options to the underlying |
|
|
31 qhull command |
|
|
32 */ |
|
|
33 |
|
|
34 #include <iostream> |
|
|
35 #ifdef HAVE_CONFIG_H |
|
|
36 #include <config.h> |
|
|
37 #endif |
|
|
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; |
|
|
169 octave_value_list F; |
|
|
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; |
|
|
224 retval(2) = AtInf; |
|
|
225 |
|
6945
|
226 // free long memory |
|
|
227 qh_freeqhull (! qh_ALL); |
|
6823
|
228 |
|
6945
|
229 // free short memory and memory allocator |
|
6823
|
230 int curlong, totlong; |
|
|
231 qh_memfreeshort (&curlong, &totlong); |
|
|
232 |
|
6945
|
233 if (curlong || totlong) |
|
|
234 warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong); |
|
6823
|
235 } |
|
6945
|
236 else |
|
|
237 error ("__voronoi__: qhull failed"); |
|
|
238 |
|
6823
|
239 #else |
|
|
240 error ("__voronoi__: not available in this version of Octave"); |
|
|
241 #endif |
|
6945
|
242 |
|
6823
|
243 return retval; |
|
|
244 } |
