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
convhulln.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 29. July 2000 - Kai Habel: first release
25 2002-04-22 Paul Kienzle
26 * Use warning(...) function rather than writing to cerr
27 2006-05-01 Tom Holroyd
28 * add support for consistent winding in all dimensions; output is
29 * guaranteed to be simplicial.
30 */
31 
32 #ifdef HAVE_CONFIG_H
33 #include <config.h>
34 #endif
35 
36 #include "oct-locbuf.h"
37 
38 #include "Cell.h"
39 #include "defun-dld.h"
40 #include "error.h"
41 #include "oct-obj.h"
42 #include "parse.h"
43 #include "unwind-prot.h"
44 
45 #if defined (HAVE_QHULL)
46 # include "oct-qhull.h"
47 # if defined (NEED_QHULL_VERSION)
48 char qh_version[] = "convhulln.oct 2007-07-24";
49 # endif
50 #endif
51 
52 static void
53 close_fcn (FILE *f)
54 {
55  gnulib::fclose (f);
56 }
57 
58 static bool
60 {
61  if (sizeof (octave_idx_type) > sizeof (int))
62  {
63  int maxval = std::numeric_limits<int>::max ();
64 
65  if (dim > maxval || n > maxval)
66  {
67  error ("%s: dimension too large for Qhull", who);
68  return false;
69  }
70  }
71 
72  return true;
73 }
74 
75 DEFUN_DLD (convhulln, args, nargout,
76  "-*- texinfo -*-\n\
77 @deftypefn {Loadable Function} {@var{h} =} convhulln (@var{pts})\n\
78 @deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{pts}, @var{options})\n\
79 @deftypefnx {Loadable Function} {[@var{h}, @var{v}] =} convhulln (@dots{})\n\
80 Compute the convex hull of the set of points @var{pts}.\n\
81 \n\
82 @var{pts} is a matrix of size [n, dim] containing n points in a space of\n\
83 dimension dim.\n\
84 \n\
85 The hull @var{h} is an index vector into the set of points and specifies\n\
86 which points form the enclosing hull.\n\
87 \n\
88 An optional second argument, which must be a string or cell array of strings,\n\
89 contains options passed to the underlying qhull command.\n\
90 See the documentation for the Qhull library for details\n\
91 @url{http://www.qhull.org/html/qh-quick.htm#options}.\n\
92 The default options depend on the dimension of the input:\n\
93 \n\
94 @itemize\n\
95 @item 2D, 3D, 4D: @var{options} = @code{@{\"Qt\"@}}\n\
96 \n\
97 @item 5D and higher: @var{options} = @code{@{\"Qt\", \"Qx\"@}}\n\
98 @end itemize\n\
99 \n\
100 If @var{options} is not present or @code{[]} then the default arguments are\n\
101 used. Otherwise, @var{options} replaces the default argument list.\n\
102 To append user options to the defaults it is necessary to repeat the\n\
103 default arguments in @var{options}. Use a null string to pass no arguments.\n\
104 \n\
105 If the second output @var{v} is requested the volume of the enclosing\n\
106 convex hull is calculated.\n\n\
107 @seealso{convhull, delaunayn, voronoin}\n\
108 @end deftypefn")
109 {
110  octave_value_list retval;
111 
112 #if defined (HAVE_QHULL)
113 
114  int nargin = args.length ();
115  if (nargin < 1 || nargin > 2)
116  {
117  print_usage ();
118  return retval;
119  }
120 
121  Matrix points (args(0).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, "convhulln"))
126  return retval;
127 
128  points = points.transpose ();
129 
130  std::string options;
131 
132  if (dim <= 4)
133  options = " Qt";
134  else
135  options = " Qt Qx";
136 
137  if (nargin == 2)
138  {
139  if (args(1).is_string ())
140  options = " " + args(1).string_value ();
141  else if (args(1).is_empty ())
142  ; // Use default options.
143  else if (args(1).is_cellstr ())
144  {
145  options = "";
146 
147  Array<std::string> tmp = args(1).cellstr_value ();
148 
149  for (octave_idx_type i = 0; i < tmp.numel (); i++)
150  options += " " + tmp(i);
151  }
152  else
153  {
154  error ("convhulln: OPTIONS must be a string, cell array of strings, or empty");
155  return retval;
156  }
157  }
158 
159  boolT ismalloc = false;
160 
161  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 = gnulib::fopen ("NUL", "w");
166 #else
167  FILE *outfile = gnulib::fopen ("/dev/null", "w");
168 #endif
169  FILE *errfile = stderr;
170 
171  if (outfile)
172  frame.add_fcn (close_fcn, outfile);
173  else
174  {
175  error ("convhulln: unable to create temporary file for output");
176  return retval;
177  }
178 
179  // qh_new_qhull command and points arguments are not const...
180 
181  std::string cmd = "qhull" + options;
182 
183  OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1);
184 
185  strcpy (cmd_str, cmd.c_str ());
186 
187  int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (),
188  ismalloc, cmd_str, outfile, errfile);
189  if (! exitcode)
190  {
191  bool nonsimp_seen = false;
192 
193  octave_idx_type nf = qh num_facets;
194 
195  Matrix idx (nf, dim + 1);
196 
197  facetT *facet;
198 
199  octave_idx_type i = 0;
200 
201  FORALLfacets
202  {
203  octave_idx_type j = 0;
204 
205  if (! (nonsimp_seen || facet->simplicial || qh hull_dim == 2))
206  {
207  nonsimp_seen = true;
208 
209  if (cmd.find ("QJ") != std::string::npos)
210  {
211  // Should never happen with QJ.
212  error ("convhulln: qhull failed: option 'QJ' returned non-simplicial facet");
213  return retval;
214  }
215  }
216 
217  if (dim == 3)
218  {
219  setT *vertices = qh_facet3vertex (facet);
220 
221  vertexT *vertex, **vertexp;
222 
223  FOREACHvertex_ (vertices)
224  idx(i, j++) = 1 + qh_pointid(vertex->point);
225 
226  qh_settempfree (&vertices);
227  }
228  else
229  {
230  if (facet->toporient ^ qh_ORIENTclock)
231  {
232  vertexT *vertex, **vertexp;
233 
234  FOREACHvertex_ (facet->vertices)
235  idx(i, j++) = 1 + qh_pointid(vertex->point);
236  }
237  else
238  {
239  vertexT *vertex, **vertexp;
240 
241  FOREACHvertexreverse12_ (facet->vertices)
242  idx(i, j++) = 1 + qh_pointid(vertex->point);
243  }
244  }
245  if (j < dim)
246  warning ("convhulln: facet %d only has %d vertices", i, j);
247 
248  i++;
249  }
250 
251  // Remove extra dimension if all facets were simplicial.
252 
253  if (! nonsimp_seen)
254  idx.resize (nf, dim, 0.0);
255 
256  if (nargout == 2)
257  {
258  // Calculate volume of convex hull, taken from qhull src/geom2.c.
259 
260  realT area;
261  realT dist;
262 
263  FORALLfacets
264  {
265  if (! facet->normal)
266  continue;
267 
268  if (facet->upperdelaunay && qh ATinfinity)
269  continue;
270 
271  facet->f.area = area = qh_facetarea (facet);
272  facet->isarea = True;
273 
274  if (qh DELAUNAY)
275  {
276  if (facet->upperdelaunay == qh UPPERdelaunay)
277  qh totarea += area;
278  }
279  else
280  {
281  qh totarea += area;
282  qh_distplane (qh interior_point, facet, &dist);
283  qh totvol += -dist * area/ qh hull_dim;
284  }
285  }
286 
287  retval(1) = octave_value (qh totvol);
288  }
289 
290  retval(0) = idx;
291  }
292  else
293  error ("convhulln: qhull failed");
294 
295  // Free memory from Qhull
296  qh_freeqhull (! qh_ALL);
297 
298  int curlong, totlong;
299  qh_memfreeshort (&curlong, &totlong);
300 
301  if (curlong || totlong)
302  warning ("convhulln: did not free %d bytes of long memory (%d pieces)",
303  totlong, curlong);
304 
305 #else
306  error ("convhulln: not available in this version of Octave");
307 #endif
308 
309  return retval;
310 }
311 
312 /*
313 %!testif HAVE_QHULL
314 %! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
315 %! [h, v] = convhulln (cube, "Qt");
316 %! assert (size (h), [12 3]);
317 %! h = sortrows (sort (h, 2), [1:3]);
318 %! assert (h, [1 2 4; 1 2 6; 1 4 8; 1 5 6; 1 5 8; 2 3 4; 2 3 7; 2 6 7; 3 4 7; 4 7 8; 5 6 7; 5 7 8]);
319 %! assert (v, 1, 10*eps);
320 %! [h2, v2] = convhulln (cube); % Test defaut option = "Qt"
321 %! assert (size (h2), size (h));
322 %! h2 = sortrows (sort (h2, 2), [1:3]);
323 %! assert (h2, h);
324 %! assert (v2, v, 10*eps);
325 
326 %!testif HAVE_QHULL
327 %! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
328 %! [h, v] = convhulln (cube, "QJ");
329 %! assert (size (h), [12 3]);
330 %! assert (sortrows (sort (h, 2), [1:3]), [1 2 4; 1 2 5; 1 4 5; 2 3 4; 2 3 6; 2 5 6; 3 4 8; 3 6 7; 3 7 8; 4 5 8; 5 6 8; 6 7 8]);
331 %! assert (v, 1.0, 1e6*eps);
332 
333 %!testif HAVE_QHULL
334 %! tetrahedron = [1 1 1;-1 -1 1;-1 1 -1;1 -1 -1];
335 %! [h, v] = convhulln (tetrahedron);
336 %! h = sortrows (sort (h, 2), [1 2 3]);
337 %! assert (h, [1 2 3;1 2 4; 1 3 4; 2 3 4]);
338 %! assert (v, 8/3, 10*eps);
339 
340 %!testif HAVE_QHULL
341 %! triangle=[0 0; 1 1; 1 0; 1 2];
342 %! h = convhulln (triangle);
343 %! assert (size (h), [3 2]);
344 */
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:130
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
void error(const char *fmt,...)
Definition: error.cc:476
static void close_fcn(FILE *f)
Definition: convhulln.cc:53
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const double const double * f
void add_fcn(void(*fcn)(void))
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
static bool octave_qhull_dims_ok(octave_idx_type dim, octave_idx_type n, const char *who)
Definition: convhulln.cc:59
#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
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))