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;
    typedef DGtal::BigInteger               BigInteger;
    static const Dimension dimension = Space::dimension;

    struct UnitSegment {
      Dimension k;
      UnitSegment( Dimension d ) : k( d ) {}
    };
    struct UnitSegment {…};

    struct RightStrictUnitSegment {
      Dimension k;
      RightStrictUnitSegment( Dimension d ) : k( d ) {}
    };
    struct RightStrictUnitSegment {…};

    struct LeftStrictUnitSegment {
      Dimension k;
      LeftStrictUnitSegment( Dimension d ) : k( d ) {}
    };
    struct LeftStrictUnitSegment {…};

    struct UnitCell {
      std::vector<Dimension> dims;
      UnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
      UnitCell( std::initializer_list<Dimension> l ) {…}

      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;
      }
      friend std::ostream& {…}
    };
    struct UnitCell {…};
    
    struct RightStrictUnitCell {
      std::vector<Dimension> dims;
      RightStrictUnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
      RightStrictUnitCell( std::initializer_list<Dimension> l ) {…}

      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;
      }
      friend std::ostream& {…}
    };
    struct RightStrictUnitCell {…};

    struct LeftStrictUnitCell {
      std::vector<Dimension> dims;
      LeftStrictUnitCell( std::initializer_list<Dimension> l )
        : dims( l.begin(), l.end() ) {}
      LeftStrictUnitCell( std::initializer_list<Dimension> l ) {…}

      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;
      }
      friend std::ostream& {…}
    };
    struct LeftStrictUnitCell {…};


    ~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();

    BoundedLatticePolytope interiorPolytope() const;

    BoundedLatticePolytope closurePolytope() const;
    
 
    // ----------------------- 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:

 
    
    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 getKPoints( std::vector<Point>& pts, const Point& alpha_shift ) const;

    void getInteriorPoints( std::vector<Point>& pts ) const;

    void getBoundaryPoints( std::vector<Point>& pts ) const;

    template <typename PointSet>
    void insertPoints( PointSet& pts_set ) const;

    template <typename PointSet>
    void insertKPoints( PointSet& pts_set, const Point& alpha_shift ) const;
    
 
    // -------------- Enumeration services (old methods by scanning ) --------------
  public:

    
    Integer countByScanning() const;

    Integer countInteriorByScanning() const;

    Integer countBoundaryByScanning() const;

    Integer countWithinByScanning( Point low, Point hi ) const;

    Integer countUpToByScanning( Integer max ) const;

    void getPointsByScanning( std::vector<Point>& pts ) const;

    void getInteriorPointsByScanning( std::vector<Point>& pts ) const;

    void getBoundaryPointsByScanning( std::vector<Point>& pts ) const;

    template <typename PointSet>
    void insertPointsByScanning( 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
  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;
      }
      addEdgeConstraint( Polytope& , unsigned int , unsigned int , {…}

      static
      Vector crossProduct( const Vector& , const Vector& )
      {
        trace.error() << "[BoundedLatticePolytopeHelper::crossProduct]"
                      << " this method is only implemented in 3D." << std::endl;
        return Vector::zero;
      }
      Vector crossProduct( const Vector& , const Vector& ) {…}
    };
    struct BoundedLatticePolytopeSpecializer {…};
    
    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 );
              }
            }
      }
      addEdgeConstraint( Polytope& P, unsigned int i, unsigned int j, {…}

      static
      Vector crossProduct( const Vector& v1, const Vector& v2 )
      {
        return v1.crossProduct( v2 );
      }
      Vector crossProduct( const Vector& v1, const Vector& v2 ) {…}
    };
    struct BoundedLatticePolytopeSpecializer<3, TInteger> {…};
  }

  
  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)