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