convhulln.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2000-2021 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or <https://octave.org/copyright/>.
//
// 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
// <https://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.
*/
 
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif
 
#include <limits>
#include <string>
 
#include "Array.h"
#include "dMatrix.h"
#include "oct-locbuf.h"
#include "unwind-prot.h"
 
#include "defun-dld.h"
#include "error.h"
#include "errwarn.h"
#include "ov.h"
 
#if defined (HAVE_QHULL)
 
#  include "oct-qhull.h"
 
#  if defined (NEED_QHULL_VERSION)
char qh_version[] = "convhulln.oct 2007-07-24";
#  endif
 
static void
close_fcn (FILE *f)
{
  std::fclose (f);
}
 
static void
free_qhull_memory ()
{
  qh_freeqhull (! qh_ALL);
 
  int curlong, totlong;
  qh_memfreeshort (&curlong, &totlong);
 
  if (curlong || totlong)
    warning ("convhulln: did not free %d bytes of long memory (%d pieces)",
             totlong, curlong);
}
 
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;
}
 
#endif
 
DEFUN_DLD (convhulln, args, nargout,
           doc: /* -*- texinfo -*-
@deftypefn  {} {@var{h} =} convhulln (@var{pts})
@deftypefnx {} {@var{h} =} convhulln (@var{pts}, @var{options})
@deftypefnx {} {[@var{h}, @var{v}] =} convhulln (@dots{})
Compute the convex hull of the set of points @var{pts}.
 
@var{pts} is a matrix of size [n, dim] containing n points in a space of
dimension dim.
 
The hull @var{h} is an index vector into the set of points and specifies
which points form the enclosing hull.
 
An optional second argument, which must be a string or cell array of
strings, contains options passed to the underlying qhull command.  See the
documentation for the Qhull library for details
@url{http://www.qhull.org/html/qh-quick.htm#options}.
The default options depend on the dimension of the input:
 
@itemize
@item 2D, 3D, 4D: @var{options} = @code{@{"Qt"@}}
 
@item 5D and higher: @var{options} = @code{@{"Qt", "Qx"@}}
@end itemize
 
If @var{options} is not present or @code{[]} then the default arguments are
used.  Otherwise, @var{options} replaces the default argument list.
To append user options to the defaults it is necessary to repeat the
default arguments in @var{options}.  Use a null string to pass no arguments.
 
If the second output @var{v} is requested the volume of the enclosing
convex hull is calculated.
@seealso{convhull, delaunayn, voronoin}
@end deftypefn */)
{
#if defined (HAVE_QHULL)
 
  int nargin = args.length ();
 
  if (nargin < 1 || nargin > 2)
    print_usage ();
 
  octave_value_list retval;
 
  Matrix points (args(0).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, "convhulln"))
    return retval;
 
  points = points.transpose ();
 
  std::string options;
 
  if (dim <= 4)
    options = " Qt";
  else
    options = " Qt Qx";
 
  if (nargin == 2)
    {
      if (args(1).is_string ())
        options = ' ' + args(1).string_value ();
      else if (args(1).isempty ())
        ; // Use default options.
      else if (args(1).iscellstr ())
        {
          options = "";
 
          Array<std::string> tmp = args(1).cellstr_value ();
 
          for (octave_idx_type i = 0; i < tmp.numel (); i++)
            options += ' ' + tmp(i);
        }
      else
        error ("convhulln: OPTIONS must be a string, cell array of strings, or empty");
    }
 
  boolT ismalloc = false;
 
  octave::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 = std::fopen ("NUL", "w");
#else
  FILE *outfile = std::fopen ("/dev/null", "w");
#endif
  FILE *errfile = stderr;
 
  if (! outfile)
    error ("convhulln: 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" + 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);
 
  frame.add_fcn (free_qhull_memory);
 
  if (exitcode)
    error ("convhulln: qhull failed");
 
  bool nonsimp_seen = false;
 
  octave_idx_type nf = qh num_facets;
 
  Matrix idx (nf, dim + 1);
 
  facetT *facet;
 
  octave_idx_type i = 0;
 
  FORALLfacets
    {
      octave_idx_type j = 0;
 
      if (! (nonsimp_seen || facet->simplicial || qh hull_dim == 2))
        {
          nonsimp_seen = true;
 
          if (cmd.find ("QJ") != std::string::npos)
            // Should never happen with QJ.
            error ("convhulln: qhull failed: option 'QJ' returned non-simplicial facet");
        }
 
      if (dim == 3)
        {
          setT *vertices = qh_facet3vertex (facet);
 
          vertexT *vertex, **vertexp;
 
          FOREACHvertex_ (vertices)
            idx(i, j++) = 1 + qh_pointid(vertex->point);
 
          qh_settempfree (&vertices);
        }
      else
        {
          if (facet->toporient ^ qh_ORIENTclock)
            {
              vertexT *vertex, **vertexp;
 
              FOREACHvertex_ (facet->vertices)
                idx(i, j++) = 1 + qh_pointid(vertex->point);
            }
          else
            {
              vertexT *vertex, **vertexp;
 
              FOREACHvertexreverse12_ (facet->vertices)
                idx(i, j++) = 1 + qh_pointid(vertex->point);
            }
        }
      if (j < dim)
        warning ("convhulln: facet %" OCTAVE_IDX_TYPE_FORMAT
                 " only has %" OCTAVE_IDX_TYPE_FORMAT
                 " vertices", i, j);
 
      i++;
    }
 
  // Remove extra dimension if all facets were simplicial.
 
  if (! nonsimp_seen)
    idx.resize (nf, dim, 0.0);
 
  if (nargout == 2)
    {
      // Calculate volume of convex hull, taken from qhull src/geom2.c.
 
      realT area;
      realT dist;
 
      FORALLfacets
        {
          if (! facet->normal)
            continue;
 
          if (facet->upperdelaunay && qh ATinfinity)
            continue;
 
          facet->f.area = area = qh_facetarea (facet);
          facet->isarea = True;
 
          if (qh DELAUNAY)
            {
              if (facet->upperdelaunay == qh UPPERdelaunay)
                qh totarea += area;
            }
          else
            {
              qh totarea += area;
              qh_distplane (qh interior_point, facet, &dist);
              qh totvol += -dist * area/ qh hull_dim;
            }
        }
 
      retval(1) = octave_value (qh totvol);
    }
 
  retval(0) = idx;
 
  return retval;
 
#else
 
  octave_unused_parameter (args);
  octave_unused_parameter (nargout);
 
  err_disabled_feature ("convhulln", "Qhull");
 
#endif
}
 
/*
%!testif HAVE_QHULL
%! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
%! [h, v] = convhulln (cube, "Qt");
%! assert (size (h), [12 3]);
%! h = sortrows (sort (h, 2), [1:3]);
%! assert (h, [1 2 4; 1 2 6; 1 4 8; 1 5 6; 1 5 8; 2 3 4; 2 3 7; 2 6 7; 3 4 7; 4 7 8; 5 6 7; 5 7 8]);
%! assert (v, 1, 10*eps);
%! [h2, v2] = convhulln (cube);  # Test default option = "Qt"
%! assert (size (h2), size (h));
%! h2 = sortrows (sort (h2, 2), [1:3]);
%! assert (h2, h);
%! assert (v2, v, 10*eps);
 
%!testif HAVE_QHULL
%! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
%! [h, v] = convhulln (cube, "QJ");
%! assert (size (h), [12 3]);
%! assert (sortrows (sort (h, 2), [1:3]), [1 2 4; 1 2 5; 1 4 5; 2 3 4; 2 3 6; 2 5 6; 3 4 8; 3 6 7; 3 7 8; 4 5 8; 5 6 8; 6 7 8]);
%! assert (v, 1.0, 1e6*eps);
 
%!testif HAVE_QHULL
%! tetrahedron = [1 1 1;-1 -1 1;-1 1 -1;1 -1 -1];
%! [h, v] = convhulln (tetrahedron);
%! h = sortrows (sort (h, 2), [1 2 3]);
%! assert (h, [1 2 3;1 2 4; 1 3 4; 2 3 4]);
%! assert (v, 8/3, 10*eps);
 
%!testif HAVE_QHULL
%! triangle = [0 0; 1 1; 1 0; 1 2];
%! h = convhulln (triangle);
%! assert (size (h), [3 2]);
*/