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
__delaunayn__.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  16. July 2000 - Kai Habel: first release
25 
26  25. September 2002 - Changes by Rafael Laboissiere <[email protected]>
27 
28  * Added Qbb option to normalize the input and avoid crashes in Octave.
29  * delaunayn accepts now a second (optional) argument that must be a string
30  containing extra options to the qhull command.
31  * Fixed doc string. The dimension of the result matrix is [m, dim+1], and
32  not [n, dim-1].
33 
34  6. June 2006: Changes by Alexander Barth <[email protected]>
35 
36  * triangulate non-simplicial facets
37  * allow options to be specified as cell array of strings
38  * change the default options (for compatibility with matlab)
39 */
40 
41 #ifdef HAVE_CONFIG_H
42 #include <config.h>
43 #endif
44 
45 #include <iostream>
46 #include <string>
47 
48 #include "oct-locbuf.h"
49 
50 #include "Cell.h"
51 #include "defun-dld.h"
52 #include "error.h"
53 #include "oct-obj.h"
54 #include "unwind-prot.h"
55 
56 #if defined (HAVE_QHULL)
57 # include "oct-qhull.h"
58 # if defined (NEED_QHULL_VERSION)
59 char qh_version[] = "__delaunayn__.oct 2007-08-21";
60 # endif
61 #endif
62 
63 static void
64 close_fcn (FILE *f)
65 {
66  gnulib::fclose (f);
67 }
68 
69 static bool
71 {
72  if (sizeof (octave_idx_type) > sizeof (int))
73  {
74  int maxval = std::numeric_limits<int>::max ();
75 
76  if (dim > maxval || n > maxval)
77  {
78  error ("%s: dimension too large for Qhull", who);
79  return false;
80  }
81  }
82 
83  return true;
84 }
85 
86 DEFUN_DLD (__delaunayn__, args, ,
87  "-*- texinfo -*-\n\
88 @deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts})\n\
89 @deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts}, @var{options})\n\
90 Undocumented internal function.\n\
91 @end deftypefn")
92 
93 {
94  octave_value_list retval;
95 
96 #if defined (HAVE_QHULL)
97 
98  retval(0) = 0.0;
99 
100  int nargin = args.length ();
101  if (nargin < 1 || nargin > 2)
102  {
103  print_usage ();
104  return retval;
105  }
106 
107  Matrix p (args(0).matrix_value ());
108  const octave_idx_type dim = p.columns ();
109  const octave_idx_type n = p.rows ();
110 
111  if (! octave_qhull_dims_ok (dim, n, "__delaynayn__"))
112  return retval;
113 
114  // Default options
115  std::string options;
116  if (dim <= 3)
117  options = "Qt Qbb Qc Qz";
118  else
119  options = "Qt Qbb Qc Qx";
120 
121  if (nargin == 2)
122  {
123  if (args(1).is_string ())
124  options = args(1).string_value ();
125  else if (args(1).is_empty ())
126  ; // Use default options
127  else if (args(1).is_cellstr ())
128  {
129  options = "";
130  Array<std::string> tmp = args(1).cellstr_value ();
131 
132  for (octave_idx_type i = 0; i < tmp.numel (); i++)
133  options += tmp(i) + " ";
134  }
135  else
136  {
137  error ("__delaunayn__: OPTIONS argument must be a string, cell array of strings, or empty");
138  return retval;
139  }
140  }
141 
142  if (n > dim + 1)
143  {
144  p = p.transpose ();
145  double *pt_array = p.fortran_vec ();
146  boolT ismalloc = false;
147 
148  // Qhull flags argument is not const char*
149  OCTAVE_LOCAL_BUFFER (char, flags, 9 + options.length ());
150 
151  sprintf (flags, "qhull d %s", options.c_str ());
152 
153  unwind_protect frame;
154 
155  // Replace the outfile pointer with stdout for debugging information.
156 #if defined (OCTAVE_HAVE_WINDOWS_FILESYSTEM) && ! defined (OCTAVE_HAVE_POSIX_FILESYSTEM)
157  FILE *outfile = gnulib::fopen ("NUL", "w");
158 #else
159  FILE *outfile = gnulib::fopen ("/dev/null", "w");
160 #endif
161  FILE *errfile = stderr;
162 
163  if (outfile)
164  frame.add_fcn (close_fcn, outfile);
165  else
166  {
167  error ("__delaunayn__: unable to create temporary file for output");
168  return retval;
169  }
170 
171  int exitcode = qh_new_qhull (dim, n, pt_array,
172  ismalloc, flags, outfile, errfile);
173  if (! exitcode)
174  {
175  // triangulate non-simplicial facets
176  qh_triangulate ();
177 
178  facetT *facet;
179  vertexT *vertex, **vertexp;
180  octave_idx_type nf = 0;
181  octave_idx_type i = 0;
182 
183  FORALLfacets
184  {
185  if (! facet->upperdelaunay)
186  nf++;
187 
188  // Double check. Non-simplicial facets will cause segfault below
189  if (! facet->simplicial)
190  {
191  error ("__delaunayn__: Qhull returned non-simplicial facets -- try delaunayn with different options");
192  exitcode = 1;
193  break;
194  }
195  }
196 
197  if (! exitcode)
198  {
199  Matrix simpl (nf, dim+1);
200 
201  FORALLfacets
202  {
203  if (! facet->upperdelaunay)
204  {
205  octave_idx_type j = 0;
206 
207  FOREACHvertex_ (facet->vertices)
208  {
209  simpl(i, j++) = 1 + qh_pointid(vertex->point);
210  }
211  i++;
212  }
213  }
214 
215  retval(0) = simpl;
216  }
217  }
218  else
219  error ("__delaunayn__: qhull failed");
220 
221  // Free memory from Qhull
222  qh_freeqhull (! qh_ALL);
223 
224  int curlong, totlong;
225  qh_memfreeshort (&curlong, &totlong);
226 
227  if (curlong || totlong)
228  warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)",
229  totlong, curlong);
230  }
231  else if (n == dim + 1)
232  {
233  // one should check if nx points span a simplex
234  // I will look at this later.
235  RowVector vec (n);
236  for (octave_idx_type i = 0; i < n; i++)
237  vec(i) = i + 1.0;
238 
239  retval(0) = vec;
240  }
241 
242 #else
243  error ("__delaunayn__: not available in this version of Octave");
244 #endif
245 
246  return retval;
247 }
248 
249 /*
250 ## No test needed for internal helper function.
251 %!assert (1)
252 */
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
octave_idx_type rows(void) const
Definition: Array.h:313
static void close_fcn(FILE *f)
F77_RET_T const double const double * f
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
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
static bool octave_qhull_dims_ok(octave_idx_type dim, octave_idx_type n, const char *who)
#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