ConvexityHelper.ih

/**
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or  (at your option) any later version.
 *
 *  This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 **/
 
/**
 * @file ConvexityHelper.ih
 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
 *
 * @date 2020/01/02
 *
 * Implementation of inline methods defined in ConvexityHelper.h
 *
 * This file is part of the DGtal library.
 */
 
 
//////////////////////////////////////////////////////////////////////////////
#include <cstdlib>
#include <string>
#include <sstream>
#include "DGtal/kernel/IntegerConverter.h"
//////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
// IMPLEMENTATION of inline methods.
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
// ----------------------- Standard services ------------------------------
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSimplex
( const PointRange& input_points,
  bool remove_duplicates )
{
  PointRange X;
  if ( remove_duplicates )
    {
      std::set<Point> S;
      for ( auto&& p : input_points ) S.insert( p );
      X = PointRange( S.cbegin(), S.cend() );
    }
  else X = input_points;
  LatticePolytope P( X.cbegin(), X.cend() );
  if ( P.nbHalfSpaces() != 0 )
    return P;
  else
    return computeDegeneratedLatticePolytope( X );
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
computeDegeneratedLatticePolytope
( PointRange& input_points )
{
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  typedef typename ConvexHull::HalfSpace       ConvexHullHalfSpace;
  // input_points is a range of points with no duplicates, but which
  // seems to be not full dimensional.
  if ( input_points.size() <= 1 )
    return LatticePolytope( input_points.cbegin(), input_points.cend() );
  // At least 1-dimensional
  std::vector< Vector >  basis;
  std::vector< Integer > alpha;
  basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
  const auto n0 = basis[ 0 ].norm();
  alpha.push_back( Integer( 0 ) );
  alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
  Index i = 2;
  while ( i < input_points.size() ) {
    Vector v = input_points[ i ] - input_points[ 0 ];
    alpha.push_back( basis[ 0 ].dot( v ) );
    const auto ni = v.norm();
    const double alignment =
      fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
            - ( n0 * ni ) );
    if ( alignment > 1e-8 ) break;
    i++;
  }
  if ( i == input_points.size() ) 
    { // 1-dimensional
      Index a = 0, b = 0;
      for ( i = 1; i < input_points.size(); i++ )
        {
          if ( alpha[ i ] < alpha[ a ] ) a = i;
          if ( alpha[ i ] > alpha[ b ] ) b = i;
        }
      PointRange X( 2 );
      X[ 0 ] = input_points[ a ];
      X[ 1 ] = input_points[ b ];
      return LatticePolytope( X.cbegin(), X.cend() );
    }
  // at least 2-dimensional
  ASSERT( dimension > 1 );
  if ( dimension == 2 )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << " Weird error: found initial full dimensional simplex" << std::endl;
      return LatticePolytope();
    }
  if ( dimension >= 4 )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
                << std::endl;
      return LatticePolytope();
    }
  basis.push_back( input_points[ i ] - input_points[ 0 ] );
  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
    ::crossProduct( basis[ 0 ], basis[ 1 ] );
  if ( n == Vector::zero )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Weird numerical error, u . v != |u| |v| but u x v != 0"
                << std::endl;
      return LatticePolytope();
    }
  // Now the set of input points should be full dimensional.
  input_points.push_back( input_points[ 0 ] + n );
  // Compute domain
  Point l = input_points[ 0 ];
  Point u = input_points[ 0 ];
  for ( const auto& p : input_points ) {
    l = l.inf( p );
    u = u.sup( p );
  }
  Domain domain( l, u );
  // Compute convex hull
  ConvexHull hull;
  hull.setInput( input_points, false );
  const auto target = ConvexHull::Status::FacetsCompleted;
  IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
  bool ok_init = hull.setInitialSimplex( full_splx );
  if ( ! ok_init )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Weird error in hull.setInitialSimplex" << std::endl;
      return LatticePolytope();
    }
  bool ok_hull = hull.computeConvexHull( target );
  if ( ! ok_hull )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Weird error in hull.computeConvexHull" << std::endl;
      return LatticePolytope();
    }
  // Initialize polytope
  std::vector< ConvexHullHalfSpace > HS;
  std::vector< PolytopeHalfSpace >   PHS;
  hull.getFacetHalfSpaces( HS );
  PHS.reserve( HS.size()+2 );
  for ( auto& H : HS ) {
    Vector  N;
    Integer nu;
    for ( Dimension i = 0; i < dim; ++i )
      N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
    nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
    PHS.emplace_back( N, nu );
  }
  // Add top constraint.
  Integer nu0 = input_points[ 0 ].dot( n );
  PHS.emplace_back( n, nu0 );
  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
                          false, false );
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
computeDegeneratedConvexHullVertices
( PointRange& input_points )
{
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  typedef typename ConvexHull::HalfSpace       ConvexHullHalfSpace;
  // input_points is a range of points with no duplicates, but which
  // seems to be not full dimensional.
  if ( input_points.size() <= 1 )
    return input_points;
  // At least 1-dimensional
  std::vector< Vector >  basis;
  std::vector< Integer > alpha;
  basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
  const auto n0 = basis[ 0 ].norm();
  alpha.push_back( Integer( 0 ) );
  alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
  Index i = 2;
  while ( i < input_points.size() ) {
    Vector v = input_points[ i ] - input_points[ 0 ];
    alpha.push_back( basis[ 0 ].dot( v ) );
    const auto ni = v.norm();
    const double alignment =
      fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
            - ( n0 * ni ) );
    if ( alignment > 1e-8 ) break;
    i++;
  }
  if ( i == input_points.size() ) 
    { // 1-dimensional
      Index a = 0, b = 0;
      for ( i = 1; i < input_points.size(); i++ )
        {
          if ( alpha[ i ] < alpha[ a ] ) a = i;
          if ( alpha[ i ] > alpha[ b ] ) b = i;
        }
      PointRange X( 2 );
      X[ 0 ] = input_points[ a ];
      X[ 1 ] = input_points[ b ];
      return X;
    }
  // at least 2-dimensional
  ASSERT( dimension > 1 );
  if ( dimension == 2 )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << " Weird error: found initial full dimensional simplex" << std::endl;
      return PointRange();
    }
  if ( dimension >= 4 )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
                << std::endl;
      return PointRange();
    }
  basis.push_back( input_points[ i ] - input_points[ 0 ] );
  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
    ::crossProduct( basis[ 0 ], basis[ 1 ] );
  if ( n == Vector::zero )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << "Weird numerical error, u . v != |u| |v| but u x v != 0"
                << std::endl;
      return PointRange();
    }
  // Now the set of input points should be full dimensional.
  input_points.push_back( input_points[ 0 ] + n );
  // Compute convex hull
  ConvexHull hull;
  hull.setInput( input_points, false );
  const auto target = ConvexHull::Status::VerticesCompleted;
  IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
  bool ok_init = hull.setInitialSimplex( full_splx );
  if ( ! ok_init )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << "Weird error in hull.setInitialSimplex" << std::endl;
      return PointRange();
    }
  bool ok_hull = hull.computeConvexHull( target );
  if ( ! ok_hull )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << "Weird error in hull.computeConvexHull" << std::endl;
      return PointRange();
    }
  // Get convex hull vertices and remove top point
  PointRange X;
  hull.getVertexPositions( X );
  const std::size_t nX = X.size();
  for ( std::size_t j = 0; j < nX; j++ )
    if ( X[ j ] == input_points.back() )
      {
        X[ j ] = X.back();
        break;
      }
  X.resize( nX - 1 ); 
  return X;
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
computeLatticePolytope
( const PointRange& input_points,
  bool remove_duplicates,
  bool make_minkowski_summable )
{
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::HalfSpace       ConvexHullHalfSpace;
  typedef typename ConvexHull::Ridge           Ridge;
  if ( input_points.empty() ) return LatticePolytope();
  if ( input_points.size() <= ( dimension + 1) )
    return computeSimplex( input_points, remove_duplicates );
  // Compute domain
  Point l = input_points[ 0 ];
  Point u = input_points[ 0 ];
  for ( const auto& p : input_points ) {
    l = l.inf( p );
    u = u.sup( p );
  }
  Domain domain( l, u );
  // Compute convex hull
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  const auto target = ( make_minkowski_summable && dimension == 3 )
    ? ConvexHull::Status::VerticesCompleted
    : ConvexHull::Status::FacetsCompleted;
  bool ok = hull.computeConvexHull( target );
  if ( ! ok ) // set of points is not full dimensional
    return computeDegeneratedLatticePolytope( hull.points );
  // Initialize polytope
  std::vector< ConvexHullHalfSpace > HS;
  std::vector< PolytopeHalfSpace >   PHS;
  hull.getFacetHalfSpaces( HS );
  PHS.reserve( HS.size() );
  for ( auto& H : HS ) {
    Vector  N;
    Integer nu;
    for ( Dimension i = 0; i < dim; ++i )
      N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
    nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
    PHS.emplace_back( N, nu );
  }
  if ( make_minkowski_summable && dimension >= 4 )
    trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
                    << " Not implemented starting from dimension 4."
                    << std::endl;
  if ( make_minkowski_summable && dimension == 3 )
    {
      // Compute ridge vertices to add edge constraints.
      PointRange positions;
      std::vector< IndexRange > facet_vertices; 
      std::vector< IndexRange > ridge_vertices; 
      std::map< Ridge, Index > ridge2index;
      hull.getVertexPositions( positions );
      computeFacetAndRidgeVertices( hull, facet_vertices,
                                    ridge2index, ridge_vertices );
      for ( auto p : ridge2index ) {
        const auto   r = p.first;
        // Copy by value since PHS may be reallocated  during the iteration.
        const auto   U = PHS[ r.first  ].N; // normal of facet 1
        const auto   V = PHS[ r.second ].N; // normal of facet 2
        const auto&  S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
        ASSERT( S.size() == 2 && "Invalid ridge" );
        const auto& P0 = positions[ S[ 0 ] ];
        const auto& P1 = positions[ S[ 1 ] ];
        auto         E = P1 - P0; // edge 1, 2
        const auto UxV =
          detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
          ::crossProduct( U, V ); // parallel to E
        ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
        if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
        const auto  E1 = 
          detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
          ::crossProduct( U, E ); // edge on facet 1
        const auto  E2 = 
          detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
          ::crossProduct( E, V ); // edge on facet 2
        ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
        ASSERT( E1.dot( V ) <  0 && "Invalid E1 / V" );
        ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
        ASSERT( E2.dot( U ) <  0 && "Invalid E1 / U" );
        for ( Dimension k = 0; k < dimension; ++k ) {
          const auto W = U[ k ] * V - V[ k ] * U;
          const auto nn1 = W.dot( E1 ); 
          const auto nn2 = W.dot( E2 ); 
          if ( nn1 > 0 && nn2 > 0 ) {
            PHS.emplace_back( -W, -W.dot( P0 ) );
            ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
            ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
          }
          else if ( nn1 < 0 && nn2 < 0 ) {
            PHS.emplace_back( W, W.dot( P0 ) );
            ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
            ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
          }
        }
      }
    }
  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
                          make_minkowski_summable && ( dimension <= 3 ), true );
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
computeConvexHullVertices
( const PointRange& input_points,
  bool remove_duplicates )
{
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  PointRange positions;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  if ( !ok )
    {
      PointRange Z( input_points );
      return computeDegeneratedConvexHullVertices( Z );
    }
  hull.getVertexPositions( positions );
  return positions;
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
template < typename TSurfaceMesh >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
( TSurfaceMesh&               mesh,
  const PointRange& input_points,
  bool remove_duplicates )
{
  typedef TSurfaceMesh                         SurfaceMesh;
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  if ( !ok ) return false;
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
                      faces.cbegin(), faces.cend() );
  return true;
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
( PolygonalSurface< Point >&  polysurf,
  const PointRange& input_points,
  bool remove_duplicates )
{
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  if ( !ok ) return false;
  PointRange positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  // build polygonal surface
  polysurf.clear();
  for ( auto p : positions ) polysurf.addVertex( p );
  for ( auto f : faces )     polysurf.addPolygonalFace( f );
  return polysurf.build();
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
( ConvexCellComplex< Point >& cell_complex,
  const PointRange& input_points,
  bool remove_duplicates )
{
  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  cell_complex.clear();
  if ( ! ok ) return false;
  // Build complex, only 1 finite cell and as many faces as convex hull facets.
  // Taking care of faces for each cell (here one cell borders all faces).
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  FaceRange all_faces;
  for ( Index i = 0; i < faces.size(); i++ )
    all_faces.push_back( { i, true } );
  cell_complex.cell_faces.push_back( all_faces );
  // Vertices of this unique cell will be computed lazily on request.
  // Taking care of each face.
  for ( const auto& vtcs : faces )
    {
      // every inner face borders cell 0
      cell_complex.true_face_cell.push_back( 0 );
      // every outer face borders the infinite cell
      cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
    }
  // Taking care of vertices (in consistent order) of every face
  cell_complex.true_face_vertices.swap( faces );
  // Taking care of vertex positions
  hull.getVertexPositions( cell_complex.vertex_position );
  return true;
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
( ConvexCellComplex< Point >& cell_complex,
  const PointRange& input_points,
  bool remove_duplicates )
{
  typedef QuickHull< LatticeDelaunayKernel >   Delaunay;
  typedef typename Delaunay::Ridge             Ridge;
  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
  
  Delaunay del;
  del.setInput( input_points, remove_duplicates );
  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
  cell_complex.clear();
  if ( ! ok ) return false;
  
  // Build complex, as many maximal cells as convex hull facets.
  // convex hull facet -> cell of complex
  // convex hull ridge -> face of complex
  // (1) Get cell vertices, count ridges/faces and compute their vertices
  std::map< Ridge, Index > r2f;
  computeFacetAndRidgeVertices( del,
                                cell_complex.cell_vertices,
                                r2f,
                                cell_complex.true_face_vertices );
  // (2) assign ridges/faces to cell and conversely
  const Index  nb_r = r2f.size();
  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
  cell_complex.true_face_vertices.resize( nb_r );
  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) { 
    const auto& facet = del.facets[ cur_f ];
    FaceRange current_faces;
    for ( auto neigh_f : facet.neighbors ) {
      const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
      const bool pos = cur_f < neigh_f;
      const Index cur_r = r2f[ r ];
      cell_complex.true_face_cell [ cur_r ] = r.first;
      if ( r.second >= del.nbFiniteFacets() )
        cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
      else
        cell_complex.false_face_cell[ cur_r ] = r.second;
      current_faces.emplace_back( cur_r, pos );
    }
    cell_complex.cell_faces.push_back( current_faces );
  }
  // (3) Takes care of vertex positions
  del.getVertexPositions( cell_complex.vertex_position );
  return true;
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
( const std::vector< RealPoint >& input_points,
  Integer denominator, 
  bool remove_duplicates,
  bool make_minkowski_summable )
{
  if ( denominator < 1 )
    trace.error() << "Invalid denominator " << denominator
                  << ". Should be greater or equal to 1." << std::endl;
  typedef typename RationalPolytope::Domain     Domain;
  typedef typename RationalPolytope::HalfSpace  PolytopeHalfSpace;
  typedef QuickHull< RealConvexHullKernel >     ConvexHull;
  typedef typename ConvexHull::HalfSpace        ConvexHullHalfSpace;
  typedef typename ConvexHull::Ridge            Ridge;
  if ( input_points.empty() ) return RationalPolytope();
  // Compute convex hull
  ConvexHull hull( denominator );
  hull.setInput( input_points, remove_duplicates );
  const auto target = ( make_minkowski_summable && dimension == 3 )
    ? ConvexHull::Status::VerticesCompleted
    : ConvexHull::Status::FacetsCompleted;
  bool ok = hull.computeConvexHull( target );
  if ( ! ok ) return RationalPolytope();
  // Compute domain (as a lattice domain)
  auto l = hull.points[ 0 ];
  auto u = hull.points[ 0 ];
  for ( const auto& p : hull.points ) {
    l = l.inf( p );
    u = u.sup( p );
  }
  Domain domain( l, u );
  trace.info() << "Domain l=" << l << " u=" << u << std::endl;
  // Initialize polytope
  std::vector< ConvexHullHalfSpace > HS;
  std::vector< PolytopeHalfSpace >   PHS;
  hull.getFacetHalfSpaces( HS );
  PHS.reserve( HS.size() );
  for ( auto& H : HS ) {
    Vector  N;
    Integer nu;
    for ( Dimension i = 0; i < dim; ++i )
      N[ i ] = (Integer) H.internalNormal()[ i ];
    nu = (Integer) H.internalIntercept();
    PHS.emplace_back( N, nu );
  }
  if ( make_minkowski_summable && dimension >= 4 )
    trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
                    << " Not implemented starting from dimension 4."
                    << std::endl;
  if ( make_minkowski_summable && dimension == 3 )
    {
      // Compute ridge vertices to add edge constraints.
      PointRange positions;
      std::vector< IndexRange > facet_vertices; 
      std::vector< IndexRange > ridge_vertices; 
      std::map< Ridge, Index > ridge2index;
      hull.getVertexPositions( positions );
      computeFacetAndRidgeVertices( hull, facet_vertices,
                                    ridge2index, ridge_vertices );
      for ( auto p : ridge2index ) {
        const auto   r = p.first;
        // Copy by value since PHS may be reallocated  during the iteration.
        const auto   U = PHS[ r.first  ].N; // normal of facet 1
        const auto   V = PHS[ r.second ].N; // normal of facet 2
        const auto&  S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
        ASSERT( S.size() == 2 && "Invalid ridge" );
        const auto& P0 = positions[ S[ 0 ] ];
        const auto& P1 = positions[ S[ 1 ] ];
        auto         E = P1 - P0; // edge 1, 2
        const auto UxV =
          detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
          ::crossProduct( U, V ); // parallel to E
        ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
        if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
        const auto  E1 = 
          detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
          ::crossProduct( U, E ); // edge on facet 1
        const auto  E2 = 
          detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
          ::crossProduct( E, V ); // edge on facet 2
        ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
        ASSERT( E1.dot( V ) <  0 && "Invalid E1 / V" );
        ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
        ASSERT( E2.dot( U ) <  0 && "Invalid E1 / U" );
        for ( Dimension k = 0; k < dimension; ++k ) {
          const auto W = U[ k ] * V - V[ k ] * U;
          const auto nn1 = W.dot( E1 ); 
          const auto nn2 = W.dot( E2 ); 
          if ( nn1 > 0 && nn2 > 0 ) {
            PHS.emplace_back( -W, -W.dot( P0 ) );
            ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
            ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
          }
          else if ( nn1 < 0 && nn2 < 0 ) {
            PHS.emplace_back( W, W.dot( P0 ) );
            ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
            ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
          }
        }
      }
    }
  return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
                           make_minkowski_summable && ( dimension <= 3 ), true );
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
template < typename TSurfaceMesh >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
( TSurfaceMesh&               mesh,
  const std::vector< RealPoint >& input_points,
  double precision,
  bool remove_duplicates )
{
  typedef TSurfaceMesh                         SurfaceMesh;
  typedef QuickHull< RealConvexHullKernel >    ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  ConvexHull hull( precision );
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  if ( !ok ) return false;
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
                      faces.cbegin(), faces.cend() );
  return true;
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
( PolygonalSurface< RealPoint >&  polysurf,
  const std::vector< RealPoint >& input_points,
  double precision,
  bool remove_duplicates )
{
  typedef QuickHull< RealConvexHullKernel >    ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  ConvexHull hull( precision );
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  if ( !ok ) return false;
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  // build polygonal surface
  polysurf.clear();
  for ( auto p : positions ) polysurf.addVertex( p );
  for ( auto f : faces )     polysurf.addPolygonalFace( f );
  return polysurf.build();
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
( ConvexCellComplex< RealPoint >& cell_complex,
  const std::vector< RealPoint >& input_points,
  double precision,
  bool remove_duplicates )
{
  typedef QuickHull< RealConvexHullKernel >    ConvexHull;
  typedef typename ConvexHull::IndexRange      IndexRange;
  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
  ConvexHull hull( precision );
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  cell_complex.clear();
  if ( ! ok ) return false;
  // Build complex, only 1 finite cell and as many faces as convex hull facets.
  // Taking care of faces for each cell (here one cell borders all faces).
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  FaceRange all_faces;
  for ( Index i = 0; i < faces.size(); i++ )
    all_faces.push_back( { i, true } );
  cell_complex.cell_faces.push_back( all_faces );
  // Vertices of this unique cell will be computed lazily on request.
  // Taking care of each face.
  for ( const auto& vtcs : faces )
    {
      // every inner face borders cell 0
      cell_complex.true_face_cell.push_back( 0 );
      // every outer face borders the infinite cell
      cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
    }
  // Taking care of vertices (in consistent order) of every face
  cell_complex.true_face_vertices.swap( faces );
  // Taking care of vertex positions
  hull.getVertexPositions( cell_complex.vertex_position );
  return true;
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
bool
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
( ConvexCellComplex< RealPoint >& cell_complex,
  const std::vector< RealPoint >& input_points,
  double precision,
  bool remove_duplicates )
{
  typedef QuickHull< RealDelaunayKernel >      Delaunay;
  typedef typename Delaunay::Ridge             Ridge;
  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
  
  Delaunay del( precision );
  del.setInput( input_points, remove_duplicates );
  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
  cell_complex.clear();
  if ( ! ok ) return false;
  
  // Build complex, as many maximal cells as convex hull facets.
  // convex hull facet -> cell of complex
  // convex hull ridge -> face of complex
  // (1) Get cell vertices, count ridges/faces and compute their vertices
  std::map< Ridge, Index > r2f;
  computeFacetAndRidgeVertices( del,
                                cell_complex.cell_vertices,
                                r2f,
                                cell_complex.true_face_vertices );
  // (2) assign ridges/faces to cell and conversely
  const Index  nb_r = r2f.size();
  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
  cell_complex.true_face_vertices.resize( nb_r );
  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) { 
    const auto& facet = del.facets[ cur_f ];
    FaceRange current_faces;
    for ( auto neigh_f : facet.neighbors ) {
      const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
      const bool pos = cur_f < neigh_f;
      const Index cur_r = r2f[ r ];
      cell_complex.true_face_cell [ cur_r ] = r.first;
      if ( r.second >= del.nbFiniteFacets() )
        cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
      else
        cell_complex.false_face_cell[ cur_r ] = r.second;
      current_faces.emplace_back( cur_r, pos );
    }
    cell_complex.cell_faces.push_back( current_faces );
  }
  // (3) Takes care of vertex positions
  del.getVertexPositions( cell_complex.vertex_position );
  return true;
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
template < typename QHull >
void
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
( const QHull& hull,
  std::vector< IndexRange >& cell_vertices,
  std::map< typename QHull::Ridge, Index >& r2f,
  std::vector< IndexRange >& face_vertices )
{
  typedef typename QHull::Ridge             Ridge;
 
  ASSERT( hull.status() >= QHull::Status::VerticesCompleted
          && hull.status() <= QHull::Status::AllCompleted );
  
  // Get cell vertices and sort them
  bool ok_fv = hull.getFacetVertices( cell_vertices );
  if ( ! ok_fv )
    trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
                  << " method hull.getFacetVertices failed."
                  << " Maybe QuickHull was not computed till VerticesCompleted."
                  << std::endl;
  std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
  for ( auto& vtcs : sorted_cell_vertices )
    std::sort( vtcs.begin(), vtcs.end() );
  cell_vertices.resize( hull.nbFiniteFacets() );
  
  // Count ridges/faces and compute their vertices
  Index  nb_r = 0;
  face_vertices.clear();
  for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f )  { 
    const auto& facet = hull.facets[ cur_f ];
    for ( auto neigh_f : facet.neighbors ) {
      const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
      auto itr = r2f.find( r );
      if ( itr == r2f.end() ) {
        IndexRange result;
        std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
                               sorted_cell_vertices[ cur_f ].cend(),
                               sorted_cell_vertices[ neigh_f ].cbegin(),
                               sorted_cell_vertices[ neigh_f ].cend(), 
                               std::back_inserter( result ) );
        face_vertices.push_back( result );
        r2f[ r ] = nb_r++;
      }
    }
  }
}
 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////