DGtal 1.4.2
Loading...
Searching...
No Matches
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:
43template <typename TKSpace, typename TCellContainer>
44inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex()
45 : Parent(),
46 myTablePtr(nullptr), myPointToMaskPtr(nullptr),
47 myIsTableLoaded(false) {}
48
49// Copy constructor:
50template <typename TKSpace, typename TCellContainer>
51inline 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//-----------------------------------------------------------------------------
62template <typename TKSpace, typename TCellContainer>
63inline DGtal::VoxelComplex<TKSpace, TCellContainer> &
64DGtal::VoxelComplex<TKSpace, TCellContainer>::
65operator=(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
98template <typename TKSpace, typename TCellContainer>
99template <typename TDigitalSet>
100inline 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
108template <typename TKSpace, typename TCellContainer>
109void 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
118template <typename TKSpace, typename TCellContainer>
119void DGtal::VoxelComplex<TKSpace, TCellContainer>::copySimplicityTable(
120 const Self & other)
121{
122 myTablePtr = other.myTablePtr;
123 myPointToMaskPtr = other.myPointToMaskPtr;
124 myIsTableLoaded = other.myIsTableLoaded;
125}
126
127template <typename TKSpace, typename TCellContainer>
128const typename DGtal::VoxelComplex<TKSpace,
129 TCellContainer>::ConfigMap &
130DGtal::VoxelComplex<TKSpace, TCellContainer>::table() const
131{
132 return *myTablePtr;
133}
134
135template <typename TKSpace, typename TCellContainer>
136const bool &
137DGtal::VoxelComplex<TKSpace, TCellContainer>::isTableLoaded() const
138{
139 return myIsTableLoaded;
140}
141
142template <typename TKSpace, typename TCellContainer>
143const typename DGtal::VoxelComplex<TKSpace,
144 TCellContainer>::PointToMaskMap &
145DGtal::VoxelComplex<TKSpace, TCellContainer>::pointToMask() const
146{
147 return *myPointToMaskPtr;
148}
149//---------------------------------------------------------------------------
150template <typename TKSpace, typename TCellContainer>
151inline 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//---------------------------------------------------------------------------
166template <typename TKSpace, typename TCellContainer>
167inline 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//---------------------------------------------------------------------------
187template <typename TKSpace, typename TCellContainer>
188inline void
189DGtal::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//---------------------------------------------------------------------------
199template <typename TKSpace, typename TCellContainer>
200inline void
201DGtal::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//---------------------------------------------------------------------------
209template <typename TKSpace, typename TCellContainer>
210template <typename TDigitalSet>
211inline void
212DGtal::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
222template <typename TKSpace, typename TCellContainer>
223void 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//---------------------------------------------------------------------------
238template <typename TKSpace, typename TCellContainer>
239void 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//---------------------------------------------------------------------------
257template <typename TKSpace, typename TCellContainer>
258typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique
259DGtal::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
272template <typename TKSpace, typename TCellContainer>
273std::set<typename TKSpace::Cell>
274DGtal::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
285template <typename TKSpace, typename TCellContainer>
286std::set<typename TKSpace::Cell>
287DGtal::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//---------------------------------------------------------------------------
300template <typename TKSpace, typename TCellContainer>
301bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSpel(
302 const Cell &b) const
303{
304 return (this->space().uDim(b) == this->space().DIM);
305}
306
307//---------------------------------------------------------------------------
308template <typename TKSpace, typename TCellContainer>
309typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Cell
310DGtal::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
326template <typename TKSpace, typename TCellContainer>
327std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
328 TCellContainer>::Clique>
329DGtal::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//---------------------------------------------------------------------------
488template <typename TKSpace, typename TCellContainer>
489std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
490DGtal::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
505template <typename TKSpace, typename TCellContainer>
506std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
507DGtal::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
525template <typename TKSpace, typename TCellContainer>
526std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
527DGtal::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
648template <typename TKSpace, typename TCellContainer>
649std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
650DGtal::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
738template <typename TKSpace, typename TCellContainer>
739std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
740DGtal::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)
762template <typename TKSpace, typename TCellContainer>
763std::array<
764 typename DGtal::VoxelComplex<TKSpace,
765TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
766TCellContainer>::dimension + 1
767>
768DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
769 const Parent & cubical,
770 bool verbose
771 ) const
772
773template <typename TKSpace, typename TCellContainer>
774std::array<
775 typename DGtal::VoxelComplex<TKSpace,
776TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
777TCellContainer>::dimension + 1
778>
779DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
780 bool verbose
781 ) const
782*/
783//---------------------------------------------------------------------------
784
785template <typename TKSpace, typename TCellContainer>
786std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
787DGtal::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
809template <typename TKSpace, typename TCellContainer>
810typename DGtal::VoxelComplex<TKSpace, TCellContainer>::CliqueContainer
811DGtal::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
839std::size_t total_size = 0;
840for (const auto &sub : p_critical)
841 total_size += sub.size();
842
843critical.reserve(total_size);
844for (const auto &sub : p_critical)
845 critical.insert(critical.end(), sub.begin(), sub.end());
846
847if (verbose)
848 trace.info() << " d:" << d << " ncrit: " << critical.size();
849return 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///////////////////////////////////////////////////////////////////////////////
870template <typename TKSpace, typename TCellContainer>
871bool 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 :
893template <typename TKSpace, typename TCellContainer>
894bool 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 */
914template <typename TKSpace, typename TCellContainer>
915inline 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 */
927template <typename TKSpace, typename TCellContainer>
928inline bool
929DGtal::VoxelComplex<TKSpace, TCellContainer>::isValid() const
930{
931 return true;
932}
933
934//-----------------------------------------------------------------------------
935template <typename TKSpace, typename TCellContainer>
936inline std::string
937DGtal::VoxelComplex<TKSpace, TCellContainer>::className() const
938{
939 return "VoxelComplex";
940}
941
942///////////////////////////////////////////////////////////////////////////////
943// Implementation of inline functions //
944
945template <typename TKSpace, typename TCellContainer>
946inline std::ostream &DGtal::
947operator<<(std::ostream &out,
948 const VoxelComplex<TKSpace, TCellContainer> &object)
949{
950 object.selfDisplay(out);
951 return out;
952}
953
954// //
955///////////////////////////////////////////////////////////////////////////////