DGtal  1.5.beta
testQuickHull.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/tools/QuickHull.h"
37 #include "DGtalCatch.h"
39 
40 using namespace std;
41 using namespace DGtal;
42 
43 template <typename Point>
44 std::vector< Point >
45 randomPointsInBall( int nb, int R )
46 {
47  std::vector< Point > V;
48  Point c = Point::diagonal( R );
49  double R2 = (double) R * (double) R;
50  for ( int i = 0; i < nb; ) {
51  Point p;
52  for ( DGtal::Dimension k = 0; k < Point::dimension; ++k )
53  p[ k ] = rand() % (2*R);
54  if ( ( p - c ).squaredNorm() < R2 ) { V.push_back( p ); i++; }
55  }
56  return V;
57 }
58 
60 // Functions for testing class QuickHull in 2D.
62 
63 SCENARIO( "QuickHull< ConvexHullIntegralKernel< 2 > > unit tests", "[quickhull][integral_kernel][2d]" )
64 {
65  typedef ConvexHullIntegralKernel< 2 > QHKernel;
66  typedef QuickHull< QHKernel > QHull;
67  typedef SpaceND< 2, int > Space;
68  typedef Space::Point Point;
69  typedef QHull::Index Index;
70  typedef QHull::IndexRange IndexRange;
71 
72 
73  GIVEN( "Given a star { (0,0), (-4,-1), (-3,5), (7,3), (5, -2) } " ) {
74  std::vector<Point> V
75  = { Point(0,0), Point(-4,-1), Point(-3,5), Point(7,3), Point(5, -2) };
76  QHull hull;
77  hull.setInput( V, false );
78  hull.computeConvexHull();
79  THEN( "The convex hull is valid and contains every point" ) {
80  REQUIRE( hull.check() );
81  }
82  THEN( "Its convex hull has 4 vertices" ) {
83  REQUIRE( hull.nbVertices() == 4 );
84  }
85  THEN( "Its convex hull has 4 facets" ) {
86  REQUIRE( hull.nbFacets() == 4 );
87  }
88  THEN( "Its facets form a linked list" ) {
89  std::vector< IndexRange > facets;
90  hull.getFacetVertices( facets );
91  std::vector< Index > next( hull.nbVertices(), (Index) -1 );
92  Index nb_zero_next = hull.nbVertices();
93  Index nb_two_next = 0;
94  for ( auto f : facets ) {
95  if ( next[ f[ 0 ] ] != (Index) -1 ) nb_two_next += 1;
96  else {
97  next[ f[ 0 ] ] = f[ 1 ];
98  nb_zero_next -= 1;
99  }
100  }
101  REQUIRE( nb_zero_next == 0 );
102  REQUIRE( nb_two_next == 0 );
103  }
104  }
105  GIVEN( "Given 100 random point in a ball of radius 50 " ) {
106  std::vector<Point> V = randomPointsInBall< Point >( 100, 50 );
107  QHull hull;
108  hull.setInput( V, false );
109  hull.computeConvexHull();
110  THEN( "The convex hull is valid and contains every point" ) {
111  REQUIRE( hull.check() );
112  }
113  THEN( "Its convex hull has the same number of vertices and facets" ) {
114  REQUIRE( hull.nbVertices() == hull.nbFacets() );
115  }
116  THEN( "Its convex hull has much fewer vertices than input points" ) {
117  REQUIRE( 2*hull.nbVertices() <= hull.nbPoints() );
118  }
119  THEN( "Its facets form a linked list" ) {
120  std::vector< IndexRange > facets;
121  hull.getFacetVertices( facets );
122  std::vector< Index > next( hull.nbVertices(), (Index) -1 );
123  Index nb_zero_next = hull.nbVertices();
124  Index nb_two_next = 0;
125  for ( auto f : facets ) {
126  if ( next[ f[ 0 ] ] != (Index) -1 ) nb_two_next += 1;
127  else {
128  next[ f[ 0 ] ] = f[ 1 ];
129  nb_zero_next -= 1;
130  }
131  }
132  REQUIRE( nb_zero_next == 0 );
133  REQUIRE( nb_two_next == 0 );
134  }
135  }
136 }
137 
138 
140 // Functions for testing class QuickHull in 3D.
142 
143 SCENARIO( "QuickHull< ConvexHullIntegralKernel< 3 > > unit tests", "[quickhull][integral_kernel][3d]" )
144 {
145  typedef ConvexHullIntegralKernel< 3 > QHKernel;
146  typedef QuickHull< QHKernel > QHull;
147  typedef SpaceND< 3, int > Space;
148  typedef Space::Point Point;
149  typedef QHull::IndexRange IndexRange;
150 
151 
152  GIVEN( "Given an octahedron" ) {
153  const int R = 5;
154  std::vector<Point> V = { Point( 0,0,0 ) };
155  for ( Dimension k = 0; k < 3; ++k ) {
156  V.push_back( Point::base( k, R ) );
157  V.push_back( Point::base( k, -R ) );
158  }
159  QHull hull;
160  hull.setInput( V, false );
161  hull.setInitialSimplex( IndexRange { 0, 1, 3, 5 } );
162  hull.computeConvexHull();
163  THEN( "The convex hull is valid and contains every point" ) {
164  REQUIRE( hull.check() );
165  }
166  THEN( "Its convex hull has 6 vertices" ) {
167  REQUIRE( hull.nbVertices() == 6 );
168  }
169  THEN( "Its convex hull has 8 facets" ) {
170  REQUIRE( hull.nbFacets() == 8 );
171  }
172  }
173  GIVEN( "Given 100 random point in a ball of radius 50 " ) {
174  std::vector<Point> V = randomPointsInBall< Point >( 100, 50 );
175  QHull hull;
176  hull.setInput( V, false );
177  hull.computeConvexHull();
178  THEN( "The convex hull is valid and contains every point" ) {
179  REQUIRE( hull.check() );
180  }
181  THEN( "Its convex hull has more facets than vertices" ) {
182  REQUIRE( hull.nbVertices() < hull.nbFacets() );
183  }
184  THEN( "Its convex hull has fewer vertices than input points" ) {
185  REQUIRE( hull.nbVertices() < hull.nbPoints() );
186  }
187  }
188 }
189 
190 
192 // Functions for testing class QuickHull in 4D.
194 
195 SCENARIO( "QuickHull< ConvexHullIntegralKernel< 4 > > unit tests", "[quickhull][integral_kernel][4d]" )
196 {
197  typedef ConvexHullIntegralKernel< 4 > QHKernel;
198  typedef QuickHull< QHKernel > QHull;
199  typedef SpaceND< 4, int > Space;
200  typedef Space::Point Point;
201 
202 
203  GIVEN( "Given an octahedron" ) {
204  const int R = 5;
205  std::vector<Point> V = { Point( 0,0,0,0 ) };
206  for ( Dimension k = 0; k < 4; ++k ) {
207  V.push_back( Point::base( k, R ) );
208  V.push_back( Point::base( k, -R ) );
209  }
210  QHull hull;
211  hull.setInput( V, false );
212  hull.computeConvexHull();
213  THEN( "The convex hull is valid and contains every point" ) {
214  REQUIRE( hull.check() );
215  }
216  THEN( "Its convex hull has 8 vertices" ) {
217  REQUIRE( hull.nbVertices() == 8 );
218  }
219  THEN( "Its convex hull has 16 facets" ) {
220  REQUIRE( hull.nbFacets() == 16 );
221  }
222  }
223  GIVEN( "Given 100 random point in a ball of radius 50 " ) {
224  std::vector<Point> V = randomPointsInBall< Point >( 100, 50 );
225  QHull hull;
226  hull.setInput( V, false );
227  hull.computeConvexHull();
228  THEN( "The convex hull is valid and contains every point" ) {
229  REQUIRE( hull.check() );
230  }
231  THEN( "Its convex hull has fewer vertices than input points" ) {
232  REQUIRE( hull.nbVertices() < hull.nbPoints() );
233  }
234  }
235 }
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: a geometric kernel to compute the convex hull of digital points with integer-only arithmetic.
Aim: Implements the quickhull algorithm by Barber et al. , a famous arbitrary dimensional convex hull...
Definition: QuickHull.h:140
MyPointD Point
Definition: testClone2.cpp:383
GIVEN("A cubical complex with random 3-cells")
std::vector< Point > randomPointsInBall(int nb, int R)
SCENARIO("QuickHull< ConvexHullIntegralKernel< 2 > > unit tests", "[quickhull][integral_kernel][2d]")
REQUIRE(domain.isInside(aPoint))