DGtal  1.5.beta
BoundedRationalPolytope.h
1 
17 #pragma once
18 
31 #if defined(BoundedRationalPolytope_RECURSES)
32 #error Recursive header files inclusion detected in BoundedRationalPolytope.h
33 #else // defined(BoundedRationalPolytope_RECURSES)
35 #define BoundedRationalPolytope_RECURSES
36 
37 #if !defined BoundedRationalPolytope_h
39 #define BoundedRationalPolytope_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <list>
45 #include <vector>
46 #include <string>
47 #include "DGtal/base/Common.h"
48 #include "DGtal/kernel/CSpace.h"
49 #include "DGtal/kernel/domains/HyperRectDomain.h"
50 #include "DGtal/arithmetic/IntegerComputer.h"
51 #include "DGtal/arithmetic/ClosedIntegerHalfPlane.h"
53 
54 namespace DGtal
55 {
56 
58  // template class BoundedRationalPolytope
73  template < typename TSpace >
75  {
77 
78  public:
80  typedef TSpace Space;
81  typedef typename Space::Integer Integer;
82  typedef typename Space::Point Point;
83  typedef typename Space::Vector Vector;
84  typedef std::vector<Vector> InequalityMatrix;
85  typedef std::vector<Integer> InequalityVector;
88 #ifdef WITH_BIGINTEGER
90 #else
91  typedef DGtal::int64_t BigInteger;
92 #endif
93  static const Dimension dimension = Space::dimension;
94 
99  struct UnitSegment {
101  UnitSegment( Dimension d ) : k( d ) {}
102  };
103 
109  struct UnitCell {
110  std::vector<Dimension> dims;
111  UnitCell( std::initializer_list<Dimension> l )
112  : dims( l.begin(), l.end() ) {}
113 
120  friend std::ostream&
121  operator<< ( std::ostream & out,
122  const UnitCell & object )
123  {
124  out << "{";
125  for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ];
126  out << "}";
127  return out;
128  }
129  };
130 
131 
134  struct Rational {
135  Integer p; // numerator
136  Integer q; // denominator
140  inline Rational( Integer a, Integer b ) : p( a ), q( b ) {}
141  };
142 
145 
150 
155 
160  BoundedRationalPolytope ( const Self & other ) = default;
161 
162 
173  BoundedRationalPolytope( std::initializer_list<Point> l );
174 
188  template <typename PointIterator>
189  BoundedRationalPolytope( Integer d, PointIterator itB, PointIterator itE );
190 
217  template <typename HalfSpaceIterator>
219  const Domain& domain,
220  HalfSpaceIterator itB, HalfSpaceIterator itE,
221  bool valid_edge_constraints = false,
222  bool check_duplicate_constraints = false );
223 
250  template <typename HalfSpaceIterator>
251  void init( Integer d, const Domain& domain,
252  HalfSpaceIterator itB, HalfSpaceIterator itE,
253  bool valid_edge_constraints = false,
254  bool check_duplicate_constraints = false );
255 
273  template <typename PointIterator>
274  bool init( Integer d, PointIterator itB, PointIterator itE );
275 
281  Self & operator= ( const Self & other ) = default;
282 
284  void clear();
285 
287 
288  // ----------------------- Accessor services ------------------------------
289  public:
292 
294  const Domain& getDomain() const;
295 
298  const Domain& getLatticeDomain() const;
299 
302  const Domain& getRationalDomain() const;
303 
305  unsigned int nbHalfSpaces() const;
306 
311 
317  const Vector& getA( unsigned int i ) const;
318 
324  Integer getB( unsigned int i ) const;
325 
331  bool isLarge( unsigned int i ) const;
332 
334  const InequalityMatrix& getA() const;
335 
337  const InequalityVector& getB() const;
341  const std::vector<bool>& getI() const;
342 
347  bool canBeSummed() const;
348 
350 
351  // ----------------------- Check point services ------------------------------
352  public:
353 
356 
361  bool isInside( const Point& p ) const;
362 
369  bool isDomainPointInside( const Point& p ) const;
370 
375  bool isInterior( const Point& p ) const;
376 
381  bool isBoundary( const Point& p ) const;
382 
384 
385  // ----------------------- Modification services ------------------------------
386  public:
387 
390 
394 
407  unsigned int cut( Dimension k, bool pos, Integer b, bool large = true );
408 
426  unsigned int cut( const Vector& a, Integer b, bool large = true,
427  bool valid_edge_constraint = false );
428 
445  unsigned int cut( const HalfSpace & hs, bool large = true,
446  bool valid_edge_constraint = false );
447 
452  void swap( BoundedRationalPolytope & other );
453 
454 
461 
468 
476 
484 
486 
487  // ----------------------- Enumeration services ------------------------------
488  public:
489 
492 
500  Integer count() const;
501 
513 
525 
536  Integer countWithin( Point low, Point hi ) const;
537 
555 
564  void getPoints( std::vector<Point>& pts ) const;
565 
574  void getInteriorPoints( std::vector<Point>& pts ) const;
575 
584  void getBoundaryPoints( std::vector<Point>& pts ) const;
585 
596  template <typename PointSet>
597  void insertPoints( PointSet& pts_set ) const;
598 
600 
601 
602  // ----------------------- Interface --------------------------------------
603  public:
606 
611  void selfDisplay ( std::ostream & out ) const;
612 
619  bool isValid() const;
620 
624  std::string className() const;
625 
627 
628  // ------------------------- Protected Datas ------------------------------
629  protected:
630  // Denominator for constraints, i.e. \f$ q A x \le B \f$.
632  // The matrix A in the polytope representation \f$ q A x \le B \f$.
634  // The vector B in the polytope representation \f$ q A x \le B \f$.
636  // Tight bounded box
638  // Lattice bounded box (i.e. floor( rationalD/q ))
640  // Are inequalities large ?
641  std::vector<bool> I;
642  // Indicates if Minkowski sums with segments will be valid
644 
645  // ------------------------- Private Datas --------------------------------
646  private:
647 
648 
649  // ------------------------- Internals ------------------------------------
650  private:
658 
665 
672 
680 
687 
688  }; // end of class BoundedRationalPolytope
689 
690  namespace detail {
699  template <DGtal::Dimension N, typename TInteger>
701  typedef TInteger Integer;
703  typedef typename Space::Point Point;
704  typedef typename Space::Vector Vector;
707 
720  static void
721  addEdgeConstraint( Polytope& , unsigned int , unsigned int ,
722  const std::vector<Point>& )
723  {
724  trace.error() << "[BoundedRationalPolytopeHelper::addEdgeConstraint]"
725  << " this method is only implemented in 3D." << std::endl;
726  }
727 
730  static
731  Vector crossProduct( const Vector& , const Vector& )
732  {
733  trace.error() << "[BoundedRationalPolytopeHelper::crossProduct]"
734  << " this method is only implemented in 3D." << std::endl;
735  return Vector::zero;
736  }
737  };
738 
746  template <typename TInteger>
748  typedef TInteger Integer;
750  typedef typename Space::Point Point;
751  typedef typename Space::Vector Vector;
754 
764  static void
765  addEdgeConstraint( Polytope& P, unsigned int i, unsigned int j,
766  const std::vector<Point>& pts )
767  {
768  Vector ab = pts[ i ] - pts[ j ];
769  for ( int s = 0; s < 2; s++ )
770  for ( Dimension k = 0; k < dimension; ++k )
771  {
772  Vector n = ab.crossProduct( Point::base( k, (s == 0) ? 1 : -1 ) );
773  Integer b = n.dot( pts[ i ] );
774  std::size_t nb_in = 0;
775  for ( auto p : pts ) {
776  Integer v = n.dot( p );
777  if ( v < b ) nb_in++;
778  }
779  if ( nb_in == dimension - 1 ) {
780  P.cut( n, b, true, true );
781  }
782  }
783  }
788  static
789  Vector crossProduct( const Vector& v1, const Vector& v2 )
790  {
791  return v1.crossProduct( v2 );
792  }
793  };
794  }
795 
798 
805  template <typename TSpace>
806  std::ostream&
807  operator<< ( std::ostream & out,
808  const BoundedRationalPolytope<TSpace> & object );
809 
810 
816  template <typename TSpace>
820 
826  template <typename TSpace>
830 
831 
839  template <typename TSpace>
843 
851  template <typename TSpace>
855 
857 
858 } // namespace DGtal
859 
860 
862 // Includes inline functions.
863 #include "BoundedRationalPolytope.ih"
864 
865 // //
867 
868 #endif // !defined BoundedRationalPolytope_h
869 
870 #undef BoundedRationalPolytope_RECURSES
871 #endif // else defined(BoundedRationalPolytope_RECURSES)
Aim: Represents an nD rational polytope, i.e. a convex polyhedron bounded by vertices with rational c...
BoundedRationalPolytope< TSpace > Self
unsigned int cut(const Vector &a, Integer b, bool large=true, bool valid_edge_constraint=false)
Domain computeRationalDomain(const Domain &d)
Integer countUpTo(Integer max) const
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
const InequalityMatrix & getA() const
bool isLarge(unsigned int i) const
const InequalityVector & getB() const
bool internalInitFromSegment3D(Point a, Point b)
void swap(BoundedRationalPolytope &other)
bool init(Integer d, PointIterator itB, PointIterator itE)
void selfDisplay(std::ostream &out) const
Self & operator+=(UnitCell c)
BoundedRationalPolytope(std::initializer_list< Point > l)
Self & operator*=(Integer t)
void insertPoints(PointSet &pts_set) const
Self & operator=(const Self &other)=default
bool isBoundary(const Point &p) const
Integer getB(unsigned int i) const
void init(Integer d, const Domain &domain, HalfSpaceIterator itB, HalfSpaceIterator itE, bool valid_edge_constraints=false, bool check_duplicate_constraints=false)
bool isDomainPointInside(const Point &p) const
const std::vector< bool > & getI() const
void getInteriorPoints(std::vector< Point > &pts) const
Self & operator*=(Rational r)
void clear()
Clears the polytope.
const Domain & getRationalDomain() const
std::string className() const
unsigned int cut(Dimension k, bool pos, Integer b, bool large=true)
void getBoundaryPoints(std::vector< Point > &pts) const
bool internalInitFromTriangle3D(Point a, Point b, Point c)
bool isInside(const Point &p) const
unsigned int cut(const HalfSpace &hs, bool large=true, bool valid_edge_constraint=false)
unsigned int nbHalfSpaces() const
void getPoints(std::vector< Point > &pts) const
const Domain & getLatticeDomain() const
bool internalInitFromSegment2D(Point a, Point b)
BoundedRationalPolytope interiorPolytope() const
Integer countWithin(Point low, Point hi) const
std::vector< Integer > InequalityVector
Self & operator+=(UnitSegment s)
const Vector & getA(unsigned int i) const
BoundedRationalPolytope(const Self &other)=default
BoundedRationalPolytope(Integer d, const Domain &domain, HalfSpaceIterator itB, HalfSpaceIterator itE, bool valid_edge_constraints=false, bool check_duplicate_constraints=false)
Domain computeLatticeDomain(const Domain &d)
ClosedIntegerHalfPlane< Space > HalfSpace
bool isInterior(const Point &p) const
BoundedRationalPolytope(Integer d, PointIterator itB, PointIterator itE)
const Domain & getDomain() const
static Self base(Dimension k, Component val=1)
auto crossProduct(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::crossProduct(*this, v))
Cross product with a PointVector.
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
TInteger Integer
Arithmetic ring induced by (+,-,*) and Integer numbers.
Definition: SpaceND.h:102
static const Dimension dimension
static constants to store the dimension.
Definition: SpaceND.h:132
std::ostream & error()
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
KForm< Calculus, order, duality > operator*(const typename Calculus::Scalar &scalar, const KForm< Calculus, order, duality > &form)
mpz_class BigInteger
Multi-precision integer with GMP implementation.
Definition: BasicTypes.h:79
Circulator< TIterator > operator+(typename IteratorCirculatorTraits< TIterator >::Difference d, Circulator< TIterator > &object)
Definition: Circulator.h:453
UnitCell(std::initializer_list< Dimension > l)
friend std::ostream & operator<<(std::ostream &out, const UnitCell &object)
Aim: A half-space specified by a vector N and a constant c. The half-space is the set .
Aim: Defines the concept describing a digital space, ie a cartesian product of integer lines.
Definition: CSpace.h:106
static Vector crossProduct(const Vector &v1, const Vector &v2)
static void addEdgeConstraint(Polytope &P, unsigned int i, unsigned int j, const std::vector< Point > &pts)
Aim: It is just a helper class for BoundedRationalPolytope to add dimension specific static methods.
static void addEdgeConstraint(Polytope &, unsigned int, unsigned int, const std::vector< Point > &)
static Vector crossProduct(const Vector &, const Vector &)
int max(int a, int b)
Domain domain