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"
43 using namespace DGtal;
50 SCENARIO(
"DigitalConvexity< Z2 > unit tests",
"[digital_convexity][2d]" )
56 DConvexity dconv(
Point( -5, -5 ),
Point( 10, 10 ) );
58 GIVEN(
"Given a fat simplex { (0,0), (4,-1), (2,5) } " ) {
60 auto vertex_cover = dconv.
makeCellCover( V.begin(), V.end() );
61 auto fat_simplex = dconv.
makeSimplex ( V.begin(), V.end() );
64 auto point_cover = dconv.
makeCellCover( inside_pts.begin(), inside_pts.end() );
65 THEN(
"The fat simplex is not degenerated." ) {
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 );
73 THEN(
"Its vertex cover is a subset of its point cover" ) {
74 REQUIRE( vertex_cover.subset( point_cover ) );
76 THEN(
"Its point cover is a subset of its simplex cover" ) {
77 REQUIRE( point_cover.subset( simplex_cover ) );
79 THEN(
"Being fat, its simplex cover is equal to its point cover" ) {
80 REQUIRE( simplex_cover.subset( point_cover ) );
83 GIVEN(
"Given a thin simplex { (0,0), (4,3), (7,5) } " ) {
85 auto vertex_cover = dconv.
makeCellCover( V.begin(), V.end() );
86 auto thin_simplex = dconv.
makeSimplex ( V.begin(), V.end() );
89 auto point_cover = dconv.
makeCellCover( inside_pts.begin(), inside_pts.end() );
90 THEN(
"The thin simplex is not degenerated." ) {
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 );
98 THEN(
"Its vertex cover is a subset of its point cover" ) {
99 REQUIRE( vertex_cover.subset( point_cover ) );
101 THEN(
"Its point cover is a subset of its simplex cover" ) {
102 REQUIRE( point_cover.subset( simplex_cover ) );
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 ) );
113 SCENARIO(
"DigitalConvexity< Z2 > fully convex triangles",
"[convex_simplices][2d]" )
122 DConvexity dconv(
Point( -1, -1 ),
Point( 5, 5 ) );
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;
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;
141 THEN(
"All 2737 invalid triangles are degenerated " ) {
143 REQUIRE( nb_notsimplex == nb_degenerated );
144 REQUIRE( nb_degenerated == 2737 );
146 THEN(
"There are 12888 valid triangles" ) {
147 REQUIRE( nb_unitary + nb_common == 12888 );
149 THEN(
"There are fewer (1920) unitary triangles than common triangles (10968)" ) {
152 REQUIRE( nb_unitary < nb_common );
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 );
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;
168 if ( ! ( ( a < b ) && ( b < c ) ) )
continue;
171 bool cvx0 = dconv.
isKConvex( triangle, 0 );
172 bool cvx1 = dconv.
isKConvex( triangle, 1 );
173 bool cvx2 = dconv.
isKConvex( triangle, 2 );
178 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
180 THEN(
"All valid triangles are 0-convex." ) {
183 THEN(
"There are less 1-convex and 2-convex than 0-convex." ) {
187 THEN(
"When the triangle is 0-convex and 1-convex, then it is 2-convex." ) {
193 SCENARIO(
"DigitalConvexity< Z3 > fully convex tetrahedra",
"[convex_simplices][3d]" )
202 DConvexity dconv(
Point( -1, -1, -1 ),
Point( 4, 4, 4 ) );
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;
215 if ( ! ( ( a < b ) && ( b < c ) && ( c < d ) ) )
continue;
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;
223 THEN(
"All 4228 invalid tetrahedra are degenerated " ) {
225 REQUIRE( nb_notsimplex == nb_degenerated );
226 REQUIRE( nb_degenerated == 4228 );
228 THEN(
"There are 35483 valid tetrahedra" ) {
229 REQUIRE( nb_unitary + nb_common == 35483 );
231 THEN(
"There are fewer (2515) unitary triangles than common triangles (32968)" ) {
234 REQUIRE( nb_unitary < nb_common );
236 THEN(
"The total number of triangles (unitary, common, degenerated) is 39711" ) {
237 REQUIRE( nb_unitary + nb_common + nb_degenerated == 39711 );
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 )
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 );
260 std::vector< Point > 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;
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;
286 THEN(
"All valid tetrahedra are 0-convex." ) {
289 THEN(
"There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
294 THEN(
"When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
300 THEN(
"All methods for computing full convexity agree." ) {
307 SCENARIO(
"DigitalConvexity< Z3 > rational fully convex tetrahedra",
"[convex_simplices][3d][rational]" )
313 DConvexity dconv(
Point( -1, -1, -1 ),
Point( 10, 10, 10 ) );
314 WHEN(
"Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
315 const unsigned int nb = 50;
316 unsigned int nbsimplex= 0;
317 unsigned int nb0 = 0;
318 unsigned int nb1 = 0;
319 unsigned int nb2 = 0;
320 unsigned int nb3 = 0;
321 unsigned int nb012_not3 = 0;
322 unsigned int nbf = 0;
323 unsigned int nb0123 = 0;
324 for (
unsigned int i = 0; i < nb; ++i )
326 Point a( 2*(rand() % 10), rand() % 20, 2*(rand() % 10) );
327 Point b( rand() % 20, 2*(rand() % 10), 2*(rand() % 10) );
328 Point c( 2*(rand() % 10), 2*(rand() % 10), rand() % 20 );
329 Point d( 2*(rand() % 10), 2*(rand() % 10), 2*(rand() % 10) );
343 nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
344 nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
346 THEN(
"All valid tetrahedra are 0-convex." ) {
349 THEN(
"There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
354 THEN(
"When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
369 SCENARIO(
"DigitalConvexity< Z2 > rational fully convex tetrahedra",
"[convex_simplices][2d][rational]" )
375 DConvexity dconv(
Point( -1, -1 ),
Point( 10, 10 ) );
376 WHEN(
"Computing many triangle in domain (0,0)-(9,9)." ) {
377 const unsigned int nb = 50;
378 unsigned int nbsimplex= 0;
379 unsigned int nb0 = 0;
380 unsigned int nb1 = 0;
381 unsigned int nb2 = 0;
382 unsigned int nb01_not2 = 0;
383 unsigned int nbf = 0;
384 unsigned int nb012 = 0;
385 for (
unsigned int i = 0; i < nb; ++i )
387 Point a( 2*(rand() % 10), rand() % 20 );
388 Point b( rand() % 20 , 2*(rand() % 10) );
389 Point c( 2*(rand() % 10), 2*(rand() % 10) );
401 nb012 += ( cvx0 && cvx1 && cvx2 ) ? 1 : 0;
402 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
404 THEN(
"All valid tetrahedra are 0-convex." ) {
407 THEN(
"There are less 1-convex, 2-convex than 0-convex." ) {
411 THEN(
"When the tetrahedron is 0-convex, and 1-convex, then it is 2-convex." ) {
424 SCENARIO(
"DigitalConvexity< Z3 > full subconvexity of segments and triangles",
"[subconvexity][3d]" )
433 DConvexity dconv(
Point( -6, -6, -6 ),
Point( 6, 6, 6 ) );
435 WHEN(
"Computing many tetrahedra" ) {
436 const unsigned int nb = 50;
437 unsigned int nb_fulldim = 0;
438 unsigned int nb_ok_seg = 0;
439 unsigned int nb_ok_tri = 0;
440 for (
unsigned int l = 0; l < nb; ++l )
442 const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
443 const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
444 const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
445 const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
446 const std::vector<Point> pts = { a, b, c, d };
448 nb_fulldim += fulldim ? 1 : 0;
449 if ( ! fulldim )
continue;
450 auto simplex = dconv.
makeSimplex( pts.cbegin(), pts.cend() );
453 unsigned int nb_subconvex = 0;
454 unsigned int nb_total = 0;
455 for (
unsigned int i = 0; i < 4; i++ )
456 for (
unsigned int j = i+1; j < 4; j++ )
458 auto segment = dconv.
makeSimplex( { pts[ i ], pts[ j ] } );
460 nb_subconvex += ok ? 1 : 0;
463 trace.
info() <<
"****** SEGMENT NOT SUBCONVEX ******" << std::endl;
464 trace.
info() <<
"splx v =" << a << b << c << d << std::endl;
465 trace.
info() <<
"simplex=" << simplex << std::endl;
466 trace.
info() <<
"seg v =" << pts[ i ] << pts[ j ] << std::endl;
467 trace.
info() <<
"segment=" << segment << std::endl;
470 nb_ok_seg += ( nb_subconvex == nb_total ) ? 1 : 0;
473 unsigned int nb_subconvex = 0;
474 unsigned int nb_total = 0;
475 for (
unsigned int i = 0; i < 4; i++ )
476 for (
unsigned int j = i+1; j < 4; j++ )
477 for (
unsigned int k = j+1; k < 4; k++ )
479 auto triangle = dconv.
makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
481 nb_subconvex += ok ? 1 : 0;
484 trace.
info() <<
"****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
485 trace.
info() <<
"splx v =" << a << b << c << d << std::endl;
486 trace.
info() <<
"simplex=" << simplex << std::endl;
487 trace.
info() <<
"tri v =" << pts[ i ] << pts[ j ]
488 << pts[ k ] << std::endl;
489 trace.
info() <<
"triangle=" << triangle << std::endl;
492 nb_ok_tri += ( nb_subconvex == nb_total ) ? 1 : 0;
495 THEN(
"At least half the tetrahedra are full dimensional." ) {
496 REQUIRE( nb_fulldim >= nb / 2 );
498 THEN(
"All segments of a tetrahedron should be subconvex to it." ) {
499 REQUIRE( nb_ok_seg == nb_fulldim );
501 THEN(
"All triangles of a tetrahedron should be subconvex to it." ) {
502 REQUIRE( nb_ok_tri == nb_fulldim );
507 SCENARIO(
"DigitalConvexity< Z3 > full convexity of polyhedra",
"[full_convexity][3d]" )
516 DConvexity dconv(
Point( -36, -36, -36 ),
Point( 36, 36, 36 ) );
518 const unsigned int nb = 10;
519 unsigned int nbfg = 0;
520 unsigned int nbffast = 0;
521 unsigned int nbfenv = 0;
523 std::vector< PointRange > XX;
524 for (
unsigned int i = 0; i < nb; ++i )
526 unsigned int k = 100;
528 for (
unsigned int j = 0; j < k; ++ j )
529 X[ j ] =
Point( rand() % 10, rand() % 10, rand() % 10 );
537 for (
const auto& X : XX )
540 nbfg += fcvx ? 1 : 0;
542 double t1 = c.stopClock();
544 for (
const auto& X : XX )
547 nbffast += fcvx ? 1 : 0;
549 double t2 = c.stopClock();
551 for (
const auto& X : XX )
553 auto card = dconv.
envelope( X ).size();
554 bool fcvx = card == X.size();
555 nbfenv += fcvx ? 1 : 0;
557 double t3 = c.stopClock();
558 WHEN(
"Computing many polytopes." ) {
559 THEN(
"All three methods agree on full convexity results" ) {
570 SCENARIO(
"DigitalConvexity< Z4 > full convexity of polyhedra",
"[full_convexity][4d]" )
579 DConvexity dconv(
Point( -36, -36, -36, -36 ),
Point( 36, 36, 36, 36 ) );
581 const unsigned int nb = 4;
582 unsigned int nbfg = 0;
583 unsigned int nbffast = 0;
584 unsigned int nbfenv = 0;
586 std::vector< PointRange > XX;
587 for (
unsigned int i = 0; i < nb; ++i )
589 unsigned int k = 100;
591 for (
unsigned int j = 0; j < k; ++ j )
592 X[ j ] =
Point( rand() % 8, rand() % 8, rand() % 8, rand() % 8 );
600 for (
const auto& X : XX )
603 nbfg += fcvx ? 1 : 0;
605 double t1 = c.stopClock();
607 for (
const auto& X : XX )
610 nbffast += fcvx ? 1 : 0;
612 double t2 = c.stopClock();
614 for (
const auto& X : XX )
616 auto card = dconv.
envelope( X ).size();
617 bool fcvx = card == X.size();
618 nbfenv += fcvx ? 1 : 0;
620 double t3 = c.stopClock();
621 WHEN(
"Computing many polytopes." ) {
622 THEN(
"All three methods agree on full convexity results" ) {
632 SCENARIO(
"DigitalConvexity< Z2 > sub-convexity of polyhedra",
"[full_subconvexity][2d]" )
638 DConvexity dconv(
Point( -36, -36 ),
Point( 36, 36 ) );
640 std::vector< Point > X( k );
641 X[ 0 ] =
Point( 0,0 );
642 X[ 1 ] =
Point( 7,-2 );
643 X[ 2 ] =
Point( 3,6 );
644 X[ 3 ] =
Point( 5, 5 );
645 X[ 4 ] =
Point( 2, 3 );
646 X[ 5 ] =
Point( -1, 1 );
650 REQUIRE( CG.nbCells() ==
L.size() );
651 for (
unsigned int i = 0; i < k; i++ )
652 for (
unsigned int j = i+1; j < k; j++ )
654 std::vector< Point > Z { X[ i ], X[ j ] };
660 REQUIRE( tangent_old == tangent_new );
661 REQUIRE( tangent_ab_old == tangent_ab_new );
662 REQUIRE( tangent_new == tangent_ab_new );
666 SCENARIO(
"DigitalConvexity< Z3 > sub-convexity of polyhedra",
"[full_subconvexity][3d]" )
672 DConvexity dconv(
Point( -36, -36, -36 ),
Point( 36, 36, 36 ) );
673 std::vector< Point > X( 5 );
674 X[ 0 ] =
Point( 0,0,0 );
675 X[ 1 ] =
Point( 0,5,1 );
676 X[ 2 ] =
Point( 2,1,6 );
677 X[ 3 ] =
Point( 6,1,1 );
678 X[ 4 ] =
Point( -2,-2,-3 );
682 std::vector< Point > Y;
684 REQUIRE( CG.nbCells() ==
L.size() );
686 unsigned int nb_ok = 0;
687 unsigned int nb_tgt= 0;
688 for (
int i = 0; i < 100; i++ )
690 Point a( rand() % 6, rand() % 6, rand() % 6 );
691 Point b( rand() % 6, rand() % 6, rand() % 6 );
695 nb_tgt += tangent_ab_new ? 1 : 0;
696 nb_ok += ( tangent_ab_old == tangent_ab_new ) ? 1 : 0;
704 SCENARIO(
"DigitalConvexity< Z3 > envelope",
"[envelope][3d]" )
710 DConvexity dconv(
Point( -36, -36, -36 ),
Point( 36, 36, 36 ) );
712 WHEN(
"Computing the envelope Z of a digital set X with direct algorithm" ) {
713 THEN(
"Z contains X" ){
714 for (
int k = 0; k < 5; k++ )
716 int n = 3 + ( rand() % 7 );
718 for (
int i = 0; i < n; i++ )
719 S.insert(
Point( rand() % 10, rand() % 10, rand() % 10 ) );
720 std::vector< Point > X( S.cbegin(), S.cend() );
722 auto Z = dconv.
envelope( X, DConvexity::EnvelopeAlgorithm::DIRECT );
725 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
726 X.cbegin(), X.cend() );
727 REQUIRE( X.size() <= Z.size() );
730 THEN(
"Z is fully convex" ){
731 for (
int k = 0; k < 5; k++ )
733 int n = 3 + ( rand() % 7 );
735 for (
int i = 0; i < n; i++ )
736 S.insert(
Point( rand() % 10, rand() % 10, rand() % 10 ) );
737 std::vector< Point > X( S.cbegin(), S.cend() );
745 WHEN(
"Computing the envelope Z of a digital set X with LatticeSet algorithm" ) {
746 THEN(
"Z contains X" ){
747 for (
int k = 0; k < 5; k++ )
749 int n = 3 + ( rand() % 7 );
751 for (
int i = 0; i < n; i++ )
752 S.insert(
Point( rand() % 10, rand() % 10, rand() % 10 ) );
753 std::vector< Point > X( S.cbegin(), S.cend() );
754 auto Z = dconv.
envelope( X, DConvexity::EnvelopeAlgorithm::LATTICE_SET );
756 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
757 X.cbegin(), X.cend() );
758 REQUIRE( X.size() <= Z.size() );
761 THEN(
"Z is fully convex" ){
762 for (
int k = 0; k < 5; k++ )
764 int n = 3 + ( rand() % 7 );
766 for (
int i = 0; i < n; i++ )
767 S.insert(
Point( rand() % 10, rand() % 10, rand() % 10 ) );
768 std::vector< Point > X( S.cbegin(), S.cend() );
778 SCENARIO(
"DigitalConvexity< Z2 > envelope",
"[envelope][2d]" )
784 DConvexity dconv(
Point( -360, -360 ),
Point( 360, 360 ) );
786 WHEN(
"Computing the envelope Z of two points" ) {
787 THEN(
"it requires at most one iteration" ){
788 for (
int k = 0; k < 10; k++ )
790 std::vector< Point > X;
791 X.push_back(
Point( rand() % 100, rand() % 100 ) );
792 X.push_back(
Point( rand() % 100, rand() % 100 ) );
793 std::sort( X.begin(), X.end() );
801 SCENARIO(
"DigitalConvexity< Z2 > relative envelope",
"[rel_envelope][2d]" )
807 DConvexity dconv(
Point( -360, -360 ),
Point( 360, 360 ) );
809 std::vector< Point > X {
Point( -10, -7 ),
Point( 10, 7 ) };
810 std::vector< Point > Y {
Point( -11, -6 ),
Point( 9, 8 ) };
811 std::sort( X.begin(), X.end() );
812 std::sort( Y.begin(), Y.end() );
817 WHEN(
"Computing the envelope of X relative to Y and Y relative to X" ) {
820 THEN(
"Both sets are fully convex" ){
824 THEN(
"There are inclusion rules between sets" ){
827 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
828 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
829 REQUIRE( std::includes( X.cbegin(), X.cend(),
830 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
833 WHEN(
"Computing the envelope of X relative to Y specified by a predicate" ) {
834 auto PredY = [] (
Point p )
835 {
return ( -4 <= p.dot(
Point( 2,5 ) ) ) && ( p.dot(
Point( 2,5 ) ) < 9 ); };
837 THEN(
"It is fully convex and included in Y" ){
842 for (
auto p : FC_X_rel_Y )
844 nb_in += PredY( p ) ? 1 : 0;
852 SCENARIO(
"DigitalConvexity< Z3 > relative envelope",
"[rel_envelope][3d]" )
858 DConvexity dconv(
Point( -360, -360, -360 ),
Point( 360, 360, 360 ) );
860 std::vector< Point > X {
Point( -61, -20, -8 ),
Point( 43, 25, 9 ) };
861 std::vector< Point > Y {
Point( -50, -27, -10 ),
Point( 40, 37, 17 ) };
862 std::sort( X.begin(), X.end() );
863 std::sort( Y.begin(), Y.end() );
868 WHEN(
"Computing the envelope of X relative to Y and Y relative to X" ) {
871 THEN(
"Both sets are fully convex" ){
875 THEN(
"There are inclusion rules between sets" ){
878 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
879 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
880 REQUIRE( std::includes( X.cbegin(), X.cend(),
881 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
887 SCENARIO(
"DigitalConvexity< Z3 > full subconvexity of triangles",
"[subconvexity][3d]" )
896 DConvexity dconv(
Point( -6, -6, -6 ),
Point( 6, 6, 6 ) );
899 WHEN(
"Computing many tetrahedra" ) {
900 const unsigned int nb = 50;
901 unsigned int nb_fulldim = 0;
902 unsigned int nb_ok_tri1 = 0;
903 unsigned int nb_ok_tri2 = 0;
904 for (
unsigned int l = 0; l < nb; ++l )
906 const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
907 const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
908 const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
909 const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
910 const std::vector<Point> pts = { a, b, c, d };
912 nb_fulldim += fulldim ? 1 : 0;
913 if ( ! fulldim )
continue;
914 auto simplex = dconv.
makeSimplex( pts.cbegin(), pts.cend() );
918 unsigned int nb_subconvex1 = 0;
919 unsigned int nb_subconvex2 = 0;
920 unsigned int nb_total = 0;
921 for (
unsigned int i = 0; i < 4; i++ )
922 for (
unsigned int j = i+1; j < 4; j++ )
923 for (
unsigned int k = j+1; k < 4; k++ )
925 auto tri1 = dconv.
makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
928 nb_subconvex1 += ok1 ? 1 : 0;
929 nb_subconvex2 += ok2 ? 1 : 0;
932 trace.
info() <<
"****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
933 trace.
info() <<
"splx v =" << a << b << c << d << std::endl;
934 trace.
info() <<
"simplex=" << simplex << std::endl;
935 trace.
info() <<
"tri v =" << pts[ i ] << pts[ j ]
936 << pts[ k ] << std::endl;
937 trace.
info() <<
"tri1=" << tri1 << std::endl;
940 trace.
info() <<
"****** TRIANGLE 3D NOT SUBCONVEX ****" << std::endl;
941 trace.
info() <<
"splx v =" << a << b << c << d << std::endl;
942 trace.
info() <<
"simplex=" << simplex << std::endl;
943 trace.
info() <<
"tri v =" << pts[ i ] << pts[ j ]
944 << pts[ k ] << std::endl;
947 nb_ok_tri1 += ( nb_subconvex1 == nb_total ) ? 1 : 0;
948 nb_ok_tri2 += ( nb_subconvex2 == nb_total ) ? 1 : 0;
951 THEN(
"All triangles of a tetrahedron should be subconvex to it." ) {
952 REQUIRE( nb_ok_tri1 == nb_fulldim );
954 THEN(
"All 3D triangles of a tetrahedron should be subconvex to it." ) {
955 REQUIRE( nb_ok_tri2 == nb_fulldim );
961 SCENARIO(
"DigitalConvexity< Z3 > full subconvexity of points and triangles",
"[subconvexity][3d]" )
971 DConvexity dconv (
Point( -21, -21, -21 ),
Point( 21, 21, 21 ) );
973 WHEN(
"Computing many tetrahedra" ) {
974 const unsigned int nb = 50;
975 unsigned int nb_total = 0;
976 unsigned int nb_ok_tri = 0;
977 unsigned int nb_subconvex1 = 0;
978 unsigned int nb_subconvex2 = 0;
979 for (
unsigned int l = 0; l < nb; ++l )
981 const Point a { (rand() % 10 - 10), (rand() % 10 - 10), (rand() % 10 - 10) };
982 const Point b { (rand() % 20 ), (rand() % 20 - 10), (rand() % 20 - 10) };
983 const Point c { (rand() % 20 - 10), (rand() % 20 ), (rand() % 20 - 10) };
984 const Point d { (rand() % 20 - 10), (rand() % 20 - 10), (rand() % 20 ) };
985 const std::vector<Point> pts = { a, b, c, d };
987 if ( ! fulldim )
continue;
988 auto simplex = dconv.
makeSimplex( pts.cbegin(), pts.cend() );
992 for (
unsigned int i = 0; i < 100; i++ )
994 const Point p { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
995 const Point q { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
996 const Point r { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
998 if ( n == Vector::zero )
continue;
1002 nb_subconvex1 += ok1 ? 1 : 0;
1003 nb_subconvex2 += ok2 ? 1 : 0;
1005 std::cout <<
"***** FULL SUBCONVEXITY ERROR ON TRIANGLE ****" << std::endl;
1006 std::cout <<
"splx v =" << a << b << c << d << std::endl;
1007 std::cout <<
"simplex=" << simplex << std::endl;
1008 std::cout <<
"tri v =" << p << q << r << std::endl;
1009 std::cout <<
"tri1=" << tri1 << std::endl;
1010 std::cout <<
"tri1 is fully subconvex: " << ( ok1 ?
"YES" :
"NO" ) << std::endl;
1011 std::cout <<
"3 points are fully subconvex: " << ( ok2 ?
"YES" :
"NO" ) << std::endl;
1013 nb_ok_tri += ( ok1 == ok2 ) ? 1 : 0;
1018 THEN(
"The number of triangles and point triplets subconvex to it should be equal." ) {
1019 REQUIRE( nb_subconvex1 == nb_subconvex2 );
1021 THEN(
"Full subconvexity should agree on every subset." ) {
1022 REQUIRE( nb_ok_tri == nb_total );
void getPoints(std::vector< Point > &pts) const
static PointRange insidePoints(const LatticePolytope &polytope)
static RationalPolytope makeRationalSimplex(Integer d, PointIterator itB, PointIterator itE)
bool isKConvex(const LatticePolytope &P, const Dimension k) const
static bool isSimplexFullDimensional(PointIterator itB, PointIterator itE)
static LatticePolytope makeSimplex(PointIterator itB, PointIterator itE)
bool isFullyConvexFast(const PointRange &X) const
bool isFullyConvex(const PointRange &X, bool convex0=false) const
PointRange envelope(const PointRange &Z, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
PointRange relativeEnvelope(const PointRange &Z, const PointRange &Y, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
bool isFullySubconvex(const PointRange &Y, const LatticeSet &StarX) const
LatticeSet StarCvxH(const PointRange &X, Dimension axis=dimension) const
Size depthLastEnvelope() const
LatticePolytope makePolytope(const PointRange &X, bool make_minkowski_summable=false) const
CellGeometry makeCellCover(PointIterator itB, PointIterator itE, Dimension i=0, Dimension k=KSpace::dimension) const
static SimplexType simplexType(PointIterator itB, PointIterator itE)
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::vector< Point > PointRange
DigitalPlane::Point Vector
DGtal is the top-level namespace which contains all DGtal functions and types.
auto crossProduct(PointVector< 3, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< 3, RightEuclideanRing, RightContainer > const &rhs) -> decltype(DGtal::constructFromArithmeticConversion(lhs, rhs))
Cross product of two 3D Points/Vectors.
GIVEN("A cubical complex with random 3-cells")
SCENARIO("DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]")
HyperRectDomain< Space > Domain
REQUIRE(domain.isInside(aPoint))