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