GenericLatticeConvexHull.h

 
#pragma once
 
#if defined(GenericLatticeConvexHull_RECURSES)
#error Recursive header files inclusion detected in GenericLatticeConvexHull.h
#else // defined(GenericLatticeConvexHull_RECURSES)
#define GenericLatticeConvexHull_RECURSES
 
#if !defined GenericLatticeConvexHull_h
#define GenericLatticeConvexHull_h
 
// Inclusions
#include <iostream>
#include <string>
#include <vector>
#include <queue>
#include <set>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/geometry/tools/AffineGeometry.h"
#include "DGtal/geometry/tools/AffineBasis.h"
#include "DGtal/geometry/tools/QuickHull.h"
#include "DGtal/geometry/tools/QuickHullKernels.h"
#include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
 
namespace DGtal
{
  // Forward declaration.
  template < Dimension dim,
             typename TCoordinateInteger,
             typename TInternalInteger >
  struct GenericLatticeConvexHull;
 
  namespace detail {
    template < Dimension dim,
               typename TCoordinateInteger,
               typename TInternalInteger,
               Dimension K >
    struct GenericLatticeConvexHullComputers
    {
      typedef ConvexHullIntegralKernel
      < K,TCoordinateInteger,TInternalInteger > Kernel;
      typedef detail::GenericLatticeConvexHullComputers
      < dim, TCoordinateInteger, TInternalInteger, K-1> LowerKernels;
      typedef std::size_t                Size;
      typedef typename Kernel::CoordinatePoint Point;
      typedef typename Point::Coordinate Integer;
      typedef QuickHull< Kernel >        QHull;
      typedef SpaceND< K, Integer >      Space;
      typedef BoundedLatticePolytope< Space > LatticePolytope;
      typedef GenericLatticeConvexHull< dim,
                                        TCoordinateInteger,
                                        TInternalInteger > Computer;
      typedef typename Computer::OutputPoint OutputPoint;
      static const Dimension dimension = K;
 
      GenericLatticeConvexHullComputers( Computer* ptrGenQHull )
        : ptr_gen_qhull( ptrGenQHull ), lower_kernels( ptrGenQHull ),
          hull( Kernel(), ptrGenQHull->debug_level )
      {
        clear();
      }
 
      void clear()
      {
        hull.clear(); 
        proj_points.clear();
        polytope.clear();
        lower_kernels.clear();
      }
 
      template <typename TInputPoint>
      bool compute( const std::vector< Size >& I,
                    const std::vector< TInputPoint >& X )
      {
        typedef TInputPoint InputPoint;
        typedef AffineGeometry< InputPoint > Affine;
        typedef AffineBasis< InputPoint >    Basis;
        hull.clear();
        polytope.clear();
        if ( (I.size()-1) != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.compute( I, X );
          }
        ptr_gen_qhull->affine_dimension = dimension;
        auto& points    = ptr_gen_qhull->points;
        auto& ppoints   = ptr_gen_qhull->projected_points;
        auto& positions = ptr_gen_qhull->positions;
        auto& v2p       = ptr_gen_qhull->vertex2point;
        auto& facets    = ptr_gen_qhull->facets;
        auto& dilation  = ptr_gen_qhull->projected_dilation;
        auto& aff_basis = ptr_gen_qhull->affine_basis;
        
        Basis basis;
        if ( dimension != ptr_gen_qhull->dimension )
          { // Build the affine basis spanning the convex hull affine space.
            if ( dimension == ptr_gen_qhull->dimension-1 )
              { // codimension 1
                // builds orthogonal vector
                InputPoint normal;
                functions::getOrthogonalVector( normal, X, I );
                basis = Basis( X[0], normal, Basis::Type::ECHELON_REDUCED );
              }
            else
              { // Generic case
                basis = Basis( X, Basis::Type::SHORTEST_ECHELON_REDUCED );
              }
          }
        // Build projected points on affine basis
        dilation  = basis.projectPoints( proj_points, X );
          
        // Compute convex hull using quickhull.
        bool ok_input = hull.setInput( proj_points, false );
        bool ok_hull  = hull.computeConvexHull( QHull::Status::VerticesCompleted );
        if ( ! ok_hull || ! ok_input )
          {
            trace.error() << "[GenericLatticeConvexHullComputers::compute]"
                          << " Error in quick hull computation.\n"
                          << "qhull=" << hull << "\n";
            return false;
          }
        points.resize( X.size() );
        for ( Size i = 0; i < points.size(); i++ )
          points[ i ] = Affine::transform( X[ i ] );
        hull.getVertex2Point( v2p );
        hull.getFacetVertices( facets );
        positions.resize( v2p.size() );
        for ( Size i = 0; i < positions.size(); i++ )
          positions[ i ] = X[ v2p[ i ] ];
        ppoints.resize( proj_points.size() );
        for ( Size i = 0; i < ppoints.size(); i++ )
          {
            ppoints[ i ] = OutputPoint::zero;
            for ( Dimension j = 0; j < Point::dimension; j++ )
              ppoints[ i ][ j ] = proj_points[ i ][ j ];
          }
        aff_basis = AffineBasis<OutputPoint>( basis.origin(),
                                              basis.basis(),
                                              Basis::Type::SHORTEST_ECHELON_REDUCED,
                                              true );
        return true;
      }
 
      bool makePolytope()
      {
        typedef typename LatticePolytope::Domain     Domain;
        typedef typename LatticePolytope::HalfSpace  PolytopeHalfSpace;
        typedef typename QHull::HalfSpace            ConvexHullHalfSpace;
        typedef AffineGeometry< Point > Affine;
        if ( ptr_gen_qhull->affine_dimension != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.makePolytope();
          }
        // If polytope is already initialized returns.
        if ( polytope.nbHalfSpaces() > 0 ) return true;
        // Compute domain
        Point l = proj_points[ 0 ];
        Point u = proj_points[ 0 ];
        for ( std::size_t i = 1; i < proj_points.size(); i++ )
          {
            const auto& p = proj_points[ i ];
            l = l.inf( p );
            u = u.sup( p );
          }
        Domain domain( l, u );
        
        // Initialize polytope
        std::vector< ConvexHullHalfSpace > HS;
        std::vector< PolytopeHalfSpace >   PHS;
        hull.getFacetHalfSpaces( HS );
        PHS.reserve( HS.size() );
        for ( auto& H : HS ) {
          Point  N;
          Integer nu;
          for ( Dimension i = 0; i < dimension; ++i )
            N[ i ] = IntegerConverter< dimension, Integer >
              ::cast( H.internalNormal()[ i ]  );
          Point   Ns = Affine::simplifiedVector( N );
          Integer g  = N.norm1() / Ns.norm1();
          nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
          PHS.emplace_back( Ns, nu / g );
        }
        polytope = LatticePolytope( domain, PHS.cbegin(), PHS.cend(), false, true );
        return polytope.isValid();
      }
 
      Integer count()
      {
        if ( ptr_gen_qhull->affine_dimension != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.count();
          }
        // If polytope is not initialized returns error.
        if ( polytope.nbHalfSpaces() == 0 ) return -1;
        return polytope.count();
      }
 
      Integer countInterior()
      {
        if ( ptr_gen_qhull->affine_dimension != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.countInterior();
          }
        // If polytope is not initialized returns error.
        if ( polytope.nbHalfSpaces() == 0 ) return -1;
        return polytope.countInterior();
      }
      
      Integer countBoundary()
      {
        if ( ptr_gen_qhull->affine_dimension != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.countBoundary();
          }
        // If polytope is not initialized returns error.
        if ( polytope.nbHalfSpaces() == 0 ) return -1;
        return polytope.countBoundary();
      }
      
      Integer countUpTo( Integer max )
      {
        if ( ptr_gen_qhull->affine_dimension != dimension )
          { // This kernel is not adapted => go to lower dimension
            return lower_kernels.countUpTo( max );
          }
        // If polytope is not initialized returns error.
        if ( polytope.nbHalfSpaces() == 0 ) return -1;
        return polytope.countUpTo( max );
      }
 
      
      Computer*            ptr_gen_qhull; 
      LowerKernels         lower_kernels; 
      QHull                hull; 
      std::vector< Point > proj_points; 
      LatticePolytope      polytope; 
    };
 
    
    template < Dimension dim,
               typename TCoordinateInteger,
               typename TInternalInteger >
    struct GenericLatticeConvexHullComputers< dim, TCoordinateInteger, TInternalInteger, 1>
    {
      typedef ConvexHullIntegralKernel
      < 1,TCoordinateInteger,TInternalInteger > Kernel;
      typedef Kernel                     Type;
      typedef std::size_t                Size;
      typedef typename Kernel::CoordinatePoint Point;
      typedef typename Point::Coordinate Integer;
      typedef SpaceND<1, Integer>        Space;
      typedef GenericLatticeConvexHull< dim,
                                        TCoordinateInteger,
                                        TInternalInteger > Computer;
      typedef typename Computer::OutputPoint OutputPoint;
      static const Dimension             dimension = 1;
 
      GenericLatticeConvexHullComputers( Computer* ptrGenQHull )
        : ptr_gen_qhull( ptrGenQHull )
      {
        clear();
      }
 
      void clear()
      {
        proj_points.clear();
      }
      
      template <typename TInputPoint>
      bool compute( const std::vector< Size >& I,
                    const std::vector< TInputPoint >& X )
      {
        typedef TInputPoint InputPoint;
        typedef AffineGeometry< InputPoint > Affine;
        typedef AffineBasis< InputPoint >    Basis;
        typedef std::size_t                  Index;
 
        auto& aff_dim   = ptr_gen_qhull->affine_dimension;
        auto& points    = ptr_gen_qhull->points;
        auto& ppoints   = ptr_gen_qhull->projected_points;
        auto& positions = ptr_gen_qhull->positions;
        auto& v2p       = ptr_gen_qhull->vertex2point;
        auto& facets    = ptr_gen_qhull->facets;
        auto& dilation  = ptr_gen_qhull->projected_dilation;
        auto& aff_basis = ptr_gen_qhull->affine_basis;
        facets.clear(); // no facets
        
        if ( (I.size()-1) != dimension )
          { // This kernel is not adapted => lower dimension is either
            // 0, ie. 1 point, or -1, ie. 0 points.
            if ( ! X.empty() )
              {
                aff_dim = 0;
                points.resize( X.size() );
                for ( Size i = 0; i < points.size(); i++ )
                  points[ i ] = Affine::transform( X[ i ] );
                ppoints = points;
                positions.resize( 1 );
                positions[ 0 ] = X[ 0 ];
                v2p.resize( 1 );
                v2p[ 0 ]   = 0;
                nb_in_hull = 1;
              }
            else
              {
                aff_dim = -1;
                points.clear();
                ppoints.clear();
                positions.clear();
                v2p.clear();
                nb_in_hull = 0;
              }
            return true;
          }
        // Generic 1D case.
        aff_dim = dimension;
        Basis basis;
        if ( dimension != ptr_gen_qhull->dimension )
          {
            // Build the affine basis spanning the convex hull affine space.
            basis = Basis( X, Basis::Type::SHORTEST_ECHELON_REDUCED );
          }
        // Build projected points on affine basis
        dilation  = basis.projectPoints( proj_points, X );
        // Compute convex hull by looking at extremal points
        Index left  = 0;
        Index right = 0;
        for ( Index i = 1; i < proj_points.size(); i++ )
          {
            if ( proj_points[ i ][ 0 ] < proj_points[ left ][ 0 ] )
              left = i;
            else if ( proj_points[ i ][ 0 ] > proj_points[ right ][ 0 ] )
              right = i;
          }
        points.resize( X.size() );
        for ( Size i = 0; i < points.size(); i++ )
          points[ i ] = Affine::transform( X[ i ] );
        ppoints.resize( proj_points.size() );
        for ( Size i = 0; i < ppoints.size(); i++ )
          {
            ppoints[ i ]      = OutputPoint::zero;
            ppoints[ i ][ 0 ] = proj_points[ i ][ 0 ];
          }
        v2p.resize( 2 );
        v2p[ 0 ] = left;
        v2p[ 1 ] = right;
        positions.resize( 2 );
        positions[ 0 ] = X[ v2p[ 0 ] ];
        positions[ 1 ] = X[ v2p[ 1 ] ];
        aff_basis = AffineBasis<OutputPoint>( basis.origin(),
                                              basis.basis(),
                                              Basis::Type::SHORTEST_ECHELON_REDUCED,
                                              true );
        auto    dx  = Affine::transform( points[ v2p[ 1 ] ] - points[ v2p[ 0 ] ] );
        auto    sdx = Affine::simplifiedVector( dx );
        Integer n   = dx.normInfinity() / sdx.normInfinity();
        nb_in_hull  = n+1;
        return true;
      }
 
      bool makePolytope()
      {
        return true;
      }
 
      Integer count() const
      {
        return nb_in_hull;
      }
 
      Integer countInterior()
      {
        return nb_in_hull >= 2 ? nb_in_hull - 2 : 0;
      }
      
      Integer countBoundary()
      {
        return nb_in_hull >= 2 ? 2 : nb_in_hull;
      }
      
      Integer countUpTo( Integer max )
      {
        return nb_in_hull < max ? nb_in_hull : max;
      }
 
      Computer*            ptr_gen_qhull; 
      std::vector< Point > proj_points; 
      Integer              nb_in_hull; 
    };
  }
  
  // template class GenericLatticeConvexHull
 
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct GenericLatticeConvexHull
  {
    typedef ConvexHullIntegralKernel< dim,
                                      TCoordinateInteger,
                                      TInternalInteger > Kernel;
    typedef typename Kernel::CoordinatePoint     Point;
    typedef typename Kernel::CoordinateScalar    Integer;
    typedef typename Kernel::InternalScalar      InternalInteger;
    typedef Point                                OutputPoint;
    typedef std::size_t                Index;
    typedef std::size_t                Size;
    typedef std::vector< Index >       IndexRange;
    typedef detail::GenericLatticeConvexHullComputers
    < dim, TCoordinateInteger, TInternalInteger, dim > GenericComputers;
 
    static const Size  dimension  = Point::dimension;
 
    
    // ----------------------- standard services --------------------------
  public:
  
    GenericLatticeConvexHull( const Kernel& K_ = Kernel(), int dbg = 0 )
      : kernel( K_ ), debug_level( dbg ), generic_computers( this )
    {
      clear();
    }
 
    void clear()
    {
      generic_computers.clear();
      points.clear();
      projected_points.clear();
      affine_dimension  = -1;
      positions.clear();
      facets.clear();
      vertex2point.clear();
      affine_basis       = AffineBasis< OutputPoint >();
      projected_dilation = 1;
      polytope_computed  = false;
    }
    
    // -------------------------- Convex hull services ----------------------------
  public:
 
    
    template < typename InputPoint >
    bool compute( const std::vector< InputPoint >& input_points )
    {
      // Determine affine dimension of set of input points.
      typedef AffineGeometry< InputPoint > Affine;
      std::vector< Size > indices = Affine::affineSubset( input_points );
      bool ok = generic_computers.compute( indices, input_points );
      if ( ( ! ok ) || ( debug_level >= 1 ) )
        {
          std::cout << "Generic Convex hull #V=" << positions.size()
                    << " #F=" << facets.size() << "\n";
          for ( Size i = 0; i < facets.size(); i++ )
            {
              std::cout << "F_" << i << " = (";
              for ( auto v : facets[ i ] ) std::cout << " " << v;
              std::cout << " )\n";
            }
          for ( Size i = 0; i < positions.size(); i++ )
            std::cout << "V_" << i
                      << " pi(x)=" << projected_points[ vertex2point[ i ] ]
                      << " -> x=" << positions[ i ] << "\n";
        }
      return ok;
    }
 
    Integer count()
    {
      if ( ! polytope_computed )
        polytope_computed = generic_computers.makePolytope();
      if ( ! polytope_computed ) return -1; 
      return generic_computers.count();
    }
 
    Integer countInterior()
    {
      if ( ! polytope_computed )
        polytope_computed = generic_computers.makePolytope();
      if ( ! polytope_computed ) return -1; 
      return generic_computers.countInterior();
      }
    
    Integer countBoundary()
    {
      if ( ! polytope_computed )
        polytope_computed = generic_computers.makePolytope();
      if ( ! polytope_computed ) return -1; 
      return generic_computers.countBoundary();
    }
    
    Integer countUpTo( Integer max )
    {
      if ( ! polytope_computed )
        polytope_computed = generic_computers.makePolytope();
      if ( ! polytope_computed ) return -1; 
      return generic_computers.countUpTo( max );
    }
 
    
    
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const
    {
      out << "[GenericLatticeConvexHull"
          << " dim=" << dimension
          << " #in=" << points.size()
          << " aff_dim=" << affine_dimension
          << " #V=" << positions.size()
          << " #F=" << facets.size()
          << "]";
    }
  
    bool isValid() const
    {
      return affine_dimension >= 0;
    }
 
    // ------------------------ public datas --------------------------
  public:
  
  public:
    Kernel kernel;
    int debug_level; 
    GenericComputers generic_computers;
 
    std::vector< OutputPoint > points;
    std::vector< OutputPoint > projected_points;
    int64_t                    affine_dimension;
    std::vector< OutputPoint > positions;
    std::vector< IndexRange >  facets;
    IndexRange                 vertex2point;
    Integer                    projected_dilation;
    AffineBasis< OutputPoint > affine_basis;
    bool                       polytope_computed { false };
    
  };
 
  template < Dimension dim,
             typename TCoordinateInteger,
             typename TInternalInteger >
  std::ostream&
  operator<< ( std::ostream & out,
               const GenericLatticeConvexHull< dim, TCoordinateInteger, TInternalInteger > & object )
  {
    object.selfDisplay( out );
    return out;
  }
  
} // namespace DGtal
 
// Includes inline functions.
//                                                                           //
 
#endif // !defined GenericLatticeConvexHull_h
 
#undef GenericLatticeConvexHull_RECURSES
#endif // else defined(GenericLatticeConvexHull_RECURSES)