Mercurial > octave
view libinterp/dldfcn/__voronoi__.cc @ 20939:b17fda023ca6
maint: Use new C++ archetype in more files.
Place input validation first in files.
Move declaration of retval down in function to be closer to point of usage.
Eliminate else clause after if () error.
Use "return ovl()" where it makes sense.
* find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, givens.cc,
graphics.cc, help.cc, hess.cc, hex2num.cc, input.cc, kron.cc, load-path.cc,
load-save.cc, lookup.cc, mappers.cc, matrix_type.cc, mgorth.cc, nproc.cc,
ordschur.cc, pager.cc, pinv.cc, pr-output.cc, profiler.cc, psi.cc, quad.cc,
rcond.cc, regexp.cc, schur.cc, sighandlers.cc, sparse.cc, str2double.cc,
strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc, symtab.cc,
syscalls.cc, sysdep.cc, time.cc, toplev.cc, tril.cc, tsearch.cc, typecast.cc,
urlwrite.cc, utils.cc, variables.cc, __delaunayn__.cc, __eigs__.cc,
__glpk__.cc, __magick_read__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc,
audiodevinfo.cc, audioread.cc, chol.cc, colamd.cc, dmperm.cc, fftw.cc, qr.cc,
symbfact.cc, symrcm.cc, ov-bool-mat.cc, ov-cell.cc, ov-class.cc,
ov-classdef.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-flt-re-mat.cc,
ov-java.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-re-mat.cc, ov-struct.cc,
ov-typeinfo.cc, ov-usr-fcn.cc, ov.cc, octave.cc:
Use new C++ archetype in more files.
author | Rik <rik@octave.org> |
---|---|
date | Fri, 18 Dec 2015 15:37:22 -0800 |
parents | 2bd3b13e2c8e |
children | 48b2ad5ee801 |
line wrap: on
line source
/* Copyright (C) 2000-2015 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/>. */ /* 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 */ #ifdef HAVE_CONFIG_H #include <config.h> #endif #include <cstdio> #include <list> #include "oct-locbuf.h" #include "lo-ieee.h" #include "Cell.h" #include "defun-dld.h" #include "error.h" #include "oct-obj.h" #include "unwind-prot.h" #if defined (HAVE_QHULL) # include "oct-qhull.h" # if defined (NEED_QHULL_VERSION) char qh_version[] = "__voronoi__.oct 2007-07-24"; # endif #endif static void close_fcn (FILE *f) { gnulib::fclose (f); } static bool octave_qhull_dims_ok (octave_idx_type dim, octave_idx_type n, const char *who) { if (sizeof (octave_idx_type) > sizeof (int)) { int maxval = std::numeric_limits<int>::max (); if (dim > maxval || n > maxval) error ("%s: dimension too large for Qhull", who); } return true; } DEFUN_DLD (__voronoi__, args, , "-*- texinfo -*-\n\ @deftypefn {} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts})\n\ @deftypefnx {} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts}, @var{options})\n\ @deftypefnx {} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\ Undocumented internal function.\n\ @end deftypefn") { #if defined (HAVE_QHULL) int nargin = args.length (); if (nargin < 2 || nargin > 3) print_usage (); std::string caller = args(0).xstring_value ("__voronoi__: CALLER must be a string"); octave_value_list retval; Matrix points = args(1).matrix_value (); const octave_idx_type dim = points.columns (); const octave_idx_type num_points = points.rows (); if (! octave_qhull_dims_ok (dim, num_points, "__voronoi__")) return ovl (0.0); points = points.transpose (); std::string options; if (dim <= 3) options = " Qbb"; else options = " Qbb Qx"; if (nargin == 3) { octave_value opt_arg = args(2); if (opt_arg.is_string ()) options = " " + opt_arg.string_value (); else if (opt_arg.is_empty ()) ; // Use default options. else if (opt_arg.is_cellstr ()) { options = ""; Array<std::string> tmp = opt_arg.cellstr_value (); for (octave_idx_type i = 0; i < tmp.numel (); i++) options += " " + tmp(i); } else error ("%s: OPTIONS must be a string, cell array of strings, or empty", caller.c_str ()); } boolT ismalloc = false; unwind_protect frame; // Replace the outfile pointer with stdout for debugging information. #if defined (OCTAVE_HAVE_WINDOWS_FILESYSTEM) && ! defined (OCTAVE_HAVE_POSIX_FILESYSTEM) FILE *outfile = gnulib::fopen ("NUL", "w"); #else FILE *outfile = gnulib::fopen ("/dev/null", "w"); #endif FILE *errfile = stderr; if (! outfile) error ("__voronoi__: unable to create temporary file for output"); frame.add_fcn (close_fcn, outfile); // qh_new_qhull command and points arguments are not const... std::string cmd = "qhull v" + options; OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1); strcpy (cmd_str, cmd.c_str ()); int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (), ismalloc, cmd_str, outfile, errfile); if (exitcode) error ("%s: qhull failed", caller.c_str ()); // Calling findgood_all provides the number of Voronoi vertices // (sets qh num_good). qh_findgood_all (qh facet_list); octave_idx_type num_voronoi_regions = qh num_vertices - qh_setsize (qh del_vertices); octave_idx_type num_voronoi_vertices = qh num_good; // Find the voronoi centers for all facets. qh_setvoronoi_all (); facetT *facet; vertexT *vertex; octave_idx_type k; // Find the number of Voronoi vertices for each Voronoi cell and // store them in NI so we can use them later to set the dimensions // of the RowVector objects used to collect them. FORALLfacets { facet->seen = false; } OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, num_voronoi_regions); for (octave_idx_type i = 0; i < num_voronoi_regions; i++) ni[i] = 0; k = 0; FORALLvertices { if (qh hull_dim == 3) qh_order_vertexneighbors (vertex); bool infinity_seen = false; facetT *neighbor, **neighborp; FOREACHneighbor_ (vertex) { if (neighbor->upperdelaunay) { if (! infinity_seen) { infinity_seen = true; ni[k]++; } } else { neighbor->seen = true; ni[k]++; } } k++; } // If Qhull finds fewer regions than points, we will pad the end // of the at_inf and C arrays so that they always contain at least // as many elements as the given points array. // FIXME: is it possible (or does it make sense) for // num_voronoi_regions to ever be larger than num_points? octave_idx_type nr = (num_points > num_voronoi_regions ? num_points : num_voronoi_regions); boolMatrix at_inf (nr, 1, false); // The list of Voronoi vertices. The first element is always // Inf. Matrix F (num_voronoi_vertices+1, dim); for (octave_idx_type d = 0; d < dim; d++) F(0,d) = octave_Inf; // The cell array of vectors of indices into F that represent the // vertices of the Voronoi regions (cells). Cell C (nr, 1); // Now loop through the list of vertices again and store the // coordinates of the Voronoi vertices and the lists of indices // for the cells. FORALLfacets { facet->seen = false; } octave_idx_type i = 0; k = 0; FORALLvertices { if (qh hull_dim == 3) qh_order_vertexneighbors (vertex); bool infinity_seen = false; octave_idx_type idx = qh_pointid (vertex->point); octave_idx_type num_vertices = ni[k++]; // Qhull seems to sometimes produces regions with a single // vertex. Is that a bug? How can a region have just one // vertex? Let's skip it. if (num_vertices == 1) continue; RowVector facet_list (num_vertices); octave_idx_type m = 0; facetT *neighbor, **neighborp; FOREACHneighbor_(vertex) { if (neighbor->upperdelaunay) { if (! infinity_seen) { infinity_seen = true; facet_list(m++) = 1; at_inf(idx) = true; } } else { if (! neighbor->seen) { i++; for (octave_idx_type d = 0; d < dim; d++) F(i,d) = neighbor->center[d]; neighbor->seen = true; neighbor->visitid = i; } facet_list(m++) = neighbor->visitid + 1; } } C(idx) = facet_list; } retval = ovl (F, C, at_inf); // Free memory from Qhull qh_freeqhull (! qh_ALL); int curlong, totlong; qh_memfreeshort (&curlong, &totlong); if (curlong || totlong) warning ("%s: qhull did not free %d bytes of long memory (%d pieces)", caller.c_str (), totlong, curlong); return retval; #else error ("%s: not available in this version of Octave", caller.c_str ()); #endif } /* ## No test needed for internal helper function. %!assert (1) */