Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
BoundedRationalPolytope.h
1
16
17#pragma once
18
30
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;
89
91
96 struct UnitSegment {
98 UnitSegment( Dimension d ) : k( d ) {}
99 };
100
106 struct UnitCell {
107 std::vector<Dimension> dims;
108 UnitCell( std::initializer_list<Dimension> l )
109 : dims( l.begin(), l.end() ) {}
110
117 friend std::ostream&
118 operator<< ( std::ostream & out,
119 const UnitCell & object )
120 {
121 out << "{";
122 for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ];
123 out << "}";
124 return out;
125 }
126 };
127
128
131 struct Rational {
132 Integer p; // numerator
133 Integer q; // denominator
137 inline Rational( Integer a, Integer b ) : p( a ), q( b ) {}
138 };
139
142
147
152
157 BoundedRationalPolytope ( const Self & other ) = default;
158
159
170 BoundedRationalPolytope( std::initializer_list<Point> l );
171
185 template <typename PointIterator>
186 BoundedRationalPolytope( Integer d, PointIterator itB, PointIterator itE );
187
214 template <typename HalfSpaceIterator>
216 const Domain& domain,
217 HalfSpaceIterator itB, HalfSpaceIterator itE,
218 bool valid_edge_constraints = false,
219 bool check_duplicate_constraints = false );
220
247 template <typename HalfSpaceIterator>
248 void init( Integer d, const Domain& domain,
249 HalfSpaceIterator itB, HalfSpaceIterator itE,
250 bool valid_edge_constraints = false,
251 bool check_duplicate_constraints = false );
252
270 template <typename PointIterator>
271 bool init( Integer d, PointIterator itB, PointIterator itE );
272
278 Self & operator= ( const Self & other ) = default;
279
281 void clear();
282
284
285 // ----------------------- Accessor services ------------------------------
286 public:
289
291 const Domain& getDomain() const;
292
295 const Domain& getLatticeDomain() const;
296
299 const Domain& getRationalDomain() const;
300
302 unsigned int nbHalfSpaces() const;
303
308
314 const Vector& getA( unsigned int i ) const;
315
321 Integer getB( unsigned int i ) const;
322
328 bool isLarge( unsigned int i ) const;
329
331 const InequalityMatrix& getA() const;
332
334 const InequalityVector& getB() const;
338 const std::vector<bool>& getI() const;
339
344 bool canBeSummed() const;
345
347
348 // ----------------------- Check point services ------------------------------
349 public:
350
353
358 bool isInside( const Point& p ) const;
359
366 bool isDomainPointInside( const Point& p ) const;
367
372 bool isInterior( const Point& p ) const;
373
378 bool isBoundary( const Point& p ) const;
379
381
382 // ----------------------- Modification services ------------------------------
383 public:
384
387
391
404 unsigned int cut( Dimension k, bool pos, Integer b, bool large = true );
405
423 unsigned int cut( const Vector& a, Integer b, bool large = true,
424 bool valid_edge_constraint = false );
425
442 unsigned int cut( const HalfSpace & hs, bool large = true,
443 bool valid_edge_constraint = false );
444
450
451
458
464 Self& operator*=( Rational r );
465
472 Self& operator+=( UnitSegment s );
473
480 Self& operator+=( UnitCell c );
481
483
484 // ----------------------- Enumeration services ------------------------------
485 public:
486
489
497 Integer count() const;
498
510
522
533 Integer countWithin( Point low, Point hi ) const;
534
552
561 void getPoints( std::vector<Point>& pts ) const;
562
571 void getInteriorPoints( std::vector<Point>& pts ) const;
572
581 void getBoundaryPoints( std::vector<Point>& pts ) const;
582
593 template <typename PointSet>
594 void insertPoints( PointSet& pts_set ) const;
595
597
598
599 // ----------------------- Interface --------------------------------------
600 public:
603
608 void selfDisplay ( std::ostream & out ) const;
609
616 bool isValid() const;
617
621 std::string className() const;
622
624
625 // ------------------------- Protected Datas ------------------------------
626 protected:
627 // Denominator for constraints, i.e. \f$ q A x \le B \f$.
629 // The matrix A in the polytope representation \f$ q A x \le B \f$.
631 // The vector B in the polytope representation \f$ q A x \le B \f$.
633 // Tight bounded box
635 // Lattice bounded box (i.e. floor( rationalD/q ))
637 // Are inequalities large ?
638 std::vector<bool> I;
639 // Indicates if Minkowski sums with segments will be valid
641
642 // ------------------------- Private Datas --------------------------------
643 private:
644
645
646 // ------------------------- Internals ------------------------------------
647 private:
655
662
669
677
684
685 }; // end of class BoundedRationalPolytope
686
687 namespace detail {
696 template <DGtal::Dimension N, typename TInteger>
698 typedef TInteger Integer;
700 typedef typename Space::Point Point;
701 typedef typename Space::Vector Vector;
704
717 static void
718 addEdgeConstraint( Polytope& , unsigned int , unsigned int ,
719 const std::vector<Point>& )
720 {
721 trace.error() << "[BoundedRationalPolytopeHelper::addEdgeConstraint]"
722 << " this method is only implemented in 3D." << std::endl;
723 }
724
727 static
728 Vector crossProduct( const Vector& , const Vector& )
729 {
730 trace.error() << "[BoundedRationalPolytopeHelper::crossProduct]"
731 << " this method is only implemented in 3D." << std::endl;
732 return Vector::zero;
733 }
734 };
735
743 template <typename TInteger>
745 typedef TInteger Integer;
747 typedef typename Space::Point Point;
748 typedef typename Space::Vector Vector;
751
761 static void
762 addEdgeConstraint( Polytope& P, unsigned int i, unsigned int j,
763 const std::vector<Point>& pts )
764 {
765 Vector ab = pts[ i ] - pts[ j ];
766 for ( int s = 0; s < 2; s++ )
767 for ( Dimension k = 0; k < dimension; ++k )
768 {
769 Vector n = ab.crossProduct( Point::base( k, (s == 0) ? 1 : -1 ) );
770 Integer b = n.dot( pts[ i ] );
771 std::size_t nb_in = 0;
772 for ( auto p : pts ) {
773 Integer v = n.dot( p );
774 if ( v < b ) nb_in++;
775 }
776 if ( nb_in == dimension - 1 ) {
777 P.cut( n, b, true, true );
778 }
779 }
780 }
781
785 static
786 Vector crossProduct( const Vector& v1, const Vector& v2 )
787 {
788 return v1.crossProduct( v2 );
789 }
790 };
791 }
792
795
802 template <typename TSpace>
803 std::ostream&
804 operator<< ( std::ostream & out,
805 const BoundedRationalPolytope<TSpace> & object );
806
807
813 template <typename TSpace>
817
823 template <typename TSpace>
827
828
836 template <typename TSpace>
840
848 template <typename TSpace>
852
854
855} // namespace DGtal
856
857
859// Includes inline functions.
860#include "BoundedRationalPolytope.ih"
861
862// //
864
865#endif // !defined BoundedRationalPolytope_h
866
867#undef BoundedRationalPolytope_RECURSES
868#endif // else defined(BoundedRationalPolytope_RECURSES)
Self & operator+=(UnitCell 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 >))
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)
ClosedIntegerHalfPlane< Space > HalfSpace
bool isInterior(const Point &p) const
BoundedRationalPolytope(Integer d, PointIterator itB, PointIterator itE)
Self & operator*=(Integer t)
static Self base(Dimension k, Component val=1)
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
auto crossProduct(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::crossProduct(*this, v))
Cross product with a PointVector.
PointVector< dim, Integer > Point
Definition SpaceND.h:110
static const Dimension dimension
Definition SpaceND.h:132
PointVector< dim, Integer > Vector
Definition SpaceND.h:113
detail namespace gathers internal classes and functions.
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)
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:119
Trace trace
boost::multiprecision::number< boost::multiprecision::cpp_int_backend<>, boost::multiprecision::et_off > BigInteger
Definition BasicTypes.h:75
UnitCell(std::initializer_list< Dimension > l)
friend std::ostream & operator<<(std::ostream &out, const UnitCell &object)
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