DGtal  1.2.0
testBoundedRationalPolytope.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include <vector>
33 #include <algorithm>
34 #include "DGtal/base/Common.h"
35 #include "DGtal/kernel/SpaceND.h"
36 #include "DGtal/geometry/volumes/BoundedRationalPolytope.h"
37 #include "DGtalCatch.h"
39 
40 using namespace std;
41 using namespace DGtal;
42 
43 
45 // Functions for testing class BoundedRationalPolytope.
47 
48 SCENARIO( "BoundedRationalPolytope< Z2 > unit tests", "[rational_polytope][2d]" )
49 {
50  typedef SpaceND<2,int> Space;
51  typedef Space::Point Point;
52  typedef Space::Vector Vector;
53  typedef BoundedRationalPolytope< Space > Polytope;
54 
55  GIVEN( "A triangle P at (0,0), (5/2,0), (0,7/2)" ) {
56  Point a( 0, 0 );
57  Point b( 5, 0 );
58  Point c( 0, 7 );
59  Polytope P { Point(2,2), a, b, c };
60  THEN( "It contains more than 3 integer points" ) {
61  REQUIRE( P.count() > 3 );
62  }
63  THEN( "It contains more points than its area" ) {
64  REQUIRE( P.count() > (5*7/8) );
65  }
66  THEN( "It satisfies #In(P) <= #Int(P) + #Bd(P)" ) {
67  auto nb = P.count();
68  auto nb_int = P.countInterior();
69  auto nb_bd = P.countBoundary();
70  CAPTURE( nb );
71  CAPTURE( nb_int );
72  CAPTURE( nb_bd );
73  REQUIRE( nb <= nb_int + nb_bd );
74  }
75  WHEN( "Cut by some half-space" ) {
76  Polytope Q = P;
77  Q.cut( Vector( -1, 1 ), 3 );
78  THEN( "It contains less points" ) {
79  REQUIRE( Q.count() < P.count() );
80  }
81  }
82  WHEN( "Multiplied by 4 as Q = 4 * P, it becomes a lattice polytope" ) {
83  Polytope Q = 4 * P;
84  THEN( "It satisfies Pick's formula, ie 2*Area(P) = 2*#Int(P) + #Bd(P) - 2" ) {
85  Polytope IntQ = Q.interiorPolytope();
86  auto nb_int = IntQ.count();
87  auto nb_bd = Q.count() - IntQ.count();
88  auto area2 = (5*2)*(7*2);
89  CAPTURE( Q );
90  CAPTURE( nb_int );
91  CAPTURE( nb_bd );
92  CAPTURE( area2 );
93  REQUIRE( area2 == 2*nb_int + nb_bd - 2 );
94  }
95  }
96  }
97  GIVEN( "A closed segment S at (4/2,0/2), (-8/2,-4/2)" ) {
98  Point a( 4, 0 );
99  Point b( -8, -4 );
100  Polytope P { Point(2,2), a, b };
101  THEN( "Its interior is empty #Int(P) == 0" ) {
102  auto nb_int = P.countInterior();
103  REQUIRE( nb_int == 0 );
104  }
105  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P) == #Bd(P) == 3" ) {
106  auto nb = P.count();
107  auto nb_int = P.countInterior();
108  auto nb_bd = P.countBoundary();
109  CAPTURE( nb );
110  CAPTURE( nb_int );
111  CAPTURE( nb_bd );
112  std::vector<Point> Ppts;
113  P.getPoints( Ppts );
114  CAPTURE( P );
115  CAPTURE( Ppts );
116  REQUIRE( nb_bd == 3 );
117  REQUIRE( nb == nb_int + nb_bd );
118  }
119  }
120  GIVEN( "A thin triangle P at (4/4,2/4), (2/4,4/4), (9/4,9/4)" ) {
121  Point a( 4, 2 );
122  Point b( 2, 4 );
123  Point c( 9, 9 );
124  Polytope P { Point(4,4), a, b, c };
125  THEN( "It contains 2 integer points" ) {
126  REQUIRE( P.count() == 2 );
127  }
128  THEN( "Its boundary is empty" ) {
129  REQUIRE( P.countBoundary() == 0 );
130  }
131  WHEN( "Multiplied by 4 as Q = 4 * P, it becomes a lattice polytope" ) {
132  Polytope Q = 4 * P;
133  THEN( "It satisfies Pick's formula, ie 2*Area(P) = 2*#Int(P) + #Bd(P) - 2" ) {
134  Polytope IntQ = Q.interiorPolytope();
135  auto nb_int = IntQ.count();
136  auto nb_bd = Q.count() - IntQ.count();
137  auto area2 = 24;
138  CAPTURE( Q );
139  CAPTURE( nb_int );
140  CAPTURE( nb_bd );
141  CAPTURE( area2 );
142  REQUIRE( area2 == 2*nb_int + nb_bd - 2 );
143  }
144  }
145  WHEN( "Multiplied by 10/3 as Q = 10/3 * P, it is a rational polytope" ) {
146  Polytope Q = Polytope::Rational( 10, 3 ) * P;
147  THEN( "It has a denominator 3 * gcd(4,10) == 6" ) {
148  REQUIRE( Q.denominator() == 6 );
149  }
150  THEN( "#( 3P Cap Z2 ) <= #( Q Cap Z2 ) < #( 4P Cap Z2 )" ) {
151  Polytope R = 3 * P;
152  Polytope S = 4 * P;
153  auto nbQ = Q.count();
154  auto nbR = R.count();
155  auto nbS = S.count();
156  CAPTURE( nbR );
157  CAPTURE( nbQ );
158  CAPTURE( nbS );
159  REQUIRE( nbR <= nbQ );
160  REQUIRE( nbQ <= nbS );
161  }
162  THEN( "6/5*Q is a polytope equal to 4*P." ) {
163  Polytope R = Polytope::Rational( 6, 5 ) * Q;
164  Polytope S = 4 * P;
165  std::vector<Point> Rpts, Spts;
166  R.getPoints( Rpts );
167  S.getPoints( Spts );
168  CAPTURE( Rpts );
169  CAPTURE( Spts );
170  REQUIRE( std::equal( Rpts.cbegin(), Rpts.cend(), Spts.cbegin() ) );
171  }
172  }
173  }
174 } // SCENARIO( "BoundedRationalPolytope< Z2 > unit tests", "[rational_polytope][2d]" )
175 
176 SCENARIO( "BoundedRationalPolytope< Z3 > unit tests", "[rational_polytope][3d]" )
177 {
178  typedef SpaceND<3,int> Space;
179  typedef Space::Point Point;
180  typedef Space::Vector Vector;
181  typedef BoundedRationalPolytope< Space > Polytope;
182 
183  GIVEN( "A tetrahedron P at (0,0,0), (5/2,0,0), (0,7/2,0), (0,0,3/2)" ) {
184  Point a( 0, 0, 0 );
185  Point b( 5, 0, 0 );
186  Point c( 0, 7, 0 );
187  Point d( 0, 0, 3 );
188  Polytope P { Point(2,2), a, b, c, d };
189  THEN( "It contains more than 3 integer points" ) {
190  REQUIRE( P.count() > 3 );
191  }
192  THEN( "It contains more points than its volume" ) {
193  REQUIRE( P.count() > (5*7*3/8/6) );
194  }
195  THEN( "It satisfies #In(P) <= #Int(P) + #Bd(P)" ) {
196  auto nb = P.count();
197  auto nb_int = P.countInterior();
198  auto nb_bd = P.countBoundary();
199  CAPTURE( nb );
200  CAPTURE( nb_int );
201  CAPTURE( nb_bd );
202  REQUIRE( nb <= nb_int + nb_bd );
203  }
204  WHEN( "Multiplied by 4 as Q = 2 * P, it becomes a lattice polytope" ) {
205  Polytope Q = 2 * P;
206  THEN( "Its denominator is 1" ) {
207  REQUIRE( Q.denominator() == 1 );
208  }
209  THEN( "It has more interior and boundary points than P" ) {
210  auto nb = P.count();
211  auto nb_int = P.countInterior();
212  auto nb_bd = P.countBoundary();
213  Polytope IntQ = Q.interiorPolytope();
214  auto nb_intQ = IntQ.count();
215  auto nbQ = Q.count();
216  auto nb_bdQ = nbQ - nb_intQ;
217  CAPTURE( Q );
218  REQUIRE( nb_int < nb_intQ );
219  REQUIRE( nb_bd < nb_bdQ );
220  REQUIRE( nb < nbQ );
221  }
222  }
223  WHEN( "Considering an increasing series f_i = 3/2, 2, 9/4, 3, 11/3, the number of inside points of f_i * P is increasing" ) {
224  Polytope Q1 = Polytope::Rational( 3 , 2 ) * P;
225  Polytope Q2 = Polytope::Rational( 2 , 1 ) * P;
226  Polytope Q3 = Polytope::Rational( 9 , 4 ) * P;
227  Polytope Q4 = Polytope::Rational( 3 , 1 ) * P;
228  Polytope Q5 = Polytope::Rational( 11, 3 ) * P;
229  auto nb = P .count();
230  auto nb1 = Q1.count();
231  auto nb2 = Q2.count();
232  auto nb3 = Q3.count();
233  auto nb4 = Q4.count();
234  auto nb5 = Q5.count();
235  REQUIRE( nb < nb1 );
236  REQUIRE( nb1 < nb2 );
237  REQUIRE( nb2 < nb3 );
238  REQUIRE( nb3 < nb4 );
239  REQUIRE( nb4 < nb5 );
240  }
241  }
242 }
243 
244 // //
Aim: Represents an nD rational polytope, i.e. a convex polyhedron bounded by vertices with rational c...
DGtal is the top-level namespace which contains all DGtal functions and types.
SCENARIO("BoundedRationalPolytope< Z2 > unit tests", "[rational_polytope][2d]")
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
REQUIRE(domain.isInside(aPoint))