DGtal  1.5.beta
BoundedRationalPolytope.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file BoundedRationalPolytope.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2020/04/28
23  *
24  * Implementation of inline methods defined in BoundedRationalPolytope.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include "DGtal/math/linalg/SimpleMatrix.h"
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // ----------------------- Standard services ------------------------------
41 
42 //-----------------------------------------------------------------------------
43 template <typename TSpace>
44 DGtal::BoundedRationalPolytope<TSpace>::
45 BoundedRationalPolytope()
46  : q( NumberTraits<Integer>::ZERO ),
47  rationalD( Point::zero, Point::zero ),
48  latticeD( Point::zero, Point::zero ),
49  myValidEdgeConstraints( false )
50 {}
51 
52 //-----------------------------------------------------------------------------
53 template <typename TSpace>
54 void
55 DGtal::BoundedRationalPolytope<TSpace>::
56 clear()
57 {
58  A.clear();
59  B.clear();
60  I.clear();
61  myValidEdgeConstraints = false;
62  rationalD = Domain( Point::zero, Point::zero );
63  latticeD = Domain( Point::zero, Point::zero );
64  q = 0;
65 }
66 
67 //-----------------------------------------------------------------------------
68 template <typename TSpace>
69 DGtal::BoundedRationalPolytope<TSpace>::
70 BoundedRationalPolytope( std::initializer_list<Point> l )
71 {
72  myValidEdgeConstraints = false;
73  auto it = l.begin();
74  if ( it != l.end() )
75  {
76  Integer denom = (*it++)[ 0 ];
77  init( denom, it, l.end() );
78  }
79 }
80 
81 //-----------------------------------------------------------------------------
82 template <typename TSpace>
83 template <typename PointIterator>
84 DGtal::BoundedRationalPolytope<TSpace>::
85 BoundedRationalPolytope( Integer denom, PointIterator itB, PointIterator itE )
86 {
87  myValidEdgeConstraints = false;
88  init( denom, itB, itE );
89 }
90 
91 //-----------------------------------------------------------------------------
92 template <typename TSpace>
93 template <typename HalfSpaceIterator>
94 DGtal::BoundedRationalPolytope<TSpace>::
95 BoundedRationalPolytope( Integer denom,
96  const Domain& domain,
97  HalfSpaceIterator itB, HalfSpaceIterator itE,
98  bool valid_edge_constraints,
99  bool check_duplicate_constraints )
100  : myValidEdgeConstraints( valid_edge_constraints )
101 {
102  init( denom, domain, itB, itE,
103  valid_edge_constraints, check_duplicate_constraints );
104 }
105 
106 //-----------------------------------------------------------------------------
107 template <typename TSpace>
108 template <typename HalfSpaceIterator>
109 void
110 DGtal::BoundedRationalPolytope<TSpace>::
111 init( Integer denom,
112  const Domain& domain,
113  HalfSpaceIterator itB, HalfSpaceIterator itE,
114  bool valid_edge_constraints,
115  bool check_duplicate_constraints )
116 {
117  clear();
118  q = denom;
119  myValidEdgeConstraints = valid_edge_constraints;
120  const Dimension d = dimension;
121  const Point lo = domain.lowerBound();
122  const Point hi = domain.upperBound();
123  rationalD = Domain( lo, hi );
124  latticeD = computeLatticeDomain( rationalD );
125  // Add constraints related to sup/inf in x.
126  for ( Dimension s = 0; s < d; ++s )
127  {
128  Vector z = Vector::zero;
129  z[ s ] = q;
130  A.push_back( z );
131  B.push_back( hi[ s ] );
132  z[ s ] = -q;
133  A.push_back( z );
134  B.push_back( -lo[ s ] );
135  }
136  // Add other halfplanes
137  Integer nb_hp = 2*d;
138  if ( check_duplicate_constraints )
139  {
140  // Add other halfplanes
141  for ( auto it = itB; it != itE; ++it )
142  {
143  // Checks that is not inside.
144  const auto a = it->N;
145  const auto b = it->c;
146  const auto itAE = A.begin()+2*d;
147  const auto itF = std::find( A.begin(), itAE, a );
148  if ( itF == itAE )
149  {
150  A.push_back( q * a );
151  B.push_back( b );
152  }
153  else
154  {
155  const auto k = itF - A.begin();
156  B[ k ] = std::min( B[ k ], b );
157  }
158  }
159  }
160  else
161  { // Add other halfplanes
162  for ( auto it = itB; it != itE; ++it )
163  {
164  A.push_back( q * it->N );
165  B.push_back( it->c );
166  ++nb_hp;
167  }
168  }
169  I = std::vector<bool>( nb_hp, true ); // inequalities are large
170 }
171 
172 //-----------------------------------------------------------------------------
173 template <typename TSpace>
174 bool
175 DGtal::BoundedRationalPolytope<TSpace>::
176 internalInitFromTriangle3D( Point a, Point b, Point c )
177 {
178  Vector ab = b - a;
179  Vector bc = c - b;
180  Vector ca = a - c;
181  Vector n = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
182  crossProduct( ab, bc );
183  if ( n == Vector::zero ) { clear(); return false; }
184  A.push_back( q * n );
185  B.push_back( a.dot( n ) );
186  A.push_back( -q * n );
187  B.push_back( -a.dot( n ) );
188  Vector abn = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
189  crossProduct( ab, n );
190  A.push_back( q * abn );
191  B.push_back( a.dot( abn ) );
192  Vector bcn = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
193  crossProduct( bc, n );
194  A.push_back( q * bcn );
195  B.push_back( b.dot( bcn ) );
196  Vector can = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
197  crossProduct( ca, n );
198  A.push_back( q * can );
199  B.push_back( c.dot( can ) );
200  I = std::vector<bool>( 2*3+5, true ); // inequalities are large
201  return true;
202 }
203 
204 //-----------------------------------------------------------------------------
205 template <typename TSpace>
206 bool
207 DGtal::BoundedRationalPolytope<TSpace>::
208 internalInitFromSegment3D( Point a, Point b )
209 {
210  Vector ab = b - a;
211  if ( ab == Vector::zero ) return true; // domain and constraints already computed
212  Vector t( 1, 0, 0 );
213  Vector n = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
214  crossProduct( ab, t );
215  if ( n == Vector::zero )
216  {
217  t = Vector( 0, 1, 0 );
218  n = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
219  crossProduct( ab, t );
220  }
221  A.push_back( q * n );
222  B.push_back( a.dot( n ) );
223  A.push_back( -q * n );
224  B.push_back( -a.dot( n ) );
225  Vector w = detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::
226  crossProduct( ab, n );
227  A.push_back( q * w );
228  B.push_back( a.dot( w ) );
229  A.push_back( -q * w );
230  B.push_back( -a.dot( w ) );
231  I = std::vector<bool>( 2*3+4, true ); // inequalities are large
232  return true;
233 }
234 
235 //-----------------------------------------------------------------------------
236 template <typename TSpace>
237 bool
238 DGtal::BoundedRationalPolytope<TSpace>::
239 internalInitFromSegment2D( Point a, Point b )
240 {
241  Vector ab = b - a;
242  if ( ab == Vector::zero ) return true; // domain and constraints already computed
243  Vector n( -ab[ 1 ], ab[ 0 ] );
244  A.push_back( q * n );
245  B.push_back( a.dot( n ) );
246  A.push_back( -q * n );
247  B.push_back( -a.dot( n ) );
248  I = std::vector<bool>( 2*2+2, true ); // inequalities are large
249  return true;
250 }
251 
252 //-----------------------------------------------------------------------------
253 template <typename TSpace>
254 typename DGtal::BoundedRationalPolytope<TSpace>::Domain
255 DGtal::BoundedRationalPolytope<TSpace>::
256 computeLatticeDomain( const Domain& d )
257 {
258  Point lo, hi;
259  for ( Dimension i = 0; i < Space::dimension; i++ )
260  {
261  lo[ i ] = d.lowerBound()[ i ] / q;
262  hi[ i ] = d.upperBound()[ i ] / q;
263  }
264  return Domain( lo, hi );
265 }
266 //-----------------------------------------------------------------------------
267 template <typename TSpace>
268 typename DGtal::BoundedRationalPolytope<TSpace>::Domain
269 DGtal::BoundedRationalPolytope<TSpace>::
270 computeRationalDomain( const Domain& d )
271 {
272  Point lo, hi;
273  for ( Dimension i = 0; i < Space::dimension; i++ )
274  {
275  lo[ i ] = d.lowerBound()[ i ] * q;
276  hi[ i ] = d.upperBound()[ i ] * q;
277  }
278  return Domain( lo, hi );
279 }
280 
281 //-----------------------------------------------------------------------------
282 template <typename TSpace>
283 template <typename PointIterator>
284 bool
285 DGtal::BoundedRationalPolytope<TSpace>::
286 init( Integer denom, PointIterator itB, PointIterator itE )
287 {
288  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
289  clear();
290  q = denom;
291  const Dimension d = dimension;
292  std::vector<Point> pts;
293  for ( ; itB != itE; ++itB ) pts.push_back( *itB );
294  Point lo = pts[ 0 ];
295  Point hi = pts[ 0 ];
296  for ( Dimension s = 1; s < pts.size(); ++s )
297  {
298  lo = lo.inf( pts[ s ] );
299  hi = hi.sup( pts[ s ] );
300  }
301  // Add constraints related to sup/inf in x.
302  for ( Dimension s = 0; s < d; ++s )
303  {
304  Vector z = Vector::zero;
305  z[ s ] = q;
306  A.push_back( z );
307  B.push_back( hi[ s ] );
308  z[ s ] = -q;
309  A.push_back( z );
310  B.push_back( -lo[ s ] );
311  }
312  rationalD = Domain( lo, hi );
313  latticeD = computeLatticeDomain( rationalD );
314  if ( pts.size() != d+1 )
315  { // Some degenerated cases are taken into account.
316  myValidEdgeConstraints = true;
317  if ( d == 3 ) {
318  if ( pts.size() == 3 )
319  return internalInitFromTriangle3D( pts[ 0 ], pts[ 1 ], pts[ 2 ] );
320  else if ( pts.size() == 2 )
321  return internalInitFromSegment3D( pts[ 0 ], pts[ 1 ] );
322  } else if ( d == 2 ) {
323  if ( pts.size() == 2 )
324  return internalInitFromSegment2D( pts[ 0 ], pts[ 1 ] );
325  }
326  I = std::vector<bool>( 2*2, true ); // inequalities are large
327  if ( pts.size() == 1 ) return true;
328  clear();
329  return false;
330  }
331  // Build Matrix A and Vector b through cofactors
332  I = std::vector<bool>( 3*d+1, true ); // inequalities are large
333  Vector a;
334  Integer b;
335  for ( Dimension s = 0; s <= d; ++s )
336  {
337  // Build matrix v composed of p_i and vectors p_k - p_i for i and k != p
338  Matrix V;
339  Dimension p = (s+1) % (d+1);
340  for ( Dimension j = 0; j < d; ++j )
341  V.setComponent( 0, j, pts[ p ][ j ] - pts[ s ][ j ] );
342  for ( Dimension k = 1; k < d; ++k )
343  {
344  Dimension l = (p+k) % (d+1);
345  for ( Dimension j = 0; j < d; ++j )
346  V.setComponent( k, j, pts[ l ][ j ] - pts[ p ][ j ] );
347  }
348  b = V.determinant();
349  if ( b == 0 )
350  {
351  clear();
352  return false;
353  }
354  // Form vector [b, 0, ..., 0]
355  Vector z = Vector::zero;
356  z[ 0 ] = 1;
357  a = V.cofactor().transpose() * z;
358  b += a.dot( pts[ s ] );
359  // Check sign
360  if ( a.dot( pts[ s ] ) > b ) { a *= (Integer) -1; b *= (Integer) -1; }
361  A.push_back( q * a );
362  B.push_back( b );
363  }
364  myValidEdgeConstraints = true;
365  if ( dimension >= 3 )
366  { // One should add edges
367  for ( unsigned int i = 0; i < pts.size(); ++i )
368  for ( unsigned int j = i+1; j < pts.size(); ++j ) {
369  detail::BoundedRationalPolytopeSpecializer< dimension, Integer >::addEdgeConstraint
370  ( *this, i, j, pts );
371  }
372  if ( dimension >= 4 )
373  { // Not implemented yet
374  myValidEdgeConstraints = false;
375  }
376  }
377  return true;
378 }
379 
380 
381 //-----------------------------------------------------------------------------
382 template <typename TSpace>
383 DGtal::BoundedRationalPolytope<TSpace>
384 DGtal::BoundedRationalPolytope<TSpace>::
385 interiorPolytope() const
386 {
387  BoundedRationalPolytope P( *this );
388  P.I = std::vector<bool>( P.A.size(), false );
389  // for ( auto it = P.I.begin(), itE = P.I.end(); it != itE; ++it )
390  // *it = false;
391  return P;
392 }
393 
394 //-----------------------------------------------------------------------------
395 template <typename TSpace>
396 unsigned int
397 DGtal::BoundedRationalPolytope<TSpace>::
398 cut( Dimension k, bool pos, Integer b, bool large )
399 {
400  ASSERT( k < dimension );
401  auto i = 2*k + (pos ? 0 : 1);
402  B[ i ] = std::min( B[ i ], b );
403  I[ i ] = large;
404  Point L = rationalD.lowerBound();
405  Point U = rationalD.upperBound();
406  if ( pos ) U[ k ] = B[ i ];
407  else L[ k ] = -B[ i ];
408  rationalD = Domain( L, U );
409  latticeD = computeLatticeDomain( rationalD );
410  return k;
411 }
412 
413 //-----------------------------------------------------------------------------
414 template <typename TSpace>
415 unsigned int
416 DGtal::BoundedRationalPolytope<TSpace>::
417 cut( const Vector& a, Integer b, bool large, bool valid_edge_constraint )
418 {
419  // Checks that is not inside.
420  auto it = std::find( A.begin(), A.end(), a );
421  if ( it == A.end() )
422  {
423  A.push_back( q * a );
424  B.push_back( b );
425  I.push_back( large );
426  myValidEdgeConstraints = myValidEdgeConstraints && valid_edge_constraint; // a cut might invalidate an edge constraint
427  return static_cast<unsigned int>(A.size() - 1);
428  }
429  else
430  {
431  auto k = it - A.begin();
432  B[ k ] = std::min( B[ k ], b );
433  I[ k ] = large;
434  myValidEdgeConstraints = myValidEdgeConstraints && valid_edge_constraint; // a cut might invalidate an edge constraint
435  return static_cast<unsigned int>(k);
436  }
437 }
438 //-----------------------------------------------------------------------------
439 template <typename TSpace>
440 unsigned int
441 DGtal::BoundedRationalPolytope<TSpace>::
442 cut( const HalfSpace& hs, bool large, bool valid_edge_constraint )
443 {
444  auto a = hs.N;
445  auto b = hs.c;
446  return cut( a, b, large, valid_edge_constraint );
447 }
448 
449 //-----------------------------------------------------------------------------
450 template <typename TSpace>
451 void
452 DGtal::BoundedRationalPolytope<TSpace>::
453 swap( BoundedRationalPolytope & other )
454 {
455  A.swap( other.A );
456  B.swap( other.B );
457  I.swap( other.I );
458  std::swap( rationalD, other.rationalD );
459  std::swap( latticeD, other.latticeD );
460  std::swap( myValidEdgeConstraints, other.myValidEdgeConstraints );
461  std::swap( q, other.q );
462 }
463 
464 //-----------------------------------------------------------------------------
465 template <typename TSpace>
466 bool
467 DGtal::BoundedRationalPolytope<TSpace>::
468 isInside( const Point& p ) const
469 {
470  ASSERT( isValid() );
471  for ( Dimension i = 0; i < A.size(); ++i )
472  {
473  bool in_half_space =
474  I[ i ]
475  ? A[ i ].dot( p ) <= B[ i ]
476  : A[ i ].dot( p ) < B[ i ];
477  if ( ! in_half_space ) return false;
478  }
479  return true;
480 }
481 
482 //-----------------------------------------------------------------------------
483 template <typename TSpace>
484 bool
485 DGtal::BoundedRationalPolytope<TSpace>::
486 isDomainPointInside( const Point& p ) const
487 {
488  ASSERT( isValid() );
489  for ( Dimension i = 2*dimension; i < A.size(); ++i )
490  {
491  bool in_half_space =
492  I[ i ]
493  ? A[ i ].dot( p ) <= B[ i ]
494  : A[ i ].dot( p ) < B[ i ];
495  if ( ! in_half_space ) return false;
496  }
497  return true;
498 }
499 
500 //-----------------------------------------------------------------------------
501 template <typename TSpace>
502 bool
503 DGtal::BoundedRationalPolytope<TSpace>::
504 isInterior( const Point& p ) const
505 {
506  ASSERT( isValid() );
507  for ( Dimension i = 0; i < A.size(); ++i )
508  {
509  bool in_half_space = A[ i ].dot( p ) < B[ i ];
510  if ( ! in_half_space ) return false;
511  }
512  return true;
513 }
514 
515 //-----------------------------------------------------------------------------
516 template <typename TSpace>
517 bool
518 DGtal::BoundedRationalPolytope<TSpace>::
519 isBoundary( const Point& p ) const
520 {
521  ASSERT( isValid() );
522  bool is_boundary = false;
523  for ( Dimension i = 0; i < A.size(); ++i )
524  {
525  auto Ai_dot_p = A[ i ].dot( p );
526  if ( Ai_dot_p == B[ i ] ) is_boundary = true;
527  if ( Ai_dot_p > B[ i ] ) return false;
528  }
529  return is_boundary;
530 }
531 
532 //-----------------------------------------------------------------------------
533 template <typename TSpace>
534 typename DGtal::BoundedRationalPolytope<TSpace>::Self&
535 DGtal::BoundedRationalPolytope<TSpace>::
536 operator*=( Integer t )
537 {
538  const Integer g = IntegerComputer< Integer >::staticGcd( q, t );
539  const Integer f = t / g;
540  for ( Integer& b : B ) b *= f;
541  for ( Vector& a : A ) a /= g;
542  rationalD = Domain( rationalD.lowerBound() * f, rationalD.upperBound() * f );
543  q /= g;
544  latticeD = computeLatticeDomain( rationalD );
545  return *this;
546 }
547 
548 //-----------------------------------------------------------------------------
549 template <typename TSpace>
550 typename DGtal::BoundedRationalPolytope<TSpace>::Self&
551 DGtal::BoundedRationalPolytope<TSpace>::
552 operator*=( Rational r )
553 {
554  const Integer g = IntegerComputer< Integer >::staticGcd( q * r.q, r.p );
555  const Integer f = r.p / g;
556  for ( Integer& b : B ) { b *= f; }
557  for ( Vector& a : A ) { a *= r.q; a /= g; }
558  rationalD = Domain( rationalD.lowerBound() * f, rationalD.upperBound() * f );
559  q *= r.q;
560  q /= g;
561  latticeD = computeLatticeDomain( rationalD );
562  return *this;
563 }
564 
565 //-----------------------------------------------------------------------------
566 template <typename TSpace>
567 typename DGtal::BoundedRationalPolytope<TSpace>::Self&
568 DGtal::BoundedRationalPolytope<TSpace>::
569 operator+=( UnitSegment s )
570 {
571  for ( Dimension i = 0; i < A.size(); ++i )
572  {
573  if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
574  B[ i ] += A[ i ][ s.k ];
575  }
576  Vector z = Vector::zero;
577  z[ s.k ] = NumberTraits<Integer>::ONE;
578  rationalD = Domain( rationalD.lowerBound(), rationalD.upperBound() + q * z );
579  latticeD = Domain( latticeD.lowerBound(), latticeD.upperBound() + z );
580  return *this;
581 }
582 
583 
584 //-----------------------------------------------------------------------------
585 template <typename TSpace>
586 typename DGtal::BoundedRationalPolytope<TSpace>::Self&
587 DGtal::BoundedRationalPolytope<TSpace>::
588 operator+=( UnitCell c )
589 {
590  for ( Dimension i = 0; i < c.dims.size(); ++i )
591  *this += UnitSegment( c.dims[ i ] );
592  return *this;
593 }
594 
595 
596 //-----------------------------------------------------------------------------
597 template <typename TSpace>
598 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
599 DGtal::BoundedRationalPolytope<TSpace>::
600 count() const
601 {
602  Integer nb = 0;
603  for ( const Point & p : latticeD )
604  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
605  return nb;
606 }
607 
608 //-----------------------------------------------------------------------------
609 template <typename TSpace>
610 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
611 DGtal::BoundedRationalPolytope<TSpace>::
612 countInterior() const
613 {
614  Integer nb = 0;
615  for ( const Point & p : latticeD )
616  nb += isInterior( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
617  return nb;
618 }
619 //-----------------------------------------------------------------------------
620 template <typename TSpace>
621 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
622 DGtal::BoundedRationalPolytope<TSpace>::
623 countBoundary() const
624 {
625  Integer nb = 0;
626  for ( const Point & p : latticeD )
627  nb += isBoundary( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
628  return nb;
629 }
630 //-----------------------------------------------------------------------------
631 template <typename TSpace>
632 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
633 DGtal::BoundedRationalPolytope<TSpace>::
634 countWithin( Point lo, Point hi ) const
635 {
636  Integer nb = 0;
637  Domain D1( lo.sup( latticeD.lowerBound() ), hi.inf( latticeD.upperBound() ) );
638  for ( const Point & p : D1 )
639  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
640  return nb;
641 }
642 //-----------------------------------------------------------------------------
643 template <typename TSpace>
644 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
645 DGtal::BoundedRationalPolytope<TSpace>::
646 countUpTo( Integer max) const
647 {
648  Integer nb = 0;
649  for ( const Point & p : latticeD ) {
650  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
651  if ( nb >= max ) return max;
652  }
653  return nb;
654 }
655 //-----------------------------------------------------------------------------
656 template <typename TSpace>
657 void
658 DGtal::BoundedRationalPolytope<TSpace>::
659 getPoints( std::vector<Point>& pts ) const
660 {
661  pts.clear();
662  for ( const Point & p : latticeD )
663  if ( isDomainPointInside( p ) ) pts.push_back( p );
664 }
665 //-----------------------------------------------------------------------------
666 template <typename TSpace>
667 template <typename PointSet>
668 void
669 DGtal::BoundedRationalPolytope<TSpace>::
670 insertPoints( PointSet& pts_set ) const
671 {
672  for ( const Point & p : latticeD )
673  if ( isDomainPointInside( p ) ) pts_set.insert( p );
674 }
675 //-----------------------------------------------------------------------------
676 template <typename TSpace>
677 void
678 DGtal::BoundedRationalPolytope<TSpace>::
679 getInteriorPoints( std::vector<Point>& pts ) const
680 {
681  pts.clear();
682  for ( const Point & p : latticeD )
683  if ( isInterior( p ) ) pts.push_back( p );
684 }
685 //-----------------------------------------------------------------------------
686 template <typename TSpace>
687 void
688 DGtal::BoundedRationalPolytope<TSpace>::
689 getBoundaryPoints( std::vector<Point>& pts ) const
690 {
691  pts.clear();
692  for ( const Point & p : latticeD )
693  if ( isBoundary( p ) ) pts.push_back( p );
694 }
695 
696 //-----------------------------------------------------------------------------
697 template <typename TSpace>
698 const typename DGtal::BoundedRationalPolytope<TSpace>::Domain&
699 DGtal::BoundedRationalPolytope<TSpace>::getDomain() const
700 {
701  return latticeD;
702 }
703 
704 //-----------------------------------------------------------------------------
705 template <typename TSpace>
706 const typename DGtal::BoundedRationalPolytope<TSpace>::Domain&
707 DGtal::BoundedRationalPolytope<TSpace>::getLatticeDomain() const
708 {
709  return latticeD;
710 }
711 
712 //-----------------------------------------------------------------------------
713 template <typename TSpace>
714 const typename DGtal::BoundedRationalPolytope<TSpace>::Domain&
715 DGtal::BoundedRationalPolytope<TSpace>::getRationalDomain() const
716 {
717  return rationalD;
718 }
719 
720 //-----------------------------------------------------------------------------
721 template <typename TSpace>
722 unsigned int
723 DGtal::BoundedRationalPolytope<TSpace>::nbHalfSpaces() const
724 {
725  return A.size();
726 }
727 
728 //-----------------------------------------------------------------------------
729 template <typename TSpace>
730 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
731 DGtal::BoundedRationalPolytope<TSpace>::denominator() const
732 {
733  return q;
734 }
735 
736 
737 //-----------------------------------------------------------------------------
738 template <typename TSpace>
739 const typename DGtal::BoundedRationalPolytope<TSpace>::Vector&
740 DGtal::BoundedRationalPolytope<TSpace>::getA( unsigned int i ) const
741 {
742  ASSERT( i < nbHalfSpaces() );
743  return A[ i ];
744 }
745 
746 //-----------------------------------------------------------------------------
747 template <typename TSpace>
748 typename DGtal::BoundedRationalPolytope<TSpace>::Integer
749 DGtal::BoundedRationalPolytope<TSpace>::getB( unsigned int i ) const
750 {
751  ASSERT( i < nbHalfSpaces() );
752  return B[ i ];
753 }
754 
755 //-----------------------------------------------------------------------------
756 template <typename TSpace>
757 bool
758 DGtal::BoundedRationalPolytope<TSpace>::isLarge( unsigned int i ) const
759 {
760  ASSERT( i < nbHalfSpaces() );
761  return I[ i ];
762 }
763 
764 //-----------------------------------------------------------------------------
765 template <typename TSpace>
766 const typename DGtal::BoundedRationalPolytope<TSpace>::InequalityMatrix&
767 DGtal::BoundedRationalPolytope<TSpace>::getA() const
768 {
769  return A;
770 }
771 
772 //-----------------------------------------------------------------------------
773 template <typename TSpace>
774 const typename DGtal::BoundedRationalPolytope<TSpace>::InequalityVector&
775 DGtal::BoundedRationalPolytope<TSpace>::getB() const
776 {
777  return B;
778 }
779 
780 //-----------------------------------------------------------------------------
781 template <typename TSpace>
782 const std::vector<bool>&
783 DGtal::BoundedRationalPolytope<TSpace>::getI() const
784 {
785  return I;
786 }
787 
788 //-----------------------------------------------------------------------------
789 template <typename TSpace>
790 bool
791 DGtal::BoundedRationalPolytope<TSpace>::canBeSummed() const
792 {
793  return myValidEdgeConstraints;
794 }
795 
796 ///////////////////////////////////////////////////////////////////////////////
797 // Interface - public :
798 
799 /**
800  * Writes/Displays the object on an output stream.
801  * @param out the output stream where the object is written.
802  */
803 template <typename TSpace>
804 inline
805 void
806 DGtal::BoundedRationalPolytope<TSpace>::selfDisplay ( std::ostream & out ) const
807 {
808  out << "[BoundedRationalPolytope<" << Space::dimension << "> A.rows=" << A.size()
809  << " valid_edge_constraints=" << myValidEdgeConstraints
810  << " denom=" << q << "]" << std::endl;
811  for ( Dimension i = 0; i < A.size(); ++i )
812  {
813  out << " [";
814  for ( Dimension j = 0; j < dimension; ++j )
815  out << " " << A[ i ][ j ];
816  out << " ] . x <= " << B[ i ] << std::endl;
817  }
818 }
819 
820 /**
821  * Checks the validity/consistency of the object.
822  * @return 'true' if the object is valid, 'false' otherwise.
823  */
824 template <typename TSpace>
825 inline
826 bool
827 DGtal::BoundedRationalPolytope<TSpace>::isValid() const
828 {
829  return q > 0 && ! rationalD.isEmpty();
830 }
831 //-----------------------------------------------------------------------------
832 template <typename TSpace>
833 inline
834 std::string
835 DGtal::BoundedRationalPolytope<TSpace>::className
836 () const
837 {
838  return "BoundedRationalPolytope";
839 }
840 
841 
842 
843 ///////////////////////////////////////////////////////////////////////////////
844 // Implementation of inline functions //
845 
846 //-----------------------------------------------------------------------------
847 template <typename TSpace>
848 inline
849 std::ostream&
850 DGtal::operator<< ( std::ostream & out,
851  const BoundedRationalPolytope<TSpace> & object )
852 {
853  object.selfDisplay( out );
854  return out;
855 }
856 //-----------------------------------------------------------------------------
857 template <typename TSpace>
858 DGtal::BoundedRationalPolytope<TSpace>
859 DGtal::operator* ( typename BoundedRationalPolytope<TSpace>::Integer t,
860  const BoundedRationalPolytope<TSpace> & P )
861 {
862  BoundedRationalPolytope<TSpace> Q = P;
863  Q *= t;
864  return Q;
865 }
866 //-----------------------------------------------------------------------------
867 template <typename TSpace>
868 DGtal::BoundedRationalPolytope<TSpace>
869 DGtal::operator* ( typename BoundedRationalPolytope<TSpace>::Rational r,
870  const BoundedRationalPolytope<TSpace> & P )
871 {
872  BoundedRationalPolytope<TSpace> Q = P;
873  Q *= r;
874  return Q;
875 }
876 //-----------------------------------------------------------------------------
877 template <typename TSpace>
878 DGtal::BoundedRationalPolytope<TSpace>
879 DGtal::operator+ ( const BoundedRationalPolytope<TSpace> & P,
880  typename BoundedRationalPolytope<TSpace>::UnitSegment s )
881 {
882  BoundedRationalPolytope<TSpace> Q = P;
883  Q += s;
884  return Q;
885 }
886 //-----------------------------------------------------------------------------
887 template <typename TSpace>
888 DGtal::BoundedRationalPolytope<TSpace>
889 DGtal::operator+ ( const BoundedRationalPolytope<TSpace> & P,
890  typename BoundedRationalPolytope<TSpace>::UnitCell c )
891 {
892  BoundedRationalPolytope<TSpace> Q = P;
893  Q += c;
894  return Q;
895 }
896 
897 // //
898 ///////////////////////////////////////////////////////////////////////////////