OpenVDB  13.0.0
LevelSetUtil.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 
4 /// @file tools/LevelSetUtil.h
5 ///
6 /// @brief Miscellaneous utility methods that operate primarily
7 /// or exclusively on level set grids.
8 ///
9 /// @author Mihai Alden
10 
11 #ifndef OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
12 #define OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
13 
14 #include "MeshToVolume.h" // for traceExteriorBoundaries
15 #include "SignedFloodFill.h" // for signedFloodFillWithValues
16 
17 #include <openvdb/Types.h>
18 #include <openvdb/Grid.h>
19 #include <openvdb/openvdb.h>
20 #include <openvdb/points/PointDataGrid.h>
21 #include <tbb/blocked_range.h>
22 #include <tbb/parallel_for.h>
23 #include <tbb/parallel_reduce.h>
24 #include <tbb/parallel_sort.h>
25 #include <algorithm>
26 #include <cmath>
27 #include <cstdlib>
28 #include <deque>
29 #include <limits>
30 #include <memory>
31 #include <set>
32 #include <vector>
33 
34 
35 namespace openvdb {
36 OPENVDB_USE_VERSION_NAMESPACE
37 namespace OPENVDB_VERSION_NAME {
38 namespace tools {
39 
40 // MS Visual C++ requires this extra level of indirection in order to compile
41 // THIS MUST EXIST IN AN UNNAMED NAMESPACE IN ORDER TO COMPILE ON WINDOWS
42 namespace {
43 
44 template<typename GridType>
45 inline typename GridType::ValueType lsutilGridMax()
46 {
47  return std::numeric_limits<typename GridType::ValueType>::max();
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
78 sdfToFogVolume(
79  GridType& grid,
80  typename GridType::ValueType cutoffDistance = lsutilGridMax<GridType>());
81 
82 
83 /// @brief Threaded method to construct a boolean mask that represents interior regions
84 /// in a signed distance field.
85 ///
86 /// @return A shared pointer to either a boolean grid or tree with the same tree
87 /// configuration and potentially transform as the input @c volume and whose active
88 /// and @c true values correspond to the interior of the input signed distance field.
89 ///
90 /// @param volume Signed distance field / level set volume.
91 /// @param isovalue Threshold below which values are considered part of the
92 /// interior region.
93 template<class GridOrTreeType>
94 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
95 sdfInteriorMask(
96  const GridOrTreeType& volume,
97  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>());
98 
99 
100 /// @brief Extracts the interior regions of a signed distance field and topologically enclosed
101 /// (watertight) regions of value greater than the @a isovalue (cavities) that can arise
102 /// as the result of CSG union operations between different shapes where at least one of
103 /// the shapes has a concavity that is capped.
104 ///
105 /// For example the enclosed region of a capped bottle would include the walls and
106 /// the interior cavity.
107 ///
108 /// @return A shared pointer to either a boolean grid or tree with the same tree configuration
109 /// and potentially transform as the input @c volume and whose active and @c true values
110 /// correspond to the interior and enclosed regions in the input signed distance field.
111 ///
112 /// @param volume Signed distance field / level set volume.
113 /// @param isovalue Threshold below which values are considered part of the interior region.
114 /// @param fillMask Optional boolean tree, when provided enclosed cavity regions that are not
115 /// completely filled by this mask are ignored.
116 ///
117 /// For instance if the fill mask does not completely fill the bottle in the
118 /// previous example only the walls and cap are returned and the interior
119 /// cavity will be ignored.
120 template<typename GridOrTreeType>
121 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
122 extractEnclosedRegion(
123  const GridOrTreeType& volume,
124  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>(),
125  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
126  fillMask = nullptr);
127 
128 
129 /// @brief Return a mask of the voxels that intersect the implicit surface with
130 /// the given @a isovalue.
131 ///
132 /// @param volume Signed distance field / level set volume.
133 /// @param isovalue The crossing point that is considered the surface.
134 template<typename GridOrTreeType>
135 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
136 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue);
137 
138 
139 /// @brief Return a mask for each connected component of the given grid's active voxels.
140 ///
141 /// @param volume Input grid or tree
142 /// @param masks Output set of disjoint active topology masks sorted in descending order
143 /// based on the active voxel count.
144 template<typename GridOrTreeType>
145 void
146 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
147  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks);
148 
149 
150 /// @brief Separates disjoint active topology components into distinct grids or trees.
151 ///
152 /// @details Supports volumes with active tiles.
153 ///
154 /// @param volume Input grid or tree
155 /// @param segments Output set of disjoint active topology components sorted in
156 /// descending order based on the active voxel count.
157 template<typename GridOrTreeType>
158 void
159 segmentActiveVoxels(const GridOrTreeType& volume,
160  std::vector<typename GridOrTreeType::Ptr>& segments);
161 
162 
163 /// @brief Separates disjoint SDF surfaces into distinct grids or trees.
164 ///
165 /// @details Supports asymmetric interior / exterior narrowband widths and
166 /// SDF volumes with dense interior regions.
167 ///
168 /// @param volume Input signed distance field / level set volume
169 /// @param segments Output set of disjoint SDF surfaces found in @a volume sorted in
170 /// descending order based on the surface intersecting voxel count.
171 template<typename GridOrTreeType>
172 void
173 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
174 
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 ////////////////////////////////////////////////////////////////////////////////
178 
179 // Internal utility objects and implementation details
180 
181 /// @cond OPENVDB_DOCS_INTERNAL
182 
183 namespace level_set_util_internal {
184 
185 
186 template<typename LeafNodeType>
187 struct MaskInteriorVoxels {
188 
189  using ValueType = typename LeafNodeType::ValueType;
190  using BoolLeafNodeType = tree::LeafNode<bool, LeafNodeType::LOG2DIM>;
191 
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  {
2239  tree::ValueAccessor<const TreeType> acc(tree);
2240 
2241  typename TreeType::ValueAllIter it(*newTree);
2242  it.setMaxDepth(TreeType::ValueAllIter::LEAF_DEPTH - 2);
2243 
2244  for ( ; it; ++it) {
2245  if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
2246  it.setValue(ValueType(1.0));
2247  it.setActiveState(true);
2248  }
2249  }
2250  }
2251 
2252  // Insert missing root level tiles. (The new tree is constructed from the remaining leafnodes
2253  // and will therefore not contain any root level tiles that may exist in the original tree.)
2254  {
2255  typename TreeType::ValueAllIter it(tree);
2256  it.setMaxDepth(TreeType::ValueAllIter::ROOT_DEPTH);
2257  for ( ; it; ++it) {
2258  if (it.getValue() < ValueType(0.0)) {
2259  newTree->addTile(TreeType::ValueAllIter::ROOT_LEVEL, it.getCoord(),
2260  ValueType(1.0), true);
2261  }
2262  }
2263  }
2264 
2265  grid.setTree(newTree);
2266  grid.setGridClass(GRID_FOG_VOLUME);
2267 }
2268 
2269 
2270 ////////////////////////////////////////
2271 
2272 
2273 template <class GridOrTreeType>
2274 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2275 sdfInteriorMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2276 {
2277  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2278  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2279 
2280  using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2281  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(tree, isovalue);
2282 
2283  return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2284  volume, mask);
2285 }
2286 
2287 
2288 template<typename GridOrTreeType>
2289 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2290 extractEnclosedRegion(const GridOrTreeType& volume,
2291  typename GridOrTreeType::ValueType isovalue,
2292  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
2293  fillMask)
2294 {
2295  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2296  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2297 
2298  using CharTreePtrType = typename TreeType::template ValueConverter<char>::Type::Ptr;
2299  CharTreePtrType regionMask = level_set_util_internal::computeEnclosedRegionMask(
2300  tree, isovalue, fillMask);
2301 
2302  using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2303  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(*regionMask, 0);
2304 
2305  return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2306  volume, mask);
2307 }
2308 
2309 
2310 ////////////////////////////////////////
2311 
2312 
2313 template<typename GridOrTreeType>
2314 typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2315 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2316 {
2317  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2318  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2319 
2320  std::vector<const typename TreeType::LeafNodeType*> nodes;
2321  tree.getNodes(nodes);
2322 
2323  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2324  typename BoolTreeType::Ptr mask(new BoolTreeType(false));
2325 
2326  level_set_util_internal::MaskIsovalueCrossingVoxels<TreeType> op(tree, nodes, *mask, isovalue);
2327  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2328 
2329  return level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2330  volume, mask);
2331 }
2332 
2333 
2334 ////////////////////////////////////////
2335 
2336 
2337 template<typename GridOrTreeType>
2338 void
2339 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
2340  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks)
2341 {
2342  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2343  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2344  using BoolTreePtrType = typename BoolTreeType::Ptr;
2345  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2346  using BoolGridOrTreePtrType = typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr;
2347 
2348  using NodeMaskSegmentType = level_set_util_internal::NodeMaskSegment<BoolLeafNodeType>;
2349  using NodeMaskSegmentPtrType = typename NodeMaskSegmentType::Ptr;
2350  using NodeMaskSegmentPtrVector = typename std::vector<NodeMaskSegmentPtrType>;
2351  using NodeMaskSegmentRawPtrVector = typename std::vector<NodeMaskSegmentType*>;
2352 
2353  /////
2354 
2355  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2356 
2357  BoolTreeType topologyMask(tree, false, TopologyCopy());
2358 
2359  // prune out any inactive leaf nodes or inactive tiles
2360  tools::pruneInactive(topologyMask);
2361 
2362  if (topologyMask.hasActiveTiles()) {
2363  topologyMask.voxelizeActiveTiles();
2364  }
2365 
2366  std::vector<BoolLeafNodeType*> leafnodes;
2367  topologyMask.getNodes(leafnodes);
2368 
2369  if (leafnodes.empty()) return;
2370 
2371  // 1. Split node masks into disjoint segments
2372  // Note: The LeafNode origin coord is modified to record the 'leafnodes' array offset.
2373 
2374  std::unique_ptr<NodeMaskSegmentPtrVector[]> nodeSegmentArray(
2375  new NodeMaskSegmentPtrVector[leafnodes.size()]);
2376 
2377  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2378  level_set_util_internal::SegmentNodeMask<BoolLeafNodeType>(
2379  leafnodes, nodeSegmentArray.get()));
2380 
2381 
2382  // 2. Compute segment connectivity
2383 
2384  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2385  level_set_util_internal::ConnectNodeMaskSegments<BoolTreeType, BoolLeafNodeType>(
2386  topologyMask, nodeSegmentArray.get()));
2387 
2388  topologyMask.clear();
2389 
2390  size_t nodeSegmentCount = 0;
2391  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2392  nodeSegmentCount += nodeSegmentArray[n].size();
2393  }
2394 
2395  // 3. Group connected segments
2396 
2397  std::deque<NodeMaskSegmentRawPtrVector> nodeSegmentGroups;
2398 
2399  NodeMaskSegmentType* nextSegment = nodeSegmentArray[0][0].get();
2400  while (nextSegment) {
2401 
2402  nodeSegmentGroups.push_back(NodeMaskSegmentRawPtrVector());
2403 
2404  std::vector<NodeMaskSegmentType*>& segmentGroup = nodeSegmentGroups.back();
2405  segmentGroup.reserve(nodeSegmentCount);
2406 
2407  std::deque<NodeMaskSegmentType*> segmentQueue;
2408  segmentQueue.push_back(nextSegment);
2409  nextSegment = nullptr;
2410 
2411  while (!segmentQueue.empty()) {
2412 
2413  NodeMaskSegmentType* segment = segmentQueue.back();
2414  segmentQueue.pop_back();
2415 
2416  if (segment->visited) continue;
2417  segment->visited = true;
2418 
2419  segmentGroup.push_back(segment);
2420 
2421  // queue connected segments
2422  std::vector<NodeMaskSegmentType*>& connections = segment->connections;
2423  for (size_t n = 0, N = connections.size(); n < N; ++n) {
2424  if (!connections[n]->visited) segmentQueue.push_back(connections[n]);
2425  }
2426  }
2427 
2428  // find first unvisited segment
2429  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2430  NodeMaskSegmentPtrVector& nodeSegments = nodeSegmentArray[n];
2431  for (size_t i = 0, I = nodeSegments.size(); i < I; ++i) {
2432  if (!nodeSegments[i]->visited) nextSegment = nodeSegments[i].get();
2433  }
2434  }
2435  }
2436 
2437  // 4. Mask segment groups
2438 
2439  if (nodeSegmentGroups.size() == 1) {
2440 
2441  BoolTreePtrType mask(new BoolTreeType(tree, false, TopologyCopy()));
2442 
2443  tools::pruneInactive(*mask);
2444 
2445  if (mask->hasActiveTiles()) {
2446  mask->voxelizeActiveTiles();
2447  }
2448 
2449  masks.push_back(
2450  level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2451  volume, mask));
2452 
2453  } else if (nodeSegmentGroups.size() > 1) {
2454 
2455  for (size_t n = 0, N = nodeSegmentGroups.size(); n < N; ++n) {
2456 
2457  NodeMaskSegmentRawPtrVector& segmentGroup = nodeSegmentGroups[n];
2458 
2459  level_set_util_internal::MaskSegmentGroup<BoolTreeType> op(segmentGroup);
2460  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, segmentGroup.size()), op);
2461 
2462  masks.push_back(
2463  level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::constructMask(
2464  volume, op.mask()));
2465  }
2466  }
2467 
2468  // 5. Sort segments in descending order based on the active voxel count.
2469 
2470  if (masks.size() > 1) {
2471  const size_t segmentCount = masks.size();
2472 
2473  std::unique_ptr<size_t[]> segmentOrderArray(new size_t[segmentCount]);
2474  std::unique_ptr<size_t[]> voxelCountArray(new size_t[segmentCount]);
2475 
2476  for (size_t n = 0; n < segmentCount; ++n) {
2477  segmentOrderArray[n] = n;
2478  }
2479 
2480  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
2481  level_set_util_internal::ComputeActiveVoxelCount<BoolGridOrTreePtrType>(
2482  masks, voxelCountArray.get()));
2483 
2484  size_t *begin = segmentOrderArray.get();
2485  tbb::parallel_sort(begin, begin + masks.size(), level_set_util_internal::GreaterCount(
2486  voxelCountArray.get()));
2487 
2488  std::vector<BoolGridOrTreePtrType> orderedMasks;
2489  orderedMasks.reserve(masks.size());
2490 
2491  for (size_t n = 0; n < segmentCount; ++n) {
2492  orderedMasks.push_back(masks[segmentOrderArray[n]]);
2493  }
2494 
2495  masks.swap(orderedMasks);
2496  }
2497 
2498 } // extractActiveVoxelSegmentMasks()
2499 
2500 
2501 template<typename GridOrTreeType>
2502 void
2503 segmentActiveVoxels(const GridOrTreeType& volume,
2504  std::vector<typename GridOrTreeType::Ptr>& segments)
2505 {
2506  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2507  using TreePtrType = typename TreeType::Ptr;
2508  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2509  using BoolTreePtrType = typename BoolTreeType::Ptr;
2510 
2511  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2512 
2513  // 1. Segment active topology mask
2514  std::vector<BoolTreePtrType> maskSegmentArray;
2515  extractActiveVoxelSegmentMasks(inputTree, maskSegmentArray);
2516 
2517  // 2. Export segments
2518 
2519  const size_t numSegments = std::max(size_t(1), maskSegmentArray.size());
2520  std::vector<TreePtrType> outputSegmentArray(numSegments);
2521 
2522  if (maskSegmentArray.empty()) {
2523  // if no active voxels in the original volume, copy just the background
2524  // value of the input tree
2525  outputSegmentArray[0] = TreePtrType(new TreeType(inputTree.background()));
2526  } else if (numSegments == 1) {
2527  // if there's only one segment with active voxels, copy the input tree
2528  TreePtrType segment(new TreeType(inputTree));
2529  // however, if the leaf counts do not match due to the pruning of inactive leaf
2530  // nodes in the mask, do a topology intersection to drop these inactive leafs
2531  if (segment->leafCount() != inputTree.leafCount()) {
2532  segment->topologyIntersection(*maskSegmentArray[0]);
2533  }
2534  outputSegmentArray[0] = segment;
2535  } else {
2536  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2537  tbb::parallel_for(segmentRange,
2538  level_set_util_internal::MaskedCopy<TreeType>(inputTree, outputSegmentArray,
2539  maskSegmentArray));
2540  }
2541 
2542  for (auto& segment : outputSegmentArray) {
2543  segments.push_back(
2544  level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::construct(
2545  volume, segment));
2546  }
2547 }
2548 
2549 
2550 template<typename GridOrTreeType>
2551 void
2552 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments)
2553 {
2554  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2555  using TreePtrType = typename TreeType::Ptr;
2556  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2557  using BoolTreePtrType = typename BoolTreeType::Ptr;
2558 
2559  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2560 
2561  // 1. Mask zero crossing voxels
2562  BoolTreePtrType mask = extractIsosurfaceMask(inputTree, lsutilGridZero<GridOrTreeType>());
2563 
2564  // 2. Segment the zero crossing mask
2565  std::vector<BoolTreePtrType> maskSegmentArray;
2566  extractActiveVoxelSegmentMasks(*mask, maskSegmentArray);
2567 
2568  const size_t numSegments = std::max(size_t(1), maskSegmentArray.size());
2569  std::vector<TreePtrType> outputSegmentArray(numSegments);
2570 
2571  if (maskSegmentArray.empty()) {
2572  // if no active voxels in the original volume, copy just the background
2573  // value of the input tree
2574  outputSegmentArray[0] = TreePtrType(new TreeType(inputTree.background()));
2575  } else {
2576  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2577 
2578  // 3. Expand zero crossing mask to capture sdf narrow band
2579  tbb::parallel_for(segmentRange,
2580  level_set_util_internal::ExpandNarrowbandMask<TreeType>(inputTree, maskSegmentArray));
2581 
2582  // 4. Export sdf segments
2583 
2584  tbb::parallel_for(segmentRange, level_set_util_internal::MaskedCopy<TreeType>(
2585  inputTree, outputSegmentArray, maskSegmentArray));
2586 
2587  tbb::parallel_for(segmentRange,
2588  level_set_util_internal::FloodFillSign<TreeType>(inputTree, outputSegmentArray));
2589  }
2590 
2591  for (auto& segment : outputSegmentArray) {
2592  segments.push_back(
2593  level_set_util_internal::GridOrTreeConstructor<GridOrTreeType>::construct(
2594  volume, segment));
2595  }
2596 }
2597 
2598 
2599 ////////////////////////////////////////
2600 
2601 
2602 // Explicit Template Instantiation
2603 
2604 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
2605 
2606 #ifdef OPENVDB_INSTANTIATE_LEVELSETUTIL
2607 #include <openvdb/util/ExplicitInstantiation.h>
2608 #endif
2609 
2610 #define _FUNCTION(TreeT) \
2611  void sdfToFogVolume(Grid<TreeT>&, TreeT::ValueType)
2612 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2613 #undef _FUNCTION
2614 
2615 #define _FUNCTION(TreeT) \
2616  TreeT::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const TreeT&, TreeT::ValueType)
2617 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2618 #undef _FUNCTION
2619 
2620 #define _FUNCTION(TreeT) \
2621  Grid<TreeT>::ValueConverter<bool>::Type::Ptr sdfInteriorMask(const Grid<TreeT>&, TreeT::ValueType)
2622 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2623 #undef _FUNCTION
2624 
2625 #define _FUNCTION(TreeT) \
2626  TreeT::ValueConverter<bool>::Type::Ptr extractEnclosedRegion(\
2627  const TreeT&, TreeT::ValueType, \
2628  const TreeAdapter<TreeT>::TreeType::ValueConverter<bool>::Type*)
2629 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2630 #undef _FUNCTION
2631 
2632 #define _FUNCTION(TreeT) \
2633  Grid<TreeT>::ValueConverter<bool>::Type::Ptr extractEnclosedRegion(\
2634  const Grid<TreeT>&, TreeT::ValueType, \
2635  const TreeAdapter<Grid<TreeT>>::TreeType::ValueConverter<bool>::Type*)
2636 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2637 #undef _FUNCTION
2638 
2639 #define _FUNCTION(TreeT) \
2640  TreeT::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const TreeT&, TreeT::ValueType)
2641 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2642 #undef _FUNCTION
2643 
2644 #define _FUNCTION(TreeT) \
2645  Grid<TreeT>::ValueConverter<bool>::Type::Ptr extractIsosurfaceMask(const Grid<TreeT>&, TreeT::ValueType)
2646 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2647 #undef _FUNCTION
2648 
2649 #define _FUNCTION(TreeT) \
2650  void extractActiveVoxelSegmentMasks(\
2651  const TreeT&, std::vector<TreeT::ValueConverter<bool>::Type::Ptr>&)
2652 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2653 #undef _FUNCTION
2654 
2655 #define _FUNCTION(TreeT) \
2656  void extractActiveVoxelSegmentMasks(\
2657  const Grid<TreeT>&, std::vector<Grid<TreeT>::ValueConverter<bool>::Type::Ptr>&)
2658 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2659 #undef _FUNCTION
2660 
2661 #define _FUNCTION(TreeT) \
2662  void segmentActiveVoxels(const TreeT&, std::vector<TreeT::Ptr>&)
2663 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2664 #undef _FUNCTION
2665 
2666 #define _FUNCTION(TreeT) \
2667  void segmentActiveVoxels(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
2668 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
2669 #undef _FUNCTION
2670 
2671 #define _FUNCTION(TreeT) \
2672  void segmentSDF(const TreeT&, std::vector<TreeT::Ptr>&)
2673 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2674 #undef _FUNCTION
2675 
2676 #define _FUNCTION(TreeT) \
2677  void segmentSDF(const Grid<TreeT>&, std::vector<Grid<TreeT>::Ptr>&)
2678 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
2679 #undef _FUNCTION
2680 
2681 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
2682 
2683 
2684 } // namespace tools
2685 } // namespace OPENVDB_VERSION_NAME
2686 } // namespace openvdb
2687 
2688 #endif // OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
openvdb::v13_0::tools::sdfInteriorMask
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
openvdb::v13_0::GRID_FOG_VOLUME
Definition: Types.h:523
openvdb::v13_0::tree::ValueAccessorImpl
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
openvdb::v13_0::tools::traceExteriorBoundaries
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition: MeshToVolume.h:3049
std
Definition: Coord.h:590
openvdb::v13_0::tools::extractEnclosedRegion
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
openvdb::v13_0::tools::pruneInactive
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
openvdb::v13_0::tools::extractIsosurfaceMask
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
openvdb::v13_0::tools::composite::min
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
openvdb::v13_0::tools::extractActiveVoxelSegmentMasks
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&#39;s active voxels.
Definition: LevelSetUtil.h:2339
openvdb::v13_0::tools::segmentSDF
void segmentSDF(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint SDF surfaces into distinct grids or trees.
Definition: LevelSetUtil.h:2552
openvdb::v13_0::tree::ValueAccessorImpl::getValue
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:455
openvdb_houdini::Grid
openvdb::GridBase Grid
Definition: Utils.h:43
openvdb::v13_0::tools::segmentActiveVoxels
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
OPENVDB_ALL_TREE_INSTANTIATE
#define OPENVDB_ALL_TREE_INSTANTIATE(Function)
Definition: version.h.in:166
SignedFloodFill.h
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
openvdb::v13_0::TopologyCopy
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:750
ROOT_LEVEL
#define ROOT_LEVEL
Definition: CNanoVDB.h:56
MeshToVolume.h
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
nanovdb::util::empty
bool empty(const char *str)
tests if a c-string str is empty, that is its first value is &#39;\0&#39;
Definition: Util.h:144
openvdb
Definition: Exceptions.h:13
openvdb::v13_0::tools::sdfToFogVolume
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
openvdb::v13_0::tools::op
OutGridT XformOp & op
Definition: ValueTransformer.h:139
openvdb::v13_0::TreeAdapter
This adapter allows code that is templated on a Tree type to accept either a Tree type or a Grid type...
Definition: Grid.h:1055
nanovdb::GridType
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:214
openvdb::v13_0::TreeAdapter::TreeType
_TreeType TreeType
Definition: Grid.h:1057
openvdb::v13_0::tools::composite::max
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
openvdb::v13_0::tools::signedFloodFillWithValues
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
PointDataGrid.h
Attribute-owned data structure for points. Point attributes are stored in leaf nodes and ordered by v...
OPENVDB_REAL_TREE_INSTANTIATE
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:162
openvdb::v13_0::tools::signedFloodFill
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
openvdb::v13_0::Grid::Ptr
SharedPtr< Grid > Ptr
Definition: Grid.h:573
openvdb::v13_0::string::split
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
OPENVDB_VERSION_NAME
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
openvdb::v13_0::Index
Index32 Index
Definition: Types.h:34
OPENVDB_USE_VERSION_NAMESPACE
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
openvdb::v13_0::TreeAdapter::tree
static NonConstTreeType & tree(NonConstTreeType &t)
Definition: Grid.h:1072
Generated by   doxygen 1.8.11