Mercurial > octave-nkf

comparison libinterp/dldfcn/__voronoi__.cc @ 15195:2fc554ffbc28

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
split libinterp from src * libinterp: New directory. Move all files from src directory here except Makefile.am, main.cc, main-cli.cc, mkoctfile.in.cc, mkoctfilr.in.sh, octave-config.in.cc, octave-config.in.sh. * libinterp/Makefile.am: New file, extracted from src/Makefile.am. * src/Makefile.am: Delete everything except targets and definitions needed to build and link main and utility programs. * Makefile.am (SUBDIRS): Include libinterp in the list. * autogen.sh: Run config-module.sh in libinterp/dldfcn directory, not src/dldfcn directory. * configure.ac (AC_CONFIG_SRCDIR): Use libinterp/octave.cc, not src/octave.cc. (DL_LDFLAGS, LIBOCTINTERP): Use libinterp, not src. (AC_CONFIG_FILES): Include libinterp/Makefile in the list. * find-docstring-files.sh: Look in libinterp, not src. * gui/src/Makefile.am (liboctgui_la_CPPFLAGS): Find header files in libinterp, not src.
author John W. Eaton <jwe@octave.org>
date Sat, 18 Aug 2012 16:23:39 -0400
parents src/dldfcn/__voronoi__.cc@000587f92082
children dda043ccad7c
comparison
equal deleted inserted replaced
15194:0f0b795044c3 15195:2fc554ffbc28
1 /*
2
3 Copyright (C) 2000-2012 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 #ifdef HAVE_CONFIG_H
34 #include <config.h>
35 #endif
36
37 #include <cstdio>
38
39 #include <list>
40
41 #include "lo-ieee.h"
42
43 #include "Cell.h"
44 #include "defun-dld.h"
45 #include "error.h"
46 #include "oct-obj.h"
47 #include "unwind-prot.h"
48
49 #if defined (HAVE_QHULL)
50 # include "oct-qhull.h"
51 # if defined (NEED_QHULL_VERSION)
52 char qh_version[] = "__voronoi__.oct 2007-07-24";
53 # endif
54 #endif
55
56 static void
57 close_fcn (FILE *f)
58 {
59 gnulib::fclose (f);
60 }
61
62 DEFUN_DLD (__voronoi__, args, ,
63 "-*- texinfo -*-\n\
64 @deftypefn {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts})\n\
65 @deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts}, @var{options})\n\
66 @deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\
67 Undocumented internal function.\n\
68 @end deftypefn")
69 {
70 octave_value_list retval;
71
72 std::string caller = args(0).string_value ();
73
74 #if defined (HAVE_QHULL)
75
76 retval(0) = 0.0;
77
78 int nargin = args.length ();
79 if (nargin < 2 || nargin > 3)
80 {
81 print_usage ();
82 return retval;
83 }
84
85 Matrix points = args(1).matrix_value ();
86 const octave_idx_type dim = points.columns ();
87 const octave_idx_type num_points = points.rows ();
88
89 points = points.transpose ();
90
91 std::string options;
92
93 if (dim <= 4)
94 options = " Qbb";
95 else
96 options = " Qbb Qx";
97
98 if (nargin == 3)
99 {
100 octave_value opt_arg = args(2);
101
102 if (opt_arg.is_string ())
103 options = " " + opt_arg.string_value ();
104 else if (opt_arg.is_empty ())
105 ; // Use default options.
106 else if (opt_arg.is_cellstr ())
107 {
108 options = "";
109
110 Array<std::string> tmp = opt_arg.cellstr_value ();
111
112 for (octave_idx_type i = 0; i < tmp.numel (); i++)
113 options += " " + tmp(i);
114 }
115 else
116 {
117 error ("%s: OPTIONS must be a string, cell array of strings, or empty",
118 caller.c_str ());
119 return retval;
120 }
121 }
122
123 boolT ismalloc = false;
124
125 unwind_protect frame;
126
127 // Replace the outfile pointer with stdout for debugging information.
128 #if defined (OCTAVE_HAVE_WINDOWS_FILESYSTEM) && ! defined (OCTAVE_HAVE_POSIX_FILESYSTEM)
129 FILE *outfile = gnulib::fopen ("NUL", "w");
130 #else
131 FILE *outfile = gnulib::fopen ("/dev/null", "w");
132 #endif
133 FILE *errfile = stderr;
134
135 if (outfile)
136 frame.add_fcn (close_fcn, outfile);
137 else
138 {
139 error ("__voronoi__: unable to create temporary file for output");
140 return retval;
141 }
142
143 // qh_new_qhull command and points arguments are not const...
144
145 std::string cmd = "qhull v" + options;
146
147 OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1);
148
149 strcpy (cmd_str, cmd.c_str ());
150
151 int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (),
152 ismalloc, cmd_str, outfile, errfile);
153 if (! exitcode)
154 {
155 // Calling findgood_all provides the number of Voronoi vertices
156 // (sets qh num_good).
157
158 qh_findgood_all (qh facet_list);
159
160 octave_idx_type num_voronoi_regions
161 = qh num_vertices - qh_setsize (qh del_vertices);
162
163 octave_idx_type num_voronoi_vertices = qh num_good;
164
165 // Find the voronoi centers for all facets.
166
167 qh_setvoronoi_all ();
168
169 facetT *facet;
170 vertexT *vertex;
171 octave_idx_type k;
172
173 // Find the number of Voronoi vertices for each Voronoi cell and
174 // store them in NI so we can use them later to set the dimensions
175 // of the RowVector objects used to collect them.
176
177 FORALLfacets
178 {
179 facet->seen = false;
180 }
181
182 OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, num_voronoi_regions);
183 for (octave_idx_type i = 0; i < num_voronoi_regions; i++)
184 ni[i] = 0;
185
186 k = 0;
187
188 FORALLvertices
189 {
190 if (qh hull_dim == 3)
191 qh_order_vertexneighbors (vertex);
192
193 bool infinity_seen = false;
194
195 facetT *neighbor, **neighborp;
196
197 FOREACHneighbor_ (vertex)
198 {
199 if (neighbor->upperdelaunay)
200 {
201 if (! infinity_seen)
202 {
203 infinity_seen = true;
204 ni[k]++;
205 }
206 }
207 else
208 {
209 neighbor->seen = true;
210 ni[k]++;
211 }
212 }
213
214 k++;
215 }
216
217 // If Qhull finds fewer regions than points, we will pad the end
218 // of the at_inf and C arrays so that they always contain at least
219 // as many elements as the given points array.
220
221 // FIXME -- is it possible (or does it make sense) for
222 // num_voronoi_regions to ever be larger than num_points?
223
224 octave_idx_type nr = (num_points > num_voronoi_regions
225 ? num_points : num_voronoi_regions);
226
227 boolMatrix at_inf (nr, 1, false);
228
229 // The list of Voronoi vertices. The first element is always
230 // Inf.
231 Matrix F (num_voronoi_vertices+1, dim);
232
233 for (octave_idx_type d = 0; d < dim; d++)
234 F(0,d) = octave_Inf;
235
236 // The cell array of vectors of indices into F that represent the
237 // vertices of the Voronoi regions (cells).
238
239 Cell C (nr, 1);
240
241 // Now loop through the list of vertices again and store the
242 // coordinates of the Voronoi vertices and the lists of indices
243 // for the cells.
244
245 FORALLfacets
246 {
247 facet->seen = false;
248 }
249
250 octave_idx_type i = 0;
251 k = 0;
252
253 FORALLvertices
254 {
255 if (qh hull_dim == 3)
256 qh_order_vertexneighbors (vertex);
257
258 bool infinity_seen = false;
259
260 octave_idx_type idx = qh_pointid (vertex->point);
261
262 octave_idx_type num_vertices = ni[k++];
263
264 // Qhull seems to sometimes produces regions with a single
265 // vertex. Is that a bug? How can a region have just one
266 // vertex? Let's skip it.
267
268 if (num_vertices == 1)
269 continue;
270
271 RowVector facet_list (num_vertices);
272
273 octave_idx_type m = 0;
274
275 facetT *neighbor, **neighborp;
276
277 FOREACHneighbor_(vertex)
278 {
279 if (neighbor->upperdelaunay)
280 {
281 if (! infinity_seen)
282 {
283 infinity_seen = true;
284 facet_list(m++) = 1;
285 at_inf(idx) = true;
286 }
287 }
288 else
289 {
290 if (! neighbor->seen)
291 {
292 i++;
293 for (octave_idx_type d = 0; d < dim; d++)
294 F(i,d) = neighbor->center[d];
295
296 neighbor->seen = true;
297 neighbor->visitid = i;
298 }
299
300 facet_list(m++) = neighbor->visitid + 1;
301 }
302 }
303
304 C(idx) = facet_list;
305 }
306
307 retval(2) = at_inf;
308 retval(1) = C;
309 retval(0) = F;
310 }
311 else
312 error ("%s: qhull failed", caller.c_str ());
313
314 // Free memory from Qhull
315 qh_freeqhull (! qh_ALL);
316
317 int curlong, totlong;
318 qh_memfreeshort (&curlong, &totlong);
319
320 if (curlong || totlong)
321 warning ("%s: qhull did not free %d bytes of long memory (%d pieces)",
322 caller.c_str (), totlong, curlong);
323
324 #else
325 error ("%s: not available in this version of Octave", caller.c_str ());
326 #endif
327
328 return retval;
329 }
330
331 /*
332 ## No test needed for internal helper function.
333 %!assert (1)
334 */