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