/*
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
.
*/
/*
20. Augiust 2000 - Kai Habel: first release
*/
/*
2003-12-14 Rafael Laboissiere
Added optional second argument to pass options to the underlying
qhull command
*/
#ifdef HAVE_CONFIG_H
#include
#endif
#include
#include
#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::max ();
if (dim > maxval || n > maxval)
{
error ("%s: dimension too large for Qhull", who);
return false;
}
}
return true;
}
DEFUN_DLD (__voronoi__, args, ,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts})\n\
@deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts}, @var{options})\n\
@deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\
Undocumented internal function.\n\
@end deftypefn")
{
octave_value_list retval;
std::string caller = args(0).string_value ("__voronoi__: CALLER must be a string");
#if defined (HAVE_QHULL)
retval(0) = 0.0;
int nargin = args.length ();
if (nargin < 2 || nargin > 3)
{
print_usage ();
return 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 retval;
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 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 ());
return retval;
}
}
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)
frame.add_fcn (close_fcn, outfile);
else
{
error ("__voronoi__: unable to create temporary file for output");
return retval;
}
// 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)
{
// 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(2) = at_inf;
retval(1) = C;
retval(0) = F;
}
else
error ("%s: qhull failed", caller.c_str ());
// 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);
#else
error ("%s: not available in this version of Octave", caller.c_str ());
#endif
return retval;
}
/*
## No test needed for internal helper function.
%!assert (1)
*/