DGtal  1.4.beta
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 {
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 
113 SCENARIO( "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 }
193 SCENARIO( "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 
307 SCENARIO( "DigitalConvexity< Z3 > rational fully convex tetrahedra", "[convex_simplices][3d][rational]" )
308 {
309  typedef KhalimskySpaceND<3,int> KSpace;
310  typedef KSpace::Point Point;
311  typedef DigitalConvexity< KSpace > DConvexity;
312 
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 )
325  {
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) );
330  if ( ! dconv.isSimplexFullDimensional( { a, b, c, d } ) ) continue;
331  auto tetra = dconv.makeRationalSimplex( { Point(2,2,2), a, b, c, d } );
332  bool cvx0 = dconv.isKConvex( tetra, 0 );
333  bool cvx1 = dconv.isKConvex( tetra, 1 );
334  bool cvx2 = dconv.isKConvex( tetra, 2 );
335  bool cvx3 = dconv.isKConvex( tetra, 3 );
336  bool cvxf = dconv.isFullyConvex( tetra );
337  nbsimplex += 1;
338  nb0 += cvx0 ? 1 : 0;
339  nb1 += cvx1 ? 1 : 0;
340  nb2 += cvx2 ? 1 : 0;
341  nb3 += cvx3 ? 1 : 0;
342  nbf += cvxf ? 1 : 0;
343  nb0123 += ( cvx0 && cvx1 && cvx2 && cvx3 ) ? 1 : 0;
344  nb012_not3+= ( cvx0 && cvx1 && cvx2 && ! cvx3 ) ? 1 : 0;
345  }
346  THEN( "All valid tetrahedra are 0-convex." ) {
347  REQUIRE( nb0 == nbsimplex );
348  }
349  THEN( "There are less 1-convex, 2-convex and 3-convex than 0-convex." ) {
350  REQUIRE( nb1 < nb0 );
351  REQUIRE( nb2 < nb0 );
352  REQUIRE( nb3 < nb0 );
353  }
354  THEN( "When the tetrahedron is 0-convex, 1-convex and 2-convex, then it is 3-convex." ) {
355  CAPTURE( nb0 );
356  CAPTURE( nb1 );
357  CAPTURE( nb2 );
358  CAPTURE( nb3 );
359  CAPTURE( nb012_not3 );
360  CAPTURE( nbf );
361  REQUIRE( nb2 <= nb3 );
362  REQUIRE( nb012_not3 == 0 );
363  REQUIRE( nbf == nb0123 );
364  }
365  }
366 }
367 
368 
369 SCENARIO( "DigitalConvexity< Z2 > rational fully convex tetrahedra", "[convex_simplices][2d][rational]" )
370 {
371  typedef KhalimskySpaceND<2,int> KSpace;
372  typedef KSpace::Point Point;
373  typedef DigitalConvexity< KSpace > DConvexity;
374 
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 )
386  {
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) );
390  if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
391  auto tetra = dconv.makeRationalSimplex( { Point(2,2), a, b, c } );
392  bool cvx0 = dconv.isKConvex( tetra, 0 );
393  bool cvx1 = dconv.isKConvex( tetra, 1 );
394  bool cvx2 = dconv.isKConvex( tetra, 2 );
395  bool cvxf = dconv.isFullyConvex( tetra );
396  nbsimplex += 1;
397  nb0 += cvx0 ? 1 : 0;
398  nb1 += cvx1 ? 1 : 0;
399  nb2 += cvx2 ? 1 : 0;
400  nbf += cvxf ? 1 : 0;
401  nb012 += ( cvx0 && cvx1 && cvx2 ) ? 1 : 0;
402  nb01_not2 += ( cvx0 && cvx1 && ! cvx2 ) ? 1 : 0;
403  }
404  THEN( "All valid tetrahedra are 0-convex." ) {
405  REQUIRE( nb0 == nbsimplex );
406  }
407  THEN( "There are less 1-convex, 2-convex than 0-convex." ) {
408  REQUIRE( nb1 < nb0 );
409  REQUIRE( nb2 < nb0 );
410  }
411  THEN( "When the tetrahedron is 0-convex, and 1-convex, then it is 2-convex." ) {
412  CAPTURE( nb0 );
413  CAPTURE( nb1 );
414  CAPTURE( nb2 );
415  CAPTURE( nb01_not2 );
416  CAPTURE( nbf );
417  REQUIRE( nb1 <= nb2 );
418  REQUIRE( nb01_not2 == 0 );
419  REQUIRE( nbf == nb012 );
420  }
421  }
422 }
423 
424 SCENARIO( "DigitalConvexity< Z3 > full subconvexity of segments and triangles", "[subconvexity][3d]" )
425 {
426  typedef KhalimskySpaceND<3,int> KSpace;
427  typedef KSpace::Point Point;
428  typedef KSpace::Space Space;
429  typedef HyperRectDomain< Space > Domain;
430  typedef DigitalConvexity< KSpace > DConvexity;
431 
432  Domain domain( Point( -5, -5, -5 ), Point( 5, 5, 5 ) );
433  DConvexity dconv( Point( -6, -6, -6 ), Point( 6, 6, 6 ) );
434 
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 )
441  {
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 };
447  const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend());
448  nb_fulldim += fulldim ? 1 : 0;
449  if ( ! fulldim ) continue;
450  auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() );
451  auto cover = dconv.makeCellCover( simplex, 0, 3 );
452  {
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++ )
457  {
458  auto segment = dconv.makeSimplex( { pts[ i ], pts[ j ] } );
459  bool ok = dconv.isFullySubconvex( segment, cover );
460  nb_subconvex += ok ? 1 : 0;
461  nb_total += 1;
462  if ( ! ok ) {
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;
468  }
469  }
470  nb_ok_seg += ( nb_subconvex == nb_total ) ? 1 : 0;
471  }
472  {
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++ )
478  {
479  auto triangle = dconv.makeSimplex({ pts[ i ], pts[ j ], pts[ k ] });
480  bool ok = dconv.isFullySubconvex( triangle, cover );
481  nb_subconvex += ok ? 1 : 0;
482  nb_total += 1;
483  if ( ! ok ) {
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;
490  }
491  }
492  nb_ok_tri += ( nb_subconvex == nb_total ) ? 1 : 0;
493  }
494  }
495  THEN( "At least half the tetrahedra are full dimensional." ) {
496  REQUIRE( nb_fulldim >= nb / 2 );
497  }
498  THEN( "All segments of a tetrahedron should be subconvex to it." ) {
499  REQUIRE( nb_ok_seg == nb_fulldim );
500  }
501  THEN( "All triangles of a tetrahedron should be subconvex to it." ) {
502  REQUIRE( nb_ok_tri == nb_fulldim );
503  }
504  }
505 }
506 
507 SCENARIO( "DigitalConvexity< Z3 > full convexity of polyhedra", "[full_convexity][3d]" )
508 {
509  typedef KhalimskySpaceND<3,int> KSpace;
510  typedef KSpace::Point Point;
511  typedef KSpace::Space Space;
512  typedef HyperRectDomain< Space > Domain;
513  typedef DigitalConvexity< KSpace > DConvexity;
514 
515  Domain domain( Point( -35, -35, -35 ), Point( 35, 35, 35 ) );
516  DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
517 
518  const unsigned int nb = 10;
519  unsigned int nbfg = 0;
520  unsigned int nbffast = 0;
521  unsigned int nbfenv = 0;
522  typedef std::vector< Point > PointRange;
523  std::vector< PointRange > XX;
524  for ( unsigned int i = 0; i < nb; ++i )
525  {
526  unsigned int k = 100;
527  PointRange X( k );
528  for ( unsigned int j = 0; j < k; ++ j )
529  X[ j ] = Point( rand() % 10, rand() % 10, rand() % 10 );
530  auto P = dconv.makePolytope( X );
531  PointRange Y;
532  P.getPoints( Y );
533  XX.push_back( Y );
534  }
535  Clock c;
536  c.startClock();
537  for ( const auto& X : XX )
538  {
539  bool fcvx = dconv.isFullyConvex( X, false );
540  nbfg += fcvx ? 1 : 0;
541  }
542  double t1 = c.stopClock();
543  c.startClock();
544  for ( const auto& X : XX )
545  {
546  bool fcvx = dconv.isFullyConvexFast( X );
547  nbffast += fcvx ? 1 : 0;
548  }
549  double t2 = c.stopClock();
550  c.startClock();
551  for ( const auto& X : XX )
552  {
553  auto card = dconv.envelope( X ).size();
554  bool fcvx = card == X.size();
555  nbfenv += fcvx ? 1 : 0;
556  }
557  double t3 = c.stopClock();
558  WHEN( "Computing many polytopes." ) {
559  THEN( "All three methods agree on full convexity results" ) {
560  CAPTURE( t1 );
561  CAPTURE( t2 );
562  CAPTURE( t3 );
563  REQUIRE( nbfg == nbffast );
564  REQUIRE( nbfg == nbfenv );
565  }
566  }
567 }
568 
569 
570 SCENARIO( "DigitalConvexity< Z4 > full convexity of polyhedra", "[full_convexity][4d]" )
571 {
572  typedef KhalimskySpaceND<4,int> KSpace;
573  typedef KSpace::Point Point;
574  typedef KSpace::Space Space;
575  typedef HyperRectDomain< Space > Domain;
576  typedef DigitalConvexity< KSpace > DConvexity;
577 
578  Domain domain( Point( -35, -35, -35, -35 ), Point( 35, 35, 35, 35 ) );
579  DConvexity dconv( Point( -36, -36, -36, -36 ), Point( 36, 36, 36, 36 ) );
580 
581  const unsigned int nb = 4;
582  unsigned int nbfg = 0;
583  unsigned int nbffast = 0;
584  unsigned int nbfenv = 0;
585  typedef std::vector< Point > PointRange;
586  std::vector< PointRange > XX;
587  for ( unsigned int i = 0; i < nb; ++i )
588  {
589  unsigned int k = 100;
590  PointRange X( k );
591  for ( unsigned int j = 0; j < k; ++ j )
592  X[ j ] = Point( rand() % 8, rand() % 8, rand() % 8, rand() % 8 );
593  auto P = dconv.makePolytope( X );
594  PointRange Y;
595  P.getPoints( Y );
596  XX.push_back( Y );
597  }
598  Clock c;
599  c.startClock();
600  for ( const auto& X : XX )
601  {
602  bool fcvx = dconv.isFullyConvex( X, false );
603  nbfg += fcvx ? 1 : 0;
604  }
605  double t1 = c.stopClock();
606  c.startClock();
607  for ( const auto& X : XX )
608  {
609  bool fcvx = dconv.isFullyConvexFast( X );
610  nbffast += fcvx ? 1 : 0;
611  }
612  double t2 = c.stopClock();
613  c.startClock();
614  for ( const auto& X : XX )
615  {
616  auto card = dconv.envelope( X ).size();
617  bool fcvx = card == X.size();
618  nbfenv += fcvx ? 1 : 0;
619  }
620  double t3 = c.stopClock();
621  WHEN( "Computing many polytopes." ) {
622  THEN( "All three methods agree on full convexity results" ) {
623  CAPTURE( t1 );
624  CAPTURE( t2 );
625  CAPTURE( t3 );
626  REQUIRE( nbfg == nbffast );
627  REQUIRE( nbfg == nbfenv );
628  }
629  }
630 }
631 
632 SCENARIO( "DigitalConvexity< Z2 > sub-convexity of polyhedra", "[full_subconvexity][2d]" )
633 {
634  typedef KhalimskySpaceND<2,int> KSpace;
635  typedef KSpace::Point Point;
636  typedef DigitalConvexity< KSpace > DConvexity;
637 
638  DConvexity dconv( Point( -36, -36 ), Point( 36, 36 ) );
639  unsigned int k = 6;
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 );
647  auto P = dconv.makePolytope( X, true );
648  auto CG = dconv.makeCellCover( P, 0, 2 );
649  auto L = dconv.StarCvxH( X, 0 );
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++ )
653  {
654  std::vector< Point > Z { X[ i ], X[ j ] };
655  const auto Q = dconv.makePolytope( Z );
656  bool tangent_old = dconv.isFullySubconvex( Q, CG );
657  bool tangent_new = dconv.isFullySubconvex( Z, L );
658  bool tangent_ab_old = dconv.isFullySubconvex( X[ i ], X[ j ], CG );
659  bool tangent_ab_new = dconv.isFullySubconvex( X[ i ], X[ j ], L );
660  REQUIRE( tangent_old == tangent_new );
661  REQUIRE( tangent_ab_old == tangent_ab_new );
662  REQUIRE( tangent_new == tangent_ab_new );
663  }
664 }
665 
666 SCENARIO( "DigitalConvexity< Z3 > sub-convexity of polyhedra", "[full_subconvexity][3d]" )
667 {
668  typedef KhalimskySpaceND<3,int> KSpace;
669  typedef KSpace::Point Point;
670  typedef DigitalConvexity< KSpace > DConvexity;
671 
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 );
679  auto P = dconv.makePolytope( X, true );
680  auto CG = dconv.makeCellCover( P, 0, 3 );
681  auto L = dconv.StarCvxH( X, 0 );
682  std::vector< Point > Y;
683  P.getPoints( Y );
684  REQUIRE( CG.nbCells() == L.size() );
685  unsigned int nb = 0;
686  unsigned int nb_ok = 0;
687  unsigned int nb_tgt= 0;
688  for ( int i = 0; i < 100; i++ )
689  {
690  Point a( rand() % 6, rand() % 6, rand() % 6 );
691  Point b( rand() % 6, rand() % 6, rand() % 6 );
692  // Point b( rand() % 20 - 10, rand() % 20 - 10, rand() % 20 - 10 );
693  bool tangent_ab_old = dconv.isFullySubconvex( a, b, CG );
694  bool tangent_ab_new = dconv.isFullySubconvex( a, b, L );
695  nb_tgt += tangent_ab_new ? 1 : 0;
696  nb_ok += ( tangent_ab_old == tangent_ab_new ) ? 1 : 0;
697  nb += 1;
698  }
699  REQUIRE( nb == nb_ok );
700  REQUIRE( 0 < nb_tgt );
701  REQUIRE( nb_tgt < 100 );
702 }
703 
704 SCENARIO( "DigitalConvexity< Z3 > envelope", "[envelope][3d]" )
705 {
706  typedef KhalimskySpaceND<3,int> KSpace;
707  typedef KSpace::Point Point;
708  typedef DigitalConvexity< KSpace > DConvexity;
709 
710  DConvexity dconv( Point( -36, -36, -36 ), Point( 36, 36, 36 ) );
711 
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++ )
715  {
716  int n = 3 + ( rand() % 7 );
717  std::set< Point > S;
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() );
721  CAPTURE( X );
722  auto Z = dconv.envelope( X, DConvexity::EnvelopeAlgorithm::DIRECT );
723  CAPTURE( Z );
724  CAPTURE( dconv.depthLastEnvelope() );
725  bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
726  X.cbegin(), X.cend() );
727  REQUIRE( X.size() <= Z.size() );
728  REQUIRE( Z_includes_X );
729  }
730  THEN( "Z is fully convex" ){
731  for ( int k = 0; k < 5; k++ )
732  {
733  int n = 3 + ( rand() % 7 );
734  std::set< Point > S;
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() );
738  auto Z = dconv.envelope( X );
739  CAPTURE( dconv.depthLastEnvelope() );
740  REQUIRE( dconv.isFullyConvex( Z ) );
741  }
742  }
743  }
744  }
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++ )
748  {
749  int n = 3 + ( rand() % 7 );
750  std::set< Point > S;
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 );
755  CAPTURE( dconv.depthLastEnvelope() );
756  bool Z_includes_X = std::includes( Z.cbegin(), Z.cend(),
757  X.cbegin(), X.cend() );
758  REQUIRE( X.size() <= Z.size() );
759  REQUIRE( Z_includes_X );
760  }
761  THEN( "Z is fully convex" ){
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 );
770  CAPTURE( dconv.depthLastEnvelope() );
771  REQUIRE( dconv.isFullyConvex( Z ) );
772  }
773  }
774  }
775  }
776 }
777 
778 SCENARIO( "DigitalConvexity< Z2 > envelope", "[envelope][2d]" )
779 {
780  typedef KhalimskySpaceND<2,int> KSpace;
781  typedef KSpace::Point Point;
782  typedef DigitalConvexity< KSpace > DConvexity;
783 
784  DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
785 
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++ )
789  {
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() );
794  auto Z = dconv.envelope( X );
795  REQUIRE( dconv.depthLastEnvelope() <= 1 );
796  }
797  }
798  }
799 }
800 
801 SCENARIO( "DigitalConvexity< Z2 > relative envelope", "[rel_envelope][2d]" )
802 {
803  typedef KhalimskySpaceND<2,int> KSpace;
804  typedef KSpace::Point Point;
805  typedef DigitalConvexity< KSpace > DConvexity;
806 
807  DConvexity dconv( Point( -360, -360 ), Point( 360, 360 ) );
808 
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() );
813  X = dconv.envelope( X );
814  Y = dconv.envelope( Y );
815  REQUIRE( dconv.isFullyConvex( X ) );
816  REQUIRE( dconv.isFullyConvex( Y ) );
817  WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
818  auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
819  auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
820  THEN( "Both sets are fully convex" ){
821  REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
822  REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
823  }
824  THEN( "There are inclusion rules between sets" ){
825  CAPTURE( FC_X_rel_Y );
826  CAPTURE( FC_Y_rel_X );
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() ) );
831  }
832  }
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 ); };
836  auto FC_X_rel_Y = dconv.relativeEnvelope( X, PredY );
837  THEN( "It is fully convex and included in Y" ){
838  CAPTURE( FC_X_rel_Y );
839  REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
840  int nb = 0;
841  int nb_in = 0;
842  for ( auto p : FC_X_rel_Y )
843  {
844  nb_in += PredY( p ) ? 1 : 0;
845  nb += 1;
846  }
847  REQUIRE( nb == nb_in );
848  }
849  }
850 }
851 
852 SCENARIO( "DigitalConvexity< Z3 > relative envelope", "[rel_envelope][3d]" )
853 {
854  typedef KhalimskySpaceND<3,int> KSpace;
855  typedef KSpace::Point Point;
856  typedef DigitalConvexity< KSpace > DConvexity;
857 
858  DConvexity dconv( Point( -360, -360, -360 ), Point( 360, 360, 360 ) );
859 
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() );
864  X = dconv.envelope( X );
865  Y = dconv.envelope( Y );
866  REQUIRE( dconv.isFullyConvex( X ) );
867  REQUIRE( dconv.isFullyConvex( Y ) );
868  WHEN( "Computing the envelope of X relative to Y and Y relative to X" ) {
869  auto FC_X_rel_Y = dconv.relativeEnvelope( X, Y );
870  auto FC_Y_rel_X = dconv.relativeEnvelope( Y, X );
871  THEN( "Both sets are fully convex" ){
872  REQUIRE( dconv.isFullyConvex( FC_X_rel_Y ) );
873  REQUIRE( dconv.isFullyConvex( FC_Y_rel_X ) );
874  }
875  THEN( "There are inclusion rules between sets" ){
876  CAPTURE( FC_X_rel_Y );
877  CAPTURE( FC_Y_rel_X );
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() ) );
882  }
883  }
884 }
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.
DGtal::trace
Trace trace
Definition: Common.h:153
Point
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE
CAPTURE(thicknessHV)
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition: testCubicalComplex.cpp:70
SCENARIO
SCENARIO("DigitalConvexity< Z2 > unit tests", "[digital_convexity][2d]")
Definition: testDigitalConvexity.cpp:50
domain
Domain domain
Definition: testProjection.cpp:88
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
REQUIRE
REQUIRE(domain.isInside(aPoint))