DGtal  1.0.0
VoxelComplex.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file VoxelComplex.ih
19  * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
20  * Insitute of Fundamental Sciences, Massey University, New Zealand.
21  *
22  * @date 2018/01/01
23  *
24  * Implementation of inline methods defined in VoxelComplex.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 //////////////////////////////////////////////////////////////////////////////
30 #include <DGtal/graph/ObjectBoostGraphInterface.h>
31 #include <DGtal/topology/CubicalComplexFunctions.h>
32 #include <DGtal/topology/NeighborhoodConfigurations.h>
33 #include <boost/graph/connected_components.hpp>
34 #include <boost/graph/filtered_graph.hpp>
35 #include <boost/property_map/property_map.hpp>
36 #include <iostream>
37 #ifdef WITH_OPENMP
38 // #include <experimental/algorithm>
39 #include <omp.h>
40 #endif
41 //////////////////////////////////////////////////////////////////////////////
42 // Default constructor:
43 template <typename TKSpace, typename TCellContainer>
44 inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex()
45  : Parent(),
46  myTablePtr(nullptr), myPointToMaskPtr(nullptr),
47  myIsTableLoaded(false) {}
48 
49 // Copy constructor:
50 template <typename TKSpace, typename TCellContainer>
51 inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex(
52  const VoxelComplex &other)
53  : Parent(other),
54  myTablePtr(other.myTablePtr),
55  myPointToMaskPtr(other.myPointToMaskPtr),
56  myIsTableLoaded(other.myIsTableLoaded) {}
57 
58 ///////////////////////////////////////////////////////////////////////////////
59 // IMPLEMENTATION of inline methods.
60 ///////////////////////////////////////////////////////////////////////////////
61 //-----------------------------------------------------------------------------
62 template <typename TKSpace, typename TCellContainer>
63 inline DGtal::VoxelComplex<TKSpace, TCellContainer> &
64 DGtal::VoxelComplex<TKSpace, TCellContainer>::
65 operator=(const Self &other) {
66  if (this != &other) {
67  this->myKSpace = other.myKSpace;
68  this->myCells = other.myCells;
69  myTablePtr = other.myTablePtr;
70  myPointToMaskPtr = other.myPointToMaskPtr;
71  myIsTableLoaded = other.myIsTableLoaded;
72  }
73  return *this;
74 }
75 //---------------------------------------------------------------------------
76 ///////////////////////////////////////////////////////////////////////////////
77 // Voxel methods :
78 ///////////////////////////////////////////////////////////////////////////////
79 
80 ///////////////////////////////////////////////////////////////////////////////
81 // Interface - Voxel :
82 //---------------------------------------------------------------------------
83 // template <typename TKSpace, typename TCellContainer>
84 // const typename DGtal::VoxelComplex<TKSpace,
85 // TCellContainer>::Object::Point &
86 // DGtal::VoxelComplex<TKSpace, TCellContainer>::objPointFromVoxel(
87 // const Cell &voxel) const {
88 // ASSERT(isSpel(voxel) == true);
89 // ASSERT(this->belongs(voxel));
90 // auto &ks = this->space();
91 // return *(myObject.pointSet().find(ks.uCoords(voxel)));
92 // }
93 
94 
95 //-----------------------------------------------------------------------------
96 
97 template <typename TKSpace, typename TCellContainer>
98 template <typename TDigitalSet>
99 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::construct(
100  const TDigitalSet &input_set) {
101  Parent::construct(input_set);
102 }
103 
104 template <typename TKSpace, typename TCellContainer>
105 template <typename TDigitalSet>
106 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::construct(
107  const TDigitalSet &input_set,
108  const Alias<ConfigMap> input_table) {
109  Parent::construct(input_set);
110  setSimplicityTable(input_table);
111 }
112 
113 template <typename TKSpace, typename TCellContainer>
114 void DGtal::VoxelComplex<TKSpace, TCellContainer>::setSimplicityTable(
115  const Alias<ConfigMap> input_table) {
116  this->myTablePtr = input_table;
117  this->myPointToMaskPtr =
118  functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
119  this->myIsTableLoaded = true;
120 }
121 
122 template <typename TKSpace, typename TCellContainer>
123 void DGtal::VoxelComplex<TKSpace, TCellContainer>::copySimplicityTable(
124  const Self & other) {
125  myTablePtr = other.myTablePtr;
126  myPointToMaskPtr = other.myPointToMaskPtr;
127  myIsTableLoaded = other.myIsTableLoaded;
128 }
129 
130 template <typename TKSpace, typename TCellContainer>
131 const typename DGtal::VoxelComplex<TKSpace,
132  TCellContainer>::ConfigMap &
133 DGtal::VoxelComplex<TKSpace, TCellContainer>::table() const {
134  return *myTablePtr;
135 }
136 
137 template <typename TKSpace, typename TCellContainer>
138 const bool &
139 DGtal::VoxelComplex<TKSpace, TCellContainer>::isTableLoaded() const {
140  return myIsTableLoaded;
141 }
142 
143 template <typename TKSpace, typename TCellContainer>
144 const typename DGtal::VoxelComplex<TKSpace,
145  TCellContainer>::PointToMaskMap &
146 DGtal::VoxelComplex<TKSpace, TCellContainer>::pointToMask() const {
147  return *myPointToMaskPtr;
148 }
149 //---------------------------------------------------------------------------
150 template <typename TKSpace, typename TCellContainer>
151 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::voxelClose(
152  const Cell &kcell) {
153  auto &ks = this->space();
154  ASSERT(ks.uDim(kcell) == 3);
155  Dimension l = 2;
156  auto direct_faces = ks.uLowerIncident(kcell);
157  for (typename Cells::const_iterator cells_it = direct_faces.begin(),
158  cells_it_end = direct_faces.end();
159  cells_it != cells_it_end; ++cells_it) {
160  this->insertCell(l, *cells_it);
161  }
162  cellsClose(l, direct_faces);
163 }
164 //---------------------------------------------------------------------------
165 template <typename TKSpace, typename TCellContainer>
166 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::cellsClose(
167  Dimension k, const Cells &cells) {
168  if (k <= 0)
169  return;
170  if (cells.size() == 0)
171  return;
172  auto &ks = this->space();
173  Dimension l = k - 1;
174  for (auto const &kcell : cells) {
175  auto direct_faces = ks.uLowerIncident(kcell);
176  for (typename Cells::const_iterator cells_it = direct_faces.begin(),
177  cells_it_end = direct_faces.end();
178  cells_it != cells_it_end; ++cells_it) {
179  this->insertCell(l, *cells_it);
180  }
181  cellsClose(l, direct_faces);
182  }
183 }
184 //---------------------------------------------------------------------------
185 template <typename TKSpace, typename TCellContainer>
186 inline void
187 DGtal::VoxelComplex<TKSpace, TCellContainer>::insertVoxelCell(
188  const Cell &kcell, const bool &close_it, const Data &data) {
189  auto &ks = this->space();
190  ASSERT(ks.uDim(kcell) == 3);
191  this->insertCell(3, kcell, data);
192  if (close_it)
193  voxelClose(kcell);
194 }
195 
196 //---------------------------------------------------------------------------
197 template <typename TKSpace, typename TCellContainer>
198 inline void
199 DGtal::VoxelComplex<TKSpace, TCellContainer>::insertVoxelPoint(
200  const Point &dig_point, const bool &close_it, const Data &data) {
201  auto &ks = this->space();
202  insertVoxelCell(ks.uSpel(dig_point), close_it, data);
203 }
204 
205 //---------------------------------------------------------------------------
206 template <typename TKSpace, typename TCellContainer>
207 template <typename TDigitalSet>
208 inline void
209 DGtal::VoxelComplex<TKSpace, TCellContainer>::dumpVoxels(TDigitalSet & in_out_set) {
210  auto &ks = this->space();
211  for (auto it = this->begin(3), itE = this->end(3) ; it != itE ; ++it ){
212  in_out_set.insertNew(ks.uCoords(it->first));
213  }
214 }
215 //---------------------------------------------------------------------------
216 template <typename TKSpace, typename TCellContainer>
217 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::clear() {
218  for (Dimension d = 0; d <= dimension; ++d)
219  this->Parent::clear(d);
220 }
221 
222 //---------------------------------------------------------------------------
223 
224 template <typename TKSpace, typename TCellContainer>
225 void DGtal::VoxelComplex<TKSpace, TCellContainer>::pointelsFromCell(
226  std::set<Cell> &pointels_out, const Cell &input_cell) const {
227  auto input_dim = this->space().uDim(input_cell);
228  if (input_dim == 0) {
229  pointels_out.emplace(input_cell);
230  return;
231  } else {
232  auto ufaces = this->space().uFaces(input_cell);
233  for (auto &&f : ufaces)
234  this->pointelsFromCell(pointels_out, f);
235  }
236 }
237 
238 //---------------------------------------------------------------------------
239 template <typename TKSpace, typename TCellContainer>
240 void DGtal::VoxelComplex<TKSpace, TCellContainer>::spelsFromCell(
241  std::set<Cell> &spels_out, const Cell &input_cell) const {
242  auto input_dim = this->space().uDim(input_cell);
243  if (input_dim == this->dimension) {
244  if (this->belongs(input_cell))
245  spels_out.emplace(input_cell);
246  return;
247  }
248  auto co_faces = this->space().uCoFaces(input_cell);
249  for (auto &&f : co_faces) {
250  auto f_dim = this->space().uDim(f);
251  if (f_dim >= input_dim)
252  this->spelsFromCell(spels_out, f);
253  }
254 }
255 
256 //---------------------------------------------------------------------------
257 template <typename TKSpace, typename TCellContainer>
258 typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique
259 DGtal::VoxelComplex<TKSpace, TCellContainer>::Kneighborhood(
260  const Cell &input_cell) const {
261 
262  auto spels_out = neighborhoodVoxels(input_cell);
263 
264  auto &ks = this->space();
265  Clique clique(ks);
266  for (const auto &v : spels_out)
267  clique.insertCell(v);
268 
269  return clique;
270 }
271 
272 template <typename TKSpace, typename TCellContainer>
273 std::set<typename TKSpace::Cell>
274 DGtal::VoxelComplex<TKSpace, TCellContainer>::neighborhoodVoxels(
275  const Cell &input_spel) const {
276 
277  std::set<Cell> pointels_out;
278  std::set<Cell> spels_out;
279  pointelsFromCell(pointels_out, input_spel);
280  for (const auto &p : pointels_out)
281  spelsFromCell(spels_out, p);
282  return spels_out;
283 }
284 
285 template <typename TKSpace, typename TCellContainer>
286 std::set<typename TKSpace::Cell>
287 DGtal::VoxelComplex<TKSpace, TCellContainer>::properNeighborhoodVoxels(
288  const Cell &input_spel) const {
289 
290  auto spels_out = neighborhoodVoxels(input_spel);
291  auto search = spels_out.find(input_spel);
292  if (search != spels_out.end()) {
293  spels_out.erase(search);
294  }
295  return spels_out;
296 }
297 
298 //---------------------------------------------------------------------------
299 template <typename TKSpace, typename TCellContainer>
300 bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSpel(
301  const Cell &b) const {
302  return (this->space().uDim(b) == this->space().DIM);
303 }
304 
305 //---------------------------------------------------------------------------
306 template <typename TKSpace, typename TCellContainer>
307 typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Cell
308 DGtal::VoxelComplex<TKSpace,
309  TCellContainer>::surfelBetweenAdjacentSpels(const Cell &A,
310  const Cell &B)
311  const {
312  ASSERT(isSpel(A) == true);
313  ASSERT(isSpel(B) == true);
314  auto &ks = this->space();
315  // Digital coordinates
316  auto &&orientation_BA = ks.uCoords(B) - ks.uCoords(A);
317  ASSERT(orientation_BA.norm1() == 1);
318  // Khalimsky Coordinates
319  return ks.uCell(A.preCell().coordinates + orientation_BA);
320 }
321 //---------------------------------------------------------------------------
322 ///////////////////////////////////////////////////////////////////////////////
323 // Cliques
324 template <typename TKSpace, typename TCellContainer>
325 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
326  TCellContainer>::Clique>
327 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(
328  const typename KSpace::Point &A,
329  const typename KSpace::Point &B,
330  bool verbose) const {
331  auto &ks = this->space();
332  auto orientation_vector = B - A;
333  ASSERT(orientation_vector.norm1() == 1);
334  int direction{-1};
335  for (auto i = 0; i != ks.DIM; ++i)
336  if (orientation_vector[i] == 1 || orientation_vector[i] == -1) {
337  direction = i;
338  break;
339  }
340 
341  Point right, up;
342  if (direction == 0) {
343  right = {0, 0, 1};
344  up = {0, 1, 0};
345  } else if (direction == 1) {
346  right = {1, 0, 0};
347  up = {0, 0, 1};
348  } else {
349  right = {0, 1, 0};
350  up = {1, 0, 0};
351  }
352 
353  Cell x0 = ks.uSpel(A + right);
354  Cell x4 = ks.uSpel(A - right);
355  Cell x2 = ks.uSpel(A + up);
356  Cell x6 = ks.uSpel(A - up);
357 
358  Cell x1 = ks.uSpel(A + up + right);
359  Cell x3 = ks.uSpel(A + up - right);
360  Cell x7 = ks.uSpel(A - up + right);
361  Cell x5 = ks.uSpel(A - up - right);
362 
363  Cell y0 = ks.uSpel(B + right);
364  Cell y4 = ks.uSpel(B - right);
365  Cell y2 = ks.uSpel(B + up);
366  Cell y6 = ks.uSpel(B - up);
367 
368  Cell y1 = ks.uSpel(B + up + right);
369  Cell y3 = ks.uSpel(B + up - right);
370  Cell y7 = ks.uSpel(B - up + right);
371  Cell y5 = ks.uSpel(B - up - right);
372 
373  auto bx0 = this->belongs(KSpace::DIM, x0);
374  auto bx1 = this->belongs(KSpace::DIM, x1);
375  auto bx2 = this->belongs(KSpace::DIM, x2);
376  auto bx3 = this->belongs(KSpace::DIM, x3);
377  auto bx4 = this->belongs(KSpace::DIM, x4);
378  auto bx5 = this->belongs(KSpace::DIM, x5);
379  auto bx6 = this->belongs(KSpace::DIM, x6);
380  auto bx7 = this->belongs(KSpace::DIM, x7);
381 
382  auto by0 = this->belongs(KSpace::DIM, y0);
383  auto by1 = this->belongs(KSpace::DIM, y1);
384  auto by2 = this->belongs(KSpace::DIM, y2);
385  auto by3 = this->belongs(KSpace::DIM, y3);
386  auto by4 = this->belongs(KSpace::DIM, y4);
387  auto by5 = this->belongs(KSpace::DIM, y5);
388  auto by6 = this->belongs(KSpace::DIM, y6);
389  auto by7 = this->belongs(KSpace::DIM, y7);
390 
391  Clique k2_crit(ks);
392  if (bx0)
393  k2_crit.insertCell(x0);
394  if (bx1)
395  k2_crit.insertCell(x1);
396  if (bx2)
397  k2_crit.insertCell(x2);
398  if (bx3)
399  k2_crit.insertCell(x3);
400  if (bx4)
401  k2_crit.insertCell(x4);
402  if (bx5)
403  k2_crit.insertCell(x5);
404  if (bx6)
405  k2_crit.insertCell(x6);
406  if (bx7)
407  k2_crit.insertCell(x7);
408 
409  if (by0)
410  k2_crit.insertCell(y0);
411  if (by1)
412  k2_crit.insertCell(y1);
413  if (by2)
414  k2_crit.insertCell(y2);
415  if (by3)
416  k2_crit.insertCell(y3);
417  if (by4)
418  k2_crit.insertCell(y4);
419  if (by5)
420  k2_crit.insertCell(y5);
421  if (by6)
422  k2_crit.insertCell(y6);
423  if (by7)
424  k2_crit.insertCell(y7);
425  // Note that input spels A,B are ommited.
426 
427  /////////////////////////////////
428  // Critical Clique Conditions:
429  using namespace DGtal::functions;
430  // Intersection of k2-neighborhood with the object:
431  // (i) k2_clique must be empty or NOT 0-connected
432  bool is_empty{k2_crit.nbCells(KSpace::DIM) == 0};
433 
434  // Check connectedness using object if not empty
435  bool is_disconnected{false};
436  if (!is_empty) {
437  using DigitalTopology = DGtal::Z3i::DT26_6;
438  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
439  DGtal::Z3i::Domain,
440  std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
441  using NewObject = DGtal::Object<DigitalTopology, DigitalSet>;
442  auto new_obj = objectFromSpels<NewObject, KSpace, CellContainer>(k2_crit);
443  auto con = new_obj->computeConnectedness();
444  is_disconnected = (con == DISCONNECTED);
445  }
446 
447  bool conditionI = is_empty || is_disconnected;
448 
449  // (ii) Xi or Yi belongs to this for i={0,2,4,6}
450  std::vector<bool> bb(4);
451  bb[0] = bx0 || by0;
452  bb[1] = bx2 || by2;
453  bb[2] = bx4 || by4;
454  bb[3] = bx6 || by6;
455 
456  bool conditionII = bb[0] && bb[1] && bb[2] && bb[3];
457  // is_critical if any condition is true.
458  bool is_critical = conditionI || conditionII;
459 
460  if (verbose) {
461  trace.beginBlock(" K2 critical conditions ");
462  trace.info() << " conditionI = " << conditionI
463  << " : is_empty || is_disconnected = " << is_empty
464  << " && " << is_disconnected << std::endl;
465  trace.info() << " conditionII = " << conditionII << " : " << bb[0]
466  << " && " << bb[1] << " && " << bb[2] << " && " << bb[3]
467  << std::endl;
468  trace.info() << " is_critical = " << is_critical
469  << " conditionI || conditionII : " << conditionI << " || "
470  << conditionII << std::endl;
471  trace.endBlock();
472  }
473 
474  // Return the clique (A,B), not the mask k2_crit
475  Clique k2_clique(ks);
476  Cell Ac = ks.uSpel(A);
477  Cell Bc = ks.uSpel(B);
478  k2_clique.insertCell(Ac);
479  k2_clique.insertCell(Bc);
480  return std::make_pair(is_critical, k2_clique);
481 }
482 
483 //---------------------------------------------------------------------------
484 template <typename TKSpace, typename TCellContainer>
485 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
486  TCellContainer>::Clique>
487 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(const Cell &A,
488  const Cell &B,
489  bool verbose) const {
490  // Precondition:
491  // A and B are contiguous spels.
492  ASSERT(isSpel(A) == true);
493  ASSERT(isSpel(B) == true);
494  auto &ks = this->space();
495  auto B_coord = ks.uCoords(B);
496  auto A_coord = ks.uCoords(A);
497  return this->K_2(A_coord, B_coord, verbose);
498 }
499 //---------------------------------------------------------------------------
500 
501 template <typename TKSpace, typename TCellContainer>
502 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
503  TCellContainer>::Clique>
504 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_2(const Cell &face2,
505  bool verbose) const {
506  auto &ks = this->space();
507  ASSERT(ks.uIsSurfel(face2));
508  auto co_faces = ks.uCoFaces(face2);
509  ASSERT(co_faces.size() == 2);
510  auto &cf0 = co_faces[0];
511  auto &cf1 = co_faces[1];
512  // spels must belong to complex.
513  if (this->belongs(cf0) && this->belongs(cf1))
514  return this->K_2(cf0, cf1, verbose);
515  else
516  return std::make_pair(false, Clique(ks));
517 }
518 //---------------------------------------------------------------------------
519 
520 template <typename TKSpace, typename TCellContainer>
521 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
522  TCellContainer>::Clique>
523 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_1(const Cell &face1,
524  bool verbose) const {
525  auto &ks = this->space();
526  ASSERT(ks.uDim(face1) == 1);
527  // Get 2 orth dirs in orient_orth
528  std::vector<Point> dirs_orth;
529  for (auto q = ks.uOrthDirs(face1); q != 0; ++q) {
530  Dimension dir = *q;
531  Point positive_orth{0, 0, 0};
532  for (Dimension i = 0; i != ks.DIM; ++i)
533  if (i == dir)
534  positive_orth[i] = 1;
535 
536  dirs_orth.push_back(positive_orth);
537  }
538 
539  auto &kface = face1.preCell().coordinates;
540  Point a{kface + dirs_orth[0] + dirs_orth[1]};
541  Point b{kface + dirs_orth[0] - dirs_orth[1]};
542  Point c{kface - dirs_orth[0] + dirs_orth[1]};
543  Point d{kface - dirs_orth[0] - dirs_orth[1]};
544 
545  Cell A = ks.uCell(a);
546  Cell B = ks.uCell(b);
547  Cell C = ks.uCell(c);
548  Cell D = ks.uCell(d);
549 
550  // Now we need the other spels forming the mask
551  // Get the direction (positive) linel spans.
552  Point dir_parallel{0, 0, 0};
553  for (auto q = ks.uDirs(face1); q != 0; ++q) {
554  Dimension dir = *q;
555  for (Dimension i = 0; i != ks.DIM; ++i)
556  if (i == dir)
557  dir_parallel[i] = 1;
558  }
559  // Note that C, B are interchangeable. Same in A,D. Same between X and Y
560  // sets Changed notation from paper: XA=X0, XB=X1, XC=X2, XD=X3
561  // X
562  Point xa{a + 2 * dir_parallel};
563  Point xb{b + 2 * dir_parallel};
564  Point xc{c + 2 * dir_parallel};
565  Point xd{d + 2 * dir_parallel};
566  // Y
567  Point ya{a - 2 * dir_parallel};
568  Point yb{b - 2 * dir_parallel};
569  Point yc{c - 2 * dir_parallel};
570  Point yd{d - 2 * dir_parallel};
571 
572  // Cell of the mask from KCoords
573  Cell XA = ks.uCell(xa);
574  Cell XB = ks.uCell(xb);
575  Cell XC = ks.uCell(xc);
576  Cell XD = ks.uCell(xd);
577 
578  Cell YA = ks.uCell(ya);
579  Cell YB = ks.uCell(yb);
580  Cell YC = ks.uCell(yc);
581  Cell YD = ks.uCell(yd);
582 
583  /////////////////////////////////
584  // Critical Clique Conditions:
585 
586  /** is_critical = ConditionI && ConditionII
587  * (i) ConditionI:
588  * At least one the sets {A,D},{B,C} is subset of this complex
589  */
590  /** (ii) ConditionII = B1 OR B2
591  * B1) (U & *this != empty) AND (V & *this != empty)
592  * B2) (U & *this == empty) AND (V & *this == empty)
593  */
594 
595  const bool A1{this->belongs(KSpace::DIM, A) &&
596  this->belongs(KSpace::DIM, D)};
597  const bool A2{this->belongs(KSpace::DIM, B) &&
598  this->belongs(KSpace::DIM, C)};
599  const bool conditionI{A1 || A2};
600 
601  const bool u_not_empty{
602  this->belongs(KSpace::DIM, XA) || this->belongs(KSpace::DIM, XB) ||
603  this->belongs(KSpace::DIM, XC) || this->belongs(KSpace::DIM, XD)};
604 
605  const bool v_not_empty{
606  this->belongs(KSpace::DIM, YA) || this->belongs(KSpace::DIM, YB) ||
607  this->belongs(KSpace::DIM, YC) || this->belongs(KSpace::DIM, YD)};
608 
609  const bool B1{u_not_empty && v_not_empty};
610  const bool B2{!u_not_empty && !v_not_empty};
611  const bool conditionII{B1 || B2};
612 
613  const bool is_critical{conditionI && conditionII};
614 
615  if (verbose) {
616  trace.beginBlock(" K1 critical conditions ");
617  trace.info() << "input linel: " << face1 << std::endl;
618  trace.info() << "is_critical = " << is_critical
619  << " conditionI || condition II " << conditionI << " || "
620  << conditionII << std::endl;
621  trace.info() << " conditionI = " << conditionI << " = A1 || A2 : " << A1
622  << " || " << A2 << std::endl;
623  trace.info() << " conditionII = " << conditionII
624  << " = B1 || B2 : " << B1 << " || " << B2 << std::endl;
625  trace.endBlock();
626  }
627 
628  // out clique is the intersection between mask and object
629  Clique k1(ks);
630  if (this->belongs(KSpace::DIM, A))
631  k1.insert(A);
632  if (this->belongs(KSpace::DIM, B))
633  k1.insert(B);
634  if (this->belongs(KSpace::DIM, C))
635  k1.insert(C);
636  if (this->belongs(KSpace::DIM, D))
637  k1.insert(D);
638  return std::make_pair(is_critical, k1);
639 }
640 //---------------------------------------------------------------------------
641 
642 template <typename TKSpace, typename TCellContainer>
643 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
644  TCellContainer>::Clique>
645 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_0(const Cell &face0,
646  bool verbose) const {
647  auto &ks = this->space();
648  ASSERT(ks.uDim(face0) == 0);
649  auto &kface = face0.preCell().coordinates;
650  Point z{0, 0, 1};
651  Point y{0, 1, 0};
652  Point x{1, 0, 0};
653 
654  Point a{kface + x - y - z};
655  Point b{kface - x - y - z};
656  Point c{kface + x - y + z};
657  Point d{kface - x - y + z};
658 
659  Point e{kface + x + y - z};
660  Point f{kface - x + y - z};
661  Point g{kface + x + y + z};
662  Point h{kface - x + y + z};
663 
664  Cell A{ks.uCell(a)};
665  Cell B{ks.uCell(b)};
666  Cell C{ks.uCell(c)};
667  Cell D{ks.uCell(d)};
668 
669  Cell E{ks.uCell(e)};
670  Cell F{ks.uCell(f)};
671  Cell G{ks.uCell(g)};
672  Cell H{ks.uCell(h)};
673 
674  /////////////////////////////////
675  // Critical Clique Conditions:
676  /** is_critical = B1 || B2 || B3 || B4
677  * where:
678  * B1 = isSubset{A, H}
679  * B2 = isSubset{B, G}
680  * B3 = isSubset{C, F}
681  * B4 = isSubset{D, E}
682  * @note that the subsets define the 4 longest diagonals between the 8
683  * pixels.
684  */
685  const bool bA = this->belongs(KSpace::DIM, A);
686  const bool bB = this->belongs(KSpace::DIM, B);
687  const bool bC = this->belongs(KSpace::DIM, C);
688  const bool bD = this->belongs(KSpace::DIM, D);
689  const bool bE = this->belongs(KSpace::DIM, E);
690  const bool bF = this->belongs(KSpace::DIM, F);
691  const bool bG = this->belongs(KSpace::DIM, G);
692  const bool bH = this->belongs(KSpace::DIM, H);
693 
694  const bool B1{bA && bH};
695  const bool B2{bB && bG};
696  const bool B3{bC && bF};
697  const bool B4{bD && bE};
698  const bool is_critical{B1 || B2 || B3 || B4};
699 
700  if (verbose) {
701  trace.beginBlock(" K0 critical conditions ");
702  trace.info() << "input pointel: " << face0 << std::endl;
703  trace.info() << "is_critical = B1 || B2 || B3 || B4 " << std::endl;
704  trace.info() << is_critical << " = " << B1 << " || " << B2 << " || "
705  << B3 << " || " << B4 << std::endl;
706  trace.endBlock();
707  }
708  // out clique is the intersection between mask and object
709  Clique k0_out(ks);
710  if (bA)
711  k0_out.insert(A);
712  if (bB)
713  k0_out.insert(B);
714  if (bC)
715  k0_out.insert(C);
716  if (bD)
717  k0_out.insert(D);
718  if (bE)
719  k0_out.insert(E);
720  if (bF)
721  k0_out.insert(F);
722  if (bG)
723  k0_out.insert(G);
724  if (bH)
725  k0_out.insert(H);
726  using namespace DGtal::functions;
727 
728  return std::make_pair(is_critical, k0_out);
729 }
730 //---------------------------------------------------------------------------
731 
732 template <typename TKSpace, typename TCellContainer>
733 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
734  TCellContainer>::Clique>
735 DGtal::VoxelComplex<TKSpace, TCellContainer>::K_3(const Cell &voxel,
736  bool verbose) const {
737  auto &ks = this->space();
738  ASSERT(ks.uDim(voxel) == 3);
739  const bool is_critical = !isSimple(voxel);
740 
741  if (verbose) {
742  trace.beginBlock(" K3 critical conditions ");
743  trace.info() << "input voxel: " << voxel << std::endl;
744  trace.info() << "is_critical = " << is_critical << std::endl;
745  trace.endBlock();
746  }
747 
748  Clique clique(ks);
749  clique.insertCell(voxel);
750  return std::make_pair(is_critical, clique);
751 }
752 //---------------------------------------------------------------------------
753 
754 /* BUG workaround: MSVC compiler error C2244.
755  * It doesn't see the definition of these declarations (Moved to header)
756 template <typename TKSpace, typename TCellContainer>
757 std::array<
758  typename DGtal::VoxelComplex<TKSpace,
759 TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
760 TCellContainer>::dimension + 1
761 >
762 DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
763  const Parent & cubical,
764  bool verbose
765  ) const
766 
767 template <typename TKSpace, typename TCellContainer>
768 std::array<
769  typename DGtal::VoxelComplex<TKSpace,
770 TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
771 TCellContainer>::dimension + 1
772 >
773 DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
774  bool verbose
775  ) const
776 */
777 //---------------------------------------------------------------------------
778 
779 template <typename TKSpace, typename TCellContainer>
780 std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
781  TCellContainer>::Clique>
782 DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliquePair(
783  const Dimension d, const CellMapConstIterator &cellMapIterator) const {
784  auto &it = cellMapIterator;
785  auto &cell = it->first;
786  // auto &cell_data = it->second;
787  auto clique_p = std::make_pair(false, Clique(this->space()));
788  if (d == 0)
789  clique_p = K_0(cell);
790  else if (d == 1)
791  clique_p = K_1(cell);
792  else if (d == 2)
793  clique_p = K_2(cell);
794  else if (d == 3)
795  clique_p = K_3(cell);
796  else
797  throw std::runtime_error("Wrong dimension: " + std::to_string(d));
798 
799  return clique_p;
800 }
801 //---------------------------------------------------------------------------
802 
803 template <typename TKSpace, typename TCellContainer>
804 typename DGtal::VoxelComplex<TKSpace, TCellContainer>::CliqueContainer
805 DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliquesForD(
806  const Dimension d, const Parent &cubical, bool verbose) const {
807 #ifdef WITH_OPENMP
808 
809  ASSERT(dimension >= 0 && dimension <= 3);
810  CliqueContainer critical;
811 
812  auto nthreads = omp_get_num_procs();
813  omp_set_num_threads(nthreads);
814  std::vector<CliqueContainer> p_critical;
815  p_critical.resize(nthreads);
816 #pragma omp parallel
817  {
818 #pragma omp single nowait
819  {for (auto it = cubical.begin(d), itE = cubical.end(d); it != itE; ++it)
820 #pragma omp task firstprivate(it)
821  {auto clique_p = criticalCliquePair(d, it);
822  auto &is_critical = clique_p.first;
823  auto &clique = clique_p.second;
824  // Push
825  auto th = omp_get_thread_num();
826  if (is_critical)
827  p_critical[th].push_back(clique);
828 }
829 } // cell loop
830 #pragma omp taskwait
831 }
832 // Merge
833 std::size_t total_size = 0;
834 for (const auto &sub : p_critical)
835  total_size += sub.size();
836 
837 critical.reserve(total_size);
838 for (const auto &sub : p_critical)
839  critical.insert(critical.end(), sub.begin(), sub.end());
840 
841 if (verbose)
842  trace.info() << " d:" << d << " ncrit: " << critical.size();
843 return critical;
844 
845 #else
846 
847  ASSERT(dimension >= 0 && dimension <= 3);
848  CliqueContainer critical;
849  for (auto it = cubical.begin(d), itE = cubical.end(d); it != itE; ++it) {
850  auto clique_p = criticalCliquePair(d, it);
851  auto &is_critical = clique_p.first;
852  auto &clique = clique_p.second;
853  if (is_critical)
854  critical.push_back(clique);
855  } // cell loop
856  if (verbose)
857  trace.info() << " d:" << d << " ncrit: " << critical.size();
858  return critical;
859 
860 #endif
861 }
862 //---------------------------------------------------------------------------
863 ///////////////////////////////////////////////////////////////////////////////
864 template <typename TKSpace, typename TCellContainer>
865 bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSimpleByThinning(
866  const Cell &input_spel) const {
867  // x = input_spel ; X = VoxelComplex ~ occupancy of khalimsky space
868  // a) Get the neighborhood (voxels) of input_spel intersected
869  // with the voxel complex. -- N^{*}(x) intersection X --
870  ASSERT(this->space().uDim(input_spel) == 3);
871  auto spels_out = this->properNeighborhoodVoxels(input_spel);
872  auto &ks = this->space();
873  Clique clique(ks);
874  for (const auto &v : spels_out)
875  clique.insertCell(v);
876  clique.close();
877  // b) Apply a thinning on the result of a)
878  typename Parent::DefaultCellMapIteratorPriority default_priority;
879  bool clique_is_closed = true;
880  functions::collapse( clique, spels_out.begin(), spels_out.end(), default_priority, false /* spels_out is not closed */, clique_is_closed, false /*verbose*/);
881  // c) If the result is a single pointel, it is reducible
882  return clique.size() == 1;
883 }
884 
885 // Object wrappers :
886 template <typename TKSpace, typename TCellContainer>
887 bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSimple(
888  const Cell &input_cell) const {
889  ASSERT(isSpel(input_cell) == true);
890 
891  if (myIsTableLoaded) {
892  auto conf = functions::getSpelNeighborhoodConfigurationOccupancy<Self>(
893  *this, this->space().uCoords(input_cell), this->pointToMask());
894  return (*myTablePtr)[conf];
895  } else
896  return isSimpleByThinning(input_cell);
897 }
898 //---------------------------------------------------------------------------
899 ///////////////////////////////////////////////////////////////////////////////
900 // Interface - public :
901 
902 /**
903  * Writes/Displays the object on an output stream.
904  * @param out the output stream where the object is written.
905  */
906 template <typename TKSpace, typename TCellContainer>
907 inline void DGtal::VoxelComplex<TKSpace, TCellContainer>::selfDisplay(
908  std::ostream &out) const {
909  out << "[VoxelComplex dim=" << this->dim() << " chi=" << this->euler();
910  out << " isTableLoaded? " << ((isTableLoaded()) ? "True" : "False");
911 }
912 //---------------------------------------------------------------------------
913 
914 /**
915  * Checks the validity/consistency of the object.
916  * @return 'true' if the object is valid, 'false' otherwise.
917  */
918 template <typename TKSpace, typename TCellContainer>
919 inline bool
920 DGtal::VoxelComplex<TKSpace, TCellContainer>::isValid() const {
921  return true;
922 }
923 
924 //-----------------------------------------------------------------------------
925 template <typename TKSpace, typename TCellContainer>
926 inline std::string
927 DGtal::VoxelComplex<TKSpace, TCellContainer>::className() const {
928  return "VoxelComplex";
929 }
930 
931 ///////////////////////////////////////////////////////////////////////////////
932 // Implementation of inline functions //
933 
934 template <typename TKSpace, typename TCellContainer>
935 inline std::ostream &DGtal::
936 operator<<(std::ostream &out,
937  const VoxelComplex<TKSpace, TCellContainer> &object) {
938  object.selfDisplay(out);
939  return out;
940 }
941 
942 // //
943 ///////////////////////////////////////////////////////////////////////////////