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