DGtal  1.3.beta
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 (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 
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 (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 
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 (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 
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 (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 
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 (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 
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 (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 
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
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  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 
736 TEST_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 }
778 TEST_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 //
810 TEST_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
884 struct Fixture_X {
886  // type aliases
888  using Point = DGtal::Z3i::Point;
889  using Domain = DGtal::Z3i::Domain;
890  using KSpace = DGtal::Z3i::KSpace;
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 
977 TEST_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 
1004 TEST_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 
1059 TEST_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);
DGtal::DigitalTopology::BackgroundAdjacency
TBackgroundAdjacency BackgroundAdjacency
Definition: DigitalTopology.h:100
DGtal::KhalimskySpaceND::upperBound
const Point & upperBound() const
Return the upper bound for digital points in this space.
DGtal::DigitalTopology
Aim: Represents a digital topology as a couple of adjacency relations.
Definition: DigitalTopology.h:95
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::Z3i::DT26_6
DigitalTopology< Adj26, Adj6 > DT26_6
Definition: StdDefs.h:167
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::KhalimskySpaceND::Cells
AnyCellCollection< Cell > Cells
Definition: KhalimskySpaceND.h:439
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::Z3i::Point
Space::Point Point
Definition: StdDefs.h:168
DGtal::Z3i::KSpace
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::KhalimskySpaceND::lowerBound
const Point & lowerBound() const
Return the lower bound for digital points in this space.
DGtal::VoxelComplex< KSpace >
REQUIRE
REQUIRE(domain.isInside(aPoint))
CAPTURE
CAPTURE(thicknessHV)
DGtal::functions
functions namespace gathers all DGtal functionsxs.
TEST_CASE_METHOD
TEST_CASE_METHOD(Fixture_complex_diamond, "insertVoxel", "[insert][close]")
Definition: testVoxelComplex.cpp:120
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::functions::skelIsthmus
bool skelIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
DGtal::functions::skelWithTable
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)
DGtal::Z3i::Domain
HyperRectDomain< Space > Domain
Definition: StdDefs.h:172
DGtal::HyperRectDomain::end
const ConstIterator & end() const
Definition: HyperRectDomain.h:201
DGtal::functions::oneIsthmus
bool oneIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
DGtal::functions::isZeroSurface
bool isZeroSurface(const TObject &small_obj)
DGtal::ExactPredicateLpSeparableMetric
Aim: implements separable l_p metrics with exact predicates.
Definition: ExactPredicateLpSeparableMetric.h:87
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::DigitalTopology::ForegroundAdjacency
TForegroundAdjacency ForegroundAdjacency
Definition: DigitalTopology.h:99
SECTION
SECTION("Testing constant forward iterators")
Definition: testSimpleRandomAccessRangeFromPoint.cpp:66
isEqual
bool isEqual(Container1 &c1, Container2 &c2)
Definition: testIndexedListWithBlocks.cpp:45
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
DGtal::functions::twoIsthmus
bool twoIsthmus(const TComplex &vc, const typename TComplex::Cell &cell)
DGtal::Object
Aim: An object (or digital object) represents a set in some digital space associated with a digital t...
Definition: Object.h:119
DGtal::DistanceTransformation
Aim: Implementation of the linear in time distance transformation for separable metrics.
Definition: DistanceTransformation.h:98
DGtal::functions::isOneSurface
bool isOneSurface(const TObject &small_obj)
DGtal::HyperRectDomain::begin
const ConstIterator & begin() const
Definition: HyperRectDomain.h:176
Point
MyPointD Point
Definition: testClone2.cpp:383
DigitalSet
Z2i::DigitalSet DigitalSet
Definition: testVoronoiMapComplete.cpp:41
DGtal::DigitalSetByAssociativeContainer
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Definition: DigitalSetByAssociativeContainer.h:89
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64