DGtal 1.3.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 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 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+=( 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//-----------------------------------------------------------------------------
527template <typename TSpace>
528typename DGtal::BoundedLatticePolytope<TSpace>::Self&
529DGtal::BoundedLatticePolytope<TSpace>::
530operator+=( 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//-----------------------------------------------------------------------------
546template <typename TSpace>
547typename DGtal::BoundedLatticePolytope<TSpace>::Self&
548DGtal::BoundedLatticePolytope<TSpace>::
549operator+=( 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//-----------------------------------------------------------------------------
557template <typename TSpace>
558typename DGtal::BoundedLatticePolytope<TSpace>::Self&
559DGtal::BoundedLatticePolytope<TSpace>::
560operator+=( 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//-----------------------------------------------------------------------------
568template <typename TSpace>
569typename DGtal::BoundedLatticePolytope<TSpace>::Self&
570DGtal::BoundedLatticePolytope<TSpace>::
571operator+=( 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//-----------------------------------------------------------------------------
579template <typename TSpace>
580typename DGtal::BoundedLatticePolytope<TSpace>::Integer
581DGtal::BoundedLatticePolytope<TSpace>::
582count() const
583{
584 BoundedLatticePolytopeCounter<Space> C( *this );
585 return C.countAlongAxis( C.longestAxis() );
586}
587
588//-----------------------------------------------------------------------------
589template <typename TSpace>
590typename DGtal::BoundedLatticePolytope<TSpace>::Integer
591DGtal::BoundedLatticePolytope<TSpace>::
592countInterior() const
593{
594 BoundedLatticePolytopeCounter<Space> C( *this );
595 return C.countInteriorAlongAxis( C.longestAxis() );
596}
597//-----------------------------------------------------------------------------
598template <typename TSpace>
599typename DGtal::BoundedLatticePolytope<TSpace>::Integer
600DGtal::BoundedLatticePolytope<TSpace>::
601countBoundary() const
602{
603 const auto clP = closurePolytope();
604 return clP.count() - countInterior();
605}
606//-----------------------------------------------------------------------------
607template <typename TSpace>
608typename DGtal::BoundedLatticePolytope<TSpace>::Integer
609DGtal::BoundedLatticePolytope<TSpace>::
610countWithin( 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 I = C.intersectionIntervalAlongAxis( p, b );
628 nb += I.second - I.first;
629 }
630 return nb;
631}
632//-----------------------------------------------------------------------------
633template <typename TSpace>
634typename DGtal::BoundedLatticePolytope<TSpace>::Integer
635DGtal::BoundedLatticePolytope<TSpace>::
636countUpTo( 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 I = C.intersectionIntervalAlongAxis( p, a );
648 nb += I.second - I.first;
649 if ( nb >= max ) return max;
650 }
651 return nb;
652}
653//-----------------------------------------------------------------------------
654template <typename TSpace>
655void
656DGtal::BoundedLatticePolytope<TSpace>::
657getPoints( std::vector<Point>& pts ) const
658{
659 pts.clear();
660 BoundedLatticePolytopeCounter<Space> C( *this );
661 C.getPointsAlongAxis( pts, C.longestAxis() );
662}
663//-----------------------------------------------------------------------------
664template <typename TSpace>
665void
666DGtal::BoundedLatticePolytope<TSpace>::
667getKPoints( 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 I = C.intersectionIntervalAlongAxis( p, a );
680 Point q( 2*p - alpha_shift );
681 q[ a ] = 2*I.first - alpha_shift[ a ];
682 for ( Integer x = I.first; x < I.second; x++ )
683 {
684 pts.push_back( q );
685 q[ a ] += 2;
686 }
687 }
688}
689//-----------------------------------------------------------------------------
690template <typename TSpace>
691template <typename PointSet>
692void
693DGtal::BoundedLatticePolytope<TSpace>::
694insertPoints( PointSet& pts_set ) const
695{
696 std::vector<Point> pts;
697 getPoints( pts );
698 pts_set.insert( pts.cbegin(), pts.cend() );
699}
700//-----------------------------------------------------------------------------
701template <typename TSpace>
702template <typename PointSet>
703void
704DGtal::BoundedLatticePolytope<TSpace>::
705insertKPoints( 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 I = C.intersectionIntervalAlongAxis( p, a );
717 Point q( 2*p - alpha_shift );
718 q[ a ] = 2*I.first - alpha_shift[ a ];
719 for ( Integer x = I.first; x < I.second; x++ )
720 {
721 pts_set.insert( q );
722 q[ a ] += 2;
723 }
724 }
725}
726//-----------------------------------------------------------------------------
727template <typename TSpace>
728void
729DGtal::BoundedLatticePolytope<TSpace>::
730getInteriorPoints( std::vector<Point>& pts ) const
731{
732 pts.clear();
733 BoundedLatticePolytopeCounter<Space> C( *this );
734 C.getInteriorPointsAlongAxis( pts, C.longestAxis() );
735}
736//-----------------------------------------------------------------------------
737template <typename TSpace>
738void
739DGtal::BoundedLatticePolytope<TSpace>::
740getBoundaryPoints( 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 I = C .interiorIntersectionIntervalAlongAxis( p, a );
754 auto clI = clC.intersectionIntervalAlongAxis( p, a );
755 auto nbI = I.second - I.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 < I.first )
773 {
774 q[ a ] = clI.first;
775 pts.push_back( q );
776 }
777 if ( clI.second > I.second )
778 {
779 q[ a ] = I.second;
780 pts.push_back( q );
781 }
782 }
783 }
784}
785
786
787//-----------------------------------------------------------------------------
788template <typename TSpace>
789typename DGtal::BoundedLatticePolytope<TSpace>::Integer
790DGtal::BoundedLatticePolytope<TSpace>::
791countByScanning() 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//-----------------------------------------------------------------------------
800template <typename TSpace>
801typename DGtal::BoundedLatticePolytope<TSpace>::Integer
802DGtal::BoundedLatticePolytope<TSpace>::
803countInteriorByScanning() 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//-----------------------------------------------------------------------------
811template <typename TSpace>
812typename DGtal::BoundedLatticePolytope<TSpace>::Integer
813DGtal::BoundedLatticePolytope<TSpace>::
814countBoundaryByScanning() 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//-----------------------------------------------------------------------------
822template <typename TSpace>
823typename DGtal::BoundedLatticePolytope<TSpace>::Integer
824DGtal::BoundedLatticePolytope<TSpace>::
825countWithinByScanning( 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//-----------------------------------------------------------------------------
834template <typename TSpace>
835typename DGtal::BoundedLatticePolytope<TSpace>::Integer
836DGtal::BoundedLatticePolytope<TSpace>::
837countUpToByScanning( 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//-----------------------------------------------------------------------------
847template <typename TSpace>
848void
849DGtal::BoundedLatticePolytope<TSpace>::
850getPointsByScanning( 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//-----------------------------------------------------------------------------
857template <typename TSpace>
858template <typename PointSet>
859void
860DGtal::BoundedLatticePolytope<TSpace>::
861insertPointsByScanning( PointSet& pts_set ) const
862{
863 for ( const Point & p : D )
864 if ( isDomainPointInside( p ) ) pts_set.insert( p );
865}
866//-----------------------------------------------------------------------------
867template <typename TSpace>
868void
869DGtal::BoundedLatticePolytope<TSpace>::
870getInteriorPointsByScanning( 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//-----------------------------------------------------------------------------
877template <typename TSpace>
878void
879DGtal::BoundedLatticePolytope<TSpace>::
880getBoundaryPointsByScanning( 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//-----------------------------------------------------------------------------
888template <typename TSpace>
889const typename DGtal::BoundedLatticePolytope<TSpace>::Domain&
890DGtal::BoundedLatticePolytope<TSpace>::getDomain() const
891{
892 return D;
893}
894
895//-----------------------------------------------------------------------------
896template <typename TSpace>
897unsigned int
898DGtal::BoundedLatticePolytope<TSpace>::nbHalfSpaces() const
899{
900 return A.size();
901}
902
903//-----------------------------------------------------------------------------
904template <typename TSpace>
905const typename DGtal::BoundedLatticePolytope<TSpace>::Vector&
906DGtal::BoundedLatticePolytope<TSpace>::getA( unsigned int i ) const
907{
908 ASSERT( i < nbHalfSpaces() );
909 return A[ i ];
910}
911
912//-----------------------------------------------------------------------------
913template <typename TSpace>
914typename DGtal::BoundedLatticePolytope<TSpace>::Integer
915DGtal::BoundedLatticePolytope<TSpace>::getB( unsigned int i ) const
916{
917 ASSERT( i < nbHalfSpaces() );
918 return B[ i ];
919}
920
921//-----------------------------------------------------------------------------
922template <typename TSpace>
923bool
924DGtal::BoundedLatticePolytope<TSpace>::isLarge( unsigned int i ) const
925{
926 ASSERT( i < nbHalfSpaces() );
927 return I[ i ];
928}
929
930//-----------------------------------------------------------------------------
931template <typename TSpace>
932const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityMatrix&
933DGtal::BoundedLatticePolytope<TSpace>::getA() const
934{
935 return A;
936}
937
938//-----------------------------------------------------------------------------
939template <typename TSpace>
940const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityVector&
941DGtal::BoundedLatticePolytope<TSpace>::getB() const
942{
943 return B;
944}
945
946//-----------------------------------------------------------------------------
947template <typename TSpace>
948const std::vector<bool>&
949DGtal::BoundedLatticePolytope<TSpace>::getI() const
950{
951 return I;
952}
953
954//-----------------------------------------------------------------------------
955template <typename TSpace>
956bool
957DGtal::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 */
969template <typename TSpace>
970inline
971void
972DGtal::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 */
990template <typename TSpace>
991inline
992bool
993DGtal::BoundedLatticePolytope<TSpace>::isValid() const
994{
995 return ! D.isEmpty();
996}
997//-----------------------------------------------------------------------------
998template <typename TSpace>
999inline
1000std::string
1001DGtal::BoundedLatticePolytope<TSpace>::className
1002() const
1003{
1004 return "BoundedLatticePolytope";
1005}
1006
1007
1008
1009///////////////////////////////////////////////////////////////////////////////
1010// Implementation of inline functions //
1011
1012//-----------------------------------------------------------------------------
1013template <typename TSpace>
1014inline
1015std::ostream&
1016DGtal::operator<< ( std::ostream & out,
1017 const BoundedLatticePolytope<TSpace> & object )
1018{
1019 object.selfDisplay( out );
1020 return out;
1021}
1022//-----------------------------------------------------------------------------
1023template <typename TSpace>
1024DGtal::BoundedLatticePolytope<TSpace>
1025DGtal::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//-----------------------------------------------------------------------------
1033template <typename TSpace>
1034DGtal::BoundedLatticePolytope<TSpace>
1035DGtal::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//-----------------------------------------------------------------------------
1043template <typename TSpace>
1044DGtal::BoundedLatticePolytope<TSpace>
1045DGtal::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//-----------------------------------------------------------------------------
1053template <typename TSpace>
1054DGtal::BoundedLatticePolytope<TSpace>
1055DGtal::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//-----------------------------------------------------------------------------
1063template <typename TSpace>
1064DGtal::BoundedLatticePolytope<TSpace>
1065DGtal::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//-----------------------------------------------------------------------------
1073template <typename TSpace>
1074DGtal::BoundedLatticePolytope<TSpace>
1075DGtal::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//-----------------------------------------------------------------------------
1083template <typename TSpace>
1084DGtal::BoundedLatticePolytope<TSpace>
1085DGtal::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///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:40 for DGtal by doxygen 1.9.5