BoundedLatticePolytope.h

 
#pragma once
 
#if defined(BoundedLatticePolytope_RECURSES)
#error Recursive header files inclusion detected in BoundedLatticePolytope.h
#else // defined(BoundedLatticePolytope_RECURSES)
#define BoundedLatticePolytope_RECURSES
 
#if !defined BoundedLatticePolytope_h
#define BoundedLatticePolytope_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 BoundedLatticePolytope
  template < typename TSpace >
  class BoundedLatticePolytope 
  {
    BOOST_CONCEPT_ASSERT(( concepts::CSpace< TSpace > ));
 
  public:
    typedef BoundedLatticePolytope<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 RightStrictUnitSegment {
      Dimension k;
      RightStrictUnitSegment( Dimension d ) : k( d ) {}
    };
 
    struct LeftStrictUnitSegment {
      Dimension k;
      LeftStrictUnitSegment( 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 RightStrictUnitCell {
      std::vector<Dimension> dims;
      RightStrictUnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
 
      friend std::ostream&
      operator<< ( std::ostream & out, 
                   const RightStrictUnitCell & object )
      {
        out << "{";
        for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ];
        out << "}";
        return out;
      }
    };
 
    struct LeftStrictUnitCell {
      std::vector<Dimension> dims;
      LeftStrictUnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
 
      friend std::ostream&
      operator<< ( std::ostream & out, 
                   const LeftStrictUnitCell & object )
      {
        out << "{";
        for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ];
        out << "}";
        return out;
      }
    };
 
 
    ~BoundedLatticePolytope() = default;
 
    BoundedLatticePolytope() = default;
 
    BoundedLatticePolytope ( const Self & other ) = default;
 
    
    BoundedLatticePolytope( std::initializer_list<Point> l );
    
    template <typename PointIterator>
    BoundedLatticePolytope( PointIterator itB, PointIterator itE );
 
    template <typename HalfSpaceIterator>
    BoundedLatticePolytope( const Domain& domain,
                            HalfSpaceIterator itB, HalfSpaceIterator itE,
                            bool valid_edge_constraints = false,
                            bool check_duplicate_constraints = false );
 
    template <typename HalfSpaceIterator>
    void init( const Domain& domain,
               HalfSpaceIterator itB, HalfSpaceIterator itE,
               bool valid_edge_constraints = false,
               bool check_duplicate_constraints = false );
 
    
    template <typename PointIterator>
    bool init( PointIterator itB, PointIterator itE );
    
    Self & operator= ( const Self & other ) = default;
 
    void clear();
    
 
    // ----------------------- Accessor services ------------------------------
  public:
    
    const Domain& getDomain() const;
 
    unsigned int nbHalfSpaces() 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:
 
 
    BoundedLatticePolytope 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( BoundedLatticePolytope & other );
 
 
    Self& operator*=( Integer t );
 
    Self& operator+=( UnitSegment s );
 
    Self& operator+=( UnitCell c );
 
    Self& operator+=( RightStrictUnitSegment s );
 
    Self& operator+=( RightStrictUnitCell c );
 
    Self& operator+=( LeftStrictUnitSegment s );
 
    Self& operator+=( LeftStrictUnitCell 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:
    InequalityMatrix  A;
    InequalityVector  B;
    Domain            D;
    std::vector<bool> I;
    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 );
 
  }; // end of class BoundedLatticePolytope
 
  namespace detail {
    template <DGtal::Dimension N, typename TInteger>
    struct BoundedLatticePolytopeSpecializer {
      typedef TInteger                        Integer;
      typedef SpaceND< N, Integer>            Space;
      typedef typename Space::Point           Point;
      typedef typename Space::Vector          Vector;
      typedef BoundedLatticePolytope< Space > Polytope;
      static const Dimension dimension = Space::dimension;
 
      static void 
      addEdgeConstraint( Polytope& , unsigned int , unsigned int ,
                         const std::vector<Point>& )
      {
        trace.error() << "[BoundedLatticePolytopeHelper::addEdgeConstraint]"
                      << " this method is only implemented in 3D." << std::endl;
      }
 
      static
      Vector crossProduct( const Vector& , const Vector& )
      {
        trace.error() << "[BoundedLatticePolytopeHelper::crossProduct]"
                      << " this method is only implemented in 3D." << std::endl;
        return Vector::zero;
      }
    };
    
    template <typename TInteger>
    struct BoundedLatticePolytopeSpecializer<3, TInteger> {
      typedef TInteger                        Integer;
      typedef SpaceND< 3, Integer>            Space;
      typedef typename Space::Point           Point;
      typedef typename Space::Vector          Vector;
      typedef BoundedLatticePolytope< 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 == pts.size() - 2 ) {
                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 BoundedLatticePolytope<TSpace> & object );
 
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator* ( typename BoundedLatticePolytope<TSpace>::Integer t, 
              const BoundedLatticePolytope<TSpace> & P );
    
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::UnitSegment s );
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::UnitCell c );
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::RightStrictUnitSegment s );
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::RightStrictUnitCell c );
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::LeftStrictUnitSegment s );
 
  template <typename TSpace>
  BoundedLatticePolytope<TSpace>
  operator+ ( const BoundedLatticePolytope<TSpace> & P,
              typename BoundedLatticePolytope<TSpace>::LeftStrictUnitCell c );
 
  
} // namespace DGtal
 
 
// Includes inline functions.
#include "BoundedLatticePolytope.ih"
 
//                                                                           //
 
#endif // !defined BoundedLatticePolytope_h
 
#undef BoundedLatticePolytope_RECURSES
#endif // else defined(BoundedLatticePolytope_RECURSES)