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