DGtal  1.2.0
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 "DGtalCatch.h"
40 
41 using namespace std;
42 using namespace DGtal;
43 
44 
46 // Functions for testing class CellGeometry.
48 
49 SCENARIO( "CellGeometry< Z2 > unit tests", "[cell_geometry][2d]" )
50 {
52  typedef KSpace::Point Point;
53  typedef KSpace::Vector Vector;
54  typedef KSpace::Integer Integer;
55  typedef CellGeometry< KSpace > CGeometry;
56 
57  KSpace K;
58  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
59  GIVEN( "Some points (0,0), Point(2,1), Point(1,3)" ) {
60  std::vector< Point > V = { Point(0,0), Point(2,1), Point(1,3) };
61  CGeometry geometry( K, 0, 2, false );
62  geometry.addCellsTouchingPoints( V.begin(), V.end() );
63  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
64  REQUIRE( V.size() == geometry.computeNbCells( 0 ) );
65  REQUIRE( V.size() < geometry.computeNbCells( 1 ) );
66  REQUIRE( V.size() < geometry.computeNbCells( 2 ) );
67  }
68  THEN( "Its cells form an open complex with euler characteristic 3" ) {
69  REQUIRE( geometry.computeEuler() == 3 );
70  }
71  }
72 }
73 
74 SCENARIO( "CellGeometry< Z3 > unit tests", "[cell_geometry][3d]" )
75 {
77  typedef KSpace::Space Space;
78  typedef KSpace::Point Point;
79  typedef KSpace::Vector Vector;
80  typedef KSpace::Integer Integer;
81  typedef CellGeometry< KSpace > CGeometry;
83 
84  KSpace K;
85  const int N = 10;
86  K.init( Point( -1, -1, -1 ), Point( N, N, N ), true );
87  Domain D( Point( 0, 0, 0 ), Point( N-1, N-1, N-1 ) );
88 
89  GIVEN( "Some a block of points (0,0,0)-(N-1,N-1,N-1)" ) {
90  CGeometry geometry( K, 0, 3, false );
91  geometry.addCellsTouchingPoints( D.begin(), D.end() );
92  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
93  REQUIRE( D.size() == geometry.computeNbCells( 0 ) );
94  REQUIRE( D.size() < geometry.computeNbCells( 1 ) );
95  REQUIRE( D.size() < geometry.computeNbCells( 2 ) );
96  REQUIRE( D.size() < geometry.computeNbCells( 3 ) );
97  }
98  THEN( "Its cells form an open complex with euler characteristic -1" ) {
99  REQUIRE( geometry.computeEuler() == -1 );
100  }
101  }
102 
103  GIVEN( "Some a block of points (0,0,0)-(N-1,N-1,N-1)" ) {
104  CGeometry geometry( K, 0, 2, false );
105  geometry.addCellsTouchingPoints( D.begin(), D.end() );
106  THEN( "Its cell geometry contains more 1-cells and 2-cells than points" ) {
107  REQUIRE( D.size() == geometry.computeNbCells( 0 ) );
108  REQUIRE( D.size() < geometry.computeNbCells( 1 ) );
109  REQUIRE( D.size() < geometry.computeNbCells( 2 ) );
110  }
111  }
112 }
113 
114 SCENARIO( "CellGeometry< Z2 > intersections", "[cell_geometry][2d]" )
115 {
117  typedef KSpace::Point Point;
118  typedef KSpace::Vector Vector;
119  typedef KSpace::Integer Integer;
120  typedef KSpace::Space Space;
121  typedef CellGeometry< KSpace > CGeometry;
122  typedef BoundedLatticePolytope< Space > Polytope;
123 
124  GIVEN( "A simplex P={ Point(0,0), Point(4,2), Point(-1,4) }" ) {
125  KSpace K;
126  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
127  std::vector< Point > V = { Point(0,0), Point(4,2), Point(-1,4) };
128  Polytope P( V.begin(), V.end() );
129  CGeometry intersected_cover( K, 0, 2, false );
130  intersected_cover.addCellsTouchingPolytope( P );
131  CGeometry touched_cover( K, 0, 2, false );
132  touched_cover.addCellsTouchingPoints( V.begin(), V.end() );
133  CGeometry touched_points_cover( K, 0, 2, false );
134  touched_points_cover.addCellsTouchingPolytopePoints( P );
135  // trace.info() << "Polytope P=" << P << std::endl;
136  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
137  REQUIRE( intersected_cover.computeEuler() == 1 );
138  }
139  THEN( "Its convex hull intersects more cells than its vertices touch." ) {
140  REQUIRE( touched_cover.computeNbCells( 0 )
141  < intersected_cover.computeNbCells( 0 ) );
142  REQUIRE( touched_cover.computeNbCells( 1 )
143  < intersected_cover.computeNbCells( 1 ) );
144  REQUIRE( touched_cover.computeNbCells( 2 )
145  < intersected_cover.computeNbCells( 2 ) );
146  }
147  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
148  REQUIRE( touched_points_cover.computeNbCells( 0 )
149  <= intersected_cover.computeNbCells( 0 ) );
150  REQUIRE( touched_points_cover.computeNbCells( 1 )
151  <= intersected_cover.computeNbCells( 1 ) );
152  REQUIRE( touched_points_cover.computeNbCells( 2 )
153  <= intersected_cover.computeNbCells( 2 ) );
154  }
155  }
156 }
157 
158 SCENARIO( "CellGeometry< Z3 > intersections", "[cell_geometry][3d]" )
159 {
161  typedef KSpace::Point Point;
162  typedef KSpace::Vector Vector;
163  typedef KSpace::Integer Integer;
164  typedef KSpace::Space Space;
165  typedef CellGeometry< KSpace > CGeometry;
166  typedef BoundedLatticePolytope< Space > Polytope;
167 
168  GIVEN( "A simplex P={ Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) }" ) {
169  KSpace K;
170  K.init( Point( -5, -5, -5 ), Point( 10, 10, 10 ), true );
171  CGeometry intersected_cover( K, 0, 3, false );
172  Polytope P = { Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) };
173  intersected_cover.addCellsTouchingPolytope( P );
174  std::vector< Point > V = { Point(0,0,0), Point(4,2,1), Point(-1,4,1), Point(3,3,5) };
175  CGeometry touched_cover( K, 0, 3, false );
176  touched_cover.addCellsTouchingPoints( V.begin(), V.end() );
177  CGeometry touched_points_cover( K, 0, 3, false );
178  touched_points_cover.addCellsTouchingPolytopePoints( P );
179  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
180  REQUIRE( intersected_cover.computeEuler() == -1 );
181  }
182  THEN( "Its convex hull intersects more cells than its vertices touch." ) {
183  REQUIRE( touched_cover.computeNbCells( 0 )
184  < intersected_cover.computeNbCells( 0 ) );
185  REQUIRE( touched_cover.computeNbCells( 1 )
186  < intersected_cover.computeNbCells( 1 ) );
187  REQUIRE( touched_cover.computeNbCells( 2 )
188  < intersected_cover.computeNbCells( 2 ) );
189  REQUIRE( touched_cover.computeNbCells( 3 )
190  < intersected_cover.computeNbCells( 3 ) );
191  }
192  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
193  REQUIRE( touched_points_cover.computeNbCells( 0 )
194  <= intersected_cover.computeNbCells( 0 ) );
195  REQUIRE( touched_points_cover.computeNbCells( 1 )
196  <= intersected_cover.computeNbCells( 1 ) );
197  REQUIRE( touched_points_cover.computeNbCells( 2 )
198  <= intersected_cover.computeNbCells( 2 ) );
199  REQUIRE( touched_points_cover.computeNbCells( 3 )
200  <= intersected_cover.computeNbCells( 3 ) );
201  }
202  THEN( "The cells touched by its inside points is a subset of the cells its convex hull intersects." ) {
203  REQUIRE( touched_points_cover.subset( intersected_cover, 0 ) );
204  REQUIRE( touched_points_cover.subset( intersected_cover, 1 ) );
205  REQUIRE( touched_points_cover.subset( intersected_cover, 2 ) );
206  REQUIRE( touched_points_cover.subset( intersected_cover, 3 ) );
207  }
208  }
209 }
210 
211 SCENARIO( "CellGeometry< Z2 > rational intersections",
212  "[cell_geometry][2d][rational]" )
213 {
215  typedef KSpace::Point Point;
216  typedef KSpace::Vector Vector;
217  typedef KSpace::Integer Integer;
218  typedef KSpace::Space Space;
219  typedef CellGeometry< KSpace > CGeometry;
220  typedef BoundedRationalPolytope< Space > Polytope;
221 
222  GIVEN( "A rational simplex P={ Point(0/4,0/4), Point(17/4,8/4), Point(-5/4,15/4) }" ) {
223  KSpace K;
224  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
225  std::vector< Point > V = { Point(0,0), Point(17,8), Point(-5,15) };
226  Polytope P( 4, V.begin(), V.end() );
227  CGeometry intersected_cover( K, 0, 2, false );
228  intersected_cover.addCellsTouchingPolytope( P );
229  CGeometry touched_points_cover( K, 0, 2, false );
230  touched_points_cover.addCellsTouchingPolytopePoints( P );
231  // trace.info() << "Polytope P=" << P << std::endl;
232  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
233  REQUIRE( intersected_cover.computeEuler() == 1 );
234  }
235  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
236  REQUIRE( touched_points_cover.computeNbCells( 0 )
237  <= intersected_cover.computeNbCells( 0 ) );
238  REQUIRE( touched_points_cover.computeNbCells( 1 )
239  <= intersected_cover.computeNbCells( 1 ) );
240  REQUIRE( touched_points_cover.computeNbCells( 2 )
241  <= intersected_cover.computeNbCells( 2 ) );
242  }
243  }
244  GIVEN( "A thin rational simplex P={ Point(6/4,6/4), Point(17/4,8/4), Point(-5/4,15/4) }" ) {
245  KSpace K;
246  K.init( Point( -5, -5 ), Point( 10, 10 ), true );
247  std::vector< Point > V = { Point(6,6), Point(17,8), Point(-5,15) };
248  Polytope P( 4, V.begin(), V.end() );
249  CGeometry intersected_cover( K, 0, 2, false );
250  intersected_cover.addCellsTouchingPolytope( P );
251  CGeometry touched_points_cover( K, 0, 2, false );
252  touched_points_cover.addCellsTouchingPolytopePoints( P );
253  // trace.info() << "Polytope P=" << P << std::endl;
254  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
255  REQUIRE( intersected_cover.computeEuler() == 1 );
256  }
257  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
258  REQUIRE( touched_points_cover.computeNbCells( 0 )
259  <= intersected_cover.computeNbCells( 0 ) );
260  REQUIRE( touched_points_cover.computeNbCells( 1 )
261  <= intersected_cover.computeNbCells( 1 ) );
262  REQUIRE( touched_points_cover.computeNbCells( 2 )
263  <= intersected_cover.computeNbCells( 2 ) );
264  }
265  }
266 } // SCENARIO( "CellGeometry< Z2 > rational intersections","[cell_geometry][2d][rational]" )
267 
268 
269 SCENARIO( "CellGeometry< Z3 > rational intersections",
270  "[cell_geometry][3d]{rational]" )
271 {
273  typedef KSpace::Point Point;
274  typedef KSpace::Vector Vector;
275  typedef KSpace::Integer Integer;
276  typedef KSpace::Space Space;
277  typedef CellGeometry< KSpace > CGeometry;
278  typedef BoundedRationalPolytope< Space > Polytope;
279 
280  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) }" ) {
281  KSpace K;
282  K.init( Point( -5, -5, -5 ), Point( 10, 10, 10 ), true );
283  CGeometry intersected_cover( K, 0, 3, false );
284  Polytope P = { Point(2,2,2),
285  Point(1,0,-1), Point(7,3,1), Point(-2,9,3), Point(6,7,10) };
286  intersected_cover.addCellsTouchingPolytope( P );
287  CGeometry touched_points_cover( K, 0, 3, false );
288  touched_points_cover.addCellsTouchingPolytopePoints( P );
289  THEN( "The cells intersected by its convex hull form an open and simply connected complex." ) {
290  REQUIRE( intersected_cover.computeEuler() == -1 );
291  }
292  THEN( "Its convex hull intersects at least as many cells as its inside points touch." ) {
293  REQUIRE( touched_points_cover.computeNbCells( 0 )
294  <= intersected_cover.computeNbCells( 0 ) );
295  REQUIRE( touched_points_cover.computeNbCells( 1 )
296  <= intersected_cover.computeNbCells( 1 ) );
297  REQUIRE( touched_points_cover.computeNbCells( 2 )
298  <= intersected_cover.computeNbCells( 2 ) );
299  REQUIRE( touched_points_cover.computeNbCells( 3 )
300  <= intersected_cover.computeNbCells( 3 ) );
301  }
302  THEN( "The cells touched by its inside points is a subset of the cells its convex hull intersects." ) {
303  REQUIRE( touched_points_cover.subset( intersected_cover, 0 ) );
304  REQUIRE( touched_points_cover.subset( intersected_cover, 1 ) );
305  REQUIRE( touched_points_cover.subset( intersected_cover, 2 ) );
306  REQUIRE( touched_points_cover.subset( intersected_cover, 3 ) );
307  }
308  }
309 } // SCENARIO( "CellGeometry< Z3 > rational intersections", "[cell_geometry][3d]{rational]" )
310 
311 
312 
313 // //
Aim: Represents an nD lattice polytope, i.e. a convex polyhedron bounded with vertices with integer c...
Aim: Represents an nD rational polytope, i.e. a convex polyhedron bounded by vertices with rational c...
Aim: Computes and stores sets of cells and provides methods to compute intersections of lattice and r...
Definition: CellGeometry.h:75
const ConstIterator & end() const
const ConstIterator & begin() const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
TInteger Integer
Arithmetic ring induced by (+,-,*) and Integer numbers.
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
DGtal is the top-level namespace which contains all DGtal functions and types.
SCENARIO("CellGeometry< Z2 > unit tests", "[cell_geometry][2d]")
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
KSpace K
GIVEN("A cubical complex with random 3-cells")
HyperRectDomain< Space > Domain
REQUIRE(domain.isInside(aPoint))