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