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