DGtal  1.0.0
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 
52 using namespace std;
53 using namespace DGtal;
54 
56 // Fixture for object from a diamond set
58 struct Fixture_complex_diamond {
60  // type aliases
62  using Point = DGtal::Z3i::Point;
63  using Domain = DGtal::Z3i::Domain;
64  using KSpace = DGtal::Z3i::KSpace;
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 
120 TEST_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 
180 TEST_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 
248 TEST_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 (auto &&n : std::vector<int>(10))
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 
288 TEST_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 (auto &&n : std::vector<int>(10))
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 
313 TEST_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 (auto &&n : std::vector<int>(10))
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 
341 TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_2", "[clique]") {
342  auto &vc = complex_fixture;
343  auto itc = vc.begin(3);
344  for (auto &&n : std::vector<int>(10))
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 
380 TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_1", "[clique]") {
381  auto &vc = complex_fixture;
382  auto itc = vc.begin(3);
383  for (auto &&n : std::vector<int>(10))
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 
402 TEST_CASE_METHOD(Fixture_complex_diamond, "Cliques Masks K_0", "[clique]") {
403  auto &vc = complex_fixture;
404  auto itc = vc.begin(3);
405  for (auto &&n : std::vector<int>(10))
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 
425 TEST_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
448 struct Fixture_complex_fig4 {
450  // type aliases
452  using Point = DGtal::Z3i::Point;
453  using Domain = DGtal::Z3i::Domain;
454  using KSpace = DGtal::Z3i::KSpace;
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 };
524 TEST_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 */
544 TEST_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
563  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
565  std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
569  auto domain = DGtal::Z3i::Domain(ks_fixture.lowerBound(), ks_fixture.upperBound());
570  DigitalTopology topo(
571  adjF, adjB, DGtal::DigitalTopologyProperties::JORDAN_DT);
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
613 struct Fixture_isthmus {
615  // type aliases
617  using Point = DGtal::Z3i::Point;
618  using Domain = DGtal::Z3i::Domain;
619  using KSpace = DGtal::Z3i::KSpace;
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 
716 TEST_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  Point x(-1, 4, 0);
722  set_fixture.erase(x);
723  // Delete one point and reconstruct complex.
724  /* obj_fixture.pointSet().erase(x); */
725  vc.clear();
726  vc.construct(set_fixture);
727  SECTION("with skelUltimate") {
728  auto vc_new = asymetricThinningScheme<FixtureComplex>(
729  vc, selectFirst<FixtureComplex>, skelUltimate<FixtureComplex>);
730  CHECK(vc_new.nbCells(3) == 2);
731  }
732 }
733 
734 
735 TEST_CASE_METHOD(Fixture_isthmus, "Check isthmus", "[isthmus][function]") {
736  auto &vc = complex_fixture;
737  using namespace DGtal::functions;
738  auto &ks = vc.space();
739  Point x(-1, 4, 0); // oneIsthmus
740  Point y(-1, 1, 0); // twoIsthmus
741 
742  SECTION("checking one isthmus") {
743  std::set<Point> isthmus;
744  for (const auto &p : set_fixture) {
745  if (oneIsthmus(vc, ks.uSpel(p)))
746  isthmus.insert(p);
747  }
748  CHECK(isthmus.size() == 1);
749  auto xit = isthmus.find(x);
750  CHECK(*xit == x);
751  }
752 
753  SECTION("checking 2 isthmus") {
754  std::set<Point> isthmus;
755  for (const auto &p : set_fixture) {
756  if (twoIsthmus(vc, ks.uSpel(p)))
757  isthmus.insert(p);
758  }
759  CHECK(isthmus.size() == 1);
760  auto yit = isthmus.find(y);
761  CHECK(*yit == y);
762  }
763 
764  SECTION("checking 1 and 2 isthmus") {
765  std::set<Point> isthmus;
766  for (const auto &p : set_fixture) {
767  if (skelIsthmus(vc, ks.uSpel(p)))
768  isthmus.insert(p);
769  }
770  CHECK(isthmus.size() == 2);
771  auto xit = isthmus.find(x);
772  auto yit = isthmus.find(y);
773  CHECK(xit != isthmus.end());
774  CHECK(yit != isthmus.end());
775  }
776 }
777 TEST_CASE_METHOD(Fixture_isthmus, "Thin complex", "[isthmus][thin][function]") {
778  using namespace DGtal::functions;
779  auto &vc = complex_fixture;
780  auto &ks = vc.space();
781  SECTION("with skelUltimate") {
782  auto vc_new = asymetricThinningScheme<FixtureComplex>(
783  vc, selectFirst<FixtureComplex>, skelUltimate<FixtureComplex>);
784  CHECK(vc_new.nbCells(3) == 1);
785  }
786  SECTION("with skelEnd") {
787  auto vc_new = asymetricThinningScheme<FixtureComplex>(
788  vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>);
789  CHECK(vc_new.nbCells(3) == 5);
790  }
791  SECTION("with oneIsthmus") {
792  auto vc_new = asymetricThinningScheme<FixtureComplex>(
793  vc, selectRandom<FixtureComplex>, oneIsthmus<FixtureComplex>);
794  CHECK(vc_new.nbCells(3) == 3);
795  }
796  SECTION("with twoIsthmus") {
797  auto vc_new = asymetricThinningScheme<FixtureComplex>(
798  vc, selectRandom<FixtureComplex>, twoIsthmus<FixtureComplex>);
799  CHECK(vc_new.nbCells(3) == 1);
800  }
801  SECTION("with skelIsthmus") {
802  auto vc_new = asymetricThinningScheme<FixtureComplex>(
803  vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>);
804  CHECK(vc_new.nbCells(3) == 3);
805  }
806 }
807 //
808 TEST_CASE_METHOD(Fixture_isthmus, "Persistence thin",
809  "[persistence][isthmus][thin][function]") {
810  using namespace DGtal::functions;
811  auto &vc = complex_fixture;
812  auto &ks = vc.space();
813  SECTION("with skelUltimate") {
814  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
815  vc, selectRandom<FixtureComplex>, skelUltimate<FixtureComplex>, 0);
816  CHECK(vc_new.nbCells(3) == 1);
817  }
818  SECTION("with skelEnd") {
819  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
820  vc, selectRandom<FixtureComplex>, skelEnd<FixtureComplex>, 0);
821  CHECK(vc_new.nbCells(3) == 5);
822  }
823  SECTION("with oneIsthmus") {
824  // Not using LUT skel function (slow):
825  //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
826  // vc, selectRandom<FixtureComplex>, oneIsthmus<FixtureComplex>, 0);
827  //
828  // with LUT:
829  auto table = *functions::loadTable(isthmusicity::tableOneIsthmus);
830  auto pointToMaskMap =
831  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
832  auto oneIsthmusTable =
833  [&table, &pointToMaskMap](const FixtureComplex &fc,
834  const FixtureComplex::Cell &c) {
835  return skelWithTable(table, pointToMaskMap, fc, c);
836  };
837  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
838  vc, selectRandom<FixtureComplex>, oneIsthmusTable, 0);
839  CHECK(vc_new.nbCells(3) == 3);
840  }
841  SECTION("with twoIsthmus") {
842  // Not using LUT skel function (slow):
843  //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
844  // vc, selectRandom<FixtureComplex>, twoIsthmus<FixtureComplex>, 0);
845  //
846  // with LUT:
847  auto table = *functions::loadTable(isthmusicity::tableTwoIsthmus);
848  auto pointToMaskMap =
849  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
850  auto twoIsthmusTable =
851  [&table, &pointToMaskMap](const FixtureComplex &fc,
852  const FixtureComplex::Cell &c) {
853  return skelWithTable(table, pointToMaskMap, fc, c);
854  };
855  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
856  vc, selectRandom<FixtureComplex>, twoIsthmusTable, 0);
857  CHECK(vc_new.nbCells(3) == 1);
858  }
859  SECTION("with skelIsthmus") {
860  // Not using LUT skel function (slow):
861  //auto vc_new = persistenceAsymetricThinningScheme< FixtureComplex >(
862  // vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>, 0);
863  //
864  // with LUT:
865  auto table = *functions::loadTable(isthmusicity::tableIsthmus);
866  auto pointToMaskMap =
867  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
868  auto isthmusTable = [&table,
869  &pointToMaskMap](const FixtureComplex &fc,
870  const FixtureComplex::Cell &c) {
871  return skelWithTable(table, pointToMaskMap, fc, c);
872  };
873  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
874  vc, selectRandom<FixtureComplex>, isthmusTable, 0);
875  CHECK(vc_new.nbCells(3) == 3);
876  }
877 }
878 
880 // Fixture for a thick X
881 struct Fixture_X {
883  // type aliases
885  using Point = DGtal::Z3i::Point;
886  using Domain = DGtal::Z3i::Domain;
887  using KSpace = DGtal::Z3i::KSpace;
888 
889  using FixtureDigitalTopology = DGtal::Z3i::DT26_6;
890  using FixtureDigitalSet = DGtal::DigitalSetByAssociativeContainer<
892  std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
893  using FixtureMap = std::unordered_map<KSpace::Cell, CubicalCellData>;
894  using FixtureComplex = DGtal::VoxelComplex<KSpace>;
895 
897  // fixture data
898  FixtureComplex complex_fixture;
899  FixtureDigitalSet set_fixture;
900  KSpace ks_fixture; // needed because ConstAlias in CC constructor.
902 
904  // Constructor
906  Fixture_X() : complex_fixture(ks_fixture), set_fixture(create_set()) {
907  create_complex_from_set(set_fixture);
908  }
909 
911  // Function members
913  FixtureDigitalSet create_set() {
914  using namespace DGtal;
915 
916  Point p1(-30, -30, -30);
917  Point p2(30, 30, 30);
918  Domain domain(p1, p2);
919 
920  FixtureDigitalSet a_set(domain);
921  std::vector<Point> center_set;
922  center_set.reserve(9);
923 
924  Point c00(0, 0, 0);
925  center_set.push_back(c00);
926  Point c01x(-1, 0, 0);
927  center_set.push_back(c01x);
928  Point c10x(1, 0, 0);
929  center_set.push_back(c10x);
930  Point c02x(-2, 0, 0);
931  center_set.push_back(c02x);
932  Point c20x(2, 0, 0);
933  center_set.push_back(c20x);
934 
935  Point c01y(0, -1, 0);
936  center_set.push_back(c01y);
937  Point c10y(0, 1, 0);
938  center_set.push_back(c10y);
939  Point c02y(0, -2, 0);
940  center_set.push_back(c02y);
941  Point c20y(0, 2, 0);
942  center_set.push_back(c20y);
943 
944  Point z_pos(0, 0, 1);
945  int branch_length(10);
946  std::vector<Point> diagonals;
947  diagonals.reserve(6);
948  for (const auto &p : center_set) {
949  diagonals.clear();
950  for (int l = 0; l <= branch_length; ++l) {
951  diagonals.push_back({l, l, 0});
952  diagonals.push_back({l, -l, 0});
953  diagonals.push_back({-l, l, 0});
954  diagonals.push_back({-l, -l, 0});
955  for (int z = -2; z <= 2; ++z)
956  for (const auto &d : diagonals)
957  a_set.insert(p + d + (z * z_pos));
958  }
959  }
960 
961  return a_set;
962  }
963 
964  FixtureComplex &create_complex_from_set(FixtureDigitalSet &input_set) {
965 
966  ks_fixture.init(input_set.domain().lowerBound(),
967  input_set.domain().upperBound(), true);
968  complex_fixture = FixtureComplex(ks_fixture);
969  complex_fixture.construct(input_set);
970  return complex_fixture;
971  }
972 };
973 
974 TEST_CASE_METHOD(Fixture_X, "X Thin",
975  "[x][persistence][isthmus][thin][function]") {
976  using namespace DGtal::functions;
977  auto &vc = complex_fixture;
978  auto &ks = vc.space();
979  bool verbose = true;
980  SECTION(
981  "persistence value of 1 is equivalent to the assymetric algorithm") {
982  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
983  vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>, 1,
984  verbose);
985 
986  auto vc_assymetric = asymetricThinningScheme<FixtureComplex>(
987  vc, selectFirst<FixtureComplex>, skelEnd<FixtureComplex>, true);
988  // 10% tolerance
989  CHECK(vc_assymetric.nbCells(3) ==
990  Approx(vc_new.nbCells(3)).epsilon(0.1));
991  }
992  SECTION("persistence value: infinite is equivalent to ultimate skeleton") {
993  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
994  vc, selectRandom<FixtureComplex>, skelEnd<FixtureComplex>, 100,
995  verbose);
996  REQUIRE(vc_new.nbCells(3) == 1);
997  }
998 }
999 
1000 TEST_CASE_METHOD(Fixture_X, "X Thin with Isthmus, and tables",
1001  "[x][isthmus][thin][function][table]") {
1002  using namespace DGtal::functions;
1003  auto &vc = complex_fixture;
1004  auto &ks = vc.space();
1005  bool verbose = true;
1006 
1007  // Disabled to reduce test time
1008  // SECTION("Compute with skelIsthmus") {
1009  // trace.beginBlock("Fixture_X skelIsthmus");
1010  // auto vc_new = asymetricThinningScheme<FixtureComplex>(
1011  // vc, selectRandom<FixtureComplex>, skelIsthmus<FixtureComplex>,
1012  // verbose);
1013  // trace.endBlock();
1014  // }
1015  SECTION("Compute with skelWithTable (isIsthmus)") {
1016  trace.beginBlock("Fixture_X skelIsthmus with table");
1017  auto table = *functions::loadTable(isthmusicity::tableIsthmus);
1018  auto pointToMaskMap =
1019  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1020  auto skelWithTableIsthmus =
1021  [&table, &pointToMaskMap](const FixtureComplex &fc,
1022  const FixtureComplex::Cell &c) {
1023  return skelWithTable(table, pointToMaskMap, fc, c);
1024  };
1025 
1026  auto vc_new = asymetricThinningScheme<FixtureComplex>(
1027  vc, selectRandom<FixtureComplex>, skelWithTableIsthmus, verbose);
1028 
1029  trace.endBlock();
1030  }
1031 
1032  SECTION("Compute with skelWithTable (isIsthmus) and empty Object") {
1033  trace.beginBlock("Fixture_X skelIsthmus with table (empty objectSet)");
1034  vc.setSimplicityTable(
1035  functions::loadTable(simplicity::tableSimple26_6));
1036  vc.clear();
1037  auto table = *functions::loadTable(isthmusicity::tableIsthmus);
1038  auto pointToMaskMap =
1039  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1040  auto skelWithTableIsthmus =
1041  [&table, &pointToMaskMap](const FixtureComplex &fc,
1042  const FixtureComplex::Cell &c) {
1043  return skelWithTable(table, pointToMaskMap, fc, c);
1044  };
1045 
1046  auto vc_new = asymetricThinningScheme<FixtureComplex>(
1047  vc, selectRandom<FixtureComplex>, skelWithTableIsthmus, verbose);
1048 
1049  trace.endBlock();
1050  }
1051 }
1052 
1054 TEST_CASE_METHOD(Fixture_X, "X DistanceMap", "[x][distance][thin]") {
1055  using namespace DGtal::functions;
1056  auto &vc = complex_fixture;
1057  auto &ks = vc.space();
1058  using Predicate = Z3i::DigitalSet;
1061  bool verbose = true;
1062  SECTION("Fixture_X Distance Map") {
1063  trace.beginBlock("With a Distance Map");
1064  L3Metric l3;
1065  DT dt(set_fixture.domain(), set_fixture, l3);
1066  // Create wrap around selectMaxValue to use the thinning.
1067  auto selectDistMax = [&dt](const FixtureComplex::Clique &clique) {
1068  return selectMaxValue<DT, FixtureComplex>(dt, clique);
1069  };
1070  SECTION("asymetricThinning"){
1071  auto vc_new = asymetricThinningScheme<FixtureComplex>(
1072  vc, selectDistMax, skelEnd<FixtureComplex>, verbose);
1073  }
1074  SECTION("persistenceThinning"){
1075  auto table = *functions::loadTable(isthmusicity::tableOneIsthmus);
1076  auto pointToMaskMap =
1077  *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
1078  auto oneIsthmusTable =
1079  [&table, &pointToMaskMap](const FixtureComplex &fc,
1080  const FixtureComplex::Cell &c) {
1081  return skelWithTable(table, pointToMaskMap, fc, c);
1082  };
1083  int persistence = 3;
1084  auto vc_new = persistenceAsymetricThinningScheme<FixtureComplex>(
1085  vc, selectDistMax, oneIsthmusTable,
1086  persistence, verbose);
1087  // SECTION( "visualize the thining" ){
1088  // int argc(1);
1089  // char** argv(nullptr);
1090  // QApplication app(argc, argv);
1091  // Viewer3D<> viewer(ks_fixture);
1092  // viewer.show();
1093  //
1094  // viewer.setFillColor(Color(200, 200, 200, 100));
1095  // for ( auto it = vc_new.begin(3); it!= vc_new.end(3); ++it )
1096  // viewer << it->first;
1097  //
1098  // // All kspace voxels
1099  // viewer.setFillColor(Color(40, 40, 40, 10));
1100  // for ( auto it = vc.begin(3); it!= vc.end(3); ++it )
1101  // viewer << it->first;
1102  //
1103  // viewer << Viewer3D<>::updateDisplay;
1104  // app.exec();
1105  // }
1106  }
1107 
1108  }
1109 
1110  trace.endBlock();
1111 }
1112 
1113 // REQUIRE(vc_new.nbCells(3) == 38);
Aim: An object (or digital object) represents a set in some digital space associated with a digital t...
Definition: Object.h:119
void beginBlock(const std::string &keyword="")
const ConstIterator & end() const
CAPTURE(thicknessHV)
Aim: Implementation of the linear in time distance transformation for separable metrics.
Trace trace
Definition: Common.h:144
HyperRectDomain< Space > Domain
Definition: StdDefs.h:172
const Point & upperBound() const
Return the upper bound for digital points in this space.
TForegroundAdjacency ForegroundAdjacency
bool twoIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
double endBlock()
bool oneIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
Domain domain
bool isOneSurface(const TObject &small_obj)
bool skelIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
Aim: implements separable l_p metrics with exact predicates.
REQUIRE(domain.isInside(aPoint))
TBackgroundAdjacency BackgroundAdjacency
Aim: Represents a digital topology as a couple of adjacency relations.
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
bool isZeroSurface(const TObject &small_obj)
const ConstIterator & begin() const
DigitalTopology< Adj26, Adj6 > DT26_6
Definition: StdDefs.h:167
TEST_CASE_METHOD(Fixture_complex_diamond, "insertVoxel", "[insert][close]")
const Point & lowerBound() const
Return the lower bound for digital points in this space.
DGtal is the top-level namespace which contains all DGtal functions and types.
MyPointD Point
Definition: testClone2.cpp:383
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
bool isEqual(Container1 &c1, Container2 &c2)
Space::Point Point
Definition: StdDefs.h:168
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
functions namespace gathers all DGtal functionsxs.
Definition: SetFunctions.h:768
bool skelWithTable(const boost::dynamic_bitset<> &table, const std::unordered_map< typename TComplex::Point, unsigned int > &pointToMaskMap, const TComplex &vc, const typename TComplex::Cell &cell)
KSpace::Cell Cell
SECTION("Testing constant forward iterators")
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
AnyCellCollection< Cell > Cells