DGtal 2.1.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
44template <typename TSpace>
45void
46DGtal::BoundedLatticePolytope<TSpace>::
47clear()
48{
49 A.clear();
50 B.clear();
51 I.clear();
52 myValidEdgeConstraints = false;
53 D = Domain( Point::zero, Point::zero );
54}
55
56//-----------------------------------------------------------------------------
57template <typename TSpace>
58DGtal::BoundedLatticePolytope<TSpace>::
59BoundedLatticePolytope( std::initializer_list<Point> l )
60{
61 myValidEdgeConstraints = false;
62 init( l.begin(), l.end() );
63}
64
65//-----------------------------------------------------------------------------
66template <typename TSpace>
67template <typename PointIterator>
68DGtal::BoundedLatticePolytope<TSpace>::
69BoundedLatticePolytope( PointIterator itB, PointIterator itE )
70{
71 myValidEdgeConstraints = false;
72 init( itB, itE );
73}
74
75//-----------------------------------------------------------------------------
76template <typename TSpace>
77template <typename HalfSpaceIterator>
78DGtal::BoundedLatticePolytope<TSpace>::
79BoundedLatticePolytope( 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//-----------------------------------------------------------------------------
90template <typename TSpace>
91template <typename HalfSpaceIterator>
92void
93DGtal::BoundedLatticePolytope<TSpace>::
94init( 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//-----------------------------------------------------------------------------
153template <typename TSpace>
154bool
155DGtal::BoundedLatticePolytope<TSpace>::
156internalInitFromTriangle3D( 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//-----------------------------------------------------------------------------
181template <typename TSpace>
182bool
183DGtal::BoundedLatticePolytope<TSpace>::
184internalInitFromSegment3D( 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//-----------------------------------------------------------------------------
206template <typename TSpace>
207bool
208DGtal::BoundedLatticePolytope<TSpace>::
209internalInitFromSegment2D( 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//-----------------------------------------------------------------------------
223template <typename TSpace>
224template <typename PointIterator>
225bool
226DGtal::BoundedLatticePolytope<TSpace>::
227init( 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//-----------------------------------------------------------------------------
320template <typename TSpace>
321DGtal::BoundedLatticePolytope<TSpace>
322DGtal::BoundedLatticePolytope<TSpace>::
323interiorPolytope() 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//-----------------------------------------------------------------------------
332template <typename TSpace>
333DGtal::BoundedLatticePolytope<TSpace>
334DGtal::BoundedLatticePolytope<TSpace>::
335closurePolytope() 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//-----------------------------------------------------------------------------
344template <typename TSpace>
345unsigned int
346DGtal::BoundedLatticePolytope<TSpace>::
347cut( 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//-----------------------------------------------------------------------------
362template <typename TSpace>
363unsigned int
364DGtal::BoundedLatticePolytope<TSpace>::
365cut( 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//-----------------------------------------------------------------------------
387template <typename TSpace>
388unsigned int
389DGtal::BoundedLatticePolytope<TSpace>::
390cut( 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//-----------------------------------------------------------------------------
398template <typename TSpace>
399void
400DGtal::BoundedLatticePolytope<TSpace>::
401swap( 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//-----------------------------------------------------------------------------
411template <typename TSpace>
412bool
413DGtal::BoundedLatticePolytope<TSpace>::
414isInside( 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//-----------------------------------------------------------------------------
429template <typename TSpace>
430bool
431DGtal::BoundedLatticePolytope<TSpace>::
432isDomainPointInside( 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//-----------------------------------------------------------------------------
447template <typename TSpace>
448bool
449DGtal::BoundedLatticePolytope<TSpace>::
450isInterior( 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//-----------------------------------------------------------------------------
462template <typename TSpace>
463bool
464DGtal::BoundedLatticePolytope<TSpace>::
465isBoundary( 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//-----------------------------------------------------------------------------
479template <typename TSpace>
480typename DGtal::BoundedLatticePolytope<TSpace>::Self&
481DGtal::BoundedLatticePolytope<TSpace>::
482operator*=( 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//-----------------------------------------------------------------------------
490template <typename TSpace>
491typename DGtal::BoundedLatticePolytope<TSpace>::Self&
492DGtal::BoundedLatticePolytope<TSpace>::
493operator+=( 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//-----------------------------------------------------------------------------
507template <typename TSpace>
508typename DGtal::BoundedLatticePolytope<TSpace>::Self&
509DGtal::BoundedLatticePolytope<TSpace>::
510operator+=( StrictUnitSegment s )
511{
512 for ( Dimension i = 0; i < A.size(); ++i )
513 {
514 if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
515 {
516 B[ i ] += A[ i ][ s.k ];
517 I[ i ] = false;
518 }
519 else if ( A[ i ][ s.k ] < NumberTraits<Integer>::ZERO )
520 I[ i ] = false;
521 }
522 Vector z = Vector::zero;
523 z[ s.k ] = NumberTraits<Integer>::ONE;
524 D = Domain( D.lowerBound(), D.upperBound() + z );
525 return *this;
526}
527
528//-----------------------------------------------------------------------------
529template <typename TSpace>
530typename DGtal::BoundedLatticePolytope<TSpace>::Self&
531DGtal::BoundedLatticePolytope<TSpace>::
532operator+=( LeftStrictUnitSegment s )
533{
534 I[ 2*s.k + 1 ] = false;
535 for ( Dimension i = 0; i < A.size(); ++i )
536 {
537 if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
538 B[ i ] += A[ i ][ s.k ];
539 if ( A[ i ][ s.k ] < NumberTraits<Integer>::ZERO )
540 I[ i ] = false;
541 }
542 Vector z = Vector::zero;
543 z[ s.k ] = NumberTraits<Integer>::ONE;
544 D = Domain( D.lowerBound() + z, D.upperBound() + z );
545 return *this;
546}
547
548//-----------------------------------------------------------------------------
549template <typename TSpace>
550typename DGtal::BoundedLatticePolytope<TSpace>::Self&
551DGtal::BoundedLatticePolytope<TSpace>::
552operator+=( RightStrictUnitSegment s )
553{
554 I[ 2*s.k ] = false;
555 for ( Dimension i = 0; i < A.size(); ++i )
556 {
557 if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO ) {
558 B[ i ] += A[ i ][ s.k ];
559 I[ i ] = false;
560 }
561 }
562 Vector z = Vector::zero;
563 z[ s.k ] = NumberTraits<Integer>::ONE;
564 return *this;
565}
566
567//-----------------------------------------------------------------------------
568template <typename TSpace>
569typename DGtal::BoundedLatticePolytope<TSpace>::Self&
570DGtal::BoundedLatticePolytope<TSpace>::
571operator+=( UnitCell c )
572{
573 for ( Dimension i = 0; i < c.dims.size(); ++i )
574 *this += UnitSegment( c.dims[ i ] );
575 return *this;
576}
577
578//-----------------------------------------------------------------------------
579template <typename TSpace>
580typename DGtal::BoundedLatticePolytope<TSpace>::Self&
581DGtal::BoundedLatticePolytope<TSpace>::
582operator+=( StrictUnitCell c )
583{
584 for ( Dimension i = 0; i < c.dims.size(); ++i )
585 *this += StrictUnitSegment( c.dims[ i ] );
586 return *this;
587}
588
589//-----------------------------------------------------------------------------
590template <typename TSpace>
591typename DGtal::BoundedLatticePolytope<TSpace>::Self&
592DGtal::BoundedLatticePolytope<TSpace>::
593operator+=( RightStrictUnitCell c )
594{
595 for ( Dimension i = 0; i < c.dims.size(); ++i )
596 *this += RightStrictUnitSegment( c.dims[ i ] );
597 return *this;
598}
599
600//-----------------------------------------------------------------------------
601template <typename TSpace>
602typename DGtal::BoundedLatticePolytope<TSpace>::Self&
603DGtal::BoundedLatticePolytope<TSpace>::
604operator+=( LeftStrictUnitCell c )
605{
606 for ( Dimension i = 0; i < c.dims.size(); ++i )
607 *this += LeftStrictUnitSegment( c.dims[ i ] );
608 return *this;
609}
610
611//-----------------------------------------------------------------------------
612template <typename TSpace>
613typename DGtal::BoundedLatticePolytope<TSpace>::Integer
614DGtal::BoundedLatticePolytope<TSpace>::
615count() const
616{
617 BoundedLatticePolytopeCounter<Space> C( *this );
618 return C.countAlongAxis( C.longestAxis() );
619}
620
621//-----------------------------------------------------------------------------
622template <typename TSpace>
623typename DGtal::BoundedLatticePolytope<TSpace>::Integer
624DGtal::BoundedLatticePolytope<TSpace>::
625countInterior() const
626{
627 BoundedLatticePolytopeCounter<Space> C( *this );
628 return C.countInteriorAlongAxis( C.longestAxis() );
629}
630//-----------------------------------------------------------------------------
631template <typename TSpace>
632typename DGtal::BoundedLatticePolytope<TSpace>::Integer
633DGtal::BoundedLatticePolytope<TSpace>::
634countBoundary() const
635{
636 const auto clP = closurePolytope();
637 return clP.count() - countInterior();
638}
639//-----------------------------------------------------------------------------
640template <typename TSpace>
641typename DGtal::BoundedLatticePolytope<TSpace>::Integer
642DGtal::BoundedLatticePolytope<TSpace>::
643countWithin( Point lo, Point hi ) const
644{
645 BoundedLatticePolytopeCounter<Space> C( *this );
646 Dimension b = 0;
647 lo = lo.sup( D.lowerBound() );
648 hi = hi.inf( D.upperBound() );
649 auto b_size = hi[ 0 ] - lo[ 0 ];
650 for ( Dimension a = 1; a < dimension; a++ )
651 {
652 const auto a_size = hi[ a ] - lo[ a ];
653 if ( b_size < a_size ) { b = a; b_size = a_size; }
654 }
655 hi[ b ] = lo[ b ];
656 Integer nb = 0;
657 Domain localD( lo, hi );
658 for ( auto&& p : localD )
659 {
660 auto II = C.intersectionIntervalAlongAxis( p, b );
661 nb += II.second - II.first;
662 }
663 return nb;
664}
665//-----------------------------------------------------------------------------
666template <typename TSpace>
667typename DGtal::BoundedLatticePolytope<TSpace>::Integer
668DGtal::BoundedLatticePolytope<TSpace>::
669countUpTo( Integer max) const
670{
671 BoundedLatticePolytopeCounter<Space> C( *this );
672 Dimension a = C.longestAxis();
673 Integer nb = 0;
674 Point lo = D.lowerBound();
675 Point hi = D.upperBound();
676 hi[ a ] = lo[ a ];
677 Domain localD( lo, hi );
678 for ( auto&& p : localD )
679 {
680 auto II = C.intersectionIntervalAlongAxis( p, a );
681 nb += II.second - II.first;
682 if ( nb >= max ) return max;
683 }
684 return nb;
685}
686//-----------------------------------------------------------------------------
687template <typename TSpace>
688void
689DGtal::BoundedLatticePolytope<TSpace>::
690getPoints( std::vector<Point>& pts ) const
691{
692 pts.clear();
693 BoundedLatticePolytopeCounter<Space> C( *this );
694 C.getPointsAlongAxis( pts, C.longestAxis() );
695}
696//-----------------------------------------------------------------------------
697template <typename TSpace>
698void
699DGtal::BoundedLatticePolytope<TSpace>::
700getKPoints( std::vector<Point>& pts, const Point& alpha_shift ) const
701{
702 pts.clear();
703 BoundedLatticePolytopeCounter<Space> C( *this );
704 Dimension a = C.longestAxis();
705 Integer nb = 0;
706 Point lo = D.lowerBound();
707 Point hi = D.upperBound();
708 hi[ a ] = lo[ a ];
709 Domain localD( lo, hi );
710 for ( auto&& p : localD )
711 {
712 auto II = C.intersectionIntervalAlongAxis( p, a );
713 Point q( 2*p - alpha_shift );
714 q[ a ] = 2*II.first - alpha_shift[ a ];
715 for ( Integer x = II.first; x < II.second; x++ )
716 {
717 pts.push_back( q );
718 q[ a ] += 2;
719 }
720 }
721}
722//-----------------------------------------------------------------------------
723template <typename TSpace>
724template <typename PointSet>
725void
726DGtal::BoundedLatticePolytope<TSpace>::
727insertPoints( PointSet& pts_set ) const
728{
729 std::vector<Point> pts;
730 getPoints( pts );
731 pts_set.insert( pts.cbegin(), pts.cend() );
732}
733//-----------------------------------------------------------------------------
734template <typename TSpace>
735template <typename PointSet>
736void
737DGtal::BoundedLatticePolytope<TSpace>::
738insertKPoints( PointSet& pts_set, const Point& alpha_shift ) const
739{
740 BoundedLatticePolytopeCounter<Space> C( *this );
741 Dimension a = C.longestAxis();
742 Integer nb = 0;
743 Point lo = D.lowerBound();
744 Point hi = D.upperBound();
745 hi[ a ] = lo[ a ];
746 Domain localD( lo, hi );
747 for ( auto&& p : localD )
748 {
749 auto II = C.intersectionIntervalAlongAxis( p, a );
750 Point q( 2*p - alpha_shift );
751 q[ a ] = 2*II.first - alpha_shift[ a ];
752 for ( Integer x = II.first; x < II.second; x++ )
753 {
754 pts_set.insert( q );
755 q[ a ] += 2;
756 }
757 }
758}
759//-----------------------------------------------------------------------------
760template <typename TSpace>
761void
762DGtal::BoundedLatticePolytope<TSpace>::
763getInteriorPoints( std::vector<Point>& pts ) const
764{
765 pts.clear();
766 BoundedLatticePolytopeCounter<Space> C( *this );
767 C.getInteriorPointsAlongAxis( pts, C.longestAxis() );
768}
769//-----------------------------------------------------------------------------
770template <typename TSpace>
771void
772DGtal::BoundedLatticePolytope<TSpace>::
773getBoundaryPoints( std::vector<Point>& pts ) const
774{
775 const auto clP = closurePolytope();
776 BoundedLatticePolytopeCounter<Space> C ( *this );
777 BoundedLatticePolytopeCounter<Space> clC( clP );
778 pts.clear();
779 const Dimension a = clC.longestAxis();
780 Point lo = clP.getDomain().lowerBound();
781 Point hi = clP.getDomain().upperBound();
782 hi[ a ] = lo[ a ];
783 Domain localD( lo, hi );
784 for ( auto&& p : localD )
785 {
786 auto II = C .interiorIntersectionIntervalAlongAxis( p, a );
787 auto clI = clC.intersectionIntervalAlongAxis( p, a );
788 auto nbI = II.second - II.first;
789 auto nbclI = clI.second - clI.first;
790 if ( nbI == nbclI ) continue;
791 if ( nbI > nbclI )
792 trace.error() << "BoundedLatticePolytope::getBoundaryPoints: bad count"
793 << std::endl;
794 Point q = p;
795 if ( nbI == 0 )
796 {
797 for ( Integer x = clI.first; x != clI.second; x++ )
798 {
799 q[ a ] = x;
800 pts.push_back( q );
801 }
802 }
803 else
804 {
805 if ( clI.first < II.first )
806 {
807 q[ a ] = clI.first;
808 pts.push_back( q );
809 }
810 if ( clI.second > II.second )
811 {
812 q[ a ] = II.second;
813 pts.push_back( q );
814 }
815 }
816 }
817}
818
819
820//-----------------------------------------------------------------------------
821template <typename TSpace>
822typename DGtal::BoundedLatticePolytope<TSpace>::Integer
823DGtal::BoundedLatticePolytope<TSpace>::
824countByScanning() const
825{
826 Integer nb = 0;
827 for ( const Point & p : D )
828 nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
829 return nb;
830}
831
832//-----------------------------------------------------------------------------
833template <typename TSpace>
834typename DGtal::BoundedLatticePolytope<TSpace>::Integer
835DGtal::BoundedLatticePolytope<TSpace>::
836countInteriorByScanning() const
837{
838 Integer nb = 0;
839 for ( const Point & p : D )
840 nb += isInterior( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
841 return nb;
842}
843//-----------------------------------------------------------------------------
844template <typename TSpace>
845typename DGtal::BoundedLatticePolytope<TSpace>::Integer
846DGtal::BoundedLatticePolytope<TSpace>::
847countBoundaryByScanning() const
848{
849 Integer nb = 0;
850 for ( const Point & p : D )
851 nb += isBoundary( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
852 return nb;
853}
854//-----------------------------------------------------------------------------
855template <typename TSpace>
856typename DGtal::BoundedLatticePolytope<TSpace>::Integer
857DGtal::BoundedLatticePolytope<TSpace>::
858countWithinByScanning( Point lo, Point hi ) const
859{
860 Integer nb = 0;
861 Domain D1( lo.sup( D.lowerBound() ), hi.inf( D.upperBound() ) );
862 for ( const Point & p : D1 )
863 nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
864 return nb;
865}
866//-----------------------------------------------------------------------------
867template <typename TSpace>
868typename DGtal::BoundedLatticePolytope<TSpace>::Integer
869DGtal::BoundedLatticePolytope<TSpace>::
870countUpToByScanning( Integer max) const
871{
872 Integer nb = 0;
873 for ( const Point & p : D ) {
874 nb += isDomainPointInside( p ) ? NumberTraits<Integer>::ONE : NumberTraits<Integer>::ZERO;
875 if ( nb >= max ) return max;
876 }
877 return nb;
878}
879//-----------------------------------------------------------------------------
880template <typename TSpace>
881void
882DGtal::BoundedLatticePolytope<TSpace>::
883getPointsByScanning( std::vector<Point>& pts ) const
884{
885 pts.clear();
886 for ( const Point & p : D )
887 if ( isDomainPointInside( p ) ) pts.push_back( p );
888}
889//-----------------------------------------------------------------------------
890template <typename TSpace>
891template <typename PointSet>
892void
893DGtal::BoundedLatticePolytope<TSpace>::
894insertPointsByScanning( PointSet& pts_set ) const
895{
896 for ( const Point & p : D )
897 if ( isDomainPointInside( p ) ) pts_set.insert( p );
898}
899//-----------------------------------------------------------------------------
900template <typename TSpace>
901void
902DGtal::BoundedLatticePolytope<TSpace>::
903getInteriorPointsByScanning( std::vector<Point>& pts ) const
904{
905 pts.clear();
906 for ( const Point & p : D )
907 if ( isInterior( p ) ) pts.push_back( p );
908}
909//-----------------------------------------------------------------------------
910template <typename TSpace>
911void
912DGtal::BoundedLatticePolytope<TSpace>::
913getBoundaryPointsByScanning( std::vector<Point>& pts ) const
914{
915 pts.clear();
916 for ( const Point & p : D )
917 if ( isBoundary( p ) ) pts.push_back( p );
918}
919
920//-----------------------------------------------------------------------------
921template <typename TSpace>
922const typename DGtal::BoundedLatticePolytope<TSpace>::Domain&
923DGtal::BoundedLatticePolytope<TSpace>::getDomain() const
924{
925 return D;
926}
927
928//-----------------------------------------------------------------------------
929template <typename TSpace>
930unsigned int
931DGtal::BoundedLatticePolytope<TSpace>::nbHalfSpaces() const
932{
933 return static_cast<unsigned int>(A.size());
934}
935
936//-----------------------------------------------------------------------------
937template <typename TSpace>
938const typename DGtal::BoundedLatticePolytope<TSpace>::Vector&
939DGtal::BoundedLatticePolytope<TSpace>::getA( unsigned int i ) const
940{
941 ASSERT( i < nbHalfSpaces() );
942 return A[ i ];
943}
944
945//-----------------------------------------------------------------------------
946template <typename TSpace>
947typename DGtal::BoundedLatticePolytope<TSpace>::Integer
948DGtal::BoundedLatticePolytope<TSpace>::getB( unsigned int i ) const
949{
950 ASSERT( i < nbHalfSpaces() );
951 return B[ i ];
952}
953
954//-----------------------------------------------------------------------------
955template <typename TSpace>
956bool
957DGtal::BoundedLatticePolytope<TSpace>::isLarge( unsigned int i ) const
958{
959 ASSERT( i < nbHalfSpaces() );
960 return I[ i ];
961}
962
963//-----------------------------------------------------------------------------
964template <typename TSpace>
965void
966DGtal::BoundedLatticePolytope<TSpace>::setLarge( unsigned int i )
967{
968 ASSERT( i < nbHalfSpaces() );
969 I[ i ] = true;
970}
971
972//-----------------------------------------------------------------------------
973template <typename TSpace>
974void
975DGtal::BoundedLatticePolytope<TSpace>::setStrict( unsigned int i )
976{
977 ASSERT( i < nbHalfSpaces() );
978 I[ i ] = false;
979}
980
981
982//-----------------------------------------------------------------------------
983template <typename TSpace>
984const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityMatrix&
985DGtal::BoundedLatticePolytope<TSpace>::getA() const
986{
987 return A;
988}
989
990//-----------------------------------------------------------------------------
991template <typename TSpace>
992const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityVector&
993DGtal::BoundedLatticePolytope<TSpace>::getB() const
994{
995 return B;
996}
997
998//-----------------------------------------------------------------------------
999template <typename TSpace>
1000const std::vector<bool>&
1001DGtal::BoundedLatticePolytope<TSpace>::getI() const
1002{
1003 return I;
1004}
1005
1006//-----------------------------------------------------------------------------
1007template <typename TSpace>
1008bool
1009DGtal::BoundedLatticePolytope<TSpace>::canBeSummed() const
1010{
1011 return myValidEdgeConstraints;
1012}
1013
1014///////////////////////////////////////////////////////////////////////////////
1015// Interface - public :
1016
1017/**
1018 * Writes/Displays the object on an output stream.
1019 * @param out the output stream where the object is written.
1020 */
1021template <typename TSpace>
1022inline
1023void
1024DGtal::BoundedLatticePolytope<TSpace>::selfDisplay ( std::ostream & out ) const
1025{
1026 out << "[BoundedLatticePolytope<" << Space::dimension << "> A.rows=" << A.size()
1027 << " valid_edge_constraints=" << myValidEdgeConstraints
1028 << "]" << std::endl;
1029 for ( Dimension i = 0; i < A.size(); ++i )
1030 {
1031 out << " [";
1032 for ( Dimension j = 0; j < dimension; ++j )
1033 out << " " << A[ i ][ j ];
1034 out << " ] . x " << ( isLarge( i ) ? "<=" : "<" ) << " " << B[ i ] << std::endl;
1035 }
1036}
1037
1038/**
1039 * Checks the validity/consistency of the object.
1040 * @return 'true' if the object is valid, 'false' otherwise.
1041 */
1042template <typename TSpace>
1043inline
1044bool
1045DGtal::BoundedLatticePolytope<TSpace>::isValid() const
1046{
1047 return ! D.isEmpty();
1048}
1049//-----------------------------------------------------------------------------
1050template <typename TSpace>
1051inline
1052std::string
1053DGtal::BoundedLatticePolytope<TSpace>::className
1054() const
1055{
1056 return "BoundedLatticePolytope";
1057}
1058
1059
1060
1061///////////////////////////////////////////////////////////////////////////////
1062// Implementation of inline functions //
1063
1064//-----------------------------------------------------------------------------
1065template <typename TSpace>
1066inline
1067std::ostream&
1068DGtal::operator<< ( std::ostream & out,
1069 const BoundedLatticePolytope<TSpace> & object )
1070{
1071 object.selfDisplay( out );
1072 return out;
1073}
1074//-----------------------------------------------------------------------------
1075template <typename TSpace>
1076DGtal::BoundedLatticePolytope<TSpace>
1077DGtal::operator* ( typename BoundedLatticePolytope<TSpace>::Integer t,
1078 const BoundedLatticePolytope<TSpace> & P )
1079{
1080 BoundedLatticePolytope<TSpace> Q = P;
1081 Q *= t;
1082 return Q;
1083}
1084//-----------------------------------------------------------------------------
1085template <typename TSpace>
1086DGtal::BoundedLatticePolytope<TSpace>
1087DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1088 typename BoundedLatticePolytope<TSpace>::UnitSegment s )
1089{
1090 BoundedLatticePolytope<TSpace> Q = P;
1091 Q += s;
1092 return Q;
1093}
1094//-----------------------------------------------------------------------------
1095template <typename TSpace>
1096DGtal::BoundedLatticePolytope<TSpace>
1097DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1098 typename BoundedLatticePolytope<TSpace>::UnitCell c )
1099{
1100 BoundedLatticePolytope<TSpace> Q = P;
1101 Q += c;
1102 return Q;
1103}
1104//-----------------------------------------------------------------------------
1105template <typename TSpace>
1106DGtal::BoundedLatticePolytope<TSpace>
1107DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1108 typename BoundedLatticePolytope<TSpace>::StrictUnitSegment s )
1109{
1110 BoundedLatticePolytope<TSpace> Q = P;
1111 Q += s;
1112 return Q;
1113}
1114//-----------------------------------------------------------------------------
1115template <typename TSpace>
1116DGtal::BoundedLatticePolytope<TSpace>
1117DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1118 typename BoundedLatticePolytope<TSpace>::StrictUnitCell c )
1119{
1120 BoundedLatticePolytope<TSpace> Q = P;
1121 Q += c;
1122 return Q;
1123}
1124//-----------------------------------------------------------------------------
1125template <typename TSpace>
1126DGtal::BoundedLatticePolytope<TSpace>
1127DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1128 typename BoundedLatticePolytope<TSpace>::RightStrictUnitSegment s )
1129{
1130 BoundedLatticePolytope<TSpace> Q = P;
1131 Q += s;
1132 return Q;
1133}
1134//-----------------------------------------------------------------------------
1135template <typename TSpace>
1136DGtal::BoundedLatticePolytope<TSpace>
1137DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1138 typename BoundedLatticePolytope<TSpace>::RightStrictUnitCell c )
1139{
1140 BoundedLatticePolytope<TSpace> Q = P;
1141 Q += c;
1142 return Q;
1143}
1144//-----------------------------------------------------------------------------
1145template <typename TSpace>
1146DGtal::BoundedLatticePolytope<TSpace>
1147DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1148 typename BoundedLatticePolytope<TSpace>::LeftStrictUnitSegment s )
1149{
1150 BoundedLatticePolytope<TSpace> Q = P;
1151 Q += s;
1152 return Q;
1153}
1154//-----------------------------------------------------------------------------
1155template <typename TSpace>
1156DGtal::BoundedLatticePolytope<TSpace>
1157DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
1158 typename BoundedLatticePolytope<TSpace>::LeftStrictUnitCell c )
1159{
1160 BoundedLatticePolytope<TSpace> Q = P;
1161 Q += c;
1162 return Q;
1163}
1164
1165// //
1166///////////////////////////////////////////////////////////////////////////////