DGtal 1.4.0
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 <experimental/algorithm>
40#include <omp.h>
41#endif
42//////////////////////////////////////////////////////////////////////////////
43// Default constructor:
44template <typename TKSpace, typename TCellContainer>
45inline DGtal::VoxelComplex<TKSpace, TCellContainer>::VoxelComplex()
46 : Parent(),
47 myTablePtr(nullptr), myPointToMaskPtr(nullptr),
48 myIsTableLoaded(false) {}
49
50// Copy constructor:
51template <typename TKSpace, typename TCellContainer>
52inline 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//-----------------------------------------------------------------------------
63template <typename TKSpace, typename TCellContainer>
64inline DGtal::VoxelComplex<TKSpace, TCellContainer> &
65DGtal::VoxelComplex<TKSpace, TCellContainer>::
66operator=(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
99template <typename TKSpace, typename TCellContainer>
100template <typename TDigitalSet>
101inline 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
109template <typename TKSpace, typename TCellContainer>
110void 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
119template <typename TKSpace, typename TCellContainer>
120void DGtal::VoxelComplex<TKSpace, TCellContainer>::copySimplicityTable(
121 const Self & other)
122{
123 myTablePtr = other.myTablePtr;
124 myPointToMaskPtr = other.myPointToMaskPtr;
125 myIsTableLoaded = other.myIsTableLoaded;
126}
127
128template <typename TKSpace, typename TCellContainer>
129const typename DGtal::VoxelComplex<TKSpace,
130 TCellContainer>::ConfigMap &
131DGtal::VoxelComplex<TKSpace, TCellContainer>::table() const
132{
133 return *myTablePtr;
134}
135
136template <typename TKSpace, typename TCellContainer>
137const bool &
138DGtal::VoxelComplex<TKSpace, TCellContainer>::isTableLoaded() const
139{
140 return myIsTableLoaded;
141}
142
143template <typename TKSpace, typename TCellContainer>
144const typename DGtal::VoxelComplex<TKSpace,
145 TCellContainer>::PointToMaskMap &
146DGtal::VoxelComplex<TKSpace, TCellContainer>::pointToMask() const
147{
148 return *myPointToMaskPtr;
149}
150//---------------------------------------------------------------------------
151template <typename TKSpace, typename TCellContainer>
152inline 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//---------------------------------------------------------------------------
167template <typename TKSpace, typename TCellContainer>
168inline 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//---------------------------------------------------------------------------
188template <typename TKSpace, typename TCellContainer>
189inline void
190DGtal::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//---------------------------------------------------------------------------
200template <typename TKSpace, typename TCellContainer>
201inline void
202DGtal::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//---------------------------------------------------------------------------
210template <typename TKSpace, typename TCellContainer>
211template <typename TDigitalSet>
212inline void
213DGtal::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
223template <typename TKSpace, typename TCellContainer>
224void 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//---------------------------------------------------------------------------
239template <typename TKSpace, typename TCellContainer>
240void 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//---------------------------------------------------------------------------
258template <typename TKSpace, typename TCellContainer>
259typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique
260DGtal::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
273template <typename TKSpace, typename TCellContainer>
274std::set<typename TKSpace::Cell>
275DGtal::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
286template <typename TKSpace, typename TCellContainer>
287std::set<typename TKSpace::Cell>
288DGtal::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//---------------------------------------------------------------------------
301template <typename TKSpace, typename TCellContainer>
302bool DGtal::VoxelComplex<TKSpace, TCellContainer>::isSpel(
303 const Cell &b) const
304{
305 return (this->space().uDim(b) == this->space().DIM);
306}
307
308//---------------------------------------------------------------------------
309template <typename TKSpace, typename TCellContainer>
310typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Cell
311DGtal::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
327template <typename TKSpace, typename TCellContainer>
328std::pair<bool, typename DGtal::VoxelComplex<TKSpace,
329 TCellContainer>::Clique>
330DGtal::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//---------------------------------------------------------------------------
489template <typename TKSpace, typename TCellContainer>
490std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
491DGtal::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
506template <typename TKSpace, typename TCellContainer>
507std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
508DGtal::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
526template <typename TKSpace, typename TCellContainer>
527std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
528DGtal::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
649template <typename TKSpace, typename TCellContainer>
650std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
651DGtal::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
739template <typename TKSpace, typename TCellContainer>
740std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
741DGtal::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)
763template <typename TKSpace, typename TCellContainer>
764std::array<
765 typename DGtal::VoxelComplex<TKSpace,
766TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
767TCellContainer>::dimension + 1
768>
769DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
770 const Parent & cubical,
771 bool verbose
772 ) const
773
774template <typename TKSpace, typename TCellContainer>
775std::array<
776 typename DGtal::VoxelComplex<TKSpace,
777TCellContainer>::CliqueContainer, DGtal::VoxelComplex<TKSpace,
778TCellContainer>::dimension + 1
779>
780DGtal::VoxelComplex<TKSpace, TCellContainer>::criticalCliques(
781 bool verbose
782 ) const
783*/
784//---------------------------------------------------------------------------
785
786template <typename TKSpace, typename TCellContainer>
787std::pair<bool, typename DGtal::VoxelComplex<TKSpace, TCellContainer>::Clique>
788DGtal::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
810template <typename TKSpace, typename TCellContainer>
811typename DGtal::VoxelComplex<TKSpace, TCellContainer>::CliqueContainer
812DGtal::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
840std::size_t total_size = 0;
841for (const auto &sub : p_critical)
842 total_size += sub.size();
843
844critical.reserve(total_size);
845for (const auto &sub : p_critical)
846 critical.insert(critical.end(), sub.begin(), sub.end());
847
848if (verbose)
849 trace.info() << " d:" << d << " ncrit: " << critical.size();
850return 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///////////////////////////////////////////////////////////////////////////////
871template <typename TKSpace, typename TCellContainer>
872bool 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 :
894template <typename TKSpace, typename TCellContainer>
895bool 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 */
915template <typename TKSpace, typename TCellContainer>
916inline 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 */
928template <typename TKSpace, typename TCellContainer>
929inline bool
930DGtal::VoxelComplex<TKSpace, TCellContainer>::isValid() const
931{
932 return true;
933}
934
935//-----------------------------------------------------------------------------
936template <typename TKSpace, typename TCellContainer>
937inline std::string
938DGtal::VoxelComplex<TKSpace, TCellContainer>::className() const
939{
940 return "VoxelComplex";
941}
942
943///////////////////////////////////////////////////////////////////////////////
944// Implementation of inline functions //
945
946template <typename TKSpace, typename TCellContainer>
947inline std::ostream &DGtal::
948operator<<(std::ostream &out,
949 const VoxelComplex<TKSpace, TCellContainer> &object)
950{
951 object.selfDisplay(out);
952 return out;
953}
954
955// //
956///////////////////////////////////////////////////////////////////////////////