Mercurial > octave-nkf
view src/DLD-FUNCTIONS/convhulln.cc @ 7301:89d546610556 ss-2-9-19
[project @ 2007-12-12 03:56:59 by jwe]
author | jwe |
---|---|
date | Wed, 12 Dec 2007 03:57:00 +0000 |
parents | a1dbe9d80eee |
children | b166043585a8 |
line wrap: on
line source
/* Copyright (C) 2000, 2007 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 3 of the License, 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, see <http://www.gnu.org/licenses/>. */ /* 29. July 2000 - Kai Habel: first release 2002-04-22 Paul Kienzle * Use warning(...) function rather than writing to cerr 2006-05-01 Tom Holroyd * add support for consistent winding in all dimensions; output is * guaranteed to be simplicial. */ #include <sstream> #ifdef HAVE_CONFIG_H #include <config.h> #endif #include "oct.h" #include "Cell.h" #ifdef HAVE_QHULL #if defined(HAVE__SNPRINTF) && !defined(HAVE_SNPRINTF) #define snprintf _snprintf #endif extern "C" { #include <qhull/qhull_a.h> } #ifdef NEED_QHULL_VERSION char qh_version[] = "convhulln.oct 2007-07-24"; #endif #endif DEFUN_DLD (convhulln, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{H} =} convhulln (@var{p})\n\ @deftypefnx {Loadable Function} {@var{H} =} convhulln (@var{p}, @var{opt})\n\ Returns an index vector to the points of the enclosing convex hull.\n\ The input matrix of size [n, dim] contains n points of dimension dim.\n\n\ If a second optional argument is given, it must be a string or cell array\n\ of strings containing options for the underlying qhull command. (See\n\ the Qhull documentation for the available options.) The default options\n\ are \"s Qci Tcv\".\n\n\ @seealso{convhull, delaunayn}\n\ @end deftypefn") { octave_value_list retval; #ifdef HAVE_QHULL std::string options; int nargin = args.length (); if (nargin < 1 || nargin > 2) { print_usage (); return retval; } if (nargin == 2) { if (args (1).is_string ()) options = args(1).string_value (); else if (args(1).is_cell ()) { Cell c = args(1).cell_value (); options = ""; for (octave_idx_type i = 0; i < c.numel (); i++) { if (! c.elem(i).is_string ()) { error ("convhulln: second argument must be a string or cell array of strings"); return retval; } options = options + c.elem(i).string_value() + " "; } } else { error ("convhulln: second argument must be a string or cell array of strings"); return retval; } } else // turn on some consistency checks options = "s Qci Tcv"; Matrix p (args(0).matrix_value ()); const octave_idx_type dim = p.columns (); const octave_idx_type n = p.rows (); p = p.transpose (); double *pt_array = p.fortran_vec (); boolT ismalloc = False; OCTAVE_LOCAL_BUFFER (char, flags, 250); // hmm, lots of options for qhull here // QJ guarantees that the output will be triangles snprintf (flags, 250, "qhull QJ %s", options.c_str ()); if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, NULL, stderr)) { // If you want some debugging information replace the NULL // pointer with stdout vertexT *vertex, **vertexp; facetT *facet; setT *vertices; unsigned int nf = qh num_facets; Matrix idx (nf, dim); octave_idx_type j, i = 0; FORALLfacets { j = 0; if (! facet->simplicial) // should never happen with QJ error ("convhulln: non-simplicial facet"); if (dim == 3) { vertices = qh_facet3vertex (facet); FOREACHvertex_ (vertices) idx(i, j++) = 1 + qh_pointid(vertex->point); qh_settempfree (&vertices); } else { if (facet->toporient ^ qh_ORIENTclock) { FOREACHvertex_ (facet->vertices) idx(i, j++) = 1 + qh_pointid(vertex->point); } else { FOREACHvertexreverse12_ (facet->vertices) idx(i, j++) = 1 + qh_pointid(vertex->point); } } if (j < dim) // likewise but less fatal warning ("facet %d only has %d vertices", i, j); i++; } retval(0) = octave_value (idx); } // free long memory qh_freeqhull (! qh_ALL); // free short memory and memory allocator int curlong, totlong; qh_memfreeshort (&curlong, &totlong); if (curlong || totlong) warning ("convhulln: did not free %d bytes of long memory (%d pieces)", totlong, curlong); #else error ("convhulln: not available in this version of Octave"); #endif return retval; }