BoundedRationalPolytope.h

 
#pragma once
 
#if defined(BoundedRationalPolytope_RECURSES)
#error Recursive header files inclusion detected in BoundedRationalPolytope.h
#else // defined(BoundedRationalPolytope_RECURSES)
#define BoundedRationalPolytope_RECURSES
 
#if !defined BoundedRationalPolytope_h
#define BoundedRationalPolytope_h
 
// Inclusions
#include <iostream>
#include <list>
#include <vector>
#include <string>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/CSpace.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/arithmetic/IntegerComputer.h"
#include "DGtal/arithmetic/ClosedIntegerHalfPlane.h"
 
namespace DGtal
{
 
  // template class BoundedRationalPolytope
  template < typename TSpace >
  class BoundedRationalPolytope 
  {
    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
 
  public:
    typedef BoundedRationalPolytope<TSpace>  Self;
    typedef TSpace                          Space;
    typedef typename Space::Integer         Integer;
    typedef typename Space::Point           Point;
    typedef typename Space::Vector          Vector;
    typedef std::vector<Vector>             InequalityMatrix;
    typedef std::vector<Integer>            InequalityVector;
    typedef HyperRectDomain< Space >        Domain; 
    typedef ClosedIntegerHalfPlane< Space > HalfSpace;
#ifdef WITH_BIGINTEGER
    typedef DGtal::BigInteger               BigInteger;
#else
    typedef DGtal::int64_t                  BigInteger;
#endif
    static const Dimension dimension = Space::dimension;
 
    struct UnitSegment {
      Dimension k;
      UnitSegment( Dimension d ) : k( d ) {}
    };
 
    struct UnitCell {
      std::vector<Dimension> dims;
      UnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
 
      friend std::ostream&
      operator<< ( std::ostream & out, 
                   const UnitCell & object )
      {
        out << "{";
        for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ];
        out << "}";
        return out;
      }
    };
    
 
    struct Rational {
      Integer p; // numerator
      Integer q; // denominator
      inline Rational( Integer a, Integer b ) : p( a ), q( b ) {}
    };
    
 
    ~BoundedRationalPolytope() = default;
 
    BoundedRationalPolytope();
 
    BoundedRationalPolytope ( const Self & other ) = default;
 
    
    BoundedRationalPolytope( std::initializer_list<Point> l );
    
    template <typename PointIterator>
    BoundedRationalPolytope( Integer d, PointIterator itB, PointIterator itE );
 
    template <typename HalfSpaceIterator>
    BoundedRationalPolytope( Integer d,
                             const Domain& domain,
                             HalfSpaceIterator itB, HalfSpaceIterator itE,
                             bool valid_edge_constraints = false,
                             bool check_duplicate_constraints = false );
 
    template <typename HalfSpaceIterator>
    void init( Integer d, const Domain& domain,
               HalfSpaceIterator itB, HalfSpaceIterator itE,
               bool valid_edge_constraints = false,
               bool check_duplicate_constraints = false );
    
    template <typename PointIterator>
    bool init( Integer d, PointIterator itB, PointIterator itE );
    
    Self & operator= ( const Self & other ) = default;
 
    void clear();
    
 
    // ----------------------- Accessor services ------------------------------
  public:
    
    const Domain& getDomain() const;
 
    const Domain& getLatticeDomain() const;
 
    const Domain& getRationalDomain() const;
 
    unsigned int nbHalfSpaces() const;
 
    Integer denominator() const;
    
    const Vector& getA( unsigned int i ) const;
 
    Integer getB( unsigned int i ) const;
 
    bool isLarge( unsigned int i ) const;
 
    const InequalityMatrix& getA() const;
    
    const InequalityVector& getB() const;
    const std::vector<bool>& getI() const;
 
    bool canBeSummed() const;
    
    
    // ----------------------- Check point services ------------------------------
  public:
 
 
    bool isInside( const Point& p ) const;
 
    bool isDomainPointInside( const Point& p ) const;
 
    bool isInterior( const Point& p ) const;
 
    bool isBoundary( const Point& p ) const;
 
    
    // ----------------------- Modification services ------------------------------
  public:
 
 
    BoundedRationalPolytope interiorPolytope() const;
      
    unsigned int cut( Dimension k, bool pos, Integer b, bool large = true );
 
    unsigned int cut( const Vector& a, Integer b, bool large = true,
                      bool valid_edge_constraint = false );
 
    unsigned int cut( const HalfSpace & hs, bool large = true,
                      bool valid_edge_constraint = false );
    
    void swap( BoundedRationalPolytope & other );
 
 
    Self& operator*=( Integer t );
 
    Self& operator*=( Rational r );
 
    Self& operator+=( UnitSegment s );
 
    Self& operator+=( UnitCell c );
 
 
    // ----------------------- Enumeration services ------------------------------
  public:
 
    
    Integer count() const;
 
    Integer countInterior() const;
 
    Integer countBoundary() const;
 
    Integer countWithin( Point low, Point hi ) const;
 
    Integer countUpTo( Integer max ) const;
 
    void getPoints( std::vector<Point>& pts ) const;
 
    void getInteriorPoints( std::vector<Point>& pts ) const;
 
    void getBoundaryPoints( std::vector<Point>& pts ) const;
 
    template <typename PointSet>
    void insertPoints( PointSet& pts_set ) const;
 
    
    
    // ----------------------- Interface --------------------------------------
  public:
 
    void selfDisplay ( std::ostream & out ) const;
 
    bool isValid() const;
 
    std::string className() const;
 
 
    // ------------------------- Protected Datas ------------------------------
  protected:
    // Denominator for constraints, i.e. \f$ q A x \le B \f$.
    Integer           q;
    // The matrix A in the polytope representation \f$ q A x \le B \f$.
    InequalityMatrix  A;
    // The vector B in the polytope representation \f$ q A x \le B \f$.
    InequalityVector  B;
    // Tight bounded box
    Domain            rationalD;
    // Lattice bounded box (i.e. floor( rationalD/q ))
    Domain            latticeD;
    // Are inequalities large ?
    std::vector<bool> I;
    // Indicates if Minkowski sums with segments will be valid
    bool myValidEdgeConstraints;
 
    // ------------------------- Private Datas --------------------------------
  private:
 
 
    // ------------------------- Internals ------------------------------------
  private:
    bool internalInitFromTriangle3D( Point a, Point b, Point c );
 
    bool internalInitFromSegment3D( Point a, Point b );
 
    bool internalInitFromSegment2D( Point a, Point b );
 
    Domain computeLatticeDomain( const Domain& d );
 
    Domain computeRationalDomain( const Domain& d );
    
  }; // end of class BoundedRationalPolytope
 
  namespace detail {
    template <DGtal::Dimension N, typename TInteger>
    struct BoundedRationalPolytopeSpecializer {
      typedef TInteger                        Integer;
      typedef SpaceND< N, Integer>            Space;
      typedef typename Space::Point           Point;
      typedef typename Space::Vector          Vector;
      typedef BoundedRationalPolytope< Space > Polytope;
      static const Dimension dimension = Space::dimension;
 
      static void 
      addEdgeConstraint( Polytope& , unsigned int , unsigned int ,
                         const std::vector<Point>& )
      {
        trace.error() << "[BoundedRationalPolytopeHelper::addEdgeConstraint]"
                      << " this method is only implemented in 3D." << std::endl;
      }
 
      static
      Vector crossProduct( const Vector& , const Vector& )
      {
        trace.error() << "[BoundedRationalPolytopeHelper::crossProduct]"
                      << " this method is only implemented in 3D." << std::endl;
        return Vector::zero;
      }
    };
    
    template <typename TInteger>
    struct BoundedRationalPolytopeSpecializer<3, TInteger> {
      typedef TInteger                        Integer;
      typedef SpaceND< 3, Integer>            Space;
      typedef typename Space::Point           Point;
      typedef typename Space::Vector          Vector;
      typedef BoundedRationalPolytope< Space > Polytope;
      static const Dimension dimension = Space::dimension;
 
      static void 
      addEdgeConstraint( Polytope& P, unsigned int i, unsigned int j,
                         const std::vector<Point>& pts )
      {
        Vector ab = pts[ i ] - pts[ j ];
        for ( int s = 0; s < 2; s++ )
          for ( Dimension k = 0; k < dimension; ++k )
            {
              Vector  n = ab.crossProduct( Point::base( k, (s == 0) ? 1 : -1 ) );
              Integer b = n.dot( pts[ i ] );
              std::size_t nb_in = 0;
              for ( auto p : pts ) {
                Integer v = n.dot( p );
                if ( v < b )  nb_in++;
              }
              if ( nb_in == dimension - 1 ) {
                P.cut( n, b, true, true );
              }
            }
      }
      static
      Vector crossProduct( const Vector& v1, const Vector& v2 )
      {
        return v1.crossProduct( v2 );
      }
    };
  }
 
  
  template <typename TSpace>
  std::ostream&
  operator<< ( std::ostream & out, 
               const BoundedRationalPolytope<TSpace> & object );
 
 
  template <typename TSpace>
  BoundedRationalPolytope<TSpace>
  operator* ( typename BoundedRationalPolytope<TSpace>::Integer t, 
              const BoundedRationalPolytope<TSpace> & P );
 
  template <typename TSpace>
  BoundedRationalPolytope<TSpace>
  operator* ( typename BoundedRationalPolytope<TSpace>::Rational r, 
              const BoundedRationalPolytope<TSpace> & P );
    
 
  template <typename TSpace>
  BoundedRationalPolytope<TSpace>
  operator+ ( const BoundedRationalPolytope<TSpace> & P,
              typename BoundedRationalPolytope<TSpace>::UnitSegment s );
 
  template <typename TSpace>
  BoundedRationalPolytope<TSpace>
  operator+ ( const BoundedRationalPolytope<TSpace> & P,
              typename BoundedRationalPolytope<TSpace>::UnitCell c );
 
  
} // namespace DGtal
 
 
// Includes inline functions.
#include "BoundedRationalPolytope.ih"
 
//                                                                           //
 
#endif // !defined BoundedRationalPolytope_h
 
#undef BoundedRationalPolytope_RECURSES
#endif // else defined(BoundedRationalPolytope_RECURSES)