OpenVDB  8.1.0
FastSweeping.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
26 
27 #ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
28 #define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
29 
30 //#define BENCHMARK_FAST_SWEEPING
31 
32 #include <type_traits>// for static_assert
33 #include <cmath>
34 #include <limits>
35 #include <deque>
36 #include <unordered_map>
37 #include <utility>// for std::make_pair
38 
39 #include <tbb/parallel_for.h>
40 #include <tbb/enumerable_thread_specific.h>
41 #include <tbb/task_group.h>
42 
43 #include <openvdb/math/Math.h> // for Abs() and isExactlyEqual()
44 #include <openvdb/math/Stencils.h> // for GradStencil
46 #include "LevelSetUtil.h"
47 #include "Morphology.h"
48 
49 #include "Statistics.h"
50 #ifdef BENCHMARK_FAST_SWEEPING
51 #include <openvdb/util/CpuTimer.h>
52 #endif
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
87 template<typename GridT>
88 typename GridT::Ptr
89 fogToSdf(const GridT &fogGrid,
90  typename GridT::ValueType isoValue,
91  int nIter = 1);
92 
120 template<typename GridT>
121 typename GridT::Ptr
122 sdfToSdf(const GridT &sdfGrid,
123  typename GridT::ValueType isoValue = 0,
124  int nIter = 1);
125 
158 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
159 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
160 fogToExt(const FogGridT &fogGrid,
161  const ExtOpT &op,
162  const ExtValueT& background,
163  typename FogGridT::ValueType isoValue,
164  int nIter = 1);
165 
196 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
197 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
198 sdfToExt(const SdfGridT &sdfGrid,
199  const ExtOpT &op,
200  const ExtValueT &background,
201  typename SdfGridT::ValueType isoValue = 0,
202  int nIter = 1);
203 
238 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
239 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
240 fogToSdfAndExt(const FogGridT &fogGrid,
241  const ExtOpT &op,
242  const ExtValueT &background,
243  typename FogGridT::ValueType isoValue,
244  int nIter = 1);
245 
280 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
281 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
282 sdfToSdfAndExt(const SdfGridT &sdfGrid,
283  const ExtOpT &op,
284  const ExtValueT &background,
285  typename SdfGridT::ValueType isoValue = 0,
286  int nIter = 1);
287 
305 template<typename GridT>
306 typename GridT::Ptr
307 dilateSdf(const GridT &sdfGrid,
308  int dilation,
310  int nIter = 1);
311 
330 template<typename GridT, typename MaskTreeT>
331 typename GridT::Ptr
332 maskSdf(const GridT &sdfGrid,
333  const Grid<MaskTreeT> &mask,
334  bool ignoreActiveTiles = false,
335  int nIter = 1);
336 
349 template<typename SdfGridT, typename ExtValueT = typename SdfGridT::ValueType>
351 {
352  static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
353  "FastSweeping requires SdfGridT to have floating-point values");
354  // Defined types related to the signed disntance (or fog) grid
355  using SdfValueT = typename SdfGridT::ValueType;
356  using SdfTreeT = typename SdfGridT::TreeType;
357  using SdfAccT = tree::ValueAccessor<SdfTreeT, false>;//don't register accessors
358 
359  // define types related to the extension field
360  using ExtGridT = typename SdfGridT::template ValueConverter<ExtValueT>::Type;
361  using ExtTreeT = typename ExtGridT::TreeType;
363 
364  // define types related to the tree that masks out the active voxels to be solved for
365  using SweepMaskTreeT = typename SdfTreeT::template ValueConverter<ValueMask>::Type;
366  using SweepMaskAccT = tree::ValueAccessor<SweepMaskTreeT, false>;//don't register accessors
367 
368 public:
369 
371  FastSweeping();
372 
374  ~FastSweeping() { this->clear(); }
375 
377  FastSweeping(const FastSweeping&) = delete;
378 
381 
388  typename SdfGridT::Ptr sdfGrid() { return mSdfGrid; }
389 
396  typename ExtGridT::Ptr extGrid() { return mExtGrid; }
397 
418  bool initSdf(const SdfGridT &sdfGrid, SdfValueT isoValue, bool isInputSdf);
419 
448  template <typename ExtOpT>
449  bool initExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, SdfValueT isoValue, bool isInputSdf);
450 
466  bool initDilate(const SdfGridT &sdfGrid, int dilation, NearestNeighbors nn = NN_FACE);
467 
486  template<typename MaskTreeT>
487  bool initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles = false);
488 
501  void sweep(int nIter = 1, bool finalize = true);
502 
504  void clear();
505 
507  size_t sweepingVoxelCount() const { return mSweepingVoxelCount; }
508 
510  size_t boundaryVoxelCount() const { return mBoundaryVoxelCount; }
511 
513  bool isValid() const { return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
514 
515 private:
516 
518  void computeSweepMaskLeafOrigins();
519 
520  // Private utility classes
521  template<typename>
522  struct MaskKernel;// initialization to extend a SDF into a mask
523  template<typename>
524  struct InitExt;
525  struct InitSdf;
526  struct DilateKernel;// initialization to dilate a SDF
527  struct MinMaxKernel;
528  struct SweepingKernel;// performs the actual concurrent sparse fast sweeping
529 
530  // Define the topology (i.e. stencil) of the neighboring grid points
531  static const Coord mOffset[6];// = {{-1,0,0},{1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
532 
533  // Private member data of FastSweeping
534  typename SdfGridT::Ptr mSdfGrid;
535  typename ExtGridT::Ptr mExtGrid;
536  SweepMaskTreeT mSweepMask; // mask tree containing all non-boundary active voxels
537  std::vector<Coord> mSweepMaskLeafOrigins; // cache of leaf node origins for mask tree
538  size_t mSweepingVoxelCount, mBoundaryVoxelCount;
539 };// FastSweeping
540 
542 
543 // Static member data initialization
544 template <typename SdfGridT, typename ExtValueT>
545 const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
546  {0,-1,0},{0,1,0},
547  {0,0,-1},{0,0,1}};
548 
549 template <typename SdfGridT, typename ExtValueT>
551  : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0)
552 {
553 }
554 
555 template <typename SdfGridT, typename ExtValueT>
557 {
558  mSdfGrid.reset();
559  mExtGrid.reset();
560  mSweepMask.clear();
561  mSweepingVoxelCount = mBoundaryVoxelCount = 0;
562 }
563 
564 template <typename SdfGridT, typename ExtValueT>
566 {
567  // replace any inactive leaf nodes with tiles and voxelize any active tiles
568 
569  pruneInactive(mSweepMask);
570  mSweepMask.voxelizeActiveTiles();
571 
572  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
573  using LeafT = typename SweepMaskTreeT::LeafNodeType;
574  LeafManagerT leafManager(mSweepMask);
575 
576  mSweepMaskLeafOrigins.resize(leafManager.leafCount());
577  std::atomic<size_t> sweepingVoxelCount{0};
578  auto kernel = [&](const LeafT& leaf, size_t leafIdx) {
579  mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
580  sweepingVoxelCount += leaf.onVoxelCount();
581  };
582  leafManager.foreach(kernel, /*threaded=*/true, /*grainsize=*/1024);
583 
584  mBoundaryVoxelCount = 0;
585  mSweepingVoxelCount = sweepingVoxelCount;
586  if (mSdfGrid) {
587  const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
588  assert( totalCount >= mSweepingVoxelCount );
589  mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
590  }
591 }// FastSweeping::computeSweepMaskLeafOrigins
592 
593 template <typename SdfGridT, typename ExtValueT>
594 bool FastSweeping<SdfGridT, ExtValueT>::initSdf(const SdfGridT &fogGrid, SdfValueT isoValue, bool isInputSdf)
595 {
596  this->clear();
597  mSdfGrid = fogGrid.deepCopy();//very fast
598  InitSdf kernel(*this);
599  kernel.run(isoValue, isInputSdf);
600  return this->isValid();
601 }
602 
603 template <typename SdfGridT, typename ExtValueT>
604 template <typename OpT>
605 bool FastSweeping<SdfGridT, ExtValueT>::initExt(const SdfGridT &fogGrid, const OpT &op, const ExtValueT &background, SdfValueT isoValue, bool isInputSdf)
606 {
607  this->clear();
608  mSdfGrid = fogGrid.deepCopy();//very fast
609  mExtGrid = createGrid<ExtGridT>( background );
610  mExtGrid->topologyUnion( *mSdfGrid );//very fast
611  InitExt<OpT> kernel(*this);
612  kernel.run(isoValue, op, isInputSdf);
613  return this->isValid();
614 }
615 
616 template <typename SdfGridT, typename ExtValueT>
617 bool FastSweeping<SdfGridT, ExtValueT>::initDilate(const SdfGridT &sdfGrid, int dilate, NearestNeighbors nn)
618 {
619  this->clear();
620  mSdfGrid = sdfGrid.deepCopy();//very fast
621  DilateKernel kernel(*this);
622  kernel.run(dilate, nn);
623  return this->isValid();
624 }
625 
626 template <typename SdfGridT, typename ExtValueT>
627 template<typename MaskTreeT>
628 bool FastSweeping<SdfGridT, ExtValueT>::initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles)
629 {
630  this->clear();
631  mSdfGrid = sdfGrid.deepCopy();//very fast
632 
633  if (mSdfGrid->transform() != mask.transform()) {
634  OPENVDB_THROW(RuntimeError, "FastSweeping: Mask not aligned with the grid!");
635  }
636 
637  if (mask.getGridClass() == GRID_LEVEL_SET) {
638  using T = typename MaskTreeT::template ValueConverter<bool>::Type;
639  typename Grid<T>::Ptr tmp = sdfInteriorMask(mask);//might have active tiles
640  tmp->tree().voxelizeActiveTiles();//multi-threaded
641  MaskKernel<T> kernel(*this);
642  kernel.run(tmp->tree());//multi-threaded
643  } else {
644  if (ignoreActiveTiles || !mask.tree().hasActiveTiles()) {
645  MaskKernel<MaskTreeT> kernel(*this);
646  kernel.run(mask.tree());//multi-threaded
647  } else {
648  using T = typename MaskTreeT::template ValueConverter<ValueMask>::Type;
649  T tmp(mask.tree(), false, TopologyCopy());//multi-threaded
650  tmp.voxelizeActiveTiles(true);//multi-threaded
651  MaskKernel<T> kernel(*this);
652  kernel.run(tmp);//multi-threaded
653  }
654  }
655  return this->isValid();
656 }// FastSweeping::initMask
657 
658 template <typename SdfGridT, typename ExtValueT>
659 void FastSweeping<SdfGridT, ExtValueT>::sweep(int nIter, bool finalize)
660 {
661  if (!mSdfGrid) {
662  OPENVDB_THROW(RuntimeError, "FastSweeping::sweep called before initialization");
663  }
664  if (this->boundaryVoxelCount() == 0) {
665  OPENVDB_THROW(RuntimeError, "FastSweeping: No boundary voxels found!");
666  } else if (this->sweepingVoxelCount() == 0) {
667  OPENVDB_THROW(RuntimeError, "FastSweeping: No computing voxels found!");
668  }
669 
670  // note: SweepingKernel is non copy-constructible, so use a deque instead of a vector
671  std::deque<SweepingKernel> kernels;
672  for (int i = 0; i < 4; i++) kernels.emplace_back(*this);
673 
674  { // compute voxel slices
675 #ifdef BENCHMARK_FAST_SWEEPING
676  util::CpuTimer timer("Computing voxel slices");
677 #endif
678 
679  // Exploiting nested parallelism - all voxel slice data is precomputed
680  tbb::task_group tasks;
681  tasks.run([&] { kernels[0].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]+a[2]; });/*+++ & ---*/ });
682  tasks.run([&] { kernels[1].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]-a[2]; });/*++- & --+*/ });
683  tasks.run([&] { kernels[2].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]+a[2]; });/*+-+ & -+-*/ });
684  tasks.run([&] { kernels[3].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]-a[2]; });/*+-- & -++*/ });
685  tasks.wait();
686 
687 #ifdef BENCHMARK_FAST_SWEEPING
688  timer.stop();
689 #endif
690  }
691 
692  // perform nIter iterations of bi-directional sweeping in all directions
693  for (int i = 0; i < nIter; ++i) {
694  for (SweepingKernel& kernel : kernels) kernel.sweep();
695  }
696 
697  if (finalize) {
698 #ifdef BENCHMARK_FAST_SWEEPING
699  util::CpuTimer timer("Computing extrema values");
700 #endif
701  MinMaxKernel kernel;
702  auto e = kernel.run(*mSdfGrid);//multi-threaded
703  //auto e = extrema(mGrid->beginValueOn());// 100x slower!!!!
704 #ifdef BENCHMARK_FAST_SWEEPING
705  std::cerr << "Min = " << e.min() << " Max = " << e.max() << std::endl;
706  timer.restart("Changing asymmetric background value");
707 #endif
708  changeAsymmetricLevelSetBackground(mSdfGrid->tree(), e.max(), e.min());//multi-threaded
709 
710 #ifdef BENCHMARK_FAST_SWEEPING
711  timer.stop();
712 #endif
713  }
714 }// FastSweeping::sweep
715 
719 template <typename SdfGridT, typename ExtValueT>
720 struct FastSweeping<SdfGridT, ExtValueT>::MinMaxKernel
721 {
723  using LeafRange = typename LeafMgr::LeafRange;
724  MinMaxKernel() : mMin(std::numeric_limits<SdfValueT>::max()), mMax(-mMin) {}
725  MinMaxKernel(MinMaxKernel& other, tbb::split) : mMin(other.mMin), mMax(other.mMax) {}
726 
727  math::MinMax<SdfValueT> run(const SdfGridT &grid)
728  {
729  LeafMgr mgr(grid.tree());// super fast
730  tbb::parallel_reduce(mgr.leafRange(), *this);
731  return math::MinMax<SdfValueT>(mMin, mMax);
732  }
733 
734  void operator()(const LeafRange& r)
735  {
736  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
737  for (auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
738  const SdfValueT v = *voxelIter;
739  if (v < mMin) mMin = v;
740  if (v > mMax) mMax = v;
741  }
742  }
743  }
744 
745  void join(const MinMaxKernel& other)
746  {
747  if (other.mMin < mMin) mMin = other.mMin;
748  if (other.mMax > mMax) mMax = other.mMax;
749  }
750 
751  SdfValueT mMin, mMax;
752 };// FastSweeping::MinMaxKernel
753 
755 
757 template <typename SdfGridT, typename ExtValueT>
758 struct FastSweeping<SdfGridT, ExtValueT>::DilateKernel
759 {
762  : mParent(&parent), mBackground(parent.mSdfGrid->background())
763  {
764  }
765  DilateKernel(const DilateKernel &parent) = default;// for tbb::parallel_for
767 
768  void run(int dilation, NearestNeighbors nn)
769  {
770 #ifdef BENCHMARK_FAST_SWEEPING
771  util::CpuTimer timer("Construct LeafManager");
772 #endif
773  tree::LeafManager<SdfTreeT> mgr(mParent->mSdfGrid->tree());// super fast
774 
775 #ifdef BENCHMARK_FAST_SWEEPING
776  timer.restart("Changing background value");
777 #endif
778  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
779  changeLevelSetBackground(mgr, Unknown);//multi-threaded
780 
781  #ifdef BENCHMARK_FAST_SWEEPING
782  timer.restart("Dilating and updating mgr (parallel)");
783  //timer.restart("Dilating and updating mgr (serial)");
784 #endif
785 
786  const int delta = 5;
787  for (int i=0, d = dilation/delta; i<d; ++i) dilateActiveValues(mgr, delta, nn, IGNORE_TILES);
788  dilateActiveValues(mgr, dilation % delta, nn, IGNORE_TILES);
789  //for (int i=0, n=5, d=dilation/n; i<d; ++i) dilateActiveValues(mgr, n, nn, IGNORE_TILES);
790  //dilateVoxels(mgr, dilation, nn);
791 
792 #ifdef BENCHMARK_FAST_SWEEPING
793  timer.restart("Initializing grid and sweep mask");
794 #endif
795 
796  mParent->mSweepMask.clear();
797  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
798 
800  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
801  LeafManagerT leafManager(mParent->mSdfGrid->tree());
802 
803  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
804  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
805  const SdfValueT background = mBackground;//local copy
806  auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
807  assert(maskLeaf);
808  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
809  const SdfValueT value = *voxelIter;
810  if (math::Abs(value) < background) {// disable boundary voxels from the mask tree
811  maskLeaf->setValueOff(voxelIter.pos());
812  } else {
813  voxelIter.setValue(value > 0 ? Unknown : -Unknown);
814  }
815  }
816  };
817 
818  leafManager.foreach( kernel );
819 
820  // cache the leaf node origins for fast lookup in the sweeping kernels
821 
822  mParent->computeSweepMaskLeafOrigins();
823 
824 #ifdef BENCHMARK_FAST_SWEEPING
825  timer.stop();
826 #endif
827  }// FastSweeping::DilateKernel::run
828 
829  // Private member data of DilateKernel
831  const SdfValueT mBackground;
832 };// FastSweeping::DilateKernel
833 
835 template <typename SdfGridT, typename ExtValueT>
836 struct FastSweeping<SdfGridT, ExtValueT>::InitSdf
837 {
839  InitSdf(FastSweeping &parent): mParent(&parent),
840  mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
841  InitSdf(const InitSdf&) = default;// for tbb::parallel_for
842  InitSdf& operator=(const InitSdf&) = delete;
843 
844  void run(SdfValueT isoValue, bool isInputSdf)
845  {
846  mIsoValue = isoValue;
847  mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
848  SdfTreeT &tree = mSdfGrid->tree();//sdf
849  const bool hasActiveTiles = tree.hasActiveTiles();
850 
851  if (isInputSdf && hasActiveTiles) {
852  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
853  }
854 
855 #ifdef BENCHMARK_FAST_SWEEPING
856  util::CpuTimer timer("Initialize voxels");
857 #endif
858  mParent->mSweepMask.clear();
859  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
860 
861  {// Process all voxels
862  tree::LeafManager<SdfTreeT> mgr(tree, 1);// we need one auxiliary buffer
863  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
864  mgr.swapLeafBuffer(1);//swap voxel values
865  }
866 
867 #ifdef BENCHMARK_FAST_SWEEPING
868  timer.restart("Initialize tiles - new");
869 #endif
870  // Process all tiles
871  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree);
872  mgr.foreachBottomUp(*this);//multi-threaded
873  tree.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
874  if (hasActiveTiles) tree.voxelizeActiveTiles();//multi-threaded
875 
876  // cache the leaf node origins for fast lookup in the sweeping kernels
877 
878  mParent->computeSweepMaskLeafOrigins();
879  }// FastSweeping::InitSdf::run
880 
881  void operator()(const LeafRange& r) const
882  {
883  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
884  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
885  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
886  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
887  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
888  SdfValueT* sdf = leafIter.buffer(1).data();
889  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
890  const SdfValueT value = *voxelIter;
891  const bool isAbove = value > isoValue;
892  if (!voxelIter.isValueOn()) {// inactive voxels
893  sdf[voxelIter.pos()] = isAbove ? above : -above;
894  } else {// active voxels
895  const Coord ijk = voxelIter.getCoord();
896  stencil.moveTo(ijk, value);
897  const auto mask = stencil.intersectionMask( isoValue );
898  if (mask.none()) {// most common case
899  sdf[voxelIter.pos()] = isAbove ? above : -above;
900  } else {// compute distance to iso-surface
901  // disable boundary voxels from the mask tree
902  sweepMaskAcc.setValueOff(ijk);
903  const SdfValueT delta = value - isoValue;//offset relative to iso-value
904  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
905  sdf[voxelIter.pos()] = 0;
906  } else {//voxel is neighboring the iso-surface
907  SdfValueT sum = 0;
908  for (int i=0; i<6;) {
909  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
910  if (mask.test(i++)) d = math::Abs(delta/(value-stencil.getValue(i)));
911  if (mask.test(i++)) {
912  d2 = math::Abs(delta/(value-stencil.getValue(i)));
913  if (d2 < d) d = d2;
914  }
915  if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
916  }
917  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum) : -h / math::Sqrt(sum);
918  }// voxel is neighboring the iso-surface
919  }// intersecting voxels
920  }// active voxels
921  }// loop over voxels
922  }// loop over leaf nodes
923  }// FastSweeping::InitSdf::operator(const LeafRange&)
924 
925  template<typename RootOrInternalNodeT>
926  void operator()(const RootOrInternalNodeT& node) const
927  {
928  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
929  for (auto it = node.cbeginValueAll(); it; ++it) {
930  SdfValueT& v = const_cast<SdfValueT&>(*it);
931  v = v > isoValue ? above : -above;
932  }//loop over all tiles
933  }// FastSweeping::InitSdf::operator()(const RootOrInternalNodeT&)
934 
935  // Public member data
937  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
938  SdfValueT mIsoValue;
939  SdfValueT mAboveSign;//sign of distance values above the iso-value
940 };// FastSweeping::InitSdf
941 
942 
944 template <typename SdfGridT, typename ExtValueT>
945 template <typename OpT>
946 struct FastSweeping<SdfGridT, ExtValueT>::InitExt
947 {
948  using LeafRange = typename tree::LeafManager<SdfTreeT>::LeafRange;
949  using OpPoolT = tbb::enumerable_thread_specific<OpT>;
950  InitExt(FastSweeping &parent) : mParent(&parent),
951  mOpPool(nullptr), mSdfGrid(parent.mSdfGrid.get()),
952  mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
953  InitExt(const InitExt&) = default;// for tbb::parallel_for
954  InitExt& operator=(const InitExt&) = delete;
955  void run(SdfValueT isoValue, const OpT &opPrototype, bool isInputSdf)
956  {
957  static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::value, "Invalid return type of functor");
958  if (!mExtGrid) {
959  OPENVDB_THROW(RuntimeError, "FastSweeping::InitExt expected an extension grid!");
960  }
961 
962  mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
963  mIsoValue = isoValue;
964  auto &tree1 = mSdfGrid->tree();
965  auto &tree2 = mExtGrid->tree();
966  const bool hasActiveTiles = tree1.hasActiveTiles();//very fast
967 
968  if (isInputSdf && hasActiveTiles) {
969  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
970  }
971 
972 #ifdef BENCHMARK_FAST_SWEEPING
973  util::CpuTimer timer("Initialize voxels");
974 #endif
975 
976  mParent->mSweepMask.clear();
977  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
978 
979  {// Process all voxels
980  // Define thread-local operators
981  OpPoolT opPool(opPrototype);
982  mOpPool = &opPool;
983 
984  tree::LeafManager<SdfTreeT> mgr(tree1, 1);// we need one auxiliary buffer
985  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
986  mgr.swapLeafBuffer(1);//swap out auxiliary buffer
987  }
988 
989 #ifdef BENCHMARK_FAST_SWEEPING
990  timer.restart("Initialize tiles");
991 #endif
992  {// Process all tiles
993  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
994  mgr.foreachBottomUp(*this);//multi-threaded
995  tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
996  if (hasActiveTiles) {
997 #ifdef BENCHMARK_FAST_SWEEPING
998  timer.restart("Voxelizing active tiles");
999 #endif
1000  tree1.voxelizeActiveTiles();//multi-threaded
1001  tree2.voxelizeActiveTiles();//multi-threaded
1002  }
1003  }
1004 
1005  // cache the leaf node origins for fast lookup in the sweeping kernels
1006 
1007  mParent->computeSweepMaskLeafOrigins();
1008 
1009 #ifdef BENCHMARK_FAST_SWEEPING
1010  timer.stop();
1011 #endif
1012  }// FastSweeping::InitExt::run
1013 
1014  // int implements an update if minD needs to be updated
1015  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1016  void sumHelper(ExtT& sum2, ExtT ext, bool update, const SdfT& /* d2 */) const { if (update) sum2 = ext; }// int implementation
1017 
1018  // non-int implements a weighted sum
1019  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1020  void sumHelper(ExtT& sum2, ExtT ext, bool /* update */, const SdfT& d2) const { sum2 += static_cast<ExtValueT>(d2 * ext); }// non-int implementation
1021 
1022  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1023  ExtT extValHelper(ExtT extSum, const SdfT& /* sdfSum */) const { return extSum; }// int implementation
1024 
1025  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1026  ExtT extValHelper(ExtT extSum, const SdfT& sdfSum) const {return ExtT((SdfT(1) / sdfSum) * extSum); }// non-int implementation
1027 
1028  void operator()(const LeafRange& r) const
1029  {
1030  ExtAccT acc(mExtGrid->tree());
1031  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1032  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1033  const math::Transform& xform = mExtGrid->transform();
1034  typename OpPoolT::reference op = mOpPool->local();
1035  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1036  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1037  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1038  SdfValueT *sdf = leafIter.buffer(1).data();
1039  ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();//should be safe!
1040  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1041  const SdfValueT value = *voxelIter;
1042  const bool isAbove = value > isoValue;
1043  if (!voxelIter.isValueOn()) {// inactive voxels
1044  sdf[voxelIter.pos()] = isAbove ? above : -above;
1045  } else {// active voxels
1046  const Coord ijk = voxelIter.getCoord();
1047  stencil.moveTo(ijk, value);
1048  const auto mask = stencil.intersectionMask( isoValue );
1049  if (mask.none()) {// no zero-crossing neighbors, most common case
1050  sdf[voxelIter.pos()] = isAbove ? above : -above;
1051  // the ext grid already has its active values set to the bakground value
1052  } else {// compute distance to iso-surface
1053  // disable boundary voxels from the mask tree
1054  sweepMaskAcc.setValueOff(ijk);
1055  const SdfValueT delta = value - isoValue;//offset relative to iso-value
1056  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1057  sdf[voxelIter.pos()] = 0;
1058  ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1059  } else {//voxel is neighboring the iso-surface
1060  SdfValueT sum1 = 0;
1061  ExtValueT sum2 = zeroVal<ExtValueT>();
1062  // minD is used to update sum2 in the integer case,
1063  // where we choose the value of the extension grid corresponding to
1064  // the smallest value of d in the 6 direction neighboring the current
1065  // voxel.
1066  SdfValueT minD = std::numeric_limits<SdfValueT>::max();
1067  for (int n=0, i=0; i<6;) {
1068  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1069  if (mask.test(i++)) {
1070  d = math::Abs(delta/(value-stencil.getValue(i)));
1071  n = i - 1;
1072  }
1073  if (mask.test(i++)) {
1074  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1075  if (d2 < d) {
1076  d = d2;
1077  n = i - 1;
1078  }
1079  }
1081  d2 = 1/(d*d);
1082  sum1 += d2;
1083  const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1084  static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1085  static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1086  // If current d is less than minD, update minD
1087  sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1088  if (d < minD) minD = d;
1089  }
1090  }//look over six cases
1091  ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1092  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum1) : -h / math::Sqrt(sum1);
1093  }// voxel is neighboring the iso-surface
1094  }// intersecting voxels
1095  }// active voxels
1096  }// loop over voxels
1097  }// loop over leaf nodes
1098  }// FastSweeping::InitExt::operator(const LeafRange& r)
1099 
1100  template<typename RootOrInternalNodeT>
1101  void operator()(const RootOrInternalNodeT& node) const
1102  {
1103  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1104  for (auto it = node.cbeginValueAll(); it; ++it) {
1105  SdfValueT& v = const_cast<SdfValueT&>(*it);
1106  v = v > isoValue ? above : -above;
1107  }//loop over all tiles
1108  }
1109  // Public member data
1110  FastSweeping *mParent;
1111  OpPoolT *mOpPool;
1112  SdfGridT *mSdfGrid;
1113  ExtGridT *mExtGrid;
1114  SdfValueT mIsoValue;
1115  SdfValueT mAboveSign;//sign of distance values above the iso-value
1116 };// FastSweeping::InitExt
1117 
1119 template <typename SdfGridT, typename ExtValueT>
1120 template <typename MaskTreeT>
1121 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1122 {
1123  using LeafRange = typename tree::LeafManager<const MaskTreeT>::LeafRange;
1124  MaskKernel(FastSweeping &parent) : mParent(&parent),
1125  mSdfGrid(parent.mSdfGrid.get()) {}
1126  MaskKernel(const MaskKernel &parent) = default;// for tbb::parallel_for
1127  MaskKernel& operator=(const MaskKernel&) = delete;
1128 
1129  void run(const MaskTreeT &mask)
1130  {
1131 #ifdef BENCHMARK_FAST_SWEEPING
1132  util::CpuTimer timer;
1133 #endif
1134  auto &lsTree = mSdfGrid->tree();
1135 
1136  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1137 
1138 #ifdef BENCHMARK_FAST_SWEEPING
1139  timer.restart("Changing background value");
1140 #endif
1141  changeLevelSetBackground(lsTree, Unknown);//multi-threaded
1142 
1143 #ifdef BENCHMARK_FAST_SWEEPING
1144  timer.restart("Union with mask");//multi-threaded
1145 #endif
1146  lsTree.topologyUnion(mask);//multi-threaded
1147 
1148  // ignore active tiles since the input grid is assumed to be a level set
1149  tree::LeafManager<const MaskTreeT> mgr(mask);// super fast
1150 
1151 #ifdef BENCHMARK_FAST_SWEEPING
1152  timer.restart("Initializing grid and sweep mask");
1153 #endif
1154 
1155  mParent->mSweepMask.clear();
1156  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1157 
1158  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1159  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1160  LeafManagerT leafManager(mParent->mSweepMask);
1161 
1162  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
1163  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1164  SdfAccT acc(mSdfGrid->tree());
1165  // The following hack is safe due to the topoloyUnion in
1166  // init and the fact that SdfValueT is known to be a floating point!
1167  SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1168  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {// mask voxels
1169  if (math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1170  // disable boundary voxels from the mask tree
1171  voxelIter.setValue(false);
1172  }
1173  }
1174  };
1175  leafManager.foreach( kernel );
1176 
1177  // cache the leaf node origins for fast lookup in the sweeping kernels
1178  mParent->computeSweepMaskLeafOrigins();
1179 
1180 #ifdef BENCHMARK_FAST_SWEEPING
1181  timer.stop();
1182 #endif
1183  }// FastSweeping::MaskKernel::run
1184 
1185  // Private member data of MaskKernel
1186  FastSweeping *mParent;
1187  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1188 };// FastSweeping::MaskKernel
1189 
1191 template <typename SdfGridT, typename ExtValueT>
1192 struct FastSweeping<SdfGridT, ExtValueT>::SweepingKernel
1193 {
1194  SweepingKernel(FastSweeping &parent) : mParent(&parent) {}
1195  SweepingKernel(const SweepingKernel&) = delete;
1197 
1199  template<typename HashOp>
1200  void computeVoxelSlices(HashOp hash)
1201  {
1202 #ifdef BENCHMARK_FAST_SWEEPING
1203  util::CpuTimer timer;
1204 #endif
1205 
1206  // mask of the active voxels to be solved for, i.e. excluding boundary voxels
1207  const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1208 
1209  using LeafManagerT = typename tree::LeafManager<const SweepMaskTreeT>;
1210  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1211  LeafManagerT leafManager(maskTree);
1212 
1213  // compute the leaf node slices that have active voxels in them
1214  // the sliding window of the has keys is -14 to 21 (based on an 8x8x8 leaf node
1215  // and the extrema hash values i-j-k and i+j+k), but we use a larger mask window here to
1216  // easily accomodate any leaf dimension. The mask offset is used to be able to
1217  // store this in a fixed-size byte array
1218  constexpr int maskOffset = LeafT::DIM * 3;
1219  constexpr int maskRange = maskOffset * 2;
1220 
1221  // mark each possible slice in each leaf node that has one or more active voxels in it
1222  std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1223  auto kernel1 = [&](const LeafT& leaf, size_t leafIdx) {
1224  const size_t leafOffset = leafIdx * maskRange;
1225  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1226  const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1227  leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1228  }
1229  };
1230  leafManager.foreach( kernel1 );
1231 
1232  // compute the voxel slice map using a thread-local-storage hash map
1233  // the key of the hash map is the slice index of the voxel coord (ijk.x() + ijk.y() + ijk.z())
1234  // the values are an array of indices for every leaf that has active voxels with this slice index
1235  using ThreadLocalMap = std::unordered_map</*voxelSliceKey=*/int64_t, /*leafIdx=*/std::deque<size_t>>;
1236  tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1237  auto kernel2 = [&](const LeafT& leaf, size_t leafIdx) {
1238  ThreadLocalMap& map = pool.local();
1239  const Coord& origin = leaf.origin();
1240  const int64_t leafKey = hash(origin);
1241  const size_t leafOffset = leafIdx * maskRange;
1242  for (int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1243  if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1244  const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1245  map[voxelSliceKey].emplace_back(leafIdx);
1246  }
1247  }
1248  };
1249  leafManager.foreach( kernel2 );
1250 
1251  // combine into a single ordered map keyed by the voxel slice key
1252  // note that this is now stored in a map ordered by voxel slice key,
1253  // so sweep slices can be processed in order
1254  for (auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1255  const ThreadLocalMap& map = *poolIt;
1256  for (const auto& it : map) {
1257  for (const size_t leafIdx : it.second) {
1258  mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1259  }
1260  }
1261  }
1262 
1263  // extract the voxel slice keys for random access into the map
1264  mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1265  for (const auto& it : mVoxelSliceMap) {
1266  mVoxelSliceKeys.push_back(it.first);
1267  }
1268 
1269  // allocate the node masks in parallel, as the map is populated in serial
1270  auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1271  for (size_t i = range.begin(); i < range.end(); i++) {
1272  const int64_t key = mVoxelSliceKeys[i];
1273  for (auto& it : mVoxelSliceMap[key]) {
1274  it.second = std::make_unique<NodeMaskT>();
1275  }
1276  }
1277  };
1278  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1279 
1280  // each voxel slice contains a leafIdx-nodeMask pair,
1281  // this routine populates these node masks to select only the active voxels
1282  // from the mask tree that have the same voxel slice key
1283  // TODO: a small optimization here would be to union this leaf node mask with
1284  // a pre-computed one for this particular slice pattern
1285  auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1286  for (size_t i = range.begin(); i < range.end(); i++) {
1287  const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1288  LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1289  for (LeafSlice& leafSlice : leafSliceArray) {
1290  const size_t leafIdx = leafSlice.first;
1291  NodeMaskPtrT& nodeMask = leafSlice.second;
1292  const LeafT& leaf = leafManager.leaf(leafIdx);
1293  const Coord& origin = leaf.origin();
1294  const int64_t leafKey = hash(origin);
1295  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1296  const Index voxelIdx = voxelIter.pos();
1297  const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1298  const int64_t key = leafKey + hash(ijk);
1299  if (key == voxelSliceKey) {
1300  nodeMask->setOn(voxelIdx);
1301  }
1302  }
1303  }
1304  }
1305  };
1306  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1307  }// FastSweeping::SweepingKernel::computeVoxelSlices
1308 
1309  // Private struct for nearest neighbor grid points (very memory light!)
1310  struct NN {
1311  SdfValueT v;
1312  int n;
1313  inline static Coord ijk(const Coord &p, int i) { return p + FastSweeping::mOffset[i]; }
1314  NN() : v(), n() {}
1315  NN(const SdfAccT &a, const Coord &p, int i) : v(math::Abs(a.getValue(ijk(p,i)))), n(i) {}
1316  inline Coord operator()(const Coord &p) const { return ijk(p, n); }
1317  inline bool operator<(const NN &rhs) const { return v < rhs.v; }
1318  inline operator bool() const { return v < SdfValueT(1000); }
1319  };// NN
1320 
1322  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1323  ExtT twoNghbr(const NN& d1, const NN& d2, const SdfT& /* w */, const ExtT& v1, const ExtT& v2) const { return d1.v < d2.v ? v1 : v2; }// int implementation
1324 
1326  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1327  ExtT twoNghbr(const NN& d1, const NN& d2, const SdfT& w, const ExtT& v1, const ExtT& v2) const { return ExtT(w*(d1.v*v1 + d2.v*v2)); }// non-int implementation
1328 
1330  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1331  ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& /* w */, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1332  math::Vec3<SdfT> d(d1.v, d2.v, d3.v);
1333  math::Vec3<ExtT> v(v1, v2, v3);
1334  return v[static_cast<int>(math::MinIndex(d))];
1335  }// int implementation
1336 
1338  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1339  ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& w, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1340  return ExtT(w*(d1.v*v1 + d2.v*v2 + d3.v*v3));
1341  }// non-int implementation
1342 
1343  void sweep()
1344  {
1345  typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() : nullptr;
1346 
1347  const SdfValueT h = static_cast<SdfValueT>(mParent->mSdfGrid->voxelSize()[0]);
1348  const SdfValueT sqrt2h = math::Sqrt(SdfValueT(2))*h;
1349 
1350  const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1351 
1352  int64_t voxelSliceIndex(0);
1353 
1354  auto kernel = [&](const tbb::blocked_range<size_t>& range) {
1355  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
1356 
1357  SdfAccT acc1(mParent->mSdfGrid->tree());
1358  auto acc2 = std::unique_ptr<ExtAccT>(tree2 ? new ExtAccT(*tree2) : nullptr);
1359  SdfValueT absV, sign, update, D;
1360  NN d1, d2, d3;//distance values and coordinates of closest neighbor points
1361 
1362  const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1363 
1364  // Solves Goudonov's scheme: [x-d1]^2 + [x-d2]^2 + [x-d3]^2 = h^2
1365  // where [X] = (X>0?X:0) and ai=min(di+1,di-1)
1366  for (size_t i = range.begin(); i < range.end(); ++i) {
1367 
1368  // iterate over all leafs in the slice and extract the leaf
1369  // and node mask for each slice pattern
1370 
1371  const LeafSlice& leafSlice = leafSliceArray[i];
1372  const size_t leafIdx = leafSlice.first;
1373  const NodeMaskPtrT& nodeMask = leafSlice.second;
1374 
1375  const Coord& origin = leafNodeOrigins[leafIdx];
1376 
1377  Coord ijk;
1378  for (auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1379 
1380  // Get coordinate of center point of the FD stencil
1381  ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1382 
1383  // Find the closes neighbors in the three axial directions
1384  d1 = std::min(NN(acc1, ijk, 0), NN(acc1, ijk, 1));
1385  d2 = std::min(NN(acc1, ijk, 2), NN(acc1, ijk, 3));
1386  d3 = std::min(NN(acc1, ijk, 4), NN(acc1, ijk, 5));
1387 
1388  if (!(d1 || d2 || d3)) continue;//no valid neighbors
1389 
1390  // Get the center point of the FD stencil (assumed to be an active voxel)
1391  // Note this const_cast is normally unsafe but by design we know the tree
1392  // to be static, of floating-point type and containing active voxels only!
1393  SdfValueT &value = const_cast<SdfValueT&>(acc1.getValue(ijk));
1394 
1395  // Extract the sign
1396  sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1397 
1398  // Absolute value
1399  absV = math::Abs(value);
1400 
1401  // sort values so d1 <= d2 <= d3
1402  if (d2 < d1) std::swap(d1, d2);
1403  if (d3 < d2) std::swap(d2, d3);
1404  if (d2 < d1) std::swap(d1, d2);
1405 
1406  // Test if there is a solution depending on ONE of the neighboring voxels
1407  // if d2 - d1 >= h => d2 >= d1 + h then:
1408  // (x-d1)^2=h^2 => x = d1 + h
1409  update = d1.v + h;
1410  if (update <= d2.v) {
1411  if (update < absV) {
1412  value = sign * update;
1413  if (acc2) acc2->setValue(ijk, acc2->getValue(d1(ijk)));//update ext?
1414  }//update sdf?
1415  continue;
1416  }// one neighbor case
1417 
1418  // Test if there is a solution depending on TWO of the neighboring voxels
1419  // (x-d1)^2 + (x-d2)^2 = h^2
1420  //D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1421  //if (D >= SdfValueT(0)) {// non-negative discriminant
1422  if (d2.v <= sqrt2h + d1.v) {
1423  D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1424  update = SdfValueT(0.5) * (d1.v + d2.v + std::sqrt(D));
1425  if (update > d2.v && update <= d3.v) {
1426  if (update < absV) {
1427  value = sign * update;
1428  if (acc2) {
1429  d1.v -= update;
1430  d2.v -= update;
1431  // affine combination of two neighboring extension values
1432  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v);
1433  const ExtValueT v1 = acc2->getValue(d1(ijk));
1434  const ExtValueT v2 = acc2->getValue(d2(ijk));
1435  const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1436  acc2->setValue(ijk, extVal);
1437  }//update ext?
1438  }//update sdf?
1439  continue;
1440  }//test for two neighbor case
1441  }//test for non-negative determinant
1442 
1443  // Test if there is a solution depending on THREE of the neighboring voxels
1444  // (x-d1)^2 + (x-d2)^2 + (x-d3)^2 = h^2
1445  // 3x^2 - 2(d1 + d2 + d3)x + d1^2 + d2^2 + d3^2 = h^2
1446  // ax^2 + bx + c=0, a=3, b=-2(d1+d2+d3), c=d1^2 + d2^2 + d3^2 - h^2
1447  const SdfValueT d123 = d1.v + d2.v + d3.v;
1448  D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.v*d2.v + d3.v*d3.v - h * h);
1449  if (D >= SdfValueT(0)) {// non-negative discriminant
1450  update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));//always passes test
1451  //if (update > d3.v) {//disabled due to round-off errors
1452  if (update < absV) {
1453  value = sign * update;
1454  if (acc2) {
1455  d1.v -= update;
1456  d2.v -= update;
1457  d3.v -= update;
1458  // affine combination of three neighboring extension values
1459  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v+d3.v);
1460  const ExtValueT v1 = acc2->getValue(d1(ijk));
1461  const ExtValueT v2 = acc2->getValue(d2(ijk));
1462  const ExtValueT v3 = acc2->getValue(d3(ijk));
1463  const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1464  acc2->setValue(ijk, extVal);
1465  }//update ext?
1466  }//update sdf?
1467  }//test for non-negative determinant
1468  }//loop over coordinates
1469  }
1470  };
1471 
1472 #ifdef BENCHMARK_FAST_SWEEPING
1473  util::CpuTimer timer("Forward sweep");
1474 #endif
1475 
1476  for (size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1477  voxelSliceIndex = mVoxelSliceKeys[i];
1478  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1479  }
1480 
1481 #ifdef BENCHMARK_FAST_SWEEPING
1482  timer.restart("Backward sweeps");
1483 #endif
1484  for (size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1485  voxelSliceIndex = mVoxelSliceKeys[i-1];
1486  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1487  }
1488 
1489 #ifdef BENCHMARK_FAST_SWEEPING
1490  timer.stop();
1491 #endif
1492  }// FastSweeping::SweepingKernel::sweep
1493 
1494 private:
1495  using NodeMaskT = typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1496  using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1497  // using a unique ptr for the NodeMask allows for parallel allocation,
1498  // but makes this class not copy-constructible
1499  using LeafSlice = std::pair</*leafIdx=*/size_t, /*leafMask=*/NodeMaskPtrT>;
1500  using LeafSliceArray = std::deque<LeafSlice>;
1501  using VoxelSliceMap = std::map</*voxelSliceKey=*/int64_t, LeafSliceArray>;
1502 
1503  // Private member data of SweepingKernel
1504  FastSweeping *mParent;
1505  VoxelSliceMap mVoxelSliceMap;
1506  std::vector<int64_t> mVoxelSliceKeys;
1507 };// FastSweeping::SweepingKernel
1508 
1510 
1511 template<typename GridT>
1512 typename GridT::Ptr
1513 fogToSdf(const GridT &fogGrid,
1514  typename GridT::ValueType isoValue,
1515  int nIter)
1516 {
1518  if (fs.initSdf(fogGrid, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1519  return fs.sdfGrid();
1520 }
1521 
1522 template<typename GridT>
1523 typename GridT::Ptr
1524 sdfToSdf(const GridT &sdfGrid,
1525  typename GridT::ValueType isoValue,
1526  int nIter)
1527 {
1529  if (fs.initSdf(sdfGrid, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1530  return fs.sdfGrid();
1531 }
1532 
1533 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1534 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1535 fogToExt(const FogGridT &fogGrid,
1536  const ExtOpT &op,
1537  const ExtValueT& background,
1538  typename FogGridT::ValueType isoValue,
1539  int nIter)
1540 {
1542  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1543  return fs.extGrid();
1544 }
1545 
1546 template<typename SdfGridT, typename OpT, typename ExtValueT>
1547 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1548 sdfToExt(const SdfGridT &sdfGrid,
1549  const OpT &op,
1550  const ExtValueT &background,
1551  typename SdfGridT::ValueType isoValue,
1552  int nIter)
1553 {
1555  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1556  return fs.extGrid();
1557 }
1558 
1559 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1560 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1561 fogToSdfAndExt(const FogGridT &fogGrid,
1562  const ExtOpT &op,
1563  const ExtValueT &background,
1564  typename FogGridT::ValueType isoValue,
1565  int nIter)
1566 {
1568  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1569  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1570 }
1571 
1572 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
1573 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1574 sdfToSdfAndExt(const SdfGridT &sdfGrid,
1575  const ExtOpT &op,
1576  const ExtValueT &background,
1577  typename SdfGridT::ValueType isoValue,
1578  int nIter)
1579 {
1581  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1582  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1583 }
1584 
1585 template<typename GridT>
1586 typename GridT::Ptr
1587 dilateSdf(const GridT &sdfGrid,
1588  int dilation,
1589  NearestNeighbors nn,
1590  int nIter)
1591 {
1593  if (fs.initDilate(sdfGrid, dilation, nn)) fs.sweep(nIter);
1594  return fs.sdfGrid();
1595 }
1596 
1597 template<typename GridT, typename MaskTreeT>
1598 typename GridT::Ptr
1599 maskSdf(const GridT &sdfGrid,
1600  const Grid<MaskTreeT> &mask,
1601  bool ignoreActiveTiles,
1602  int nIter)
1603 {
1605  if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1606  return fs.sdfGrid();
1607 }
1608 
1609 } // namespace tools
1610 } // namespace OPENVDB_VERSION_NAME
1611 } // namespace openvdb
1612 
1613 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean,...
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:415
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:577
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:917
SharedPtr< Grid > Ptr
Definition: Grid.h:579
Definition: openvdb/Exceptions.h:63
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
std::bitset< 6 > intersectionMask(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true a bit-mask where the 6 bits indicates if the center of the stencil intersects the iso-con...
Definition: Stencils.h:188
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:26
Definition: Stencils.h:1232
Templated class to compute the minimum and maximum values.
Definition: Stats.h:31
Definition: Vec3.h:24
Computes signed distance values from an initial iso-surface and optionally performs velocty extension...
Definition: FastSweeping.h:351
bool initMask(const SdfGridT &sdfGrid, const Grid< MaskTreeT > &mask, bool ignoreActiveTiles=false)
Initializer used for the extamnsion of an exsiting signed distance field into the active values of an...
Definition: FastSweeping.h:628
size_t boundaryVoxelCount() const
Return the number of voxels that defined the boundary condition.
Definition: FastSweeping.h:510
FastSweeping()
Constructor.
Definition: FastSweeping.h:550
bool isValid() const
Return true if there are voxels and boundaries to solve for.
Definition: FastSweeping.h:513
size_t sweepingVoxelCount() const
Return the number of voxels that will be solved for.
Definition: FastSweeping.h:507
~FastSweeping()
Destructor.
Definition: FastSweeping.h:374
bool initDilate(const SdfGridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE)
Initializer used when dilating an exsiting signed distance field.
Definition: FastSweeping.h:617
FastSweeping & operator=(const FastSweeping &)=delete
Disallow copy assignment.
bool initSdf(const SdfGridT &sdfGrid, SdfValueT isoValue, bool isInputSdf)
Initializer for input grids that are either a signed distance field or a scalar fog volume.
Definition: FastSweeping.h:594
SdfGridT::Ptr sdfGrid()
Returns a shared pointer to the signed distance field computed by this class.
Definition: FastSweeping.h:388
void sweep(int nIter=1, bool finalize=true)
Perform nIter iterations of the fast sweeping algorithm.
Definition: FastSweeping.h:659
ExtGridT::Ptr extGrid()
Returns a shared pointer to the extension field computed by this class.
Definition: FastSweeping.h:396
void clear()
Clears all the grids and counters so initializtion can be called again.
Definition: FastSweeping.h:556
FastSweeping(const FastSweeping &)=delete
Disallow copy construction.
bool initExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, SdfValueT isoValue, bool isInputSdf)
Initializer used whenever velocity extension is performed in addition to the computation of signed di...
Definition: LeafManager.h:102
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node's buffer with the nth corresponding auxiliary buffer, where n = bufferIdx.
Definition: LeafManager.h:359
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays,...
Definition: NodeManager.h:531
Definition: ValueAccessor.h:183
void setValueOff(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as inactive.
Definition: ValueAccessor.h:266
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:219
Simple timer for basic profiling.
Definition: CpuTimer.h:67
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
double restart()
Re-start timer.
Definition: CpuTimer.h:150
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance.
Definition: Math.h:350
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:938
Vec3< double > Vec3d
Definition: Vec3.h:668
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
GridT::Ptr fogToSdf(const GridT &fogGrid, typename GridT::ValueType isoValue, int nIter=1)
Converts a scalar fog volume into a signed distance function. Active input voxels with scalar values ...
Definition: FastSweeping.h:1513
NearestNeighbors
Voxel topology of nearest neighbors.
Definition: Morphology.h:56
@ NN_FACE
Definition: Morphology.h:56
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2270
GridT::Ptr sdfToSdf(const GridT &sdfGrid, typename GridT::ValueType isoValue=0, int nIter=1)
Given an existing approximate SDF it solves the Eikonal equation for all its active voxels....
Definition: FastSweeping.h:1524
GridT::Ptr dilateSdf(const GridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE, int nIter=1)
Dilates an existing signed distance field by a specified number of voxels.
Definition: FastSweeping.h:1587
std::pair< typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr > sdfToSdfAndExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue=0, int nIter=1)
Computes the signed distance field and the extension of a scalar field, defined by the specified func...
Definition: FastSweeping.h:1574
void changeAsymmetricLevelSetBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=32)
Replace the background values in all the nodes of a floating-point tree containing a possibly asymmet...
Definition: ChangeBackground.h:217
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1049
GridT::Ptr maskSdf(const GridT &sdfGrid, const Grid< MaskTreeT > &mask, bool ignoreActiveTiles=false, int nIter=1)
Fills mask by extending an existing signed distance field into the active values of this input ree of...
Definition: FastSweeping.h:1599
SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr sdfToExt(const SdfGridT &sdfGrid, const OpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue, int nIter)
Definition: FastSweeping.h:1548
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3705
std::pair< typename FogGridT::Ptr, typename FogGridT::template ValueConverter< ExtValueT >::Type::Ptr > fogToSdfAndExt(const FogGridT &fogGrid, const ExtOpT &op, const ExtValueT &background, typename FogGridT::ValueType isoValue, int nIter=1)
Computes the signed distance field and the extension of a scalar field, defined by the specified func...
Definition: FastSweeping.h:1561
@ IGNORE_TILES
Definition: Morphology.h:78
FogGridT::template ValueConverter< ExtValueT >::Type::Ptr fogToExt(const FogGridT &fogGrid, const ExtOpT &op, const ExtValueT &background, typename FogGridT::ValueType isoValue, int nIter=1)
Computes the extension of a field, defined by the specified functor, off an iso-surface from an input...
Definition: FastSweeping.h:1535
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 changeLevelSetBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &halfWidth, bool threaded=true, size_t grainSize=32)
Replace the background value in all the nodes of a floating-point tree containing a symmetric narrow-...
Definition: ChangeBackground.h:233
SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr sdfToExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue=0, int nIter=1)
Computes the extension of a scalar field, defined by the specified functor, off an iso-surface from a...
Index32 Index
Definition: openvdb/Types.h:50
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:333
math::Vec3< Real > Vec3R
Definition: openvdb/Types.h:68
Definition: openvdb/Exceptions.h:13
Definition: Coord.h:587
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
Private class of FastSweeping to perform multi-threaded initialization.
Definition: FastSweeping.h:759
DilateKernel & operator=(const DilateKernel &)=delete
const SdfValueT mBackground
Definition: FastSweeping.h:831
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:760
void run(int dilation, NearestNeighbors nn)
Definition: FastSweeping.h:768
DilateKernel(const DilateKernel &parent)=default
DilateKernel(FastSweeping &parent)
Definition: FastSweeping.h:761
FastSweeping * mParent
Definition: FastSweeping.h:830
Definition: FastSweeping.h:837
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:838
void operator()(const RootOrInternalNodeT &node) const
Definition: FastSweeping.h:926
void operator()(const LeafRange &r) const
Definition: FastSweeping.h:881
void run(SdfValueT isoValue, bool isInputSdf)
Definition: FastSweeping.h:844
SdfGridT * mSdfGrid
Definition: FastSweeping.h:937
SdfValueT mAboveSign
Definition: FastSweeping.h:939
SdfValueT mIsoValue
Definition: FastSweeping.h:938
InitSdf & operator=(const InitSdf &)=delete
InitSdf(FastSweeping &parent)
Definition: FastSweeping.h:839
FastSweeping * mParent
Definition: FastSweeping.h:936
SdfValueT mMin
Definition: FastSweeping.h:751
MinMaxKernel(MinMaxKernel &other, tbb::split)
Definition: FastSweeping.h:725
void operator()(const LeafRange &r)
Definition: FastSweeping.h:734
void join(const MinMaxKernel &other)
Definition: FastSweeping.h:745
typename LeafMgr::LeafRange LeafRange
Definition: FastSweeping.h:723
math::MinMax< SdfValueT > run(const SdfGridT &grid)
Definition: FastSweeping.h:727
SdfValueT mMax
Definition: FastSweeping.h:751
MinMaxKernel()
Definition: FastSweeping.h:724
SdfValueT v
Definition: FastSweeping.h:1311
Coord operator()(const Coord &p) const
Definition: FastSweeping.h:1316
bool operator<(const NN &rhs) const
Definition: FastSweeping.h:1317
NN(const SdfAccT &a, const Coord &p, int i)
Definition: FastSweeping.h:1315
static Coord ijk(const Coord &p, int i)
Definition: FastSweeping.h:1313
Private class of FastSweeping to perform concurrent fast sweeping in two directions.
Definition: FastSweeping.h:1193
ExtT threeNghbr(const NN &d1, const NN &d2, const NN &d3, const SdfT &, const ExtT &v1, const ExtT &v2, const ExtT &v3) const
Definition: FastSweeping.h:1331
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &w, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1327
void sweep()
Definition: FastSweeping.h:1343
SweepingKernel(const SweepingKernel &)=delete
SweepingKernel(FastSweeping &parent)
Definition: FastSweeping.h:1194
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1323
SweepingKernel & operator=(const SweepingKernel &)=delete
ExtT threeNghbr(const NN &d1, const NN &d2, const NN &d3, const SdfT &w, const ExtT &v1, const ExtT &v2, const ExtT &v3) const
Definition: FastSweeping.h:1339
void computeVoxelSlices(HashOp hash)
Main method that performs concurrent bi-directional sweeps.
Definition: FastSweeping.h:1200
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178