DGtal 1.3.0
Loading...
Searching...
No Matches
testDigitalConvexity.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/topology/KhalimskySpaceND.h"
37#include "DGtal/geometry/volumes/CellGeometry.h"
38#include "DGtal/geometry/volumes/DigitalConvexity.h"
39#include "DGtalCatch.h"
41
42using namespace std;
43using namespace DGtal;
44
45
47// Functions for testing class DigitalConvexity.
49
50SCENARIO( "DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]" )
51{
52 typedef KhalimskySpaceND<2,int> KSpace;
53 typedef KSpace::Point Point;
54 typedef DigitalConvexity< KSpace > DConvexity;
55
56 DConvexity dconv( Point( -5, -5 ), Point( 10, 10 ) );
57
58 GIVEN( "Given a fat simplex { (0,0), (4,-1), (2,5) } " ) {
59 std::vector<Point> V = { Point(0,0), Point(4,-1), Point(2,5) };
60 auto vertex_cover = dconv.makeCellCover( V.begin(), V.end() );
61 auto fat_simplex = dconv.makeSimplex ( V.begin(), V.end() );
62 auto inside_pts = dconv.insidePoints ( fat_simplex );
63 auto simplex_cover = dconv.makeCellCover( fat_simplex );
64 auto point_cover = dconv.makeCellCover( inside_pts.begin(), inside_pts.end() );
65 THEN( "The fat simplex is not degenerated." ) {
66 REQUIRE( dconv.isSimplexFullDimensional( V.begin(), V.end() ) );
67 }
68 THEN( "Its vertex cover contains 3 0-cells, 12 1-cells, 12 2-cells" ) {
69 REQUIRE( vertex_cover.computeNbCells( 0 ) == 3 );
70 REQUIRE( vertex_cover.computeNbCells( 1 ) == 12 );
71 REQUIRE( vertex_cover.computeNbCells( 2 ) == 12 );
72 }
73 THEN( "Its vertex cover is a subset of its point cover" ) {
74 REQUIRE( vertex_cover.subset( point_cover ) );
75 }
76 THEN( "Its point cover is a subset of its simplex cover" ) {
77 REQUIRE( point_cover.subset( simplex_cover ) );
78 }
79 THEN( "Being fat, its simplex cover is equal to its point cover" ) {
80 REQUIRE( simplex_cover.subset( point_cover ) );
81 }
82 }
83 GIVEN( "Given a thin simplex { (0,0), (4,3), (7,5) } " ) {
84 std::vector<Point> V = { Point(0,0), Point(4,3), Point(7,5) };
85 auto vertex_cover = dconv.makeCellCover( V.begin(), V.end() );
86 auto thin_simplex = dconv.makeSimplex ( V.begin(), V.end() );
87 auto inside_pts = dconv.insidePoints ( thin_simplex );
88 auto simplex_cover = dconv.makeCellCover( thin_simplex );
89 auto point_cover = dconv.makeCellCover( inside_pts.begin(), inside_pts.end() );
90 THEN( "The thin simplex is not degenerated." ) {
91 REQUIRE( dconv.isSimplexFullDimensional( V.begin(), V.end() ) );
92 }
93 THEN( "Its vertex cover contains 3 0-cells, 12 1-cells, 12 2-cells" ) {
94 REQUIRE( vertex_cover.computeNbCells( 0 ) == 3 );
95 REQUIRE( vertex_cover.computeNbCells( 1 ) == 12 );
96 REQUIRE( vertex_cover.computeNbCells( 2 ) == 12 );
97 }
98 THEN( "Its vertex cover is a subset of its point cover" ) {
99 REQUIRE( vertex_cover.subset( point_cover ) );
100 }
101 THEN( "Its point cover is a subset of its simplex cover" ) {
102 REQUIRE( point_cover.subset( simplex_cover ) );
103 }
104 THEN( "Being thin, its simplex cover is not equal to its point cover for 1<=dim<=2" ) {
105 REQUIRE( ! simplex_cover.subset( point_cover ) );
106 REQUIRE( simplex_cover.subset( point_cover, 0 ) );
107 REQUIRE( ! simplex_cover.subset( point_cover, 1 ) );
108 REQUIRE( ! simplex_cover.subset( point_cover, 2 ) );
109 }
110 }
111}
112
113SCENARIO( "DigitalConvexity< Z2 > fully convex triangles", "[convex_simplices][2d]" )
114{
115 typedef KhalimskySpaceND<2,int> KSpace;
116 typedef KSpace::Point Point;
117 typedef KSpace::Space Space;
118 typedef HyperRectDomain< Space > Domain;
119 typedef DigitalConvexity< KSpace > DConvexity;
120
121 Domain domain( Point( 0, 0 ), Point( 4, 4 ) );
122 DConvexity dconv( Point( -1, -1 ), Point( 5, 5 ) );
123
124 WHEN( "Computing all triangles in domain (0,0)-(4,4)." ) {
125 unsigned int nb_notsimplex = 0;
126 unsigned int nb_invalid = 0;
127 unsigned int nb_degenerated= 0;
128 unsigned int nb_common = 0;
129 unsigned int nb_unitary = 0;
130 for ( auto a : domain )
131 for ( auto b : domain )
132 for ( auto c : domain )
133 {
134 nb_notsimplex += ! dconv.isSimplexFullDimensional( { a, b, c } ) ? 1 : 0;
135 auto tri_type = dconv.simplexType( { a, b, c } );
136 nb_degenerated += tri_type == DConvexity::SimplexType::DEGENERATED ? 1 : 0;
137 nb_invalid += tri_type == DConvexity::SimplexType::INVALID ? 1 : 0;
138 nb_unitary += tri_type == DConvexity::SimplexType::UNITARY ? 1 : 0;
139 nb_common += tri_type == DConvexity::SimplexType::COMMON ? 1 : 0;
140 }
141 THEN( "All 2737 invalid triangles are degenerated " ) {
142 REQUIRE( nb_invalid == 0 );
143 REQUIRE( nb_notsimplex == nb_degenerated );
144 REQUIRE( nb_degenerated == 2737 );
145 }
146 THEN( "There are 12888 valid triangles" ) {
147 REQUIRE( nb_unitary + nb_common == 12888 );
148 }
149 THEN( "There are fewer (1920) unitary triangles than common triangles (10968)" ) {
150 REQUIRE( nb_unitary == 1920 );
151 REQUIRE( nb_common == 10968 );
152 REQUIRE( nb_unitary < nb_common );
153 }
154 THEN( "The total number of triangles (unitary, common, degenerated) is (domain size)^3, i.e. 5^6" ) {
155 REQUIRE( nb_unitary + nb_common + nb_degenerated == 5*5*5*5*5*5 );
156 }
157 }
158 WHEN( "Computing all triangles in domain (0,0)-(4,4)." ) {
159 unsigned int nbsimplex= 0;
160 unsigned int nb0 = 0;
161 unsigned int nb1 = 0;
162 unsigned int nb2 = 0;
163 unsigned int nb01_not2 = 0;
164 for ( auto a : domain )
165 for ( auto b : domain )
166 for ( auto c : domain )
167 {
168 if ( ! ( ( a < b ) && ( b < c ) ) ) continue;
169 if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
170 auto triangle = dconv.makeSimplex( { a, b, c } );
171 bool cvx0 = dconv.isKConvex( triangle, 0 );
172 bool cvx1 = dconv.isKConvex( triangle, 1 );
173 bool cvx2 = dconv.isKConvex( triangle, 2 );
174 nbsimplex += 1;
175 nb0 += cvx0 ? 1 : 0;
176 nb1 += cvx1 ? 1 : 0;
177 nb2 += cvx2 ? 1 : 0;
178 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
179 }
180 THEN( "All valid triangles are 0-convex." ) {
181 REQUIRE( nb0 == nbsimplex );
182 }
183 THEN( "There are less 1-convex and 2-convex than 0-convex." ) {
184 REQUIRE( nb1 < nb0 );
185 REQUIRE( nb2 < nb0 );
186 }
187 THEN( "When the triangle is 0-convex and 1-convex, then it is 2-convex." ) {
188 REQUIRE( nb1 <= nb2 );
189 REQUIRE( nb01_not2 == 0 );
190 }
191 }
192}
193SCENARIO( "DigitalConvexity< Z3 > fully convex tetrahedra", "[convex_simplices][3d]" )
194{
195 typedef KhalimskySpaceND<3,int> KSpace;
196 typedef KSpace::Point Point;
197 typedef KSpace::Space Space;
198 typedef HyperRectDomain< Space > Domain;
199 typedef DigitalConvexity< KSpace > DConvexity;
200
201 Domain domain( Point( 0, 0, 0 ), Point( 3, 3, 3 ) );
202 DConvexity dconv( Point( -1, -1, -1 ), Point( 4, 4, 4 ) );
203
204 WHEN( "Computing all lexicographically ordered tetrahedra anchored at (0,0,0) in domain (0,0,0)-(3,3,3)." ) {
205 unsigned int nb_notsimplex = 0;
206 unsigned int nb_invalid = 0;
207 unsigned int nb_degenerated= 0;
208 unsigned int nb_common = 0;
209 unsigned int nb_unitary = 0;
210 Point a(0,0,0);
211 for ( auto b : domain )
212 for ( auto c : domain )
213 for ( auto d : domain )
214 {
215 if ( ! ( ( a < b ) && ( b < c ) && ( c < d ) ) ) continue;
216 nb_notsimplex += ! dconv.isSimplexFullDimensional( {a,b,c,d} ) ? 1 : 0;
217 auto tri_type = dconv.simplexType( { a, b, c, d } );
218 nb_degenerated += tri_type == DConvexity::SimplexType::DEGENERATED ? 1 : 0;
219 nb_invalid += tri_type == DConvexity::SimplexType::INVALID ? 1 : 0;
220 nb_unitary += tri_type == DConvexity::SimplexType::UNITARY ? 1 : 0;
221 nb_common += tri_type == DConvexity::SimplexType::COMMON ? 1 : 0;
222 }
223 THEN( "All 4228 invalid tetrahedra are degenerated " ) {
224 REQUIRE( nb_invalid == 0 );
225 REQUIRE( nb_notsimplex == nb_degenerated );
226 REQUIRE( nb_degenerated == 4228 );
227 }
228 THEN( "There are 35483 valid tetrahedra" ) {
229 REQUIRE( nb_unitary + nb_common == 35483 );
230 }
231 THEN( "There are fewer (2515) unitary triangles than common triangles (32968)" ) {
232 REQUIRE( nb_unitary == 2515 );
233 REQUIRE( nb_common == 32968 );
234 REQUIRE( nb_unitary < nb_common );
235 }
236 THEN( "The total number of triangles (unitary, common, degenerated) is 39711" ) {
237 REQUIRE( nb_unitary + nb_common + nb_degenerated == 39711 );
238 }
239 }
240 WHEN( "Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
241 const unsigned int nb = 50;
242 unsigned int nbsimplex= 0;
243 unsigned int nb0 = 0;
244 unsigned int nb1 = 0;
245 unsigned int nb2 = 0;
246 unsigned int nb3 = 0;
247 unsigned int nb012_not3 = 0;
248 unsigned int nbf = 0;
249 unsigned int nbfg = 0;
250 unsigned int nbffast = 0;
251 unsigned int nb0123 = 0;
252 for ( unsigned int i = 0; i < nb; ++i )
253 {
254 Point a( rand() % 5, rand() % 5, rand() % 5 );
255 Point b( rand() % 5, rand() % 5, rand() % 5 );
256 Point c( rand() % 5, rand() % 5, rand() % 5 );
257 Point d( rand() % 5, rand() % 5, rand() % 5 );
258 if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
259 auto tetra = dconv.makeSimplex( { a, b, c, d } );
260 std::vector< Point > X;
261 tetra.getPoints( X );
262 bool cvx0 = dconv.isKConvex( tetra, 0 );
263 bool cvx1 = dconv.isKConvex( tetra, 1 );
264 bool cvx2 = dconv.isKConvex( tetra, 2 );
265 bool cvx3 = dconv.isKConvex( tetra, 3 );
266 bool cvxf = dconv.isFullyConvex( tetra );
267 bool cvxfg = dconv.isFullyConvex( X, false );
268 bool cvxffast = dconv.isFullyConvexFast( X );
269 if ( cvxf != cvxfg || cvxf != cvxffast) {
270 std::cout << "[" << cvx0 << cvx1 << cvx2 << cvx3 << "] "
271 << "[" << cvxf << "] [" << cvxfg
272 << "] [" << cvxffast << "]"
273 << a << b << c << d << std::endl;
274 }
275 nbsimplex += 1;
276 nb0 += cvx0 ? 1 : 0;
277 nb1 += cvx1 ? 1 : 0;
278 nb2 += cvx2 ? 1 : 0;
279 nb3 += cvx3 ? 1 : 0;
280 nbf += cvxf ? 1 : 0;
281 nbfg += cvxfg ? 1 : 0;
282 nbffast += cvxffast ? 1 : 0;
283 nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
284 nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
285 }
286 THEN( "All valid tetrahedra are 0-convex." ) {
287 REQUIRE( nb0 == nbsimplex );
288 }
289 THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
290 REQUIRE( nb1 < nb0 );
291 REQUIRE( nb2 < nb0 );
292 REQUIRE( nb3 < nb0 );
293 }
294 THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
295 REQUIRE( nb1 <= nb3 );
296 REQUIRE( nb2 <= nb3 );
297 REQUIRE( nb012_not3 == 0 );
298 REQUIRE( nbf == nb0123 );
299 }
300 THEN( "All methods for computing full convexity agree." ) {
301 REQUIRE( nbf == nbfg );
302 REQUIRE( nbf == nbffast );
303 }
304 }
305}
306
307SCENARIO( "DigitalConvexity< Z3 > rational fully convex tetrahedra", "[convex_simplices][3d][rational]" )
308{
309 typedef KhalimskySpaceND<3,int> KSpace;
310 typedef KSpace::Point Point;
311 typedef KSpace::Space Space;
312 typedef DigitalConvexity< KSpace > DConvexity;
313
314 DConvexity dconv( Point( -1, -1, -1 ), Point( 10, 10, 10 ) );
315 WHEN( "Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
316 const unsigned int nb = 50;
317 unsigned int nbsimplex= 0;
318 unsigned int nb0 = 0;
319 unsigned int nb1 = 0;
320 unsigned int nb2 = 0;
321 unsigned int nb3 = 0;
322 unsigned int nb012_not3 = 0;
323 unsigned int nbf = 0;
324 unsigned int nb0123 = 0;
325 for ( unsigned int i = 0; i < nb; ++i )
326 {
327 Point a( 2*(rand() % 10), rand() % 20, 2*(rand() % 10) );
328 Point b( rand() % 20, 2*(rand() % 10), 2*(rand() % 10) );
329 Point c( 2*(rand() % 10), 2*(rand() % 10), rand() % 20 );
330 Point d( 2*(rand() % 10), 2*(rand() % 10), 2*(rand() % 10) );
331 if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
332 auto tetra = dconv.makeRationalSimplex( { Point(2,2,2), a, b, c, d } );
333 bool cvx0 = dconv.isKConvex( tetra, 0 );
334 bool cvx1 = dconv.isKConvex( tetra, 1 );
335 bool cvx2 = dconv.isKConvex( tetra, 2 );
336 bool cvx3 = dconv.isKConvex( tetra, 3 );
337 bool cvxf = dconv.isFullyConvex( tetra );
338 nbsimplex += 1;
339 nb0 += cvx0 ? 1 : 0;
340 nb1 += cvx1 ? 1 : 0;
341 nb2 += cvx2 ? 1 : 0;
342 nb3 += cvx3 ? 1 : 0;
343 nbf += cvxf ? 1 : 0;
344 nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
345 nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
346 }
347 THEN( "All valid tetrahedra are 0-convex." ) {
348 REQUIRE( nb0 == nbsimplex );
349 }
350 THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
351 REQUIRE( nb1 < nb0 );
352 REQUIRE( nb2 < nb0 );
353 REQUIRE( nb3 < nb0 );
354 }
355 THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
356 CAPTURE( nb0 );
357 CAPTURE( nb1 );
358 CAPTURE( nb2 );
359 CAPTURE( nb3 );
360 CAPTURE( nb012_not3 );
361 CAPTURE( nbf );
362 REQUIRE( nb2 <= nb3 );
363 REQUIRE( nb012_not3 == 0 );
364 REQUIRE( nbf == nb0123 );
365 }
366 }
367}
368
369
370SCENARIO( "DigitalConvexity< Z2 > rational fully convex tetrahedra", "[convex_simplices][2d][rational]" )
371{
372 typedef KhalimskySpaceND<2,int> KSpace;
373 typedef KSpace::Point Point;
374 typedef KSpace::Space Space;
375 typedef DigitalConvexity< KSpace > DConvexity;
376
377 DConvexity dconv( Point( -1, -1 ), Point( 10, 10 ) );
378 WHEN( "Computing many triangle in domain (0,0)-(9,9)." ) {
379 const unsigned int nb = 50;
380 unsigned int nbsimplex= 0;
381 unsigned int nb0 = 0;
382 unsigned int nb1 = 0;
383 unsigned int nb2 = 0;
384 unsigned int nb01_not2 = 0;
385 unsigned int nbf = 0;
386 unsigned int nb012 = 0;
387 for ( unsigned int i = 0; i < nb; ++i )
388 {
389 Point a( 2*(rand() % 10), rand() % 20 );
390 Point b( rand() % 20 , 2*(rand() % 10) );
391 Point c( 2*(rand() % 10), 2*(rand() % 10) );
392 if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
393 auto tetra = dconv.makeRationalSimplex( { Point(2,2), a, b, c } );
394 bool cvx0 = dconv.isKConvex( tetra, 0 );
395 bool cvx1 = dconv.isKConvex( tetra, 1 );
396 bool cvx2 = dconv.isKConvex( tetra, 2 );
397 bool cvxf = dconv.isFullyConvex( tetra );
398 nbsimplex += 1;
399 nb0 += cvx0 ? 1 : 0;
400 nb1 += cvx1 ? 1 : 0;
401 nb2 += cvx2 ? 1 : 0;
402 nbf += cvxf ? 1 : 0;
403 nb012 += ( cvx0 && cvx1 && cvx2 ) ? 1 : 0;
404 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
405 }
406 THEN( "All valid tetrahedra are 0-convex." ) {
407 REQUIRE( nb0 == nbsimplex );
408 }
409 THEN( "There are less 1-convex, 2-convex than 0-convex." ) {
410 REQUIRE( nb1 < nb0 );
411 REQUIRE( nb2 < nb0 );
412 }
413 THEN( "When the tetrahedron is 0-convex, and 1-convex, then it is 2-convex." ) {
414 CAPTURE( nb0 );
415 CAPTURE( nb1 );
416 CAPTURE( nb2 );
417 CAPTURE( nb01_not2 );
418 CAPTURE( nbf );
419 REQUIRE( nb1 <= nb2 );
420 REQUIRE( nb01_not2 == 0 );
421 REQUIRE( nbf == nb012 );
422 }
423 }
424}
425
426SCENARIO( "DigitalConvexity< Z3 > full subconvexity of segments and triangles", "[subconvexity][3d]" )
427{
428 typedef KhalimskySpaceND<3,int> KSpace;
429 typedef KSpace::Point Point;
430 typedef KSpace::Space Space;
431 typedef HyperRectDomain< Space > Domain;
432 typedef DigitalConvexity< KSpace > DConvexity;
433
434 Domain domain( Point( -5, -5, -5 ), Point( 5, 5, 5 ) );
435 DConvexity dconv( Point( -6, -6, -6 ), Point( 6, 6, 6 ) );
436
437 WHEN( "Computing many tetrahedra" ) {
438 const unsigned int nb = 50;
439 unsigned int nb_fulldim = 0;
440 unsigned int nb_ok_seg = 0;
441 unsigned int nb_ok_tri = 0;
442 for ( unsigned int l = 0; l < nb; ++l )
443 {
444 const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
445 const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
446 const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
447 const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
448 const std::vector<Point> pts = { a, b, c, d };
449 const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
450 nb_fulldim += fulldim ? 1 : 0;
451 if ( ! fulldim ) continue;
452 auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
453 auto cover = dconv.makeCellCover( simplex, 0, 3 );
454 {
455 unsigned int nb_subconvex = 0;
456 unsigned int nb_total = 0;
457 for ( unsigned int i = 0; i < 4; i++ )
458 for ( unsigned int j = i+1; j < 4; j++ )
459 {
460 auto segment = dconv.makeSimplex( { pts[ i ], pts[ j ] } );
461 bool ok = dconv.isFullySubconvex( segment, cover );
462 nb_subconvex += ok ? 1 : 0;
463 nb_total += 1;
464 if ( ! ok ) {
465 trace.info() << "****** SEGMENT NOT SUBCONVEX ******" << std::endl;
466 trace.info() << "splx v =" << a << b << c << d << std::endl;
467 trace.info() << "simplex=" << simplex << std::endl;
468 trace.info() << "seg v =" << pts[ i ] << pts[ j ] << std::endl;
469 trace.info() << "segment=" << segment << std::endl;
470 }
471 }
472 nb_ok_seg += ( nb_subconvex == nb_total ) ? 1 : 0;
473 }
474 {
475 unsigned int nb_subconvex = 0;
476 unsigned int nb_total = 0;
477 for ( unsigned int i = 0; i < 4; i++ )
478 for ( unsigned int j = i+1; j < 4; j++ )
479 for ( unsigned int k = j+1; k < 4; k++ )
480 {
481 auto triangle = dconv.makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
482 bool ok = dconv.isFullySubconvex( triangle, cover );
483 nb_subconvex += ok ? 1 : 0;
484 nb_total += 1;
485 if ( ! ok ) {
486 trace.info() << "****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
487 trace.info() << "splx v =" << a << b << c << d << std::endl;
488 trace.info() << "simplex=" << simplex << std::endl;
489 trace.info() << "tri v =" << pts[ i ] << pts[ j ]
490 << pts[ k ] << std::endl;
491 trace.info() << "triangle=" << triangle << std::endl;
492 }
493 }
494 nb_ok_tri += ( nb_subconvex == nb_total ) ? 1 : 0;
495 }
496 }
497 THEN( "At least half the tetrahedra are full dimensional." ) {
498 REQUIRE( nb_fulldim >= nb / 2 );
499 }
500 THEN( "All segments of a tetrahedron should be subconvex to it." ) {
501 REQUIRE( nb_ok_seg == nb_fulldim );
502 }
503 THEN( "All triangles of a tetrahedron should be subconvex to it." ) {
504 REQUIRE( nb_ok_tri == nb_fulldim );
505 }
506 }
507}
508
509SCENARIO( "DigitalConvexity< Z3 > full convexity of polyhedra", "[full_convexity][3d]" )
510{
511 typedef KhalimskySpaceND<3,int> KSpace;
512 typedef KSpace::Point Point;
513 typedef KSpace::Space Space;
514 typedef HyperRectDomain< Space > Domain;
515 typedef DigitalConvexity< KSpace > DConvexity;
516
517 Domain domain( Point( -35, -35, -35 ), Point( 35, 35, 35 ) );
518 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
519
520 const unsigned int nb = 10;
521 unsigned int nbfg = 0;
522 unsigned int nbffast = 0;
523 unsigned int nbfenv = 0;
524 typedef std::vector< Point > PointRange;
525 std::vector< PointRange > XX;
526 for ( unsigned int i = 0; i < nb; ++i )
527 {
528 unsigned int k = 100;
529 PointRange X( k );
530 for ( unsigned int j = 0; j < k; ++ j )
531 X[ j ] = Point( rand() % 10, rand() % 10, rand() % 10 );
532 auto P = dconv.makePolytope( X );
533 PointRange Y;
534 P.getPoints( Y );
535 XX.push_back( Y );
536 }
537 Clock c;
538 c.startClock();
539 for ( const auto& X : XX )
540 {
541 bool fcvx = dconv.isFullyConvex( X, false );
542 nbfg += fcvx ? 1 : 0;
543 }
544 double t1 = c.stopClock();
545 c.startClock();
546 for ( const auto& X : XX )
547 {
548 bool fcvx = dconv.isFullyConvexFast( X );
549 nbffast += fcvx ? 1 : 0;
550 }
551 double t2 = c.stopClock();
552 c.startClock();
553 for ( const auto& X : XX )
554 {
555 auto card = dconv.envelope( X ).size();
556 bool fcvx = card == X.size();
557 nbfenv += fcvx ? 1 : 0;
558 }
559 double t3 = c.stopClock();
560 WHEN( "Computing many polytopes." ) {
561 THEN( "All three methods agree on full convexity results" ) {
562 CAPTURE( t1 );
563 CAPTURE( t2 );
564 CAPTURE( t3 );
565 REQUIRE( nbfg == nbffast );
566 REQUIRE( nbfg == nbfenv );
567 }
568 }
569}
570
571
572SCENARIO( "DigitalConvexity< Z4 > full convexity of polyhedra", "[full_convexity][4d]" )
573{
574 typedef KhalimskySpaceND<4,int> KSpace;
575 typedef KSpace::Point Point;
576 typedef KSpace::Space Space;
577 typedef HyperRectDomain< Space > Domain;
578 typedef DigitalConvexity< KSpace > DConvexity;
579
580 Domain domain( Point( -35, -35, -35, -35 ), Point( 35, 35, 35, 35 ) );
581 DConvexity dconv( Point( -36, -36, -36, -36 ), Point( 36, 36, 36, 36 ) );
582
583 const unsigned int nb = 4;
584 unsigned int nbfg = 0;
585 unsigned int nbffast = 0;
586 unsigned int nbfenv = 0;
587 typedef std::vector< Point > PointRange;
588 std::vector< PointRange > XX;
589 for ( unsigned int i = 0; i < nb; ++i )
590 {
591 unsigned int k = 100;
592 PointRange X( k );
593 for ( unsigned int j = 0; j < k; ++ j )
594 X[ j ] = Point( rand() % 8, rand() % 8, rand() % 8, rand() % 8 );
595 auto P = dconv.makePolytope( X );
596 PointRange Y;
597 P.getPoints( Y );
598 XX.push_back( Y );
599 }
600 Clock c;
601 c.startClock();
602 for ( const auto& X : XX )
603 {
604 bool fcvx = dconv.isFullyConvex( X, false );
605 nbfg += fcvx ? 1 : 0;
606 }
607 double t1 = c.stopClock();
608 c.startClock();
609 for ( const auto& X : XX )
610 {
611 bool fcvx = dconv.isFullyConvexFast( X );
612 nbffast += fcvx ? 1 : 0;
613 }
614 double t2 = c.stopClock();
615 c.startClock();
616 for ( const auto& X : XX )
617 {
618 auto card = dconv.envelope( X ).size();
619 bool fcvx = card == X.size();
620 nbfenv += fcvx ? 1 : 0;
621 }
622 double t3 = c.stopClock();
623 WHEN( "Computing many polytopes." ) {
624 THEN( "All three methods agree on full convexity results" ) {
625 CAPTURE( t1 );
626 CAPTURE( t2 );
627 CAPTURE( t3 );
628 REQUIRE( nbfg == nbffast );
629 REQUIRE( nbfg == nbfenv );
630 }
631 }
632}
633
634SCENARIO( "DigitalConvexity< Z2 > sub-convexity of polyhedra", "[full_subconvexity][2d]" )
635{
636 typedef KhalimskySpaceND<2,int> KSpace;
637 typedef KSpace::Point Point;
638 typedef KSpace::Space Space;
639 typedef HyperRectDomain< Space > Domain;
640 typedef DigitalConvexity< KSpace > DConvexity;
641
642 DConvexity dconv( Point( -36, -36 ), Point( 36, 36 ) );
643 unsigned int k = 6;
644 std::vector< Point > X( k );
645 X[ 0 ] = Point( 0,0 );
646 X[ 1 ] = Point( 7,-2 );
647 X[ 2 ] = Point( 3,6 );
648 X[ 3 ] = Point( 5, 5 );
649 X[ 4 ] = Point( 2, 3 );
650 X[ 5 ] = Point( -1, 1 );
651 auto P = dconv.makePolytope( X, true );
652 auto CG = dconv.makeCellCover( P, 0, 2 );
653 auto L = dconv.StarCvxH( X, 0 );
654 REQUIRE( CG.nbCells() == L.size() );
655 for ( int i = 0; i < k; i++ )
656 for ( int j = i+1; j < k; j++ )
657 {
658 std::vector< Point > Z { X[ i ], X[ j ] };
659 const auto Q = dconv.makePolytope( Z );
660 bool tangent_old = dconv.isFullySubconvex( Q, CG );
661 bool tangent_new = dconv.isFullySubconvex( Z, L );
662 bool tangent_ab_old = dconv.isFullySubconvex( X[ i ], X[ j ], CG );
663 bool tangent_ab_new = dconv.isFullySubconvex( X[ i ], X[ j ], L );
664 REQUIRE( tangent_old == tangent_new );
665 REQUIRE( tangent_ab_old == tangent_ab_new );
666 REQUIRE( tangent_new == tangent_ab_new );
667 }
668}
669
670SCENARIO( "DigitalConvexity< Z3 > sub-convexity of polyhedra", "[full_subconvexity][3d]" )
671{
672 typedef KhalimskySpaceND<3,int> KSpace;
673 typedef KSpace::Point Point;
674 typedef KSpace::Space Space;
675 typedef HyperRectDomain< Space > Domain;
676 typedef DigitalConvexity< KSpace > DConvexity;
677
678 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
679 std::vector< Point > X( 5 );
680 X[ 0 ] = Point( 0,0,0 );
681 X[ 1 ] = Point( 0,5,1 );
682 X[ 2 ] = Point( 2,1,6 );
683 X[ 3 ] = Point( 6,1,1 );
684 X[ 4 ] = Point( -2,-2,-3 );
685 auto P = dconv.makePolytope( X, true );
686 auto CG = dconv.makeCellCover( P, 0, 3 );
687 auto L = dconv.StarCvxH( X, 0 );
688 std::vector< Point > Y;
689 P.getPoints( Y );
690 REQUIRE( CG.nbCells() == L.size() );
691 unsigned int nb = 0;
692 unsigned int nb_ok = 0;
693 unsigned int nb_tgt= 0;
694 for ( int i = 0; i < 100; i++ )
695 {
696 Point a( rand() % 6, rand() % 6, rand() % 6 );
697 Point b( rand() % 6, rand() % 6, rand() % 6 );
698 // Point b( rand() % 20 - 10, rand() % 20 - 10, rand() % 20 - 10 );
699 bool tangent_ab_old = dconv.isFullySubconvex( a, b, CG );
700 bool tangent_ab_new = dconv.isFullySubconvex( a, b, L );
701 nb_tgt += tangent_ab_new ? 1 : 0;
702 nb_ok += ( tangent_ab_old == tangent_ab_new ) ? 1 : 0;
703 nb += 1;
704 }
705 REQUIRE( nb == nb_ok );
706 REQUIRE( 0 < nb_tgt );
707 REQUIRE( nb_tgt < 100 );
708}
709
710SCENARIO( "DigitalConvexity< Z3 > envelope", "[envelope][3d]" )
711{
712 typedef KhalimskySpaceND<3,int> KSpace;
713 typedef KSpace::Point Point;
714 typedef KSpace::Space Space;
715 typedef HyperRectDomain< Space > Domain;
716 typedef DigitalConvexity< KSpace > DConvexity;
717
718 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
719
720 WHEN( "Computing the envelope Z of a digital set X with direct algorithm" ) {
721 THEN( "Z contains X" ){
722 for ( int k = 0; k < 5; k++ )
723 {
724 int n = 3 + ( rand() % 7 );
725 std::set< Point > S;
726 for ( int i = 0; i < n; i++ )
727 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
728 std::vector< Point > X( S.cbegin(), S.cend() );
729 CAPTURE( X );
730 auto Z = dconv.envelope( X, DConvexity::EnvelopeAlgorithm::DIRECT );
731 CAPTURE( Z );
732 CAPTURE( dconv.depthLastEnvelope() );
733 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
734 X.cbegin(), X.cend() );
735 REQUIRE( X.size() <= Z.size() );
736 REQUIRE( Z_includes_X );
737 }
738 THEN( "Z is fully convex" ){
739 for ( int k = 0; k < 5; k++ )
740 {
741 int n = 3 + ( rand() % 7 );
742 std::set< Point > S;
743 for ( int i = 0; i < n; i++ )
744 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
745 std::vector< Point > X( S.cbegin(), S.cend() );
746 auto Z = dconv.envelope( X );
747 CAPTURE( dconv.depthLastEnvelope() );
748 REQUIRE( dconv.isFullyConvex( Z ) );
749 }
750 }
751 }
752 }
753 WHEN( "Computing the envelope Z of a digital set X with LatticeSet algorithm" ) {
754 THEN( "Z contains X" ){
755 for ( int k = 0; k < 5; k++ )
756 {
757 int n = 3 + ( rand() % 7 );
758 std::set< Point > S;
759 for ( int i = 0; i < n; i++ )
760 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
761 std::vector< Point > X( S.cbegin(), S.cend() );
762 auto Z = dconv.envelope( X, DConvexity::EnvelopeAlgorithm::LATTICE_SET );
763 CAPTURE( dconv.depthLastEnvelope() );
764 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
765 X.cbegin(), X.cend() );
766 REQUIRE( X.size() <= Z.size() );
767 REQUIRE( Z_includes_X );
768 }
769 THEN( "Z is fully convex" ){
770 for ( int k = 0; k < 5; k++ )
771 {
772 int n = 3 + ( rand() % 7 );
773 std::set< Point > S;
774 for ( int i = 0; i < n; i++ )
775 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
776 std::vector< Point > X( S.cbegin(), S.cend() );
777 auto Z = dconv.envelope( X );
778 CAPTURE( dconv.depthLastEnvelope() );
779 REQUIRE( dconv.isFullyConvex( Z ) );
780 }
781 }
782 }
783 }
784}
785
786SCENARIO( "DigitalConvexity< Z2 > envelope", "[envelope][2d]" )
787{
788 typedef KhalimskySpaceND<2,int> KSpace;
789 typedef KSpace::Point Point;
790 typedef KSpace::Space Space;
791 typedef HyperRectDomain< Space > Domain;
792 typedef DigitalConvexity< KSpace > DConvexity;
793
794 DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
795
796 WHEN( "Computing the envelope Z of two points" ) {
797 THEN( "it requires at most one iteration" ){
798 for ( int k = 0; k < 10; k++ )
799 {
800 std::vector< Point > X;
801 X.push_back( Point( rand() % 100, rand() % 100 ) );
802 X.push_back( Point( rand() % 100, rand() % 100 ) );
803 std::sort( X.begin(), X.end() );
804 auto Z = dconv.envelope( X );
805 REQUIRE( dconv.depthLastEnvelope() <= 1 );
806 }
807 }
808 }
809}
810
811SCENARIO( "DigitalConvexity< Z2 > relative envelope", "[rel_envelope][2d]" )
812{
813 typedef KhalimskySpaceND<2,int> KSpace;
814 typedef KSpace::Point Point;
815 typedef KSpace::Space Space;
816 typedef HyperRectDomain< Space > Domain;
817 typedef DigitalConvexity< KSpace > DConvexity;
818
819 DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
820
821 std::vector< Point > X { Point( -10, -7 ), Point( 10, 7 ) };
822 std::vector< Point > Y { Point( -11, -6 ), Point( 9, 8 ) };
823 std::sort( X.begin(), X.end() );
824 std::sort( Y.begin(), Y.end() );
825 X = dconv.envelope( X );
826 Y = dconv.envelope( Y );
827 REQUIRE( dconv.isFullyConvex( X ) );
828 REQUIRE( dconv.isFullyConvex( Y ) );
829 WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
830 auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
831 auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
832 THEN( "Both sets are fully convex" ){
833 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
834 REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
835 }
836 THEN( "There are inclusion rules between sets" ){
837 CAPTURE( FC_X_rel_Y );
838 CAPTURE( FC_Y_rel_X );
839 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
840 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
841 REQUIRE( std::includes( X.cbegin(), X.cend(),
842 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
843 }
844 }
845 WHEN( "Computing the envelope of X relative to Y specified by a predicate" ) {
846 auto PredY = [] ( Point p )
847 { return ( -4 <= p.dot( Point( 2,5 ) ) ) && ( p.dot( Point( 2,5 ) ) < 9 ); };
848 auto FC_X_rel_Y = dconv.relativeEnvelope( X, PredY );
849 THEN( "It is fully convex and included in Y" ){
850 CAPTURE( FC_X_rel_Y );
851 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
852 int nb = 0;
853 int nb_in = 0;
854 for ( auto p : FC_X_rel_Y )
855 {
856 nb_in += PredY( p ) ? 1 : 0;
857 nb += 1;
858 }
859 REQUIRE( nb == nb_in );
860 }
861 }
862}
863
864SCENARIO( "DigitalConvexity< Z3 > relative envelope", "[rel_envelope][3d]" )
865{
866 typedef KhalimskySpaceND<3,int> KSpace;
867 typedef KSpace::Point Point;
868 typedef KSpace::Space Space;
869 typedef HyperRectDomain< Space > Domain;
870 typedef DigitalConvexity< KSpace > DConvexity;
871
872 DConvexity dconv( Point( -360, -360, -360 ), Point( 360, 360, 360 ) );
873
874 std::vector< Point > X { Point( -61, -20, -8 ), Point( 43, 25, 9 ) };
875 std::vector< Point > Y { Point( -50, -27, -10 ), Point( 40, 37, 17 ) };
876 std::sort( X.begin(), X.end() );
877 std::sort( Y.begin(), Y.end() );
878 X = dconv.envelope( X );
879 Y = dconv.envelope( Y );
880 REQUIRE( dconv.isFullyConvex( X ) );
881 REQUIRE( dconv.isFullyConvex( Y ) );
882 WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
883 auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
884 auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
885 THEN( "Both sets are fully convex" ){
886 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
887 REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
888 }
889 THEN( "There are inclusion rules between sets" ){
890 CAPTURE( FC_X_rel_Y );
891 CAPTURE( FC_Y_rel_X );
892 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
893 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
894 REQUIRE( std::includes( X.cbegin(), X.cend(),
895 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
896 }
897 }
898}
DGtal::BoundedLatticePolytope::getPoints
void getPoints(std::vector< Point > &pts) const
DGtal::Clock::startClock
void startClock()
DGtal::DigitalConvexity::insidePoints
static PointRange insidePoints(const LatticePolytope &polytope)
DGtal::DigitalConvexity::makeRationalSimplex
static RationalPolytope makeRationalSimplex(Integer d, PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::isKConvex
bool isKConvex(const LatticePolytope &P, const Dimension k) const
DGtal::DigitalConvexity::isSimplexFullDimensional
static bool isSimplexFullDimensional(PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::makeSimplex
static LatticePolytope makeSimplex(PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::isFullyConvexFast
bool isFullyConvexFast(const PointRange &X) const
DGtal::DigitalConvexity::isFullyConvex
bool isFullyConvex(const PointRange &X, bool convex0=false) const
DGtal::DigitalConvexity::envelope
PointRange envelope(const PointRange &Z, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
DGtal::DigitalConvexity::relativeEnvelope
PointRange relativeEnvelope(const PointRange &Z, const PointRange &Y, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
DGtal::DigitalConvexity::isFullySubconvex
bool isFullySubconvex(const PointRange &Y, const LatticeSet &StarX) const
DGtal::DigitalConvexity::StarCvxH
LatticeSet StarCvxH(const PointRange &X, Dimension axis=dimension) const
DGtal::DigitalConvexity::depthLastEnvelope
Size depthLastEnvelope() const
DGtal::DigitalConvexity::makePolytope
LatticePolytope makePolytope(const PointRange &X, bool make_minkowski_summable=false) const
DGtal::DigitalConvexity::makeCellCover
CellGeometry makeCellCover(PointIterator itB, PointIterator itE, Dimension i=0, Dimension k=KSpace::dimension) const
DGtal::DigitalConvexity::simplexType
static SimplexType simplexType(PointIterator itB, PointIterator itE)
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:394
DGtal::Trace::info
std::ostream & info()
PointRange
std::vector< Point > PointRange
Definition: digitalPolyhedronBuilder3D.cpp:84
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition: ClosedIntegerHalfPlane.h:49
DGtal::trace
Trace trace
Definition: Common.h:154
std
STL namespace.
Point
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE
CAPTURE(thicknessHV)
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition: testCubicalComplex.cpp:70
domain
Domain domain
Definition: testProjection.cpp:88
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
REQUIRE
REQUIRE(domain.isInside(aPoint))
SCENARIO
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")
Definition: testUnorderedSetByBlock.cpp:58
Generated on Sun Nov 27 2022 15:15:56 for DGtal by doxygen 1.9.5