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