1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2000-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 /*
27 20. Augiust 2000 - Kai Habel: first release
28 */
29
30 /*
31 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
32 Added optional second argument to pass options to the underlying
33 qhull command
34 */
35
36 #if defined (HAVE_CONFIG_H)
37 # include "config.h"
38 #endif
39
40 #include <cstdio>
41
42 #include <limits>
43 #include <string>
44
45 #include "Array.h"
46 #include "boolMatrix.h"
47 #include "dMatrix.h"
48 #include "dRowVector.h"
49 #include "oct-locbuf.h"
50 #include "unwind-prot.h"
51
52 #include "Cell.h"
53 #include "defun-dld.h"
54 #include "error.h"
55 #include "errwarn.h"
56 #include "ov.h"
57 #include "ovl.h"
58
59 #if defined (HAVE_QHULL)
60
61 # include "oct-qhull.h"
62
63 # if defined (NEED_QHULL_VERSION)
64 char qh_version[] = "__voronoi__.oct 2007-07-24";
65 # endif
66
67 static void
close_fcn(FILE * f)68 close_fcn (FILE *f)
69 {
70 std::fclose (f);
71 }
72
73 static void
free_qhull_memory()74 free_qhull_memory ()
75 {
76 qh_freeqhull (! qh_ALL);
77
78 int curlong, totlong;
79 qh_memfreeshort (&curlong, &totlong);
80
81 if (curlong || totlong)
82 warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)",
83 totlong, curlong);
84 }
85
86 static bool
octave_qhull_dims_ok(octave_idx_type dim,octave_idx_type n,const char * who)87 octave_qhull_dims_ok (octave_idx_type dim, octave_idx_type n, const char *who)
88 {
89 if (sizeof (octave_idx_type) > sizeof (int))
90 {
91 int maxval = std::numeric_limits<int>::max ();
92
93 if (dim > maxval || n > maxval)
94 error ("%s: dimension too large for Qhull", who);
95 }
96
97 return true;
98 }
99
100 #endif
101
102 DEFUN_DLD (__voronoi__, args, ,
103 doc: /* -*- texinfo -*-
104 @deftypefn {} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts})
105 @deftypefnx {} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts}, @var{options})
106 @deftypefnx {} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})
107 Undocumented internal function.
108 @end deftypefn */)
109 {
110 #if defined (HAVE_QHULL)
111
112 int nargin = args.length ();
113
114 if (nargin < 2 || nargin > 3)
115 print_usage ();
116
117 std::string caller = args(0).xstring_value ("__voronoi__: CALLER must be a string");
118
119 octave_value_list retval;
120
121 Matrix points = args(1).matrix_value ();
122 const octave_idx_type dim = points.columns ();
123 const octave_idx_type num_points = points.rows ();
124
125 if (! octave_qhull_dims_ok (dim, num_points, "__voronoi__"))
126 return ovl (0.0);
127
128 points = points.transpose ();
129
130 std::string options;
131
132 if (dim <= 3)
133 options = " Qbb";
134 else
135 options = " Qbb Qx";
136
137 if (nargin == 3)
138 {
139 octave_value opt_arg = args(2);
140
141 if (opt_arg.is_string ())
142 options = ' ' + opt_arg.string_value ();
143 else if (opt_arg.isempty ())
144 ; // Use default options.
145 else if (opt_arg.iscellstr ())
146 {
147 options = "";
148
149 Array<std::string> tmp = opt_arg.cellstr_value ();
150
151 for (octave_idx_type i = 0; i < tmp.numel (); i++)
152 options += ' ' + tmp(i);
153 }
154 else
155 error ("%s: OPTIONS must be a string, cell array of strings, or empty",
156 caller.c_str ());
157 }
158
159 boolT ismalloc = false;
160
161 octave::unwind_protect frame;
162
163 // Replace the outfile pointer with stdout for debugging information.
164 #if defined (OCTAVE_HAVE_WINDOWS_FILESYSTEM) && ! defined (OCTAVE_HAVE_POSIX_FILESYSTEM)
165 FILE *outfile = std::fopen ("NUL", "w");
166 #else
167 FILE *outfile = std::fopen ("/dev/null", "w");
168 #endif
169 FILE *errfile = stderr;
170
171 if (! outfile)
172 error ("__voronoi__: unable to create temporary file for output");
173
174 frame.add_fcn (close_fcn, outfile);
175
176 // qh_new_qhull command and points arguments are not const...
177
178 std::string cmd = "qhull v" + options;
179
180 OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1);
181
182 strcpy (cmd_str, cmd.c_str ());
183
184 int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (),
185 ismalloc, cmd_str, outfile, errfile);
186
187 frame.add_fcn (free_qhull_memory);
188
189 if (exitcode)
190 error ("%s: qhull failed", caller.c_str ());
191
192 // Calling findgood_all provides the number of Voronoi vertices
193 // (sets qh num_good).
194
195 qh_findgood_all (qh facet_list);
196
197 octave_idx_type num_voronoi_regions
198 = qh num_vertices - qh_setsize (qh del_vertices);
199
200 octave_idx_type num_voronoi_vertices = qh num_good;
201
202 // Find the voronoi centers for all facets.
203
204 qh_setvoronoi_all ();
205
206 facetT *facet;
207 vertexT *vertex;
208 octave_idx_type k;
209
210 // Find the number of Voronoi vertices for each Voronoi cell and
211 // store them in NI so we can use them later to set the dimensions
212 // of the RowVector objects used to collect them.
213
214 FORALLfacets
215 {
216 facet->seen = false;
217 }
218
219 OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, num_voronoi_regions);
220 for (octave_idx_type i = 0; i < num_voronoi_regions; i++)
221 ni[i] = 0;
222
223 k = 0;
224
225 FORALLvertices
226 {
227 if (qh hull_dim == 3)
228 qh_order_vertexneighbors (vertex);
229
230 bool infinity_seen = false;
231
232 facetT *neighbor, **neighborp;
233
234 FOREACHneighbor_ (vertex)
235 {
236 if (neighbor->upperdelaunay)
237 {
238 if (! infinity_seen)
239 {
240 infinity_seen = true;
241 ni[k]++;
242 }
243 }
244 else
245 {
246 neighbor->seen = true;
247 ni[k]++;
248 }
249 }
250
251 k++;
252 }
253
254 // If Qhull finds fewer regions than points, we will pad the end
255 // of the at_inf and C arrays so that they always contain at least
256 // as many elements as the given points array.
257
258 // FIXME: is it possible (or does it make sense) for
259 // num_voronoi_regions to ever be larger than num_points?
260
261 octave_idx_type nr = (num_points > num_voronoi_regions
262 ? num_points : num_voronoi_regions);
263
264 boolMatrix at_inf (nr, 1, false);
265
266 // The list of Voronoi vertices. The first element is always
267 // Inf.
268 Matrix F (num_voronoi_vertices+1, dim);
269
270 for (octave_idx_type d = 0; d < dim; d++)
271 F(0,d) = octave::numeric_limits<double>::Inf ();
272
273 // The cell array of vectors of indices into F that represent the
274 // vertices of the Voronoi regions (cells).
275
276 Cell C (nr, 1);
277
278 // Now loop through the list of vertices again and store the
279 // coordinates of the Voronoi vertices and the lists of indices
280 // for the cells.
281
282 FORALLfacets
283 {
284 facet->seen = false;
285 }
286
287 octave_idx_type i = 0;
288 k = 0;
289
290 FORALLvertices
291 {
292 if (qh hull_dim == 3)
293 qh_order_vertexneighbors (vertex);
294
295 bool infinity_seen = false;
296
297 octave_idx_type idx = qh_pointid (vertex->point);
298
299 octave_idx_type num_vertices = ni[k++];
300
301 // Qhull seems to sometimes produces regions with a single
302 // vertex. Is that a bug? How can a region have just one
303 // vertex? Let's skip it.
304
305 if (num_vertices == 1)
306 continue;
307
308 RowVector facet_list (num_vertices);
309
310 octave_idx_type m = 0;
311
312 facetT *neighbor, **neighborp;
313
314 FOREACHneighbor_(vertex)
315 {
316 if (neighbor->upperdelaunay)
317 {
318 if (! infinity_seen)
319 {
320 infinity_seen = true;
321 facet_list(m++) = 1;
322 at_inf(idx) = true;
323 }
324 }
325 else
326 {
327 if (! neighbor->seen)
328 {
329 i++;
330 for (octave_idx_type d = 0; d < dim; d++)
331 F(i,d) = neighbor->center[d];
332
333 neighbor->seen = true;
334 neighbor->visitid = i;
335 }
336
337 facet_list(m++) = neighbor->visitid + 1;
338 }
339 }
340
341 C(idx) = facet_list;
342 }
343
344 retval = ovl (F, C, at_inf);
345
346 return retval;
347
348 #else
349
350 octave_unused_parameter (args);
351
352 std::string caller
353 = (args.length () > 0
354 ? args(0).xstring_value ("__voronoi__: CALLER must be a string")
355 : "__voronoi__");
356
357 err_disabled_feature (caller, "Qhull");
358
359 #endif
360 }
361
362 /*
363 ## No test needed for internal helper function.
364 %!assert (1)
365 */
366