DGtal 1.4.0
Loading...
Searching...
No Matches
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
42using namespace std;
43using namespace DGtal;
44
45
47// Functions for testing class CellGeometry.
49
50SCENARIO( "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
117SCENARIO( "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
177SCENARIO( "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
200SCENARIO( "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
238SCENARIO( "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
280SCENARIO( "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
331SCENARIO( "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
387SCENARIO( "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
Aim: Parallelepidec region of a digital space, model of a 'CDomain'.
Definition HyperRectDomain.h:100
DGtal::HyperRectDomain::begin
const ConstIterator & begin() const
Definition HyperRectDomain.h:176
DGtal::HyperRectDomain::end
const ConstIterator & end() const
Definition HyperRectDomain.h:201
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.
Definition ClosedIntegerHalfPlane.h:49
std
STL namespace.
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))
SCENARIO
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")
Definition testUnorderedSetByBlock.cpp:58