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