DGtal  1.2.0
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 
42 using namespace std;
43 using namespace DGtal;
44 
45 
47 // Functions for testing class DigitalConvexity.
49 
50 SCENARIO( "DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]" )
51 {
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 
113 SCENARIO( "DigitalConvexity< Z2 > fully convex triangles", "[convex_simplices][2d]" )
114 {
116  typedef KSpace::Point Point;
117  typedef KSpace::Space Space;
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 }
193 SCENARIO( "DigitalConvexity< Z3 > fully convex tetrahedra", "[convex_simplices][3d]" )
194 {
196  typedef KSpace::Point Point;
197  typedef KSpace::Space Space;
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 = 100;
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 nb0123 = 0;
250  for ( unsigned int i = 0; i < nb; ++i )
251  {
252  Point a( rand() % 5, rand() % 5, rand() % 5 );
253  Point b( rand() % 5, rand() % 5, rand() % 5 );
254  Point c( rand() % 5, rand() % 5, rand() % 5 );
255  Point d( rand() % 5, rand() % 5, rand() % 5 );
256  if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
257  auto tetra = dconv.makeSimplex( { a, b, c, d } );
258  bool cvx0 = dconv.isKConvex( tetra, 0 );
259  bool cvx1 = dconv.isKConvex( tetra, 1 );
260  bool cvx2 = dconv.isKConvex( tetra, 2 );
261  bool cvx3 = dconv.isKConvex( tetra, 3 );
262  bool cvxf = dconv.isFullyConvex( tetra );
263  nbsimplex += 1;
264  nb0 += cvx0 ? 1 : 0;
265  nb1 += cvx1 ? 1 : 0;
266  nb2 += cvx2 ? 1 : 0;
267  nb3 += cvx3 ? 1 : 0;
268  nbf += cvxf ? 1 : 0;
269  nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
270  nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
271  }
272  THEN( "All valid tetrahedra are 0-convex." ) {
273  REQUIRE( nb0 == nbsimplex );
274  }
275  THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
276  REQUIRE( nb1 < nb0 );
277  REQUIRE( nb2 < nb0 );
278  REQUIRE( nb3 < nb0 );
279  }
280  THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
281  REQUIRE( nb1 <= nb3 );
282  REQUIRE( nb2 <= nb3 );
283  REQUIRE( nb012_not3 == 0 );
284  REQUIRE( nbf == nb0123 );
285  }
286  }
287 }
288 
289 SCENARIO( "DigitalConvexity< Z3 > rational fully convex tetrahedra", "[convex_simplices][3d][rational]" )
290 {
292  typedef KSpace::Point Point;
293  typedef KSpace::Space Space;
294  typedef DigitalConvexity< KSpace > DConvexity;
295 
296  DConvexity dconv( Point( -1, -1, -1 ), Point( 10, 10, 10 ) );
297  WHEN( "Computing many tetrahedra in domain (0,0,0)-(4,4,4)." ) {
298  const unsigned int nb = 100;
299  unsigned int nbsimplex= 0;
300  unsigned int nb0 = 0;
301  unsigned int nb1 = 0;
302  unsigned int nb2 = 0;
303  unsigned int nb3 = 0;
304  unsigned int nb012_not3 = 0;
305  unsigned int nbf = 0;
306  unsigned int nb0123 = 0;
307  for ( unsigned int i = 0; i < nb; ++i )
308  {
309  Point a( 2*(rand() % 10), rand() % 20, 2*(rand() % 10) );
310  Point b( rand() % 20, 2*(rand() % 10), 2*(rand() % 10) );
311  Point c( 2*(rand() % 10), 2*(rand() % 10), rand() % 20 );
312  Point d( 2*(rand() % 10), 2*(rand() % 10), 2*(rand() % 10) );
313  if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
314  auto tetra = dconv.makeRationalSimplex( { Point(2,2,2), a, b, c, d } );
315  bool cvx0 = dconv.isKConvex( tetra, 0 );
316  bool cvx1 = dconv.isKConvex( tetra, 1 );
317  bool cvx2 = dconv.isKConvex( tetra, 2 );
318  bool cvx3 = dconv.isKConvex( tetra, 3 );
319  bool cvxf = dconv.isFullyConvex( tetra );
320  nbsimplex += 1;
321  nb0 += cvx0 ? 1 : 0;
322  nb1 += cvx1 ? 1 : 0;
323  nb2 += cvx2 ? 1 : 0;
324  nb3 += cvx3 ? 1 : 0;
325  nbf += cvxf ? 1 : 0;
326  nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
327  nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
328  }
329  THEN( "All valid tetrahedra are 0-convex." ) {
330  REQUIRE( nb0 == nbsimplex );
331  }
332  THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
333  REQUIRE( nb1 < nb0 );
334  REQUIRE( nb2 < nb0 );
335  REQUIRE( nb3 < nb0 );
336  }
337  THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
338  CAPTURE( nb0 );
339  CAPTURE( nb1 );
340  CAPTURE( nb2 );
341  CAPTURE( nb3 );
342  CAPTURE( nb012_not3 );
343  CAPTURE( nbf );
344  REQUIRE( nb2 <= nb3 );
345  REQUIRE( nb012_not3 == 0 );
346  REQUIRE( nbf == nb0123 );
347  }
348  }
349 }
350 
351 
352 SCENARIO( "DigitalConvexity< Z2 > rational fully convex tetrahedra", "[convex_simplices][2d][rational]" )
353 {
355  typedef KSpace::Point Point;
356  typedef KSpace::Space Space;
357  typedef DigitalConvexity< KSpace > DConvexity;
358 
359  DConvexity dconv( Point( -1, -1 ), Point( 10, 10 ) );
360  WHEN( "Computing many triangle in domain (0,0)-(9,9)." ) {
361  const unsigned int nb = 100;
362  unsigned int nbsimplex= 0;
363  unsigned int nb0 = 0;
364  unsigned int nb1 = 0;
365  unsigned int nb2 = 0;
366  unsigned int nb01_not2 = 0;
367  unsigned int nbf = 0;
368  unsigned int nb012 = 0;
369  for ( unsigned int i = 0; i < nb; ++i )
370  {
371  Point a( 2*(rand() % 10), rand() % 20 );
372  Point b( rand() % 20 , 2*(rand() % 10) );
373  Point c( 2*(rand() % 10), 2*(rand() % 10) );
374  if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
375  auto tetra = dconv.makeRationalSimplex( { Point(2,2), a, b, c } );
376  bool cvx0 = dconv.isKConvex( tetra, 0 );
377  bool cvx1 = dconv.isKConvex( tetra, 1 );
378  bool cvx2 = dconv.isKConvex( tetra, 2 );
379  bool cvxf = dconv.isFullyConvex( tetra );
380  nbsimplex += 1;
381  nb0 += cvx0 ? 1 : 0;
382  nb1 += cvx1 ? 1 : 0;
383  nb2 += cvx2 ? 1 : 0;
384  nbf += cvxf ? 1 : 0;
385  nb012 += ( cvx0 && cvx1 && cvx2 ) ? 1 : 0;
386  nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
387  }
388  THEN( "All valid tetrahedra are 0-convex." ) {
389  REQUIRE( nb0 == nbsimplex );
390  }
391  THEN( "There are less 1-convex, 2-convex than 0-convex." ) {
392  REQUIRE( nb1 < nb0 );
393  REQUIRE( nb2 < nb0 );
394  }
395  THEN( "When the tetrahedron is 0-convex, and 1-convex, then it is 2-convex." ) {
396  CAPTURE( nb0 );
397  CAPTURE( nb1 );
398  CAPTURE( nb2 );
399  CAPTURE( nb01_not2 );
400  CAPTURE( nbf );
401  REQUIRE( nb1 <= nb2 );
402  REQUIRE( nb01_not2 == 0 );
403  REQUIRE( nbf == nb012 );
404  }
405  }
406 }
407 
408 SCENARIO( "DigitalConvexity< Z3 > full subconvexity of segments and triangles", "[subconvexity][3d]" )
409 {
411  typedef KSpace::Point Point;
412  typedef KSpace::Space Space;
414  typedef DigitalConvexity< KSpace > DConvexity;
415 
416  Domain domain( Point( -5, -5, -5 ), Point( 5, 5, 5 ) );
417  DConvexity dconv( Point( -6, -6, -6 ), Point( 6, 6, 6 ) );
418 
419  WHEN( "Computing many tetrahedra" ) {
420  const unsigned int nb = 100;
421  unsigned int nb_fulldim = 0;
422  unsigned int nb_ok_seg = 0;
423  unsigned int nb_ok_tri = 0;
424  for ( unsigned int l = 0; l < nb; ++l )
425  {
426  const Point a { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
427  const Point b { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
428  const Point c { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
429  const Point d { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) };
430  const std::vector<Point> pts = { a, b, c, d };
431  const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
432  nb_fulldim += fulldim ? 1 : 0;
433  if ( ! fulldim ) continue;
434  auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
435  auto cover = dconv.makeCellCover( simplex, 0, 3 );
436  {
437  unsigned int nb_subconvex = 0;
438  unsigned int nb_total = 0;
439  for ( unsigned int i = 0; i < 4; i++ )
440  for ( unsigned int j = i+1; j < 4; j++ )
441  {
442  auto segment = dconv.makeSimplex( { pts[ i ], pts[ j ] } );
443  bool ok = dconv.isFullySubconvex( segment, cover );
444  nb_subconvex += ok ? 1 : 0;
445  nb_total += 1;
446  if ( ! ok ) {
447  trace.info() << "****** SEGMENT NOT SUBCONVEX ******" << std::endl;
448  trace.info() << "splx v =" << a << b << c << d << std::endl;
449  trace.info() << "simplex=" << simplex << std::endl;
450  trace.info() << "seg v =" << pts[ i ] << pts[ j ] << std::endl;
451  trace.info() << "segment=" << segment << std::endl;
452  }
453  }
454  nb_ok_seg += ( nb_subconvex == nb_total ) ? 1 : 0;
455  }
456  {
457  unsigned int nb_subconvex = 0;
458  unsigned int nb_total = 0;
459  for ( unsigned int i = 0; i < 4; i++ )
460  for ( unsigned int j = i+1; j < 4; j++ )
461  for ( unsigned int k = j+1; k < 4; k++ )
462  {
463  auto triangle = dconv.makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
464  bool ok = dconv.isFullySubconvex( triangle, cover );
465  nb_subconvex += ok ? 1 : 0;
466  nb_total += 1;
467  if ( ! ok ) {
468  trace.info() << "****** TRIANGLE NOT SUBCONVEX ****" << std::endl;
469  trace.info() << "splx v =" << a << b << c << d << std::endl;
470  trace.info() << "simplex=" << simplex << std::endl;
471  trace.info() << "tri v =" << pts[ i ] << pts[ j ]
472  << pts[ k ] << std::endl;
473  trace.info() << "triangle=" << triangle << std::endl;
474  }
475  }
476  nb_ok_tri += ( nb_subconvex == nb_total ) ? 1 : 0;
477  }
478  }
479  THEN( "At least half the tetrahedra are full dimensional." ) {
480  REQUIRE( nb_fulldim >= nb / 2 );
481  }
482  THEN( "All segments of a tetrahedron should be subconvex to it." ) {
483  REQUIRE( nb_ok_seg == nb_fulldim );
484  }
485  THEN( "All triangles of a tetrahedron should be subconvex to it." ) {
486  REQUIRE( nb_ok_tri == nb_fulldim );
487  }
488  }
489 }
Aim: A helper class to build polytopes from digital sets and to check digital k-convexity.
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::ostream & info()
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
SCENARIO("DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]")
Domain domain
HyperRectDomain< Space > Domain
REQUIRE(domain.isInside(aPoint))