DGtal 1.3.0
Loading...
Searching...
No Matches
testVoxelComplex.cpp
Go to the documentation of this file.
1
31#include "DGtal/base/SetFunctions.h"
32#include "DGtal/helpers/StdDefs.h"
33#include "DGtal/topology/CubicalComplexFunctions.h"
34#include "DGtal/topology/CubicalComplex.h"
35#include "DGtal/topology/KhalimskyCellHashFunctions.h"
36#include "DGtal/topology/VoxelComplex.h"
37#include "DGtal/topology/VoxelComplexFunctions.h"
38#include "DGtalCatch.h"
39#include <iostream>
40#include <unordered_map>
41
42#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
43#include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
44#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
45#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
46#include "DGtal/kernel/BasicPointPredicates.h"
47#include "DGtal/topology/NeighborhoodConfigurations.h"
48#include "DGtal/topology/tables/NeighborhoodTables.h"
49// #include <DGtal/io/viewers/Viewer3D.h>
51
52using namespace std;
53using namespace DGtal;
54
56// Fixture for object from a diamond set
58struct Fixture_complex_diamond {
60 // type aliases
65
66 using FixtureDigitalTopology = DGtal::Z3i::DT26_6;
67 using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
69 std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
70 using FixtureComplex = DGtal::VoxelComplex<KSpace>;
71
72
74 // Constructor
76 Fixture_complex_diamond() :
77 complex_fixture(ks_fixture) {
78 create_complex_from_set(create_set());
79 }
80
82 // fixture data
83 FixtureComplex complex_fixture;
84 KSpace ks_fixture; // needed because ConstAlias in CC constructor.
86
88 // Function members
90 FixtureDigitalSet create_set() {
91 using namespace DGtal;
92
93 // trace.beginBlock ( "Create Fixture_object_diamond" );
94 Point p1(-10, -10, -10);
95 Point p2(10, 10, 10);
96 Domain domain(p1, p2);
97 Point c(0, 0, 0);
98
99 // diamond of radius 4
100 FixtureDigitalSet set_fixture(domain);
101 for (auto it = domain.begin(); it != domain.end(); ++it) {
102 if ((*it - c).norm1() <= 3)
103 set_fixture.insertNew(*it);
104 }
105 set_fixture.erase(c);
106
107 return set_fixture;
108 }
109
110 FixtureComplex &create_complex_from_set(const FixtureDigitalSet &input_set) {
111
112 ks_fixture.init(input_set.domain().lowerBound(),
113 input_set.domain().upperBound(), true);
114 complex_fixture = FixtureComplex(ks_fixture);
115 complex_fixture.construct(input_set);
116 return complex_fixture;
117 }
118};
119
120TEST_CASE_METHOD(Fixture_complex_diamond, "insertVoxel", "[insert][close]") {
121 auto &vc = complex_fixture;
122 auto &ks = vc.space();
123 auto s3 = vc.nbCells(3);
124 auto s2 = vc.nbCells(2);
125 auto s1 = vc.nbCells(1);
126 auto s0 = vc.nbCells(0);
127 SECTION("insertVoxel in border") {
128 Point p{0, 0, 4};
129 auto vc_new = vc;
130 vc_new.insertVoxelPoint(p);
131 auto sa3 = vc_new.nbCells(3);
132 auto sa2 = vc_new.nbCells(2);
133 auto sa1 = vc_new.nbCells(1);
134 auto sa0 = vc_new.nbCells(0);
135 auto d3 = sa3 - s3;
136 CHECK(d3 == 1);
137 auto d2 = sa2 - s2;
138 CHECK(d2 == 5);
139 auto d1 = sa1 - s1;
140 CHECK(d1 == 8);
141 auto d0 = sa0 - s0;
142 CHECK(d0 == 4);
143 }
144 SECTION("insertVoxel isolated") {
145 Point p{5, 0, 0};
146 auto vc_new = vc;
147 vc_new.insertVoxelCell(ks.uSpel(p));
148 auto sa3 = vc_new.nbCells(3);
149 auto sa2 = vc_new.nbCells(2);
150 auto sa1 = vc_new.nbCells(1);
151 auto sa0 = vc_new.nbCells(0);
152 auto d3 = sa3 - s3;
153 CHECK(d3 == 1);
154 auto d2 = sa2 - s2;
155 CHECK(d2 == 6);
156 auto d1 = sa1 - s1;
157 CHECK(d1 == 12);
158 auto d0 = sa0 - s0;
159 CHECK(d0 == 8);
160 }
161 SECTION("insertVoxel interior") {
162 Point p{0, 2, 0};
163 auto vc_new = vc;
164 vc_new.insertVoxelCell(ks.uSpel(p));
165 auto sa3 = vc_new.nbCells(3);
166 auto sa2 = vc_new.nbCells(2);
167 auto sa1 = vc_new.nbCells(1);
168 auto sa0 = vc_new.nbCells(0);
169 auto d3 = sa3 - s3;
170 CHECK(d3 == 0);
171 auto d2 = sa2 - s2;
172 CHECK(d2 == 0);
173 auto d1 = sa1 - s1;
174 CHECK(d1 == 0);
175 auto d0 = sa0 - s0;
176 CHECK(d0 == 0);
177 }
178}
179
180TEST_CASE_METHOD(Fixture_complex_diamond, "Faces of voxel",
181 "[neighborhood][faces]") {
182 auto &vc = complex_fixture;
183 auto &ks = vc.space();
184 bool closed = true;
185 bool no_hint = false;
186 auto not_found = vc.end(3);
187
188 SECTION("Voxel is in the interior of the complex.") {
189 Point p{0, 0, 2};
190 auto cellMapIt = vc.findCell(3, ks.uSpel(p));
191 CHECK(cellMapIt != not_found);
192 auto cell = cellMapIt->first;
193 auto kfaces = ks.uFaces(cell);
194 CAPTURE(p);
195 CAPTURE(cell);
196 CHECK(kfaces.size() == 26);
197 auto vcBoundary_closed = vc.cellBoundary(cell, closed);
198 CHECK(vcBoundary_closed.size() == 26);
199 auto vcBoundary_no_hint = vc.cellBoundary(cell, no_hint);
200 CHECK(vcBoundary_no_hint.size() == 26);
201 }
202 SECTION("Voxel is in the exterior") {
203 Point p{0, 0, 5};
204 auto cellMapIt = vc.findCell(3, ks.uSpel(p));
205 CHECK(cellMapIt == not_found);
206 auto cell = ks.uSpel(p);
207 auto kfaces = ks.uFaces(cell);
208 CAPTURE(p);
209 CAPTURE(cell);
210 CHECK(kfaces.size() == 26);
211 // We know that is not closed, but checking the differences:
212 auto vcBoundary_closed = vc.cellBoundary(cell, closed);
213 CHECK(vcBoundary_closed.size() == 26);
214 // Right method to call:
215 auto vcBoundary_no_hint = vc.cellBoundary(cell, no_hint);
216 CHECK(vcBoundary_no_hint.size() == 0);
217 }
218 SECTION("Voxel is out, but surrounded by in voxels") {
219 Point p{0, 0, 0};
220 auto cellMapIt = vc.findCell(3, ks.uSpel(p));
221 CHECK(cellMapIt == not_found);
222 auto cell = ks.uSpel(p);
223 auto kfaces = ks.uFaces(cell);
224 CAPTURE(p);
225 CAPTURE(cell);
226 CHECK(kfaces.size() == 26);
227 auto vcBoundary_closed = vc.cellBoundary(cell, closed);
228 CHECK(vcBoundary_closed.size() == 26);
229 auto vcBoundary_no_hint = vc.cellBoundary(cell, no_hint);
230 CHECK(vcBoundary_no_hint.size() == 26);
231 }
232 SECTION("Voxel is in the border") {
233 Point p{0, 0, 3};
234 auto cellMapIt = vc.findCell(3, ks.uSpel(p));
235 CHECK(cellMapIt != not_found);
236 auto cell = cellMapIt->first;
237 auto kfaces = ks.uFaces(cell);
238 CAPTURE(p);
239 CAPTURE(cell);
240 CHECK(kfaces.size() == 26);
241 auto vcBoundary_closed = vc.cellBoundary(cell, closed);
242 CHECK(vcBoundary_closed.size() == 26);
243 auto vcBoundary_no_hint = vc.cellBoundary(cell, no_hint);
244 CHECK(vcBoundary_no_hint.size() == 26);
245 }
246}
247
248TEST_CASE_METHOD(Fixture_complex_diamond, "Neighbors from Object and KSpace",
249 "[neighborhood]") {
250 auto &vc = complex_fixture;
251 SECTION(" get neighbors from Kspace") {
252 {
253 size_t dim_voxel = 3;
254 auto it = vc.begin(dim_voxel);
255 for (std::size_t n = 0; n < 10; ++n)
256 ++it;
257 auto cell = it->first;
258 auto point_from_kspace_1 = cell.preCell().coordinates;
259 size_t expected_num_adjacent_voxels = 12;
260
261 SECTION("properNeighborhood from KSpace"
262 "does not output all adjacent voxels.") {
263
264 size_t expected_kspace_neighborhood = 7;
265 auto neigh_k = vc.space().uNeighborhood(cell);
266 CHECK(neigh_k.size() == expected_kspace_neighborhood);
267
268 auto propN_k = vc.space().uProperNeighborhood(cell);
269 CHECK(propN_k.size() == expected_kspace_neighborhood - 1);
270 }
271
272 SECTION("Getting associated pointels and voxels from input_cell") {
273 std::set<FixtureComplex::Cell> pointel_set;
274 vc.pointelsFromCell(pointel_set, cell);
275 std::set<FixtureComplex::Cell> voxel_set;
276 for (auto &&p : pointel_set)
277 vc.spelsFromCell(voxel_set, p);
278 SECTION("Gets desired full adjancency for voxels") {
279 CHECK(voxel_set.size() == expected_num_adjacent_voxels);
280 }
281 auto clique = vc.Kneighborhood(cell);
282 CHECK(clique.nbCells(3) == expected_num_adjacent_voxels);
283 }
284 }
285 }
286}
287
288TEST_CASE_METHOD(Fixture_complex_diamond, "Test Simplicity", "[simplicity]") {
289 auto &vc = complex_fixture;
290 size_t dim_voxel = 3;
291 auto cit = vc.begin(dim_voxel);
292 for (std::size_t n = 0; n < 10; ++n)
293 ++cit;
294 auto cell = cit->first;
295
296 SECTION("querying voxel simplicity") {
297 size_t nsimples{0};
298 std::vector<FixtureComplex::Cell> non_simple_cells;
299 for (auto it = vc.begin(dim_voxel); it != vc.end(dim_voxel); ++it) {
300 if (vc.isSimple(it->first))
301 ++nsimples;
302 else
303 non_simple_cells.emplace_back(it->first);
304 }
305
306 // Border points are simple in diamond.
307 // auto border_size = vc.object().border().size();
308 size_t border_size = 44;
309 REQUIRE(nsimples == border_size);
310 }
311}
312
313TEST_CASE_METHOD(Fixture_complex_diamond, "Test table wrappers",
314 "[table][simple]") {
315 auto &vc = complex_fixture;
316 trace.beginBlock("loadTable");
317 vc.setSimplicityTable(functions::loadTable(simplicity::tableSimple26_6));
318 trace.endBlock();
319 auto dim_voxel = 3;
320 auto cit = vc.begin(dim_voxel);
321 for (std::size_t n = 0; n < 10; ++n)
322 ++cit;
323 auto cell = cit->first;
324
325 SECTION("querying voxel simplicity") {
326 size_t nsimples{0};
327 std::vector<FixtureComplex::Cell> non_simple_cells;
328 for (auto it = vc.begin(dim_voxel); it != vc.end(dim_voxel); ++it) {
329 if (vc.isSimple(it->first))
330 ++nsimples;
331 else
332 non_simple_cells.emplace_back(it->first);
333 }
334
335 // Border points are simple in diamond.
336 size_t border_size = 44;
337 REQUIRE(nsimples == border_size);
338 }
339}
340
341TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_2", "[clique]") {
342 auto &vc = complex_fixture;
343 auto itc = vc.begin(3);
344 for (std::size_t n = 0; n < 10; ++n)
345 ++itc;
346 auto cell = itc->first;
347 Point p_cell = vc.space().uCoords(cell);
348 auto neigh6 = vc.space().uProperNeighborhood(cell);
349 REQUIRE(neigh6.size() == 6);
350 KSpace::Cells cells;
351 for(auto & n : neigh6)
352 if(vc.belongs(n)) cells.push_back(n);
353 REQUIRE(cells.size() == 2);
354
355 std::vector<std::pair<bool, FixtureComplex::Clique>> cliques_p;
356 bool verbose = true;
357 for (auto && cell_n : cells) {
358 auto cell_point = vc.space().uCoords(cell_n);
359 cliques_p.emplace_back(vc.K_2(p_cell, cell_point, verbose));
360 auto &is_critical = cliques_p.back().first;
361 auto &k2_clique = cliques_p.back().second;
362 CHECK(is_critical == true);
363 CHECK(k2_clique.nbCells(3) == 2);
364 }
365
366 SECTION(" Check same results for different K_2 interfaces") {
367 for (auto && cell_n : cells) {
368 auto co_face = vc.surfelBetweenAdjacentSpels(cell_n, cell);
369
370 auto cell_point = vc.space().uCoords(cell_n);
371 using namespace DGtal::functions;
372 CHECK(isEqual(vc.K_2(p_cell, cell_point, false).second,
373 vc.K_2(cell, cell_n, false).second) == true);
374 CHECK(isEqual(vc.K_2(p_cell, cell_point, false).second,
375 vc.K_2(co_face, false).second) == true);
376 } // endfor
377 } // then_interfaces
378} // test
379
380TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_1", "[clique]") {
381 auto &vc = complex_fixture;
382 auto itc = vc.begin(3);
383 for (std::size_t n = 0; n < 10; ++n)
384 ++itc;
385 auto cell = itc->first;
386 Point p_cell = vc.space().uCoords(cell);
387
388 SECTION("K_1 mask from a linel") {
389 auto linel =
390 vc.space().uCell(cell.preCell().coordinates + Point{0, 1, 1});
391
392 using namespace DGtal::functions;
393 auto k1p = vc.K_1(linel, true);
394 auto &is_critical = k1p.first;
395 CHECK(is_critical == true);
396 auto &k1_clique = k1p.second;
397 // numer of voxels in clique
398 CHECK(k1_clique.nbCells(3) == 3);
399 }
400} // test
401
402TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_0", "[clique]") {
403 auto &vc = complex_fixture;
404 auto itc = vc.begin(3);
405 for (std::size_t n = 0; n < 10; ++n)
406 ++itc;
407 auto cell = itc->first;
408 Point p_cell = vc.space().uCoords(cell);
409
410 SECTION("K_0 mask from a pointel") {
411 auto pointel = vc.space().uPointel(p_cell);
412 // auto pointel = vc.space().uPointel(Point{0,0,0});
413
414 using namespace DGtal::functions;
415 auto k0_mask = vc.K_0(pointel, true);
416
417 auto &is_critical = k0_mask.first;
418 CHECK(is_critical == false);
419 auto &k0_clique = k0_mask.second;
420 // numer of voxels in clique
421 CHECK(k0_clique.nbCells(3) == 1);
422 }
423} // test
424
425TEST_CASE_METHOD(Fixture_complex_diamond, "Get All Critical Cliques of diamond",
426 "[critical][clique]") {
427 auto &vc = complex_fixture;
428 using namespace DGtal::functions;
429 SECTION(" Get all criticalCliques() ") {
430 auto criticals = vc.criticalCliques();
431 CHECK(criticals.size() == 4);
432
433 CHECK(vc.nbCells(3) == 62);
434 CHECK(criticals[3].size() == 18);
435 CHECK(vc.nbCells(2) == 264);
436 CHECK(criticals[2].size() == 108);
437 CHECK(vc.nbCells(1) == 360);
438 CHECK(criticals[1].size() == 168);
439 CHECK(vc.nbCells(0) == 160);
440 CHECK(criticals[0].size() == 32);
441 }
442}
443
444
446// Fixture for complex fig 4 of Asymmetric parallel 3D thinning scheme
448struct Fixture_complex_fig4 {
450 // type aliases
452 using Point = DGtal::Z3i::Point;
455
456 using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
458 std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
459 using FixtureMap = std::unordered_map<KSpace::Cell, CubicalCellData>;
460 using FixtureComplex = DGtal::VoxelComplex<KSpace, FixtureMap>;
461
463 // fixture data
464 FixtureComplex complex_fixture;
465 KSpace ks_fixture; // needed because ConstAlias in CC constructor.
467
469 // Constructor
471 Fixture_complex_fig4() : complex_fixture(ks_fixture) {
472 create_complex_from_set(create_set());
473 }
474
476 // Function members
478 FixtureDigitalSet create_set() {
479 using namespace DGtal;
480
481 Point p1(-10, -10, -10);
482 Point p2(10, 10, 10);
483 Domain domain(p1, p2);
484
485 // 12 voxels of fig4, centered in critical voxel
486 FixtureDigitalSet fig4_set(domain);
487 Point a1(0, 0, 0);
488 fig4_set.insertNew(a1);
489 Point a2(0, -1, 0);
490 fig4_set.insertNew(a2);
491 Point a3(1, -1, 0);
492 fig4_set.insertNew(a3);
493 Point c1(0, 1, -1);
494 fig4_set.insertNew(c1);
495 Point c2(1, 1, -1);
496 fig4_set.insertNew(c2);
497 Point c3(1, 2, -1);
498 fig4_set.insertNew(c3);
499 Point c4(0, 2, -1);
500 fig4_set.insertNew(c4);
501 Point b1(0, 2, 0);
502 fig4_set.insertNew(b1);
503 Point b2(-1, 3, 0);
504 fig4_set.insertNew(b2);
505 Point l1(1, 0, -2);
506 fig4_set.insertNew(l1);
507 Point l2(2, 0, -2);
508 fig4_set.insertNew(l2);
509 Point r1(1, 3, -2);
510 fig4_set.insertNew(r1);
511
512 return fig4_set;
513 }
514
515 FixtureComplex &create_complex_from_set(const FixtureDigitalSet &input_set) {
516
517 ks_fixture.init(input_set.domain().lowerBound(),
518 input_set.domain().upperBound(), true);
519 complex_fixture = FixtureComplex(ks_fixture);
520 complex_fixture.construct(input_set);
521 return complex_fixture;
522 }
523};
524TEST_CASE_METHOD(Fixture_complex_fig4, "Get All Critical Cliques of fig4",
525 "[critical][clique]") {
526 auto &vc = complex_fixture;
527 using namespace DGtal::functions;
528 SECTION(" Get all criticalCliques() ") {
529 auto criticals = vc.criticalCliques(true);
530 CHECK(criticals.size() == 4);
531
532 CHECK(vc.nbCells(3) == 12);
533 CHECK(criticals[3].size() == 1);
534 CHECK(vc.nbCells(2) == 64);
535 CHECK(criticals[2].size() == 3);
536 CHECK(vc.nbCells(1) == 109);
537 CHECK(criticals[1].size() == 2);
538 CHECK(vc.nbCells(0) == 58);
539 CHECK(criticals[0].size() == 6);
540 }
541}
542
543/* zeroSurface and oneSurface */
544TEST_CASE_METHOD(Fixture_complex_fig4, "zeroSurface and oneSurface",
545 "[isSurface][function]") {
546 auto &vc = complex_fixture;
547 using namespace DGtal::functions;
548 vc.clear();
549 Point c(0, 0, 0);
550 {
551 vc.insertCell(3, vc.space().uSpel(c));
552 }
553 Point r(0, 1, 0);
554 {
555 vc.insertCell(3, vc.space().uSpel(r));
556 }
557 Point l(0, -1, 0);
558 {
559 vc.insertCell(3, vc.space().uSpel(l));
560 }
561 // Init an object
565 std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
569 auto domain = DGtal::Z3i::Domain(ks_fixture.lowerBound(), ks_fixture.upperBound());
570 DigitalTopology topo(
572
573 SECTION("checking zero surfaces of original set") {
574 std::set<Point> zeroSurfacesCells;
575 DigitalSet voxel_set(domain);
576 vc.dumpVoxels(voxel_set);
577 Object object(topo, voxel_set);
578 for (const auto &p : voxel_set) {
579 auto small_obj = object.properNeighborhood(p);
580 if (isZeroSurface(small_obj))
581 zeroSurfacesCells.insert(p);
582 }
583 CHECK(zeroSurfacesCells.size() == 1);
584 }
585
586 SECTION("checking one surfaces of original set") {
587 Point u(-1, 0, 0);
588 {
589 vc.insertCell(3, vc.space().uSpel(u));
590 }
591 Point d(1, 0, 0);
592 {
593 vc.insertCell(3, vc.space().uSpel(d));
594 }
595
596 std::set<Point> oneSurfacesCells;
597 DigitalSet voxel_set(domain);
598 vc.dumpVoxels(voxel_set);
599 Object object(topo, voxel_set);
600
601 for (const auto &p : voxel_set) {
602 auto small_obj = object.properNeighborhood(p);
603 if (isOneSurface(small_obj))
604 oneSurfacesCells.insert(p);
605 }
606 CHECK(oneSurfacesCells.size() == 1);
607 }
608}
609
611// Fixture for complex fig 4 of Asymmetric parallel 3D thinning scheme
613struct Fixture_isthmus {
615 // type aliases
617 using Point = DGtal::Z3i::Point;
620
621 using FixtureDigitalTopology = DGtal::Z3i::DT26_6;
622 using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
624 std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
625 using FixtureMap = std::unordered_map<KSpace::Cell, CubicalCellData>;
626 using FixtureComplex = DGtal::VoxelComplex<KSpace, FixtureMap>;
627
629 // fixture data
630 FixtureComplex complex_fixture;
631 FixtureDigitalSet set_fixture;
632 KSpace ks_fixture; // needed because ConstAlias in CC constructor.
634
636 // Constructor
638 Fixture_isthmus() : complex_fixture(ks_fixture), set_fixture(create_set()) {
639 create_complex_from_set(set_fixture);
640 }
641
643 // Function members
645 FixtureDigitalSet create_set() {
646 using namespace DGtal;
647
648 Point p1(-10, -10, -10);
649 Point p2(10, 10, 10);
650 Domain domain(p1, p2);
651
652 FixtureDigitalSet fig6_set(domain);
653
654 Point c00(0, 0, 0);
655 fig6_set.insertNew(c00);
656 Point c01(-1, 0, 0);
657 fig6_set.insertNew(c01);
658 Point c02(-2, 0, 0);
659 fig6_set.insertNew(c02);
660 Point c03(-2, 0, 1);
661 fig6_set.insertNew(c03);
662
663 Point c10(0, 1, 0);
664 fig6_set.insertNew(c10);
665 Point y(-1, 1, 0);
666 fig6_set.insertNew(y);
667 Point c12(-2, 1, 0);
668 fig6_set.insertNew(c12);
669 Point c13(-2, 1, 1);
670 fig6_set.insertNew(c13);
671
672 Point c20(0, 2, 0);
673 fig6_set.insertNew(c20);
674 Point c21(-1, 2, 0);
675 fig6_set.insertNew(c21);
676 Point c22(-1, 2, 1);
677 fig6_set.insertNew(c22);
678 Point c33(-1, 2, -1);
679 fig6_set.insertNew(c33);
680
681 Point c30(-1, 3, 0);
682 fig6_set.insertNew(c30);
683 Point c31(-1, 3, 1);
684 fig6_set.insertNew(c31);
685
686 Point x(-1, 4, 0);
687 fig6_set.insertNew(x);
688
689 Point c50(-1, 5, 0);
690 fig6_set.insertNew(c50);
691 Point c51(-1, 5, -1);
692 fig6_set.insertNew(c51);
693 Point c52(-2, 5, 0);
694 fig6_set.insertNew(c52);
695
696 Point c60(-1, 6, 0);
697 fig6_set.insertNew(c60);
698 Point c61(-2, 6, 0);
699 fig6_set.insertNew(c61);
700 Point c62(-2, 6, -1);
701 fig6_set.insertNew(c62);
702
703 return fig6_set;
704 }
705
706 FixtureComplex &create_complex_from_set(const FixtureDigitalSet &input_set) {
707
708 ks_fixture.init(input_set.domain().lowerBound(),
709 input_set.domain().upperBound(), true);
710 complex_fixture = FixtureComplex(ks_fixture);
711 complex_fixture.construct(input_set);
712 return complex_fixture;
713 }
714};
715
716TEST_CASE_METHOD(Fixture_isthmus, "Thin disconnected complex",
717 "[isthmus][thin][function]") {
718 using namespace DGtal::functions;
719 auto &vc = complex_fixture;
720 auto &ks = vc.space();
721 boost::ignore_unused_variable_warning(ks);
722 Point x(-1, 4, 0);
723 set_fixture.erase(x);
724 // Delete one point and reconstruct complex.
725 /* obj_fixture.pointSet().erase(x); */
726 vc.clear();
727 vc.construct(set_fixture);
728 SECTION("with skelUltimate") {
729 auto vc_new = asymetricThinningScheme<FixtureComplex>(
730 vc, selectFirst<FixtureComplex>, skelUltimate<FixtureComplex>);
731 CHECK(vc_new.nbCells(3) == 2);
732 }
733}
734
735
736TEST_CASE_METHOD(Fixture_isthmus, "Check isthmus", "[isthmus][function]") {
737 auto &vc = complex_fixture;
738 using namespace DGtal::functions;
739 auto &ks = vc.space();
740 Point x(-1, 4, 0); // oneIsthmus
741 Point y(-1, 1, 0); // twoIsthmus
742
743 SECTION("checking one isthmus") {
744 std::set<Point> isthmus;
745 for (const auto &p : set_fixture) {
746 if (oneIsthmus(vc, ks.uSpel(p)))
747 isthmus.insert(p);
748 }
749 CHECK(isthmus.size() == 1);
750 auto xit = isthmus.find(x);
751 CHECK(*xit == x);
752 }
753
754 SECTION("checking 2 isthmus") {
755 std::set<Point> isthmus;
756 for (const auto &p : set_fixture) {
757 if (twoIsthmus(vc, ks.uSpel(p)))
758 isthmus.insert(p);
759 }
760 CHECK(isthmus.size() == 1);
761 auto yit = isthmus.find(y);
762 CHECK(*yit == y);
763 }
764
765 SECTION("checking 1 and 2 isthmus") {
766 std::set<Point> isthmus;
767 for (const auto &p : set_fixture) {
768 if (skelIsthmus(vc, ks.uSpel(p)))
769 isthmus.insert(p);
770 }
771 CHECK(isthmus.size() == 2);
772 auto xit = isthmus.find(x);
773 auto yit = isthmus.find(y);
774 CHECK(xit != isthmus.end());
775 CHECK(yit != isthmus.end());
776 }
777}
778TEST_CASE_METHOD(Fixture_isthmus, "Thin complex", "[isthmus][thin][function]") {
779 using namespace DGtal::functions;
780 auto &vc = complex_fixture;
781 auto &ks = vc.space();
782 boost::ignore_unused_variable_warning(ks);
783 SECTION("with skelUltimate") {
784 auto vc_new = asymetricThinningScheme<FixtureComplex>(
785 vc, selectFirst<FixtureComplex>, skelUltimate<FixtureComplex>);
786 CHECK(vc_new.nbCells(3) == 1);
787 }
788 SECTION("with skelEnd") {
789 auto vc_new = asymetricThinningScheme<FixtureComplex>(
790 vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>);
791 CHECK(vc_new.nbCells(3) == 5);
792 }
793 SECTION("with oneIsthmus") {
794 auto vc_new = asymetricThinningScheme<FixtureComplex>(
795 vc, selectRandom<FixtureComplex>, oneIsthmus<FixtureComplex>);
796 CHECK(vc_new.nbCells(3) == 3);
797 }
798 SECTION("with twoIsthmus") {
799 auto vc_new = asymetricThinningScheme<FixtureComplex>(
800 vc, selectRandom<FixtureComplex>, twoIsthmus<FixtureComplex>);
801 CHECK(vc_new.nbCells(3) == 1);
802 }
803 SECTION("with skelIsthmus") {
804 auto vc_new = asymetricThinningScheme<FixtureComplex>(
805 vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>);
806 CHECK(vc_new.nbCells(3) == 3);
807 }
808}
809//
810TEST_CASE_METHOD(Fixture_isthmus, "Persistence thin",
811 "[persistence][isthmus][thin][function]") {
812 using namespace DGtal::functions;
813 auto &vc = complex_fixture;
814 auto &ks = vc.space();
815 boost::ignore_unused_variable_warning(ks);
816 SECTION("with skelUltimate") {
817 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
818 vc, selectRandom<FixtureComplex>, skelUltimate<FixtureComplex>, 0);
819 CHECK(vc_new.nbCells(3) == 1);
820 }
821 SECTION("with skelEnd") {
822 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
823 vc, selectRandom<FixtureComplex>, skelEnd<FixtureComplex>, 0);
824 CHECK(vc_new.nbCells(3) == 5);
825 }
826 SECTION("with oneIsthmus") {
827 // Not using LUT skel function (slow):
828 //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
829 // vc, selectRandom<FixtureComplex>, oneIsthmus<FixtureComplex>, 0);
830 //
831 // with LUT:
832 auto table = *functions::loadTable(isthmusicity::tableOneIsthmus);
833 auto pointToMaskMap =
834 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
835 auto oneIsthmusTable =
836 [&table, &pointToMaskMap](const FixtureComplex &fc,
837 const FixtureComplex::Cell &c) {
838 return skelWithTable(table, pointToMaskMap, fc, c);
839 };
840 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
841 vc, selectRandom<FixtureComplex>, oneIsthmusTable, 0);
842 CHECK(vc_new.nbCells(3) == 3);
843 }
844 SECTION("with twoIsthmus") {
845 // Not using LUT skel function (slow):
846 //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
847 // vc, selectRandom<FixtureComplex>, twoIsthmus<FixtureComplex>, 0);
848 //
849 // with LUT:
850 auto table = *functions::loadTable(isthmusicity::tableTwoIsthmus);
851 auto pointToMaskMap =
852 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
853 auto twoIsthmusTable =
854 [&table, &pointToMaskMap](const FixtureComplex &fc,
855 const FixtureComplex::Cell &c) {
856 return skelWithTable(table, pointToMaskMap, fc, c);
857 };
858 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
859 vc, selectRandom<FixtureComplex>, twoIsthmusTable, 0);
860 CHECK(vc_new.nbCells(3) == 1);
861 }
862 SECTION("with skelIsthmus") {
863 // Not using LUT skel function (slow):
864 //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
865 // vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>, 0);
866 //
867 // with LUT:
868 auto table = *functions::loadTable(isthmusicity::tableIsthmus);
869 auto pointToMaskMap =
870 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
871 auto isthmusTable = [&table,
872 &pointToMaskMap](const FixtureComplex &fc,
873 const FixtureComplex::Cell &c) {
874 return skelWithTable(table, pointToMaskMap, fc, c);
875 };
876 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
877 vc, selectRandom<FixtureComplex>, isthmusTable, 0);
878 CHECK(vc_new.nbCells(3) == 3);
879 }
880}
881
883// Fixture for an X
884struct Fixture_X {
886 // type aliases
888 using Point = DGtal::Z3i::Point;
891
892 using FixtureDigitalTopology = DGtal::Z3i::DT26_6;
893 using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
895 std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
896 using FixtureMap = std::unordered_map<KSpace::Cell, CubicalCellData>;
897 using FixtureComplex = DGtal::VoxelComplex<KSpace>;
898
900 // fixture data
901 FixtureComplex complex_fixture;
902 FixtureDigitalSet set_fixture;
903 KSpace ks_fixture; // needed because ConstAlias in CC constructor.
905
907 // Constructor
909 Fixture_X() : complex_fixture(ks_fixture), set_fixture(create_set()) {
910 create_complex_from_set(set_fixture);
911 }
912
914 // Function members
916 FixtureDigitalSet create_set() {
917 using namespace DGtal;
918
919 Point p1(-16, -16, -16);
920 Point p2(16, 16, 16);
921 Domain domain(p1, p2);
922
923 FixtureDigitalSet a_set(domain);
924 std::vector<Point> center_set;
925 center_set.reserve(9);
926
927 Point c00(0, 0, 0);
928 center_set.push_back(c00);
929 Point c01x(-1, 0, 0);
930 center_set.push_back(c01x);
931 Point c10x(1, 0, 0);
932 center_set.push_back(c10x);
933 Point c02x(-2, 0, 0);
934 center_set.push_back(c02x);
935 Point c20x(2, 0, 0);
936 center_set.push_back(c20x);
937
938 Point c01y(0, -1, 0);
939 center_set.push_back(c01y);
940 Point c10y(0, 1, 0);
941 center_set.push_back(c10y);
942 Point c02y(0, -2, 0);
943 center_set.push_back(c02y);
944 Point c20y(0, 2, 0);
945 center_set.push_back(c20y);
946
947 Point z_pos(0, 0, 1);
948 int branch_length(4);
949 std::vector<Point> diagonals;
950 diagonals.reserve(6);
951 for (const auto &p : center_set) {
952 diagonals.clear();
953 for (int l = 0; l <= branch_length; ++l) {
954 diagonals.push_back({l, l, 0});
955 diagonals.push_back({l, -l, 0});
956 diagonals.push_back({-l, l, 0});
957 diagonals.push_back({-l, -l, 0});
958 for (int z = -1; z <= 1; ++z)
959 for (const auto &d : diagonals)
960 a_set.insert(p + d + (z * z_pos));
961 }
962 }
963
964 return a_set;
965 }
966
967 FixtureComplex &create_complex_from_set(FixtureDigitalSet &input_set) {
968
969 ks_fixture.init(input_set.domain().lowerBound(),
970 input_set.domain().upperBound(), true);
971 complex_fixture = FixtureComplex(ks_fixture);
972 complex_fixture.construct(input_set);
973 return complex_fixture;
974 }
975};
976
977TEST_CASE_METHOD(Fixture_X, "X Thin",
978 "[x][persistence][isthmus][thin][function]") {
979 using namespace DGtal::functions;
980 auto &vc = complex_fixture;
981 auto &ks = vc.space();
982 boost::ignore_unused_variable_warning(ks);
983 bool verbose = true;
984 SECTION(
985 "persistence value of 1 is equivalent to the assymetric algorithm") {
986 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
987 vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>, 1,
988 verbose);
989
990 auto vc_assymetric = asymetricThinningScheme<FixtureComplex>(
991 vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>, true);
992 // 10% tolerance
993 CHECK(vc_assymetric.nbCells(3) ==
994 Approx(vc_new.nbCells(3)).epsilon(0.1));
995 }
996 SECTION("persistence value: infinite is equivalent to ultimate skeleton") {
997 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
998 vc, selectRandom<FixtureComplex>, skelEnd<FixtureComplex>, 100,
999 verbose);
1000 REQUIRE(vc_new.nbCells(3) == 1);
1001 }
1002}
1003
1004TEST_CASE_METHOD(Fixture_X, "X Thin with Isthmus, and tables",
1005 "[x][isthmus][thin][function][table]") {
1006 using namespace DGtal::functions;
1007 auto &vc = complex_fixture;
1008 auto &ks = vc.space();
1009 boost::ignore_unused_variable_warning(ks);
1010 bool verbose = true;
1011
1012 // Disabled to reduce test time
1013 // SECTION("Compute with skelIsthmus") {
1014 // trace.beginBlock("Fixture_X skelIsthmus");
1015 // auto vc_new = asymetricThinningScheme<FixtureComplex>(
1016 // vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>,
1017 // verbose);
1018 // trace.endBlock();
1019 // }
1020 SECTION("Compute with skelWithTable (isIsthmus)") {
1021 trace.beginBlock("Fixture_X skelIsthmus with table");
1022 auto table = *functions::loadTable(isthmusicity::tableIsthmus);
1023 auto pointToMaskMap =
1024 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1025 auto skelWithTableIsthmus =
1026 [&table, &pointToMaskMap](const FixtureComplex &fc,
1027 const FixtureComplex::Cell &c) {
1028 return skelWithTable(table, pointToMaskMap, fc, c);
1029 };
1030
1031 auto vc_new = asymetricThinningScheme<FixtureComplex>(
1032 vc, selectRandom<FixtureComplex>, skelWithTableIsthmus, verbose);
1033
1034 trace.endBlock();
1035 }
1036
1037 SECTION("Compute with skelWithTable (isIsthmus) and empty Object") {
1038 trace.beginBlock("Fixture_X skelIsthmus with table (empty objectSet)");
1039 vc.setSimplicityTable(
1040 functions::loadTable(simplicity::tableSimple26_6));
1041 vc.clear();
1042 auto table = *functions::loadTable(isthmusicity::tableIsthmus);
1043 auto pointToMaskMap =
1044 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1045 auto skelWithTableIsthmus =
1046 [&table, &pointToMaskMap](const FixtureComplex &fc,
1047 const FixtureComplex::Cell &c) {
1048 return skelWithTable(table, pointToMaskMap, fc, c);
1049 };
1050
1051 auto vc_new = asymetricThinningScheme<FixtureComplex>(
1052 vc, selectRandom<FixtureComplex>, skelWithTableIsthmus, verbose);
1053
1054 trace.endBlock();
1055 }
1056}
1057
1059TEST_CASE_METHOD(Fixture_X, "X DistanceMap", "[x][distance][thin]") {
1060 using namespace DGtal::functions;
1061 auto &vc = complex_fixture;
1062 auto &ks = vc.space();
1063 boost::ignore_unused_variable_warning(ks);
1064 using Predicate = Z3i::DigitalSet;
1067 bool verbose = true;
1068 SECTION("Fixture_X Distance Map") {
1069 trace.beginBlock("With a Distance Map");
1070 L3Metric l3;
1071 DT dt(set_fixture.domain(), set_fixture, l3);
1072 // Create wrap around selectMaxValue to use the thinning.
1073 auto selectDistMax = [&dt](const FixtureComplex::Clique &clique) {
1074 return selectMaxValue<DT, FixtureComplex>(dt, clique);
1075 };
1076 SECTION("asymetricThinning"){
1077 auto vc_new = asymetricThinningScheme<FixtureComplex>(
1078 vc, selectDistMax, skelEnd<FixtureComplex>, verbose);
1079 }
1080 SECTION("persistenceThinning"){
1081 auto table = *functions::loadTable(isthmusicity::tableOneIsthmus);
1082 auto pointToMaskMap =
1083 *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1084 auto oneIsthmusTable =
1085 [&table, &pointToMaskMap](const FixtureComplex &fc,
1086 const FixtureComplex::Cell &c) {
1087 return skelWithTable(table, pointToMaskMap, fc, c);
1088 };
1089 int persistence = 3;
1090 auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
1091 vc, selectDistMax, oneIsthmusTable,
1092 persistence, verbose);
1093 // SECTION( "visualize the thining" ){
1094 // int argc(1);
1095 // char** argv(nullptr);
1096 // QApplication app(argc, argv);
1097 // Viewer3D<> viewer(ks_fixture);
1098 // viewer.show();
1099 //
1100 // viewer.setFillColor(Color(200, 200, 200, 100));
1101 // for ( auto it = vc_new.begin(3); it!= vc_new.end(3); ++it )
1102 // viewer << it->first;
1103 //
1104 // // All kspace voxels
1105 // viewer.setFillColor(Color(40, 40, 40, 10));
1106 // for ( auto it = vc.begin(3); it!= vc.end(3); ++it )
1107 // viewer << it->first;
1108 //
1109 // viewer << Viewer3D<>::updateDisplay;
1110 // app.exec();
1111 // }
1112 }
1113
1114 }
1115
1116 trace.endBlock();
1117}
1118
1119// REQUIRE(vc_new.nbCells(3) == 38);
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Aim: Represents a digital topology as a couple of adjacency relations.
TForegroundAdjacency ForegroundAdjacency
TBackgroundAdjacency BackgroundAdjacency
Aim: Implementation of the linear in time distance transformation for separable metrics.
Aim: implements separable l_p metrics with exact predicates.
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,...
const Point & lowerBound() const
Return the lower bound for digital points in this space.
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
const Point & upperBound() const
Return the upper bound for digital points in this space.
AnyCellCollection< Cell > Cells
Aim: An object (or digital object) represents a set in some digital space associated with a digital t...
Definition: Object.h:120
void beginBlock(const std::string &keyword="")
double endBlock()
This class represents a voxel complex living in some Khalimsky space. Voxel complexes are derived fro...
Definition: VoxelComplex.h:91
HyperRectDomain< Space > Domain
Definition: StdDefs.h:172
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
Space::Point Point
Definition: StdDefs.h:168
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
Definition: StdDefs.h:173
DigitalTopology< Adj26, Adj6 > DT26_6
Definition: StdDefs.h:167
functions namespace gathers all DGtal functionsxs.
DGtal::CountedPtr< boost::dynamic_bitset<> > loadTable(const std::string &input_filename, const unsigned int known_size, const bool compressed=true)
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
bool isEqual(Container1 &c1, Container2 &c2)
Domain domain
SECTION("Testing constant forward iterators")
REQUIRE(domain.isInside(aPoint))
TEST_CASE_METHOD(Fixture_complex_diamond, "insertVoxel", "[insert][close]")