Mercurial > octave-nkf
view src/DLD-FUNCTIONS/__voronoi__.cc @ 6880:e00a8f661f06
[project @ 2007-09-07 21:39:55 by dbateman]
author | dbateman |
---|---|
date | Fri, 07 Sep 2007 21:39:56 +0000 |
parents | 9fddcc586065 |
children | 6bbf56a9718a |
line wrap: on
line source
/* 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 #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; } const char *options; 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; OCTAVE_LOCAL_BUFFER(char, flags, 250); // hmm lot's of options for qhull here sprintf(flags,"qhull v Fv T0 %s",options); // If you want some debugging information replace the NULL // pointer with outfile. FILE *outfile = stdout; FILE *errfile = stderr; 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; }