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 "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::HalfSpace       ConvexHullHalfSpace;
 
  // Rewrite computeDegeneratedLatticePolytope using AffineSubset
  typedef PointVector<dim, InternalInteger>   InternalPoint;
  typedef AffineGeometry<InternalPoint>        Affine;
  auto  indices   = Affine::affineSubset( input_points );
  auto  ref_basis = Affine::affineBasis ( input_points, indices );
  auto  ref       = ref_basis.first;
  auto& basis     = ref_basis.second;
  if ( basis.size() >= dimension )
    {
      trace.error() << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                    << " Weird error: found initial full dimensional simplex" << std::endl;
      return LatticePolytope();
    }  
  if ( basis.size() == 0 )
    { // Polytope is degenerated to one point (nD)
      PointRange X = { ref };
      return LatticePolytope( X.cbegin(), X.cend() );
    }
  if ( basis.size() == 1 )
    { // Polytope is a 1-dimensional straight segment (nD).
      // Find the two extremal points.
      std::vector< InternalInteger > alpha( input_points.size() );
      for ( std::size_t i = 0; i < input_points.size(); i++ )
        alpha[ i ] = basis[ 0 ].dot( input_points[ i ] - ref );      
      Index a = 0, b = 0;
      for ( std::size_t i = 1; i < input_points.size(); i++ )
        {
          if ( alpha[ i ] < alpha[ a ] ) a = i;
          if ( alpha[ i ] > alpha[ b ] ) b = i;
        }
      PointRange X = { input_points[ a ], input_points[ b ] };
      return LatticePolytope( X.cbegin(), X.cend() );
    }
  if ( basis.size()+1 != dimension )
    {
      trace.error() << "[ConvexityHelper::computeDegeneratedLatticePolytope] "
                    << "computation of a " << basis.size() << "-dimensional polytope"
                    << " in " << dimension << "D is not implemented yet." << std::endl;
      return LatticePolytope();
    }
  // Polytope is d-1 dimensional (nD)
  // Compute an independent vector to the basis (one of the canonic vector)
  Point e = Affine::template independentVector<Point>( basis );
  // Compute an orthogonal vector to the basis
  Affine::completeBasis( basis, true );
  InternalPoint n = basis.back(); // used to compute the two half-spaces that sandwich the polytope
  
  // 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 );
 
  // We complete the affine subset so that it is full dimensional.
  // We could use `ref + n` but it induces bigger integers.
  // Now the set of input points should be full dimensional.
  input_points.push_back( ref + e );
  
  // Compute convex hull
  ConvexHull hull;
  hull.setInput( input_points, false );
  const auto target = ConvexHull::Status::FacetsCompleted;
  // the array `indices` provides the initial full dimensional simplex,
  // if we had the last point.
  indices.push_back( input_points.size() - 1 );
  bool ok_init = hull.setInitialSimplex( indices );
  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 ii = 0; ii < dim; ++ii )
      N[ ii ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ ii ] );
    nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
    PHS.emplace_back( N, nu );
  }
  // Add top constraint.
  Vector tn = IntegerConverter< dimension, Integer >::cast( n );
  Integer nu0 = input_points[ 0 ].dot( tn );
  PHS.emplace_back( tn, nu0 );
  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
                          false, false );
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DTriangle
( const Point& a, const Point& b, const Point& c,
  bool make_minkowski_summable )
{
  if constexpr( dim != 3 ) return LatticePolytope();
  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  // Compute domain
  const Vector ab   = b - a;
  const Vector bc   = c - b;
  const Vector ca   = a - c;
  const Vector n    = Op::crossProduct( ab, bc );
  if ( n == Vector::zero )
    return computeDegeneratedTriangle( a, b, c );
  const Point  low  = a.inf( b ).inf( c );
  const Point  high = a.sup( b ).sup( c );
  // Initialize polytope
  std::vector< PolytopeHalfSpace >   PHS;
  PHS.reserve( make_minkowski_summable ? 11 : 5 );
  const Integer n_a = n.dot( a );
  const Vector  u   = Op::crossProduct( ab, n );
  const Vector  v   = Op::crossProduct( bc, n );
  const Vector  w   = Op::crossProduct( ca, n );
  PHS.emplace_back(  n,  n_a  );
  PHS.emplace_back( -n, -n_a );
  if ( ! make_minkowski_summable )
    { // It is enough to have one constraint per edge.
      PHS.emplace_back(  u, u.dot( a ) );
      PHS.emplace_back(  v, v.dot( b ) );
      PHS.emplace_back(  w, w.dot( c ) );
    }
  else
    { // Compute additionnal constraints on edges so that the
      // Minkowski sum with axis-aligned edges is valid.
      for ( Integer d = -1; d <= 1; d += 2 )
        for ( Dimension k = 0; k < dim; k++ )
          {
            const Vector  i     = Vector::base( k, d );
            const Vector  eab   = Op::crossProduct( ab, i );
            const Integer eab_a = eab.dot( a );
            if ( eab.dot( c ) < eab_a ) // c must be below plane (a,eab)
              PHS.emplace_back( eab, eab_a );
            const Vector  ebc   = Op::crossProduct( bc, i );
            const Integer ebc_b = ebc.dot( b );
            if ( ebc.dot( a ) < ebc_b ) // a must be below plane (b,ebc)
              PHS.emplace_back( ebc, ebc_b );
            const Vector  eca   = Op::crossProduct( ca, i );
            const Integer eca_c = eca.dot( c );
            if ( eca.dot( b ) < eca_c ) // b must be below plane (c,eca)
              PHS.emplace_back( eca, eca_c );
          }
    }
  return LatticePolytope( Domain( low, high ),
                          PHS.cbegin(), PHS.cend(),
                          make_minkowski_summable, false );
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DOpenTriangle
( const Point& a, const Point& b, const Point& c,
  bool make_minkowski_summable )
{
  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
  const Vector ab   = b - a;
  const Vector bc   = c - b;
  const Vector n    = Op::crossProduct( ab, bc );
  if ( n == Vector::zero )
    return LatticePolytope(); // empty polytope
  auto P = compute3DTriangle( a, b, c, make_minkowski_summable );
  // Compute domain
  for ( auto k = 0; k < P.nbHalfSpaces(); k++ )
    {
      if ( Op::crossProduct( P.getA( k ), n ) != Vector::zero )
        P.setStrict( k );
    }
  return P;
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
computeDegeneratedTriangle
( const Point& a, const Point& b, const Point& c )
{
  if ( a == b ) return computeSegment( a, c );
  if ( ( a == c ) || ( b == c ) ) return computeSegment( a, b );
  // The three points are distinct, hence aligned. One is in-between the two others.
  const Point  low  = a.inf( b ).inf( c );
  const Point  high = a.sup( b ).sup( c );
  for ( Dimension k = 0; k < dim; k++ )
    {
      const auto lk = low [ k ];
      const auto hk = high[ k ];      
      if ( ( a[ k ] != lk ) && ( a[ k ] != hk ) ) return computeSegment( b, c );
      if ( ( b[ k ] != lk ) && ( b[ k ] != hk ) ) return computeSegment( a, c );
      if ( ( c[ k ] != lk ) && ( c[ k ] != hk ) ) return computeSegment( a, b );
    }
  trace.error() << "[ConvexityHelper::computeSegmentFromDegeneratedTriangle] "
                << "Should never arrive here." << std::endl;
  return computeSegment( a, a );
}
 
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeOpenSegment
( const Point& a, const Point& b )
{
 
  if constexpr( dim != 3 ) return LatticePolytope();
  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  // Compute domain
  const Point  low  = a.inf( b );
  const Point  high = a.sup( b );
  const Vector ab   = b - a;
  bool  degenerate  = ( ab == Vector::zero );
  if ( degenerate ) return LatticePolytope();
  // Initialize polytope
  std::vector< PolytopeHalfSpace >   PHS;
  PHS.reserve( 4*dim );
  // Compute additionnal constraints on domain boundary to make it open.
  for ( Dimension k = 0; k < dim; k++ )
    {
      const Vector bp = Vector::base( k,  1 );
      PHS.emplace_back( bp, high[ k ] );
      const Vector bn = Vector::base( k, -1 );
      PHS.emplace_back( bn, -low[ k ] );
    }      
  // Compute additionnal constraints on edges so that the
  // Minkowski sum with axis-aligned edges is valid.
  for ( Integer d = -1; d <= 1; d += 2 )
    for ( Dimension k = 0; k < dim; k++ )
      {
        const Vector  i   = Vector::base( k, d );
        const Vector  e   = Op::crossProduct( ab, i );
        if ( e != Vector::zero )
          {
            const Integer e_a = e.dot( a );
            PHS.emplace_back( e, e_a );
          }
      }
  auto P = LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
  // Fix < inequalities
  for ( auto k = 2 * dim; k < 4 * dim; k++ )
    {
      auto V = P.getA( k );
      if ( V.dot( a ) != V.dot( b ) ) P.setStrict( k );
    }
  return P;
 
}
 
//-----------------------------------------------------------------------------
template < int dim, typename TInteger, typename TInternalInteger >
typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSegment
( const Point& a, const Point& b )
{
  if constexpr( dim != 3 ) return LatticePolytope();
  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
  typedef typename LatticePolytope::Domain     Domain;
  typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
  // Compute domain
  const Point  low  = a.inf( b );
  const Point  high = a.sup( b );
  const Vector ab   = b - a;
  bool  degenerate  = ( ab == Vector::zero );
  // Initialize polytope
  std::vector< PolytopeHalfSpace >   PHS;
  if ( degenerate ) 
    return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
  PHS.reserve( 2*dim );
  // Compute additionnal constraints on edges so that the
  // Minkowski sum with axis-aligned edges is valid.
  for ( Integer d = -1; d <= 1; d += 2 )
    for ( Dimension k = 0; k < dim; k++ )
      {
        const Vector  i   = Vector::base( k, d );
        const Vector  e   = Op::crossProduct( ab, i );
        if ( e != Vector::zero )
          {
            const Integer e_a = e.dot( a );
            PHS.emplace_back( e, e_a );
          }
      }
  return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, 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 QuickHull< LatticeConvexHullKernel > ConvexHull;
 
  // Rewrite computeDegeneratedLatticePolytope using AffineSubset
  typedef PointVector<dim, TInternalInteger>   InternalPoint;
  typedef AffineGeometry<InternalPoint>        Affine;
  auto  indices   = Affine::affineSubset( input_points );
  auto  ref_basis = Affine::affineBasis ( input_points, indices );
  auto  ref       = ref_basis.first;
  auto& basis     = ref_basis.second;
  if ( basis.size() >= dimension )
    {
      trace.error() << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                    << " Weird error: found initial full dimensional simplex" << std::endl;
      return PointRange();
    }  
  if ( basis.size() == 0 )
    { // Polytope is degenerated to one point (nD)
      PointRange X = { ref };
      return X;
    }
  if ( basis.size() == 1 )
    { // Polytope is a 1-dimensional straight segment (nD).
      // Find the two extremal points.
      std::vector< Integer > alpha( input_points.size() );
      for ( std::size_t i = 0; i < input_points.size(); i++ )
        alpha[ i ] = basis[ 0 ].dot( input_points[ i ] - ref );      
      Index a = 0, b = 0;
      for ( std::size_t i = 1; i < input_points.size(); i++ )
        {
          if ( alpha[ i ] < alpha[ a ] ) a = i;
          if ( alpha[ i ] > alpha[ b ] ) b = i;
        }
      PointRange X = { input_points[ a ], input_points[ b ] };
      return X;
    }
  if ( basis.size()+1 != dimension )
    {
      trace.error() << "[ConvexityHelper::computeDegeneratedLatticePolytope] "
                    << "computation of a " << basis.size() << "-dimensional polytope"
                    << " in " << dimension << "D is not implemented yet." << std::endl;
      return PointRange();
    }
  // Polytope is d-1 dimensional (nD)
 
  // old code
  // typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
  // // 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();
  //   }
 
  // Compute an independent vector to the basis (one of the canonic vector)
  Point e = Affine::template independentVector<Point>( basis );
  if ( e == Point::zero )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
                << "Weird numerical error, the basis is not full dimensional and there is no independent canonic vector !"
                << std::endl;
      return PointRange();
    }
  // We complete the affine subset so that it is full dimensional.
  // Now the set of input points should be full dimensional.
  input_points.push_back( ref + e );
  // Compute convex hull
  ConvexHull hull;
  hull.setInput( input_points, false );
  const auto target = ConvexHull::Status::VerticesCompleted; //< to get vertices
  // the array `indices` provides the initial full dimensional simplex,
  // if we had the last point.
  indices.push_back( input_points.size() - 1 );
  bool ok_init = hull.setInitialSimplex( indices );
  if ( ! ok_init )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Weird error in hull.setInitialSimplex" << std::endl;
      return PointRange();
    }
  bool ok_hull = hull.computeConvexHull( target );
  if ( ! ok_hull )
    {
      std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
                << "Weird error in hull.computeConvexHull" << std::endl;
      return PointRange();
    }
 
  // old code
  // 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 ( std::size_t i = 1; i < input_points.size(); i++ )
    {
      const auto& p = input_points[ i ];
      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;
  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;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  std::cout << "[ConvexityHelper::computeConvexHullBoundary(SurfMesh)] status=" << (int) hull.status() << "\n";
  if ( !ok ) return false;
  std::vector< RealPoint > positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  for ( auto p : positions ) std::cout << " " << p;
  for ( auto f : faces ) {
    std::cout << " (";
    for ( auto i : f ) std::cout << " " << i;
    std::cout << " )";
  }
  std::cout << "\n";
  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;
  ConvexHull hull;
  hull.setInput( input_points, remove_duplicates );
  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
  std::cout << "[ConvexityHelper::computeConvexHullBoundary(PolySurf)] status=" << (int) hull.status() << "\n";
  if ( !ok ) return false;
  PointRange positions;
  hull.getVertexPositions( positions );
  std::vector< IndexRange > faces;
  hull.getFacetVertices( faces );
  for ( auto p : positions ) std::cout << " " << p;
  for ( auto f : faces ) {
    std::cout << " (";
    for ( auto i : f ) std::cout << " " << i;
    std::cout << " )";
  }
  std::cout << "\n";
  // 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 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 ( Index i = 0; i < faces.size(); i++ )
    {
      // 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;
  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;
  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 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 ( Index i = 0; i < faces.size(); i++ )
    {
      // 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++;
      }
    }
  }
}
 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////