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