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