DGtal  1.5.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 #include "DGtal/geometry/volumes/BoundedLatticePolytopeCounter.h"
34 //////////////////////////////////////////////////////////////////////////////
35 
36 ///////////////////////////////////////////////////////////////////////////////
37 // IMPLEMENTATION of inline methods.
38 ///////////////////////////////////////////////////////////////////////////////
39 
40 ///////////////////////////////////////////////////////////////////////////////
41 // ----------------------- Standard services ------------------------------
42 
43 //-----------------------------------------------------------------------------
44 template <typename TSpace>
45 void
46 DGtal::BoundedLatticePolytope<TSpace>::
47 clear()
48 {
49  A.clear();
50  B.clear();
51  I.clear();
52  myValidEdgeConstraints = false;
53  D = Domain( Point::zero, Point::zero );
54 }
55 
56 //-----------------------------------------------------------------------------
57 template <typename TSpace>
58 DGtal::BoundedLatticePolytope<TSpace>::
59 BoundedLatticePolytope( std::initializer_list<Point> l )
60 {
61  myValidEdgeConstraints = false;
62  init( l.begin(), l.end() );
63 }
64 
65 //-----------------------------------------------------------------------------
66 template <typename TSpace>
67 template <typename PointIterator>
68 DGtal::BoundedLatticePolytope<TSpace>::
69 BoundedLatticePolytope( PointIterator itB, PointIterator itE )
70 {
71  myValidEdgeConstraints = false;
72  init( itB, itE );
73 }
74 
75 //-----------------------------------------------------------------------------
76 template <typename TSpace>
77 template <typename HalfSpaceIterator>
78 DGtal::BoundedLatticePolytope<TSpace>::
79 BoundedLatticePolytope( const Domain& domain,
80  HalfSpaceIterator itB, HalfSpaceIterator itE,
81  bool valid_edge_constraints,
82  bool check_duplicate_constraints )
83  : myValidEdgeConstraints( valid_edge_constraints )
84 {
85  init( domain, itB, itE,
86  valid_edge_constraints, check_duplicate_constraints );
87 }
88 
89 //-----------------------------------------------------------------------------
90 template <typename TSpace>
91 template <typename HalfSpaceIterator>
92 void
93 DGtal::BoundedLatticePolytope<TSpace>::
94 init( const Domain& domain,
95  HalfSpaceIterator itB, HalfSpaceIterator itE,
96  bool valid_edge_constraints,
97  bool check_duplicate_constraints )
98 {
99  clear();
100  myValidEdgeConstraints = valid_edge_constraints;
101  const Dimension d = dimension;
102  const Point lo = domain.lowerBound();
103  const Point hi = domain.upperBound();
104  D = Domain( lo, hi );
105  // Add constraints related to sup/inf in x.
106  for ( Dimension s = 0; s < d; ++s )
107  {
108  Vector z = Vector::zero;
109  z[ s ] = NumberTraits<Integer>::ONE;
110  A.push_back( z );
111  B.push_back( hi[ s ] );
112  z[ s ] = -NumberTraits<Integer>::ONE;
113  A.push_back( z );
114  B.push_back( -lo[ s ] );
115  }
116  Integer nb_hp = 2*d;
117  if ( check_duplicate_constraints )
118  {
119  // Add other halfplanes
120  for ( auto it = itB; it != itE; ++it )
121  {
122  // Checks that is not inside.
123  const auto a = it->N;
124  const auto b = it->c;
125  const auto itAE = A.begin()+2*d;
126  const auto itF = std::find( A.begin(), itAE , a );
127  if ( itF == itAE )
128  {
129  A.push_back( a );
130  B.push_back( b );
131  ++nb_hp;
132  }
133  else
134  {
135  const auto k = itF - A.begin();
136  B[ k ] = std::min( B[ k ], b );
137  }
138  }
139  }
140  else
141  { // Add other halfplanes
142  for ( auto it = itB; it != itE; ++it )
143  {
144  A.push_back( it->N );
145  B.push_back( it->c );
146  ++nb_hp;
147  }
148  }
149  I = std::vector<bool>( nb_hp, true ); // inequalities are large
150 }
151 
152 //-----------------------------------------------------------------------------
153 template <typename TSpace>
154 bool
155 DGtal::BoundedLatticePolytope<TSpace>::
156 internalInitFromTriangle3D( Point a, Point b, Point c )
157 {
158  Vector ab = b - a;
159  Vector bc = c - b;
160  Vector ca = a - c;
161  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::
162  crossProduct( ab, bc );
163  if ( n == Vector::zero ) { clear(); return false; }
164  A.push_back( n );
165  B.push_back( a.dot( A.back() ) );
166  Vector mn = -1 * n;
167  A.push_back( mn );
168  B.push_back( a.dot( A.back() ) );
169  I = std::vector<bool>( B.size(), true ); // inequalities are large
170  std::vector<Point> pts = { a, b, c };
171  detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::
172  addEdgeConstraint( *this, 0, 1, pts );
173  detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::
174  addEdgeConstraint( *this, 1, 2, pts );
175  detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::
176  addEdgeConstraint( *this, 2, 0, pts );
177  return true;
178 }
179 
180 //-----------------------------------------------------------------------------
181 template <typename TSpace>
182 bool
183 DGtal::BoundedLatticePolytope<TSpace>::
184 internalInitFromSegment3D( Point a, Point b )
185 {
186  Vector ab = b - a;
187  if ( ab == Vector::zero ) return true; // domain and constraints already computed
188  Dimension nb = 0;
189  for ( Dimension k = 0; k < 3; ++k )
190  {
191  const auto t = Vector::base( k );
192  Vector w = detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::
193  crossProduct( ab, t );
194  if ( w == Vector::zero ) continue;
195  A.push_back( w );
196  B.push_back( a.dot( w ) );
197  A.push_back( -1 * w );
198  B.push_back( a.dot( A.back() ) );
199  nb += 2;
200  }
201  I = std::vector<bool>( 2 * 3 + nb, true ); // inequalities are large
202  return true;
203 }
204 
205 //-----------------------------------------------------------------------------
206 template <typename TSpace>
207 bool
208 DGtal::BoundedLatticePolytope<TSpace>::
209 internalInitFromSegment2D( Point a, Point b )
210 {
211  Vector ab = b - a;
212  if ( ab == Vector::zero ) return true; // domain and constraints already computed
213  Vector n( -ab[ 1 ], ab[ 0 ] );
214  A.push_back( n );
215  B.push_back( a.dot( n ) );
216  A.push_back( -1 * n );
217  B.push_back( a.dot( A.back() ) );
218  I = std::vector<bool>( 2*2+2, true ); // inequalities are large
219  return true;
220 }
221 
222 //-----------------------------------------------------------------------------
223 template <typename TSpace>
224 template <typename PointIterator>
225 bool
226 DGtal::BoundedLatticePolytope<TSpace>::
227 init( PointIterator itB, PointIterator itE )
228 {
229  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
230  clear();
231  const Dimension d = dimension;
232  std::vector<Point> pts;
233  for ( ; itB != itE; ++itB ) pts.push_back( *itB );
234  Point lo = pts[ 0 ];
235  Point hi = pts[ 0 ];
236  for ( Dimension s = 1; s < pts.size(); ++s )
237  {
238  lo = lo.inf( pts[ s ] );
239  hi = hi.sup( pts[ s ] );
240  }
241  // Add constraints related to sup/inf in x.
242  for ( Dimension s = 0; s < d; ++s )
243  {
244  Vector z = Vector::zero;
245  z[ s ] = NumberTraits<Integer>::ONE;
246  A.push_back( z );
247  B.push_back( hi[ s ] );
248  z[ s ] = -NumberTraits<Integer>::ONE;
249  A.push_back( z );
250  B.push_back( -lo[ s ] );
251  }
252  D = Domain( lo, hi );
253  if ( pts.size() != d+1 )
254  { // Some degenerated cases are taken into account.
255  myValidEdgeConstraints = true;
256  if ( d == 3 ) {
257  if ( pts.size() == 3 )
258  return internalInitFromTriangle3D( pts[ 0 ], pts[ 1 ], pts[ 2 ] );
259  else if ( pts.size() == 2 )
260  return internalInitFromSegment3D( pts[ 0 ], pts[ 1 ] );
261  } else if ( d == 2 ) {
262  if ( pts.size() == 2 )
263  return internalInitFromSegment2D( pts[ 0 ], pts[ 1 ] );
264  }
265  I = std::vector<bool>( 2*2, true ); // inequalities are large
266  if ( pts.size() == 1 ) return true;
267  clear();
268  return false;
269  }
270  // Build Matrix A and Vector b through cofactors
271  I = std::vector<bool>( 3*d+1, true ); // inequalities are large
272  Vector a;
273  Integer b;
274  for ( Dimension s = 0; s <= d; ++s )
275  {
276  // Build matrix v composed of p_i and vectors p_k - p_i for i and k != p
277  Matrix V;
278  Dimension p = (s+1) % (d+1);
279  for ( Dimension j = 0; j < d; ++j )
280  V.setComponent( 0, j, pts[ p ][ j ] - pts[ s ][ j ] );
281  for ( Dimension k = 1; k < d; ++k )
282  {
283  Dimension l = (p+k) % (d+1);
284  for ( Dimension j = 0; j < d; ++j )
285  V.setComponent( k, j, pts[ l ][ j ] - pts[ p ][ j ] );
286  }
287  b = V.determinant();
288  if ( b == 0 )
289  {
290  clear();
291  D = Domain();
292  return false;
293  }
294  // Form vector [b, 0, ..., 0]
295  Vector z = Vector::zero;
296  z[ 0 ] = 1;
297  a = V.cofactor().transpose() * z;
298  b += a.dot( pts[ s ] );
299  // Check sign
300  if ( a.dot( pts[ s ] ) > b ) { a *= (Integer) -1; b *= (Integer) -1; }
301  A.push_back( a );
302  B.push_back( b );
303  }
304  myValidEdgeConstraints = ( dimension == 2 );
305  if ( dimension == 3 )
306  { // One should add edges
307  for ( unsigned int i = 0; i < pts.size(); ++i )
308  for ( unsigned int j = i+1; j < pts.size(); ++j ) {
309  detail::BoundedLatticePolytopeSpecializer< dimension, Integer >::addEdgeConstraint
310  ( *this, i, j, pts );
311  }
312  myValidEdgeConstraints = true;
313  }
314  // not implemented for dimension > 3
315  return true;
316 }
317 
318 
319 //-----------------------------------------------------------------------------
320 template <typename TSpace>
321 DGtal::BoundedLatticePolytope<TSpace>
322 DGtal::BoundedLatticePolytope<TSpace>::
323 interiorPolytope() const
324 {
325  BoundedLatticePolytope P( *this );
326  for ( auto it = P.I.begin(), itE = P.I.end(); it != itE; ++it )
327  *it = false;
328  return P;
329 }
330 
331 //-----------------------------------------------------------------------------
332 template <typename TSpace>
333 DGtal::BoundedLatticePolytope<TSpace>
334 DGtal::BoundedLatticePolytope<TSpace>::
335 closurePolytope() const
336 {
337  BoundedLatticePolytope P( *this );
338  for ( auto it = P.I.begin(), itE = P.I.end(); it != itE; ++it )
339  *it = true;
340  return P;
341 }
342 
343 //-----------------------------------------------------------------------------
344 template <typename TSpace>
345 unsigned int
346 DGtal::BoundedLatticePolytope<TSpace>::
347 cut( Dimension k, bool pos, Integer b, bool large )
348 {
349  ASSERT( k < dimension );
350  auto i = 2*k + (pos ? 0 : 1);
351  B[ i ] = std::min( B[ i ], b );
352  I[ i ] = large;
353  Point L = D.lowerBound();
354  Point U = D.upperBound();
355  if ( pos ) U[ k ] = B[ i ];
356  else L[ k ] = -B[ i ];
357  D = Domain( L, U );
358  return k;
359 }
360 
361 //-----------------------------------------------------------------------------
362 template <typename TSpace>
363 unsigned int
364 DGtal::BoundedLatticePolytope<TSpace>::
365 cut( const Vector& a, Integer b, bool large, bool valid_edge_constraint )
366 {
367  // Checks that is not inside.
368  auto it = std::find( A.begin(), A.end(), a );
369  if ( it == A.end() )
370  {
371  A.push_back( a );
372  B.push_back( b );
373  I.push_back( large );
374  myValidEdgeConstraints = myValidEdgeConstraints && valid_edge_constraint; // a cut might invalidate an edge constraint
375  return (unsigned int)(A.size() - 1);
376  }
377  else
378  {
379  auto k = it - A.begin();
380  B[ k ] = std::min( B[ k ], b );
381  I[ k ] = large;
382  myValidEdgeConstraints = myValidEdgeConstraints && valid_edge_constraint; // a cut might invalidate an edge constraint
383  return (unsigned int)k;
384  }
385 }
386 //-----------------------------------------------------------------------------
387 template <typename TSpace>
388 unsigned int
389 DGtal::BoundedLatticePolytope<TSpace>::
390 cut( const HalfSpace& hs, bool large, bool valid_edge_constraint )
391 {
392  auto a = hs.N;
393  auto b = hs.c;
394  return cut( a, b, large, valid_edge_constraint );
395 }
396 
397 //-----------------------------------------------------------------------------
398 template <typename TSpace>
399 void
400 DGtal::BoundedLatticePolytope<TSpace>::
401 swap( BoundedLatticePolytope & other )
402 {
403  A.swap( other.A );
404  B.swap( other.B );
405  I.swap( other.I );
406  std::swap( D, other.D );
407  std::swap( myValidEdgeConstraints, other.myValidEdgeConstraints );
408 }
409 
410 //-----------------------------------------------------------------------------
411 template <typename TSpace>
412 bool
413 DGtal::BoundedLatticePolytope<TSpace>::
414 isInside( const Point& p ) const
415 {
416  ASSERT( isValid() );
417  for ( Dimension i = 0; i < A.size(); ++i )
418  {
419  bool in_half_space =
420  I[ i ]
421  ? A[ i ].dot( p ) <= B[ i ]
422  : A[ i ].dot( p ) < B[ i ];
423  if ( ! in_half_space ) return false;
424  }
425  return true;
426 }
427 
428 //-----------------------------------------------------------------------------
429 template <typename TSpace>
430 bool
431 DGtal::BoundedLatticePolytope<TSpace>::
432 isDomainPointInside( const Point& p ) const
433 {
434  ASSERT( isValid() );
435  for ( Dimension i = 2*dimension; i < A.size(); ++i )
436  {
437  bool in_half_space =
438  I[ i ]
439  ? A[ i ].dot( p ) <= B[ i ]
440  : A[ i ].dot( p ) < B[ i ];
441  if ( ! in_half_space ) return false;
442  }
443  return true;
444 }
445 
446 //-----------------------------------------------------------------------------
447 template <typename TSpace>
448 bool
449 DGtal::BoundedLatticePolytope<TSpace>::
450 isInterior( const Point& p ) const
451 {
452  ASSERT( isValid() );
453  for ( Dimension i = 0; i < A.size(); ++i )
454  {
455  bool in_half_space = A[ i ].dot( p ) < B[ i ];
456  if ( ! in_half_space ) return false;
457  }
458  return true;
459 }
460 
461 //-----------------------------------------------------------------------------
462 template <typename TSpace>
463 bool
464 DGtal::BoundedLatticePolytope<TSpace>::
465 isBoundary( const Point& p ) const
466 {
467  ASSERT( isValid() );
468  bool is_boundary = false;
469  for ( Dimension i = 0; i < A.size(); ++i )
470  {
471  auto Ai_dot_p = A[ i ].dot( p );
472  if ( Ai_dot_p == B[ i ] ) is_boundary = true;
473  if ( Ai_dot_p > B[ i ] ) return false;
474  }
475  return is_boundary;
476 }
477 
478 //-----------------------------------------------------------------------------
479 template <typename TSpace>
480 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
481 DGtal::BoundedLatticePolytope<TSpace>::
482 operator*=( Integer t )
483 {
484  for ( Integer& b : B ) b *= t;
485  D = Domain( D.lowerBound() * t, D.upperBound() * t );
486  return *this;
487 }
488 
489 //-----------------------------------------------------------------------------
490 template <typename TSpace>
491 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
492 DGtal::BoundedLatticePolytope<TSpace>::
493 operator+=( UnitSegment s )
494 {
495  for ( Dimension i = 0; i < A.size(); ++i )
496  {
497  if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
498  B[ i ] += A[ i ][ s.k ];
499  }
500  Vector z = Vector::zero;
501  z[ s.k ] = NumberTraits<Integer>::ONE;
502  D = Domain( D.lowerBound(), D.upperBound() + z );
503  return *this;
504 }
505 
506 //-----------------------------------------------------------------------------
507 template <typename TSpace>
508 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
509 DGtal::BoundedLatticePolytope<TSpace>::
510 operator+=( LeftStrictUnitSegment s )
511 {
512  I[ 2*s.k + 1 ] = false;
513  for ( Dimension i = 0; i < A.size(); ++i )
514  {
515  if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
516  B[ i ] += A[ i ][ s.k ];
517  if ( A[ i ][ s.k ] < NumberTraits<Integer>::ZERO )
518  I[ i ] = false;
519  }
520  Vector z = Vector::zero;
521  z[ s.k ] = NumberTraits<Integer>::ONE;
522  D = Domain( D.lowerBound() + z, D.upperBound() + z );
523  return *this;
524 }
525 
526 //-----------------------------------------------------------------------------
527 template <typename TSpace>
528 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
529 DGtal::BoundedLatticePolytope<TSpace>::
530 operator+=( RightStrictUnitSegment s )
531 {
532  I[ 2*s.k ] = false;
533  for ( Dimension i = 0; i < A.size(); ++i )
534  {
535  if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO ) {
536  B[ i ] += A[ i ][ s.k ];
537  I[ i ] = false;
538  }
539  }
540  Vector z = Vector::zero;
541  z[ s.k ] = NumberTraits<Integer>::ONE;
542  return *this;
543 }
544 
545 //-----------------------------------------------------------------------------
546 template <typename TSpace>
547 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
548 DGtal::BoundedLatticePolytope<TSpace>::
549 operator+=( UnitCell c )
550 {
551  for ( Dimension i = 0; i < c.dims.size(); ++i )
552  *this += UnitSegment( 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+=( RightStrictUnitCell c )
561 {
562  for ( Dimension i = 0; i < c.dims.size(); ++i )
563  *this += RightStrictUnitSegment( c.dims[ i ] );
564  return *this;
565 }
566 
567 //-----------------------------------------------------------------------------
568 template <typename TSpace>
569 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
570 DGtal::BoundedLatticePolytope<TSpace>::
571 operator+=( LeftStrictUnitCell c )
572 {
573  for ( Dimension i = 0; i < c.dims.size(); ++i )
574  *this += LeftStrictUnitSegment( c.dims[ i ] );
575  return *this;
576 }
577 
578 //-----------------------------------------------------------------------------
579 template <typename TSpace>
580 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
581 DGtal::BoundedLatticePolytope<TSpace>::
582 count() const
583 {
584  BoundedLatticePolytopeCounter<Space> C( *this );
585  return C.countAlongAxis( C.longestAxis() );
586 }
587 
588 //-----------------------------------------------------------------------------
589 template <typename TSpace>
590 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
591 DGtal::BoundedLatticePolytope<TSpace>::
592 countInterior() const
593 {
594  BoundedLatticePolytopeCounter<Space> C( *this );
595  return C.countInteriorAlongAxis( C.longestAxis() );
596 }
597 //-----------------------------------------------------------------------------
598 template <typename TSpace>
599 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
600 DGtal::BoundedLatticePolytope<TSpace>::
601 countBoundary() const
602 {
603  const auto clP = closurePolytope();
604  return clP.count() - countInterior();
605 }
606 //-----------------------------------------------------------------------------
607 template <typename TSpace>
608 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
609 DGtal::BoundedLatticePolytope<TSpace>::
610 countWithin( Point lo, Point hi ) const
611 {
612  BoundedLatticePolytopeCounter<Space> C( *this );
613  Dimension b = 0;
614  lo = lo.sup( D.lowerBound() );
615  hi = hi.inf( D.upperBound() );
616  auto b_size = hi[ 0 ] - lo[ 0 ];
617  for ( Dimension a = 1; a < dimension; a++ )
618  {
619  const auto a_size = hi[ a ] - lo[ a ];
620  if ( b_size < a_size ) { b = a; b_size = a_size; }
621  }
622  hi[ b ] = lo[ b ];
623  Integer nb = 0;
624  Domain localD( lo, hi );
625  for ( auto&& p : localD )
626  {
627  auto II = C.intersectionIntervalAlongAxis( p, b );
628  nb += II.second - II.first;
629  }
630  return nb;
631 }
632 //-----------------------------------------------------------------------------
633 template <typename TSpace>
634 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
635 DGtal::BoundedLatticePolytope<TSpace>::
636 countUpTo( Integer max) const
637 {
638  BoundedLatticePolytopeCounter<Space> C( *this );
639  Dimension a = C.longestAxis();
640  Integer nb = 0;
641  Point lo = D.lowerBound();
642  Point hi = D.upperBound();
643  hi[ a ] = lo[ a ];
644  Domain localD( lo, hi );
645  for ( auto&& p : localD )
646  {
647  auto II = C.intersectionIntervalAlongAxis( p, a );
648  nb += II.second - II.first;
649  if ( nb >= max ) return max;
650  }
651  return nb;
652 }
653 //-----------------------------------------------------------------------------
654 template <typename TSpace>
655 void
656 DGtal::BoundedLatticePolytope<TSpace>::
657 getPoints( std::vector<Point>& pts ) const
658 {
659  pts.clear();
660  BoundedLatticePolytopeCounter<Space> C( *this );
661  C.getPointsAlongAxis( pts, C.longestAxis() );
662 }
663 //-----------------------------------------------------------------------------
664 template <typename TSpace>
665 void
666 DGtal::BoundedLatticePolytope<TSpace>::
667 getKPoints( std::vector<Point>& pts, const Point& alpha_shift ) const
668 {
669  pts.clear();
670  BoundedLatticePolytopeCounter<Space> C( *this );
671  Dimension a = C.longestAxis();
672  Integer nb = 0;
673  Point lo = D.lowerBound();
674  Point hi = D.upperBound();
675  hi[ a ] = lo[ a ];
676  Domain localD( lo, hi );
677  for ( auto&& p : localD )
678  {
679  auto II = C.intersectionIntervalAlongAxis( p, a );
680  Point q( 2*p - alpha_shift );
681  q[ a ] = 2*II.first - alpha_shift[ a ];
682  for ( Integer x = II.first; x < II.second; x++ )
683  {
684  pts.push_back( q );
685  q[ a ] += 2;
686  }
687  }
688 }
689 //-----------------------------------------------------------------------------
690 template <typename TSpace>
691 template <typename PointSet>
692 void
693 DGtal::BoundedLatticePolytope<TSpace>::
694 insertPoints( PointSet& pts_set ) const
695 {
696  std::vector<Point> pts;
697  getPoints( pts );
698  pts_set.insert( pts.cbegin(), pts.cend() );
699 }
700 //-----------------------------------------------------------------------------
701 template <typename TSpace>
702 template <typename PointSet>
703 void
704 DGtal::BoundedLatticePolytope<TSpace>::
705 insertKPoints( PointSet& pts_set, const Point& alpha_shift ) const
706 {
707  BoundedLatticePolytopeCounter<Space> C( *this );
708  Dimension a = C.longestAxis();
709  Integer nb = 0;
710  Point lo = D.lowerBound();
711  Point hi = D.upperBound();
712  hi[ a ] = lo[ a ];
713  Domain localD( lo, hi );
714  for ( auto&& p : localD )
715  {
716  auto II = C.intersectionIntervalAlongAxis( p, a );
717  Point q( 2*p - alpha_shift );
718  q[ a ] = 2*II.first - alpha_shift[ a ];
719  for ( Integer x = II.first; x < II.second; x++ )
720  {
721  pts_set.insert( q );
722  q[ a ] += 2;
723  }
724  }
725 }
726 //-----------------------------------------------------------------------------
727 template <typename TSpace>
728 void
729 DGtal::BoundedLatticePolytope<TSpace>::
730 getInteriorPoints( std::vector<Point>& pts ) const
731 {
732  pts.clear();
733  BoundedLatticePolytopeCounter<Space> C( *this );
734  C.getInteriorPointsAlongAxis( pts, C.longestAxis() );
735 }
736 //-----------------------------------------------------------------------------
737 template <typename TSpace>
738 void
739 DGtal::BoundedLatticePolytope<TSpace>::
740 getBoundaryPoints( std::vector<Point>& pts ) const
741 {
742  const auto clP = closurePolytope();
743  BoundedLatticePolytopeCounter<Space> C ( *this );
744  BoundedLatticePolytopeCounter<Space> clC( clP );
745  pts.clear();
746  const Dimension a = clC.longestAxis();
747  Point lo = clP.getDomain().lowerBound();
748  Point hi = clP.getDomain().upperBound();
749  hi[ a ] = lo[ a ];
750  Domain localD( lo, hi );
751  for ( auto&& p : localD )
752  {
753  auto II = C .interiorIntersectionIntervalAlongAxis( p, a );
754  auto clI = clC.intersectionIntervalAlongAxis( p, a );
755  auto nbI = II.second - II.first;
756  auto nbclI = clI.second - clI.first;
757  if ( nbI == nbclI ) continue;
758  if ( nbI > nbclI )
759  trace.error() << "BoundedLatticePolytope::getBoundaryPoints: bad count"
760  << std::endl;
761  Point q = p;
762  if ( nbI == 0 )
763  {
764  for ( Integer x = clI.first; x != clI.second; x++ )
765  {
766  q[ a ] = x;
767  pts.push_back( q );
768  }
769  }
770  else
771  {
772  if ( clI.first < II.first )
773  {
774  q[ a ] = clI.first;
775  pts.push_back( q );
776  }
777  if ( clI.second > II.second )
778  {
779  q[ a ] = II.second;
780  pts.push_back( q );
781  }
782  }
783  }
784 }
785 
786 
787 //-----------------------------------------------------------------------------
788 template <typename TSpace>
789 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
790 DGtal::BoundedLatticePolytope<TSpace>::
791 countByScanning() const
792 {
793  Integer nb = 0;
794  for ( const Point & p : D )
795  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
796  return nb;
797 }
798 
799 //-----------------------------------------------------------------------------
800 template <typename TSpace>
801 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
802 DGtal::BoundedLatticePolytope<TSpace>::
803 countInteriorByScanning() const
804 {
805  Integer nb = 0;
806  for ( const Point & p : D )
807  nb += isInterior( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
808  return nb;
809 }
810 //-----------------------------------------------------------------------------
811 template <typename TSpace>
812 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
813 DGtal::BoundedLatticePolytope<TSpace>::
814 countBoundaryByScanning() const
815 {
816  Integer nb = 0;
817  for ( const Point & p : D )
818  nb += isBoundary( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
819  return nb;
820 }
821 //-----------------------------------------------------------------------------
822 template <typename TSpace>
823 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
824 DGtal::BoundedLatticePolytope<TSpace>::
825 countWithinByScanning( Point lo, Point hi ) const
826 {
827  Integer nb = 0;
828  Domain D1( lo.sup( D.lowerBound() ), hi.inf( D.upperBound() ) );
829  for ( const Point & p : D1 )
830  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
831  return nb;
832 }
833 //-----------------------------------------------------------------------------
834 template <typename TSpace>
835 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
836 DGtal::BoundedLatticePolytope<TSpace>::
837 countUpToByScanning( Integer max) const
838 {
839  Integer nb = 0;
840  for ( const Point & p : D ) {
841  nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
842  if ( nb >= max ) return max;
843  }
844  return nb;
845 }
846 //-----------------------------------------------------------------------------
847 template <typename TSpace>
848 void
849 DGtal::BoundedLatticePolytope<TSpace>::
850 getPointsByScanning( std::vector<Point>& pts ) const
851 {
852  pts.clear();
853  for ( const Point & p : D )
854  if ( isDomainPointInside( p ) ) pts.push_back( p );
855 }
856 //-----------------------------------------------------------------------------
857 template <typename TSpace>
858 template <typename PointSet>
859 void
860 DGtal::BoundedLatticePolytope<TSpace>::
861 insertPointsByScanning( PointSet& pts_set ) const
862 {
863  for ( const Point & p : D )
864  if ( isDomainPointInside( p ) ) pts_set.insert( p );
865 }
866 //-----------------------------------------------------------------------------
867 template <typename TSpace>
868 void
869 DGtal::BoundedLatticePolytope<TSpace>::
870 getInteriorPointsByScanning( std::vector<Point>& pts ) const
871 {
872  pts.clear();
873  for ( const Point & p : D )
874  if ( isInterior( p ) ) pts.push_back( p );
875 }
876 //-----------------------------------------------------------------------------
877 template <typename TSpace>
878 void
879 DGtal::BoundedLatticePolytope<TSpace>::
880 getBoundaryPointsByScanning( std::vector<Point>& pts ) const
881 {
882  pts.clear();
883  for ( const Point & p : D )
884  if ( isBoundary( p ) ) pts.push_back( p );
885 }
886 
887 //-----------------------------------------------------------------------------
888 template <typename TSpace>
889 const typename DGtal::BoundedLatticePolytope<TSpace>::Domain&
890 DGtal::BoundedLatticePolytope<TSpace>::getDomain() const
891 {
892  return D;
893 }
894 
895 //-----------------------------------------------------------------------------
896 template <typename TSpace>
897 unsigned int
898 DGtal::BoundedLatticePolytope<TSpace>::nbHalfSpaces() const
899 {
900  return static_cast<unsigned int>(A.size());
901 }
902 
903 //-----------------------------------------------------------------------------
904 template <typename TSpace>
905 const typename DGtal::BoundedLatticePolytope<TSpace>::Vector&
906 DGtal::BoundedLatticePolytope<TSpace>::getA( unsigned int i ) const
907 {
908  ASSERT( i < nbHalfSpaces() );
909  return A[ i ];
910 }
911 
912 //-----------------------------------------------------------------------------
913 template <typename TSpace>
914 typename DGtal::BoundedLatticePolytope<TSpace>::Integer
915 DGtal::BoundedLatticePolytope<TSpace>::getB( unsigned int i ) const
916 {
917  ASSERT( i < nbHalfSpaces() );
918  return B[ i ];
919 }
920 
921 //-----------------------------------------------------------------------------
922 template <typename TSpace>
923 bool
924 DGtal::BoundedLatticePolytope<TSpace>::isLarge( unsigned int i ) const
925 {
926  ASSERT( i < nbHalfSpaces() );
927  return I[ i ];
928 }
929 
930 //-----------------------------------------------------------------------------
931 template <typename TSpace>
932 const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityMatrix&
933 DGtal::BoundedLatticePolytope<TSpace>::getA() const
934 {
935  return A;
936 }
937 
938 //-----------------------------------------------------------------------------
939 template <typename TSpace>
940 const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityVector&
941 DGtal::BoundedLatticePolytope<TSpace>::getB() const
942 {
943  return B;
944 }
945 
946 //-----------------------------------------------------------------------------
947 template <typename TSpace>
948 const std::vector<bool>&
949 DGtal::BoundedLatticePolytope<TSpace>::getI() const
950 {
951  return I;
952 }
953 
954 //-----------------------------------------------------------------------------
955 template <typename TSpace>
956 bool
957 DGtal::BoundedLatticePolytope<TSpace>::canBeSummed() const
958 {
959  return myValidEdgeConstraints;
960 }
961 
962 ///////////////////////////////////////////////////////////////////////////////
963 // Interface - public :
964 
965 /**
966  * Writes/Displays the object on an output stream.
967  * @param out the output stream where the object is written.
968  */
969 template <typename TSpace>
970 inline
971 void
972 DGtal::BoundedLatticePolytope<TSpace>::selfDisplay ( std::ostream & out ) const
973 {
974  out << "[BoundedLatticePolytope<" << Space::dimension << "> A.rows=" << A.size()
975  << " valid_edge_constraints=" << myValidEdgeConstraints
976  << "]" << std::endl;
977  for ( Dimension i = 0; i < A.size(); ++i )
978  {
979  out << " [";
980  for ( Dimension j = 0; j < dimension; ++j )
981  out << " " << A[ i ][ j ];
982  out << " ] . x <= " << B[ i ] << std::endl;
983  }
984 }
985 
986 /**
987  * Checks the validity/consistency of the object.
988  * @return 'true' if the object is valid, 'false' otherwise.
989  */
990 template <typename TSpace>
991 inline
992 bool
993 DGtal::BoundedLatticePolytope<TSpace>::isValid() const
994 {
995  return ! D.isEmpty();
996 }
997 //-----------------------------------------------------------------------------
998 template <typename TSpace>
999 inline
1000 std::string
1001 DGtal::BoundedLatticePolytope<TSpace>::className
1002 () const
1003 {
1004  return "BoundedLatticePolytope";
1005 }
1006 
1007 
1008 
1009 ///////////////////////////////////////////////////////////////////////////////
1010 // Implementation of inline functions //
1011 
1012 //-----------------------------------------------------------------------------
1013 template <typename TSpace>
1014 inline
1015 std::ostream&
1016 DGtal::operator<< ( std::ostream & out,
1017  const BoundedLatticePolytope<TSpace> & object )
1018 {
1019  object.selfDisplay( out );
1020  return out;
1021 }
1022 //-----------------------------------------------------------------------------
1023 template <typename TSpace>
1024 DGtal::BoundedLatticePolytope<TSpace>
1025 DGtal::operator* ( typename BoundedLatticePolytope<TSpace>::Integer t,
1026  const BoundedLatticePolytope<TSpace> & P )
1027 {
1028  BoundedLatticePolytope<TSpace> Q = P;
1029  Q *= t;
1030  return Q;
1031 }
1032 //-----------------------------------------------------------------------------
1033 template <typename TSpace>
1034 DGtal::BoundedLatticePolytope<TSpace>
1035 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1036  typename BoundedLatticePolytope<TSpace>::UnitSegment s )
1037 {
1038  BoundedLatticePolytope<TSpace> Q = P;
1039  Q += s;
1040  return Q;
1041 }
1042 //-----------------------------------------------------------------------------
1043 template <typename TSpace>
1044 DGtal::BoundedLatticePolytope<TSpace>
1045 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1046  typename BoundedLatticePolytope<TSpace>::UnitCell c )
1047 {
1048  BoundedLatticePolytope<TSpace> Q = P;
1049  Q += c;
1050  return Q;
1051 }
1052 //-----------------------------------------------------------------------------
1053 template <typename TSpace>
1054 DGtal::BoundedLatticePolytope<TSpace>
1055 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1056  typename BoundedLatticePolytope<TSpace>::RightStrictUnitSegment s )
1057 {
1058  BoundedLatticePolytope<TSpace> Q = P;
1059  Q += s;
1060  return Q;
1061 }
1062 //-----------------------------------------------------------------------------
1063 template <typename TSpace>
1064 DGtal::BoundedLatticePolytope<TSpace>
1065 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1066  typename BoundedLatticePolytope<TSpace>::RightStrictUnitCell c )
1067 {
1068  BoundedLatticePolytope<TSpace> Q = P;
1069  Q += c;
1070  return Q;
1071 }
1072 //-----------------------------------------------------------------------------
1073 template <typename TSpace>
1074 DGtal::BoundedLatticePolytope<TSpace>
1075 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1076  typename BoundedLatticePolytope<TSpace>::LeftStrictUnitSegment s )
1077 {
1078  BoundedLatticePolytope<TSpace> Q = P;
1079  Q += s;
1080  return Q;
1081 }
1082 //-----------------------------------------------------------------------------
1083 template <typename TSpace>
1084 DGtal::BoundedLatticePolytope<TSpace>
1085 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1086  typename BoundedLatticePolytope<TSpace>::LeftStrictUnitCell c )
1087 {
1088  BoundedLatticePolytope<TSpace> Q = P;
1089  Q += c;
1090  return Q;
1091 }
1092 
1093 // //
1094 ///////////////////////////////////////////////////////////////////////////////