OpenVDB  10.0.0
LevelSetUtil.h
Go to the documentation of this file.
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>
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 {
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>
45 inline typename GridType::ValueType lsutilGridMax()
46 {
48 }
49 
50 template<typename GridType>
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
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
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
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 
192  MaskInteriorVoxels(
193  ValueType isovalue, const LeafNodeType ** nodes, BoolLeafNodeType ** maskNodes)
194  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
195  {
196  }
197 
198  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  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 
238  MaskInteriorTiles(ValueType isovalue, const TreeType& tree, InternalNodeType ** maskNodes)
239  : mTree(&tree), mMaskNodes(maskNodes), mIsovalue(isovalue) { }
240 
241  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 
266  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 
275  PopulateTree(PopulateTree& rhs, tbb::split)
276  : mNewTree(rhs.mNewTree.background())
277  , mTreePt(&mNewTree)
278  , mNodes(rhs.mNodes)
279  , mNodeIndexMap(rhs.mNodeIndexMap)
280  {
281  }
282 
283  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 
300  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 
317  LabelBoundaryVoxels(
318  ValueType isovalue, const LeafNodeType ** nodes, CharLeafNodeType ** maskNodes)
319  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
320  {
321  }
322 
323  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  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 
362  FlipRegionSign(LeafNodeType ** nodes) : mNodes(nodes) { }
363 
364  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 
382  FindMinVoxelValue(LeafNodeType const * const * const leafnodes)
383  : minValue(std::numeric_limits<ValueType>::max())
384  , mNodes(leafnodes)
385  {
386  }
387 
388  FindMinVoxelValue(FindMinVoxelValue& rhs, tbb::split)
389  : minValue(std::numeric_limits<ValueType>::max())
390  , mNodes(rhs.mNodes)
391  {
392  }
393 
394  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 
403  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 
416  FindMinTileValue(InternalNodeType const * const * const nodes)
417  : minValue(std::numeric_limits<ValueType>::max())
418  , mNodes(nodes)
419  {
420  }
421 
422  FindMinTileValue(FindMinTileValue& rhs, tbb::split)
423  : minValue(std::numeric_limits<ValueType>::max())
424  , mNodes(rhs.mNodes)
425  {
426  }
427 
428  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 
437  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 
450  SDFVoxelsToFogVolume(LeafNodeType ** nodes, ValueType cutoffDistance)
451  : mNodes(nodes), mWeight(ValueType(1.0) / cutoffDistance)
452  {
453  }
454 
455  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 
483  SDFTilesToFogVolume(const TreeType& tree, InternalNodeType ** nodes)
484  : mTree(&tree), mNodes(nodes) { }
485 
486  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 
515  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 
525  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.
617  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 
673  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 
726  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 
783  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>
839  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>
874  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>
908  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
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
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 
1167  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 
1180  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 
1189  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  delete maskNodePt;
1258  }
1259 
1260  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 
1284  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
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 
1376  SegmentNodeMask(std::vector<NodeType*>& nodes, NodeMaskSegmentVector* nodeMaskArray)
1377  : mNodes(!nodes.empty() ? &nodes.front() : nullptr)
1378  , mNodeMaskArray(nodeMaskArray)
1379  {
1380  }
1381 
1382  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 
1406  ConnectNodeMaskSegments(const TreeType& tree, NodeMaskSegmentVector* nodeMaskArray)
1407  : mTree(&tree)
1408  , mNodeMaskArray(nodeMaskArray)
1409  {
1410  }
1411 
1412  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 
1537  static inline size_t getNodeOffset(const NodeType& node) {
1538  return static_cast<size_t>(node.origin()[0]);
1539  }
1540 
1541  static inline NodeMaskSegmentType*
1542  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
1557  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 
1577  MaskSegmentGroup(const std::vector<NodeMaskSegmentType*>& segments)
1578  : mSegments(!segments.empty() ? &segments.front() : nullptr)
1579  , mTree(new TreeType(false))
1580  {
1581  }
1582 
1583  MaskSegmentGroup(const MaskSegmentGroup& rhs, tbb::split)
1584  : mSegments(rhs.mSegments)
1585  , mTree(new TreeType(false))
1586  {
1587  }
1588 
1589  TreeTypePtr& mask() { return mTree; }
1590 
1591  void join(MaskSegmentGroup& rhs) { mTree->merge(*rhs.mTree); }
1592 
1593  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 
1625  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 
1634  ExpandLeafNodeRegion(const ExpandLeafNodeRegion& rhs, tbb::split)
1635  : mDistTree(rhs.mDistTree)
1636  , mMaskTree(rhs.mMaskTree)
1637  , mMaskNodes(rhs.mMaskNodes)
1638  , mNewMaskTree(false)
1639  {
1640  }
1641 
1642  BoolTreeType& newMaskTree() { return mNewMaskTree; }
1643 
1644  void join(ExpandLeafNodeRegion& rhs) { mNewMaskTree.merge(rhs.mNewMaskTree); }
1645 
1646  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
1739  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 
1768  FillLeafNodeVoxels(const TreeType& tree, std::vector<BoolLeafNodeType*>& maskNodes)
1769  : mTree(&tree), mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
1770  {
1771  }
1772 
1773  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 
1869  ExpandNarrowbandMask(const TreeType& tree, std::vector<BoolTreeTypePtr>& segments)
1870  : mTree(&tree), mSegments(!segments.empty() ? &segments.front() : nullptr)
1871  {
1872  }
1873 
1874  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 
1923  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 
1954  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 
1981  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 
1989  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 {
2009  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 
2016  Copy(const Copy& rhs, tbb::split)
2017  : mInputTree(rhs.mInputTree)
2018  , mMaskNodes(rhs.mMaskNodes)
2019  , mOutputTreePtr(new TreeType(mInputTree->background()))
2020  {
2021  }
2022 
2023  TreeTypePtr& outputTree() { return mOutputTreePtr; }
2024 
2025  void join(Copy& rhs) { mOutputTreePtr->merge(*rhs.mOutputTreePtr); }
2026 
2027  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 {
2078  ComputeActiveVoxelCount(std::vector<VolumePtrType>& segments, size_t *countArray)
2079  : mSegments(!segments.empty() ? &segments.front() : nullptr)
2080  , mCountArray(countArray)
2081  {
2082  }
2083 
2084  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 {
2097  GreaterCount(const size_t *countArray) : mCountArray(countArray) {}
2098 
2099  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 
2116  static BoolTreePtrType constructMask(const TreeType&, BoolTreePtrType& maskTree)
2117  { return maskTree; }
2118  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  {
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
2608 #endif
2609 
2610 #define _FUNCTION(TreeT) \
2611  void sdfToFogVolume(Grid<TreeT>&, TreeT::ValueType)
2613 #undef _FUNCTION
2614 
2615 #define _FUNCTION(TreeT) \
2616  TreeT::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const TreeT&, TreeT::ValueType)
2618 #undef _FUNCTION
2619 
2620 #define _FUNCTION(TreeT) \
2621  Grid<TreeT>::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const Grid<TreeT>&, TreeT::ValueType)
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*)
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*)
2637 #undef _FUNCTION
2638 
2639 #define _FUNCTION(TreeT) \
2640  TreeT::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const TreeT&, TreeT::ValueType)
2642 #undef _FUNCTION
2643 
2644 #define _FUNCTION(TreeT) \
2645  Grid<TreeT>::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const Grid<TreeT>&, TreeT::ValueType)
2647 #undef _FUNCTION
2648 
2649 #define _FUNCTION(TreeT) \
2650  void extractActiveVoxelSegmentMasks(\
2651  const TreeT&, std::vector<TreeT::ValueConverter<bool>::Type::Ptr>&)
2653 #undef _FUNCTION
2654 
2655 #define _FUNCTION(TreeT) \
2656  void extractActiveVoxelSegmentMasks(\
2657  const Grid<TreeT>&, std::vector<Grid<TreeT>::ValueConverter<bool>::Type::Ptr>&)
2659 #undef _FUNCTION
2660 
2661 #define _FUNCTION(TreeT) \
2662  void segmentActiveVoxels(const TreeT&, std::vector<TreeT::Ptr>&)
2664 #undef _FUNCTION
2665 
2666 #define _FUNCTION(TreeT) \
2667  void segmentActiveVoxels(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
2669 #undef _FUNCTION
2670 
2671 #define _FUNCTION(TreeT) \
2672  void segmentSDF(const TreeT&, std::vector<TreeT::Ptr>&)
2674 #undef _FUNCTION
2675 
2676 #define _FUNCTION(TreeT) \
2677  void segmentSDF(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
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
#define ROOT_LEVEL
Definition: CNanoVDB.h:53
ValueT value
Definition: GridBuilder.h:1290
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
Attribute-owned data structure for points. Point attributes are stored in leaf nodes and ordered by v...
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
Definition: ValueAccessor.h:191
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:235
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:243
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:253
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractEnclosedRegion(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >(), const typename TreeAdapter< GridOrTreeType >::TreeType::template ValueConverter< bool >::Type *fillMask=nullptr)
Extracts the interior regions of a signed distance field and topologically enclosed (watertight) regi...
Definition: LevelSetUtil.h:2290
void segmentActiveVoxels(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint active topology components into distinct grids or trees.
Definition: LevelSetUtil.h:2503
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractIsosurfaceMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue)
Return a mask of the voxels that intersect the implicit surface with the given isovalue.
Definition: LevelSetUtil.h:2315
void extractActiveVoxelSegmentMasks(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::template ValueConverter< bool >::Type::Ptr > &masks)
Return a mask for each connected component of the given grid's active voxels.
Definition: LevelSetUtil.h:2339
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2275
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition: MeshToVolume.h:3025
OPENVDB_DOCS_INTERNAL void sdfToFogVolume(GridType &grid, typename GridType::ValueType cutoffDistance)
Threaded method to convert a sparse level set/SDF into a sparse fog volume.
Definition: LevelSetUtil.h:2159
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:267
void segmentSDF(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint SDF surfaces into distinct grids or trees.
Definition: LevelSetUtil.h:2552
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:355
Index32 Index
Definition: Types.h:54
@ GRID_FOG_VOLUME
Definition: Types.h:417
openvdb::GridBase Grid
Definition: Utils.h:34
Definition: Exceptions.h:13
Definition: Coord.h:587
This adapter allows code that is templated on a Tree type to accept either a Tree type or a Grid type...
Definition: Grid.h:1060
static TreeType & tree(TreeType &t)
Definition: Grid.h:1076
_TreeType TreeType
Definition: Grid.h:1061
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157
#define OPENVDB_ALL_TREE_INSTANTIATE(Function)
Definition: version.h.in:161