DGtal  1.5.beta
testConvexityHelper.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/geometry/volumes/ConvexityHelper.h"
37 #include "DGtal/shapes/SurfaceMesh.h"
38 #include "DGtalCatch.h"
40
41 using namespace std;
42 using namespace DGtal;
43
44
46 // Functions for testing class ConvexityHelper in 2D.
48
49 SCENARIO( "ConvexityHelper< 2 > unit tests",
50  "[convexity_helper][lattice_polytope][2d]" )
51 {
52  typedef ConvexityHelper< 2 > Helper;
53  typedef Helper::Point Point;
54  typedef ConvexCellComplex< Point > CvxCellComplex;
55  GIVEN( "Given a star { (0,0), (-4,-1), (-3,5), (7,3), (5, -2) } " ) {
56  std::vector<Point> V
57  = { Point(0,0), Point(-4,-1), Point(-3,5), Point(7,3), Point(5, -2) };
58  WHEN( "Computing its lattice polytope" ){
59  const auto P = Helper::computeLatticePolytope( V, false, true );
60  CAPTURE( P );
61  THEN( "The polytope is valid and has 4 non trivial facets" ) {
62  REQUIRE( P.nbHalfSpaces() - 4 == 4 );
63  }
64  THEN( "The polytope is Minkowski summable" ) {
65  REQUIRE( P.canBeSummed() );
66  }
67  THEN( "The polytope contains the input points" ) {
68  REQUIRE( P.isInside( V[ 0 ] ) );
69  REQUIRE( P.isInside( V[ 1 ] ) );
70  REQUIRE( P.isInside( V[ 2 ] ) );
71  REQUIRE( P.isInside( V[ 3 ] ) );
72  REQUIRE( P.isInside( V[ 4 ] ) );
73  }
74  THEN( "The polytope contains 58 points" ) {
75  REQUIRE( P.count() == 58 );
76  }
77  THEN( "The interior of the polytope contains 53 points" ) {
78  REQUIRE( P.countInterior() == 53 );
79  }
80  THEN( "The boundary of the polytope contains 5 points" ) {
81  REQUIRE( P.countBoundary() == 5 );
82  }
83  }
84  }
85  GIVEN( "Given a square with an additional outside vertex " ) {
86  std::vector<Point> V
87  = { Point(-10,-10), Point(10,-10), Point(-10,10), Point(10,10),
88  Point(0,18) };
89  WHEN( "Computing its Delaunay cell complex" ){
90  CvxCellComplex complex;
91  bool ok = Helper::computeDelaunayCellComplex( complex, V, false );
92  CAPTURE( complex );
93  THEN( "The complex has 2 cells, 6 faces, 5 vertices" ) {
94  REQUIRE( ok );
95  REQUIRE( complex.nbCells() == 2 );
96  REQUIRE( complex.nbFaces() == 6 );
97  REQUIRE( complex.nbVertices() == 5 );
98  }
99  THEN( "The faces of cells are finite" ) {
100  bool ok_finite = true;
101  for ( CvxCellComplex::Size c = 0; c < complex.nbCells(); ++c ) {
102  const auto faces = complex.cellFaces( c );
103  for ( auto f : faces )
104  ok_finite = ok_finite && ! complex.isInfinite( complex.faceCell( f ) );
105  }
106  REQUIRE( ok_finite );
107  }
108  THEN( "The opposite of faces of cells are infinite except two" ) {
109  int nb_finite = 0;
110  for ( CvxCellComplex::Size c = 0; c < complex.nbCells(); ++c ) {
111  const auto faces = complex.cellFaces( c );
112  for ( auto f : faces ) {
113  const auto opp_f = complex.opposite( f );
114  nb_finite += complex.isInfinite( complex.faceCell( opp_f ) ) ? 0 : 1;
115  }
116  }
117  REQUIRE( nb_finite == 2 );
118  }
119  }
120  }
121  GIVEN( "Given a degenerated polytope { (0,0), (3,-1), (9,-3), (-6,2) } " ) {
122  std::vector<Point> V
123  = { Point(0,0), Point(3,-1), Point(9,-3), Point(-6,2) };
124  WHEN( "Computing its lattice polytope" ){
125  const auto P = Helper::computeLatticePolytope( V, false, true );
126  CAPTURE( P );
127  THEN( "The polytope is valid and has 2 non trivial facets" ) {
128  REQUIRE( P.nbHalfSpaces() - 4 == 2 );
129  }
130  THEN( "The polytope contains 6 points" ) {
131  REQUIRE( P.count() == 6 );
132  }
133  THEN( "The polytope contains no interior points" ) {
134  REQUIRE( P.countInterior() == 0 );
135  }
136  }
137  WHEN( "Computing the vertices of its convex hull" ){
138  auto X = Helper::computeConvexHullVertices( V, false );
139  std::sort( X.begin(), X.end() );
140  CAPTURE( X );
141  THEN( "The polytope is a segment defined by two points" ) {
142  REQUIRE( X.size() == 2 );
143  REQUIRE( X[ 0 ] == Point(-6, 2) );
144  REQUIRE( X[ 1 ] == Point( 9,-3) );
145  }
146  }
147  }
148  GIVEN( "Given a degenerated simplex { (4,0), (7,2), (-5,-6) } " ) {
149  std::vector<Point> V
150  = { Point(4,0), Point(7,2), Point(-5,-6) };
151  WHEN( "Computing its lattice polytope" ){
152  const auto P = Helper::computeLatticePolytope( V, false, true );
153  CAPTURE( P );
154  THEN( "The polytope is valid and has 2 non trivial facets" ) {
155  REQUIRE( P.nbHalfSpaces() - 4 == 2 );
156  }
157  THEN( "The polytope contains 5 points" ) {
158  REQUIRE( P.count() == 5 );
159  }
160  THEN( "The polytope contains no interior points" ) {
161  REQUIRE( P.countInterior() == 0 );
162  }
163  }
164  WHEN( "Computing the vertices of its convex hull" ){
165  auto X = Helper::computeConvexHullVertices( V, false );
166  std::sort( X.begin(), X.end() );
167  CAPTURE( X );
168  THEN( "The polytope is a segment defined by two points" ) {
169  REQUIRE( X.size() == 2 );
170  REQUIRE( X[ 0 ] == Point(-5,-6) );
171  REQUIRE( X[ 1 ] == Point( 7, 2) );
172  }
173  }
174  }
175 }
176
178 // Functions for testing class ConvexityHelper in 3D.
180
181 SCENARIO( "ConvexityHelper< 3 > unit tests",
182  "[convexity_helper][3d]" )
183 {
184  typedef ConvexityHelper< 3 > Helper;
185  typedef Helper::Point Point;
189  typedef PolygonalSurface< Point > LatticePolySurf;
190  typedef ConvexCellComplex< Point > CvxCellComplex;
191  GIVEN( "Given an octahedron star { (0,0,0), (-2,0,0), (2,0,0), (0,-2,0), (0,2,0), (0,0,-2), (0,0,2) } " ) {
192  std::vector<Point> V
193  = { Point(0,0,0), Point(-2,0,0), Point(2,0,0), Point(0,-2,0), Point(0,2,0),
194  Point(0,0,-2), Point(0,0,2) };
195  WHEN( "Computing its lattice polytope" ){
196  const auto P = Helper::computeLatticePolytope( V, false, true );
197  CAPTURE( P );
198  THEN( "The polytope is valid and has 8 non trivial facets plus 12 edge constraints" ) {
199  REQUIRE( P.nbHalfSpaces() - 6 == 20 );
200  }
201  THEN( "The polytope is Minkowski summable" ) {
202  REQUIRE( P.canBeSummed() );
203  }
204  THEN( "The polytope contains the input points" ) {
205  REQUIRE( P.isInside( V[ 0 ] ) );
206  REQUIRE( P.isInside( V[ 1 ] ) );
207  REQUIRE( P.isInside( V[ 2 ] ) );
208  REQUIRE( P.isInside( V[ 3 ] ) );
209  REQUIRE( P.isInside( V[ 4 ] ) );
210  REQUIRE( P.isInside( V[ 5 ] ) );
211  REQUIRE( P.isInside( V[ 6 ] ) );
212  }
213  THEN( "The polytope contains 25 points" ) {
214  REQUIRE( P.count() == 25 );
215  }
216  THEN( "The interior of the polytope contains 7 points" ) {
217  REQUIRE( P.countInterior() == 7 );
218  }
219  THEN( "The boundary of the polytope contains 18 points" ) {
220  REQUIRE( P.countBoundary() == 18 );
221  }
222  }
223  WHEN( "Computing the boundary of its convex hull as a SurfaceMesh" ){
224  SMesh smesh;
225  bool ok = Helper::computeConvexHullBoundary( smesh, V, false );
226  CAPTURE( smesh );
227  THEN( "The surface mesh is valid and has 6 vertices, 12 edges and 8 faces" ) {
228  REQUIRE( ok );
229  REQUIRE( smesh.nbVertices() == 6 );
230  REQUIRE( smesh.nbEdges() == 12 );
231  REQUIRE( smesh.nbFaces() == 8 );
232  }
233  THEN( "The surface mesh has the topology of a sphere" ) {
234  REQUIRE( smesh.Euler() == 2 );
235  REQUIRE( smesh.computeManifoldBoundaryEdges().size() == 0 );
236  REQUIRE( smesh.computeNonManifoldEdges().size() == 0 );
237  }
238  }
239  WHEN( "Computing the boundary of its convex hull as a lattice PolygonalSurface" ){
240  LatticePolySurf lpsurf;
241  bool ok = Helper::computeConvexHullBoundary( lpsurf, V, false );
242  CAPTURE( lpsurf );
243  THEN( "The polygonal surface is valid and has 6 vertices, 12 edges and 8 faces" ) {
244  REQUIRE( ok );
245  REQUIRE( lpsurf.nbVertices() == 6 );
246  REQUIRE( lpsurf.nbEdges() == 12 );
247  REQUIRE( lpsurf.nbFaces() == 8 );
248  REQUIRE( lpsurf.nbArcs() == 24 );
249  }
250  THEN( "The polygonal surface has the topology of a sphere and no boundary" ) {
251  REQUIRE( lpsurf.Euler() == 2 );
252  REQUIRE( lpsurf.allBoundaryArcs().size() == 0 );
253  REQUIRE( lpsurf.allBoundaryVertices().size() == 0 );
254  }
255  }
256  WHEN( "Computing its convex hull as a ConvexCellComplex" ){
257  CvxCellComplex complex;
258  bool ok = Helper::computeConvexHullCellComplex( complex, V, false );
259  CAPTURE( complex );
260  THEN( "The convex cell complex is valid and has 6 vertices, 8 faces and 1 finite cell" ) {
261  REQUIRE( ok );
262  REQUIRE( complex.nbVertices() == 6 );
263  REQUIRE( complex.nbFaces() == 8 );
264  REQUIRE( complex.nbCells() == 1 );
265  }
266  }
267  WHEN( "Computing the vertices of its convex hull" ){
268  const auto X = Helper::computeConvexHullVertices( V, false );
269  CAPTURE( X );
270  THEN( "The polytope has 6 vertices" ) {
271  REQUIRE( X.size() == 6 );
272  }
273  }
274  }
275  GIVEN( "Given a cube with an additional outside vertex " ) {
276  std::vector<Point> V
277  = { Point(-10,-10,-10), Point(10,-10,-10), Point(-10,10,-10), Point(10,10,-10),
278  Point(-10,-10,10), Point(10,-10,10), Point(-10,10,10), Point(10,10,10),
279  Point(0,0,18) };
280  WHEN( "Computing its Delaunay cell complex" ){
281  CvxCellComplex complex;
282  bool ok = Helper::computeDelaunayCellComplex( complex, V, false );
283  CAPTURE( complex );
284  THEN( "The complex has 2 cells, 10 faces, 9 vertices" ) {
285  REQUIRE( ok );
286  REQUIRE( complex.nbCells() == 2 );
287  REQUIRE( complex.nbFaces() == 10 );
288  REQUIRE( complex.nbVertices() == 9 );
289  }
290  THEN( "The faces of cells are finite" ) {
291  bool ok_finite = true;
292  for ( CvxCellComplex::Size c = 0; c < complex.nbCells(); ++c ) {
293  const auto faces = complex.cellFaces( c );
294  for ( auto f : faces )
295  ok_finite = ok_finite && ! complex.isInfinite( complex.faceCell( f ) );
296  }
297  REQUIRE( ok_finite );
298  }
299  THEN( "The opposite of faces of cells are infinite except two" ) {
300  int nb_finite = 0;
301  for ( CvxCellComplex::Size c = 0; c < complex.nbCells(); ++c ) {
302  const auto faces = complex.cellFaces( c );
303  for ( auto f : faces ) {
304  const auto opp_f = complex.opposite( f );
305  nb_finite += complex.isInfinite( complex.faceCell( opp_f ) ) ? 0 : 1;
306  }
307  }
308  REQUIRE( nb_finite == 2 );
309  }
310  }
311  WHEN( "Computing the vertices of its convex hull" ){
312  const auto X = Helper::computeConvexHullVertices( V, false );
313  CAPTURE( X );
314  THEN( "The polytope has 9 vertices" ) {
315  REQUIRE( X.size() == 9 );
316  }
317  }
318  }
319  GIVEN( "Given a degenerated 1d polytope { (0,0,1), (3,-1,2), (9,-3,4), (-6,2,-1) } " ) {
320  std::vector<Point> V
321  = { Point(0,0,1), Point(3,-1,2), Point(9,-3,4), Point(-6,2,-1) };
322  WHEN( "Computing its lattice polytope" ){
323  const auto P = Helper::computeLatticePolytope( V, false, true );
324  CAPTURE( P );
325  THEN( "The polytope is valid and has 6 non trivial facets" ) {
326  REQUIRE( P.nbHalfSpaces() - 6 == 6 );
327  }
328  THEN( "The polytope contains 6 points" ) {
329  REQUIRE( P.count() == 6 );
330  }
331  THEN( "The polytope contains no interior points" ) {
332  REQUIRE( P.countInterior() == 0 );
333  }
334  }
335  WHEN( "Computing the vertices of its convex hull" ){
336  auto X = Helper::computeConvexHullVertices( V, false );
337  std::sort( X.begin(), X.end() );
338  CAPTURE( X );
339  THEN( "The polytope is a segment defined by two points" ) {
340  REQUIRE( X.size() == 2 );
341  REQUIRE( X[ 0 ] == Point(-6, 2,-1) );
342  REQUIRE( X[ 1 ] == Point( 9,-3, 4) );
343  }
344  }
345  }
346  GIVEN( "Given a degenerated 1d simplex { (1,0,-1), Point(4,-1,-2), Point(10,-3,-4) } " ) {
347  std::vector<Point> V
348  = { Point(1,0,-1), Point(4,-1,-2), Point(10,-3,-4) };
349  WHEN( "Computing its lattice polytope" ){
350  const auto P = Helper::computeLatticePolytope( V, false, true );
351  CAPTURE( P );
352  THEN( "The polytope is valid and has 6 non trivial facets" ) {
353  REQUIRE( P.nbHalfSpaces() - 6 == 6 );
354  }
355  THEN( "The polytope contains 4 points" ) {
356  REQUIRE( P.count() == 4 );
357  }
358  THEN( "The polytope contains no interior points" ) {
359  REQUIRE( P.countInterior() == 0 );
360  }
361  }
362  WHEN( "Computing the vertices of its convex hull" ){
363  auto X = Helper::computeConvexHullVertices( V, false );
364  std::sort( X.begin(), X.end() );
365  CAPTURE( X );
366  THEN( "The polytope is a segment defined by two points" ) {
367  REQUIRE( X.size() == 2 );
368  REQUIRE( X[ 0 ] == Point( 1, 0,-1) );
369  REQUIRE( X[ 1 ] == Point(10,-3,-4) );
370  }
371  }
372  }
373  GIVEN( "Given a degenerated 2d polytope { (2,1,0), (1,0,1), (1,2,1), (0,1,2), (0,3,2) } " ) {
374  std::vector<Point> V
375  = { Point(2,1,0), Point(1,0,1), Point(1,2,1), Point(0,1,2), Point(0,3,2) };
376  WHEN( "Computing its lattice polytope" ){
377  const auto P = Helper::computeLatticePolytope( V, false, true );
378  CAPTURE( P );
379  THEN( "The polytope is valid and has more than 6 non trivial facets" ) {
380  REQUIRE( P.nbHalfSpaces() - 6 == 6 );
381  }
382  THEN( "The polytope contains 7 points" ) {
383  REQUIRE( P.count() == 7 );
384  }
385  THEN( "The polytope contains no interior points" ) {
386  REQUIRE( P.countInterior() == 0 );
387  }
388  }
389  WHEN( "Computing the vertices of its convex hull" ){
390  auto X = Helper::computeConvexHullVertices( V, false );
391  std::sort( X.begin(), X.end() );
392  CAPTURE( X );
393  THEN( "The polytope is a quad" ) {
394  REQUIRE( X.size() == 4 );
395  REQUIRE( X[ 0 ] == Point( 0, 1, 2) );
396  REQUIRE( X[ 1 ] == Point( 0, 3, 2) );
397  REQUIRE( X[ 2 ] == Point( 1, 0, 1) );
398  REQUIRE( X[ 3 ] == Point( 2, 1, 0) );
399  }
400  }
401  }
402  GIVEN( "Given a degenerated 2d simplex { (2,1,0), (1,0,1), (1,5,1), (0,3,2) } " ) {
403  std::vector<Point> V
404  = { Point(2,1,0), Point(1,0,1), Point(1,5,1), Point(0,3,2) };
405  WHEN( "Computing its lattice polytope" ){
406  const auto P = Helper::computeLatticePolytope( V, false, true );
407  CAPTURE( P );
408  THEN( "The polytope is valid and has more than 6 non trivial facets" ) {
409  REQUIRE( P.nbHalfSpaces() - 6 == 6 );
410  }
411  THEN( "The polytope contains 8 points" ) {
412  REQUIRE( P.count() == 8 );
413  }
414  THEN( "The polytope contains no interior points" ) {
415  REQUIRE( P.countInterior() == 0 );
416  }
417  }
418  WHEN( "Computing the vertices of its convex hull" ){
419  auto X = Helper::computeConvexHullVertices( V, false );
420  std::sort( X.begin(), X.end() );
421  CAPTURE( X );
422  THEN( "The polytope is a quad" ) {
423  REQUIRE( X.size() == 4 );
424  REQUIRE( X[ 0 ] == Point( 0, 3, 2) );
425  REQUIRE( X[ 1 ] == Point( 1, 0, 1) );
426  REQUIRE( X[ 2 ] == Point( 1, 5, 1) );
427  REQUIRE( X[ 3 ] == Point( 2, 1, 0) );
428  }
429  }
430  }
431 }
432
433
435 // Functions for testing class ConvexityHelper in 3D.
437
438 SCENARIO( "ConvexityHelper< 3 > triangle tests",
439  "[convexity_helper][triangle][3d]" )
440 {
441  typedef ConvexityHelper< 3 > Helper;
442  typedef Helper::Point Point;
443  typedef Helper::Vector Vector;
444  GIVEN( "Given non degenerated triplets of points" ) {
445  Helper::LatticePolytope P;
446  Helper::LatticePolytope T;
447  Point a, b, c;
448  Vector n;
449  int nb_total = 0;
450  int nb_ok = 0;
451  int nbi_ok = 0;
452  int nb_P = 0, nb_T = 0, nbi_P = 0, nbi_T = 0;
453  for ( int i = 0; i < 20; i++ )
454  {
455  do {
456  a = Point( rand() % 10, rand() % 10, rand() % 10 );
457  b = Point( rand() % 10, rand() % 10, rand() % 10 );
458  c = Point( rand() % 10, rand() % 10, rand() % 10 );
459  n = (b-a).crossProduct(c-a);
460  } while ( n == Vector::zero );
461  std::vector<Point> V = { a, b, c };
462  P = Helper::computeLatticePolytope( V, false, false );
463  T = Helper::compute3DTriangle( a, b, c );
464  nb_P = P.count();
465  nb_T = T.count();
466  nbi_P = P.countInterior();
467  nbi_T = T.countInterior();
468  nb_total += 1;
469  nb_ok += ( nb_P == nb_T ) ? 1 : 0;
470  nbi_ok += ( nbi_P == 0 && nbi_T == 0 ) ? 1 : 0;
471  if ( ( nb_ok != nb_total ) || ( nbi_ok != nb_total ) ) break;
472  }
473  CAPTURE( a ); CAPTURE( b ); CAPTURE( c ); CAPTURE( n );
474  CAPTURE( P ); CAPTURE( T );
475  CAPTURE( nb_P ); CAPTURE( nb_T ); CAPTURE( nbi_P ); CAPTURE( nbi_T );
476  WHEN( "Computing their tightiest polytope or triangle" ) {
477  THEN( "They have the same number of inside points" ) {
478  REQUIRE( nb_ok == nb_total );
479  }
480  THEN( "They do not have interior points" ) {
481  REQUIRE( nbi_ok == nb_total );
482  }
483  }
484  }
485
486  GIVEN( "Given non degenerated triplets of points" ) {
487  typedef Helper::LatticePolytope::UnitSegment UnitSegment;
488  Helper::LatticePolytope P;
489  Helper::LatticePolytope T;
490  Point a, b, c;
491  Vector n;
492  int nb_total = 0;
493  int nb_ok = 0;
494  int nbi_ok = 0;
495  int nb_P = 0, nb_T = 0, nbi_P = 0, nbi_T = 0;
496  for ( int i = 0; i < 20; i++ )
497  {
498  do {
499  a = Point( rand() % 10, rand() % 10, rand() % 10 );
500  b = Point( rand() % 10, rand() % 10, rand() % 10 );
501  c = Point( rand() % 10, rand() % 10, rand() % 10 );
502  n = (b-a).crossProduct(c-a);
503  } while ( n == Vector::zero );
504  std::vector<Point> V = { a, b, c };
505  P = Helper::computeLatticePolytope( V, false, true );
506  T = Helper::compute3DTriangle( a, b, c, true );
507  for ( Dimension k = 0; k < 3; k++ )
508  {
509  P += UnitSegment( k );
510  T += UnitSegment( k );
511  }
512  nb_P = P.count();
513  nb_T = T.count();
514  nbi_P = P.countInterior();
515  nbi_T = T.countInterior();
516  nb_total += 1;
517  nb_ok += ( nb_P == nb_T ) ? 1 : 0;
518  nbi_ok += ( nbi_P == nbi_T ) ? 1 : 0;
519  if ( ( nb_ok != nb_total ) || ( nbi_ok != nb_total ) ) break;
520  }
521  CAPTURE( a ); CAPTURE( b ); CAPTURE( c ); CAPTURE( n );
522  CAPTURE( P ); CAPTURE( T );
523  CAPTURE( nb_P ); CAPTURE( nb_T ); CAPTURE( nbi_P ); CAPTURE( nbi_T );
524  WHEN( "Computing their tightiest polytope or triangle, dilated by a cube" ) {
525  THEN( "They have the same number of inside points" ) {
526  REQUIRE( nb_ok == nb_total );
527  }
528  THEN( "They have the same number of interior points" ) {
529  REQUIRE( nbi_ok == nb_total );
530  }
531  }
532  }
533 }
534
535
536 SCENARIO( "ConvexityHelper< 3 > degenerated triangle tests",
537  "[convexity_helper][triangle][degenerate][3d]" )
538 {
539  typedef ConvexityHelper< 3 > Helper;
540  typedef Helper::Point Point;
541  typedef Helper::Vector Vector;
542  GIVEN( "Given degenerated triplets of points" ) {
543  Helper::LatticePolytope P;
544  Helper::LatticePolytope T;
545  Point a, b, c;
546  Vector n;
547  int nb_total = 0;
548  int nb_ok = 0;
549  int nbi_ok = 0;
550  int nb_P = 0, nb_T = 0, nbi_P = 0, nbi_T = 0;
551  for ( int i = 0; i < 20; i++ )
552  {
553  do {
554  a = Point( rand() % 10, rand() % 10, rand() % 10 );
555  b = Point( rand() % 10, rand() % 10, rand() % 10 );
556  c = Point( rand() % 10, rand() % 10, rand() % 10 );
557  n = (b-a).crossProduct(c-a);
558  } while ( n != Vector::zero );
559  std::vector<Point> V = { a, b, c };
560  P = Helper::computeLatticePolytope( V, true, false );
561  T = Helper::compute3DTriangle( a, b, c );
562  nb_P = P.count();
563  nb_T = T.count();
564  nbi_P = P.countInterior();
565  nbi_T = T.countInterior();
566  nb_total += 1;
567  nb_ok += ( nb_P == nb_T ) ? 1 : 0;
568  nbi_ok += ( nbi_P == 0 && nbi_T == 0 ) ? 1 : 0;
569  if ( ( nb_ok != nb_total ) || ( nbi_ok != nb_total ) ) break;
570  }
571  CAPTURE( a ); CAPTURE( b ); CAPTURE( c ); CAPTURE( n );
572  CAPTURE( P ); CAPTURE( T );
573  CAPTURE( nb_P ); CAPTURE( nb_T ); CAPTURE( nbi_P ); CAPTURE( nbi_T );
574  WHEN( "Computing their tightiest polytope or triangle" ) {
575  THEN( "They have the same number of inside points" ) {
576  REQUIRE( nb_ok == nb_total );
577  }
578  THEN( "They do not have interior points" ) {
579  REQUIRE( nbi_ok == nb_total );
580  }
581  }
582  }
583
584  GIVEN( "Given degenerated triplets of points" ) {
585  typedef Helper::LatticePolytope::UnitSegment UnitSegment;
586  Helper::LatticePolytope P;
587  Helper::LatticePolytope T;
588  Point a, b, c;
589  Vector n;
590  int nb_total = 0;
591  int nb_ok = 0;
592  int nbi_ok = 0;
593  int nb_P = 0, nb_T = 0, nbi_P = 0, nbi_T = 0;
594  for ( int i = 0; i < 20; i++ )
595  {
596  do {
597  a = Point( rand() % 10, rand() % 10, rand() % 10 );
598  b = Point( rand() % 10, rand() % 10, rand() % 10 );
599  c = Point( rand() % 10, rand() % 10, rand() % 10 );
600  n = (b-a).crossProduct(c-a);
601  } while ( n != Vector::zero );
602  std::vector<Point> V = { a, b, c };
603  P = Helper::computeLatticePolytope( V, true, true );
604  T = Helper::compute3DTriangle( a, b, c, true );
605  for ( Dimension k = 0; k < 3; k++ )
606  {
607  P += UnitSegment( k );
608  T += UnitSegment( k );
609  }
610  nb_P = P.count();
611  nb_T = T.count();
612  nbi_P = P.countInterior();
613  nbi_T = T.countInterior();
614  nb_total += 1;
615  nb_ok += ( nb_P == nb_T ) ? 1 : 0;
616  nbi_ok += ( nbi_P == nbi_T ) ? 1 : 0;
617  if ( ( nb_ok != nb_total ) || ( nbi_ok != nb_total ) ) break;
618  }
619  CAPTURE( a ); CAPTURE( b ); CAPTURE( c ); CAPTURE( n );
620  CAPTURE( P ); CAPTURE( T );
621  CAPTURE( nb_P ); CAPTURE( nb_T ); CAPTURE( nbi_P ); CAPTURE( nbi_T );
622  WHEN( "Computing their tightiest polytope or triangle, dilated by a cube" ) {
623  THEN( "They have the same number of inside points" ) {
624  REQUIRE( nb_ok == nb_total );
625  }
626  THEN( "They have the same number of interior points" ) {
627  REQUIRE( nbi_ok == nb_total );
628  }
629  }
630  }
631 }
Aim: Represents a polygon mesh, i.e. a 2-dimensional combinatorial surface whose faces are (topologic...
Space::RealVector RealVector
DigitalPlane::Point Vector
SurfaceMesh< RealPoint, RealVector > SMesh
DGtal is the top-level namespace which contains all DGtal functions and types.
auto crossProduct(PointVector< 3, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< 3, RightEuclideanRing, RightContainer > const &rhs) -> decltype(DGtal::constructFromArithmeticConversion(lhs, rhs))
Cross product of two 3D Points/Vectors.
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: represents a d-dimensional complex in a d-dimensional space with the following properties and re...
Aim: Provides a set of functions to facilitate the computation of convex hulls and polytopes,...
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
long Euler() const
Definition: SurfaceMesh.h:302
Size nbFaces() const
Definition: SurfaceMesh.h:296
Size nbEdges() const
Definition: SurfaceMesh.h:292
Edges computeManifoldBoundaryEdges() const
Size nbVertices() const
Definition: SurfaceMesh.h:288
Edges computeNonManifoldEdges() const
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
SCENARIO("ConvexityHelper< 2 > unit tests", "[convexity_helper][lattice_polytope][2d]")
GIVEN("A cubical complex with random 3-cells")
HalfEdgeDataStructure::Size Size
REQUIRE(domain.isInside(aPoint))
PointVector< 3, double > RealPoint