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