DGtal  1.5.beta
testCellGeometry.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 CellGeometry.
49 
50 SCENARIO( "CellGeometry< Z2 > segment tests", "[cell_geometry][2d][segment]" )
51 {
52  typedef KhalimskySpaceND<2,int> KSpace;
53  typedef KSpace::Point Point;
54  typedef CellGeometry< KSpace > CGeometry;
55  typedef DigitalConvexity< KSpace > DConvexity;
56 
57  KSpace K;
58  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
59  DConvexity dconv( K );
60  GIVEN( "two points (1,1), Point(3,-2)" ) {
61  CGeometry geometry( K, 0, 2, false );
62  geometry.addCellsTouchingSegment( Point(1,1), Point(3,-2) );
63  THEN( "Its cell geometry contains 2 0-cells, 11 1-cells, 10 2-cells" ) {
64  auto C0 = geometry.getKPoints( 0 );
65  auto C1 = geometry.getKPoints( 1 );
66  auto C2 = geometry.getKPoints( 2 );
67  CAPTURE( C0 );
68  CAPTURE( C1 );
69  CAPTURE( C2 );
70  REQUIRE( C0.size() == 2 );
71  REQUIRE( C1.size() == 11 );
72  REQUIRE( C2.size() == 10 );
73  }
74  }
75  WHEN( "Computing random segments in domain (-5,-5)-(10,10)." ) {
76  unsigned int nb = 20;
77  for ( unsigned int i = 0; i < nb; ++i )
78  {
79  Point a( rand() % 13 - 4, rand() % 13 - 4 );
80  Point b( rand() % 13 - 4, rand() % 13 - 4 );
81  if ( a == b ) continue;
82  CGeometry segm_geometry( K, 0, 2, false );
83  segm_geometry.addCellsTouchingSegment( a, b );
84  CGeometry splx_geometry( K, 0, 2, false );
85  auto splx = dconv.makeSimplex( { a, b } );
86  splx_geometry.addCellsTouchingPolytope( splx );
87  auto segm_0 = segm_geometry.getKPoints( 0 );
88  auto splx_0 = splx_geometry.getKPoints( 0 );
89  auto segm_1 = segm_geometry.getKPoints( 1 );
90  auto splx_1 = splx_geometry.getKPoints( 1 );
91  auto segm_2 = segm_geometry.getKPoints( 2 );
92  auto splx_2 = splx_geometry.getKPoints( 2 );
93  THEN( "Generic addCellsTouchingPolytope and specialized addCellsTouchingSegment should provide the same result" ) {
94  REQUIRE( segm_0.size() == splx_0.size() );
95  REQUIRE( segm_1.size() == splx_1.size() );
96  REQUIRE( segm_2.size() == splx_2.size() );
97  std::sort( segm_0.begin(), segm_0.end() );
98  std::sort( segm_1.begin(), segm_1.end() );
99  std::sort( segm_2.begin(), segm_2.end() );
100  std::sort( splx_0.begin(), splx_0.end() );
101  std::sort( splx_1.begin(), splx_1.end() );
102  std::sort( splx_2.begin(), splx_2.end() );
103  CAPTURE( segm_0 );
104  CAPTURE( segm_1 );
105  CAPTURE( segm_2 );
106  CAPTURE( splx_0 );
107  CAPTURE( splx_1 );
108  CAPTURE( splx_2 );
109  REQUIRE( std::equal( segm_0.cbegin(), segm_0.cend(), splx_0.cbegin() ) );
110  REQUIRE( std::equal( segm_1.cbegin(), segm_1.cend(), splx_1.cbegin() ) );
111  REQUIRE( std::equal( segm_2.cbegin(), segm_2.cend(), splx_2.cbegin() ) );
112  }
113  }
114  }
115 }
116 
117 SCENARIO( "CellGeometry< Z3 > segment tests", "[cell_geometry][3d][segment]" )
118 {
119  typedef KhalimskySpaceND<3,int> KSpace;
120  typedef KSpace::Point Point;
121  typedef CellGeometry< KSpace > CGeometry;
122  typedef DigitalConvexity< KSpace > DConvexity;
123 
124  KSpace K;
125  K.init( Point( -5, -5, -5 ), Point( 10, 10, 10 ), true );
126  DConvexity dconv( K );
127  WHEN( "Computing random segments in domain (-5,-5,-5)-(10,10,10)." ) {
128  unsigned int nb = 20;
129  for ( unsigned int i = 0; i < nb; ++i )
130  {
131  Point a( rand() % 13 - 4, rand() % 13 - 4, rand() % 13 - 4 );
132  Point b( rand() % 13 - 4, rand() % 13 - 4, rand() % 13 - 4 );
133  if ( a == b ) continue;
134  CGeometry segm_geometry( K, 0, 3, false );
135  segm_geometry.addCellsTouchingSegment( a, b );
136  CGeometry splx_geometry( K, 0, 3, false );
137  auto splx = dconv.makeSimplex( { a, b } );
138  splx_geometry.addCellsTouchingPolytope( splx );
139  auto segm_0 = segm_geometry.getKPoints( 0 );
140  auto splx_0 = splx_geometry.getKPoints( 0 );
141  auto segm_1 = segm_geometry.getKPoints( 1 );
142  auto splx_1 = splx_geometry.getKPoints( 1 );
143  auto segm_2 = segm_geometry.getKPoints( 2 );
144  auto splx_2 = splx_geometry.getKPoints( 2 );
145  auto segm_3 = segm_geometry.getKPoints( 3 );
146  auto splx_3 = splx_geometry.getKPoints( 3 );
147  THEN( "Generic addCellsTouchingPolytope and specialized addCellsTouchingSegment should provide the same result" ) {
148  REQUIRE( segm_0.size() == splx_0.size() );
149  REQUIRE( segm_1.size() == splx_1.size() );
150  REQUIRE( segm_2.size() == splx_2.size() );
151  REQUIRE( segm_3.size() == splx_3.size() );
152  std::sort( segm_0.begin(), segm_0.end() );
153  std::sort( segm_1.begin(), segm_1.end() );
154  std::sort( segm_2.begin(), segm_2.end() );
155  std::sort( segm_3.begin(), segm_3.end() );
156  std::sort( splx_0.begin(), splx_0.end() );
157  std::sort( splx_1.begin(), splx_1.end() );
158  std::sort( splx_2.begin(), splx_2.end() );
159  std::sort( splx_3.begin(), splx_3.end() );
160  CAPTURE( segm_0 );
161  CAPTURE( segm_1 );
162  CAPTURE( segm_2 );
163  CAPTURE( segm_3 );
164  CAPTURE( splx_0 );
165  CAPTURE( splx_1 );
166  CAPTURE( splx_2 );
167  CAPTURE( splx_3 );
168  REQUIRE( std::equal( segm_0.cbegin(), segm_0.cend(), splx_0.cbegin() ) );
169  REQUIRE( std::equal( segm_1.cbegin(), segm_1.cend(), splx_1.cbegin() ) );
170  REQUIRE( std::equal( segm_2.cbegin(), segm_2.cend(), splx_2.cbegin() ) );
171  REQUIRE( std::equal( segm_3.cbegin(), segm_3.cend(), splx_3.cbegin() ) );
172  }
173  }
174  }
175 }
176 
177 SCENARIO( "CellGeometry< Z2 > unit tests", "[cell_geometry][2d]" )
178 {
179  typedef KhalimskySpaceND<2,int> KSpace;
180  typedef KSpace::Point Point;
181  typedef CellGeometry< KSpace > CGeometry;
182 
183  KSpace K;
184  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
185  GIVEN( "Some points (0,0), Point(2,1), Point(1,3)" ) {
186  std::vector< Point > V = { Point(0,0), Point(2,1), Point(1,3) };
187  CGeometry geometry( K, 0, 2, false );
188  geometry.addCellsTouchingPoints( V.begin(), V.end() );
189  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
190  REQUIRE( V.size() == geometry.computeNbCells( 0 ) );
191  REQUIRE( V.size() < geometry.computeNbCells( 1 ) );
192  REQUIRE( V.size() < geometry.computeNbCells( 2 ) );
193  }
194  THEN( "Its cells form an open complex with euler characteristic 3" ) {
195  REQUIRE( geometry.computeEuler() == 3 );
196  }
197  }
198 }
199 
200 SCENARIO( "CellGeometry< Z3 > unit tests", "[cell_geometry][3d]" )
201 {
202  typedef KhalimskySpaceND<3,int> KSpace;
203  typedef KSpace::Space Space;
204  typedef KSpace::Point Point;
205  typedef CellGeometry< KSpace > CGeometry;
206  typedef HyperRectDomain<Space> Domain;
207 
208  KSpace K;
209  const int N = 10;
210  K.init( Point( -1, -1, -1 ), Point( N, N, N ), true );
211  Domain D( Point( 0, 0, 0 ), Point( N-1, N-1, N-1 ) );
212 
213  GIVEN( "Some a block of points (0,0,0)-(N-1,N-1,N-1)" ) {
214  CGeometry geometry( K, 0, 3, false );
215  geometry.addCellsTouchingPoints( D.begin(), D.end() );
216  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
217  REQUIRE( D.size() == geometry.computeNbCells( 0 ) );
218  REQUIRE( D.size() < geometry.computeNbCells( 1 ) );
219  REQUIRE( D.size() < geometry.computeNbCells( 2 ) );
220  REQUIRE( D.size() < geometry.computeNbCells( 3 ) );
221  }
222  THEN( "Its cells form an open complex with euler characteristic -1" ) {
223  REQUIRE( geometry.computeEuler() == -1 );
224  }
225  }
226 
227  GIVEN( "Some a block of points (0,0,0)-(N-1,N-1,N-1)" ) {
228  CGeometry geometry( K, 0, 2, false );
229  geometry.addCellsTouchingPoints( D.begin(), D.end() );
230  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
231  REQUIRE( D.size() == geometry.computeNbCells( 0 ) );
232  REQUIRE( D.size() < geometry.computeNbCells( 1 ) );
233  REQUIRE( D.size() < geometry.computeNbCells( 2 ) );
234  }
235  }
236 }
237 
238 SCENARIO( "CellGeometry< Z2 > intersections", "[cell_geometry][2d]" )
239 {
240  typedef KhalimskySpaceND<2,int> KSpace;
241  typedef KSpace::Point Point;
242  typedef KSpace::Space Space;
243  typedef CellGeometry< KSpace > CGeometry;
244  typedef BoundedLatticePolytope< Space > Polytope;
245 
246  GIVEN( "A simplex P={ Point(0,0), Point(4,2), Point(-1,4) }" ) {
247  KSpace K;
248  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
249  std::vector< Point > V = { Point(0,0), Point(4,2), Point(-1,4) };
250  Polytope P( V.begin(), V.end() );
251  CGeometry intersected_cover( K, 0, 2, false );
252  intersected_cover.addCellsTouchingPolytope( P );
253  CGeometry touched_cover( K, 0, 2, false );
254  touched_cover.addCellsTouchingPoints( V.begin(), V.end() );
255  CGeometry touched_points_cover( K, 0, 2, false );
256  touched_points_cover.addCellsTouchingPolytopePoints( P );
257  // trace.info() << "Polytope P=" << P << std::endl;
258  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
259  REQUIRE( intersected_cover.computeEuler() == 1 );
260  }
261  THEN( "Its convex hull intersects more cells than its vertices touch." ) {
262  REQUIRE( touched_cover.computeNbCells( 0 )
263  < intersected_cover.computeNbCells( 0 ) );
264  REQUIRE( touched_cover.computeNbCells( 1 )
265  < intersected_cover.computeNbCells( 1 ) );
266  REQUIRE( touched_cover.computeNbCells( 2 )
267  < intersected_cover.computeNbCells( 2 ) );
268  }
269  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
270  REQUIRE( touched_points_cover.computeNbCells( 0 )
271  <= intersected_cover.computeNbCells( 0 ) );
272  REQUIRE( touched_points_cover.computeNbCells( 1 )
273  <= intersected_cover.computeNbCells( 1 ) );
274  REQUIRE( touched_points_cover.computeNbCells( 2 )
275  <= intersected_cover.computeNbCells( 2 ) );
276  }
277  }
278 }
279 
280 SCENARIO( "CellGeometry< Z3 > intersections", "[cell_geometry][3d]" )
281 {
282  typedef KhalimskySpaceND<3,int> KSpace;
283  typedef KSpace::Point Point;
284  typedef KSpace::Space Space;
285  typedef CellGeometry< KSpace > CGeometry;
286  typedef BoundedLatticePolytope< Space > Polytope;
287 
288  GIVEN( "A simplex P={ Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) }" ) {
289  KSpace K;
290  K.init( Point( -5, -5, -5 ), Point( 10, 10, 10 ), true );
291  CGeometry intersected_cover( K, 0, 3, false );
292  Polytope P = { Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) };
293  intersected_cover.addCellsTouchingPolytope( P );
294  std::vector< Point > V = { Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) };
295  CGeometry touched_cover( K, 0, 3, false );
296  touched_cover.addCellsTouchingPoints( V.begin(), V.end() );
297  CGeometry touched_points_cover( K, 0, 3, false );
298  touched_points_cover.addCellsTouchingPolytopePoints( P );
299  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
300  REQUIRE( intersected_cover.computeEuler() == -1 );
301  }
302  THEN( "Its convex hull intersects more cells than its vertices touch." ) {
303  REQUIRE( touched_cover.computeNbCells( 0 )
304  < intersected_cover.computeNbCells( 0 ) );
305  REQUIRE( touched_cover.computeNbCells( 1 )
306  < intersected_cover.computeNbCells( 1 ) );
307  REQUIRE( touched_cover.computeNbCells( 2 )
308  < intersected_cover.computeNbCells( 2 ) );
309  REQUIRE( touched_cover.computeNbCells( 3 )
310  < intersected_cover.computeNbCells( 3 ) );
311  }
312  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
313  REQUIRE( touched_points_cover.computeNbCells( 0 )
314  <= intersected_cover.computeNbCells( 0 ) );
315  REQUIRE( touched_points_cover.computeNbCells( 1 )
316  <= intersected_cover.computeNbCells( 1 ) );
317  REQUIRE( touched_points_cover.computeNbCells( 2 )
318  <= intersected_cover.computeNbCells( 2 ) );
319  REQUIRE( touched_points_cover.computeNbCells( 3 )
320  <= intersected_cover.computeNbCells( 3 ) );
321  }
322  THEN( "The cells touched by its inside points is a subset of the cells its convex hull intersects." ) {
323  REQUIRE( touched_points_cover.subset( intersected_cover, 0 ) );
324  REQUIRE( touched_points_cover.subset( intersected_cover, 1 ) );
325  REQUIRE( touched_points_cover.subset( intersected_cover, 2 ) );
326  REQUIRE( touched_points_cover.subset( intersected_cover, 3 ) );
327  }
328  }
329 }
330 
331 SCENARIO( "CellGeometry< Z2 > rational intersections",
332  "[cell_geometry][2d][rational]" )
333 {
334  typedef KhalimskySpaceND<2,int> KSpace;
335  typedef KSpace::Point Point;
336  typedef KSpace::Space Space;
337  typedef CellGeometry< KSpace > CGeometry;
338  typedef BoundedRationalPolytope< Space > Polytope;
339 
340  GIVEN( "A rational simplex P={ Point(0/4,0/4), Point(17/4,8/4), Point(-5/4,15/4) }" ) {
341  KSpace K;
342  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
343  std::vector< Point > V = { Point(0,0), Point(17,8), Point(-5,15) };
344  Polytope P( 4, V.begin(), V.end() );
345  CGeometry intersected_cover( K, 0, 2, false );
346  intersected_cover.addCellsTouchingPolytope( P );
347  CGeometry touched_points_cover( K, 0, 2, false );
348  touched_points_cover.addCellsTouchingPolytopePoints( P );
349  // trace.info() << "Polytope P=" << P << std::endl;
350  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
351  REQUIRE( intersected_cover.computeEuler() == 1 );
352  }
353  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
354  REQUIRE( touched_points_cover.computeNbCells( 0 )
355  <= intersected_cover.computeNbCells( 0 ) );
356  REQUIRE( touched_points_cover.computeNbCells( 1 )
357  <= intersected_cover.computeNbCells( 1 ) );
358  REQUIRE( touched_points_cover.computeNbCells( 2 )
359  <= intersected_cover.computeNbCells( 2 ) );
360  }
361  }
362  GIVEN( "A thin rational simplex P={ Point(6/4,6/4), Point(17/4,8/4), Point(-5/4,15/4) }" ) {
363  KSpace K;
364  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
365  std::vector< Point > V = { Point(6,6), Point(17,8), Point(-5,15) };
366  Polytope P( 4, V.begin(), V.end() );
367  CGeometry intersected_cover( K, 0, 2, false );
368  intersected_cover.addCellsTouchingPolytope( P );
369  CGeometry touched_points_cover( K, 0, 2, false );
370  touched_points_cover.addCellsTouchingPolytopePoints( P );
371  // trace.info() << "Polytope P=" << P << std::endl;
372  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
373  REQUIRE( intersected_cover.computeEuler() == 1 );
374  }
375  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
376  REQUIRE( touched_points_cover.computeNbCells( 0 )
377  <= intersected_cover.computeNbCells( 0 ) );
378  REQUIRE( touched_points_cover.computeNbCells( 1 )
379  <= intersected_cover.computeNbCells( 1 ) );
380  REQUIRE( touched_points_cover.computeNbCells( 2 )
381  <= intersected_cover.computeNbCells( 2 ) );
382  }
383  }
384 } // SCENARIO( "CellGeometry< Z2 > rational intersections","[cell_geometry][2d][rational]" )
385 
386 
387 SCENARIO( "CellGeometry< Z3 > rational intersections",
388  "[cell_geometry][3d]{rational]" )
389 {
390  typedef KhalimskySpaceND<3,int> KSpace;
391  typedef KSpace::Point Point;
392  typedef KSpace::Space Space;
393  typedef CellGeometry< KSpace > CGeometry;
394  typedef BoundedRationalPolytope< Space > Polytope;
395 
396  GIVEN( "A simplex P={ Point(1/2,0/2,-1/2), Point(7/2,3/2,1/2), Point(-2/2,9/2,3/2), Point(6/2,7/2,10/2) }" ) {
397  KSpace K;
398  K.init( Point( -5, -5, -5 ), Point( 10, 10, 10 ), true );
399  CGeometry intersected_cover( K, 0, 3, false );
400  Polytope P = { Point(2,2,2),
401  Point(1,0,-1), Point(7,3,1), Point(-2,9,3), Point(6,7,10) };
402  intersected_cover.addCellsTouchingPolytope( P );
403  CGeometry touched_points_cover( K, 0, 3, false );
404  touched_points_cover.addCellsTouchingPolytopePoints( P );
405  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
406  REQUIRE( intersected_cover.computeEuler() == -1 );
407  }
408  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
409  REQUIRE( touched_points_cover.computeNbCells( 0 )
410  <= intersected_cover.computeNbCells( 0 ) );
411  REQUIRE( touched_points_cover.computeNbCells( 1 )
412  <= intersected_cover.computeNbCells( 1 ) );
413  REQUIRE( touched_points_cover.computeNbCells( 2 )
414  <= intersected_cover.computeNbCells( 2 ) );
415  REQUIRE( touched_points_cover.computeNbCells( 3 )
416  <= intersected_cover.computeNbCells( 3 ) );
417  }
418  THEN( "The cells touched by its inside points is a subset of the cells its convex hull intersects." ) {
419  REQUIRE( touched_points_cover.subset( intersected_cover, 0 ) );
420  REQUIRE( touched_points_cover.subset( intersected_cover, 1 ) );
421  REQUIRE( touched_points_cover.subset( intersected_cover, 2 ) );
422  REQUIRE( touched_points_cover.subset( intersected_cover, 3 ) );
423  }
424  }
425 } // SCENARIO( "CellGeometry< Z3 > rational intersections", "[cell_geometry][3d]{rational]" )
426 
427 
428 
429 // //
DGtal::BoundedLatticePolytope
Aim: Represents an nD lattice polytope, i.e. a convex polyhedron bounded with vertices with integer c...
Definition: BoundedLatticePolytope.h:74
DGtal::BoundedLatticePolytope::getKPoints
void getKPoints(std::vector< Point > &pts, const Point &alpha_shift) const
DGtal::BoundedRationalPolytope
Aim: Represents an nD rational polytope, i.e. a convex polyhedron bounded by vertices with rational c...
Definition: BoundedRationalPolytope.h:75
DGtal::DigitalConvexity::makeSimplex
static LatticePolytope makeSimplex(PointIterator itB, PointIterator itE)
DGtal::HyperRectDomain::end
const ConstIterator & end() const
Definition: HyperRectDomain.h:201
DGtal::HyperRectDomain::begin
const ConstIterator & begin() const
Definition: HyperRectDomain.h:176
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:394
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
SCENARIO
SCENARIO("CellGeometry< Z2 > segment tests", "[cell_geometry][2d][segment]")
Definition: testCellGeometry.cpp:50
Point
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE
CAPTURE(thicknessHV)
K
KSpace K
Definition: testCubicalComplex.cpp:62
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition: testCubicalComplex.cpp:70
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
REQUIRE
REQUIRE(domain.isInside(aPoint))