DGtal::QuickHull< TKernel > Struct Template Reference

Aim: Implements the quickhull algorithm by Barber et al. [9], a famous arbitrary dimensional convex hull computation algorithm. It relies on dedicated geometric kernels for computing and comparing facet geometries. More...

#include <DGtal/geometry/tools/QuickHull.h>

Inheritance diagram for DGtal::QuickHull< TKernel >:

Data Structures
struct	Facet

Public Types
enum	{ UNASSIGNED = (Index) -1 }
	Label for points that are not assigned to any facet. More...
enum class	Status {   Uninitialized = 0 , InputInitialized = 1 , SimplexCompleted = 2 , FacetsCompleted = 3 ,   VerticesCompleted = 4 , AllCompleted = 5 , NotFullDimensional = 10 , InvalidRidge = 11 ,   InvalidConvexHull = 12 }
	Represents the status of a QuickHull object. More...
typedef TKernel	Kernel
typedef Kernel::CoordinatePoint	Point
typedef Kernel::CoordinateVector	Vector
typedef Kernel::CoordinateScalar	Scalar
typedef Kernel::InternalScalar	InternalScalar
typedef std::size_t	Index
typedef std::size_t	Size
typedef std::vector< Index >	IndexRange
typedef Kernel::HalfSpace	HalfSpace
typedef Kernel::CombinatorialPlaneSimplex	CombinatorialPlaneSimplex
typedef std::pair< Index, Index >	Ridge

Public Member Functions
	BOOST_STATIC_ASSERT ((Point::dimension==Vector::dimension))
Standard services (construction, initialization, accessors)
	QuickHull (const Kernel &K=Kernel(), int dbg=0)
Status	status () const
void	clear ()
	Clears the object.
Size	memory () const
Size	nbPoints () const
Size	nbFacets () const
Size	nbVertices () const
Size	nbFiniteFacets () const
Size	nbInfiniteFacets () const
Initialization services
template<typename InputPoint>
bool	setInput (const std::vector< InputPoint > &input_points, bool remove_duplicates=true)
bool	setInitialSimplex (const IndexRange &full_splx)
Convex hull services
bool	computeConvexHull (Status target=Status::VerticesCompleted)
bool	computeInitialSimplex ()
bool	computeFacets ()
bool	computeVertices ()
Output services
template<typename OutputPoint>
bool	getVertexPositions (std::vector< OutputPoint > &vertex_positions)
bool	getVertex2Point (IndexRange &vertex_to_point)
bool	getPoint2Vertex (IndexRange &point_to_vertex)
bool	getFacetVertices (std::vector< IndexRange > &facet_vertices) const
bool	getFacetHalfSpaces (std::vector< HalfSpace > &facet_halfspaces)
Check hull services
bool	check ()
bool	checkFacets ()
bool	checkHull ()
Interface
void	selfDisplay (std::ostream &out) const
bool	isValid () const

Data Fields
public datas
Kernel	kernel
	Kernel that is duplicated for computing facet geometries.
int	debug_level
	debug_level from 0:no to 2
std::vector< Point >	points
	the set of points, indexed as in the array.
IndexRange	input2comp
IndexRange	comp2input
IndexRange	processed_points
	Points already processed (and within the convex hull).
std::vector< Index >	assignment
	assignment of points to facets
std::vector< Facet >	facets
	the current set of facets.
std::set< Index >	deleted_facets
	set of deleted facets
IndexRange	p2v
	point index -> vertex index (or UNASSIGNED)
IndexRange	v2p
	vertex index -> point index
Size	nb_finite_facets
	Number of finite facets.
Size	nb_infinite_facets
	Number of infinite facets (!= 0 only for specific kernels)
std::vector< double >	timings
	Timings of the different phases: 0: init, 1: facets, 2: vertices.
std::vector< Size >	facet_counter
	Counts the number of facets with a given number of vertices.

Static Public Attributes
static const Size	dimension = Point::dimension

Protected Attributes
protected datas
Status	myStatus

protected services
InternalScalar	height (const Facet &F, const Point &p) const
bool	above (const Facet &F, const Point &p) const
bool	aboveOrOn (const Facet &F, const Point &p) const
bool	on (const Facet &F, const Point &p) const
void	cleanFacets ()
void	renumberInfiniteFacets ()
bool	processFacet (std::queue< Index > &Q)
bool	checkFacet (Index f) const
Index	newFacet ()
void	deleteFacet (Index f)
void	makeNeighbors (const Index if1, const Index if2)
void	unmakeNeighbors (const Index if1, const Index if2)
Index	mergeParallelFacets (const Index if1, const Index if2)
bool	areFacetsParallel (const Index if1, const Index if2) const
bool	areFacetsNeighbor (const Index if1, const Index if2) const
Facet	makeFacet (const IndexRange &base, Index below) const
IndexRange	pointsOnRidge (const Ridge &R) const
IndexRange	orientedFacetPoints (Index f) const
void	filterVerticesOnFacet (const Index f)
IndexRange	pickInitialSimplex () const
bool	computeSimplexConfiguration (const IndexRange &full_simplex)
static IndexRange	pickIntegers (const Size d, const Size n)

Detailed Description

template<typename TKernel>
struct DGtal::QuickHull< TKernel >

Aim: Implements the quickhull algorithm by Barber et al. [9], a famous arbitrary dimensional convex hull computation algorithm. It relies on dedicated geometric kernels for computing and comparing facet geometries.

Description of template class 'QuickHull'

You can use it to build convex hulls of points with integral coordinate (using kernel ConvexHullIntegralKernel), convex hulls with rational coordinates (using kernel ConvexHullRationalKernel) or to build Delaunay triangulations (using kernels DelaunayIntegralKernel and DelaunayRationalKernel).

See also

QuickHull algorithm in arbitrary dimension for convex hull and Delaunay cell complex computation

Below is a complete example that computes the convex hull of points randomly defined in a ball, builds a 3D mesh out of it and output it as an OBJ file.

#include "DGtal/base/Common.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/io/writers/SurfaceMeshWriter.h"
#include "QuickHull.h"
 
using namespace DGtal::Z3i;
int main( int argc, char* argv[] )
{
  int nb = argc > 1 ? atoi( argv[ 1 ] ) : 100; // nb points
  int R  = argc > 2 ? atoi( argv[ 2 ] ) : 10;  // x-radius of ellipsoid
  // (0) typedefs
  typedef DGtal::ConvexHullIntegralKernel< 3 > Kernel3D;
  typedef DGtal::QuickHull< Kernel3D >         QuickHull3D;
  // (1) create range of random points in ball
  std::vector< Point > V;
  const auto R2 = R
  for ( int i = 0; i < nb; ) {
    Point p( rand() % (2*R+1) - R, rand() % (2*R+1) - R, rand() % (2*R+1) - R );
    if ( p.squaredNorm() < R2 ) { V.push_back( p ); i++; }
  }
  // (2) compute convex hull
  QuickHull3D hull;
  hull.setInput( V );
  hull.computeConvexHull();
  std::cout << "#points=" << hull.nbPoints()
            << " #vertices=" << hull.nbVertices()
            << " #facets=" << hull.nbFacets() << std::endl;
  // (3) build mesh
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< std::vector< std::size_t > > facets;
  hull.getFacetVertices( facets );
  typedef DGtal::SurfaceMesh< RealPoint, RealVector> SMesh;
  SMesh mesh( positions.cbegin(), positions.cend(), facets.cbegin(), facets.cend() );
  // (4) output result as OBJ file
  std::ofstream out( "qhull.obj" );
  DGtal::SurfaceMeshWriter< Space::RealPoint, Space::RealVector >
    ::writeOBJ( out, mesh );
  out.close();
  return 0;
} 

Note

In opposition with the usual QuickHull implementation, this class uses a kernel that can be chosen in order to provide exact computations. This is the case for lattice points.

In opposition with CGAL 3D convex hull package, or with the arbitrary dimensional dD Triangulation package, this algorithm does not build a simplicial convex hull. Facets may not be trangles or simplices in higher dimensions.

This version is generally more than twice faster than CGAL convex_hull_3 for the usual CGAL kernels Cartesian and Exact_predicates_inexact_constructions_kernel.

However this implementation is not tailored for incremental dynamic convex hull computations.

Template Parameters

TKernel

any type of QuickHull kernel, like ConvexHullIntegralKernel.

Examples

geometry/tools/checkLatticeBallQuickHull.cpp, and geometry/tools/exampleLatticeBallDelaunay2D.cpp.

Definition at line 139 of file QuickHull.h.

Member Typedef Documentation

CombinatorialPlaneSimplex

template<typename TKernel>

typedef Kernel::CombinatorialPlaneSimplex DGtal::QuickHull< TKernel >::CombinatorialPlaneSimplex

Definition at line 151 of file QuickHull.h.

HalfSpace

template<typename TKernel>

typedef Kernel::HalfSpace DGtal::QuickHull< TKernel >::HalfSpace

Definition at line 150 of file QuickHull.h.

Index

template<typename TKernel>

typedef std::size_t DGtal::QuickHull< TKernel >::Index

Definition at line 146 of file QuickHull.h.

IndexRange

template<typename TKernel>

typedef std::vector< Index > DGtal::QuickHull< TKernel >::IndexRange

Definition at line 149 of file QuickHull.h.

InternalScalar

template<typename TKernel>

typedef Kernel::InternalScalar DGtal::QuickHull< TKernel >::InternalScalar

Definition at line 145 of file QuickHull.h.

Kernel

template<typename TKernel>

typedef TKernel DGtal::QuickHull< TKernel >::Kernel

Definition at line 141 of file QuickHull.h.

Point

template<typename TKernel>

typedef Kernel::CoordinatePoint DGtal::QuickHull< TKernel >::Point

Definition at line 142 of file QuickHull.h.

Ridge

template<typename TKernel>

typedef std::pair< Index, Index > DGtal::QuickHull< TKernel >::Ridge

A ridge for point p is a pair of facets, such that p is visible from first facet, but not from second facet.

Definition at line 236 of file QuickHull.h.

Scalar

template<typename TKernel>

typedef Kernel::CoordinateScalar DGtal::QuickHull< TKernel >::Scalar

Definition at line 144 of file QuickHull.h.

Size

template<typename TKernel>

typedef std::size_t DGtal::QuickHull< TKernel >::Size

Definition at line 147 of file QuickHull.h.

Vector

template<typename TKernel>

typedef Kernel::CoordinateVector DGtal::QuickHull< TKernel >::Vector

Definition at line 143 of file QuickHull.h.

Member Enumeration Documentation

anonymous enum

template<typename TKernel>

anonymous enum

Label for points that are not assigned to any facet.

Enumerator
UNASSIGNED

Definition at line 155 of file QuickHull.h.

{ UNASSIGNED = (Index) -1 };

Status

template<typename TKernel>

enum class DGtal::QuickHull::Status

strong

Represents the status of a QuickHull object.

Enumerator
Uninitialized	QuickHull is empty and has just been instantiated.
InputInitialized	A range of input points has been given to QuickHull.
SimplexCompleted	An initial full-dimensional simplex has been found. QuickHull core algorithm can start.
FacetsCompleted	All facets of the convex hull are identified.
VerticesCompleted	All vertices of the convex hull are determined.
AllCompleted	Same as VerticesCompleted.
NotFullDimensional	Error: the initial simplex is not full dimensional.
InvalidRidge	Error: some ridge is not consistent (probably integer overflow).
InvalidConvexHull	Error: the convex hull does not contain all its points (probably integer overflow).

Definition at line 239 of file QuickHull.h.

                      {
      Uninitialized     =  0, 
      InputInitialized  =  1, 
      SimplexCompleted  =  2, 
      FacetsCompleted   =  3, 
      VerticesCompleted =  4, 
      AllCompleted      =  5, 
      NotFullDimensional= 10, 
      InvalidRidge      = 11, 
      InvalidConvexHull = 12  
    };

Constructor & Destructor Documentation

QuickHull()

template<typename TKernel>

DGtal::QuickHull< TKernel >::QuickHull ( const Kernel & K = Kernel(), int dbg = 0 )

inline

Default constructor

Parameters

[in]	K	a kernel for computing facet geometries.
[in]	dbg	the trace level, from 0 (no) to 3 (very verbose).

Definition at line 259 of file QuickHull.h.

      : kernel( K ), debug_level( dbg ), myStatus( Status::Uninitialized )
    {}

Member Function Documentation

above()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::above ( const Facet & F, const Point & p ) const

inlineprotected

Parameters

F	any valid facet
p	any point

Returns

'true' iff p is above F.

Definition at line 856 of file QuickHull.h.

    { return kernel.above( F.H, p ); }

Referenced by DGtal::QuickHull< Kernel3D >::checkFacet(), DGtal::QuickHull< Kernel3D >::checkHull(), DGtal::QuickHull< Kernel3D >::computeSimplexConfiguration(), and DGtal::QuickHull< Kernel3D >::processFacet().

aboveOrOn()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::aboveOrOn ( const Facet & F, const Point & p ) const

inlineprotected

Parameters

F	any valid facet
p	any point

Returns

'true' iff p is above F or on F.

Definition at line 862 of file QuickHull.h.

    { return kernel.aboveOrOn( F.H, p ); }

Referenced by DGtal::QuickHull< Kernel3D >::checkFacet(), and DGtal::QuickHull< Kernel3D >::processFacet().

areFacetsNeighbor()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::areFacetsNeighbor ( const Index if1, const Index if2 ) const

inlineprotected

Checks if two facets are neighbors by looking at the points on their boundary.

Parameters

if1	a valid facet index
if2	a valid facet index

Returns

'true' if the two facets have enough common points to be direct neighbors.

Definition at line 1267 of file QuickHull.h.

    {
      const Facet& f1 = facets[ if1 ];
      const Facet& f2 = facets[ if2 ];
      Index nb = 0;
      for ( Index i1 = 0, i2 = 0; i1 < f1.on_set.size() && i2 < f2.on_set.size(); )
        {
          if ( f1.on_set[ i1 ] == f2.on_set[ i2 ] )     { nb++; i1++; i2++; }
          else if ( f1.on_set[ i1 ] < f2.on_set[ i2 ] ) i1++;
          else                                          i2++;
        }
      return nb >= ( dimension - 1 );
    }

Referenced by DGtal::QuickHull< Kernel3D >::checkFacet(), and DGtal::QuickHull< Kernel3D >::processFacet().

areFacetsParallel()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::areFacetsParallel ( const Index if1, const Index if2 ) const

inlineprotected

Checks if two distinct facets are parallel (i.e. should be merged).

Parameters

if1	a valid facet index
if2	a valid facet index

Returns

'true' if the two facets are parallel.

Definition at line 1248 of file QuickHull.h.

    {
      ASSERT( if1 != if2 );
      const Facet& f1 = facets[ if1 ];
      const Facet& f2 = facets[ if2 ];
      if ( kernel.equal( f1.H, f2.H ) ) return true;
      // Need to check if one N is a multiple of the other.
      for ( auto&& v : f1.on_set )
        if ( ! on( f2, points[ v ] ) ) return false;
      return true;
    }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

BOOST_STATIC_ASSERT()

template<typename TKernel>

DGtal::QuickHull< TKernel >::BOOST_STATIC_ASSERT

(

(Point::dimension==Vector::dimension)

)

check()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::check ( )

inline

Global validity check of the convex hull after processing.

Note

Be careful, this function is much slower than computing the convex hull. It can take a long time since its complexity is \( O(n f ) \), where n is the number of input points and f the number of facets.

Definition at line 715 of file QuickHull.h.

    {
      bool ok = true;
      if ( status() < Status::FacetsCompleted || status() > Status::AllCompleted ) {
        trace.warning() << "[Quickhull::check] invalid status="
                        << (int)status() << std::endl;
        ok = false;
      }
      if ( processed_points.size() != points.size() ) {
        trace.warning() << "[Quickhull::check] not all points processed: "
                        << processed_points.size() << " / " <<  points.size()
                        << std::endl;
        ok = false;
      }
      if ( ! checkHull() ) {
        trace.warning() << "[Quickhull::check] Check hull is invalid. "
                        << std::endl;
        ok = false;
      }
      if ( ! checkFacets() ) {
        trace.warning() << "[Quickhull::check] Check facets is invalid. "
                        << std::endl;
        ok = false;
      }
      return ok;
    }

checkFacet()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::checkFacet ( Index f ) const

inlineprotected

Returns

true if the facet is valid

Definition at line 1136 of file QuickHull.h.

    {
      const Facet& F = facets[ f ];
      bool ok = F.on_set.size() >= dimension;
      for ( auto v : F.on_set )
        if ( ! on( F, points[ v ] ) ) {
          trace.error() << "[QuickHull::checkFacet( " << f
                    << ") Invalid 'on' vertex " << v << std::endl;
          ok = false;
        }
      if ( F.neighbors.size() < dimension ) {
        trace.error() << "[QuickHull::checkFacet( " << f
                  << ") Not enough neighbors " << F.neighbors.size() << std::endl;
        ok = false;
      }
      for ( auto nf : F.neighbors )
        if ( ! areFacetsNeighbor( f, nf ) ) {
          trace.error() << "[QuickHull::checkFacet( " << f
                    << ") Invalid neighbor " << nf << std::endl;
          ok = false;
        }
      if ( aboveOrOn( F, points[ F.below ] ) ) {
        trace.error() << "[QuickHull::checkFacet( " << f
                  << ") Bad below point " << F.below << std::endl;
        ok = false;
      }
      for ( auto ov : F.outside_set )
        if ( ! above( F, points[ ov ] ) ) {
          trace.error() << "[QuickHull::checkFacet( " << f
                    << ") Bad outside point " << ov << std::endl;
          ok = false;
        }
      for ( auto && v : F.on_set ) {
        Size      n = 0;
        for ( auto&& N : facets[ f ].neighbors )
          if ( on( facets[ N ], points[ v ] ) ) n += 1;
        if ( n < dimension-1 ) {
          trace.error() << "[QuickHull::checkFacet( " << f << ") 'on' point " << v
                    << " is a vertex of " << n << " facets "
                    << "(should be >= " << dimension-1 << ")" << std::endl;
          ok = false;
        }
      }
      return ok;
    }

Referenced by DGtal::QuickHull< Kernel3D >::checkFacets(), and DGtal::QuickHull< Kernel3D >::processFacet().

checkFacets()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::checkFacets ( )

inline

Returns

'true' if all facets are consistent with their neighors

Definition at line 743 of file QuickHull.h.

    {
      Size   nb = 0;
      Size nbok = 0;
      for ( Index f = 0; f < facets.size(); ++f )
        if ( deleted_facets.count( f ) == 0 ) {
          bool ok = checkFacet( f );
          nbok   += ok ? 1 : 0;
          nb     += 1;
        }
      return nb == nbok;
    }

Referenced by DGtal::QuickHull< Kernel3D >::check(), DGtal::QuickHull< Kernel3D >::computeSimplexConfiguration(), and DGtal::QuickHull< Kernel3D >::processFacet().

checkHull()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::checkHull ( )

inline

Returns

'true' iff the current set of facets encloses all the points

Note

Be careful, this function is much slower than computing the convex hull. It can take a long time since its complexity is \( O(n f ) \), where n is the number of input points and f the number of facets.

Definition at line 762 of file QuickHull.h.

    {
      Index nbok = 0;
      // for ( Index v = 0; v < points.size(); ++v ) {
      for ( auto v : processed_points ) {
        bool ok = true;
        for ( Index f = 0; f < facets.size(); ++f )
          if ( deleted_facets.count( f ) == 0 ) {
            if ( above( facets[ f ], points[ v ] ) ) {
              ok = false;
              trace.error() << "- bad vertex " << v << " " << points[ v ]
                            << " dist="
                            << ( kernel.height( facets[ f ].H, points[ v ] ) )
                            << std::endl;
              break;
            }
          }
        nbok += ok ? 1 : 0;
      }
      if ( debug_level >= 2 ) {
        trace.info() << nbok << "/"
                  << processed_points.size() << " vertices inside convex hull."
                  << std::endl;
      }
      if ( nbok != processed_points.size() ) myStatus = Status::InvalidConvexHull;
      return nbok == processed_points.size();
    }

Referenced by DGtal::QuickHull< Kernel3D >::check(), DGtal::QuickHull< Kernel3D >::computeSimplexConfiguration(), and DGtal::QuickHull< Kernel3D >::processFacet().

cleanFacets()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::cleanFacets ( )

inlineprotected

Cleans and renumber the facets so that no one belongs to deleted_facets.

Definition at line 873 of file QuickHull.h.

    {
      if ( deleted_facets.empty() ) return;
      IndexRange renumbering( facets.size() );
      Index i = 0;
      Index j = 0;
      for ( auto& l : renumbering ) {
        if ( ! deleted_facets.count( j ) ) l = i++;
        else l = UNASSIGNED;
        j++;
      }
      const Index nf = facets.size() - deleted_facets.size();
      deleted_facets.clear();
      for ( Index f = 0; f < facets.size(); f++ )
        if ( ( renumbering[ f ] != UNASSIGNED ) && ( f != renumbering[ f ] ) )
          facets[ renumbering[ f ] ] = facets[ f ];
      facets.resize( nf );
      for ( auto& F : facets ) {
        for ( auto& N : F.neighbors ) {
          if ( renumbering[ N ] == UNASSIGNED )
            trace.error() << "Invalid deleted neighboring facet." << std::endl;
          else N = renumbering[ N ];
        }
      }
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeFacets().

clear()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::clear ( )

inline

Clears the object.

Definition at line 270 of file QuickHull.h.

    {
      myStatus = Status::Uninitialized;
      points.clear();
      processed_points.clear();
      input2comp.clear();
      comp2input.clear();
      assignment.clear();
      facets.clear();
      deleted_facets.clear();
      p2v.clear();
      v2p.clear();
      timings.clear();
    }

Referenced by DGtal::QuickHull< Kernel3D >::setInput().

computeConvexHull()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::computeConvexHull ( Status target = Status::VerticesCompleted )

inline

Computes the convex hull of the given range of points until a specified target:

• Status::SimplexCompleted: an initial full dimensional simplex is determined.
• Status::FacetsCompleted: all the facets of the convex hull are extracted (core of quickhull). You can stop here if you only need an H-representation of the output polytope.
• Status::VerticesCompleted: all the vertices of the hull are determined and are consistently ordered on each facet. You need to stop here if you need a V-representation of the output polytope.

Precondition

status() must be at least Status::InputInitialized

Parameters

[in]

target

the computation target in Status::SimplexCompleted, Status::FacetsCompleted, Status::VerticesCompleted.

Returns

'true' if the computation target has been successfully achieved, 'false' if the achieved status is not the one specified.

Definition at line 460 of file QuickHull.h.

    {
      if ( target < Status::InputInitialized || target > Status::AllCompleted )
        return false;
      Clock tic;
      if ( status() == Status::InputInitialized )
        { // Initialization
          tic.startClock();
          bool ok1 = computeInitialSimplex();
          timings.push_back( tic.stopClock() );
          if ( ! ok1 )              return false;
          if ( status() == target ) return true;
        }
      if ( status() == Status::SimplexCompleted )
        { // Computes facets
          tic.startClock();
          bool ok2 = computeFacets();
          timings.push_back( tic.stopClock() );
          if ( ! ok2 )              return false;
          if ( status() == target ) return true;
        }
      if ( status() == Status::FacetsCompleted )
        { // Computes vertices
          tic.startClock();
          bool ok3 = computeVertices();
          timings.push_back( tic.stopClock() );
          if ( ! ok3 )              return false;
          if ( status() == target ) return true;
        }
      if ( target == Status::AllCompleted
           && status() == Status::VerticesCompleted )
        { // for now, Status::VerticesCompleted and
          // Status::AllCompleted are the same.
          myStatus = Status::AllCompleted;
          return true;
        }
      return false;
    }

Referenced by main().

computeFacets()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::computeFacets ( )

inline

Computes the facets of the convex hull using Quickhull algorithm. If everything went well, the status is Status::FacetsCompleted afterwards.

Precondition

the status shoud be Status::SimplexCompleted (computeInitialSimplex should have been called).

Returns

'true' except if the status is not Initialized when called.

Definition at line 523 of file QuickHull.h.

    {
      if ( status() != Status::SimplexCompleted ) return false;
      std::queue< Index > Q;
      for ( Index fi = 0; fi < facets.size(); ++fi )
        Q.push( fi );
      Index n = 0;
      while ( processFacet( Q ) ) {
        if ( debug_level >= 1 )
          trace.info() << "---- Iteration " << n++ << " #Q=" << Q.size() << std::endl;
      }
      cleanFacets();
      if ( debug_level >= 2 ) {
        trace.info() << ".... #facets=" << facets.size()
                  << " #deleted=" << deleted_facets.size() << std::endl;
      }
      myStatus = Status::FacetsCompleted;
      return true;
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeConvexHull().

computeInitialSimplex()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::computeInitialSimplex ( )

inline

Computes the initial full dimensional simplex from the input data.

Returns

'true' iff the input data contains d+1 points in general position, the object has then the status SimplexCompleted, otherwise returns 'false' and the status is NotFullDimensional.

Definition at line 504 of file QuickHull.h.

    {
      const auto full_simplex = pickInitialSimplex();
      if ( full_simplex.empty() ) {
        myStatus = Status::NotFullDimensional;
        return false;
      }
      return computeSimplexConfiguration( full_simplex );
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeConvexHull().

computeSimplexConfiguration()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::computeSimplexConfiguration ( const IndexRange & full_simplex )

inlineprotected

Computes the initial configuration induced by the given full dimensional simplex. Follows Status::InputInitialized and and terminate with Status::SimplexCompleted if everything went well.

Returns

'true' iff the input data contains d+1 points in general position, the object has then the status SimplexCompleted, otherwise returns 'false' and the status is NotFullDimensional.

Definition at line 1429 of file QuickHull.h.

    {
      assignment = std::vector< Index >( points.size(), UNASSIGNED );
      facets.resize( dimension + 1 );
      deleted_facets.clear();
      for ( Index j = 0; j < full_simplex.size(); ++j )
        {
          IndexRange lsimplex( dimension );
          IndexRange isimplex( dimension );
          Index s = 0;
          for ( Index i = 0; i <= dimension; i++ )
            if ( i != j ) {
              lsimplex[ s ] = i;
              isimplex[ s ] = full_simplex[ i ];
              s++;
            }
          facets[ j ] = makeFacet( isimplex, full_simplex[ j ] );
          facets[ j ].neighbors = lsimplex;
          for ( auto&& v : isimplex ) facets[ j ].on_set.push_back( v );
          std::sort( facets[ j ].on_set.begin(), facets[ j ].on_set.end() );
        }
      // List of unassigned vertices
      for ( Index fi = 0; fi < facets.size(); ++fi ) {
        Facet& f = facets[ fi ];
        for ( Index v = 0; v < points.size(); v++ )
          {
            if ( assignment[ v ] == UNASSIGNED && above( f, points[ v ] ) ) {
              f.outside_set.push_back( v );
              assignment[ v ] = fi;
            } 
          }
      }
      for ( Index v = 0; v < points.size(); v++ )
        if ( assignment[ v ] == UNASSIGNED )
          processed_points.push_back( v );
      
      // Display some information
      if ( debug_level >= 2 ) {
        for ( auto&& f : facets ) f.display( trace.info() );
      }
      myStatus = Status::SimplexCompleted;
      if ( debug_level >= 1 ) {
        trace.info() << "[CHECK INVARIANT] " << processed_points.size()
                     << " / " << points.size() << " points processed." << std::endl;
        bool okh = checkHull(); 
        if ( ! okh )
          trace.error() << "[computeInitialSimplex] Invalid convex hull" << std::endl;
        bool okf = checkFacets();
        if ( ! okf )
          trace.error() << "[computeInitialSimplex] Invalid facets" << std::endl;
        if ( ! ( okh && okf ) ) myStatus = Status::InvalidConvexHull;
      }
      return status() == Status::SimplexCompleted;
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeInitialSimplex(), and DGtal::QuickHull< Kernel3D >::setInitialSimplex().

computeVertices()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::computeVertices ( )

inline

Computes the vertices of the convex hull once the facets have been computed. It computes for each facet its vertices and reorder them so that, taken in order, their orientation matches the orientation of the facet. If everything went well, the status is Status::VerticesCompleted afterwards.

Precondition

the status shoud be Status::FacetsCompleted (computeFacets should have been called).

Returns

'true' except if the status is not FacetsCompleted when called.

Definition at line 553 of file QuickHull.h.

    {
      static const int MAX_NB_VPF = 10 * dimension;
      if ( status() != Status::FacetsCompleted ) return false;
 
      // Renumber infinite facets in case of Delaunay triangulation computation.
      renumberInfiniteFacets();
      
      // Builds the maps v2p: vertex -> point, and p2v : point -> vertex.
      facet_counter = IndexRange( MAX_NB_VPF, 0 );
      v2p.clear();
      p2v.resize( points.size() );
      std::vector< IndexRange > p2f( points.size() );
      for ( Index f = 0; f < facets.size(); ++f ) {
        for ( auto&& p : facets[ f ].on_set ) p2f[ p ].push_back( f ); 
      }
 
      // vertices belong to at least d facets
      Index v = 0;
      for ( Index p = 0; p < points.size(); ++p ) {
        const auto nbf = p2f[ p ].size();
        facet_counter[ std::min( (int) nbf, MAX_NB_VPF-1 ) ] += 1;
        if ( nbf >= dimension ) {
          v2p.push_back( p );
          p2v[ p ] = v++;
        }
        else p2v[ p ] = UNASSIGNED;
      }
 
      // Display debug informations
      if ( debug_level >= 1 ) {
        trace.info() << "#vertices=" << v2p.size() << " #facets=" << facets.size()
                  << std::endl;
        trace.info() << "#inc_facets/point= ";
        for ( auto n : facet_counter ) trace.info() << n << " ";
        trace.info() << std::endl;
      }
      myStatus = Status::VerticesCompleted;
      return true;
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeConvexHull().

deleteFacet()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::deleteFacet ( Index f )

inlineprotected

Deletes the given facet f.

Parameters

f	a valid facet index

Definition at line 1193 of file QuickHull.h.

    {
      for ( auto n : facets[ f ].neighbors )
        facets[ n ].subNeighbor( f );
      deleted_facets.insert( f );
      facets[ f ].clear();
    }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

filterVerticesOnFacet()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::filterVerticesOnFacet ( const Index f )

inlineprotected

Filters each vertex on the facet f to keep only the ones that are on or below the neighboring facets

Note

intended for debugging purposes.

Parameters

[in]

f

any valid facet index.

Definition at line 1352 of file QuickHull.h.

    {
      auto & on_set = facets[ f ].on_set;
      for ( Index i = 0; i < on_set.size(); )
        {
          Index     v = on_set[ i ];
          Size      n = 0;
          for ( auto&& N : facets[ f ].neighbors )
            if ( on( facets[ N ], points[ v ] ) ) n += 1;
          if ( n >= dimension-1 ) i++;
          else {
            on_set[ i ] = on_set.back();
            on_set.pop_back();
          }                        
        }
      std::sort( on_set.begin(), on_set.end() );
    }

getFacetHalfSpaces()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::getFacetHalfSpaces ( std::vector< HalfSpace > & facet_halfspaces )

inline

This methods return the halfspaces corresponding to each facet. It is useful to build a polytope.

Parameters

[out]

facet_halfspaces

the range giving for each facet its halfspace as a pair (normal, intercept).

Returns

'true' if the convex hull was computed before and status() was Status::FacetsCompleted, Status::VerticesCompleted or Status::AllCompleted, 'false' otherwise.

Definition at line 690 of file QuickHull.h.

    {
      facet_halfspaces.clear();
      if ( ! ( status() >= Status::FacetsCompleted
               && status() <= Status::AllCompleted ) ) return false;
      facet_halfspaces.reserve( nbFacets() );
      for ( Index f = 0; f < nbFacets(); ++f )  {
        facet_halfspaces.push_back( facets[ f ].H );
      }
      return true;
    }

getFacetVertices()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::getFacetVertices ( std::vector< IndexRange > & facet_vertices ) const

inline

Parameters

[out]

facet_vertices

the range giving for each facet the indices of its vertices.

Returns

'true' if the convex hull was computed before and status() was Status::VerticesCompleted or Status::AllCompleted, 'false' otherwise.

Definition at line 666 of file QuickHull.h.

    {
      facet_vertices.clear();
      if ( ! ( status() >= Status::VerticesCompleted
               && status() <= Status::AllCompleted ) ) return false;
      facet_vertices.reserve( nbFacets() );
      for ( Index f = 0; f < nbFacets(); ++f )  {
        IndexRange ofacet = orientedFacetPoints( f );
        for ( auto& v : ofacet ) v = p2v[ v ];
        facet_vertices.push_back( ofacet );
      }
      return true;
    }

Referenced by main().

getPoint2Vertex()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::getPoint2Vertex ( IndexRange & point_to_vertex )

inline

Gets for each point its index in the output range of vertices (or UNASSIGNED if the point is not part of the convex hull)

Parameters

[out]

point_to_vertex

the range giving for each point its index in the output range of vertices (or UNASSIGNED if the point is not part of the convex hull)

Returns

'true' if the convex hull was computed before and status() was Status::VerticesCompleted or Status::AllCompleted, 'false' otherwise.

Definition at line 651 of file QuickHull.h.

    {
      point_to_vertex.clear();
      if ( ! ( status() >= Status::VerticesCompleted
               && status() <= Status::AllCompleted ) ) return false;
      point_to_vertex = p2v;
      return true;
    }

getVertex2Point()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::getVertex2Point ( IndexRange & vertex_to_point )

inline

Gets for each vertex its index in the input range of points.

Parameters

[out]

vertex_to_point

the range giving for each vertex its index in the input range of points.

Returns

'true' if the convex hull was computed before and status() was Status::VerticesCompleted or Status::AllCompleted, 'false' otherwise.

Definition at line 632 of file QuickHull.h.

    {
      vertex_to_point.clear();
      if ( ! ( status() >= Status::VerticesCompleted
               && status() <= Status::AllCompleted ) ) return false;
      vertex_to_point = v2p;
      return true;
    }

getVertexPositions()

template<typename TKernel>

template<typename OutputPoint>

bool DGtal::QuickHull< TKernel >::getVertexPositions ( std::vector< OutputPoint > & vertex_positions )

inline

Gets the positions of the convex hull vertices.

Template Parameters

OutputPoint

a model of point such that the Point datatype is convertible to it.

Parameters

[out]

vertex_positions

the range of vertex positions.

Returns

'true' if the convex hull was computed before and status() was Status::VerticesCompleted or Status::AllCompleted, 'false' otherwise.

Definition at line 612 of file QuickHull.h.

    {
      vertex_positions.clear();
      if ( ! ( status() >= Status::VerticesCompleted
               && status() <= Status::AllCompleted ) ) return false;
      vertex_positions.resize( v2p.size() );
      for ( Index i = 0; i < v2p.size(); i++ ) {
        kernel.convertPointTo( points[ v2p[ i ] ], vertex_positions[ i ] );
      }
      return true;
    }

Referenced by main().

height()

template<typename TKernel>

InternalScalar DGtal::QuickHull< TKernel >::height ( const Facet & F, const Point & p ) const

inlineprotected

Parameters

F	any valid facet
p	any point

Returns

the height of p wrt F (0: on, >0: above ).

Definition at line 850 of file QuickHull.h.

    { return kernel.height( F.H, p ); }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

isValid()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::isValid ( ) const

inline

Checks the validity/consistency of the object.

Returns

'true' if the object is valid, 'false' otherwise.

Definition at line 1515 of file QuickHull.h.

    {
      return status() >= Status::Uninitialized
        && status() <= Status::AllCompleted;
    }

makeFacet()

template<typename TKernel>

Facet DGtal::QuickHull< TKernel >::makeFacet ( const IndexRange & base, Index below ) const

inlineprotected

Builds a facet from a base convex set of at least d different points and a point below.

Parameters

[in]	base	a range containing at least d distinct points in general position.
[in]	below	a point below the hyperplane containing simplex.

Returns

a facet object of normal and c parameter such as the points of simplex lies on the facet hyperplane while point below lies under.

Definition at line 1292 of file QuickHull.h.

    {
      CombinatorialPlaneSimplex simplex;
      for ( Size i = 0; i < dimension; i++ ) simplex[ i ] = base[ i ];
      auto plane = kernel.compute( points, simplex, below );
      return Facet( plane, below );
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeSimplexConfiguration(), and DGtal::QuickHull< Kernel3D >::processFacet().

makeNeighbors()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::makeNeighbors ( const Index if1, const Index if2 )

inlineprotected

Makes two distinct facets if1 and if2 as neighbors

Parameters

if1	a valid facet index
if2	a valid facet index

Definition at line 1204 of file QuickHull.h.

    {
      facets[ if1 ].addNeighbor( if2 );
      facets[ if2 ].addNeighbor( if1 );
    }

Referenced by DGtal::QuickHull< Kernel3D >::mergeParallelFacets(), and DGtal::QuickHull< Kernel3D >::processFacet().

memory()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::memory ( ) const

inline

Returns

an estimation of the current memory occupation of this object.

Definition at line 287 of file QuickHull.h.

    {
      // int debug_level;
      Size M = sizeof( kernel ) + sizeof( int );
      // std::vector< Point > points;    
      M += sizeof( std::vector< Point > )
        + points.capacity() * sizeof( Point );
      M += sizeof( std::vector< Index > )
        + processed_points.capacity() * sizeof( Index );
      M += sizeof( std::vector< Index > )
        + input2comp.capacity() * sizeof( Index );
      M += sizeof( std::vector< Index > )
        + comp2input.capacity() * sizeof( Index );
      // std::vector< Index > assignment;
      M += sizeof( std::vector< Index > )
        + assignment.capacity() * sizeof( Index );
      // std::vector< Facet > facets;
      M += sizeof( std::vector< Facet > )
        + facets.capacity() * sizeof( Facet );
      for ( const auto& f : facets ) M += f.variableMemory();
      // std::set< Index > deleted_facets;
      M += sizeof( std::set< Index > )
        + deleted_facets.size() * ( sizeof( Index ) + 2*sizeof(Index*) );
      // IndexRange p2v;
      M += sizeof( std::vector< Index > )
        + p2v.capacity() * sizeof( Index );
      // IndexRange v2p;
      M += sizeof( std::vector< Index > )
        + v2p.capacity() * sizeof( Index );
      // std::vector< double > timings;
      M += sizeof( std::vector< double > )
        + timings.capacity() * sizeof( double );
      return M;
    }

mergeParallelFacets()

template<typename TKernel>

Index DGtal::QuickHull< TKernel >::mergeParallelFacets ( const Index if1, const Index if2 )

inlineprotected

Merge the two facets and return the index of the second one, which is the one deleted.

Parameters

if1	a valid facet index
if2	a valid facet index

Precondition

the two facets should be distinct and parallel

Definition at line 1225 of file QuickHull.h.

    {
      Facet& f1 = facets[ if1 ];
      Facet& f2 = facets[ if2 ];
      std::copy( f2.outside_set.cbegin(), f2.outside_set.cend(),
                 std::back_inserter( f1.outside_set ) );
      IndexRange merge_idx;
      std::set_union( f1.on_set.cbegin(), f1.on_set.cend(),
                      f2.on_set.cbegin(), f2.on_set.cend(),
                      std::back_inserter( merge_idx ) );
      f1.on_set.swap( merge_idx );
      for ( auto && nf2 : f2.neighbors ) {
        if ( nf2 == if1 ) continue;
        facets[ nf2 ].subNeighbor( if2 );
        makeNeighbors( if1, nf2 );
      }
      return if2;
    }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

nbFacets()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nbFacets ( ) const

inline

Returns

the number of facets of the convex hull.

Precondition

status() >= Status::FacetsCompleted

Definition at line 328 of file QuickHull.h.

    {
      ASSERT( status() >= Status::FacetsCompleted
              && status() <= Status::AllCompleted );
      return facets.size();
    }

Referenced by DGtal::QuickHull< Kernel3D >::getFacetHalfSpaces(), DGtal::QuickHull< Kernel3D >::getFacetVertices(), main(), and DGtal::QuickHull< Kernel3D >::selfDisplay().

nbFiniteFacets()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nbFiniteFacets ( ) const

inline

Returns

the number of finite facets of the convex hull.

Precondition

status() >= Status::VerticesCompleted

Definition at line 346 of file QuickHull.h.

    {
      ASSERT( status() >= Status::VerticesCompleted
              && status() <= Status::AllCompleted );
      return nb_finite_facets;
    }

nbInfiniteFacets()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nbInfiniteFacets ( ) const

inline

Returns

the number of infinite facets of the convex hull.

Precondition

status() >= Status::VerticesCompleted

Definition at line 355 of file QuickHull.h.

    {
      ASSERT( status() >= Status::VerticesCompleted
              && status() <= Status::AllCompleted );
      return nb_infinite_facets;
    }

nbPoints()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nbPoints ( ) const

inline

Returns

the number of points used as input.

Definition at line 323 of file QuickHull.h.

    { return points.size(); }

Referenced by main(), and DGtal::QuickHull< Kernel3D >::selfDisplay().

nbVertices()

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nbVertices ( ) const

inline

Returns

the number of vertices of the convex hull.

Precondition

status() >= Status::VerticesCompleted

Definition at line 337 of file QuickHull.h.

    {
      ASSERT( status() >= Status::VerticesCompleted
              && status() <= Status::AllCompleted );
      return v2p.size();
    }

Referenced by main(), and DGtal::QuickHull< Kernel3D >::selfDisplay().

newFacet()

template<typename TKernel>

Index DGtal::QuickHull< TKernel >::newFacet ( )

inlineprotected

Returns

an unused facet index

Definition at line 1183 of file QuickHull.h.

    {
      // SLightly faster to postpone deletion of intermediate facets.
      const Index f = facets.size();
      facets.push_back( Facet() );
      return f;
    }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

on()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::on ( const Facet & F, const Point & p ) const

inlineprotected

Parameters

F	any valid facet
p	any point

Returns

'true' iff p lies on F.

Definition at line 868 of file QuickHull.h.

    { return kernel.on( F.H, p ); }

Referenced by DGtal::QuickHull< Kernel3D >::areFacetsParallel(), DGtal::QuickHull< Kernel3D >::checkFacet(), and DGtal::QuickHull< Kernel3D >::filterVerticesOnFacet().

orientedFacetPoints()

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::orientedFacetPoints ( Index f ) const

inlineprotected

Given a facet index f, return its points oriented consistently with respect to the normal.

Parameters

f	any valid facet index

Returns

the range of point indices that support this facet, in a consistent ordering.

Note

the order of points is consistent if, picking any d-simplex in order in this range, their associated half-spaces have all the same orientation.

Definition at line 1323 of file QuickHull.h.

    {
      const Facet& F = facets[ f ];
      IndexRange result = F.on_set;
      // Sort a facet such that its points, taken in order, have
      // always the same orientation of the facet.  More precisely,
      // facets span a `dimension-1` vector space. There are thus
      // dimension-2 fixed points, and the last ones (at least two)
      // may be reordered.
      CombinatorialPlaneSimplex splx;
      for ( Dimension k = 0; k < dimension-2; k++ )
        splx[ k ] = result[ k ];
      std::sort( result.begin()+dimension-2, result.end(),
                 [&] ( Index i, Index j )
                 {
                   splx[ dimension-2 ] = i;
                   splx[ dimension-1 ] = j;
                   const auto H = kernel.compute( points, splx );
                   const auto orient = kernel.dot( F.H, H );
                   return orient > 0;
                 } );
      return result;
    }

Referenced by DGtal::QuickHull< Kernel3D >::getFacetVertices().

pickInitialSimplex()

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::pickInitialSimplex ( ) const

inlineprotected

Returns

a full dimensional simplex as a vector of d + 1 distinct indices of input points, or an empty vector if none was found.

Definition at line 1372 of file QuickHull.h.

    {
      const Size          nb = points.size();
      if ( nb < dimension + 1 ) return IndexRange();
      IndexRange        best = pickIntegers( dimension + 1, nb );
      CombinatorialPlaneSimplex splx;
      for ( Index j = 0; j < dimension; ++j ) splx[ j ] = best[ j ];
      const auto     first_H = kernel.compute( points, splx, best.back() );
      auto       best_volume = kernel.volume ( first_H, points[ best.back() ] );
      const Size     nbtries = std::min( (Size) 10, 1 + nb / 10 );
      const Size max_nbtries = std::max( (Size) 10, 2 * nb );
      for ( Size i = 0; i < max_nbtries; i++ )
        {
          IndexRange tmp = pickIntegers( dimension + 1, nb );
          for ( Index j = 0; j < dimension; ++j ) splx[ j ] = tmp[ j ];
          const auto        tmp_H = kernel.compute( points, splx, tmp.back() );
          const auto   tmp_volume = kernel.volume ( tmp_H, points[ tmp.back() ] );
          if ( best_volume < tmp_volume ) {
            if ( debug_level >= 1 ) {
              trace.info() << "(" << i << ")"
                        << " new_volume = " << tmp_volume
                        << " > " << best_volume << std::endl;
            }
            best = tmp;
            best_volume = tmp_volume;
          }
          if ( i >= nbtries && best_volume > 0 ) return best;
        }
      return IndexRange();
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeInitialSimplex().

pickIntegers()

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::pickIntegers ( const Size d, const Size n )

inlinestaticprotected

Returns

a vector of d distinct integers in {0, 1, ..., n-1} randomly chosen.

Parameters

[in]	d	the number of returned integers
[in]	n	the range of possible integers {0, 1, ..., n-1}

Definition at line 1406 of file QuickHull.h.

    {
      IndexRange result( d );
      bool distinct = false;
      while ( ! distinct )
        {
          distinct = true;
          for ( Index i = 0; i < d; i++ ) result[ i ] = rand() % n;
          std::sort( result.begin(), result.end() );
          for ( Index i = 1; distinct && i < d; i++ )
            distinct = result[ i-1 ] != result[ i ];
        }
      return result;
    }

Referenced by DGtal::QuickHull< Kernel3D >::pickInitialSimplex().

pointsOnRidge()

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::pointsOnRidge ( const Ridge & R ) const

inlineprotected

Parameters

[in]

R

a ridge between two facets (a pair of facets).

Returns

the points lying on a ridge.

Definition at line 1302 of file QuickHull.h.

    {
      // Slightly faster to look at points instead at true geometry.
      IndexRange result;
      const Facet& f1 = facets[ R.first ];
      const Facet& f2 = facets[ R.second ];
      std::set_intersection( f1.on_set.cbegin(), f1.on_set.cend(),
                             f2.on_set.cbegin(), f2.on_set.cend(),
                             std::back_inserter( result ) );
      return result;
    }

Referenced by DGtal::QuickHull< Kernel3D >::processFacet().

processFacet()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::processFacet ( std::queue< Index > & Q )

inlineprotected

Process top facet in queue Q as in Quickhull algorithm, and updates the queue accordingly (top facet poped, new facets pushed).

Parameters

[in,out]

Q

a queue of facet index to process. which is updated by the method.

Returns

'true' if there is still work to do, 'false' when finished

Note

Core of Quickhull algorithm.

Definition at line 944 of file QuickHull.h.

    {
      // If Q empty, we are done
      if ( Q.empty() ) return false;
      Index F = Q.front();
      Q.pop();
      // If F is already deleted, proceed to next in queue.
      if ( deleted_facets.count( F ) ) return true;
      // Take car of current facet.
      const Facet& facet = facets[ F ];
      if ( debug_level >= 3 ) {
        trace.info() << "---------------------------------------------" << std::endl;
        trace.info() << "---- ACTIVE FACETS---------------------------" << std::endl;
        bool ok = true;
        for ( Index i = 0; i < facets.size(); i++ )
          if ( ! deleted_facets.count( i ) ) {
            trace.info() << "- facet " << i << " ";
            facets[ i ].display( trace.info() );
            ok = ok && checkFacet( i );
          }
        if ( ! ok ) { // should not happen.
          myStatus = Status::InvalidConvexHull;
          return false;
        }
      }
      if ( debug_level >= 2 ) {
        trace.info() << "---------------------------------------------" << std::endl;
        trace.info() << "Processing facet " << F << " ";
        facet.display( trace.info() );
      }
      if ( facet.outside_set.empty() ) return true;
      // Selects furthest vertex
      Index  furthest_v = facet.outside_set[ 0 ];
      auto   furthest_h = height( facet, points[ furthest_v ] );
      for ( Index v = 1; v < facet.outside_set.size(); v++ ) {
        auto h = height( facet, points[ facet.outside_set[ v ] ] );
        if ( h > furthest_h ) {
          furthest_h = h;
          furthest_v = facet.outside_set[ v ];
        }
      }
      const Point& p = points[ furthest_v ];
      // Extracts Visible facets V and Horizon Ridges H
      std::vector< Index > V;   // visible facets
      std::set< Index >    M;   // marked facets (are in E or were in E)
      std::queue< Index >  E;   // queue to extract visible facets
      std::vector< Ridge > H;   // visible facets
      E.push  ( F );
      M.insert( F );
      while ( ! E.empty() ) {
        Index G = E.front(); E.pop();
        V.push_back( G );
        for ( auto& N : facets[ G ].neighbors ) {
          if ( aboveOrOn( facets[ N ], p ) ) {
            if ( M.count( N ) ) continue;
            E.push( N );
          } else {
            H.push_back( { G, N } );
          }
          M.insert( N );
        }
      } // while ( ! E.empty() ) 
      if ( debug_level >= 1 ) {
        trace.info() << "#Visible=" << V.size() << " #Horizon=" << H.size()
                  << " furthest_v=" << furthest_v << std::endl;
      }
      // Create new facets
      IndexRange new_facets;
      // For each ridge R in H
      for ( Index i = 0; i < H.size(); i++ )
        {
          // Create a new facet from R and p
          IndexRange ridge = pointsOnRidge( H[ i ] );
          if ( debug_level >= 3 ) {
            trace.info() << "Ridge (" << H[i].first << "," << H[i].second << ") = {";
            for ( auto&& r : ridge ) trace.info() << " " << r;
            trace.info() << " } furthest_v=" << furthest_v << std::endl;
          }
          IndexRange base( 1 + ridge.size() );
          Index j     = 0;
          base[ j++ ] = furthest_v;
          for ( auto&& v : ridge ) base[ j++ ] = v;
          if ( j < dimension ) {
            trace.error() << "Bad ridge between " << std::endl
                      << "- facet " << H[i].first << " ";
            facets[ H[i].first ].display( trace.error() );
            trace.error() << "- facet " << H[i].second << " ";
            facets[ H[i].second ].display( trace.error() );
          }                    
          Index nf = newFacet();
          new_facets.push_back( nf );
          facets[ nf ] = makeFacet( base, facets[ H[i].first ].below );
          facets[ nf ].on_set = IndexRange { base.cbegin(), base.cend() };
          std::sort( facets[ nf ].on_set.begin(), facets[ nf ].on_set.end() );
          makeNeighbors( nf, H[ i ].second );
          if ( debug_level >= 3 ) {
            trace.info() << "* New facet " << nf << " ";
            facets[ nf ].display( trace.info() );
          }
          // Checks that the facet is not parallel to another in the Horizon
          for ( Index k = 0; k < new_facets.size() - 1; k++ )
            if ( areFacetsParallel( new_facets[ k ], nf ) ) {
              if ( debug_level >= 1 ) {
                trace.info() << "Facets " << new_facets[ k ] << " and " << nf
                          << " are parallel => merge." << std::endl;
              }
              mergeParallelFacets( new_facets[ k ], nf );
              new_facets.pop_back();
              deleteFacet( nf );
              if ( debug_level >= 3 ) {
                facets[ new_facets[ k ] ].display( trace.info() );
              }
            }
        }
      // For each new facet 
      for ( Index i = 0; i < new_facets.size(); i++ )
        { // link the new facet to its neighbors
          for ( Index j = i + 1; j < new_facets.size(); j++ )
            {
              const Index nfi = new_facets[ i ];
              const Index nfj = new_facets[ j ];
              if ( areFacetsNeighbor( nfi, nfj ) )
                makeNeighbors( nfi, nfj );
            }
        }
      // Extracts all outside points from visible facets V
      IndexRange outside_pts;
      for ( auto&& vf : V ) {
        for ( auto&& v : facets[ vf ].outside_set ) {
          if ( v != furthest_v ) {
            outside_pts.push_back( v );
            assignment[ v ] = UNASSIGNED;
          }
        }
      }
      // For each new facet F'
      for ( Index i = 0; i < new_facets.size(); i++ ) {
        Facet& Fp = facets[ new_facets[ i ] ];
        Index max_j = outside_pts.size();
        for ( Index j = 0; j < max_j; ) {
          const Index v = outside_pts[ j ];
          if ( above( Fp, points[ v ] ) ) {
            Fp.outside_set.push_back( v );
            assignment[ v ]  = new_facets[ i ];
            outside_pts[ j ] = outside_pts.back();
            outside_pts.pop_back();
            max_j--;
          } else j++;
        }
        if ( debug_level >= 3 ) {
          trace.info() << "- New facet " << new_facets[ i ] << " ";
          Fp.display( trace.info() );
        }
      }
      // Update processed points
      processed_points.push_back( furthest_v );
      for ( auto v : outside_pts ) processed_points.push_back( v );
      
      // Delete the facets in V
      for ( auto&& v : V ) {
        if ( debug_level >= 2 ) {
          trace.info() << "Delete facet " << v << " ";
          facets[ v ].display( trace.info() );
        }
        deleteFacet( v );
      }
 
      // Add new facets to queue
      for ( Index i = 0; i < new_facets.size(); i++ )
        Q.push( new_facets[ i ] );
      if ( debug_level >= 1 ) {
        trace.info() << "#facets=" << facets.size()
                  << " #deleted=" << deleted_facets.size() << std::endl;
      }
 
      // Checks that everything is ok.
      if ( debug_level >= 1 ) {
        trace.info() << "[CHECK INVARIANT] " << processed_points.size()
                     << " / " << points.size() << " points processed." << std::endl;
        bool okh = checkHull(); 
        if ( ! okh )
          trace.error() << "[computeFacet] Invalid convex hull" << std::endl;
        bool okf = checkFacets();
        if ( ! okf )
          trace.error() << "[computeFacet] Invalid facets" << std::endl;
        if ( ! ( okh && okf ) ) myStatus = Status::InvalidConvexHull;
      }
 
      return status() == Status::SimplexCompleted;
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeFacets().

renumberInfiniteFacets()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::renumberInfiniteFacets ( )

inlineprotected

Determine infinite facets and renumber them so that finite facets come first and infinite facets come after.

Definition at line 901 of file QuickHull.h.

    {
      nb_finite_facets   = facets.size();
      nb_infinite_facets = 0;
      if ( ! kernel.hasInfiniteFacets()  ) return;
      IndexRange renumbering( facets.size() );
      Index i = 0;
      Index k = facets.size();
      Index j = 0;
      for ( auto& l : renumbering ) {
        if ( ! kernel.isHalfSpaceFacetInfinite( facets[ j ].H ) ) l = i++;
        else l = --k;
        j++;
      }
      if ( i != k )
        trace.error() << "[Quickhull::renumberInfiniteFacets]"
                      << " Error renumbering infinite facets "
                      << " up finite=" << i << " low infinite=" << k << std::endl;
      std::vector< Facet > new_facets( facets.size() );
      for ( Index f = 0; f < facets.size(); f++ )
        new_facets[ renumbering[ f ] ].swap( facets[ f ] );
      facets.swap( new_facets );
      for ( auto& F : facets ) {
        for ( auto& N : F.neighbors ) {
          N = renumbering[ N ];
        }
      }
      // Assign correct number of facets.
      nb_finite_facets   = i;
      nb_infinite_facets = facets.size() - nb_finite_facets;
    }

Referenced by DGtal::QuickHull< Kernel3D >::computeVertices().

selfDisplay()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::selfDisplay ( std::ostream & out ) const

inline

Writes/Displays the object on an output stream.

Parameters

out	the output stream where the object is written.

Definition at line 1494 of file QuickHull.h.

    {
      out << "[QuickHull " << dimension << "D"
        //        << " status=" << status()
          << " #P=" << nbPoints();
      if ( status() >= Status::FacetsCompleted && status() <= Status::AllCompleted )
        out << " #F=" << nbFacets();
      if ( status() >= Status::VerticesCompleted && status() <= Status::AllCompleted )
        out << " #V=" << nbVertices();
      out << "]";
      if ( status() >= Status::FacetsCompleted && status() <= Status::AllCompleted
           && nbFacets() == 24 ) {
        for ( auto f : facets ) f.display( out );
        for ( auto v : v2p ) out << points[ v2p[ v ] ] << std::endl;
      }
    }

setInitialSimplex()

template<typename TKernel>

bool DGtal::QuickHull< TKernel >::setInitialSimplex ( const IndexRange & full_splx )

inline

Sets the initial full dimensional simplex

Precondition

status() must be Status::InputInitialized

Parameters

full_splx

a dimension+1-simplex specified as indices in the vector of input point.

Returns

'true' iff this initial simplex is full dimensional, the object has then the status SimplexCompleted, otherwise returns 'false' and the status is NotFullDimensional.

Definition at line 410 of file QuickHull.h.

    {
      if ( status() != Status::InputInitialized ) return false;
      if ( full_splx.size() != dimension + 1 )
        {
          trace.error() << "[QuickHull::setInitialSimplex]"
                        << " not a full dimensional simplex" << std::endl;
          myStatus = Status::NotFullDimensional;
          return false;
        }          
      CombinatorialPlaneSimplex splx;
      for ( Index j = 0; j < dimension; ++j )
        splx[ j ] = full_splx[ j ];
      const auto      H = kernel.compute( points, splx, full_splx.back() );
      const auto volume = kernel.volume( H, points[ full_splx.back() ] );
      if ( volume > 0 )
        return computeSimplexConfiguration( full_splx );
      myStatus = Status::NotFullDimensional;
      return false;
    }

setInput()

template<typename TKernel>

template<typename InputPoint>

bool DGtal::QuickHull< TKernel >::setInput ( const std::vector< InputPoint > & input_points, bool remove_duplicates = true )

inline

Sets the input data for the QuickHull convex hull algorithm, which is a range of points.

Template Parameters

InputPoint

a model of point that is convertible to Point datatype.

Parameters

[in]	input_points	the range of input points.
[in]	remove_duplicates	should be set to 'true' if the input data has duplicates.

Returns

'true' if the object is successfully initialized, status must be Status::InputInitialized, 'false' otherwise.

Examples

geometry/tools/exampleLatticeBallDelaunay2D.cpp.

Definition at line 382 of file QuickHull.h.

    {
      Clock tic;
      tic.startClock();
      clear();
      timings.clear();
      kernel.makeInput( points, input2comp,  comp2input,
                        input_points, remove_duplicates );
      timings.push_back( tic.stopClock() );
      if ( points.size() <= dimension ) {
        myStatus = Status::NotFullDimensional;
        return false;
      }
      myStatus = Status::InputInitialized;
      return true;
    }

Referenced by main().

status()

template<typename TKernel>

Status DGtal::QuickHull< TKernel >::status ( ) const

inline

Returns

the current status of this object, in Uninitialized, NotFullDimensional, InputInitialized, Completed, InvalidRidge, InvalidConvexHull

Definition at line 266 of file QuickHull.h.

    { return myStatus; }

Referenced by DGtal::QuickHull< Kernel3D >::check(), DGtal::QuickHull< Kernel3D >::computeConvexHull(), DGtal::QuickHull< Kernel3D >::computeFacets(), DGtal::QuickHull< Kernel3D >::computeSimplexConfiguration(), DGtal::QuickHull< Kernel3D >::computeVertices(), DGtal::QuickHull< Kernel3D >::getFacetHalfSpaces(), DGtal::QuickHull< Kernel3D >::getFacetVertices(), DGtal::QuickHull< Kernel3D >::getPoint2Vertex(), DGtal::QuickHull< Kernel3D >::getVertex2Point(), DGtal::QuickHull< Kernel3D >::getVertexPositions(), DGtal::QuickHull< Kernel3D >::isValid(), DGtal::QuickHull< Kernel3D >::nbFacets(), DGtal::QuickHull< Kernel3D >::nbFiniteFacets(), DGtal::QuickHull< Kernel3D >::nbInfiniteFacets(), DGtal::QuickHull< Kernel3D >::nbVertices(), DGtal::QuickHull< Kernel3D >::processFacet(), DGtal::QuickHull< Kernel3D >::selfDisplay(), and DGtal::QuickHull< Kernel3D >::setInitialSimplex().

unmakeNeighbors()

template<typename TKernel>

void DGtal::QuickHull< TKernel >::unmakeNeighbors ( const Index if1, const Index if2 )

inlineprotected

Makes two distinct facets if1 and if2 no more neighbors

Parameters

if1	a valid facet index
if2	a valid facet index

Definition at line 1213 of file QuickHull.h.

    {
      facets[ if1 ].subNeighbor( if2 );
      facets[ if2 ].subNeighbor( if1 );
    }

Field Documentation

assignment

template<typename TKernel>

std::vector< Index > DGtal::QuickHull< TKernel >::assignment

assignment of points to facets

Definition at line 812 of file QuickHull.h.

comp2input

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::comp2input

the injective mapping between the convex hull point range and the input range.

Definition at line 808 of file QuickHull.h.

debug_level

template<typename TKernel>

int DGtal::QuickHull< TKernel >::debug_level

debug_level from 0:no to 2

Definition at line 800 of file QuickHull.h.

deleted_facets

template<typename TKernel>

std::set< Index > DGtal::QuickHull< TKernel >::deleted_facets

set of deleted facets

Definition at line 816 of file QuickHull.h.

dimension

template<typename TKernel>

const Size DGtal::QuickHull< TKernel >::dimension = Point::dimension

static

Definition at line 152 of file QuickHull.h.

facet_counter

template<typename TKernel>

std::vector< Size > DGtal::QuickHull< TKernel >::facet_counter

Counts the number of facets with a given number of vertices.

Definition at line 828 of file QuickHull.h.

facets

template<typename TKernel>

std::vector< Facet > DGtal::QuickHull< TKernel >::facets

the current set of facets.

Definition at line 814 of file QuickHull.h.

input2comp

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::input2comp

the surjective mapping between the input range and the output range used for convex hull computation.

Definition at line 805 of file QuickHull.h.

kernel

template<typename TKernel>

Kernel DGtal::QuickHull< TKernel >::kernel

mutable

Kernel that is duplicated for computing facet geometries.

Definition at line 798 of file QuickHull.h.

myStatus

template<typename TKernel>

Status DGtal::QuickHull< TKernel >::myStatus

protected

The status of the object: Uninitialized, InputInitialized, SimplexCompleted, FacetsCompleted, VerticesCompleted, InvalidRidge, InvalidConvexHull, NotFullDimensional

Definition at line 839 of file QuickHull.h.

nb_finite_facets

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nb_finite_facets

Number of finite facets.

Definition at line 822 of file QuickHull.h.

nb_infinite_facets

template<typename TKernel>

Size DGtal::QuickHull< TKernel >::nb_infinite_facets

Number of infinite facets (!= 0 only for specific kernels)

Definition at line 824 of file QuickHull.h.

p2v

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::p2v

point index -> vertex index (or UNASSIGNED)

Definition at line 818 of file QuickHull.h.

points

template<typename TKernel>

std::vector< Point > DGtal::QuickHull< TKernel >::points

the set of points, indexed as in the array.

Definition at line 802 of file QuickHull.h.

processed_points

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::processed_points

Points already processed (and within the convex hull).

Definition at line 810 of file QuickHull.h.

timings

template<typename TKernel>

std::vector< double > DGtal::QuickHull< TKernel >::timings

Timings of the different phases: 0: init, 1: facets, 2: vertices.

Definition at line 826 of file QuickHull.h.

Referenced by main().

v2p

template<typename TKernel>

IndexRange DGtal::QuickHull< TKernel >::v2p

vertex index -> point index

Definition at line 820 of file QuickHull.h.

---

The documentation for this struct was generated from the following file:

• QuickHull.h