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