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