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