OpenVDB  8.1.1
LevelSetUtil.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
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 
58 
59 
74 template<class GridType>
75 inline void
77  GridType& grid,
78  typename GridType::ValueType cutoffDistance = lsutilGridMax<GridType>());
79 
80 
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 
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 
132 template<typename GridOrTreeType>
133 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
134 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue);
135 
136 
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 
155 template<typename GridOrTreeType>
156 inline void
157 segmentActiveVoxels(const GridOrTreeType& volume,
158  std::vector<typename GridOrTreeType::Ptr>& segments);
159 
160 
169 template<typename GridOrTreeType>
170 inline void
171 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
172 
173 
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 
272  PopulateTree(PopulateTree& rhs, tbb::split)
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 
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 
400  void join(FindMinVoxelValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
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 
434  void join(FindMinTileValue& rhs) { minValue = std::min(minValue, rhs.minValue); }
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,
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,
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,
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 
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 
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 
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 
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];
1382  nodeMaskSegmentation(node, mNodeMaskArray[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 
1580  MaskSegmentGroup(const MaskSegmentGroup& rhs, tbb::split)
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 
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 
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 
1772  tree::ValueAccessor<const TreeType> distAcc(*mTree);
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 
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 
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 
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 
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 
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,
2287  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
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 
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 
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 
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 
2578  tbb::parallel_for(segmentRange, level_set_util_internal::MaskedCopy<TreeType>(
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
ComputeActiveVoxelCount(std::vector< VolumePtrType > &segments, size_t *countArray)
Definition: LevelSetUtil.h:2075
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:2126
void setOff(Index32 n)
Set the nth bit off.
Definition: NodeMasks.h:457
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1972
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:391
typename NodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1279
typename NodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1398
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1617
void sdfToFogVolume(GridType &grid, typename GridType::ValueType cutoffDistance=lsutilGridMax< GridType >())
Threaded method to convert a sparse level set/SDF into a sparse fog volume.
Definition: LevelSetUtil.h:2154
const Coord & origin() const
Return the grid index coordinates of this node&#39;s local origin.
Definition: LeafNode.h:172
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...
Definition: ValueAccessor.h:384
typename TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1913
ValueType const mIsovalue
Definition: LevelSetUtil.h:253
TreeType::template ValueConverter< bool >::Type::Ptr computeInteriorMask(const TreeType &tree, typename TreeType::ValueType iso)
Definition: LevelSetUtil.h:1068
Index64 onVoxelCount() const
Return the number of voxels marked On.
Definition: LeafNode.h:140
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1614
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1976
SDFTilesToFogVolume(const TreeType &tree, InternalNodeType **nodes)
Definition: LevelSetUtil.h:480
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.
Definition: ValueAccessor.h:348
typename LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:186
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2270
void addTile(Index level, const Coord &xyz, const ValueType &value, bool state)
Add a tile at the specified tree level that contains voxel (x, y, z), possibly deleting existing node...
Definition: ValueAccessor.h:337
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:252
FindMinTileValue(FindMinTileValue &rhs, tbb::split)
Definition: LevelSetUtil.h:419
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1162
MaskInteriorTiles(ValueType isovalue, const TreeType &tree, InternalNodeType **maskNodes)
Definition: LevelSetUtil.h:235
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1379
This adapter allows code that is templated on a Tree type to accept either a Tree type or a Grid type...
Definition: Grid.h:1070
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:483
void join(MaskIsovalueCrossingVoxels &rhs)
Definition: LevelSetUtil.h:1257
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1570
GreaterCount(const size_t *countArray)
Definition: LevelSetUtil.h:2094
void segmentActiveVoxels(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint active topology components into distinct grids or trees.
Definition: LevelSetUtil.h:2497
FlipRegionSign(LeafNodeType **nodes)
Definition: LevelSetUtil.h:359
TreeType const *const mTree
Definition: LevelSetUtil.h:499
typename NodeMaskSegmentType::Ptr NodeMaskSegmentTypePtr
Definition: LevelSetUtil.h:1370
FillLeafNodeVoxels(const TreeType &tree, std::vector< BoolLeafNodeType * > &maskNodes)
Definition: LevelSetUtil.h:1765
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:508
void join(PopulateTree &rhs)
Definition: LevelSetUtil.h:297
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
typename TreeType::template ValueConverter< bool >::Type::Ptr BoolTreePtrType
Definition: LevelSetUtil.h:2111
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:320
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1971
typename std::vector< NodeMaskSegmentTypePtr > NodeMaskSegmentVector
Definition: LevelSetUtil.h:1401
void join(ExpandLeafNodeRegion &rhs)
Definition: LevelSetUtil.h:1641
Definition: Coord.h:587
VolumePtrType *const mSegments
Definition: LevelSetUtil.h:2087
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:266
typename BoolTreeType::Ptr BoolTreeTypePtr
Definition: LevelSetUtil.h:1864
TreeType const *const mTree
Definition: LevelSetUtil.h:1854
typename InternalNodeType::ValueType ValueType
Definition: LevelSetUtil.h:411
typename InputTreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1161
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:238
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:260
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1862
typename InputTreeType::LeafNodeType InputLeafNodeType
Definition: LevelSetUtil.h:1160
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:425
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:280
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1871
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:509
typename NodeMaskSegmentType::Ptr NodeMaskSegmentTypePtr
Definition: LevelSetUtil.h:1400
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:522
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractEnclosedRegion(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >(), const typename TreeAdapter< GridOrTreeType >::TreeType::template ValueConverter< bool >::Type *fillMask=nullptr)
Extracts the interior regions of a signed distance field and topologically enclosed (watertight) regi...
Definition: LevelSetUtil.h:2285
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1186
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:233
Tolerance for floating-point comparison.
Definition: Math.h:147
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:404
PopulateTree(PopulateTree &rhs, tbb::split)
Definition: LevelSetUtil.h:272
Index32 Index
Definition: openvdb/Types.h:50
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:510
static BoolGridPtrType constructMask(const GridType &grid, BoolTreePtrType &maskTree)
Definition: LevelSetUtil.h:2131
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1915
ExpandNarrowbandMask(const TreeType &tree, std::vector< BoolTreeTypePtr > &segments)
Definition: LevelSetUtil.h:1866
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.
Definition: ValueAccessor.h:329
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
typename LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:357
std::vector< NodeMaskSegment * > connections
Definition: LevelSetUtil.h:1283
size_t const *const mCountArray
Definition: LevelSetUtil.h:2101
std::shared_ptr< T > SharedPtr
Definition: openvdb/Types.h:110
typename LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:377
typename TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:2110
typename TreeType::RootNodeType RootNodeType
Definition: LevelSetUtil.h:1916
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:1761
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
void segmentSDF(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint SDF surfaces into distinct grids or trees.
Definition: LevelSetUtil.h:2546
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1760
InternalNodeType **const mNodes
Definition: LevelSetUtil.h:500
typename TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1970
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
LeafNodeType **const mNodes
Definition: LevelSetUtil.h:370
typename Grid< TreeType >::Ptr GridTypePtr
Definition: LevelSetUtil.h:2123
ValueType const mIsovalue
Definition: LevelSetUtil.h:226
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
typename BoolTreeType::Ptr BoolTreePtrType
Definition: LevelSetUtil.h:2127
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1618
typename LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:445
void nodeMaskSegmentation(const NodeType &node, std::vector< typename NodeMaskSegment< NodeType >::Ptr > &segments)
Definition: LevelSetUtil.h:1292
typename BoolGridType::Ptr BoolGridPtrType
Definition: LevelSetUtil.h:2129
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
void extractActiveVoxelSegmentMasks(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::template ValueConverter< bool >::Type::Ptr > &masks)
Return a mask for each connected component of the given grid&#39;s active voxels.
Definition: LevelSetUtil.h:2334
Definition: openvdb/Exceptions.h:13
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1590
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:219
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:1974
typename BoolTreeType::Ptr BoolTreeTypePtr
Definition: LevelSetUtil.h:1975
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:255
void join(MaskSegmentGroup &rhs)
Definition: LevelSetUtil.h:1588
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1951
SharedPtr< NodeMaskSegment > Ptr
Definition: LevelSetUtil.h:1278
ValueType const mWeight
Definition: LevelSetUtil.h:473
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:37
math::Transform & transform()
Return a reference to this grid&#39;s transform, which might be shared with other grids.
Definition: Grid.h:415
typename NodeChainType::template Get< 1 > InternalNodeType
Definition: LevelSetUtil.h:1918
void join(FindMinTileValue &rhs)
Definition: LevelSetUtil.h:434
_TreeType TreeType
Definition: Grid.h:1072
typename LeafNodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1615
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractIsosurfaceMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue)
Return a mask of the voxels that intersect the implicit surface with the given isovalue.
Definition: LevelSetUtil.h:2310
MaskSegmentGroup(const std::vector< NodeMaskSegmentType * > &segments)
Definition: LevelSetUtil.h:1574
NodeType *const *const mNodes
Definition: LevelSetUtil.h:1390
typename LeafNodeType::ValueType ValueType
Definition: LevelSetUtil.h:311
TreeType const *const mTree
Definition: LevelSetUtil.h:251
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:1643
void join(FindMinVoxelValue &rhs)
Definition: LevelSetUtil.h:400
MaskIsovalueCrossingVoxels(const InputTreeType &inputTree, const std::vector< const InputLeafNodeType * > &inputLeafNodes, BoolTreeType &maskTree, InputValueType iso)
Definition: LevelSetUtil.h:1164
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...
Definition: ValueAccessor.h:379
BoolTreeTypePtr *const mSegments
Definition: LevelSetUtil.h:1906
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:226
Definition: openvdb/Types.h:334
NodeMaskType mask
Definition: LevelSetUtil.h:1284
ExpandLeafNodeRegion(const ExpandLeafNodeRegion &rhs, tbb::split)
Definition: LevelSetUtil.h:1631
MaskIsovalueCrossingVoxels(MaskIsovalueCrossingVoxels &rhs, tbb::split)
Definition: LevelSetUtil.h:1177
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1613
ValueType const mIsovalue
Definition: LevelSetUtil.h:351
typename InputTreeType::ValueType InputValueType
Definition: LevelSetUtil.h:1159
TreeTypePtr & mask()
Definition: LevelSetUtil.h:1586
MaskSegmentGroup(const MaskSegmentGroup &rhs, tbb::split)
Definition: LevelSetUtil.h:1580
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
BoolLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:225
SDFVoxelsToFogVolume(LeafNodeType **nodes, ValueType cutoffDistance)
Definition: LevelSetUtil.h:447
typename TreeType::LeafNodeType LeafNodeType
Definition: LevelSetUtil.h:261
FindMinVoxelValue(FindMinVoxelValue &rhs, tbb::split)
Definition: LevelSetUtil.h:385
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:195
TreeType & tree() const
Return a reference to the tree associated with this accessor.
Definition: ValueAccessor.h:110
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:2081
FloodFillSign(const TreeType &tree, std::vector< TreeTypePtr > &segments)
Definition: LevelSetUtil.h:1920
BoolLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:1855
static TreeType & tree(TreeType &t)
Definition: Grid.h:1087
typename RootNodeType::NodeChainType NodeChainType
Definition: LevelSetUtil.h:1917
CharLeafNodeType **const mMaskNodes
Definition: LevelSetUtil.h:350
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition: MeshToVolume.h:3019
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: LevelSetUtil.h:1863
FillMaskBoundary(const TreeType &tree, ValueType isovalue, const BoolTreeType &fillMask, const BoolLeafNodeType **fillNodes, BoolLeafNodeType **newNodes)
Definition: LevelSetUtil.h:512
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:349
MaskInteriorVoxels(ValueType isovalue, const LeafNodeType **nodes, BoolLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:189
FindMinTileValue(InternalNodeType const *const *const nodes)
Definition: LevelSetUtil.h:413
PopulateTree(TreeType &tree, LeafNodeType **leafnodes, const size_t *nodexIndexMap, ValueType background)
Definition: LevelSetUtil.h:263
Ptr copy() const
Definition: Transform.h:50
int getValueDepth(const Coord &xyz) const
Definition: ValueAccessor.h:238
typename TreeType::Ptr TreeTypePtr
Definition: LevelSetUtil.h:1571
ConnectNodeMaskSegments(const TreeType &tree, NodeMaskSegmentVector *nodeMaskArray)
Definition: LevelSetUtil.h:1403
size_t *const mCountArray
Definition: LevelSetUtil.h:2088
MaskedCopy(const TreeType &tree, std::vector< TreeTypePtr > &segments, std::vector< BoolTreeTypePtr > &masks)
Definition: LevelSetUtil.h:1978
SharedPtr< Grid > Ptr
Definition: Grid.h:579
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:1914
LeafNodeType const *const *const mNodes
Definition: LevelSetUtil.h:224
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:452
InternalNodeType **const mMaskNodes
Definition: LevelSetUtil.h:252
const NodeMaskType & getValueMask() const
Definition: LeafNode.h:867
typename LeafNodeType::NodeMaskType NodeMaskType
Definition: LevelSetUtil.h:1762
NodeMaskSegmentVector *const mNodeMaskArray
Definition: LevelSetUtil.h:1391
static GridTypePtr construct(const GridType &grid, TreeTypePtr &maskTree)
Definition: LevelSetUtil.h:2137
LeafNodeType **const mNodes
Definition: LevelSetUtil.h:472
LabelBoundaryVoxels(ValueType isovalue, const LeafNodeType **nodes, CharLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:314
Negative active values are set 0, everything else is set to 1.
Definition: LevelSetUtil.h:309
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1409
FindMinVoxelValue(LeafNodeType const *const *const leafnodes)
Definition: LevelSetUtil.h:379
static TreeTypePtr construct(const TreeType &, TreeTypePtr &tree)
Definition: LevelSetUtil.h:2115
void setOrigin(const Coord &origin)
Set the grid index coordinates of this node&#39;s local origin.
Definition: LeafNode.h:169
TreeType const *const mTree
Definition: LevelSetUtil.h:1905
BoolTreeType & newMaskTree()
Definition: LevelSetUtil.h:1639
typename std::vector< NodeMaskSegmentTypePtr > NodeMaskSegmentVector
Definition: LevelSetUtil.h:1371
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1770
static BoolTreePtrType constructMask(const TreeType &, BoolTreePtrType &maskTree)
Definition: LevelSetUtil.h:2113
ExpandLeafNodeRegion(const TreeType &distTree, BoolTreeType &maskTree, std::vector< BoolLeafNodeType * > &maskNodes)
Definition: LevelSetUtil.h:1622
typename TreeType::ValueType ValueType
Definition: LevelSetUtil.h:507
InternalNodeType const *const *const mNodes
Definition: LevelSetUtil.h:438
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
SegmentNodeMask(std::vector< NodeType * > &nodes, NodeMaskSegmentVector *nodeMaskArray)
Definition: LevelSetUtil.h:1373
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:1986
void setValueOn(const Coord &xyz)
Mark the voxel at the given coordinates as active but don&#39;t change its value.
Definition: LeafNode.h:411
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:361