DGtal 1.4.2
Loading...
Searching...
No Matches
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
54namespace 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
92#endif
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
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...
Self & operator+=(UnitCell c)
ClosedIntegerHalfPlane< Space > HalfSpace
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 >))
Self & operator+=(UnitSegment s)
bool isLarge(unsigned int i) const
Self & operator=(const Self &other)=default
bool internalInitFromSegment3D(Point a, Point b)
void swap(BoundedRationalPolytope &other)
bool init(Integer d, PointIterator itB, PointIterator itE)
void selfDisplay(std::ostream &out) const
BoundedRationalPolytope(std::initializer_list< Point > l)
void insertPoints(PointSet &pts_set) const
bool isBoundary(const Point &p) const
const Domain & getLatticeDomain() const
const std::vector< bool > & getI() const
Integer getB(unsigned int i) const
const Domain & getRationalDomain() const
void init(Integer d, const Domain &domain, HalfSpaceIterator itB, HalfSpaceIterator itE, bool valid_edge_constraints=false, bool check_duplicate_constraints=false)
const InequalityMatrix & getA() const
bool isDomainPointInside(const Point &p) const
void getInteriorPoints(std::vector< Point > &pts) const
Self & operator*=(Rational r)
void clear()
Clears the polytope.
std::string className() const
const InequalityVector & getB() 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
bool internalInitFromSegment2D(Point a, Point b)
BoundedRationalPolytope interiorPolytope() const
Integer countWithin(Point low, Point hi) const
const Domain & getDomain() const
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)
bool isInterior(const Point &p) const
BoundedRationalPolytope< TSpace > Self
BoundedRationalPolytope(Integer d, PointIterator itB, PointIterator itE)
Self & operator*=(Integer t)
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 const Dimension dimension
static constants to store the dimension.
Definition SpaceND.h:132
TInteger Integer
Arithmetic ring induced by (+,-,*) and Integer numbers.
Definition SpaceND.h:102
std::ostream & error()
DGtal is the top-level namespace which contains all DGtal functions and types.
KForm< Calculus, order, duality > operator*(const typename Calculus::Scalar &scalar, const KForm< Calculus, order, duality > &form)
boost::int64_t int64_t
signed 94-bit integer.
Definition BasicTypes.h:74
Circulator< TIterator > operator+(typename IteratorCirculatorTraits< TIterator >::Difference d, Circulator< TIterator > &object)
Definition Circulator.h:453
std::ostream & operator<<(std::ostream &out, const ClosedIntegerHalfPlane< TSpace > &object)
DGtal::uint32_t Dimension
Definition Common.h:136
Trace trace
Definition Common.h:153
mpz_class BigInteger
Multi-precision integer with GMP implementation.
Definition BasicTypes.h:79
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