DGtal  1.5.beta
testBoundedLatticePolytope.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/BoundedLatticePolytope.h"
37 #include "DGtalCatch.h"
39 
40 using namespace std;
41 using namespace DGtal;
42 
43 
45 // Functions for testing class BoundedLatticePolytope.
47 
48 SCENARIO( "BoundedLatticePolytope< Z2 > unit tests", "[lattice_polytope][2d]" )
49 {
50  typedef SpaceND<2,int> Space;
51  typedef Space::Point Point;
52  typedef Space::Vector Vector;
53  typedef BoundedLatticePolytope< Space > Polytope;
54 
55  GIVEN( "A triangle P at (0,0), (5,0), (0,7)" ) {
56  Point a( 0, 0 );
57  Point b( 5, 0 );
58  Point c( 0, 7 );
59  Polytope P { a, b, c };
60  THEN( "Its domain contains its vertices" ) {
61  REQUIRE( P.isDomainPointInside( a ) );
62  REQUIRE( P.isDomainPointInside( b ) );
63  REQUIRE( P.isDomainPointInside( c ) );
64  }
65  THEN( "Its vertices lie on its boundary" ) {
66  REQUIRE( P.isBoundary( a ) );
67  REQUIRE( P.isBoundary( b ) );
68  REQUIRE( P.isBoundary( c ) );
69  REQUIRE( ! P.isInterior( a ) );
70  REQUIRE( ! P.isInterior( b ) );
71  REQUIRE( ! P.isInterior( c ) );
72  }
73  THEN( "It contains more than 3 integer points" ) {
74  REQUIRE( P.count() > 3 );
75  }
76  THEN( "It contains more points than its area" ) {
77  REQUIRE( P.count() > (5*7/2) );
78  }
79  THEN( "It satisfies Pick's formula, ie 2*Area(P) = 2*#Int(P) + #Bd(P) - 2" ) {
80  Polytope IntP = P.interiorPolytope();
81  auto nb_int = IntP.count();
82  auto nb_bd = P.count() - IntP.count();
83  auto area2 = 5*7;
84  CAPTURE( nb_int );
85  CAPTURE( nb_bd );
86  CAPTURE( area2 );
87  REQUIRE( area2 == 2*nb_int + nb_bd - 2 );
88  }
89  THEN( "It satisfies #In(P) <= #Int(P) + #Bd(P)" ) {
90  auto nb = P.count();
91  auto nb_int = P.countInterior();
92  auto nb_bd = P.countBoundary();
93  CAPTURE( nb );
94  CAPTURE( nb_int );
95  CAPTURE( nb_bd );
96  REQUIRE( nb <= nb_int + nb_bd );
97  }
98  WHEN( "Cut by some half-space" ) {
99  Polytope Q = P;
100  Q.cut( Vector( -1, 1 ), 3 );
101  THEN( "It contains less points" ) {
102  REQUIRE( Q.count() < P.count() );
103  }
104  }
105  THEN( "Its boundary points and interior points are its inside points (closed polytope)" ) {
106  std::vector<Point> inside;
107  std::vector<Point> interior;
108  std::vector<Point> boundary;
109  std::vector<Point> all;
110  P.getPoints( inside );
111  P.getInteriorPoints( interior );
112  P.getBoundaryPoints( boundary );
113  std::sort( inside.begin(), inside.end() );
114  std::sort( interior.begin(), interior.end() );
115  std::sort( boundary.begin(), boundary.end() );
116  CAPTURE( inside );
117  CAPTURE( interior );
118  CAPTURE( boundary );
119  std::set_union( interior.cbegin(), interior.cend(), boundary.cbegin(), boundary.cend(),
120  std::back_inserter( all ) );
121  REQUIRE( std::equal( inside.cbegin(), inside.cend(), all.cbegin() ) );
122  }
123  }
124  GIVEN( "A closed segment S at (4,0), (-8,-4)" ) {
125  Point a( 4, 0 );
126  Point b( -8, -4 );
127  Polytope P { a, b };
128  THEN( "Its interior is empty #Int(P) == 0" ) {
129  auto nb_int = P.countInterior();
130  REQUIRE( nb_int == 0 );
131  }
132  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P) == #Bd(P) == 5" ) {
133  auto nb = P.count();
134  auto nb_int = P.countInterior();
135  auto nb_bd = P.countBoundary();
136  CAPTURE( nb );
137  CAPTURE( nb_int );
138  CAPTURE( nb_bd );
139  std::vector<Point> Ppts;
140  P.getPoints( Ppts );
141  CAPTURE( Ppts );
142  REQUIRE( nb_bd == 5 );
143  REQUIRE( nb == nb_int + nb_bd );
144  }
145  }
146 }
147 
148 SCENARIO( "BoundedLatticePolytope< Z3 > unit tests", "[lattice_polytope][3d]" )
149 {
150  typedef SpaceND<3,int> Space;
151  typedef Space::Point Point;
152  typedef BoundedLatticePolytope< Space > Polytope;
153 
154  GIVEN( "A twisted simplex P at (0,0,0), (1,0,0), (0,1,0), (1,1,z)" ) {
155  Point a( 0, 0, 0 );
156  Point b( 1, 0, 0 );
157  Point c( 0, 1, 0 );
158  Point d( 1, 1, 8 );
159  Polytope P { a, b, c, d };
160  THEN( "Its domain contains its vertices" ) {
161  REQUIRE( P.isDomainPointInside( a ) );
162  REQUIRE( P.isDomainPointInside( b ) );
163  REQUIRE( P.isDomainPointInside( c ) );
164  REQUIRE( P.isDomainPointInside( d ) );
165  }
166  THEN( "Its vertices lie on its boundary" ) {
167  REQUIRE( P.isBoundary( a ) );
168  REQUIRE( P.isBoundary( b ) );
169  REQUIRE( P.isBoundary( c ) );
170  REQUIRE( P.isBoundary( d ) );
171  REQUIRE( ! P.isInterior( a ) );
172  REQUIRE( ! P.isInterior( b ) );
173  REQUIRE( ! P.isInterior( c ) );
174  REQUIRE( ! P.isInterior( d ) );
175  }
176  THEN( "It satisfies #In(P) <= #Int(P) + #Bd(P)" ) {
177  auto nb = P.count();
178  auto nb_int = P.countInterior();
179  auto nb_bd = P.countBoundary();
180  CAPTURE( nb );
181  CAPTURE( nb_int );
182  CAPTURE( nb_bd );
183  REQUIRE( nb <= nb_int + nb_bd );
184  }
185  THEN( "It contains only 4 integer points" ) {
186  REQUIRE( P.count() == 4 );
187  }
188 
189  }
190 
191  GIVEN( "A closed arbitrary simplex P at (0,0,0), (6,3,0), (0,5,10), (6,4,8)" ) {
192  Point a( 0, 0, 0 );
193  Point b( 6, 3, 0 );
194  Point c( 0, 5, 10 );
195  Point d( 6, 4, 8 );
196  Polytope P { a, b, c, d };
197 
198  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P)" ) {
199  auto nb = P.count();
200  auto nb_int = P.countInterior();
201  auto nb_bd = P.countBoundary();
202  CAPTURE( nb );
203  CAPTURE( nb_int );
204  CAPTURE( nb_bd );
205  REQUIRE( nb == nb_int + nb_bd );
206  }
207  THEN( "Its boundary points and interior points are its inside points (closed polytope)" ) {
208  std::vector<Point> inside;
209  std::vector<Point> interior;
210  std::vector<Point> boundary;
211  std::vector<Point> all;
212  P.getPoints( inside );
213  P.getInteriorPoints( interior );
214  P.getBoundaryPoints( boundary );
215  std::sort( inside.begin(), inside.end() );
216  std::sort( interior.begin(), interior.end() );
217  std::sort( boundary.begin(), boundary.end() );
218  CAPTURE( inside );
219  CAPTURE( interior );
220  CAPTURE( boundary );
221  std::set_union( interior.cbegin(), interior.cend(), boundary.cbegin(), boundary.cend(),
222  std::back_inserter( all ) );
223  REQUIRE( std::equal( inside.cbegin(), inside.cend(), all.cbegin() ) );
224  }
225  WHEN( "Cut by some axis aligned half-space (1,0,0).x <= 3" ) {
226  Polytope Q = P;
227  Q.cut( 0, true, 3 );
228  THEN( "It contains less points" ) {
229  CAPTURE( P );
230  CAPTURE( P.getDomain() );
231  CAPTURE( Q );
232  CAPTURE( Q.getDomain() );
233  auto Pnb = P.count();
234  auto Pnb_int = P.countInterior();
235  auto Pnb_bd = P.countBoundary();
236  CAPTURE( Pnb );
237  CAPTURE( Pnb_int );
238  CAPTURE( Pnb_bd );
239  auto Qnb = Q.count();
240  auto Qnb_int = Q.countInterior();
241  auto Qnb_bd = Q.countBoundary();
242  CAPTURE( Qnb );
243  CAPTURE( Qnb_int );
244  CAPTURE( Qnb_bd );
245  std::vector<Point> Qpts;
246  Q.getPoints( Qpts );
247  CAPTURE( Qpts );
248  REQUIRE( Q.count() < P.count() );
249  }
250  }
251  }
252  GIVEN( "A closed triangle P at (0,0,0), (6,3,0), (2,5,10)" ) {
253  Point a( 0, 0, 0 );
254  Point b( 6, 3, 0 );
255  Point c( 2, 5, 10 );
256  Polytope P { a, b, c };
257  THEN( "Its interior is empty #Int(P) == 0" ) {
258  auto nb_int = P.countInterior();
259  REQUIRE( nb_int == 0 );
260  }
261  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P)" ) {
262  auto nb = P.count();
263  auto nb_int = P.countInterior();
264  auto nb_bd = P.countBoundary();
265  CAPTURE( nb );
266  CAPTURE( nb_int );
267  CAPTURE( nb_bd );
268  std::vector<Point> Ppts;
269  P.getPoints( Ppts );
270  CAPTURE( Ppts );
271  REQUIRE( nb == nb_int + nb_bd );
272  }
273  }
274  GIVEN( "A closed triangle P with relatively prime coordinates (3,0,0), (0,4,0), (0,0,5)" ) {
275  Point a( 3, 0, 0 );
276  Point b( 0, 4, 0 );
277  Point c( 0, 0, 5 );
278  Polytope P { a, b, c };
279  THEN( "Its interior is empty #Int(P) == 0" ) {
280  auto nb_int = P.countInterior();
281  REQUIRE( nb_int == 0 );
282  }
283  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P) == #Bd(P) == 3" ) {
284  auto nb = P.count();
285  auto nb_int = P.countInterior();
286  auto nb_bd = P.countBoundary();
287  CAPTURE( nb );
288  CAPTURE( nb_int );
289  CAPTURE( nb_bd );
290  std::vector<Point> Ppts;
291  P.getPoints( Ppts );
292  CAPTURE( Ppts );
293  REQUIRE( nb_bd == 3 );
294  REQUIRE( nb == nb_int + nb_bd );
295  }
296  }
297  GIVEN( "A closed segment S at (0,0,0), (8,-4,2)" ) {
298  Point a( 0, 0, 0 );
299  Point b( 8, -4, 2 );
300  Polytope P { a, b };
301  THEN( "Its interior is empty #Int(P) == 0" ) {
302  auto nb_int = P.countInterior();
303  REQUIRE( nb_int == 0 );
304  }
305  THEN( "It satisfies #In(P) == #Int(P) + #Bd(P) == #Bd(P) == 3" ) {
306  auto nb = P.count();
307  auto nb_int = P.countInterior();
308  auto nb_bd = P.countBoundary();
309  CAPTURE( nb );
310  CAPTURE( nb_int );
311  CAPTURE( nb_bd );
312  std::vector<Point> Ppts;
313  P.getPoints( Ppts );
314  CAPTURE( Ppts );
315  REQUIRE( nb_bd == 3 );
316  REQUIRE( nb == nb_int + nb_bd );
317  }
318  }
319 }
320 
321 // //
Aim: Represents an nD lattice polytope, i.e. a convex polyhedron bounded with vertices with integer c...
DigitalPlane::Point Vector
DGtal is the top-level namespace which contains all DGtal functions and types.
SCENARIO("BoundedLatticePolytope< Z2 > unit tests", "[lattice_polytope][2d]")
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
REQUIRE(domain.isInside(aPoint))