DGtal 1.3.0
Loading...
Searching...
No Matches
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
40using namespace std;
41using namespace DGtal;
42
43
45// Functions for testing class BoundedLatticePolytope.
47
48SCENARIO( "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
148SCENARIO( "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...
DGtal is the top-level namespace which contains all DGtal functions and types.
STL namespace.
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
REQUIRE(domain.isInside(aPoint))
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")