DGtal 2.1.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 "DGtal/geometry/tools/AffineGeometry.h"
40#include "DGtal/geometry/tools/AffineBasis.h"
41#include "DGtalCatch.h"
43
44using namespace std;
45using namespace DGtal;
46
47
49// Functions for testing class DigitalConvexity.
51
52SCENARIO( "DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]" )
53{
54 typedef KhalimskySpaceND<2,int> KSpace;
55 typedef KSpace::Point Point;
56 typedef DigitalConvexity< KSpace > DConvexity;
57
58 DConvexity dconv( Point( -5, -5 ), Point( 10, 10 ) );
59
60 GIVEN( "Given a fat simplex { (0,0), (4,-1), (2,5) } " ) {
61 std::vector<Point> V = { Point(0,0), Point(4,-1), Point(2,5) };
62 auto vertex_cover = dconv.makeCellCover( V.begin(), V.end() );
63 auto fat_simplex = dconv.makeSimplex ( V.begin(), V.end() );
64 auto inside_pts = dconv.insidePoints ( fat_simplex );
65 auto simplex_cover = dconv.makeCellCover( fat_simplex );
66 auto point_cover = dconv.makeCellCover( inside_pts.begin(), inside_pts.end() );
67 THEN( "The fat simplex is not degenerated." ) {
68 REQUIRE( dconv.isSimplexFullDimensional( V.begin(), V.end() ) );
69 }
70 THEN( "Its vertex cover contains 3 0-cells, 12 1-cells, 12 2-cells" ) {
71 REQUIRE( vertex_cover.computeNbCells( 0 ) == 3 );
72 REQUIRE( vertex_cover.computeNbCells( 1 ) == 12 );
73 REQUIRE( vertex_cover.computeNbCells( 2 ) == 12 );
74 }
75 THEN( "Its vertex cover is a subset of its point cover" ) {
76 REQUIRE( vertex_cover.subset( point_cover ) );
77 }
78 THEN( "Its point cover is a subset of its simplex cover" ) {
79 REQUIRE( point_cover.subset( simplex_cover ) );
80 }
81 THEN( "Being fat, its simplex cover is equal to its point cover" ) {
82 REQUIRE( simplex_cover.subset( point_cover ) );
83 }
84 }
85 GIVEN( "Given a thin simplex { (0,0), (4,3), (7,5) } " ) {
86 std::vector<Point> V = { Point(0,0), Point(4,3), Point(7,5) };
87 auto vertex_cover = dconv.makeCellCover( V.begin(), V.end() );
88 auto thin_simplex = dconv.makeSimplex ( V.begin(), V.end() );
89 auto inside_pts = dconv.insidePoints ( thin_simplex );
90 auto simplex_cover = dconv.makeCellCover( thin_simplex );
91 auto point_cover = dconv.makeCellCover( inside_pts.begin(), inside_pts.end() );
92 THEN( "The thin simplex is not degenerated." ) {
93 REQUIRE( dconv.isSimplexFullDimensional( V.begin(), V.end() ) );
94 }
95 THEN( "Its vertex cover contains 3 0-cells, 12 1-cells, 12 2-cells" ) {
96 REQUIRE( vertex_cover.computeNbCells( 0 ) == 3 );
97 REQUIRE( vertex_cover.computeNbCells( 1 ) == 12 );
98 REQUIRE( vertex_cover.computeNbCells( 2 ) == 12 );
99 }
100 THEN( "Its vertex cover is a subset of its point cover" ) {
101 REQUIRE( vertex_cover.subset( point_cover ) );
102 }
103 THEN( "Its point cover is a subset of its simplex cover" ) {
104 REQUIRE( point_cover.subset( simplex_cover ) );
105 }
106 THEN( "Being thin, its simplex cover is not equal to its point cover for 1<=dim<=2" ) {
107 REQUIRE( ! simplex_cover.subset( point_cover ) );
108 REQUIRE( simplex_cover.subset( point_cover, 0 ) );
109 REQUIRE( ! simplex_cover.subset( point_cover, 1 ) );
110 REQUIRE( ! simplex_cover.subset( point_cover, 2 ) );
111 }
112 }
113}
114
115SCENARIO( "DigitalConvexity< Z2 > fully convex triangles", "[convex_simplices][2d]" )
116{
117 typedef KhalimskySpaceND<2,int> KSpace;
118 typedef KSpace::Point Point;
119 typedef KSpace::Space Space;
120 typedef HyperRectDomain< Space > Domain;
121 typedef DigitalConvexity< KSpace > DConvexity;
122
123 Domain domain( Point( 0, 0 ), Point( 4, 4 ) );
124 DConvexity dconv( Point( -1, -1 ), Point( 5, 5 ) );
125
126 WHEN( "Computing all triangles in domain (0,0)-(4,4)." ) {
127 unsigned int nb_notsimplex = 0;
128 unsigned int nb_invalid = 0;
129 unsigned int nb_degenerated= 0;
130 unsigned int nb_common = 0;
131 unsigned int nb_unitary = 0;
132 for ( auto a : domain )
133 for ( auto b : domain )
134 for ( auto c : domain )
135 {
136 nb_notsimplex += ! dconv.isSimplexFullDimensional( { a, b, c } ) ? 1 : 0;
137 auto tri_type = dconv.simplexType( { a, b, c } );
138 nb_degenerated += tri_type == DConvexity::SimplexType::DEGENERATED ? 1 : 0;
139 nb_invalid += tri_type == DConvexity::SimplexType::INVALID ? 1 : 0;
140 nb_unitary += tri_type == DConvexity::SimplexType::UNITARY ? 1 : 0;
141 nb_common += tri_type == DConvexity::SimplexType::COMMON ? 1 : 0;
142 }
143 THEN( "All 2737 invalid triangles are degenerated " ) {
144 REQUIRE( nb_invalid == 0 );
145 REQUIRE( nb_notsimplex == nb_degenerated );
146 REQUIRE( nb_degenerated == 2737 );
147 }
148 THEN( "There are 12888 valid triangles" ) {
149 REQUIRE( nb_unitary + nb_common == 12888 );
150 }
151 THEN( "There are fewer (1920) unitary triangles than common triangles (10968)" ) {
152 REQUIRE( nb_unitary == 1920 );
153 REQUIRE( nb_common == 10968 );
154 REQUIRE( nb_unitary < nb_common );
155 }
156 THEN( "The total number of triangles (unitary, common, degenerated) is (domain size)^3, i.e. 5^6" ) {
157 REQUIRE( nb_unitary + nb_common + nb_degenerated == 5*5*5*5*5*5 );
158 }
159 }
160 WHEN( "Computing all triangles in domain (0,0)-(4,4)." ) {
161 unsigned int nbsimplex= 0;
162 unsigned int nb0 = 0;
163 unsigned int nb1 = 0;
164 unsigned int nb2 = 0;
165 unsigned int nb01_not2 = 0;
166 for ( auto a : domain )
167 for ( auto b : domain )
168 for ( auto c : domain )
169 {
170 if ( ! ( ( a < b ) && ( b < c ) ) ) continue;
171 if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
172 auto triangle = dconv.makeSimplex( { a, b, c } );
173 bool cvx0 = dconv.isKConvex( triangle, 0 );
174 bool cvx1 = dconv.isKConvex( triangle, 1 );
175 bool cvx2 = dconv.isKConvex( triangle, 2 );
176 nbsimplex += 1;
177 nb0 += cvx0 ? 1 : 0;
178 nb1 += cvx1 ? 1 : 0;
179 nb2 += cvx2 ? 1 : 0;
180 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
181 }
182 THEN( "All valid triangles are 0-convex." ) {
183 REQUIRE( nb0 == nbsimplex );
184 }
185 THEN( "There are less 1-convex and 2-convex than 0-convex." ) {
186 REQUIRE( nb1 < nb0 );
187 REQUIRE( nb2 < nb0 );
188 }
189 THEN( "When the triangle is 0-convex and 1-convex, then it is 2-convex." ) {
190 REQUIRE( nb1 <= nb2 );
191 REQUIRE( nb01_not2 == 0 );
192 }
193 }
194}
195SCENARIO( "DigitalConvexity< Z3 > fully convex tetrahedra", "[convex_simplices][3d]" )
196{
197 typedef KhalimskySpaceND<3,int> KSpace;
198 typedef KSpace::Point Point;
199 typedef KSpace::Space Space;
200 typedef HyperRectDomain< Space > Domain;
201 typedef DigitalConvexity< KSpace > DConvexity;
202
203 Domain domain( Point( 0, 0, 0 ), Point( 3, 3, 3 ) );
204 DConvexity dconv( Point( -1, -1, -1 ), Point( 4, 4, 4 ) );
205
206 WHEN( "Computing all lexicographically ordered tetrahedra anchored at (0,0,0) in domain (0,0,0)-(3,3,3)." ) {
207 unsigned int nb_notsimplex = 0;
208 unsigned int nb_invalid = 0;
209 unsigned int nb_degenerated= 0;
210 unsigned int nb_common = 0;
211 unsigned int nb_unitary = 0;
212 Point a(0,0,0);
213 for ( auto b : domain )
214 for ( auto c : domain )
215 for ( auto d : domain )
216 {
217 if ( ! ( ( a < b ) && ( b < c ) && ( c < d ) ) ) continue;
218 nb_notsimplex += ! dconv.isSimplexFullDimensional( {a,b,c,d} ) ? 1 : 0;
219 auto tri_type = dconv.simplexType( { a, b, c, d } );
220 nb_degenerated += tri_type == DConvexity::SimplexType::DEGENERATED ? 1 : 0;
221 nb_invalid += tri_type == DConvexity::SimplexType::INVALID ? 1 : 0;
222 nb_unitary += tri_type == DConvexity::SimplexType::UNITARY ? 1 : 0;
223 nb_common += tri_type == DConvexity::SimplexType::COMMON ? 1 : 0;
224 }
225 THEN( "All 4228 invalid tetrahedra are degenerated " ) {
226 REQUIRE( nb_invalid == 0 );
227 REQUIRE( nb_notsimplex == nb_degenerated );
228 REQUIRE( nb_degenerated == 4228 );
229 }
230 THEN( "There are 35483 valid tetrahedra" ) {
231 REQUIRE( nb_unitary + nb_common == 35483 );
232 }
233 THEN( "There are fewer (2515) unitary triangles than common triangles (32968)" ) {
234 REQUIRE( nb_unitary == 2515 );
235 REQUIRE( nb_common == 32968 );
236 REQUIRE( nb_unitary < nb_common );
237 }
238 THEN( "The total number of triangles (unitary, common, degenerated) is 39711" ) {
239 REQUIRE( nb_unitary + nb_common + nb_degenerated == 39711 );
240 }
241 }
242 WHEN( "Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
243 const unsigned int nb = 50;
244 unsigned int nbsimplex= 0;
245 unsigned int nb0 = 0;
246 unsigned int nb1 = 0;
247 unsigned int nb2 = 0;
248 unsigned int nb3 = 0;
249 unsigned int nb012_not3 = 0;
250 unsigned int nbf = 0;
251 unsigned int nbfg = 0;
252 unsigned int nbffast = 0;
253 unsigned int nb0123 = 0;
254 unsigned int nbdim3 = 0;
255 unsigned int nbfull = 0;
256 for ( unsigned int i = 0; i < nb; ++i )
257 {
258 Point a( rand() % 5, rand() % 5, rand() % 5 );
259 Point b( rand() % 5, rand() % 5, rand() % 5 );
260 Point c( rand() % 5, rand() % 5, rand() % 5 );
261 Point d( rand() % 5, rand() % 5, rand() % 5 );
262 if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
263 nbfull += 1;
264 auto dim = functions::computeAffineDimension( std::vector<Point>{ a, b, c, d } );
265 nbdim3 += ( dim == 3 ) ? 1 : 0;
266 auto tetra = dconv.makeSimplex( { a, b, c, d } );
267 std::vector< Point > X;
268 tetra.getPoints( X );
269 bool cvx0 = dconv.isKConvex( tetra, 0 );
270 bool cvx1 = dconv.isKConvex( tetra, 1 );
271 bool cvx2 = dconv.isKConvex( tetra, 2 );
272 bool cvx3 = dconv.isKConvex( tetra, 3 );
273 bool cvxf = dconv.isFullyConvex( tetra );
274 bool cvxfg = dconv.isFullyConvex( X, false );
275 bool cvxffast = dconv.isFullyConvexFast( X );
276 if ( cvxf != cvxfg || cvxf != cvxffast) {
277 bool cvxfc = dconv.FC( X ).size() == X.size();
278 std::cout << "[0123 " << cvx0 << cvx1 << cvx2 << cvx3 << "] "
279 << "[K " << cvxf << "] [M " << cvxfg
280 << "] [MF " << cvxffast << "] [FC " << cvxfc << "] "
281 << a << b << c << d << "\n";
282 std::cout << "X=";
283 for ( auto p : X ) std::cout << " " << p;
284 std::cout << "\n";
285 }
286 nbsimplex += 1;
287 nb0 += cvx0 ? 1 : 0;
288 nb1 += cvx1 ? 1 : 0;
289 nb2 += cvx2 ? 1 : 0;
290 nb3 += cvx3 ? 1 : 0;
291 nbf += cvxf ? 1 : 0;
292 nbfg += cvxfg ? 1 : 0;
293 nbffast += cvxffast ? 1 : 0;
294 nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
295 nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
296 }
297 THEN( "All valid tetrahedra are 0-convex." ) {
298 REQUIRE( nb0 == nbsimplex );
299 }
300 THEN( "All full dimensional simplices have affine dimension 3" ) {
301 REQUIRE( nbfull == nbdim3 );
302 }
303 THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
304 REQUIRE( nb1 < nb0 );
305 REQUIRE( nb2 < nb0 );
306 REQUIRE( nb3 < nb0 );
307 }
308 THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
309 REQUIRE( nb1 <= nb3 );
310 REQUIRE( nb2 <= nb3 );
311 REQUIRE( nb012_not3 == 0 );
312 REQUIRE( nbf == nb0123 );
313 }
314 THEN( "All methods for computing full convexity agree." ) {
315 REQUIRE( nbf == nbfg );
316 REQUIRE( nbf == nbffast );
317 }
318 }
319}
320
321SCENARIO( "DigitalConvexity< Z3 > rational fully convex tetrahedra", "[convex_simplices][3d][rational]" )
322{
323 typedef KhalimskySpaceND<3,int> KSpace;
324 typedef KSpace::Point Point;
325 typedef DigitalConvexity< KSpace > DConvexity;
326
327 DConvexity dconv( Point( -10, -10, -10 ), Point( 100, 100, 100 ) );
328 WHEN( "Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
329 const unsigned int nb = 30;
330 unsigned int nbsimplex= 0;
331 unsigned int nb0 = 0;
332 unsigned int nb1 = 0;
333 unsigned int nb2 = 0;
334 unsigned int nb3 = 0;
335 unsigned int nb012_not3 = 0;
336 unsigned int nbf = 0;
337 unsigned int nb0123 = 0;
338 for ( unsigned int i = 0; i < nb; ++i )
339 {
340 Point a( 2*(rand() % 10), rand() % 20, 2*(rand() % 10) );
341 Point b( rand() % 20, 2*(rand() % 10), 2*(rand() % 10) );
342 Point c( 2*(rand() % 10), 2*(rand() % 10), rand() % 20 );
343 Point d( 2*(rand() % 10), 2*(rand() % 10), 2*(rand() % 10) );
344 if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
345 auto tetra = dconv.makeRationalSimplex( { Point(2,2,2), a, b, c, d } );
346 bool cvx0 = dconv.isKConvex( tetra, 0 );
347 bool cvx1 = dconv.isKConvex( tetra, 1 );
348 bool cvx2 = dconv.isKConvex( tetra, 2 );
349 bool cvx3 = dconv.isKConvex( tetra, 3 );
350 bool cvxf = dconv.isFullyConvex( tetra );
351 nbsimplex += 1;
352 nb0 += cvx0 ? 1 : 0;
353 nb1 += cvx1 ? 1 : 0;
354 nb2 += cvx2 ? 1 : 0;
355 nb3 += cvx3 ? 1 : 0;
356 nbf += cvxf ? 1 : 0;
357 nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
358 nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
359 }
360 THEN( "All valid tetrahedra are 0-convex." ) {
361 REQUIRE( nb0 == nbsimplex );
362 }
363 THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
364 REQUIRE( nb1 <= nb0 );
365 REQUIRE( nb2 <= nb0 );
366 REQUIRE( nb3 <= nb0 );
367 }
368 THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
369 CAPTURE( nb0 );
370 CAPTURE( nb1 );
371 CAPTURE( nb2 );
372 CAPTURE( nb3 );
373 CAPTURE( nb012_not3 );
374 CAPTURE( nbf );
375 REQUIRE( nb2 <= nb3 );
376 REQUIRE( nb012_not3 == 0 );
377 REQUIRE( nbf == nb0123 );
378 }
379 }
380}
381
382
383SCENARIO( "DigitalConvexity< Z2 > rational fully convex triangles", "[convex_simplices][2d][rational]" )
384{
385 typedef KhalimskySpaceND<2,int> KSpace;
386 typedef KSpace::Point Point;
387 typedef DigitalConvexity< KSpace > DConvexity;
388
389 DConvexity dconv( Point( -1, -1 ), Point( 30, 30 ) );
390 WHEN( "Computing many triangle in domain (0,0)-(9,9)." ) {
391 const unsigned int nb = 50;
392 unsigned int nbsimplex= 0;
393 unsigned int nb0 = 0;
394 unsigned int nb1 = 0;
395 unsigned int nb2 = 0;
396 unsigned int nb01_not2 = 0;
397 unsigned int nbf = 0;
398 unsigned int nb012 = 0;
399 for ( unsigned int i = 0; i < nb; ++i )
400 {
401 Point a( 2*(rand() % 10), rand() % 20 );
402 Point b( rand() % 20 , 2*(rand() % 10) );
403 Point c( 2*(rand() % 10), 2*(rand() % 10) );
404 if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
405 auto triangle = dconv.makeRationalSimplex( { Point(2,2), a, b, c } );
406 bool cvx0 = dconv.isKConvex( triangle, 0 );
407 bool cvx1 = dconv.isKConvex( triangle, 1 );
408 bool cvx2 = dconv.isKConvex( triangle, 2 );
409 bool cvxf = dconv.isFullyConvex( triangle );
410 nbsimplex += 1;
411 nb0 += cvx0 ? 1 : 0;
412 nb1 += cvx1 ? 1 : 0;
413 nb2 += cvx2 ? 1 : 0;
414 nbf += cvxf ? 1 : 0;
415 nb012 += ( cvx0 && cvx1 && cvx2 ) ? 1 : 0;
416 nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
417 }
418 THEN( "All valid tetrahedra are 0-convex." ) {
419 REQUIRE( nb0 == nbsimplex );
420 }
421 THEN( "There are less 1-convex, 2-convex than 0-convex." ) {
422 REQUIRE( nb1 <= nb0 );
423 REQUIRE( nb2 <= nb0 );
424 }
425 THEN( "When the tetrahedron is 0-convex, and 1-convex, then it is 2-convex." ) {
426 CAPTURE( nb0 );
427 CAPTURE( nb1 );
428 CAPTURE( nb2 );
429 CAPTURE( nb01_not2 );
430 CAPTURE( nbf );
431 REQUIRE( nb1 <= nb2 );
432 REQUIRE( nb01_not2 == 0 );
433 REQUIRE( nbf == nb012 );
434 }
435 }
436}
437
438SCENARIO( "DigitalConvexity< Z3 > full subconvexity of segments and triangles", "[subconvexity][3d]" )
439{
440 typedef KhalimskySpaceND<3,int> KSpace;
441 typedef KSpace::Point Point;
442 typedef KSpace::Space Space;
443 typedef HyperRectDomain< Space > Domain;
444 typedef DigitalConvexity< KSpace > DConvexity;
445
446 Domain domain( Point( -5, -5, -5 ), Point( 5, 5, 5 ) );
447 DConvexity dconv( Point( -6, -6, -6 ), Point( 6, 6, 6 ) );
448
449 WHEN( "Computing many tetrahedra" ) {
450 const unsigned int nb = 50;
451 unsigned int nb_fulldim = 0;
452 unsigned int nb_ok_seg = 0;
453 unsigned int nb_ok_tri = 0;
454 for ( unsigned int l = 0; l < nb; ++l )
455 {
456 const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
457 const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
458 const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
459 const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
460 const std::vector<Point> pts = { a, b, c, d };
461 const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
462 nb_fulldim += fulldim ? 1 : 0;
463 if ( ! fulldim ) continue;
464 auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
465 auto cover = dconv.makeCellCover( simplex, 0, 3 );
466 {
467 unsigned int nb_subconvex = 0;
468 unsigned int nb_total = 0;
469 for ( unsigned int i = 0; i < 4; i++ )
470 for ( unsigned int j = i+1; j < 4; j++ )
471 {
472 auto segment = dconv.makeSimplex( { pts[ i ], pts[ j ] } );
473 bool ok = dconv.isFullySubconvex( segment, cover );
474 nb_subconvex += ok ? 1 : 0;
475 nb_total += 1;
476 if ( ! ok ) {
477 trace.info() << "****** SEGMENT NOT SUBCONVEX ******" << std::endl;
478 trace.info() << "splx v =" << a << b << c << d << std::endl;
479 trace.info() << "simplex=" << simplex << std::endl;
480 trace.info() << "seg v =" << pts[ i ] << pts[ j ] << std::endl;
481 trace.info() << "segment=" << segment << std::endl;
482 }
483 }
484 nb_ok_seg += ( nb_subconvex == nb_total ) ? 1 : 0;
485 }
486 {
487 unsigned int nb_subconvex = 0;
488 unsigned int nb_total = 0;
489 for ( unsigned int i = 0; i < 4; i++ )
490 for ( unsigned int j = i+1; j < 4; j++ )
491 for ( unsigned int k = j+1; k < 4; k++ )
492 {
493 auto triangle = dconv.makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
494 bool ok = dconv.isFullySubconvex( triangle, cover );
495 nb_subconvex += ok ? 1 : 0;
496 nb_total += 1;
497 if ( ! ok ) {
498 trace.info() << "****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
499 trace.info() << "splx v =" << a << b << c << d << std::endl;
500 trace.info() << "simplex=" << simplex << std::endl;
501 trace.info() << "tri v =" << pts[ i ] << pts[ j ]
502 << pts[ k ] << std::endl;
503 trace.info() << "triangle=" << triangle << std::endl;
504 }
505 }
506 nb_ok_tri += ( nb_subconvex == nb_total ) ? 1 : 0;
507 }
508 }
509 THEN( "At least half the tetrahedra are full dimensional." ) {
510 REQUIRE( nb_fulldim >= nb / 2 );
511 }
512 THEN( "All segments of a tetrahedron should be subconvex to it." ) {
513 REQUIRE( nb_ok_seg == nb_fulldim );
514 }
515 THEN( "All triangles of a tetrahedron should be subconvex to it." ) {
516 REQUIRE( nb_ok_tri == nb_fulldim );
517 }
518 }
519}
520
521SCENARIO( "DigitalConvexity< Z3 > full convexity of polyhedra", "[full_convexity][3d]" )
522{
523 typedef KhalimskySpaceND<3,int> KSpace;
524 typedef KSpace::Point Point;
525 typedef KSpace::Space Space;
526 typedef HyperRectDomain< Space > Domain;
527 typedef DigitalConvexity< KSpace > DConvexity;
528
529 Domain domain( Point( -35, -35, -35 ), Point( 35, 35, 35 ) );
530 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
531
532 const unsigned int nb = 10;
533 unsigned int nbfg = 0;
534 unsigned int nbffast = 0;
535 unsigned int nbfenv = 0;
536 typedef std::vector< Point > PointRange;
537 std::vector< PointRange > XX;
538 for ( unsigned int i = 0; i < nb; ++i )
539 {
540 unsigned int k = 100;
541 PointRange X( k );
542 for ( unsigned int j = 0; j < k; ++ j )
543 X[ j ] = Point( rand() % 10, rand() % 10, rand() % 10 );
544 auto P = dconv.makePolytope( X );
545 PointRange Y;
546 P.getPoints( Y );
547 XX.push_back( Y );
548 }
549 Clock c;
550 c.startClock();
551 for ( const auto& X : XX )
552 {
553 bool fcvx = dconv.isFullyConvex( X, false );
554 nbfg += fcvx ? 1 : 0;
555 }
556 double t1 = c.stopClock();
557 c.startClock();
558 for ( const auto& X : XX )
559 {
560 bool fcvx = dconv.isFullyConvexFast( X );
561 nbffast += fcvx ? 1 : 0;
562 }
563 double t2 = c.stopClock();
564 c.startClock();
565 for ( const auto& X : XX )
566 {
567 auto card = dconv.envelope( X ).size();
568 bool fcvx = card == X.size();
569 nbfenv += fcvx ? 1 : 0;
570 }
571 double t3 = c.stopClock();
572 WHEN( "Computing many polytopes." ) {
573 THEN( "All three methods agree on full convexity results" ) {
574 CAPTURE( t1 );
575 CAPTURE( t2 );
576 CAPTURE( t3 );
577 REQUIRE( nbfg == nbffast );
578 REQUIRE( nbfg == nbfenv );
579 }
580 }
581}
582
583
584SCENARIO( "DigitalConvexity< Z4 > full convexity of polyhedra", "[full_convexity][4d]" )
585{
586 typedef KhalimskySpaceND<4,int> KSpace;
587 typedef KSpace::Point Point;
588 typedef KSpace::Space Space;
589 typedef HyperRectDomain< Space > Domain;
590 typedef DigitalConvexity< KSpace > DConvexity;
591
592 Domain domain( Point( -35, -35, -35, -35 ), Point( 35, 35, 35, 35 ) );
593 DConvexity dconv( Point( -36, -36, -36, -36 ), Point( 36, 36, 36, 36 ) );
594
595 const unsigned int nb = 4;
596 unsigned int nbfg = 0;
597 unsigned int nbffast = 0;
598 unsigned int nbfenv = 0;
599 typedef std::vector< Point > PointRange;
600 std::vector< PointRange > XX;
601 for ( unsigned int i = 0; i < nb; ++i )
602 {
603 unsigned int k = 100;
604 PointRange X( k );
605 for ( unsigned int j = 0; j < k; ++ j )
606 X[ j ] = Point( rand() % 8, rand() % 8, rand() % 8, rand() % 8 );
607 auto P = dconv.makePolytope( X );
608 PointRange Y;
609 P.getPoints( Y );
610 XX.push_back( Y );
611 }
612 Clock c;
613 c.startClock();
614 for ( const auto& X : XX )
615 {
616 bool fcvx = dconv.isFullyConvex( X, false );
617 nbfg += fcvx ? 1 : 0;
618 }
619 double t1 = c.stopClock();
620 c.startClock();
621 for ( const auto& X : XX )
622 {
623 bool fcvx = dconv.isFullyConvexFast( X );
624 nbffast += fcvx ? 1 : 0;
625 }
626 double t2 = c.stopClock();
627 c.startClock();
628 for ( const auto& X : XX )
629 {
630 auto card = dconv.envelope( X ).size();
631 bool fcvx = card == X.size();
632 nbfenv += fcvx ? 1 : 0;
633 }
634 double t3 = c.stopClock();
635 WHEN( "Computing many polytopes." ) {
636 THEN( "All three methods agree on full convexity results" ) {
637 CAPTURE( t1 );
638 CAPTURE( t2 );
639 CAPTURE( t3 );
640 REQUIRE( nbfg == nbffast );
641 REQUIRE( nbfg == nbfenv );
642 }
643 }
644}
645
646SCENARIO( "DigitalConvexity< Z2 > sub-convexity of polyhedra", "[full_subconvexity][2d]" )
647{
648 typedef KhalimskySpaceND<2,int> KSpace;
649 typedef KSpace::Point Point;
650 typedef DigitalConvexity< KSpace > DConvexity;
651
652 DConvexity dconv( Point( -36, -36 ), Point( 36, 36 ) );
653 unsigned int k = 6;
654 std::vector< Point > X( k );
655 X[ 0 ] = Point( 0,0 );
656 X[ 1 ] = Point( 7,-2 );
657 X[ 2 ] = Point( 3,6 );
658 X[ 3 ] = Point( 5, 5 );
659 X[ 4 ] = Point( 2, 3 );
660 X[ 5 ] = Point( -1, 1 );
661 auto P = dconv.makePolytope( X, true );
662 auto CG = dconv.makeCellCover( P, 0, 2 );
663 auto L = dconv.StarCvxH( X, 0 );
664 REQUIRE( CG.nbCells() == L.size() );
665 for ( unsigned int i = 0; i < k; i++ )
666 for ( unsigned int j = i+1; j < k; j++ )
667 {
668 std::vector< Point > Z { X[ i ], X[ j ] };
669 const auto Q = dconv.makePolytope( Z );
670 bool tangent_old = dconv.isFullySubconvex( Q, CG );
671 bool tangent_new = dconv.isFullySubconvex( Z, L );
672 bool tangent_ab_old = dconv.isFullySubconvex( X[ i ], X[ j ], CG );
673 bool tangent_ab_new = dconv.isFullySubconvex( X[ i ], X[ j ], L );
674 REQUIRE( tangent_old == tangent_new );
675 REQUIRE( tangent_ab_old == tangent_ab_new );
676 REQUIRE( tangent_new == tangent_ab_new );
677 }
678}
679
680SCENARIO( "DigitalConvexity< Z3 > sub-convexity of polyhedra", "[full_subconvexity][3d]" )
681{
682 typedef KhalimskySpaceND<3,int> KSpace;
683 typedef KSpace::Point Point;
684 typedef DigitalConvexity< KSpace > DConvexity;
685
686 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
687 std::vector< Point > X( 5 );
688 X[ 0 ] = Point( 0,0,0 );
689 X[ 1 ] = Point( 0,5,1 );
690 X[ 2 ] = Point( 2,1,6 );
691 X[ 3 ] = Point( 6,1,1 );
692 X[ 4 ] = Point( -2,-2,-3 );
693 auto P = dconv.makePolytope( X, true );
694 auto CG = dconv.makeCellCover( P, 0, 3 );
695 auto L = dconv.StarCvxH( X, 0 );
696 std::vector< Point > Y;
697 P.getPoints( Y );
698 REQUIRE( CG.nbCells() == L.size() );
699 unsigned int nb = 0;
700 unsigned int nb_ok = 0;
701 unsigned int nb_tgt= 0;
702 for ( int i = 0; i < 100; i++ )
703 {
704 Point a( rand() % 6, rand() % 6, rand() % 6 );
705 Point b( rand() % 6, rand() % 6, rand() % 6 );
706 // Point b( rand() % 20 - 10, rand() % 20 - 10, rand() % 20 - 10 );
707 bool tangent_ab_old = dconv.isFullySubconvex( a, b, CG );
708 bool tangent_ab_new = dconv.isFullySubconvex( a, b, L );
709 nb_tgt += tangent_ab_new ? 1 : 0;
710 nb_ok += ( tangent_ab_old == tangent_ab_new ) ? 1 : 0;
711 nb += 1;
712 }
713 REQUIRE( nb == nb_ok );
714 REQUIRE( 0 < nb_tgt );
715 REQUIRE( nb_tgt < 100 );
716}
717
718SCENARIO( "DigitalConvexity< Z3 > envelope", "[envelope][3d]" )
719{
720 typedef KhalimskySpaceND<3,int> KSpace;
721 typedef KSpace::Point Point;
722 typedef DigitalConvexity< KSpace > DConvexity;
723
724 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
725
726 WHEN( "Computing the envelope Z of a digital set X with direct algorithm" ) {
727 THEN( "Z contains X" ){
728 for ( int k = 0; k < 5; k++ )
729 {
730 int n = 3 + ( rand() % 7 );
731 std::set< Point > S;
732 for ( int i = 0; i < n; i++ )
733 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
734 std::vector< Point > X( S.cbegin(), S.cend() );
735 CAPTURE( X );
736 auto Z = dconv.envelope( X, DConvexity::EnvelopeAlgorithm::DIRECT );
737 CAPTURE( Z );
738 CAPTURE( dconv.depthLastEnvelope() );
739 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
740 X.cbegin(), X.cend() );
741 REQUIRE( X.size() <= Z.size() );
742 REQUIRE( Z_includes_X );
743 }
744 THEN( "Z is fully convex" ){
745 for ( int k = 0; k < 5; k++ )
746 {
747 int n = 3 + ( rand() % 7 );
748 std::set< Point > S;
749 for ( int i = 0; i < n; i++ )
750 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
751 std::vector< Point > X( S.cbegin(), S.cend() );
752 auto Z = dconv.envelope( X );
753 CAPTURE( X );
754 CAPTURE( dconv.depthLastEnvelope() );
755 REQUIRE( dconv.isFullyConvex( Z ) );
756 }
757 }
758 }
759 }
760 WHEN( "Computing the envelope Z of a digital set X with LatticeSet algorithm" ) {
761 THEN( "Z contains X" ){
762 for ( int k = 0; k < 5; k++ )
763 {
764 int n = 3 + ( rand() % 7 );
765 std::set< Point > S;
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() );
769 auto Z = dconv.envelope( X, DConvexity::EnvelopeAlgorithm::LATTICE_SET );
770 CAPTURE( dconv.depthLastEnvelope() );
771 bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
772 X.cbegin(), X.cend() );
773 REQUIRE( X.size() <= Z.size() );
774 REQUIRE( Z_includes_X );
775 }
776 THEN( "Z is fully convex" ){
777 for ( int k = 0; k < 5; k++ )
778 {
779 int n = 3 + ( rand() % 7 );
780 std::set< Point > S;
781 for ( int i = 0; i < n; i++ )
782 S.insert( Point( rand() % 10, rand() % 10, rand() % 10 ) );
783 std::vector< Point > X( S.cbegin(), S.cend() );
784 auto Z = dconv.envelope( X );
785 CAPTURE( dconv.depthLastEnvelope() );
786 REQUIRE( dconv.isFullyConvex( Z ) );
787 }
788 }
789 }
790 }
791}
792
793SCENARIO( "DigitalConvexity< Z2 > envelope", "[envelope][2d]" )
794{
795 typedef KhalimskySpaceND<2,int> KSpace;
796 typedef KSpace::Point Point;
797 typedef DigitalConvexity< KSpace > DConvexity;
798
799 DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
800
801 WHEN( "Computing the envelope Z of two points" ) {
802 THEN( "it requires at most one iteration" ){
803 for ( int k = 0; k < 10; k++ )
804 {
805 std::vector< Point > X;
806 X.push_back( Point( rand() % 100, rand() % 100 ) );
807 X.push_back( Point( rand() % 100, rand() % 100 ) );
808 std::sort( X.begin(), X.end() );
809 auto Z = dconv.envelope( X );
810 REQUIRE( dconv.depthLastEnvelope() <= 1 );
811 }
812 }
813 }
814}
815
816SCENARIO( "DigitalConvexity< Z2 > relative envelope", "[rel_envelope][2d]" )
817{
818 typedef KhalimskySpaceND<2,int> KSpace;
819 typedef KSpace::Point Point;
820 typedef DigitalConvexity< KSpace > DConvexity;
821
822 DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
823
824 std::vector< Point > X { Point( -10, -7 ), Point( 10, 7 ) };
825 std::vector< Point > Y { Point( -11, -6 ), Point( 9, 8 ) };
826 std::sort( X.begin(), X.end() );
827 std::sort( Y.begin(), Y.end() );
828 X = dconv.envelope( X );
829 Y = dconv.envelope( Y );
830 REQUIRE( dconv.isFullyConvex( X ) );
831 REQUIRE( dconv.isFullyConvex( Y ) );
832 WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
833 auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
834 auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
835 THEN( "Both sets are fully convex" ){
836 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
837 REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
838 }
839 THEN( "There are inclusion rules between sets" ){
840 CAPTURE( FC_X_rel_Y );
841 CAPTURE( FC_Y_rel_X );
842 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
843 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
844 REQUIRE( std::includes( X.cbegin(), X.cend(),
845 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
846 }
847 }
848 WHEN( "Computing the envelope of X relative to Y specified by a predicate" ) {
849 auto PredY = [] ( Point p )
850 { return ( -4 <= p.dot( Point( 2,5 ) ) ) && ( p.dot( Point( 2,5 ) ) < 9 ); };
851 auto FC_X_rel_Y = dconv.relativeEnvelope( X, PredY );
852 THEN( "It is fully convex and included in Y" ){
853 CAPTURE( FC_X_rel_Y );
854 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
855 int nb = 0;
856 int nb_in = 0;
857 for ( auto p : FC_X_rel_Y )
858 {
859 nb_in += PredY( p ) ? 1 : 0;
860 nb += 1;
861 }
862 REQUIRE( nb == nb_in );
863 }
864 }
865}
866
867SCENARIO( "DigitalConvexity< Z3 > relative envelope", "[rel_envelope][3d]" )
868{
869 typedef KhalimskySpaceND<3,int> KSpace;
870 typedef KSpace::Point Point;
871 typedef DigitalConvexity< KSpace > DConvexity;
872
873 DConvexity dconv( Point( -360, -360, -360 ), Point( 360, 360, 360 ) );
874
875 std::vector< Point > X { Point( -61, -20, -8 ), Point( 43, 25, 9 ) };
876 std::vector< Point > Y { Point( -50, -27, -10 ), Point( 40, 37, 17 ) };
877 std::sort( X.begin(), X.end() );
878 std::sort( Y.begin(), Y.end() );
879 X = dconv.envelope( X );
880 Y = dconv.envelope( Y );
881 REQUIRE( dconv.isFullyConvex( X ) );
882 REQUIRE( dconv.isFullyConvex( Y ) );
883 WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
884 auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
885 auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
886 THEN( "Both sets are fully convex" ){
887 REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
888 REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
889 }
890 THEN( "There are inclusion rules between sets" ){
891 CAPTURE( FC_X_rel_Y );
892 CAPTURE( FC_Y_rel_X );
893 REQUIRE( std::includes( Y.cbegin(), Y.cend(),
894 FC_X_rel_Y.cbegin(), FC_X_rel_Y.cend() ) );
895 REQUIRE( std::includes( X.cbegin(), X.cend(),
896 FC_Y_rel_X.cbegin(), FC_Y_rel_X.cend() ) );
897 }
898 }
899}
900
901
902SCENARIO( "DigitalConvexity< Z3 > full subconvexity of triangles", "[subconvexity][3d]" )
903{
904 typedef KhalimskySpaceND<3,int> KSpace;
905 typedef KSpace::Point Point;
906 typedef KSpace::Space Space;
907 typedef HyperRectDomain< Space > Domain;
908 typedef DigitalConvexity< KSpace > DConvexity;
909
910 Domain domain( Point( -5, -5, -5 ), Point( 5, 5, 5 ) );
911 DConvexity dconv( Point( -6, -6, -6 ), Point( 6, 6, 6 ) );
912
913
914 WHEN( "Computing many tetrahedra" ) {
915 const unsigned int nb = 50;
916 unsigned int nb_fulldim = 0;
917 unsigned int nb_ok_tri1 = 0;
918 unsigned int nb_ok_tri2 = 0;
919 for ( unsigned int l = 0; l < nb; ++l )
920 {
921 const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
922 const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
923 const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
924 const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
925 const std::vector<Point> pts = { a, b, c, d };
926 const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
927 nb_fulldim += fulldim ? 1 : 0;
928 if ( ! fulldim ) continue;
929 auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
930 auto cover = dconv.makeCellCover( simplex, 0, 3 );
931 auto ls = dconv.StarCvxH( pts );
932 {
933 unsigned int nb_subconvex1 = 0;
934 unsigned int nb_subconvex2 = 0;
935 unsigned int nb_total = 0;
936 for ( unsigned int i = 0; i < 4; i++ )
937 for ( unsigned int j = i+1; j < 4; j++ )
938 for ( unsigned int k = j+1; k < 4; k++ )
939 {
940 auto tri1 = dconv.makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
941 bool ok1 = dconv.isFullySubconvex( tri1, cover );
942 bool ok2 = dconv.isFullySubconvex( pts[ i ], pts[ j ], pts[ k ], ls );
943 nb_subconvex1 += ok1 ? 1 : 0;
944 nb_subconvex2 += ok2 ? 1 : 0;
945 nb_total += 1;
946 if ( ! ok1 ) {
947 trace.info() << "****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
948 trace.info() << "splx v =" << a << b << c << d << std::endl;
949 trace.info() << "simplex=" << simplex << std::endl;
950 trace.info() << "tri v =" << pts[ i ] << pts[ j ]
951 << pts[ k ] << std::endl;
952 trace.info() << "tri1=" << tri1 << std::endl;
953 }
954 if ( ! ok2 ) {
955 trace.info() << "****** TRIANGLE 3D NOT SUBCONVEX ****" << std::endl;
956 trace.info() << "splx v =" << a << b << c << d << std::endl;
957 trace.info() << "simplex=" << simplex << std::endl;
958 trace.info() << "tri v =" << pts[ i ] << pts[ j ]
959 << pts[ k ] << std::endl;
960 }
961 }
962 nb_ok_tri1 += ( nb_subconvex1 == nb_total ) ? 1 : 0;
963 nb_ok_tri2 += ( nb_subconvex2 == nb_total ) ? 1 : 0;
964 }
965 }
966 THEN( "All triangles of a tetrahedron should be subconvex to it." ) {
967 REQUIRE( nb_ok_tri1 == nb_fulldim );
968 }
969 THEN( "All 3D triangles of a tetrahedron should be subconvex to it." ) {
970 REQUIRE( nb_ok_tri2 == nb_fulldim );
971 }
972 }
973
974}
975
976SCENARIO( "DigitalConvexity< Z3 > full subconvexity of points and triangles", "[subconvexity][3d]" )
977{
978 typedef KhalimskySpaceND<3,int> KSpace;
979 typedef KSpace::Point Point;
980 typedef KSpace::Vector Vector;
981 typedef KSpace::Space Space;
982 typedef HyperRectDomain< Space > Domain;
983 typedef DigitalConvexity< KSpace > DConvexity;
984
985 Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) );
986 DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) );
987
988 WHEN( "Computing many tetrahedra" ) {
989 const unsigned int nb = 50;
990 unsigned int nb_total = 0;
991 unsigned int nb_ok_tri = 0;
992 unsigned int nb_subconvex1 = 0;
993 unsigned int nb_subconvex2 = 0;
994 for ( unsigned int l = 0; l < nb; ++l )
995 {
996 const Point a { (rand() % 10 - 10), (rand() % 10 - 10), (rand() % 10 - 10) };
997 const Point b { (rand() % 20 ), (rand() % 20 - 10), (rand() % 20 - 10) };
998 const Point c { (rand() % 20 - 10), (rand() % 20 ), (rand() % 20 - 10) };
999 const Point d { (rand() % 20 - 10), (rand() % 20 - 10), (rand() % 20 ) };
1000 const std::vector<Point> pts = { a, b, c, d };
1001 const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
1002 if ( ! fulldim ) continue;
1003 auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
1004 auto cover = dconv.makeCellCover( simplex, 0, 3 );
1005 auto ls = dconv.StarCvxH( pts );
1006 {
1007 for ( unsigned int i = 0; i < 100; i++ )
1008 {
1009 const Point p { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1010 const Point q { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1011 const Point r { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1012 const Vector n = ( q - p ).crossProduct( r - p );
1013 if ( n == Vector::zero ) continue;
1014 auto tri1 = dconv.makeSimplex( { p, q, r } );
1015 bool ok1 = dconv.isFullySubconvex( tri1, cover );
1016 bool ok2 = dconv.isFullySubconvex( p, q, r, ls );
1017 nb_subconvex1 += ok1 ? 1 : 0;
1018 nb_subconvex2 += ok2 ? 1 : 0;
1019 if ( ok1 != ok2 ) {
1020 std::cout << "***** FULL SUBCONVEXITY ERROR ON TRIANGLE ****" << std::endl;
1021 std::cout << "splx v =" << a << b << c << d << std::endl;
1022 std::cout << "simplex=" << simplex << std::endl;
1023 std::cout << "tri v =" << p << q << r << std::endl;
1024 std::cout << "tri1=" << tri1 << std::endl;
1025 std::cout << "tri1 is fully subconvex: " << ( ok1 ? "YES" : "NO" ) << std::endl;
1026 std::cout << "3 points are fully subconvex: " << ( ok2 ? "YES" : "NO" ) << std::endl;
1027 }
1028 nb_ok_tri += ( ok1 == ok2 ) ? 1 : 0;
1029 nb_total += 1;
1030 }
1031 }
1032 }
1033 THEN( "The number of triangles and point triplets subconvex to it should be equal." ) {
1034 REQUIRE( nb_subconvex1 == nb_subconvex2 );
1035 }
1036 THEN( "Full subconvexity should agree on every subset." ) {
1037 REQUIRE( nb_ok_tri == nb_total );
1038 }
1039 }
1040
1041}
1042
1043SCENARIO( "DigitalConvexity< Z3 > full covering of segments", "[full_cover][3d]" )
1044{
1045 typedef KhalimskySpaceND<3,int> KSpace;
1046 typedef KSpace::Point Point;
1047 // typedef KSpace::Vector Vector;
1048 typedef KSpace::Space Space;
1049 typedef HyperRectDomain< Space > Domain;
1050 typedef DigitalConvexity< KSpace > DConvexity;
1051
1052 Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) );
1053 DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) );
1054 {
1055 Point a( 0, 0, 0 );
1056 Point b( 3, 2, 1 );
1057 auto LS = dconv.CoverCvxH( a, b );
1058 auto P = LS.toPointRange();
1059 CAPTURE( P );
1060 REQUIRE( P.size() == 9 );
1061 }
1062 {
1063 Point a( 0, 0, 3 );
1064 Point b( 3, 0, 0 );
1065 auto LS = dconv.CoverCvxH( a, b );
1066 auto P = LS.toPointRange();
1067 CAPTURE( P );
1068 REQUIRE( P.size() == 7 );
1069 }
1070 {
1071 Point a( 3, 0, 0 );
1072 Point b( 0, 2, 5 );
1073 auto LS = dconv.CoverCvxH( a, b );
1074 auto P = LS.toPointRange();
1075 CAPTURE( P );
1076 REQUIRE( P.size() == 17 );
1077 }
1078}
1079
1080SCENARIO( "DigitalConvexity< Z3 > full covering of triangles", "[full_cover][3d]" )
1081{
1082 typedef KhalimskySpaceND<3,int> KSpace;
1083 typedef KSpace::Point Point;
1084// typedef KSpace::Vector Vector;
1085 typedef KSpace::Space Space;
1086 typedef HyperRectDomain< Space > Domain;
1087 typedef DigitalConvexity< KSpace > DConvexity;
1088
1089 Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) );
1090 DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) );
1091 {
1092 Point a( 0, 0, 0 );
1093 Point b( 2, 1, 0 );
1094 Point c( 2, 2, 0 );
1095 auto LS = dconv.CoverCvxH( a, b, c );
1096 auto P = LS.toPointRange();
1097 CAPTURE( P );
1098 REQUIRE( P.size() == 10 );
1099 }
1100 {
1101 Point a( 0, 0, 1 );
1102 Point b( 2, 1, 0 );
1103 Point c( 2, 2, 0 );
1104 auto LS = dconv.CoverCvxH( a, b, c );
1105 auto P = LS.toPointRange();
1106 CAPTURE( P );
1107 REQUIRE( P.size() == 10 );
1108 }
1109 {
1110 Point a( 1, 0, 0 );
1111 Point b( 0, 1, 0 );
1112 Point c( 0, 0, 1 );
1113 auto LS = dconv.CoverCvxH( a, b, c );
1114 auto P = LS.toPointRange();
1115 CAPTURE( P );
1116 REQUIRE( P.size() == 7 );
1117 }
1118}
1119
1120
1121
1122SCENARIO( "DigitalConvexity< Z3 > full subconvexity and full covering of triangles", "[subconvexity][3d][full_cover]" )
1123{
1124 typedef KhalimskySpaceND<3,int> KSpace;
1125 typedef KSpace::Point Point;
1126 typedef KSpace::Vector Vector;
1127 typedef KSpace::Space Space;
1128 typedef HyperRectDomain< Space > Domain;
1129 typedef DigitalConvexity< KSpace > DConvexity;
1130
1131 Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) );
1132 DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) );
1133
1134 WHEN( "Computing many tetrahedra" ) {
1135 const unsigned int nb = 50;
1136 unsigned int nb_total = 0;
1137 unsigned int nb_ok_tri = 0;
1138 unsigned int nb_ok_seg = 0;
1139 unsigned int nb_ko_seg = 0;
1140 unsigned int nb_subconvex = 0;
1141 unsigned int nb_covered = 0;
1142 for ( unsigned int l = 0; l < nb; ++l )
1143 {
1144 const Point a { (rand() % 10 - 10), (rand() % 10 - 10), (rand() % 10 - 10) };
1145 const Point b { (rand() % 20 ), (rand() % 20 - 10), (rand() % 20 - 10) };
1146 const Point c { (rand() % 20 - 10), (rand() % 20 ), (rand() % 20 - 10) };
1147 const Point d { (rand() % 20 - 10), (rand() % 20 - 10), (rand() % 20 ) };
1148 const std::vector<Point> pts = { a, b, c, d };
1149 const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
1150 if ( ! fulldim ) continue;
1151 auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
1152 auto cover = dconv.makeCellCover( simplex, 0, 3 );
1153 auto ls = dconv.StarCvxH( pts );
1154 {
1155 for ( unsigned int i = 0; i < 100; i++ )
1156 {
1157 const Point p { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1158 const Point q { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1159 const Point r { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
1160 const Vector n = ( q - p ).crossProduct( r - p );
1161 if ( n == Vector::zero ) continue;
1162 auto tri1 = dconv.makeSimplex( { p, q, r } );
1163 bool ok1 = dconv.isFullySubconvex( p, q, r, ls );
1164 bool ok2 = dconv.isFullyCovered( p, q, r, ls );
1165 if ( ok2 )
1166 {
1167 // check covering of segments
1168 if ( dconv.isFullyCovered( p, q, ls ) ) nb_ok_seg += 1;
1169 else nb_ko_seg += 1;
1170 if ( dconv.isFullyCovered( p, r, ls ) ) nb_ok_seg += 1;
1171 else nb_ko_seg += 1;
1172 if ( dconv.isFullyCovered( q, r, ls ) ) nb_ok_seg += 1;
1173 else nb_ko_seg += 1;
1174 }
1175
1176 nb_subconvex += ok1 ? 1 : 0;
1177 nb_covered += ok2 ? 1 : 0;
1178 bool ok = ok1 == ok2;
1179 if ( ! ok )
1180 {
1181 std::cout << "***** FULL SUBCONVEXITY ERROR ON TRIANGLE ****" << std::endl;
1182 std::cout << "splx v =" << a << b << c << d << std::endl;
1183 std::cout << "simplex=" << simplex << std::endl;
1184 std::cout << "tri v =" << p << q << r << std::endl;
1185 std::cout << "tri1=" << tri1 << std::endl;
1186 std::cout << "3 points are fully subconvex: " << ( ok1 ? "YES" : "NO" ) << std::endl;
1187 std::cout << "3 points are fully covered by star: " << ( ok2 ? "YES" : "NO" ) << std::endl;
1188 }
1189 nb_ok_tri += ( ok ) ? 1 : 0;
1190 nb_total += 1;
1191 }
1192 }
1193 }
1194 THEN( "The number of subconvex and covered triangles should be equal." ) {
1195 REQUIRE( nb_subconvex == nb_covered );
1196 }
1197 THEN( "Full subconvexity and covering should agree on every subset." ) {
1198 REQUIRE( nb_ok_tri == nb_total );
1199 }
1200 THEN( "When a triangle is fully covered, its edges are fully covered also." ) {
1201 REQUIRE( nb_ko_seg == 0 );
1202 }
1203 }
1204
1205}
1206
1207SCENARIO( "DigitalConvexity< Z3 > envelope bug", "[envelope][3d]" )
1208{
1209 typedef KhalimskySpaceND<3,int> KSpace;
1210 typedef KSpace::Point Point;
1211 typedef DigitalConvexity< KSpace > DConvexity;
1212 typedef AffineGeometry< Point > Affine;
1213 typedef AffineBasis< Point > Basis;
1214
1215 DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
1216
1217 WHEN( "Using basis B = (1, 0, -2) (1, 0, -1)" ) {
1218 std::vector< Point > b = { Point( 0, 0, 0), Point(1, 0, -2), Point(1, 0, -1) };
1219 const auto [ o, B ] = Affine::affineBasis( b );
1220 Point e = functions::computeIndependentVector( B );
1221 Basis AB( Point( 0,0 ), b, Basis::Type::ECHELON_REDUCED );
1222 bool parallel = AB.isParallel( e );
1223 const auto [ d, L, r ] = AB.decomposeVector( e );
1224 CAPTURE( B );
1225 CAPTURE( AB.basis() );
1226 CAPTURE( d );
1227 CAPTURE( L );
1228 CAPTURE( r );
1229 CAPTURE( e );
1230 REQUIRE( ! parallel );
1231 REQUIRE( r != Point::zero );
1232 }
1233
1234 WHEN( "Computing the envelope Z of a digital set X with direct algorithm" ) {
1235 std::vector< Point > X = { Point(5, 1, 9), Point(8, 1, 8), Point(9, 1, 1) };
1236 auto Z = dconv.envelope( X );
1237 THEN( "Z is fully convex" ){
1238 CAPTURE( X );
1239 CAPTURE( dconv.depthLastEnvelope() );
1240 REQUIRE( dconv.isFullyConvex( Z ) );
1241 }
1242 }
1243}
DGtal::BoundedLatticePolytope::getPoints
void getPoints(std::vector< Point > &pts) const
DGtal::Clock::startClock
void startClock()
DGtal::Clock::stopClock
double stopClock() const
DGtal::DigitalConvexity::insidePoints
static PointRange insidePoints(const LatticePolytope &polytope)
DGtal::DigitalConvexity::makeRationalSimplex
static RationalPolytope makeRationalSimplex(Integer d, PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::isFullyCovered
bool isFullyCovered(const Point &a, const Point &b, const Point &c, const LatticeSet &cells) const
DGtal::DigitalConvexity::FC
PointRange FC(const PointRange &Z, EnvelopeAlgorithm algo=EnvelopeAlgorithm::DIRECT) const
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::CoverCvxH
LatticeSet CoverCvxH(const Point &a, const Point &b, Dimension axis=dimension) 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
Vector
DigitalPlane::Point Vector
Definition examplePlaneProbingParallelepipedEstimator.cpp:45
DGtal::functions::computeAffineDimension
DGtal::int64_t computeAffineDimension(const std::vector< TPoint > &X, const double tolerance=1e-12)
Definition AffineGeometry.h:1068
DGtal::functions::computeIndependentVector
TPoint computeIndependentVector(const std::vector< TPoint > &basis, const double tolerance=1e-12)
Definition AffineGeometry.h:1202
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition ClosedIntegerHalfPlane.h:49
DGtal::crossProduct
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.
DGtal::trace
Trace trace
std
STL namespace.
DGtal::AffineBasis
Aim: Utility class to determine the affine geometry of an input set of points. It provides exact resu...
Definition AffineBasis.h:112
DGtal::AffineGeometry
Aim: Utility class to determine the affine geometry of an input set of points. It provides exact resu...
Definition AffineGeometry.h:454
Point
MyPointD Point
Definition testClone2.cpp:381
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