convhulln.cc

Go to the documentation of this file.

/*

Copyright (C) 2000-2018 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
<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 <iostream>
#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 %d only has %d 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]);
*/