Mercurial > octave-nkf
diff src/DLD-FUNCTIONS/__voronoi__.cc @ 6823:9fddcc586065
[project @ 2007-08-24 08:27:27 by dbateman]
author | dbateman |
---|---|
date | Fri, 24 Aug 2007 08:27:29 +0000 |
parents | |
children | e00a8f661f06 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/__voronoi__.cc Fri Aug 24 08:27:29 2007 +0000 @@ -0,0 +1,228 @@ +/* + +Copyright (C) 2000 Kai Habel + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. + +*/ + +/* +20. Augiust 2000 - Kai Habel: first release +*/ + +/* +2003-12-14 Rafael Laboissiere <rafael@laboissiere.net> +Added optional second argument to pass options to the underlying +qhull command +*/ + +#include <iostream> +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif +#include "oct.h" +#include "Cell.h" + +#ifdef HAVE_QHULL +extern "C" { +#include <qhull/qhull_a.h> +} + +#ifdef NEED_QHULL_VERSION +char qh_version[] = "__voronoi__.oct 2007-07-24"; +#endif + +FILE *outfile = stdout; +FILE *errfile = stderr; +char flags[250]; +const char *options; +#endif + +DEFUN_DLD (__voronoi__, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\ +@deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\ +Internal function for voronoi.\n\ +@end deftypefn") +{ + octave_value_list retval; +#ifdef HAVE_QHULL + retval(0) = 0.0; + + int nargin = args.length(); + if (nargin < 1 || nargin > 2) + { + print_usage (); + return retval; + } + + if (nargin == 2) + { + if (!args (1).is_string ()) + { + error ("__voronoi__: second argument must be a string"); + return retval; + } + options = args (1).string_value().c_str(); + } + else + options = ""; + + Matrix p(args(0).matrix_value()); + + const octave_idx_type dim = p.columns(); + const octave_idx_type np = p.rows(); + p = p.transpose(); + + double *pt_array = p.fortran_vec(); + //double pt_array[dim * np]; + //for (int i = 0; i < np; i++) + // { + // for (int j = 0; j < dim; j++) + // { + // pt_array[j+i*dim] = p(i,j); + // } + // } + + boolT ismalloc = false; + + // hmm lot's of options for qhull here + sprintf(flags,"qhull v Fv T0 %s",options); + if (!qh_new_qhull (dim, np, pt_array, ismalloc, flags, NULL, errfile)) + { + + // If you want some debugging information replace the NULL + // pointer with outfile. + facetT *facet; + vertexT *vertex; + + octave_idx_type i = 0, n = 1, k = 0, m = 0, fidx = 0, j = 0, r = 0; + OCTAVE_LOCAL_BUFFER(octave_idx_type, ni, np); + + for (i = 0; i < np; i++) + ni[i] = 0; + qh_setvoronoi_all(); + bool infinity_seen = false; + facetT *neighbor,**neighborp; + coordT *voronoi_vertex; + FORALLfacets + { + facet->seen = false; + } + FORALLvertices + { + if (qh hull_dim == 3) + qh_order_vertexneighbors(vertex); + infinity_seen = false; + FOREACHneighbor_(vertex) + { + if (!neighbor->upperdelaunay) + { + if (!neighbor->seen) + { + n++; + neighbor->seen=True; + } + ni[k]++; + } + else if (!infinity_seen) + { + infinity_seen = true; + ni[k]++; + } + } + k++; + } + + Matrix v(n, dim); + for (octave_idx_type d = 0; d < dim; d++) + v(0,d) = octave_Inf; + + boolMatrix AtInf(np, 1); + for (i = 0; i < np; i++) + AtInf(i) = false; + octave_value_list F; + k = 0; + i = 0; + FORALLfacets + { + facet->seen = false; + } + FORALLvertices + { + if (qh hull_dim == 3) + qh_order_vertexneighbors(vertex); + infinity_seen = false; + RowVector facet_list(ni[k++]); + m = 0; + FOREACHneighbor_(vertex) + { + if (neighbor->upperdelaunay) + { + if (!infinity_seen) + { + infinity_seen = true; + facet_list(m++) = 1; + AtInf(j) = true; + } + } + else + { + if (!neighbor->seen) + { + voronoi_vertex = neighbor->center; + fidx = neighbor->id; + i++; + for (octave_idx_type d = 0; d < dim; d++) + { + v(i,d) = *(voronoi_vertex++); + } + neighbor->seen = true; + neighbor->visitid = i; + } + facet_list(m++)=neighbor->visitid + 1; + } + } + F(r++) = facet_list; + j++; + } + + Cell C(r, 1); + for (i = 0; i < r; i++) + C.elem (i) = F(i); + + retval(0) = v; + retval(1) = C; + retval(2) = AtInf; + + //free long memory + qh_freeqhull(!qh_ALL); + + //free short memory and memory allocator + int curlong, totlong; + qh_memfreeshort (&curlong, &totlong); + + if (curlong || totlong) { + warning("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong); + } + } +#else + error ("__voronoi__: not available in this version of Octave"); +#endif + return retval; +}