Mercurial > hg > octave-avbm
changeset 13862:6d7e133a4bed
compatibility fixes for __voronoi__
* __voronoi__.cc (F__voronoi__): Use Fv option for Qhull, not FV.
Delete unused variable fidx. Count vertices to get size of NI array.
Skip facets that contain only one point. Always return AtInf. Use a
list of accumulate vertex lists. Pad the cell array of facet vertex
lists with empty matrices if there are fewer facets than points.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 14 Nov 2011 04:14:04 -0500 |
parents | 4d8317915db9 |
children | 7908b27de857 |
files | src/DLD-FUNCTIONS/__voronoi__.cc |
diffstat | 1 files changed, 47 insertions(+), 19 deletions(-) [+] |
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/__voronoi__.cc +++ b/src/DLD-FUNCTIONS/__voronoi__.cc @@ -36,6 +36,8 @@ #include <cstdio> +#include <list> + #include "lo-ieee.h" #include "Cell.h" @@ -53,7 +55,7 @@ #endif #endif -DEFUN_DLD (__voronoi__, args, nargout, +DEFUN_DLD (__voronoi__, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts})\n\ @deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts}, @var{options})\n\ @@ -113,7 +115,7 @@ // Qhull flags argument is not const char* OCTAVE_LOCAL_BUFFER (char, flags, 12 + options.length ()); - sprintf (flags, "qhull v FV %s", options.c_str ()); + sprintf (flags, "qhull v Fv %s", options.c_str ()); int exitcode = qh_new_qhull (dim, np, pt_array, ismalloc, flags, outfile, errfile); @@ -122,11 +124,21 @@ 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); + octave_idx_type i = 0, n = 1, k = 0, m = 0, j = 0, r = 0; - for (i = 0; i < np; i++) + // Count number of vertices for size of NI. FIXME -- does Qhull + // have a way to query this value directly? + octave_idx_type nv = 0; + FORALLvertices + { + nv++; + } + + OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, nv); + + for (i = 0; i < nv; i++) ni[i] = 0; + qh_setvoronoi_all (); bool infinity_seen = false; @@ -162,6 +174,7 @@ ni[k]++; } } + k++; } @@ -169,8 +182,8 @@ for (octave_idx_type d = 0; d < dim; d++) v(0,d) = octave_Inf; - boolMatrix AtInf (np, 1, false); - octave_value_list F (np, octave_value ()); + boolMatrix AtInf (nv, 1, false); + std::list<RowVector> F; k = 0; i = 0; @@ -182,10 +195,20 @@ FORALLvertices { if (qh hull_dim == 3) - qh_order_vertexneighbors(vertex); + qh_order_vertexneighbors (vertex); infinity_seen = false; - RowVector facet_list (ni[k++]); + + octave_idx_type n_vertices = ni[k++]; + + // Qhull seems to sometimes produce "facets" with a single + // vertex. Is that a bug? How can a facet have just one + // vertex? Let's skip it. + + if (n_vertices == 1) + continue; + + RowVector facet_list (n_vertices); m = 0; FOREACHneighbor_(vertex) @@ -204,7 +227,6 @@ if (! neighbor->seen) { voronoi_vertex = neighbor->center; - fidx = neighbor->id; i++; for (octave_idx_type d = 0; d < dim; d++) { @@ -216,19 +238,25 @@ facet_list(m++) = neighbor->visitid + 1; } } - F(r++) = facet_list; + F.push_back (facet_list); j++; } - Cell C(r, 1); - for (i = 0; i < r; i++) - C.elem (i) = F(i); + // For compatibility with Matlab, pad the cell array of vertex + // lists with empty matrices if there are fewer facets than + // points. + octave_idx_type f_len = F.size (); + Cell C(np > f_len ? np : f_len, 1); - if (nargout == 3) - { - AtInf.resize (r, 1); - retval(2) = AtInf; - } + i = 0; + for (std::list<RowVector>::const_iterator it = F.begin (); + it != F.end (); it++) + C.elem (i++) = *it; + + F.clear (); + + AtInf.resize (f_len, 1); + retval(2) = AtInf; retval(1) = C; retval(0) = v; }