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