OpenVDB  8.1.1
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 
380  FastSweeping& operator=(const FastSweeping&) = delete;
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>
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
766  DilateKernel& operator=(const DilateKernel&) = delete;
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  dilateActiveValues(mgr, dilation, nn, IGNORE_TILES);
787 
788 #ifdef BENCHMARK_FAST_SWEEPING
789  timer.restart("Initializing grid and sweep mask");
790 #endif
791 
792  mParent->mSweepMask.clear();
793  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
794 
796  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
797  LeafManagerT leafManager(mParent->mSdfGrid->tree());
798 
799  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
800  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
801  const SdfValueT background = mBackground;//local copy
802  auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
803  assert(maskLeaf);
804  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
805  const SdfValueT value = *voxelIter;
806  if (math::Abs(value) < background) {// disable boundary voxels from the mask tree
807  maskLeaf->setValueOff(voxelIter.pos());
808  } else {
809  voxelIter.setValue(value > 0 ? Unknown : -Unknown);
810  }
811  }
812  };
813 
814  leafManager.foreach( kernel );
815 
816  // cache the leaf node origins for fast lookup in the sweeping kernels
817 
818  mParent->computeSweepMaskLeafOrigins();
819 
820 #ifdef BENCHMARK_FAST_SWEEPING
821  timer.stop();
822 #endif
823  }// FastSweeping::DilateKernel::run
824 
825  // Private member data of DilateKernel
827  const SdfValueT mBackground;
828 };// FastSweeping::DilateKernel
829 
831 template <typename SdfGridT, typename ExtValueT>
832 struct FastSweeping<SdfGridT, ExtValueT>::InitSdf
833 {
835  InitSdf(FastSweeping &parent): mParent(&parent),
836  mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
837  InitSdf(const InitSdf&) = default;// for tbb::parallel_for
838  InitSdf& operator=(const InitSdf&) = delete;
839 
840  void run(SdfValueT isoValue, bool isInputSdf)
841  {
842  mIsoValue = isoValue;
843  mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
844  SdfTreeT &tree = mSdfGrid->tree();//sdf
845  const bool hasActiveTiles = tree.hasActiveTiles();
846 
847  if (isInputSdf && hasActiveTiles) {
848  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
849  }
850 
851 #ifdef BENCHMARK_FAST_SWEEPING
852  util::CpuTimer timer("Initialize voxels");
853 #endif
854  mParent->mSweepMask.clear();
855  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
856 
857  {// Process all voxels
858  tree::LeafManager<SdfTreeT> mgr(tree, 1);// we need one auxiliary buffer
859  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
860  mgr.swapLeafBuffer(1);//swap voxel values
861  }
862 
863 #ifdef BENCHMARK_FAST_SWEEPING
864  timer.restart("Initialize tiles - new");
865 #endif
866  // Process all tiles
867  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree);
868  mgr.foreachBottomUp(*this);//multi-threaded
869  tree.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
870  if (hasActiveTiles) tree.voxelizeActiveTiles();//multi-threaded
871 
872  // cache the leaf node origins for fast lookup in the sweeping kernels
873 
874  mParent->computeSweepMaskLeafOrigins();
875  }// FastSweeping::InitSdf::run
876 
877  void operator()(const LeafRange& r) const
878  {
879  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
880  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
881  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
882  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
883  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
884  SdfValueT* sdf = leafIter.buffer(1).data();
885  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
886  const SdfValueT value = *voxelIter;
887  const bool isAbove = value > isoValue;
888  if (!voxelIter.isValueOn()) {// inactive voxels
889  sdf[voxelIter.pos()] = isAbove ? above : -above;
890  } else {// active voxels
891  const Coord ijk = voxelIter.getCoord();
892  stencil.moveTo(ijk, value);
893  const auto mask = stencil.intersectionMask( isoValue );
894  if (mask.none()) {// most common case
895  sdf[voxelIter.pos()] = isAbove ? above : -above;
896  } else {// compute distance to iso-surface
897  // disable boundary voxels from the mask tree
898  sweepMaskAcc.setValueOff(ijk);
899  const SdfValueT delta = value - isoValue;//offset relative to iso-value
900  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
901  sdf[voxelIter.pos()] = 0;
902  } else {//voxel is neighboring the iso-surface
903  SdfValueT sum = 0;
904  for (int i=0; i<6;) {
905  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
906  if (mask.test(i++)) d = math::Abs(delta/(value-stencil.getValue(i)));
907  if (mask.test(i++)) {
908  d2 = math::Abs(delta/(value-stencil.getValue(i)));
909  if (d2 < d) d = d2;
910  }
911  if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
912  }
913  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum) : -h / math::Sqrt(sum);
914  }// voxel is neighboring the iso-surface
915  }// intersecting voxels
916  }// active voxels
917  }// loop over voxels
918  }// loop over leaf nodes
919  }// FastSweeping::InitSdf::operator(const LeafRange&)
920 
921  template<typename RootOrInternalNodeT>
922  void operator()(const RootOrInternalNodeT& node) const
923  {
924  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
925  for (auto it = node.cbeginValueAll(); it; ++it) {
926  SdfValueT& v = const_cast<SdfValueT&>(*it);
927  v = v > isoValue ? above : -above;
928  }//loop over all tiles
929  }// FastSweeping::InitSdf::operator()(const RootOrInternalNodeT&)
930 
931  // Public member data
933  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
934  SdfValueT mIsoValue;
935  SdfValueT mAboveSign;//sign of distance values above the iso-value
936 };// FastSweeping::InitSdf
937 
938 
940 template <typename SdfGridT, typename ExtValueT>
941 template <typename OpT>
942 struct FastSweeping<SdfGridT, ExtValueT>::InitExt
943 {
945  using OpPoolT = tbb::enumerable_thread_specific<OpT>;
946  InitExt(FastSweeping &parent) : mParent(&parent),
947  mOpPool(nullptr), mSdfGrid(parent.mSdfGrid.get()),
948  mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
949  InitExt(const InitExt&) = default;// for tbb::parallel_for
950  InitExt& operator=(const InitExt&) = delete;
951  void run(SdfValueT isoValue, const OpT &opPrototype, bool isInputSdf)
952  {
953  static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::value, "Invalid return type of functor");
954  if (!mExtGrid) {
955  OPENVDB_THROW(RuntimeError, "FastSweeping::InitExt expected an extension grid!");
956  }
957 
958  mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
959  mIsoValue = isoValue;
960  auto &tree1 = mSdfGrid->tree();
961  auto &tree2 = mExtGrid->tree();
962  const bool hasActiveTiles = tree1.hasActiveTiles();//very fast
963 
964  if (isInputSdf && hasActiveTiles) {
965  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
966  }
967 
968 #ifdef BENCHMARK_FAST_SWEEPING
969  util::CpuTimer timer("Initialize voxels");
970 #endif
971 
972  mParent->mSweepMask.clear();
973  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
974 
975  {// Process all voxels
976  // Define thread-local operators
977  OpPoolT opPool(opPrototype);
978  mOpPool = &opPool;
979 
980  tree::LeafManager<SdfTreeT> mgr(tree1, 1);// we need one auxiliary buffer
981  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
982  mgr.swapLeafBuffer(1);//swap out auxiliary buffer
983  }
984 
985 #ifdef BENCHMARK_FAST_SWEEPING
986  timer.restart("Initialize tiles");
987 #endif
988  {// Process all tiles
989  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
990  mgr.foreachBottomUp(*this);//multi-threaded
991  tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
992  if (hasActiveTiles) {
993 #ifdef BENCHMARK_FAST_SWEEPING
994  timer.restart("Voxelizing active tiles");
995 #endif
996  tree1.voxelizeActiveTiles();//multi-threaded
997  tree2.voxelizeActiveTiles();//multi-threaded
998  }
999  }
1000 
1001  // cache the leaf node origins for fast lookup in the sweeping kernels
1002 
1003  mParent->computeSweepMaskLeafOrigins();
1004 
1005 #ifdef BENCHMARK_FAST_SWEEPING
1006  timer.stop();
1007 #endif
1008  }// FastSweeping::InitExt::run
1009 
1010  // int implements an update if minD needs to be updated
1011  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1012  void sumHelper(ExtT& sum2, ExtT ext, bool update, const SdfT& /* d2 */) const { if (update) sum2 = ext; }// int implementation
1013 
1014  // non-int implements a weighted sum
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 { sum2 += static_cast<ExtValueT>(d2 * ext); }// non-int implementation
1017 
1018  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1019  ExtT extValHelper(ExtT extSum, const SdfT& /* sdfSum */) const { return extSum; }// int implementation
1020 
1021  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1022  ExtT extValHelper(ExtT extSum, const SdfT& sdfSum) const {return ExtT((SdfT(1) / sdfSum) * extSum); }// non-int implementation
1023 
1024  void operator()(const LeafRange& r) const
1025  {
1026  ExtAccT acc(mExtGrid->tree());
1027  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1028  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1029  const math::Transform& xform = mExtGrid->transform();
1030  typename OpPoolT::reference op = mOpPool->local();
1031  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1032  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1033  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1034  SdfValueT *sdf = leafIter.buffer(1).data();
1035  ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();//should be safe!
1036  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1037  const SdfValueT value = *voxelIter;
1038  const bool isAbove = value > isoValue;
1039  if (!voxelIter.isValueOn()) {// inactive voxels
1040  sdf[voxelIter.pos()] = isAbove ? above : -above;
1041  } else {// active voxels
1042  const Coord ijk = voxelIter.getCoord();
1043  stencil.moveTo(ijk, value);
1044  const auto mask = stencil.intersectionMask( isoValue );
1045  if (mask.none()) {// no zero-crossing neighbors, most common case
1046  sdf[voxelIter.pos()] = isAbove ? above : -above;
1047  // the ext grid already has its active values set to the bakground value
1048  } else {// compute distance to iso-surface
1049  // disable boundary voxels from the mask tree
1050  sweepMaskAcc.setValueOff(ijk);
1051  const SdfValueT delta = value - isoValue;//offset relative to iso-value
1052  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1053  sdf[voxelIter.pos()] = 0;
1054  ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1055  } else {//voxel is neighboring the iso-surface
1056  SdfValueT sum1 = 0;
1057  ExtValueT sum2 = zeroVal<ExtValueT>();
1058  // minD is used to update sum2 in the integer case,
1059  // where we choose the value of the extension grid corresponding to
1060  // the smallest value of d in the 6 direction neighboring the current
1061  // voxel.
1062  SdfValueT minD = std::numeric_limits<SdfValueT>::max();
1063  for (int n=0, i=0; i<6;) {
1064  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1065  if (mask.test(i++)) {
1066  d = math::Abs(delta/(value-stencil.getValue(i)));
1067  n = i - 1;
1068  }
1069  if (mask.test(i++)) {
1070  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1071  if (d2 < d) {
1072  d = d2;
1073  n = i - 1;
1074  }
1075  }
1076  if (d < std::numeric_limits<SdfValueT>::max()) {
1077  d2 = 1/(d*d);
1078  sum1 += d2;
1079  const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1080  static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1081  static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1082  // If current d is less than minD, update minD
1083  sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1084  if (d < minD) minD = d;
1085  }
1086  }//look over six cases
1087  ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1088  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum1) : -h / math::Sqrt(sum1);
1089  }// voxel is neighboring the iso-surface
1090  }// intersecting voxels
1091  }// active voxels
1092  }// loop over voxels
1093  }// loop over leaf nodes
1094  }// FastSweeping::InitExt::operator(const LeafRange& r)
1095 
1096  template<typename RootOrInternalNodeT>
1097  void operator()(const RootOrInternalNodeT& node) const
1098  {
1099  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1100  for (auto it = node.cbeginValueAll(); it; ++it) {
1101  SdfValueT& v = const_cast<SdfValueT&>(*it);
1102  v = v > isoValue ? above : -above;
1103  }//loop over all tiles
1104  }
1105  // Public member data
1106  FastSweeping *mParent;
1107  OpPoolT *mOpPool;
1108  SdfGridT *mSdfGrid;
1109  ExtGridT *mExtGrid;
1110  SdfValueT mIsoValue;
1111  SdfValueT mAboveSign;//sign of distance values above the iso-value
1112 };// FastSweeping::InitExt
1113 
1115 template <typename SdfGridT, typename ExtValueT>
1116 template <typename MaskTreeT>
1117 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1118 {
1120  MaskKernel(FastSweeping &parent) : mParent(&parent),
1121  mSdfGrid(parent.mSdfGrid.get()) {}
1122  MaskKernel(const MaskKernel &parent) = default;// for tbb::parallel_for
1123  MaskKernel& operator=(const MaskKernel&) = delete;
1124 
1125  void run(const MaskTreeT &mask)
1126  {
1127 #ifdef BENCHMARK_FAST_SWEEPING
1128  util::CpuTimer timer;
1129 #endif
1130  auto &lsTree = mSdfGrid->tree();
1131 
1132  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1133 
1134 #ifdef BENCHMARK_FAST_SWEEPING
1135  timer.restart("Changing background value");
1136 #endif
1137  changeLevelSetBackground(lsTree, Unknown);//multi-threaded
1138 
1139 #ifdef BENCHMARK_FAST_SWEEPING
1140  timer.restart("Union with mask");//multi-threaded
1141 #endif
1142  lsTree.topologyUnion(mask);//multi-threaded
1143 
1144  // ignore active tiles since the input grid is assumed to be a level set
1145  tree::LeafManager<const MaskTreeT> mgr(mask);// super fast
1146 
1147 #ifdef BENCHMARK_FAST_SWEEPING
1148  timer.restart("Initializing grid and sweep mask");
1149 #endif
1150 
1151  mParent->mSweepMask.clear();
1152  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1153 
1154  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1155  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1156  LeafManagerT leafManager(mParent->mSweepMask);
1157 
1158  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
1159  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1160  SdfAccT acc(mSdfGrid->tree());
1161  // The following hack is safe due to the topoloyUnion in
1162  // init and the fact that SdfValueT is known to be a floating point!
1163  SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1164  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {// mask voxels
1165  if (math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1166  // disable boundary voxels from the mask tree
1167  voxelIter.setValue(false);
1168  }
1169  }
1170  };
1171  leafManager.foreach( kernel );
1172 
1173  // cache the leaf node origins for fast lookup in the sweeping kernels
1174  mParent->computeSweepMaskLeafOrigins();
1175 
1176 #ifdef BENCHMARK_FAST_SWEEPING
1177  timer.stop();
1178 #endif
1179  }// FastSweeping::MaskKernel::run
1180 
1181  // Private member data of MaskKernel
1182  FastSweeping *mParent;
1183  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1184 };// FastSweeping::MaskKernel
1185 
1187 template <typename SdfGridT, typename ExtValueT>
1188 struct FastSweeping<SdfGridT, ExtValueT>::SweepingKernel
1189 {
1190  SweepingKernel(FastSweeping &parent) : mParent(&parent) {}
1191  SweepingKernel(const SweepingKernel&) = delete;
1192  SweepingKernel& operator=(const SweepingKernel&) = delete;
1193 
1195  template<typename HashOp>
1196  void computeVoxelSlices(HashOp hash)
1197  {
1198 #ifdef BENCHMARK_FAST_SWEEPING
1199  util::CpuTimer timer;
1200 #endif
1201 
1202  // mask of the active voxels to be solved for, i.e. excluding boundary voxels
1203  const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1204 
1205  using LeafManagerT = typename tree::LeafManager<const SweepMaskTreeT>;
1206  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1207  LeafManagerT leafManager(maskTree);
1208 
1209  // compute the leaf node slices that have active voxels in them
1210  // the sliding window of the has keys is -14 to 21 (based on an 8x8x8 leaf node
1211  // and the extrema hash values i-j-k and i+j+k), but we use a larger mask window here to
1212  // easily accomodate any leaf dimension. The mask offset is used to be able to
1213  // store this in a fixed-size byte array
1214  constexpr int maskOffset = LeafT::DIM * 3;
1215  constexpr int maskRange = maskOffset * 2;
1216 
1217  // mark each possible slice in each leaf node that has one or more active voxels in it
1218  std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1219  auto kernel1 = [&](const LeafT& leaf, size_t leafIdx) {
1220  const size_t leafOffset = leafIdx * maskRange;
1221  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1222  const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1223  leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1224  }
1225  };
1226  leafManager.foreach( kernel1 );
1227 
1228  // compute the voxel slice map using a thread-local-storage hash map
1229  // the key of the hash map is the slice index of the voxel coord (ijk.x() + ijk.y() + ijk.z())
1230  // the values are an array of indices for every leaf that has active voxels with this slice index
1231  using ThreadLocalMap = std::unordered_map</*voxelSliceKey=*/int64_t, /*leafIdx=*/std::deque<size_t>>;
1232  tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1233  auto kernel2 = [&](const LeafT& leaf, size_t leafIdx) {
1234  ThreadLocalMap& map = pool.local();
1235  const Coord& origin = leaf.origin();
1236  const int64_t leafKey = hash(origin);
1237  const size_t leafOffset = leafIdx * maskRange;
1238  for (int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1239  if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1240  const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1241  map[voxelSliceKey].emplace_back(leafIdx);
1242  }
1243  }
1244  };
1245  leafManager.foreach( kernel2 );
1246 
1247  // combine into a single ordered map keyed by the voxel slice key
1248  // note that this is now stored in a map ordered by voxel slice key,
1249  // so sweep slices can be processed in order
1250  for (auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1251  const ThreadLocalMap& map = *poolIt;
1252  for (const auto& it : map) {
1253  for (const size_t leafIdx : it.second) {
1254  mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1255  }
1256  }
1257  }
1258 
1259  // extract the voxel slice keys for random access into the map
1260  mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1261  for (const auto& it : mVoxelSliceMap) {
1262  mVoxelSliceKeys.push_back(it.first);
1263  }
1264 
1265  // allocate the node masks in parallel, as the map is populated in serial
1266  auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1267  for (size_t i = range.begin(); i < range.end(); i++) {
1268  const int64_t key = mVoxelSliceKeys[i];
1269  for (auto& it : mVoxelSliceMap[key]) {
1270  it.second = std::make_unique<NodeMaskT>();
1271  }
1272  }
1273  };
1274  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1275 
1276  // each voxel slice contains a leafIdx-nodeMask pair,
1277  // this routine populates these node masks to select only the active voxels
1278  // from the mask tree that have the same voxel slice key
1279  // TODO: a small optimization here would be to union this leaf node mask with
1280  // a pre-computed one for this particular slice pattern
1281  auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1282  for (size_t i = range.begin(); i < range.end(); i++) {
1283  const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1284  LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1285  for (LeafSlice& leafSlice : leafSliceArray) {
1286  const size_t leafIdx = leafSlice.first;
1287  NodeMaskPtrT& nodeMask = leafSlice.second;
1288  const LeafT& leaf = leafManager.leaf(leafIdx);
1289  const Coord& origin = leaf.origin();
1290  const int64_t leafKey = hash(origin);
1291  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1292  const Index voxelIdx = voxelIter.pos();
1293  const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1294  const int64_t key = leafKey + hash(ijk);
1295  if (key == voxelSliceKey) {
1296  nodeMask->setOn(voxelIdx);
1297  }
1298  }
1299  }
1300  }
1301  };
1302  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1303  }// FastSweeping::SweepingKernel::computeVoxelSlices
1304 
1305  // Private struct for nearest neighbor grid points (very memory light!)
1306  struct NN {
1307  SdfValueT v;
1308  int n;
1309  inline static Coord ijk(const Coord &p, int i) { return p + FastSweeping::mOffset[i]; }
1310  NN() : v(), n() {}
1311  NN(const SdfAccT &a, const Coord &p, int i) : v(math::Abs(a.getValue(ijk(p,i)))), n(i) {}
1312  inline Coord operator()(const Coord &p) const { return ijk(p, n); }
1313  inline bool operator<(const NN &rhs) const { return v < rhs.v; }
1314  inline operator bool() const { return v < SdfValueT(1000); }
1315  };// NN
1316 
1318  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1319  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
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 ExtT(w*(d1.v*v1 + d2.v*v2)); }// non-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 threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& /* w */, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1328  math::Vec3<SdfT> d(d1.v, d2.v, d3.v);
1329  math::Vec3<ExtT> v(v1, v2, v3);
1330  return v[static_cast<int>(math::MinIndex(d))];
1331  }// int implementation
1332 
1334  template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1335  ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& w, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1336  return ExtT(w*(d1.v*v1 + d2.v*v2 + d3.v*v3));
1337  }// non-int implementation
1338 
1339  void sweep()
1340  {
1341  typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() : nullptr;
1342 
1343  const SdfValueT h = static_cast<SdfValueT>(mParent->mSdfGrid->voxelSize()[0]);
1344  const SdfValueT sqrt2h = math::Sqrt(SdfValueT(2))*h;
1345 
1346  const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1347 
1348  int64_t voxelSliceIndex(0);
1349 
1350  auto kernel = [&](const tbb::blocked_range<size_t>& range) {
1351  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
1352 
1353  SdfAccT acc1(mParent->mSdfGrid->tree());
1354  auto acc2 = std::unique_ptr<ExtAccT>(tree2 ? new ExtAccT(*tree2) : nullptr);
1355  SdfValueT absV, sign, update, D;
1356  NN d1, d2, d3;//distance values and coordinates of closest neighbor points
1357 
1358  const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1359 
1360  // Solves Goudonov's scheme: [x-d1]^2 + [x-d2]^2 + [x-d3]^2 = h^2
1361  // where [X] = (X>0?X:0) and ai=min(di+1,di-1)
1362  for (size_t i = range.begin(); i < range.end(); ++i) {
1363 
1364  // iterate over all leafs in the slice and extract the leaf
1365  // and node mask for each slice pattern
1366 
1367  const LeafSlice& leafSlice = leafSliceArray[i];
1368  const size_t leafIdx = leafSlice.first;
1369  const NodeMaskPtrT& nodeMask = leafSlice.second;
1370 
1371  const Coord& origin = leafNodeOrigins[leafIdx];
1372 
1373  Coord ijk;
1374  for (auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1375 
1376  // Get coordinate of center point of the FD stencil
1377  ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1378 
1379  // Find the closes neighbors in the three axial directions
1380  d1 = std::min(NN(acc1, ijk, 0), NN(acc1, ijk, 1));
1381  d2 = std::min(NN(acc1, ijk, 2), NN(acc1, ijk, 3));
1382  d3 = std::min(NN(acc1, ijk, 4), NN(acc1, ijk, 5));
1383 
1384  if (!(d1 || d2 || d3)) continue;//no valid neighbors
1385 
1386  // Get the center point of the FD stencil (assumed to be an active voxel)
1387  // Note this const_cast is normally unsafe but by design we know the tree
1388  // to be static, of floating-point type and containing active voxels only!
1389  SdfValueT &value = const_cast<SdfValueT&>(acc1.getValue(ijk));
1390 
1391  // Extract the sign
1392  sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1393 
1394  // Absolute value
1395  absV = math::Abs(value);
1396 
1397  // sort values so d1 <= d2 <= d3
1398  if (d2 < d1) std::swap(d1, d2);
1399  if (d3 < d2) std::swap(d2, d3);
1400  if (d2 < d1) std::swap(d1, d2);
1401 
1402  // Test if there is a solution depending on ONE of the neighboring voxels
1403  // if d2 - d1 >= h => d2 >= d1 + h then:
1404  // (x-d1)^2=h^2 => x = d1 + h
1405  update = d1.v + h;
1406  if (update <= d2.v) {
1407  if (update < absV) {
1408  value = sign * update;
1409  if (acc2) acc2->setValue(ijk, acc2->getValue(d1(ijk)));//update ext?
1410  }//update sdf?
1411  continue;
1412  }// one neighbor case
1413 
1414  // Test if there is a solution depending on TWO of the neighboring voxels
1415  // (x-d1)^2 + (x-d2)^2 = h^2
1416  //D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1417  //if (D >= SdfValueT(0)) {// non-negative discriminant
1418  if (d2.v <= sqrt2h + d1.v) {
1419  D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1420  update = SdfValueT(0.5) * (d1.v + d2.v + std::sqrt(D));
1421  if (update > d2.v && update <= d3.v) {
1422  if (update < absV) {
1423  value = sign * update;
1424  if (acc2) {
1425  d1.v -= update;
1426  d2.v -= update;
1427  // affine combination of two neighboring extension values
1428  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v);
1429  const ExtValueT v1 = acc2->getValue(d1(ijk));
1430  const ExtValueT v2 = acc2->getValue(d2(ijk));
1431  const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1432  acc2->setValue(ijk, extVal);
1433  }//update ext?
1434  }//update sdf?
1435  continue;
1436  }//test for two neighbor case
1437  }//test for non-negative determinant
1438 
1439  // Test if there is a solution depending on THREE of the neighboring voxels
1440  // (x-d1)^2 + (x-d2)^2 + (x-d3)^2 = h^2
1441  // 3x^2 - 2(d1 + d2 + d3)x + d1^2 + d2^2 + d3^2 = h^2
1442  // ax^2 + bx + c=0, a=3, b=-2(d1+d2+d3), c=d1^2 + d2^2 + d3^2 - h^2
1443  const SdfValueT d123 = d1.v + d2.v + d3.v;
1444  D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.v*d2.v + d3.v*d3.v - h * h);
1445  if (D >= SdfValueT(0)) {// non-negative discriminant
1446  update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));//always passes test
1447  //if (update > d3.v) {//disabled due to round-off errors
1448  if (update < absV) {
1449  value = sign * update;
1450  if (acc2) {
1451  d1.v -= update;
1452  d2.v -= update;
1453  d3.v -= update;
1454  // affine combination of three neighboring extension values
1455  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v+d3.v);
1456  const ExtValueT v1 = acc2->getValue(d1(ijk));
1457  const ExtValueT v2 = acc2->getValue(d2(ijk));
1458  const ExtValueT v3 = acc2->getValue(d3(ijk));
1459  const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1460  acc2->setValue(ijk, extVal);
1461  }//update ext?
1462  }//update sdf?
1463  }//test for non-negative determinant
1464  }//loop over coordinates
1465  }
1466  };
1467 
1468 #ifdef BENCHMARK_FAST_SWEEPING
1469  util::CpuTimer timer("Forward sweep");
1470 #endif
1471 
1472  for (size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1473  voxelSliceIndex = mVoxelSliceKeys[i];
1474  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1475  }
1476 
1477 #ifdef BENCHMARK_FAST_SWEEPING
1478  timer.restart("Backward sweeps");
1479 #endif
1480  for (size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1481  voxelSliceIndex = mVoxelSliceKeys[i-1];
1482  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1483  }
1484 
1485 #ifdef BENCHMARK_FAST_SWEEPING
1486  timer.stop();
1487 #endif
1488  }// FastSweeping::SweepingKernel::sweep
1489 
1490 private:
1491  using NodeMaskT = typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1492  using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1493  // using a unique ptr for the NodeMask allows for parallel allocation,
1494  // but makes this class not copy-constructible
1495  using LeafSlice = std::pair</*leafIdx=*/size_t, /*leafMask=*/NodeMaskPtrT>;
1496  using LeafSliceArray = std::deque<LeafSlice>;
1497  using VoxelSliceMap = std::map</*voxelSliceKey=*/int64_t, LeafSliceArray>;
1498 
1499  // Private member data of SweepingKernel
1500  FastSweeping *mParent;
1501  VoxelSliceMap mVoxelSliceMap;
1502  std::vector<int64_t> mVoxelSliceKeys;
1503 };// FastSweeping::SweepingKernel
1504 
1506 
1507 template<typename GridT>
1508 typename GridT::Ptr
1509 fogToSdf(const GridT &fogGrid,
1510  typename GridT::ValueType isoValue,
1511  int nIter)
1512 {
1514  if (fs.initSdf(fogGrid, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1515  return fs.sdfGrid();
1516 }
1517 
1518 template<typename GridT>
1519 typename GridT::Ptr
1520 sdfToSdf(const GridT &sdfGrid,
1521  typename GridT::ValueType isoValue,
1522  int nIter)
1523 {
1525  if (fs.initSdf(sdfGrid, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1526  return fs.sdfGrid();
1527 }
1528 
1529 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1530 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1531 fogToExt(const FogGridT &fogGrid,
1532  const ExtOpT &op,
1533  const ExtValueT& background,
1534  typename FogGridT::ValueType isoValue,
1535  int nIter)
1536 {
1538  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1539  return fs.extGrid();
1540 }
1541 
1542 template<typename SdfGridT, typename OpT, typename ExtValueT>
1543 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1544 sdfToExt(const SdfGridT &sdfGrid,
1545  const OpT &op,
1546  const ExtValueT &background,
1547  typename SdfGridT::ValueType isoValue,
1548  int nIter)
1549 {
1551  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1552  return fs.extGrid();
1553 }
1554 
1555 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1556 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1557 fogToSdfAndExt(const FogGridT &fogGrid,
1558  const ExtOpT &op,
1559  const ExtValueT &background,
1560  typename FogGridT::ValueType isoValue,
1561  int nIter)
1562 {
1564  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1565  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1566 }
1567 
1568 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
1569 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1570 sdfToSdfAndExt(const SdfGridT &sdfGrid,
1571  const ExtOpT &op,
1572  const ExtValueT &background,
1573  typename SdfGridT::ValueType isoValue,
1574  int nIter)
1575 {
1577  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1578  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1579 }
1580 
1581 template<typename GridT>
1582 typename GridT::Ptr
1583 dilateSdf(const GridT &sdfGrid,
1584  int dilation,
1585  NearestNeighbors nn,
1586  int nIter)
1587 {
1589  if (fs.initDilate(sdfGrid, dilation, nn)) fs.sweep(nIter);
1590  return fs.sdfGrid();
1591 }
1592 
1593 template<typename GridT, typename MaskTreeT>
1594 typename GridT::Ptr
1595 maskSdf(const GridT &sdfGrid,
1596  const Grid<MaskTreeT> &mask,
1597  bool ignoreActiveTiles,
1598  int nIter)
1599 {
1601  if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1602  return fs.sdfGrid();
1603 }
1604 
1605 } // namespace tools
1606 } // namespace OPENVDB_VERSION_NAME
1607 } // namespace openvdb
1608 
1609 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
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:1520
Definition: ValueAccessor.h:182
SdfGridT * mSdfGrid
Definition: FastSweeping.h:933
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:1335
math::MinMax< SdfValueT > run(const SdfGridT &grid)
Definition: FastSweeping.h:727
static Coord ijk(const Coord &p, int i)
Definition: FastSweeping.h:1309
Definition: Morphology.h:56
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:1327
SweepingKernel(FastSweeping &parent)
Definition: FastSweeping.h:1190
double restart()
Re-start timer.
Definition: CpuTimer.h:150
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
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
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
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
void clear()
Clears all the grids and counters so initializtion can be called again.
Definition: FastSweeping.h:556
SdfValueT mAboveSign
Definition: FastSweeping.h:935
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:834
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:1557
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:1531
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
void join(const MinMaxKernel &other)
Definition: FastSweeping.h:745
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:1544
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean, variance, etc.) of grid values.
Definition: Coord.h:587
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays, one for each level of the tree.
Definition: NodeManager.h:30
Definition: Stencils.h:1231
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Templated class to compute the minimum and maximum values.
Definition: Stats.h:30
SdfValueT v
Definition: FastSweeping.h:1307
typename LeafMgr::LeafRange LeafRange
Definition: FastSweeping.h:723
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:108
MinMaxKernel(MinMaxKernel &other, tbb::split)
Definition: FastSweeping.h:725
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
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1319
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Index32 Index
Definition: openvdb/Types.h:50
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:1509
void operator()(const LeafRange &r)
Definition: FastSweeping.h:734
MinMaxKernel()
Definition: FastSweeping.h:724
void computeVoxelSlices(HashOp hash)
Main method that performs concurrent bi-directional sweeps.
Definition: FastSweeping.h:1196
Simple timer for basic profiling.
Definition: CpuTimer.h:66
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:938
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
bool isValid() const
Return true if there are voxels and boundaries to solve for.
Definition: FastSweeping.h:513
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
void run(int dilation, NearestNeighbors nn)
Definition: FastSweeping.h:768
Vec3< double > Vec3d
Definition: Vec3.h:668
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
void sweep()
Definition: FastSweeping.h:1339
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
size_t boundaryVoxelCount() const
Return the number of voxels that defined the boundary condition.
Definition: FastSweeping.h:510
SdfGridT::Ptr sdfGrid()
Returns a shared pointer to the signed distance field computed by this class.
Definition: FastSweeping.h:388
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:1570
Definition: openvdb/Exceptions.h:63
FastSweeping & operator=(const FastSweeping &)=delete
Disallow copy assignment.
Implementation of morphological dilation and erosion.
const SdfValueT mBackground
Definition: FastSweeping.h:827
Definition: openvdb/Exceptions.h:13
NN(const SdfAccT &a, const Coord &p, int i)
Definition: FastSweeping.h:1311
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
Private class of FastSweeping to perform concurrent fast sweeping in two directions.
Definition: FastSweeping.h:1188
Definition: Morphology.h:78
math::Transform & transform()
Return a reference to this grid&#39;s transform, which might be shared with other grids.
Definition: Grid.h:415
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
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
void operator()(const LeafRange &r) const
Definition: FastSweeping.h:877
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &w, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1323
Definition: LeafManager.h:101
void sweep(int nIter=1, bool finalize=true)
Perform nIter iterations of the fast sweeping algorithm.
Definition: FastSweeping.h:659
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:760
NearestNeighbors
Voxel topology of nearest neighbors.
Definition: Morphology.h:56
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
Definition: ValueAccessor.h:379
~FastSweeping()
Destructor.
Definition: FastSweeping.h:374
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3705
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
SdfValueT mMin
Definition: FastSweeping.h:751
Definition: openvdb/Types.h:333
ExtGridT::Ptr extGrid()
Returns a shared pointer to the extension field computed by this class.
Definition: FastSweeping.h:396
FastSweeping * mParent
Definition: FastSweeping.h:826
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...
Computes signed distance values from an initial iso-surface and optionally performs velocty extension...
Definition: FastSweeping.h:350
TreeType & tree()
Return a reference to this grid&#39;s tree, which might be shared with other grids.
Definition: Grid.h:917
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
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...
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
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
Coord operator()(const Coord &p) const
Definition: FastSweeping.h:1312
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
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node&#39;s buffer with the nth corresponding auxiliary buffer, where n = bufferIdx...
Definition: LeafManager.h:359
SharedPtr< Grid > Ptr
Definition: Grid.h:579
void run(SdfValueT isoValue, bool isInputSdf)
Definition: FastSweeping.h:840
FastSweeping * mParent
Definition: FastSweeping.h:932
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
Private class of FastSweeping to perform multi-threaded initialization.
Definition: FastSweeping.h:758
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
DilateKernel(FastSweeping &parent)
Definition: FastSweeping.h:761
bool operator<(const NN &rhs) const
Definition: FastSweeping.h:1313
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
SdfValueT mIsoValue
Definition: FastSweeping.h:934
SdfValueT mMax
Definition: FastSweeping.h:751
size_t sweepingVoxelCount() const
Return the number of voxels that will be solved for.
Definition: FastSweeping.h:507
bool initDilate(const SdfGridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE)
Initializer used when dilating an exsiting signed distance field.
Definition: FastSweeping.h:617
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:1595
InitSdf(FastSweeping &parent)
Definition: FastSweeping.h:835
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
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:1583
Definition: Transform.h:39
Definition: FastSweeping.h:832
void operator()(const RootOrInternalNodeT &node) const
Definition: FastSweeping.h:922