DGtal 1.3.0
Searching...
No Matches
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
41using namespace std;
42using namespace DGtal;
43
44
46// Functions for testing class ConvexityHelper in 2D.
48
49SCENARIO( "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
181SCENARIO( "ConvexityHelper< 3 > unit tests",
182 "[convexity_helper][3d]" )
183{
184 typedef ConvexityHelper< 3 > Helper;
185 typedef Helper::Point Point;
186 typedef Helper::RealPoint RealPoint;
187 typedef Helper::RealVector RealVector;
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}
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
Aim: Represents a polygon mesh, i.e. a 2-dimensional combinatorial surface whose faces are (topologic...
SurfaceMesh< RealPoint, RealVector > SMesh
DGtal is the top-level namespace which contains all DGtal functions and types.
STL namespace.
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:294
Size nbFaces() const
Definition: SurfaceMesh.h:288
Size nbEdges() const
Definition: SurfaceMesh.h:284
Edges computeManifoldBoundaryEdges() const
Size nbVertices() const
Definition: SurfaceMesh.h:280
Edges computeNonManifoldEdges() const
Z2i::RealPoint RealPoint
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
REQUIRE(domain.isInside(aPoint))
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")