DGtal  1.4.beta
testConvexityHelper.cpp
Go to the documentation of this file.
1 
30 #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 ( auto 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 ( auto 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 ( auto 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 ( auto 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 }
DGtal::SurfaceMesh::Euler
long Euler() const
Definition: SurfaceMesh.h:294
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
REQUIRE
REQUIRE(domain.isInside(aPoint))
DGtal::SurfaceMesh::nbVertices
Size nbVertices() const
Definition: SurfaceMesh.h:280
DGtal::ConvexCellComplex
Aim: represents a d-dimensional complex in a d-dimensional space with the following properties and re...
Definition: ConvexCellComplex.h:85
CAPTURE
CAPTURE(thicknessHV)
DGtal::SurfaceMesh::computeNonManifoldEdges
Edges computeNonManifoldEdges() const
DGtal::SurfaceMesh::nbEdges
Size nbEdges() const
Definition: SurfaceMesh.h:284
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
RealVector
Space::RealVector RealVector
Definition: digitalPolyhedronBuilder3D.cpp:82
DGtal::ConvexityHelper
Aim: Provides a set of functions to facilitate the computation of convex hulls and polytopes,...
Definition: ConvexityHelper.h:164
SCENARIO
SCENARIO("ConvexityHelper< 2 > unit tests", "[convexity_helper][lattice_polytope][2d]")
Definition: testConvexityHelper.cpp:49
SMesh
SurfaceMesh< RealPoint, RealVector > SMesh
Definition: fullConvexitySphereGeodesics.cpp:115
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition: testCubicalComplex.cpp:70
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::PolygonalSurface
Aim: Represents a polygon mesh, i.e. a 2-dimensional combinatorial surface whose faces are (topologic...
Definition: PolygonalSurface.h:86
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288