1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @file tools/LevelSetUtil.h
5 ///
6 /// @brief Miscellaneous utility methods that operate primarily
7 /// or exclusively on level set grids.
8 ///
9 /// @author Mihai Alden
10
11 #ifndef OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
12 #define OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
13
14 #include "MeshToVolume.h" // for traceExteriorBoundaries
15 #include "SignedFloodFill.h" // for signedFloodFillWithValues
16
17 #include <openvdb/Types.h>
18 #include <openvdb/Grid.h>
19 #include <openvdb/openvdb.h>
20 #include <openvdb/points/PointDataGrid.h>
21 #include <tbb/blocked_range.h>
22 #include <tbb/parallel_for.h>
23 #include <tbb/parallel_reduce.h>
24 #include <tbb/parallel_sort.h>
25 #include <algorithm>
26 #include <cmath>
27 #include <cstdlib>
28 #include <deque>
29 #include <limits>
30 #include <memory>
31 #include <set>
32 #include <vector>
33
34
35 namespace openvdb {
36 OPENVDB_USE_VERSION_NAMESPACE
37 namespace OPENVDB_VERSION_NAME {
38 namespace tools {
39
40 // MS Visual C++ requires this extra level of indirection in order to compile
41 // THIS MUST EXIST IN AN UNNAMED NAMESPACE IN ORDER TO COMPILE ON WINDOWS
42 namespace {
43
44 template<typename GridType>
lsutilGridMax()45 inline typename GridType::ValueType lsutilGridMax()
46 {
47 return std::numeric_limits<typename GridType::ValueType>::max();
48 }
49
50 template<typename GridType>
lsutilGridZero()51 inline typename GridType::ValueType lsutilGridZero()
52 {
53 return zeroVal<typename GridType::ValueType>();
54 }
55
56 } // unnamed namespace
57
58
59 ////////////////////////////////////////
60
61
62 /// @brief Threaded method to convert a sparse level set/SDF into a sparse fog volume
63 ///
64 /// @details For a level set, the active and negative-valued interior half of the
65 /// narrow band becomes a linear ramp from 0 to 1; the inactive interior becomes
66 /// active with a constant value of 1; and the exterior, including the background
67 /// and the active exterior half of the narrow band, becomes inactive with a constant
68 /// value of 0. The interior, though active, remains sparse.
69 /// @details For a generic SDF, a specified cutoff distance determines the width
70 /// of the ramp, but otherwise the result is the same as for a level set.
71 ///
72 /// @param grid level set/SDF grid to transform
73 /// @param cutoffDistance optional world space cutoff distance for the ramp
74 /// (automatically clamped if greater than the interior
75 /// narrow band width)
76 template<class GridType>
77 void
78 sdfToFogVolume(
79 GridType& grid,
80 typename GridType::ValueType cutoffDistance = lsutilGridMax<GridType>());
81
82
83 /// @brief Threaded method to construct a boolean mask that represents interior regions
84 /// in a signed distance field.
85 ///
86 /// @return A shared pointer to either a boolean grid or tree with the same tree
87 /// configuration and potentially transform as the input @c volume and whose active
88 /// and @c true values correspond to the interior of the input signed distance field.
89 ///
90 /// @param volume Signed distance field / level set volume.
91 /// @param isovalue Threshold below which values are considered part of the
92 /// interior region.
93 template<class GridOrTreeType>
94 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
95 sdfInteriorMask(
96 const GridOrTreeType& volume,
97 typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>());
98
99
100 /// @brief Extracts the interior regions of a signed distance field and topologically enclosed
101 /// (watertight) regions of value greater than the @a isovalue (cavities) that can arise
102 /// as the result of CSG union operations between different shapes where at least one of
103 /// the shapes has a concavity that is capped.
104 ///
105 /// For example the enclosed region of a capped bottle would include the walls and
106 /// the interior cavity.
107 ///
108 /// @return A shared pointer to either a boolean grid or tree with the same tree configuration
109 /// and potentially transform as the input @c volume and whose active and @c true values
110 /// correspond to the interior and enclosed regions in the input signed distance field.
111 ///
112 /// @param volume Signed distance field / level set volume.
113 /// @param isovalue Threshold below which values are considered part of the interior region.
114 /// @param fillMask Optional boolean tree, when provided enclosed cavity regions that are not
115 /// completely filled by this mask are ignored.
116 ///
117 /// For instance if the fill mask does not completely fill the bottle in the
118 /// previous example only the walls and cap are returned and the interior
119 /// cavity will be ignored.
120 template<typename GridOrTreeType>
121 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
122 extractEnclosedRegion(
123 const GridOrTreeType& volume,
124 typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>(),
125 const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
126 fillMask = nullptr);
127
128
129 /// @brief Return a mask of the voxels that intersect the implicit surface with
130 /// the given @a isovalue.
131 ///
132 /// @param volume Signed distance field / level set volume.
133 /// @param isovalue The crossing point that is considered the surface.
134 template<typename GridOrTreeType>
135 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
136 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue);
137
138
139 /// @brief Return a mask for each connected component of the given grid's active voxels.
140 ///
141 /// @param volume Input grid or tree
142 /// @param masks Output set of disjoint active topology masks sorted in descending order
143 /// based on the active voxel count.
144 template<typename GridOrTreeType>
145 void
146 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
147 std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks);
148
149
150 /// @brief Separates disjoint active topology components into distinct grids or trees.
151 ///
152 /// @details Supports volumes with active tiles.
153 ///
154 /// @param volume Input grid or tree
155 /// @param segments Output set of disjoint active topology components sorted in
156 /// descending order based on the active voxel count.
157 template<typename GridOrTreeType>
158 void
159 segmentActiveVoxels(const GridOrTreeType& volume,
160 std::vector<typename GridOrTreeType::Ptr>& segments);
161
162
163 /// @brief Separates disjoint SDF surfaces into distinct grids or trees.
164 ///
165 /// @details Supports asymmetric interior / exterior narrowband widths and
166 /// SDF volumes with dense interior regions.
167 ///
168 /// @param volume Input signed distance field / level set volume
169 /// @param segments Output set of disjoint SDF surfaces found in @a volume sorted in
170 /// descending order based on the surface intersecting voxel count.
171 template<typename GridOrTreeType>
172 void
173 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
174
175
176 ////////////////////////////////////////////////////////////////////////////////
177 ////////////////////////////////////////////////////////////////////////////////
178
179 // Internal utility objects and implementation details
180
181 /// @cond OPENVDB_DOCS_INTERNAL
182
183 namespace level_set_util_internal {
184
185
186 template<typename LeafNodeType>
187 struct MaskInteriorVoxels {
188
189 using ValueType = typename LeafNodeType::ValueType;
190 using BoolLeafNodeType = tree::LeafNode<bool, LeafNodeType::LOG2DIM>;
191
MaskInteriorVoxelsMaskInteriorVoxels192 MaskInteriorVoxels(
193 ValueType isovalue, const LeafNodeType ** nodes, BoolLeafNodeType ** maskNodes)
194 : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
195 {
196 }
197
operatorMaskInteriorVoxels198 void operator()(const tbb::blocked_range<size_t>& range) const {
199
200 BoolLeafNodeType * maskNodePt = nullptr;
201
202 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
203
204 mMaskNodes[n] = nullptr;
205 const LeafNodeType& node = *mNodes[n];
206
207 if (!maskNodePt) {
208 maskNodePt = new BoolLeafNodeType(node.origin(), false);
209 } else {
210 maskNodePt->setOrigin(node.origin());
211 }
212
213 const ValueType* values = &node.getValue(0);
214 for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
215 if (values[i] < mIsovalue) maskNodePt->setValueOn(i, true);
216 }
217
218 if (maskNodePt->onVoxelCount() > 0) {
219 mMaskNodes[n] = maskNodePt;
220 maskNodePt = nullptr;
221 }
222 }
223
224 if (maskNodePt) delete maskNodePt;
225 }
226
227 LeafNodeType const * const * const mNodes;
228 BoolLeafNodeType ** const mMaskNodes;
229 ValueType const mIsovalue;
230 }; // MaskInteriorVoxels
231
232
233 template<typename TreeType, typename InternalNodeType>
234 struct MaskInteriorTiles {
235
236 using ValueType = typename TreeType::ValueType;
237
MaskInteriorTilesMaskInteriorTiles238 MaskInteriorTiles(ValueType isovalue, const TreeType& tree, InternalNodeType ** maskNodes)
239 : mTree(&tree), mMaskNodes(maskNodes), mIsovalue(isovalue) { }
240
operatorMaskInteriorTiles241 void operator()(const tbb::blocked_range<size_t>& range) const {
242 tree::ValueAccessor<const TreeType> acc(*mTree);
243 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
244 typename InternalNodeType::ValueAllIter it = mMaskNodes[n]->beginValueAll();
245 for (; it; ++it) {
246 if (acc.getValue(it.getCoord()) < mIsovalue) {
247 it.setValue(true);
248 it.setValueOn(true);
249 }
250 }
251 }
252 }
253
254 TreeType const * const mTree;
255 InternalNodeType ** const mMaskNodes;
256 ValueType const mIsovalue;
257 }; // MaskInteriorTiles
258
259
260 template<typename TreeType>
261 struct PopulateTree {
262
263 using ValueType = typename TreeType::ValueType;
264 using LeafNodeType = typename TreeType::LeafNodeType;
265
PopulateTreePopulateTree266 PopulateTree(TreeType& tree, LeafNodeType** leafnodes,
267 const size_t * nodexIndexMap, ValueType background)
268 : mNewTree(background)
269 , mTreePt(&tree)
270 , mNodes(leafnodes)
271 , mNodeIndexMap(nodexIndexMap)
272 {
273 }
274
PopulateTreePopulateTree275 PopulateTree(PopulateTree& rhs, tbb::split)
276 : mNewTree(rhs.mNewTree.background())
277 , mTreePt(&mNewTree)
278 , mNodes(rhs.mNodes)
279 , mNodeIndexMap(rhs.mNodeIndexMap)
280 {
281 }
282
operatorPopulateTree283 void operator()(const tbb::blocked_range<size_t>& range) {
284
285 tree::ValueAccessor<TreeType> acc(*mTreePt);
286
287 if (mNodeIndexMap) {
288 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
289 for (size_t i = mNodeIndexMap[n], I = mNodeIndexMap[n + 1]; i < I; ++i) {
290 if (mNodes[i] != nullptr) acc.addLeaf(mNodes[i]);
291 }
292 }
293 } else {
294 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
295 acc.addLeaf(mNodes[n]);
296 }
297 }
298 }
299
joinPopulateTree300 void join(PopulateTree& rhs) { mTreePt->merge(*rhs.mTreePt); }
301
302 private:
303 TreeType mNewTree;
304 TreeType * const mTreePt;
305 LeafNodeType ** const mNodes;
306 size_t const * const mNodeIndexMap;
307 }; // PopulateTree
308
309
310 /// @brief Negative active values are set @c 0, everything else is set to @c 1.
311 template<typename LeafNodeType>
312 struct LabelBoundaryVoxels {
313
314 using ValueType = typename LeafNodeType::ValueType;
315 using CharLeafNodeType = tree::LeafNode<char, LeafNodeType::LOG2DIM>;
316
LabelBoundaryVoxelsLabelBoundaryVoxels317 LabelBoundaryVoxels(
318 ValueType isovalue, const LeafNodeType ** nodes, CharLeafNodeType ** maskNodes)
319 : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
320 {
321 }
322
operatorLabelBoundaryVoxels323 void operator()(const tbb::blocked_range<size_t>& range) const {
324
325 CharLeafNodeType * maskNodePt = nullptr;
326
327 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
328
329 mMaskNodes[n] = nullptr;
330 const LeafNodeType& node = *mNodes[n];
331
332 if (!maskNodePt) {
333 maskNodePt = new CharLeafNodeType(node.origin(), 1);
334 } else {
335 maskNodePt->setOrigin(node.origin());
336 }
337
338 typename LeafNodeType::ValueOnCIter it;
339 for (it = node.cbeginValueOn(); it; ++it) {
340 maskNodePt->setValueOn(it.pos(), ((*it - mIsovalue) < 0.0) ? 0 : 1);
341 }
342
343 if (maskNodePt->onVoxelCount() > 0) {
344 mMaskNodes[n] = maskNodePt;
345 maskNodePt = nullptr;
346 }
347 }
348
349 if (maskNodePt) delete maskNodePt;
350 }
351
352 LeafNodeType const * const * const mNodes;
353 CharLeafNodeType ** const mMaskNodes;
354 ValueType const mIsovalue;
355 }; // LabelBoundaryVoxels
356
357
358 template<typename LeafNodeType>
359 struct FlipRegionSign {
360 using ValueType = typename LeafNodeType::ValueType;
361
FlipRegionSignFlipRegionSign362 FlipRegionSign(LeafNodeType ** nodes) : mNodes(nodes) { }
363
operatorFlipRegionSign364 void operator()(const tbb::blocked_range<size_t>& range) const {
365 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
366 ValueType* values = const_cast<ValueType*>(&mNodes[n]->getValue(0));
367 for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
368 values[i] = values[i] < 0 ? 1 : -1;
369 }
370 }
371 }
372
373 LeafNodeType ** const mNodes;
374 }; // FlipRegionSign
375
376
377 template<typename LeafNodeType>
378 struct FindMinVoxelValue {
379
380 using ValueType = typename LeafNodeType::ValueType;
381
FindMinVoxelValueFindMinVoxelValue382 FindMinVoxelValue(LeafNodeType const * const * const leafnodes)
383 : minValue(std::numeric_limits<ValueType>::max())
384 , mNodes(leafnodes)
385 {
386 }
387
FindMinVoxelValueFindMinVoxelValue388 FindMinVoxelValue(FindMinVoxelValue& rhs, tbb::split)
389 : minValue(std::numeric_limits<ValueType>::max())
390 , mNodes(rhs.mNodes)
391 {
392 }
393
operatorFindMinVoxelValue394 void operator()(const tbb::blocked_range<size_t>& range) {
395 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
396 const ValueType* data = mNodes[n]->buffer().data();
397 for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
398 minValue = std::min(minValue, data[i]);
399 }
400 }
401 }
402
joinFindMinVoxelValue403 void join(FindMinVoxelValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
404
405 ValueType minValue;
406
407 LeafNodeType const * const * const mNodes;
408 }; // FindMinVoxelValue
409
410
411 template<typename InternalNodeType>
412 struct FindMinTileValue {
413
414 using ValueType = typename InternalNodeType::ValueType;
415
FindMinTileValueFindMinTileValue416 FindMinTileValue(InternalNodeType const * const * const nodes)
417 : minValue(std::numeric_limits<ValueType>::max())
418 , mNodes(nodes)
419 {
420 }
421
FindMinTileValueFindMinTileValue422 FindMinTileValue(FindMinTileValue& rhs, tbb::split)
423 : minValue(std::numeric_limits<ValueType>::max())
424 , mNodes(rhs.mNodes)
425 {
426 }
427
operatorFindMinTileValue428 void operator()(const tbb::blocked_range<size_t>& range) {
429 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
430 typename InternalNodeType::ValueAllCIter it = mNodes[n]->beginValueAll();
431 for (; it; ++it) {
432 minValue = std::min(minValue, *it);
433 }
434 }
435 }
436
joinFindMinTileValue437 void join(FindMinTileValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
438
439 ValueType minValue;
440
441 InternalNodeType const * const * const mNodes;
442 }; // FindMinTileValue
443
444
445 template<typename LeafNodeType>
446 struct SDFVoxelsToFogVolume {
447
448 using ValueType = typename LeafNodeType::ValueType;
449
SDFVoxelsToFogVolumeSDFVoxelsToFogVolume450 SDFVoxelsToFogVolume(LeafNodeType ** nodes, ValueType cutoffDistance)
451 : mNodes(nodes), mWeight(ValueType(1.0) / cutoffDistance)
452 {
453 }
454
operatorSDFVoxelsToFogVolume455 void operator()(const tbb::blocked_range<size_t>& range) const {
456
457 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
458
459 LeafNodeType& node = *mNodes[n];
460 node.setValuesOff();
461
462 ValueType* values = node.buffer().data();
463 for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
464 values[i] = values[i] > ValueType(0.0) ? ValueType(0.0) : values[i] * mWeight;
465 if (values[i] > ValueType(0.0)) node.setValueOn(i);
466 }
467
468 if (node.onVoxelCount() == 0) {
469 delete mNodes[n];
470 mNodes[n] = nullptr;
471 }
472 }
473 }
474
475 LeafNodeType ** const mNodes;
476 ValueType const mWeight;
477 }; // SDFVoxelsToFogVolume
478
479
480 template<typename TreeType, typename InternalNodeType>
481 struct SDFTilesToFogVolume {
482
SDFTilesToFogVolumeSDFTilesToFogVolume483 SDFTilesToFogVolume(const TreeType& tree, InternalNodeType ** nodes)
484 : mTree(&tree), mNodes(nodes) { }
485
operatorSDFTilesToFogVolume486 void operator()(const tbb::blocked_range<size_t>& range) const {
487
488 using ValueType = typename TreeType::ValueType;
489 tree::ValueAccessor<const TreeType> acc(*mTree);
490
491 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
492 typename InternalNodeType::ValueAllIter it = mNodes[n]->beginValueAll();
493 for (; it; ++it) {
494 if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
495 it.setValue(ValueType(1.0));
496 it.setValueOn(true);
497 }
498 }
499 }
500 }
501
502 TreeType const * const mTree;
503 InternalNodeType ** const mNodes;
504 }; // SDFTilesToFogVolume
505
506
507 template<typename TreeType>
508 struct FillMaskBoundary {
509
510 using ValueType = typename TreeType::ValueType;
511 using LeafNodeType = typename TreeType::LeafNodeType;
512 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
513 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
514
FillMaskBoundaryFillMaskBoundary515 FillMaskBoundary(const TreeType& tree, ValueType isovalue, const BoolTreeType& fillMask,
516 const BoolLeafNodeType ** fillNodes, BoolLeafNodeType ** newNodes)
517 : mTree(&tree)
518 , mFillMask(&fillMask)
519 , mFillNodes(fillNodes)
520 , mNewNodes(newNodes)
521 , mIsovalue(isovalue)
522 {
523 }
524
operatorFillMaskBoundary525 void operator()(const tbb::blocked_range<size_t>& range) const {
526
527 tree::ValueAccessor<const BoolTreeType> maskAcc(*mFillMask);
528 tree::ValueAccessor<const TreeType> distAcc(*mTree);
529
530 std::unique_ptr<char[]> valueMask(new char[BoolLeafNodeType::SIZE]);
531
532 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
533
534 mNewNodes[n] = nullptr;
535 const BoolLeafNodeType& node = *mFillNodes[n];
536 const Coord& origin = node.origin();
537
538 const bool denseNode = node.isDense();
539
540 // possible early out if the fill mask is dense
541 if (denseNode) {
542
543 int denseNeighbors = 0;
544
545 const BoolLeafNodeType* neighborNode =
546 maskAcc.probeConstLeaf(origin.offsetBy(-1, 0, 0));
547 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
548
549 neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(BoolLeafNodeType::DIM, 0, 0));
550 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
551
552 neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, -1, 0));
553 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
554
555 neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, BoolLeafNodeType::DIM, 0));
556 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
557
558 neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, -1));
559 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
560
561 neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, BoolLeafNodeType::DIM));
562 if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
563
564 if (denseNeighbors == 6) continue;
565 }
566
567 // rest value mask
568 memset(valueMask.get(), 0, sizeof(char) * BoolLeafNodeType::SIZE);
569
570 const typename TreeType::LeafNodeType* distNode = distAcc.probeConstLeaf(origin);
571
572 // check internal voxel neighbors
573
574 bool earlyTermination = false;
575
576 if (!denseNode) {
577 if (distNode) {
578 evalInternalNeighborsP(valueMask.get(), node, *distNode);
579 evalInternalNeighborsN(valueMask.get(), node, *distNode);
580 } else if (distAcc.getValue(origin) > mIsovalue) {
581 earlyTermination = evalInternalNeighborsP(valueMask.get(), node);
582 if (!earlyTermination) {
583 earlyTermination = evalInternalNeighborsN(valueMask.get(), node);
584 }
585 }
586 }
587
588 // check external voxel neighbors
589
590 if (!earlyTermination) {
591 evalExternalNeighborsX<true>(valueMask.get(), node, maskAcc, distAcc);
592 evalExternalNeighborsX<false>(valueMask.get(), node, maskAcc, distAcc);
593 evalExternalNeighborsY<true>(valueMask.get(), node, maskAcc, distAcc);
594 evalExternalNeighborsY<false>(valueMask.get(), node, maskAcc, distAcc);
595 evalExternalNeighborsZ<true>(valueMask.get(), node, maskAcc, distAcc);
596 evalExternalNeighborsZ<false>(valueMask.get(), node, maskAcc, distAcc);
597 }
598
599 // Export marked boundary voxels.
600
601 int numBoundaryValues = 0;
602 for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
603 numBoundaryValues += valueMask[i] == 1;
604 }
605
606 if (numBoundaryValues > 0) {
607 mNewNodes[n] = new BoolLeafNodeType(origin, false);
608 for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
609 if (valueMask[i] == 1) mNewNodes[n]->setValueOn(i);
610 }
611 }
612 }
613 }
614
615 private:
616 // Check internal voxel neighbors in positive {x, y, z} directions.
evalInternalNeighborsPFillMaskBoundary617 void evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node,
618 const LeafNodeType& distNode) const
619 {
620 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
621 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
622 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
623 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
624 for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
625 const Index pos = yPos + z;
626
627 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
628
629 if (!node.isValueOn(pos + 1) && distNode.getValue(pos + 1) > mIsovalue) {
630 valueMask[pos] = 1;
631 }
632 }
633 }
634 }
635
636 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
637 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
638 for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
639 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
640 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
641 const Index pos = yPos + z;
642
643 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
644
645 if (!node.isValueOn(pos + BoolLeafNodeType::DIM) &&
646 distNode.getValue(pos + BoolLeafNodeType::DIM) > mIsovalue) {
647 valueMask[pos] = 1;
648 }
649 }
650 }
651 }
652
653 for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
654 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
655 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
656 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
657 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
658 const Index pos = yPos + z;
659
660 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
661
662 if (!node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
663 (distNode.getValue(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)
664 > mIsovalue))
665 {
666 valueMask[pos] = 1;
667 }
668 }
669 }
670 }
671 }
672
evalInternalNeighborsPFillMaskBoundary673 bool evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node) const {
674
675 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
676 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
677 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
678 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
679 for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
680 const Index pos = yPos + z;
681
682 if (node.isValueOn(pos) && !node.isValueOn(pos + 1)) {
683 valueMask[pos] = 1;
684 return true;
685 }
686 }
687 }
688 }
689
690 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
691 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
692 for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
693 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
694 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
695 const Index pos = yPos + z;
696
697 if (node.isValueOn(pos) && !node.isValueOn(pos + BoolLeafNodeType::DIM)) {
698 valueMask[pos] = 1;
699 return true;
700 }
701 }
702 }
703 }
704
705 for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
706 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
707 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
708 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
709 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
710 const Index pos = yPos + z;
711
712 if (node.isValueOn(pos) &&
713 !node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
714 valueMask[pos] = 1;
715 return true;
716 }
717 }
718 }
719 }
720
721 return false;
722 }
723
724 // Check internal voxel neighbors in negative {x, y, z} directions.
725
evalInternalNeighborsNFillMaskBoundary726 void evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node,
727 const LeafNodeType& distNode) const
728 {
729 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
730 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
731 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
732 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
733 for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
734 const Index pos = yPos + z;
735
736 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
737
738 if (!node.isValueOn(pos - 1) && distNode.getValue(pos - 1) > mIsovalue) {
739 valueMask[pos] = 1;
740 }
741 }
742 }
743 }
744
745 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
746 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
747 for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
748 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
749 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
750 const Index pos = yPos + z;
751
752 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
753
754 if (!node.isValueOn(pos - BoolLeafNodeType::DIM) &&
755 distNode.getValue(pos - BoolLeafNodeType::DIM) > mIsovalue) {
756 valueMask[pos] = 1;
757 }
758 }
759 }
760 }
761
762 for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
763 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
764 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
765 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
766 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
767 const Index pos = yPos + z;
768
769 if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
770
771 if (!node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
772 (distNode.getValue(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)
773 > mIsovalue))
774 {
775 valueMask[pos] = 1;
776 }
777 }
778 }
779 }
780 }
781
782
evalInternalNeighborsNFillMaskBoundary783 bool evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node) const {
784
785 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
786 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
787 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
788 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
789 for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
790 const Index pos = yPos + z;
791
792 if (node.isValueOn(pos) && !node.isValueOn(pos - 1)) {
793 valueMask[pos] = 1;
794 return true;
795 }
796 }
797 }
798 }
799
800 for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
801 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
802 for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
803 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
804 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
805 const Index pos = yPos + z;
806
807 if (node.isValueOn(pos) && !node.isValueOn(pos - BoolLeafNodeType::DIM)) {
808 valueMask[pos] = 1;
809 return true;
810 }
811 }
812 }
813 }
814
815 for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
816 const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
817 for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
818 const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
819 for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
820 const Index pos = yPos + z;
821
822 if (node.isValueOn(pos) &&
823 !node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
824 valueMask[pos] = 1;
825 return true;
826 }
827 }
828 }
829 }
830
831 return false;
832 }
833
834
835 // Check external voxel neighbors
836
837 // If UpWind is true check the X+ oriented node face, else the X- oriented face.
838 template<bool UpWind>
evalExternalNeighborsXFillMaskBoundary839 void evalExternalNeighborsX(char* valueMask, const BoolLeafNodeType& node,
840 const tree::ValueAccessor<const BoolTreeType>& maskAcc,
841 const tree::ValueAccessor<const TreeType>& distAcc) const {
842
843 const Coord& origin = node.origin();
844 Coord ijk(0, 0, 0), nijk;
845 int step = -1;
846
847 if (UpWind) {
848 step = 1;
849 ijk[0] = int(BoolLeafNodeType::DIM) - 1;
850 }
851
852 const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
853
854 for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
855 const Index yPos = xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
856
857 for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
858 const Index pos = yPos + ijk[2];
859
860 if (valueMask[pos] == 0 && node.isValueOn(pos)) {
861
862 nijk = origin + ijk.offsetBy(step, 0, 0);
863
864 if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
865 valueMask[pos] = 1;
866 }
867 }
868 }
869 }
870 }
871
872 // If UpWind is true check the Y+ oriented node face, else the Y- oriented face.
873 template<bool UpWind>
evalExternalNeighborsYFillMaskBoundary874 void evalExternalNeighborsY(char* valueMask, const BoolLeafNodeType& node,
875 const tree::ValueAccessor<const BoolTreeType>& maskAcc,
876 const tree::ValueAccessor<const TreeType>& distAcc) const {
877
878 const Coord& origin = node.origin();
879 Coord ijk(0, 0, 0), nijk;
880 int step = -1;
881
882 if (UpWind) {
883 step = 1;
884 ijk[1] = int(BoolLeafNodeType::DIM) - 1;
885 }
886
887 const Index yPos = ijk[1] << int(BoolLeafNodeType::LOG2DIM);
888
889 for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
890 const Index xPos = yPos + (ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM)));
891
892 for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
893 const Index pos = xPos + ijk[2];
894
895 if (valueMask[pos] == 0 && node.isValueOn(pos)) {
896
897 nijk = origin + ijk.offsetBy(0, step, 0);
898 if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
899 valueMask[pos] = 1;
900 }
901 }
902 }
903 }
904 }
905
906 // If UpWind is true check the Z+ oriented node face, else the Z- oriented face.
907 template<bool UpWind>
evalExternalNeighborsZFillMaskBoundary908 void evalExternalNeighborsZ(char* valueMask, const BoolLeafNodeType& node,
909 const tree::ValueAccessor<const BoolTreeType>& maskAcc,
910 const tree::ValueAccessor<const TreeType>& distAcc) const {
911
912 const Coord& origin = node.origin();
913 Coord ijk(0, 0, 0), nijk;
914 int step = -1;
915
916 if (UpWind) {
917 step = 1;
918 ijk[2] = int(BoolLeafNodeType::DIM) - 1;
919 }
920
921 for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
922 const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
923
924 for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
925 const Index pos = ijk[2] + xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
926
927 if (valueMask[pos] == 0 && node.isValueOn(pos)) {
928
929 nijk = origin + ijk.offsetBy(0, 0, step);
930 if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
931 valueMask[pos] = 1;
932 }
933 }
934 }
935 }
936 }
937
938 //////////
939
940 TreeType const * const mTree;
941 BoolTreeType const * const mFillMask;
942 BoolLeafNodeType const * const * const mFillNodes;
943 BoolLeafNodeType ** const mNewNodes;
944 ValueType const mIsovalue;
945 }; // FillMaskBoundary
946
947
948 /// @brief Constructs a memory light char tree that represents the exterior region with @c +1
949 /// and the interior regions with @c -1.
950 template <class TreeType>
951 typename TreeType::template ValueConverter<char>::Type::Ptr
computeEnclosedRegionMask(const TreeType & tree,typename TreeType::ValueType isovalue,const typename TreeType::template ValueConverter<bool>::Type * fillMask)952 computeEnclosedRegionMask(const TreeType& tree, typename TreeType::ValueType isovalue,
953 const typename TreeType::template ValueConverter<bool>::Type* fillMask)
954 {
955 using LeafNodeType = typename TreeType::LeafNodeType;
956 using RootNodeType = typename TreeType::RootNodeType;
957 using NodeChainType = typename RootNodeType::NodeChainType;
958 using InternalNodeType = typename NodeChainType::template Get<1>;
959
960 using CharTreeType = typename TreeType::template ValueConverter<char>::Type;
961 using CharLeafNodeType = typename CharTreeType::LeafNodeType;
962
963 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
964 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
965
966 const TreeType* treePt = &tree;
967
968 size_t numLeafNodes = 0, numInternalNodes = 0;
969
970 std::vector<const LeafNodeType*> nodes;
971 std::vector<size_t> leafnodeCount;
972
973 {
974 // compute the prefix sum of the leafnode count in each internal node.
975 std::vector<const InternalNodeType*> internalNodes;
976 treePt->getNodes(internalNodes);
977
978 numInternalNodes = internalNodes.size();
979
980 leafnodeCount.push_back(0);
981 for (size_t n = 0; n < numInternalNodes; ++n) {
982 leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
983 }
984
985 numLeafNodes = leafnodeCount.back();
986
987 // extract all leafnodes
988 nodes.reserve(numLeafNodes);
989
990 for (size_t n = 0; n < numInternalNodes; ++n) {
991 internalNodes[n]->getNodes(nodes);
992 }
993 }
994
995 // create mask leafnodes
996 std::unique_ptr<CharLeafNodeType*[]> maskNodes(new CharLeafNodeType*[numLeafNodes]);
997
998 tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
999 LabelBoundaryVoxels<LeafNodeType>(isovalue, nodes.data(), maskNodes.get()));
1000
1001 // create mask grid
1002 typename CharTreeType::Ptr maskTree(new CharTreeType(1));
1003
1004 PopulateTree<CharTreeType> populate(*maskTree, maskNodes.get(), leafnodeCount.data(), 1);
1005 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1006
1007 // optionally evaluate the fill mask
1008
1009 std::vector<CharLeafNodeType*> extraMaskNodes;
1010
1011 if (fillMask) {
1012
1013 std::vector<const BoolLeafNodeType*> fillMaskNodes;
1014 fillMask->getNodes(fillMaskNodes);
1015
1016 std::unique_ptr<BoolLeafNodeType*[]> boundaryMaskNodes(
1017 new BoolLeafNodeType*[fillMaskNodes.size()]);
1018
1019 tbb::parallel_for(tbb::blocked_range<size_t>(0, fillMaskNodes.size()),
1020 FillMaskBoundary<TreeType>(tree, isovalue, *fillMask, fillMaskNodes.data(),
1021 boundaryMaskNodes.get()));
1022
1023 tree::ValueAccessor<CharTreeType> maskAcc(*maskTree);
1024
1025 for (size_t n = 0, N = fillMaskNodes.size(); n < N; ++n) {
1026
1027 if (boundaryMaskNodes[n] == nullptr) continue;
1028
1029 const BoolLeafNodeType& boundaryNode = *boundaryMaskNodes[n];
1030 const Coord& origin = boundaryNode.origin();
1031
1032 CharLeafNodeType* maskNodePt = maskAcc.probeLeaf(origin);
1033
1034 if (!maskNodePt) {
1035 maskNodePt = maskAcc.touchLeaf(origin);
1036 extraMaskNodes.push_back(maskNodePt);
1037 }
1038
1039 char* data = maskNodePt->buffer().data();
1040
1041 typename BoolLeafNodeType::ValueOnCIter it = boundaryNode.cbeginValueOn();
1042 for (; it; ++it) {
1043 if (data[it.pos()] != 0) data[it.pos()] = -1;
1044 }
1045
1046 delete boundaryMaskNodes[n];
1047 }
1048 }
1049
1050 // eliminate enclosed regions
1051 tools::traceExteriorBoundaries(*maskTree);
1052
1053 // flip voxel sign to negative inside and positive outside.
1054 tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1055 FlipRegionSign<CharLeafNodeType>(maskNodes.get()));
1056
1057 if (!extraMaskNodes.empty()) {
1058 tbb::parallel_for(tbb::blocked_range<size_t>(0, extraMaskNodes.size()),
1059 FlipRegionSign<CharLeafNodeType>(extraMaskNodes.data()));
1060 }
1061
1062 // propagate sign information into tile region
1063 tools::signedFloodFill(*maskTree);
1064
1065 return maskTree;
1066 } // computeEnclosedRegionMask()
1067
1068
1069 template <class TreeType>
1070 typename TreeType::template ValueConverter<bool>::Type::Ptr
computeInteriorMask(const TreeType & tree,typename TreeType::ValueType iso)1071 computeInteriorMask(const TreeType& tree, typename TreeType::ValueType iso)
1072 {
1073 using ValueType = typename TreeType::ValueType;
1074 using LeafNodeType = typename TreeType::LeafNodeType;
1075 using RootNodeType = typename TreeType::RootNodeType;
1076 using NodeChainType = typename RootNodeType::NodeChainType;
1077 using InternalNodeType = typename NodeChainType::template Get<1>;
1078
1079 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1080 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1081 using BoolRootNodeType = typename BoolTreeType::RootNodeType;
1082 using BoolNodeChainType = typename BoolRootNodeType::NodeChainType;
1083 using BoolInternalNodeType = typename BoolNodeChainType::template Get<1>;
1084
1085 /////
1086
1087 // Clamp the isovalue to the level set's background value minus epsilon.
1088 // (In a valid narrow-band level set, all voxels, including background voxels,
1089 // have values less than or equal to the background value, so an isovalue
1090 // greater than or equal to the background value would produce a mask with
1091 // effectively infinite extent.)
1092 iso = std::min(iso,
1093 static_cast<ValueType>(tree.background() - math::Tolerance<ValueType>::value()));
1094
1095 size_t numLeafNodes = 0, numInternalNodes = 0;
1096
1097 std::vector<const LeafNodeType*> nodes;
1098 std::vector<size_t> leafnodeCount;
1099
1100 {
1101 // compute the prefix sum of the leafnode count in each internal node.
1102 std::vector<const InternalNodeType*> internalNodes;
1103 tree.getNodes(internalNodes);
1104
1105 numInternalNodes = internalNodes.size();
1106
1107 leafnodeCount.push_back(0);
1108 for (size_t n = 0; n < numInternalNodes; ++n) {
1109 leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
1110 }
1111
1112 numLeafNodes = leafnodeCount.back();
1113
1114 // extract all leafnodes
1115 nodes.reserve(numLeafNodes);
1116
1117 for (size_t n = 0; n < numInternalNodes; ++n) {
1118 internalNodes[n]->getNodes(nodes);
1119 }
1120 }
1121
1122 // create mask leafnodes
1123 std::unique_ptr<BoolLeafNodeType*[]> maskNodes(new BoolLeafNodeType*[numLeafNodes]);
1124
1125 tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1126 MaskInteriorVoxels<LeafNodeType>(iso, nodes.data(), maskNodes.get()));
1127
1128
1129 // create mask grid
1130 typename BoolTreeType::Ptr maskTree(new BoolTreeType(false));
1131
1132 PopulateTree<BoolTreeType> populate(*maskTree, maskNodes.get(), leafnodeCount.data(), false);
1133 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1134
1135
1136 // evaluate tile values
1137 std::vector<BoolInternalNodeType*> internalMaskNodes;
1138 maskTree->getNodes(internalMaskNodes);
1139
1140 tbb::parallel_for(tbb::blocked_range<size_t>(0, internalMaskNodes.size()),
1141 MaskInteriorTiles<TreeType, BoolInternalNodeType>(iso, tree, internalMaskNodes.data()));
1142
1143 tree::ValueAccessor<const TreeType> acc(tree);
1144
1145 typename BoolTreeType::ValueAllIter it(*maskTree);
1146 it.setMaxDepth(BoolTreeType::ValueAllIter::LEAF_DEPTH - 2);
1147
1148 for ( ; it; ++it) {
1149 if (acc.getValue(it.getCoord()) < iso) {
1150 it.setValue(true);
1151 it.setActiveState(true);
1152 }
1153 }
1154
1155 return maskTree;
1156 } // computeInteriorMask()
1157
1158
1159 template<typename InputTreeType>
1160 struct MaskIsovalueCrossingVoxels
1161 {
1162 using InputValueType = typename InputTreeType::ValueType;
1163 using InputLeafNodeType = typename InputTreeType::LeafNodeType;
1164 using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
1165 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1166
MaskIsovalueCrossingVoxelsMaskIsovalueCrossingVoxels1167 MaskIsovalueCrossingVoxels(
1168 const InputTreeType& inputTree,
1169 const std::vector<const InputLeafNodeType*>& inputLeafNodes,
1170 BoolTreeType& maskTree,
1171 InputValueType iso)
1172 : mInputAccessor(inputTree)
1173 , mInputNodes(!inputLeafNodes.empty() ? &inputLeafNodes.front() : nullptr)
1174 , mMaskTree(false)
1175 , mMaskAccessor(maskTree)
1176 , mIsovalue(iso)
1177 {
1178 }
1179
MaskIsovalueCrossingVoxelsMaskIsovalueCrossingVoxels1180 MaskIsovalueCrossingVoxels(MaskIsovalueCrossingVoxels& rhs, tbb::split)
1181 : mInputAccessor(rhs.mInputAccessor.tree())
1182 , mInputNodes(rhs.mInputNodes)
1183 , mMaskTree(false)
1184 , mMaskAccessor(mMaskTree)
1185 , mIsovalue(rhs.mIsovalue)
1186 {
1187 }
1188
operatorMaskIsovalueCrossingVoxels1189 void operator()(const tbb::blocked_range<size_t>& range) {
1190
1191 const InputValueType iso = mIsovalue;
1192 Coord ijk(0, 0, 0);
1193
1194 BoolLeafNodeType* maskNodePt = nullptr;
1195
1196 for (size_t n = range.begin(); mInputNodes && (n != range.end()); ++n) {
1197
1198 const InputLeafNodeType& node = *mInputNodes[n];
1199
1200 if (!maskNodePt) maskNodePt = new BoolLeafNodeType(node.origin(), false);
1201 else maskNodePt->setOrigin(node.origin());
1202
1203 bool collectedData = false;
1204
1205 for (typename InputLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
1206
1207 bool isUnder = *it < iso;
1208
1209 ijk = it.getCoord();
1210
1211 ++ijk[2];
1212 bool signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +z edge
1213 --ijk[2];
1214
1215 if (!signChange) {
1216 --ijk[2];
1217 signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -z edge
1218 ++ijk[2];
1219 }
1220
1221 if (!signChange) {
1222 ++ijk[1];
1223 signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +y edge
1224 --ijk[1];
1225 }
1226
1227 if (!signChange) {
1228 --ijk[1];
1229 signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -y edge
1230 ++ijk[1];
1231 }
1232
1233 if (!signChange) {
1234 ++ijk[0];
1235 signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +x edge
1236 --ijk[0];
1237 }
1238
1239 if (!signChange) {
1240 --ijk[0];
1241 signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -x edge
1242 ++ijk[0];
1243 }
1244
1245 if (signChange) {
1246 collectedData = true;
1247 maskNodePt->setValueOn(it.pos(), true);
1248 }
1249 }
1250
1251 if (collectedData) {
1252 mMaskAccessor.addLeaf(maskNodePt);
1253 maskNodePt = nullptr;
1254 }
1255 }
1256
1257 if (maskNodePt) delete maskNodePt;
1258 }
1259
joinMaskIsovalueCrossingVoxels1260 void join(MaskIsovalueCrossingVoxels& rhs) {
1261 mMaskAccessor.tree().merge(rhs.mMaskAccessor.tree());
1262 }
1263
1264 private:
1265 tree::ValueAccessor<const InputTreeType> mInputAccessor;
1266 InputLeafNodeType const * const * const mInputNodes;
1267
1268 BoolTreeType mMaskTree;
1269 tree::ValueAccessor<BoolTreeType> mMaskAccessor;
1270
1271 InputValueType mIsovalue;
1272 }; // MaskIsovalueCrossingVoxels
1273
1274
1275 ////////////////////////////////////////
1276
1277
1278 template<typename NodeType>
1279 struct NodeMaskSegment
1280 {
1281 using Ptr = SharedPtr<NodeMaskSegment>;
1282 using NodeMaskType = typename NodeType::NodeMaskType;
1283
NodeMaskSegmentNodeMaskSegment1284 NodeMaskSegment() : connections(), mask(false), origin(0,0,0), visited(false) {}
1285
1286 std::vector<NodeMaskSegment*> connections;
1287 NodeMaskType mask;
1288 Coord origin;
1289 bool visited;
1290 }; // struct NodeMaskSegment
1291
1292
1293 template<typename NodeType>
1294 void
nodeMaskSegmentation(const NodeType & node,std::vector<typename NodeMaskSegment<NodeType>::Ptr> & segments)1295 nodeMaskSegmentation(const NodeType& node,
1296 std::vector<typename NodeMaskSegment<NodeType>::Ptr>& segments)
1297 {
1298 using NodeMaskType = typename NodeType::NodeMaskType;
1299 using NodeMaskSegmentType = NodeMaskSegment<NodeType>;
1300 using NodeMaskSegmentTypePtr = typename NodeMaskSegmentType::Ptr;
1301
1302 NodeMaskType nodeMask(node.getValueMask());
1303 std::deque<Index> indexList;
1304
1305 while (!nodeMask.isOff()) {
1306
1307 NodeMaskSegmentTypePtr segment(new NodeMaskSegmentType());
1308 segment->origin = node.origin();
1309
1310 NodeMaskType& mask = segment->mask;
1311
1312 indexList.push_back(nodeMask.findFirstOn());
1313 nodeMask.setOff(indexList.back()); // mark as visited
1314 Coord ijk(0, 0, 0);
1315
1316 while (!indexList.empty()) {
1317
1318 const Index pos = indexList.back();
1319 indexList.pop_back();
1320
1321 if (mask.isOn(pos)) continue;
1322 mask.setOn(pos);
1323
1324 ijk = NodeType::offsetToLocalCoord(pos);
1325
1326 Index npos = pos - 1;
1327 if (ijk[2] != 0 && nodeMask.isOn(npos)) {
1328 nodeMask.setOff(npos);
1329 indexList.push_back(npos);
1330 }
1331
1332 npos = pos + 1;
1333 if (ijk[2] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1334 nodeMask.setOff(npos);
1335 indexList.push_back(npos);
1336 }
1337
1338 npos = pos - NodeType::DIM;
1339 if (ijk[1] != 0 && nodeMask.isOn(npos)) {
1340 nodeMask.setOff(npos);
1341 indexList.push_back(npos);
1342 }
1343
1344 npos = pos + NodeType::DIM;
1345 if (ijk[1] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1346 nodeMask.setOff(npos);
1347 indexList.push_back(npos);
1348 }
1349
1350 npos = pos - NodeType::DIM * NodeType::DIM;
1351 if (ijk[0] != 0 && nodeMask.isOn(npos)) {
1352 nodeMask.setOff(npos);
1353 indexList.push_back(npos);
1354 }
1355
1356 npos = pos + NodeType::DIM * NodeType::DIM;
1357 if (ijk[0] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1358 nodeMask.setOff(npos);
1359 indexList.push_back(npos);
1360 }
1361
1362 }
1363
1364 segments.push_back(segment);
1365 }
1366 }
1367
1368
1369 template<typename NodeType>
1370 struct SegmentNodeMask
1371 {
1372 using NodeMaskSegmentType = NodeMaskSegment<NodeType>;
1373 using NodeMaskSegmentTypePtr = typename NodeMaskSegmentType::Ptr;
1374 using NodeMaskSegmentVector = typename std::vector<NodeMaskSegmentTypePtr>;
1375
SegmentNodeMaskSegmentNodeMask1376 SegmentNodeMask(std::vector<NodeType*>& nodes, NodeMaskSegmentVector* nodeMaskArray)
1377 : mNodes(!nodes.empty() ? &nodes.front() : nullptr)
1378 , mNodeMaskArray(nodeMaskArray)
1379 {
1380 }
1381
operatorSegmentNodeMask1382 void operator()(const tbb::blocked_range<size_t>& range) const {
1383 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1384 NodeType& node = *mNodes[n];
1385 nodeMaskSegmentation(node, mNodeMaskArray[n]);
1386
1387 // hack origin data to store array offset
1388 Coord& origin = const_cast<Coord&>(node.origin());
1389 origin[0] = static_cast<int>(n);
1390 }
1391 }
1392
1393 NodeType * const * const mNodes;
1394 NodeMaskSegmentVector * const mNodeMaskArray;
1395 }; // struct SegmentNodeMask
1396
1397
1398 template<typename TreeType, typename NodeType>
1399 struct ConnectNodeMaskSegments
1400 {
1401 using NodeMaskType = typename NodeType::NodeMaskType;
1402 using NodeMaskSegmentType = NodeMaskSegment<NodeType>;
1403 using NodeMaskSegmentTypePtr = typename NodeMaskSegmentType::Ptr;
1404 using NodeMaskSegmentVector = typename std::vector<NodeMaskSegmentTypePtr>;
1405
ConnectNodeMaskSegmentsConnectNodeMaskSegments1406 ConnectNodeMaskSegments(const TreeType& tree, NodeMaskSegmentVector* nodeMaskArray)
1407 : mTree(&tree)
1408 , mNodeMaskArray(nodeMaskArray)
1409 {
1410 }
1411
operatorConnectNodeMaskSegments1412 void operator()(const tbb::blocked_range<size_t>& range) const {
1413
1414 tree::ValueAccessor<const TreeType> acc(*mTree);
1415
1416 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1417
1418 NodeMaskSegmentVector& segments = mNodeMaskArray[n];
1419 if (segments.empty()) continue;
1420
1421 std::vector<std::set<NodeMaskSegmentType*> > connections(segments.size());
1422
1423 Coord ijk = segments[0]->origin;
1424
1425 const NodeType* node = acc.template probeConstNode<NodeType>(ijk);
1426 if (!node) continue;
1427
1428 // get neighbour nodes
1429
1430 ijk[2] += NodeType::DIM;
1431 const NodeType* nodeZUp = acc.template probeConstNode<NodeType>(ijk);
1432 ijk[2] -= (NodeType::DIM + NodeType::DIM);
1433 const NodeType* nodeZDown = acc.template probeConstNode<NodeType>(ijk);
1434 ijk[2] += NodeType::DIM;
1435
1436 ijk[1] += NodeType::DIM;
1437 const NodeType* nodeYUp = acc.template probeConstNode<NodeType>(ijk);
1438 ijk[1] -= (NodeType::DIM + NodeType::DIM);
1439 const NodeType* nodeYDown = acc.template probeConstNode<NodeType>(ijk);
1440 ijk[1] += NodeType::DIM;
1441
1442 ijk[0] += NodeType::DIM;
1443 const NodeType* nodeXUp = acc.template probeConstNode<NodeType>(ijk);
1444 ijk[0] -= (NodeType::DIM + NodeType::DIM);
1445 const NodeType* nodeXDown = acc.template probeConstNode<NodeType>(ijk);
1446 ijk[0] += NodeType::DIM;
1447
1448 const Index startPos = node->getValueMask().findFirstOn();
1449 for (Index pos = startPos; pos < NodeMaskType::SIZE; ++pos) {
1450
1451 if (!node->isValueOn(pos)) continue;
1452
1453 ijk = NodeType::offsetToLocalCoord(pos);
1454
1455 #ifdef _MSC_FULL_VER
1456 #if _MSC_FULL_VER >= 190000000 && _MSC_FULL_VER < 190024210
1457 // Visual Studio 2015 had a codegen bug that wasn't fixed until Update 3
1458 volatile Index npos = 0;
1459 #else
1460 Index npos = 0;
1461 #endif
1462 #else
1463 Index npos = 0;
1464 #endif
1465
1466 if (ijk[2] == 0) {
1467 npos = pos + (NodeType::DIM - 1);
1468 if (nodeZDown && nodeZDown->isValueOn(npos)) {
1469 NodeMaskSegmentType* nsegment =
1470 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZDown)], npos);
1471 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1472 connections[idx].insert(nsegment);
1473 }
1474 } else if (ijk[2] == (NodeType::DIM - 1)) {
1475 npos = pos - (NodeType::DIM - 1);
1476 if (nodeZUp && nodeZUp->isValueOn(npos)) {
1477 NodeMaskSegmentType* nsegment =
1478 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZUp)], npos);
1479 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1480 connections[idx].insert(nsegment);
1481 }
1482 }
1483
1484 if (ijk[1] == 0) {
1485 npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1486 if (nodeYDown && nodeYDown->isValueOn(npos)) {
1487 NodeMaskSegmentType* nsegment =
1488 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYDown)], npos);
1489 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1490 connections[idx].insert(nsegment);
1491 }
1492 } else if (ijk[1] == (NodeType::DIM - 1)) {
1493 npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1494 if (nodeYUp && nodeYUp->isValueOn(npos)) {
1495 NodeMaskSegmentType* nsegment =
1496 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYUp)], npos);
1497 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1498 connections[idx].insert(nsegment);
1499 }
1500 }
1501
1502 if (ijk[0] == 0) {
1503 npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1504 if (nodeXDown && nodeXDown->isValueOn(npos)) {
1505 NodeMaskSegmentType* nsegment =
1506 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXDown)], npos);
1507 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1508 connections[idx].insert(nsegment);
1509 }
1510 } else if (ijk[0] == (NodeType::DIM - 1)) {
1511 npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1512 if (nodeXUp && nodeXUp->isValueOn(npos)) {
1513 NodeMaskSegmentType* nsegment =
1514 findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXUp)], npos);
1515 const Index idx = findNodeMaskSegmentIndex(segments, pos);
1516 connections[idx].insert(nsegment);
1517 }
1518 }
1519 }
1520
1521 for (size_t i = 0, I = connections.size(); i < I; ++i) {
1522
1523 typename std::set<NodeMaskSegmentType*>::iterator
1524 it = connections[i].begin(), end = connections[i].end();
1525
1526 std::vector<NodeMaskSegmentType*>& segmentConnections = segments[i]->connections;
1527 segmentConnections.reserve(connections.size());
1528 for (; it != end; ++it) {
1529 segmentConnections.push_back(*it);
1530 }
1531 }
1532 } // end range loop
1533 }
1534
1535 private:
1536
getNodeOffsetConnectNodeMaskSegments1537 static inline size_t getNodeOffset(const NodeType& node) {
1538 return static_cast<size_t>(node.origin()[0]);
1539 }
1540
1541 static inline NodeMaskSegmentType*
findNodeMaskSegmentConnectNodeMaskSegments1542 findNodeMaskSegment(NodeMaskSegmentVector& segments, Index pos)
1543 {
1544 NodeMaskSegmentType* segment = nullptr;
1545
1546 for (size_t n = 0, N = segments.size(); n < N; ++n) {
1547 if (segments[n]->mask.isOn(pos)) {
1548 segment = segments[n].get();
1549 break;
1550 }
1551 }
1552
1553 return segment;
1554 }
1555
1556 static inline Index
findNodeMaskSegmentIndexConnectNodeMaskSegments1557 findNodeMaskSegmentIndex(NodeMaskSegmentVector& segments, Index pos)
1558 {
1559 for (Index n = 0, N = Index(segments.size()); n < N; ++n) {
1560 if (segments[n]->mask.isOn(pos)) return n;
1561 }
1562 return Index(-1);
1563 }
1564
1565 TreeType const * const mTree;
1566 NodeMaskSegmentVector * const mNodeMaskArray;
1567 }; // struct ConnectNodeMaskSegments
1568
1569
1570 template<typename TreeType>
1571 struct MaskSegmentGroup
1572 {
1573 using LeafNodeType = typename TreeType::LeafNodeType;
1574 using TreeTypePtr = typename TreeType::Ptr;
1575 using NodeMaskSegmentType = NodeMaskSegment<LeafNodeType>;
1576
MaskSegmentGroupMaskSegmentGroup1577 MaskSegmentGroup(const std::vector<NodeMaskSegmentType*>& segments)
1578 : mSegments(!segments.empty() ? &segments.front() : nullptr)
1579 , mTree(new TreeType(false))
1580 {
1581 }
1582
MaskSegmentGroupMaskSegmentGroup1583 MaskSegmentGroup(const MaskSegmentGroup& rhs, tbb::split)
1584 : mSegments(rhs.mSegments)
1585 , mTree(new TreeType(false))
1586 {
1587 }
1588
maskMaskSegmentGroup1589 TreeTypePtr& mask() { return mTree; }
1590
joinMaskSegmentGroup1591 void join(MaskSegmentGroup& rhs) { mTree->merge(*rhs.mTree); }
1592
operatorMaskSegmentGroup1593 void operator()(const tbb::blocked_range<size_t>& range) {
1594
1595 tree::ValueAccessor<TreeType> acc(*mTree);
1596
1597 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1598 NodeMaskSegmentType& segment = *mSegments[n];
1599 LeafNodeType* node = acc.touchLeaf(segment.origin);
1600 node->getValueMask() |= segment.mask;
1601 }
1602 }
1603
1604 private:
1605 NodeMaskSegmentType * const * const mSegments;
1606 TreeTypePtr mTree;
1607 }; // struct MaskSegmentGroup
1608
1609
1610 ////////////////////////////////////////
1611
1612
1613 template<typename TreeType>
1614 struct ExpandLeafNodeRegion
1615 {
1616 using ValueType = typename TreeType::ValueType;
1617 using LeafNodeType = typename TreeType::LeafNodeType;
1618 using NodeMaskType = typename LeafNodeType::NodeMaskType;
1619
1620 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1621 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1622
1623 /////
1624
ExpandLeafNodeRegionExpandLeafNodeRegion1625 ExpandLeafNodeRegion(const TreeType& distTree, BoolTreeType& maskTree,
1626 std::vector<BoolLeafNodeType*>& maskNodes)
1627 : mDistTree(&distTree)
1628 , mMaskTree(&maskTree)
1629 , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
1630 , mNewMaskTree(false)
1631 {
1632 }
1633
ExpandLeafNodeRegionExpandLeafNodeRegion1634 ExpandLeafNodeRegion(const ExpandLeafNodeRegion& rhs, tbb::split)
1635 : mDistTree(rhs.mDistTree)
1636 , mMaskTree(rhs.mMaskTree)
1637 , mMaskNodes(rhs.mMaskNodes)
1638 , mNewMaskTree(false)
1639 {
1640 }
1641
newMaskTreeExpandLeafNodeRegion1642 BoolTreeType& newMaskTree() { return mNewMaskTree; }
1643
joinExpandLeafNodeRegion1644 void join(ExpandLeafNodeRegion& rhs) { mNewMaskTree.merge(rhs.mNewMaskTree); }
1645
operatorExpandLeafNodeRegion1646 void operator()(const tbb::blocked_range<size_t>& range) {
1647
1648 using NodeType = LeafNodeType;
1649
1650 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1651 tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1652 tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
1653
1654 NodeMaskType maskZUp, maskZDown, maskYUp, maskYDown, maskXUp, maskXDown;
1655
1656 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1657
1658 BoolLeafNodeType& maskNode = *mMaskNodes[n];
1659 if (maskNode.isEmpty()) continue;
1660
1661 Coord ijk = maskNode.origin(), nijk;
1662
1663 const LeafNodeType* distNode = distAcc.probeConstLeaf(ijk);
1664 if (!distNode) continue;
1665
1666 const ValueType *dataZUp = nullptr, *dataZDown = nullptr,
1667 *dataYUp = nullptr, *dataYDown = nullptr,
1668 *dataXUp = nullptr, *dataXDown = nullptr;
1669
1670 ijk[2] += NodeType::DIM;
1671 getData(ijk, distAcc, maskAcc, maskZUp, dataZUp);
1672 ijk[2] -= (NodeType::DIM + NodeType::DIM);
1673 getData(ijk, distAcc, maskAcc, maskZDown, dataZDown);
1674 ijk[2] += NodeType::DIM;
1675
1676 ijk[1] += NodeType::DIM;
1677 getData(ijk, distAcc, maskAcc, maskYUp, dataYUp);
1678 ijk[1] -= (NodeType::DIM + NodeType::DIM);
1679 getData(ijk, distAcc, maskAcc, maskYDown, dataYDown);
1680 ijk[1] += NodeType::DIM;
1681
1682 ijk[0] += NodeType::DIM;
1683 getData(ijk, distAcc, maskAcc, maskXUp, dataXUp);
1684 ijk[0] -= (NodeType::DIM + NodeType::DIM);
1685 getData(ijk, distAcc, maskAcc, maskXDown, dataXDown);
1686 ijk[0] += NodeType::DIM;
1687
1688 for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
1689
1690 const Index pos = it.pos();
1691 const ValueType val = std::abs(distNode->getValue(pos));
1692
1693 ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
1694 nijk = ijk + maskNode.origin();
1695
1696 if (dataZUp && ijk[2] == (BoolLeafNodeType::DIM - 1)) {
1697 const Index npos = pos - (NodeType::DIM - 1);
1698 if (maskZUp.isOn(npos) && std::abs(dataZUp[npos]) > val) {
1699 newMaskAcc.setValueOn(nijk.offsetBy(0, 0, 1));
1700 }
1701 } else if (dataZDown && ijk[2] == 0) {
1702 const Index npos = pos + (NodeType::DIM - 1);
1703 if (maskZDown.isOn(npos) && std::abs(dataZDown[npos]) > val) {
1704 newMaskAcc.setValueOn(nijk.offsetBy(0, 0, -1));
1705 }
1706 }
1707
1708 if (dataYUp && ijk[1] == (BoolLeafNodeType::DIM - 1)) {
1709 const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1710 if (maskYUp.isOn(npos) && std::abs(dataYUp[npos]) > val) {
1711 newMaskAcc.setValueOn(nijk.offsetBy(0, 1, 0));
1712 }
1713 } else if (dataYDown && ijk[1] == 0) {
1714 const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1715 if (maskYDown.isOn(npos) && std::abs(dataYDown[npos]) > val) {
1716 newMaskAcc.setValueOn(nijk.offsetBy(0, -1, 0));
1717 }
1718 }
1719
1720 if (dataXUp && ijk[0] == (BoolLeafNodeType::DIM - 1)) {
1721 const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1722 if (maskXUp.isOn(npos) && std::abs(dataXUp[npos]) > val) {
1723 newMaskAcc.setValueOn(nijk.offsetBy(1, 0, 0));
1724 }
1725 } else if (dataXDown && ijk[0] == 0) {
1726 const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1727 if (maskXDown.isOn(npos) && std::abs(dataXDown[npos]) > val) {
1728 newMaskAcc.setValueOn(nijk.offsetBy(-1, 0, 0));
1729 }
1730 }
1731
1732 } // end value on loop
1733 } // end range loop
1734 }
1735
1736 private:
1737
1738 static inline void
getDataExpandLeafNodeRegion1739 getData(const Coord& ijk, tree::ValueAccessor<const TreeType>& distAcc,
1740 tree::ValueAccessor<const BoolTreeType>& maskAcc, NodeMaskType& mask,
1741 const ValueType*& data)
1742 {
1743 const LeafNodeType* node = distAcc.probeConstLeaf(ijk);
1744 if (node) {
1745 data = node->buffer().data();
1746 mask = node->getValueMask();
1747 const BoolLeafNodeType* maskNodePt = maskAcc.probeConstLeaf(ijk);
1748 if (maskNodePt) mask -= maskNodePt->getValueMask();
1749 }
1750 }
1751
1752 TreeType const * const mDistTree;
1753 BoolTreeType * const mMaskTree;
1754 BoolLeafNodeType ** const mMaskNodes;
1755
1756 BoolTreeType mNewMaskTree;
1757 }; // struct ExpandLeafNodeRegion
1758
1759
1760 template<typename TreeType>
1761 struct FillLeafNodeVoxels
1762 {
1763 using ValueType = typename TreeType::ValueType;
1764 using LeafNodeType = typename TreeType::LeafNodeType;
1765 using NodeMaskType = typename LeafNodeType::NodeMaskType;
1766 using BoolLeafNodeType = tree::LeafNode<bool, LeafNodeType::LOG2DIM>;
1767
FillLeafNodeVoxelsFillLeafNodeVoxels1768 FillLeafNodeVoxels(const TreeType& tree, std::vector<BoolLeafNodeType*>& maskNodes)
1769 : mTree(&tree), mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
1770 {
1771 }
1772
operatorFillLeafNodeVoxels1773 void operator()(const tbb::blocked_range<size_t>& range) const {
1774
1775 tree::ValueAccessor<const TreeType> distAcc(*mTree);
1776
1777 std::vector<Index> indexList;
1778 indexList.reserve(NodeMaskType::SIZE);
1779
1780 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1781
1782 BoolLeafNodeType& maskNode = *mMaskNodes[n];
1783
1784 const LeafNodeType * distNode = distAcc.probeConstLeaf(maskNode.origin());
1785 if (!distNode) continue;
1786
1787 NodeMaskType mask(distNode->getValueMask());
1788 NodeMaskType& narrowbandMask = maskNode.getValueMask();
1789
1790 for (Index pos = narrowbandMask.findFirstOn(); pos < NodeMaskType::SIZE; ++pos) {
1791 if (narrowbandMask.isOn(pos)) indexList.push_back(pos);
1792 }
1793
1794 mask -= narrowbandMask; // bitwise difference
1795 narrowbandMask.setOff();
1796
1797 const ValueType* data = distNode->buffer().data();
1798 Coord ijk(0, 0, 0);
1799
1800 while (!indexList.empty()) {
1801
1802 const Index pos = indexList.back();
1803 indexList.pop_back();
1804
1805 if (narrowbandMask.isOn(pos)) continue;
1806 narrowbandMask.setOn(pos);
1807
1808 const ValueType dist = std::abs(data[pos]);
1809
1810 ijk = LeafNodeType::offsetToLocalCoord(pos);
1811
1812 Index npos = pos - 1;
1813 if (ijk[2] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1814 mask.setOff(npos);
1815 indexList.push_back(npos);
1816 }
1817
1818 npos = pos + 1;
1819 if ((ijk[2] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1820 && std::abs(data[npos]) > dist)
1821 {
1822 mask.setOff(npos);
1823 indexList.push_back(npos);
1824 }
1825
1826 npos = pos - LeafNodeType::DIM;
1827 if (ijk[1] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1828 mask.setOff(npos);
1829 indexList.push_back(npos);
1830 }
1831
1832 npos = pos + LeafNodeType::DIM;
1833 if ((ijk[1] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1834 && std::abs(data[npos]) > dist)
1835 {
1836 mask.setOff(npos);
1837 indexList.push_back(npos);
1838 }
1839
1840 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1841 if (ijk[0] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1842 mask.setOff(npos);
1843 indexList.push_back(npos);
1844 }
1845
1846 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1847 if ((ijk[0] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1848 && std::abs(data[npos]) > dist)
1849 {
1850 mask.setOff(npos);
1851 indexList.push_back(npos);
1852 }
1853 } // end flood fill loop
1854 } // end range loop
1855 }
1856
1857 TreeType const * const mTree;
1858 BoolLeafNodeType ** const mMaskNodes;
1859 }; // FillLeafNodeVoxels
1860
1861
1862 template<typename TreeType>
1863 struct ExpandNarrowbandMask
1864 {
1865 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1866 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1867 using BoolTreeTypePtr = typename BoolTreeType::Ptr;
1868
ExpandNarrowbandMaskExpandNarrowbandMask1869 ExpandNarrowbandMask(const TreeType& tree, std::vector<BoolTreeTypePtr>& segments)
1870 : mTree(&tree), mSegments(!segments.empty() ? &segments.front() : nullptr)
1871 {
1872 }
1873
operatorExpandNarrowbandMask1874 void operator()(const tbb::blocked_range<size_t>& range) const {
1875
1876 const TreeType& distTree = *mTree;
1877 std::vector<BoolLeafNodeType*> nodes;
1878
1879 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1880
1881 BoolTreeType& narrowBandMask = *mSegments[n];
1882
1883 BoolTreeType candidateMask(narrowBandMask, false, TopologyCopy());
1884
1885 while (true) {
1886
1887 nodes.clear();
1888 candidateMask.getNodes(nodes);
1889 if (nodes.empty()) break;
1890
1891 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
1892
1893 tbb::parallel_for(nodeRange, FillLeafNodeVoxels<TreeType>(distTree, nodes));
1894
1895 narrowBandMask.topologyUnion(candidateMask);
1896
1897 ExpandLeafNodeRegion<TreeType> op(distTree, narrowBandMask, nodes);
1898 tbb::parallel_reduce(nodeRange, op);
1899
1900 if (op.newMaskTree().empty()) break;
1901
1902 candidateMask.clear();
1903 candidateMask.merge(op.newMaskTree());
1904 } // end expand loop
1905 } // end range loop
1906 }
1907
1908 TreeType const * const mTree;
1909 BoolTreeTypePtr * const mSegments;
1910 }; // ExpandNarrowbandMask
1911
1912
1913 template<typename TreeType>
1914 struct FloodFillSign
1915 {
1916 using TreeTypePtr = typename TreeType::Ptr;
1917 using ValueType = typename TreeType::ValueType;
1918 using LeafNodeType = typename TreeType::LeafNodeType;
1919 using RootNodeType = typename TreeType::RootNodeType;
1920 using NodeChainType = typename RootNodeType::NodeChainType;
1921 using InternalNodeType = typename NodeChainType::template Get<1>;
1922
FloodFillSignFloodFillSign1923 FloodFillSign(const TreeType& tree, std::vector<TreeTypePtr>& segments)
1924 : mTree(&tree)
1925 , mSegments(!segments.empty() ? &segments.front() : nullptr)
1926 , mMinValue(ValueType(0.0))
1927 {
1928 ValueType minSDFValue = std::numeric_limits<ValueType>::max();
1929
1930 {
1931 std::vector<const InternalNodeType*> nodes;
1932 tree.getNodes(nodes);
1933
1934 if (!nodes.empty()) {
1935 FindMinTileValue<InternalNodeType> minOp(nodes.data());
1936 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1937 minSDFValue = std::min(minSDFValue, minOp.minValue);
1938 }
1939 }
1940
1941 if (minSDFValue > ValueType(0.0)) {
1942 std::vector<const LeafNodeType*> nodes;
1943 tree.getNodes(nodes);
1944 if (!nodes.empty()) {
1945 FindMinVoxelValue<LeafNodeType> minOp(nodes.data());
1946 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1947 minSDFValue = std::min(minSDFValue, minOp.minValue);
1948 }
1949 }
1950
1951 mMinValue = minSDFValue;
1952 }
1953
operatorFloodFillSign1954 void operator()(const tbb::blocked_range<size_t>& range) const {
1955 const ValueType interiorValue = -std::abs(mMinValue);
1956 const ValueType exteriorValue = std::abs(mTree->background());
1957 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1958 tools::signedFloodFillWithValues(*mSegments[n], exteriorValue, interiorValue);
1959 }
1960 }
1961
1962 private:
1963
1964 TreeType const * const mTree;
1965 TreeTypePtr * const mSegments;
1966 ValueType mMinValue;
1967 }; // FloodFillSign
1968
1969
1970 template<typename TreeType>
1971 struct MaskedCopy
1972 {
1973 using TreeTypePtr = typename TreeType::Ptr;
1974 using ValueType = typename TreeType::ValueType;
1975 using LeafNodeType = typename TreeType::LeafNodeType;
1976
1977 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1978 using BoolTreeTypePtr = typename BoolTreeType::Ptr;
1979 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1980
MaskedCopyMaskedCopy1981 MaskedCopy(const TreeType& tree, std::vector<TreeTypePtr>& segments,
1982 std::vector<BoolTreeTypePtr>& masks)
1983 : mTree(&tree)
1984 , mSegments(!segments.empty() ? &segments.front() : nullptr)
1985 , mMasks(!masks.empty() ? &masks.front() : nullptr)
1986 {
1987 }
1988
operatorMaskedCopy1989 void operator()(const tbb::blocked_range<size_t>& range) const {
1990
1991 std::vector<const BoolLeafNodeType*> nodes;
1992
1993 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1994
1995 const BoolTreeType& mask = *mMasks[n];
1996
1997 nodes.clear();
1998 mask.getNodes(nodes);
1999
2000 Copy op(*mTree, nodes);
2001 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2002 mSegments[n] = op.outputTree();
2003 }
2004 }
2005
2006 private:
2007
2008 struct Copy {
CopyMaskedCopy::Copy2009 Copy(const TreeType& inputTree, std::vector<const BoolLeafNodeType*>& maskNodes)
2010 : mInputTree(&inputTree)
2011 , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
2012 , mOutputTreePtr(new TreeType(inputTree.background()))
2013 {
2014 }
2015
CopyMaskedCopy::Copy2016 Copy(const Copy& rhs, tbb::split)
2017 : mInputTree(rhs.mInputTree)
2018 , mMaskNodes(rhs.mMaskNodes)
2019 , mOutputTreePtr(new TreeType(mInputTree->background()))
2020 {
2021 }
2022
outputTreeMaskedCopy::Copy2023 TreeTypePtr& outputTree() { return mOutputTreePtr; }
2024
joinMaskedCopy::Copy2025 void join(Copy& rhs) { mOutputTreePtr->merge(*rhs.mOutputTreePtr); }
2026
operatorMaskedCopy::Copy2027 void operator()(const tbb::blocked_range<size_t>& range) {
2028
2029 tree::ValueAccessor<const TreeType> inputAcc(*mInputTree);
2030 tree::ValueAccessor<TreeType> outputAcc(*mOutputTreePtr);
2031
2032 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2033
2034 const BoolLeafNodeType& maskNode = *mMaskNodes[n];
2035 if (maskNode.isEmpty()) continue;
2036
2037 const Coord& ijk = maskNode.origin();
2038
2039 const LeafNodeType* inputNode = inputAcc.probeConstLeaf(ijk);
2040 if (inputNode) {
2041
2042 LeafNodeType* outputNode = outputAcc.touchLeaf(ijk);
2043
2044 for (typename BoolLeafNodeType::ValueOnCIter it = maskNode.cbeginValueOn();
2045 it; ++it)
2046 {
2047 const Index idx = it.pos();
2048 outputNode->setValueOn(idx, inputNode->getValue(idx));
2049 }
2050 } else {
2051 const int valueDepth = inputAcc.getValueDepth(ijk);
2052 if (valueDepth >= 0) {
2053 outputAcc.addTile(TreeType::RootNodeType::LEVEL - valueDepth,
2054 ijk, inputAcc.getValue(ijk), true);
2055 }
2056 }
2057 }
2058 }
2059
2060 private:
2061 TreeType const * const mInputTree;
2062 BoolLeafNodeType const * const * const mMaskNodes;
2063 TreeTypePtr mOutputTreePtr;
2064 }; // struct Copy
2065
2066 TreeType const * const mTree;
2067 TreeTypePtr * const mSegments;
2068 BoolTreeTypePtr * const mMasks;
2069 }; // MaskedCopy
2070
2071
2072 ////////////////////////////////////////
2073
2074
2075 template<typename VolumePtrType>
2076 struct ComputeActiveVoxelCount
2077 {
ComputeActiveVoxelCountComputeActiveVoxelCount2078 ComputeActiveVoxelCount(std::vector<VolumePtrType>& segments, size_t *countArray)
2079 : mSegments(!segments.empty() ? &segments.front() : nullptr)
2080 , mCountArray(countArray)
2081 {
2082 }
2083
operatorComputeActiveVoxelCount2084 void operator()(const tbb::blocked_range<size_t>& range) const {
2085 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2086 mCountArray[n] = mSegments[n]->activeVoxelCount();
2087 }
2088 }
2089
2090 VolumePtrType * const mSegments;
2091 size_t * const mCountArray;
2092 };
2093
2094
2095 struct GreaterCount
2096 {
GreaterCountGreaterCount2097 GreaterCount(const size_t *countArray) : mCountArray(countArray) {}
2098
operatorGreaterCount2099 inline bool operator() (const size_t& lhs, const size_t& rhs) const
2100 {
2101 return (mCountArray[lhs] > mCountArray[rhs]);
2102 }
2103
2104 size_t const * const mCountArray;
2105 };
2106
2107 ////////////////////////////////////////
2108
2109
2110 template<typename TreeType>
2111 struct GridOrTreeConstructor
2112 {
2113 using TreeTypePtr = typename TreeType::Ptr;
2114 using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2115
constructMaskGridOrTreeConstructor2116 static BoolTreePtrType constructMask(const TreeType&, BoolTreePtrType& maskTree)
2117 { return maskTree; }
constructGridOrTreeConstructor2118 static TreeTypePtr construct(const TreeType&, TreeTypePtr& tree) { return tree; }
2119 };
2120
2121
2122 template<typename TreeType>
2123 struct GridOrTreeConstructor<Grid<TreeType> >
2124 {
2125 using GridType = Grid<TreeType>;
2126 using GridTypePtr = typename Grid<TreeType>::Ptr;
2127 using TreeTypePtr = typename TreeType::Ptr;
2128
2129 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2130 using BoolTreePtrType = typename BoolTreeType::Ptr;
2131 using BoolGridType = Grid<BoolTreeType>;
2132 using BoolGridPtrType = typename BoolGridType::Ptr;
2133
2134 static BoolGridPtrType constructMask(const GridType& grid, BoolTreePtrType& maskTree) {
2135 BoolGridPtrType maskGrid(BoolGridType::create(maskTree));
2136 maskGrid->setTransform(grid.transform().copy());
2137 return maskGrid;
2138 }
2139
2140 static GridTypePtr construct(const GridType& grid, TreeTypePtr& maskTree) {
2141 GridTypePtr maskGrid(GridType::create(maskTree));
2142 maskGrid->setTransform(grid.transform().copy());
2143 maskGrid->insertMeta(grid);
2144 return maskGrid;
2145 }
2146 };
2147
2148
2149 } // namespace level_set_util_internal
2150
2151
2152 /// @endcond OPENVDB_DOCS_INTERNAL
2153
2154 ////////////////////////////////////////
2155
2156
2157 template <class GridType>
2158 void
2159 sdfToFogVolume(GridType& grid, typename GridType::ValueType cutoffDistance)
2160 {
2161 using ValueType = typename GridType::ValueType;
2162 using TreeType = typename GridType::TreeType;
2163 using LeafNodeType = typename TreeType::LeafNodeType;
2164 using RootNodeType = typename TreeType::RootNodeType;
2165 using NodeChainType = typename RootNodeType::NodeChainType;
2166 using InternalNodeType = typename NodeChainType::template Get<1>;
2167
2168 //////////
2169
2170 TreeType& tree = grid.tree();
2171
2172 size_t numLeafNodes = 0, numInternalNodes = 0;
2173
2174 std::vector<LeafNodeType*> nodes;
2175 std::vector<size_t> leafnodeCount;
2176
2177 {
2178 // Compute the prefix sum of the leafnode count in each internal node.
2179 std::vector<InternalNodeType*> internalNodes;
2180 tree.getNodes(internalNodes);
2181
2182 numInternalNodes = internalNodes.size();
2183
2184 leafnodeCount.push_back(0);
2185 for (size_t n = 0; n < numInternalNodes; ++n) {
2186 leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
2187 }
2188
2189 numLeafNodes = leafnodeCount.back();
2190
2191 // Steal all leafnodes (Removes them from the tree and transfers ownership.)
2192 nodes.reserve(numLeafNodes);
2193
2194 for (size_t n = 0; n < numInternalNodes; ++n) {
2195 internalNodes[n]->stealNodes(nodes, tree.background(), false);
2196 }
2197
2198 // Clamp cutoffDistance to min sdf value
2199 ValueType minSDFValue = std::numeric_limits<ValueType>::max();
2200
2201 {
2202 level_set_util_internal::FindMinTileValue<InternalNodeType> minOp(internalNodes.data());
2203 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, internalNodes.size()), minOp);
2204 minSDFValue = std::min(minSDFValue, minOp.minValue);
2205 }
2206
2207 if (minSDFValue > ValueType(0.0)) {
2208 level_set_util_internal::FindMinVoxelValue<LeafNodeType> minOp(nodes.data());
2209 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
2210 minSDFValue = std::min(minSDFValue, minOp.minValue);
2211 }
2212
2213 cutoffDistance = -std::abs(cutoffDistance);
2214 cutoffDistance = minSDFValue > cutoffDistance ? minSDFValue : cutoffDistance;
2215 }
2216
2217 // Transform voxel values and delete leafnodes that are uniformly zero after the transformation.
2218 // (Positive values are set to zero with inactive state and negative values are remapped
2219 // from zero to one with active state.)
2220 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
2221 level_set_util_internal::SDFVoxelsToFogVolume<LeafNodeType>(nodes.data(), cutoffDistance));
2222
2223 // Populate a new tree with the remaining leafnodes
2224 typename TreeType::Ptr newTree(new TreeType(ValueType(0.0)));
2225
2226 level_set_util_internal::PopulateTree<TreeType> populate(
2227 *newTree, nodes.data(), leafnodeCount.data(), 0);
2228 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
2229
2230 // Transform tile values (Negative valued tiles are set to 1.0 with active state.)
2231 std::vector<InternalNodeType*> internalNodes;
2232 newTree->getNodes(internalNodes);
2233
2234 tbb::parallel_for(tbb::blocked_range<size_t>(0, internalNodes.size()),
2235 level_set_util_internal::SDFTilesToFogVolume<TreeType, InternalNodeType>(
2236 tree, internalNodes.data()));
2237
2238 {
2239 tree::ValueAccessor<const TreeType> acc(tree);
2240
2241 typename TreeType::ValueAllIter it(*newTree);
2242 it.setMaxDepth(TreeType::ValueAllIter::LEAF_DEPTH - 2);
2243
2244 for ( ; it; ++it) {
2245 if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
2246 it.setValue(ValueType(1.0));
2247 it.setActiveState(true);
2248 }
2249 }
2250 }
2251
2252 // Insert missing root level tiles. (The new tree is constructed from the remaining leafnodes
2253 // and will therefore not contain any root level tiles that may exist in the original tree.)
2254 {
2255 typename TreeType::ValueAllIter it(tree);
2256 it.setMaxDepth(TreeType::ValueAllIter::ROOT_DEPTH);
2257 for ( ; it; ++it) {
2258 if (it.getValue() < ValueType(0.0)) {
2259 newTree->addTile(TreeType::ValueAllIter::ROOT_LEVEL, it.getCoord(),
2260 ValueType(1.0), true);
2261 }
2262 }
2263 }
2264
2265 grid.setTree(newTree);
2266 grid.setGridClass(GRID_FOG_VOLUME);
2267 }
2268
2269
2270 ////////////////////////////////////////
2271
2272
2273 template <class GridOrTreeType>
2274 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2275 sdfInteriorMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2276 {
2277 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2278 const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2279
2280 using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2281 BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(tree, isovalue);
2282
2283 return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2284 volume, mask);
2285 }
2286
2287
2288 template<typename GridOrTreeType>
2289 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2290 extractEnclosedRegion(const GridOrTreeType& volume,
2291 typename GridOrTreeType::ValueType isovalue,
2292 const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
2293 fillMask)
2294 {
2295 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2296 const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2297
2298 using CharTreePtrType = typename TreeType::template ValueConverter<char>::Type::Ptr;
2299 CharTreePtrType regionMask = level_set_util_internal::computeEnclosedRegionMask(
2300 tree, isovalue, fillMask);
2301
2302 using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2303 BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(*regionMask, 0);
2304
2305 return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2306 volume, mask);
2307 }
2308
2309
2310 ////////////////////////////////////////
2311
2312
2313 template<typename GridOrTreeType>
2314 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2315 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2316 {
2317 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2318 const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2319
2320 std::vector<const typename TreeType::LeafNodeType*> nodes;
2321 tree.getNodes(nodes);
2322
2323 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2324 typename BoolTreeType::Ptr mask(new BoolTreeType(false));
2325
2326 level_set_util_internal::MaskIsovalueCrossingVoxels<TreeType> op(tree, nodes, *mask, isovalue);
2327 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2328
2329 return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2330 volume, mask);
2331 }
2332
2333
2334 ////////////////////////////////////////
2335
2336
2337 template<typename GridOrTreeType>
2338 void
2339 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
2340 std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks)
2341 {
2342 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2343 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2344 using BoolTreePtrType = typename BoolTreeType::Ptr;
2345 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2346 using BoolGridOrTreePtrType = typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr;
2347
2348 using NodeMaskSegmentType = level_set_util_internal::NodeMaskSegment<BoolLeafNodeType>;
2349 using NodeMaskSegmentPtrType = typename NodeMaskSegmentType::Ptr;
2350 using NodeMaskSegmentPtrVector = typename std::vector<NodeMaskSegmentPtrType>;
2351 using NodeMaskSegmentRawPtrVector = typename std::vector<NodeMaskSegmentType*>;
2352
2353 /////
2354
2355 const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2356
2357 BoolTreeType topologyMask(tree, false, TopologyCopy());
2358
2359 // prune out any inactive leaf nodes or inactive tiles
2360 tools::pruneInactive(topologyMask);
2361
2362 if (topologyMask.hasActiveTiles()) {
2363 topologyMask.voxelizeActiveTiles();
2364 }
2365
2366 std::vector<BoolLeafNodeType*> leafnodes;
2367 topologyMask.getNodes(leafnodes);
2368
2369 if (leafnodes.empty()) return;
2370
2371 // 1. Split node masks into disjoint segments
2372 // Note: The LeafNode origin coord is modified to record the 'leafnodes' array offset.
2373
2374 std::unique_ptr<NodeMaskSegmentPtrVector[]> nodeSegmentArray(
2375 new NodeMaskSegmentPtrVector[leafnodes.size()]);
2376
2377 tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2378 level_set_util_internal::SegmentNodeMask<BoolLeafNodeType>(
2379 leafnodes, nodeSegmentArray.get()));
2380
2381
2382 // 2. Compute segment connectivity
2383
2384 tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2385 level_set_util_internal::ConnectNodeMaskSegments<BoolTreeType, BoolLeafNodeType>(
2386 topologyMask, nodeSegmentArray.get()));
2387
2388 topologyMask.clear();
2389
2390 size_t nodeSegmentCount = 0;
2391 for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2392 nodeSegmentCount += nodeSegmentArray[n].size();
2393 }
2394
2395 // 3. Group connected segments
2396
2397 std::deque<NodeMaskSegmentRawPtrVector> nodeSegmentGroups;
2398
2399 NodeMaskSegmentType* nextSegment = nodeSegmentArray[0][0].get();
2400 while (nextSegment) {
2401
2402 nodeSegmentGroups.push_back(NodeMaskSegmentRawPtrVector());
2403
2404 std::vector<NodeMaskSegmentType*>& segmentGroup = nodeSegmentGroups.back();
2405 segmentGroup.reserve(nodeSegmentCount);
2406
2407 std::deque<NodeMaskSegmentType*> segmentQueue;
2408 segmentQueue.push_back(nextSegment);
2409 nextSegment = nullptr;
2410
2411 while (!segmentQueue.empty()) {
2412
2413 NodeMaskSegmentType* segment = segmentQueue.back();
2414 segmentQueue.pop_back();
2415
2416 if (segment->visited) continue;
2417 segment->visited = true;
2418
2419 segmentGroup.push_back(segment);
2420
2421 // queue connected segments
2422 std::vector<NodeMaskSegmentType*>& connections = segment->connections;
2423 for (size_t n = 0, N = connections.size(); n < N; ++n) {
2424 if (!connections[n]->visited) segmentQueue.push_back(connections[n]);
2425 }
2426 }
2427
2428 // find first unvisited segment
2429 for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2430 NodeMaskSegmentPtrVector& nodeSegments = nodeSegmentArray[n];
2431 for (size_t i = 0, I = nodeSegments.size(); i < I; ++i) {
2432 if (!nodeSegments[i]->visited) nextSegment = nodeSegments[i].get();
2433 }
2434 }
2435 }
2436
2437 // 4. Mask segment groups
2438
2439 if (nodeSegmentGroups.size() == 1) {
2440
2441 BoolTreePtrType mask(new BoolTreeType(tree, false, TopologyCopy()));
2442
2443 tools::pruneInactive(*mask);
2444
2445 if (mask->hasActiveTiles()) {
2446 mask->voxelizeActiveTiles();
2447 }
2448
2449 masks.push_back(
2450 level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2451 volume, mask));
2452
2453 } else if (nodeSegmentGroups.size() > 1) {
2454
2455 for (size_t n = 0, N = nodeSegmentGroups.size(); n < N; ++n) {
2456
2457 NodeMaskSegmentRawPtrVector& segmentGroup = nodeSegmentGroups[n];
2458
2459 level_set_util_internal::MaskSegmentGroup<BoolTreeType> op(segmentGroup);
2460 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, segmentGroup.size()), op);
2461
2462 masks.push_back(
2463 level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2464 volume, op.mask()));
2465 }
2466 }
2467
2468 // 5. Sort segments in descending order based on the active voxel count.
2469
2470 if (masks.size() > 1) {
2471 const size_t segmentCount = masks.size();
2472
2473 std::unique_ptr<size_t[]> segmentOrderArray(new size_t[segmentCount]);
2474 std::unique_ptr<size_t[]> voxelCountArray(new size_t[segmentCount]);
2475
2476 for (size_t n = 0; n < segmentCount; ++n) {
2477 segmentOrderArray[n] = n;
2478 }
2479
2480 tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
2481 level_set_util_internal::ComputeActiveVoxelCount<BoolGridOrTreePtrType>(
2482 masks, voxelCountArray.get()));
2483
2484 size_t *begin = segmentOrderArray.get();
2485 tbb::parallel_sort(begin, begin + masks.size(), level_set_util_internal::GreaterCount(
2486 voxelCountArray.get()));
2487
2488 std::vector<BoolGridOrTreePtrType> orderedMasks;
2489 orderedMasks.reserve(masks.size());
2490
2491 for (size_t n = 0; n < segmentCount; ++n) {
2492 orderedMasks.push_back(masks[segmentOrderArray[n]]);
2493 }
2494
2495 masks.swap(orderedMasks);
2496 }
2497
2498 } // extractActiveVoxelSegmentMasks()
2499
2500
2501 template<typename GridOrTreeType>
2502 void
2503 segmentActiveVoxels(const GridOrTreeType& volume,
2504 std::vector<typename GridOrTreeType::Ptr>& segments)
2505 {
2506 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2507 using TreePtrType = typename TreeType::Ptr;
2508 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2509 using BoolTreePtrType = typename BoolTreeType::Ptr;
2510
2511 const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2512
2513 // 1. Segment active topology mask
2514 std::vector<BoolTreePtrType> maskSegmentArray;
2515 extractActiveVoxelSegmentMasks(inputTree, maskSegmentArray);
2516
2517 // 2. Export segments
2518
2519 const size_t numSegments = std::max(size_t(1), maskSegmentArray.size());
2520 std::vector<TreePtrType> outputSegmentArray(numSegments);
2521
2522 if (maskSegmentArray.empty()) {
2523 // if no active voxels in the original volume, copy just the background
2524 // value of the input tree
2525 outputSegmentArray[0] = TreePtrType(new TreeType(inputTree.background()));
2526 } else if (numSegments == 1) {
2527 // if there's only one segment with active voxels, copy the input tree
2528 TreePtrType segment(new TreeType(inputTree));
2529 // however, if the leaf counts do not match due to the pruning of inactive leaf
2530 // nodes in the mask, do a topology intersection to drop these inactive leafs
2531 if (segment->leafCount() != inputTree.leafCount()) {
2532 segment->topologyIntersection(*maskSegmentArray[0]);
2533 }
2534 outputSegmentArray[0] = segment;
2535 } else {
2536 const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2537 tbb::parallel_for(segmentRange,
2538 level_set_util_internal::MaskedCopy<TreeType>(inputTree, outputSegmentArray,
2539 maskSegmentArray));
2540 }
2541
2542 for (auto& segment : outputSegmentArray) {
2543 segments.push_back(
2544 level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::construct(
2545 volume, segment));
2546 }
2547 }
2548
2549
2550 template<typename GridOrTreeType>
2551 void
2552 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments)
2553 {
2554 using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2555 using TreePtrType = typename TreeType::Ptr;
2556 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2557 using BoolTreePtrType = typename BoolTreeType::Ptr;
2558
2559 const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2560
2561 // 1. Mask zero crossing voxels
2562 BoolTreePtrType mask = extractIsosurfaceMask(inputTree, lsutilGridZero<GridOrTreeType>());
2563
2564 // 2. Segment the zero crossing mask
2565 std::vector<BoolTreePtrType> maskSegmentArray;
2566 extractActiveVoxelSegmentMasks(*mask, maskSegmentArray);
2567
2568 const size_t numSegments = std::max(size_t(1), maskSegmentArray.size());
2569 std::vector<TreePtrType> outputSegmentArray(numSegments);
2570
2571 if (maskSegmentArray.empty()) {
2572 // if no active voxels in the original volume, copy just the background
2573 // value of the input tree
2574 outputSegmentArray[0] = TreePtrType(new TreeType(inputTree.background()));
2575 } else {
2576 const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2577
2578 // 3. Expand zero crossing mask to capture sdf narrow band
2579 tbb::parallel_for(segmentRange,
2580 level_set_util_internal::ExpandNarrowbandMask<TreeType>(inputTree, maskSegmentArray));
2581
2582 // 4. Export sdf segments
2583
2584 tbb::parallel_for(segmentRange, level_set_util_internal::MaskedCopy<TreeType>(
2585 inputTree, outputSegmentArray, maskSegmentArray));
2586
2587 tbb::parallel_for(segmentRange,
2588 level_set_util_internal::FloodFillSign<TreeType>(inputTree, outputSegmentArray));
2589 }
2590
2591 for (auto& segment : outputSegmentArray) {
2592 segments.push_back(
2593 level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::construct(
2594 volume, segment));
2595 }
2596 }
2597
2598
2599 ////////////////////////////////////////
2600
2601
2602 // Explicit Template Instantiation
2603
2604 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
2605
2606 #ifdef OPENVDB_INSTANTIATE_LEVELSETUTIL
2607 #include <openvdb/util/ExplicitInstantiation.h>
2608 #endif
2609
2610 #define _FUNCTION(TreeT) \
2611 void sdfToFogVolume(Grid<TreeT>&, TreeT::ValueType)
2612 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2613 #undef _FUNCTION
2614
2615 #define _FUNCTION(TreeT) \
2616 TreeT::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const TreeT&, TreeT::ValueType)
2617 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2618 #undef _FUNCTION
2619
2620 #define _FUNCTION(TreeT) \
2621 Grid<TreeT>::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const Grid<TreeT>&, TreeT::ValueType)
2622 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2623 #undef _FUNCTION
2624
2625 #define _FUNCTION(TreeT) \
2626 TreeT::ValueConverter<bool>::Type::Ptr extractEnclosedRegion(\
2627 const TreeT&, TreeT::ValueType, \
2628 const TreeAdapter<TreeT>::TreeType::ValueConverter<bool>::Type*)
2629 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2630 #undef _FUNCTION
2631
2632 #define _FUNCTION(TreeT) \
2633 Grid<TreeT>::ValueConverter<bool>::Type::Ptr extractEnclosedRegion(\
2634 const Grid<TreeT>&, TreeT::ValueType, \
2635 const TreeAdapter<Grid<TreeT>>::TreeType::ValueConverter<bool>::Type*)
2636 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2637 #undef _FUNCTION
2638
2639 #define _FUNCTION(TreeT) \
2640 TreeT::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const TreeT&, TreeT::ValueType)
2641 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2642 #undef _FUNCTION
2643
2644 #define _FUNCTION(TreeT) \
2645 Grid<TreeT>::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const Grid<TreeT>&, TreeT::ValueType)
2646 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2647 #undef _FUNCTION
2648
2649 #define _FUNCTION(TreeT) \
2650 void extractActiveVoxelSegmentMasks(\
2651 const TreeT&, std::vector<TreeT::ValueConverter<bool>::Type::Ptr>&)
2652 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2653 #undef _FUNCTION
2654
2655 #define _FUNCTION(TreeT) \
2656 void extractActiveVoxelSegmentMasks(\
2657 const Grid<TreeT>&, std::vector<Grid<TreeT>::ValueConverter<bool>::Type::Ptr>&)
2658 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2659 #undef _FUNCTION
2660
2661 #define _FUNCTION(TreeT) \
2662 void segmentActiveVoxels(const TreeT&, std::vector<TreeT::Ptr>&)
2663 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2664 #undef _FUNCTION
2665
2666 #define _FUNCTION(TreeT) \
2667 void segmentActiveVoxels(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
2668 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2669 #undef _FUNCTION
2670
2671 #define _FUNCTION(TreeT) \
2672 void segmentSDF(const TreeT&, std::vector<TreeT::Ptr>&)
2673 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2674 #undef _FUNCTION
2675
2676 #define _FUNCTION(TreeT) \
2677 void segmentSDF(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
2678 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2679 #undef _FUNCTION
2680
2681 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
2682
2683
2684 } // namespace tools
2685 } // namespace OPENVDB_VERSION_NAME
2686 } // namespace openvdb
2687
2688 #endif // OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
2689