27 #ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
28 #define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
32 #include <type_traits>
36 #include <unordered_map>
39 #include <tbb/parallel_for.h>
40 #include <tbb/enumerable_thread_specific.h>
41 #include <tbb/task_group.h>
50 #ifdef BENCHMARK_FAST_SWEEPING
87 template<
typename Gr
idT>
90 typename GridT::ValueType isoValue,
120 template<
typename Gr
idT>
123 typename GridT::ValueType isoValue = 0,
158 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
159 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
162 const ExtValueT& background,
163 typename FogGridT::ValueType isoValue,
196 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
197 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
200 const ExtValueT &background,
201 typename SdfGridT::ValueType isoValue = 0,
238 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
239 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
242 const ExtValueT &background,
243 typename FogGridT::ValueType isoValue,
280 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
281 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
284 const ExtValueT &background,
285 typename SdfGridT::ValueType isoValue = 0,
305 template<
typename Gr
idT>
330 template<
typename Gr
idT,
typename MaskTreeT>
334 bool ignoreActiveTiles =
false,
349 template<
typename SdfGr
idT,
typename ExtValueT =
typename SdfGr
idT::ValueType>
352 static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
353 "FastSweeping requires SdfGridT to have floating-point values");
355 using SdfValueT =
typename SdfGridT::ValueType;
356 using SdfTreeT =
typename SdfGridT::TreeType;
360 using ExtGridT =
typename SdfGridT::template ValueConverter<ExtValueT>::Type;
361 using ExtTreeT =
typename ExtGridT::TreeType;
365 using SweepMaskTreeT =
typename SdfTreeT::template ValueConverter<ValueMask>::Type;
388 typename SdfGridT::Ptr
sdfGrid() {
return mSdfGrid; }
396 typename ExtGridT::Ptr
extGrid() {
return mExtGrid; }
418 bool initSdf(
const SdfGridT &sdfGrid, SdfValueT isoValue,
bool isInputSdf);
448 template <
typename ExtOpT>
449 bool initExt(
const SdfGridT &sdfGrid,
const ExtOpT &op,
const ExtValueT &background, SdfValueT isoValue,
bool isInputSdf);
486 template<
typename MaskTreeT>
487 bool initMask(
const SdfGridT &sdfGrid,
const Grid<MaskTreeT> &mask,
bool ignoreActiveTiles =
false);
501 void sweep(
int nIter = 1,
bool finalize =
true);
513 bool isValid()
const {
return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
518 void computeSweepMaskLeafOrigins();
528 struct SweepingKernel;
531 static const Coord mOffset[6];
534 typename SdfGridT::Ptr mSdfGrid;
535 typename ExtGridT::Ptr mExtGrid;
536 SweepMaskTreeT mSweepMask;
537 std::vector<Coord> mSweepMaskLeafOrigins;
538 size_t mSweepingVoxelCount, mBoundaryVoxelCount;
544 template <
typename SdfGr
idT,
typename ExtValueT>
545 const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
549 template <
typename SdfGr
idT,
typename ExtValueT>
551 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0)
555 template <
typename SdfGr
idT,
typename ExtValueT>
561 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
564 template <
typename SdfGr
idT,
typename ExtValueT>
570 mSweepMask.voxelizeActiveTiles();
573 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
574 LeafManagerT leafManager(mSweepMask);
576 mSweepMaskLeafOrigins.resize(leafManager.leafCount());
577 tbb::atomic<size_t> sweepingVoxelCount = 0;
578 auto kernel = [&](
const LeafT& leaf,
size_t leafIdx) {
579 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
580 sweepingVoxelCount += leaf.onVoxelCount();
582 leafManager.foreach(kernel,
true, 1024);
584 mBoundaryVoxelCount = 0;
585 mSweepingVoxelCount = sweepingVoxelCount;
587 const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
588 assert( totalCount >= mSweepingVoxelCount );
589 mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
593 template <
typename SdfGr
idT,
typename ExtValueT>
597 mSdfGrid = fogGrid.deepCopy();
599 kernel.
run(isoValue, isInputSdf);
600 return this->isValid();
603 template <
typename SdfGr
idT,
typename ExtValueT>
604 template <
typename OpT>
608 mSdfGrid = fogGrid.deepCopy();
609 mExtGrid = createGrid<ExtGridT>( background );
610 mExtGrid->topologyUnion( *mSdfGrid );
611 InitExt<OpT> kernel(*
this);
612 kernel.run(isoValue, op, isInputSdf);
613 return this->isValid();
616 template <
typename SdfGr
idT,
typename ExtValueT>
620 mSdfGrid = sdfGrid.deepCopy();
622 kernel.
run(dilate, nn);
623 return this->isValid();
626 template <
typename SdfGr
idT,
typename ExtValueT>
627 template<
typename MaskTreeT>
631 mSdfGrid = sdfGrid.deepCopy();
633 if (mSdfGrid->transform() != mask.transform()) {
638 using T =
typename MaskTreeT::template ValueConverter<bool>::Type;
640 tmp->
tree().voxelizeActiveTiles();
641 MaskKernel<T> kernel(*
this);
642 kernel.run(tmp->
tree());
644 if (ignoreActiveTiles || !mask.
tree().hasActiveTiles()) {
645 MaskKernel<MaskTreeT> kernel(*
this);
646 kernel.run(mask.
tree());
648 using T =
typename MaskTreeT::template ValueConverter<ValueMask>::Type;
650 tmp.voxelizeActiveTiles(
true);
651 MaskKernel<T> kernel(*
this);
655 return this->isValid();
658 template <
typename SdfGr
idT,
typename ExtValueT>
664 if (this->boundaryVoxelCount() == 0) {
666 }
else if (this->sweepingVoxelCount() == 0) {
671 std::deque<SweepingKernel> kernels;
672 for (
int i = 0; i < 4; i++) kernels.emplace_back(*
this);
675 #ifdef BENCHMARK_FAST_SWEEPING
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]; }); });
687 #ifdef BENCHMARK_FAST_SWEEPING
693 for (
int i = 0; i < nIter; ++i) {
698 #ifdef BENCHMARK_FAST_SWEEPING
702 auto e = kernel.
run(*mSdfGrid);
704 #ifdef BENCHMARK_FAST_SWEEPING
705 std::cerr <<
"Min = " << e.min() <<
" Max = " << e.max() << std::endl;
706 timer.
restart(
"Changing asymmetric background value");
710 #ifdef BENCHMARK_FAST_SWEEPING
719 template <
typename SdfGr
idT,
typename ExtValueT>
730 tbb::parallel_reduce(mgr.
leafRange(), *
this);
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;
747 if (other.
mMin < mMin) mMin = other.
mMin;
748 if (other.
mMax > mMax) mMax = other.
mMax;
757 template <
typename SdfGr
idT,
typename ExtValueT>
762 : mParent(&parent), mBackground(parent.mSdfGrid->background())
770 #ifdef BENCHMARK_FAST_SWEEPING
775 #ifdef BENCHMARK_FAST_SWEEPING
776 timer.
restart(
"Changing background value");
781 #ifdef BENCHMARK_FAST_SWEEPING
782 timer.
restart(
"Dilating and updating mgr (parallel)");
792 #ifdef BENCHMARK_FAST_SWEEPING
793 timer.
restart(
"Initializing grid and sweep mask");
796 mParent->mSweepMask.clear();
797 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
800 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
801 LeafManagerT leafManager(mParent->mSdfGrid->tree());
803 auto kernel = [&](LeafT& leaf,
size_t ) {
805 const SdfValueT background = mBackground;
806 auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
808 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
809 const SdfValueT value = *voxelIter;
811 maskLeaf->setValueOff(voxelIter.pos());
813 voxelIter.setValue(value > 0 ? Unknown : -Unknown);
818 leafManager.foreach( kernel );
822 mParent->computeSweepMaskLeafOrigins();
824 #ifdef BENCHMARK_FAST_SWEEPING
835 template <
typename SdfGr
idT,
typename ExtValueT>
840 mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
844 void run(SdfValueT isoValue,
bool isInputSdf)
846 mIsoValue = isoValue;
847 mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
848 SdfTreeT &tree = mSdfGrid->tree();
849 const bool hasActiveTiles = tree.hasActiveTiles();
851 if (isInputSdf && hasActiveTiles) {
855 #ifdef BENCHMARK_FAST_SWEEPING
858 mParent->mSweepMask.clear();
859 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
863 tbb::parallel_for(mgr.
leafRange(32), *
this);
867 #ifdef BENCHMARK_FAST_SWEEPING
868 timer.
restart(
"Initialize tiles - new");
872 mgr.foreachBottomUp(*
this);
874 if (hasActiveTiles) tree.voxelizeActiveTiles();
878 mParent->computeSweepMaskLeafOrigins();
886 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
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()) {
893 sdf[voxelIter.pos()] = isAbove ? above : -above;
895 const Coord ijk = voxelIter.getCoord();
896 stencil.
moveTo(ijk, value);
899 sdf[voxelIter.pos()] = isAbove ? above : -above;
903 const SdfValueT delta = value - isoValue;
905 sdf[voxelIter.pos()] = 0;
908 for (
int i=0; i<6;) {
911 if (mask.test(i++)) {
925 template<
typename RootOrInternalNodeT>
929 for (
auto it = node.cbeginValueAll(); it; ++it) {
930 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
931 v = v > isoValue ? above : -above;
943 template <
typename SdfGr
idT,
typename ExtValueT>
944 template <
typename OpT>
948 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
950 mOpPool(
nullptr), mSdfGrid(parent.mSdfGrid.get()),
951 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
952 InitExt(
const InitExt&) =
default;
953 InitExt&
operator=(
const InitExt&) =
delete;
954 void run(SdfValueT isoValue,
const OpT &opPrototype,
bool isInputSdf)
956 static_assert(std::is_convertible<decltype(opPrototype(
Vec3d(0))),ExtValueT>::value,
"Invalid return type of functor");
961 mAboveSign = isInputSdf ? SdfValueT(1) : SdfValueT(-1);
962 mIsoValue = isoValue;
963 auto &tree1 = mSdfGrid->tree();
964 auto &tree2 = mExtGrid->tree();
965 const bool hasActiveTiles = tree1.hasActiveTiles();
967 if (isInputSdf && hasActiveTiles) {
971 #ifdef BENCHMARK_FAST_SWEEPING
975 mParent->mSweepMask.clear();
976 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
980 OpPoolT opPool(opPrototype);
984 tbb::parallel_for(mgr.leafRange(32), *
this);
985 mgr.swapLeafBuffer(1);
988 #ifdef BENCHMARK_FAST_SWEEPING
989 timer.restart(
"Initialize tiles");
992 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
993 mgr.foreachBottomUp(*
this);
995 if (hasActiveTiles) {
996 #ifdef BENCHMARK_FAST_SWEEPING
997 timer.restart(
"Voxelizing active tiles");
999 tree1.voxelizeActiveTiles();
1000 tree2.voxelizeActiveTiles();
1006 mParent->computeSweepMaskLeafOrigins();
1008 #ifdef BENCHMARK_FAST_SWEEPING
1013 void operator()(
const LeafRange& r)
const
1015 ExtAccT acc(mExtGrid->tree());
1016 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1017 math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1018 const math::Transform& xform = mExtGrid->transform();
1019 typename OpPoolT::reference op = mOpPool->local();
1021 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1022 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1023 SdfValueT *sdf = leafIter.buffer(1).data();
1024 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();
1025 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1026 const SdfValueT value = *voxelIter;
1027 const bool isAbove = value > isoValue;
1028 if (!voxelIter.isValueOn()) {
1029 sdf[voxelIter.pos()] = isAbove ? above : -above;
1031 const Coord ijk = voxelIter.getCoord();
1032 stencil.moveTo(ijk, value);
1033 const auto mask = stencil.intersectionMask( isoValue );
1035 sdf[voxelIter.pos()] = isAbove ? above : -above;
1039 sweepMaskAcc.setValueOff(ijk);
1040 const SdfValueT delta = value - isoValue;
1042 sdf[voxelIter.pos()] = 0;
1043 ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1046 ExtValueT sum2 = zeroVal<ExtValueT>();
1047 for (
int n=0, i=0; i<6;) {
1049 if (mask.test(i++)) {
1050 d =
math::Abs(delta/(value-stencil.getValue(i)));
1053 if (mask.test(i++)) {
1054 d2 =
math::Abs(delta/(value-stencil.getValue(i)));
1063 const Vec3R xyz(
static_cast<SdfValueT
>(ijk[0])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][0]),
1064 static_cast<SdfValueT
>(ijk[1])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][1]),
1065 static_cast<SdfValueT
>(ijk[2])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][2]));
1066 sum2 += d2*ExtValueT(op(xform.indexToWorld(xyz)));
1069 ext[voxelIter.pos()] = (SdfValueT(1) / sum1) * sum2;
1070 sdf[voxelIter.pos()] = isAbove ? h /
math::Sqrt(sum1) : -h / math::
Sqrt(sum1);
1078 template<
typename RootOrInternalNodeT>
1079 void operator()(
const RootOrInternalNodeT& node)
const
1082 for (
auto it = node.cbeginValueAll(); it; ++it) {
1083 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1084 v = v > isoValue ? above : -above;
1092 SdfValueT mIsoValue;
1093 SdfValueT mAboveSign;
1097 template <
typename SdfGr
idT,
typename ExtValueT>
1098 template <
typename MaskTreeT>
1099 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1101 using LeafRange =
typename tree::LeafManager<const MaskTreeT>::LeafRange;
1103 mSdfGrid(parent.mSdfGrid.get()) {}
1104 MaskKernel(
const MaskKernel &parent) =
default;
1105 MaskKernel&
operator=(
const MaskKernel&) =
delete;
1107 void run(
const MaskTreeT &mask)
1109 #ifdef BENCHMARK_FAST_SWEEPING
1110 util::CpuTimer timer;
1112 auto &lsTree = mSdfGrid->tree();
1116 #ifdef BENCHMARK_FAST_SWEEPING
1117 timer.restart(
"Changing background value");
1121 #ifdef BENCHMARK_FAST_SWEEPING
1122 timer.restart(
"Union with mask");
1124 lsTree.topologyUnion(mask);
1127 tree::LeafManager<const MaskTreeT> mgr(mask);
1129 #ifdef BENCHMARK_FAST_SWEEPING
1130 timer.restart(
"Initializing grid and sweep mask");
1133 mParent->mSweepMask.clear();
1134 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1136 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1137 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1138 LeafManagerT leafManager(mParent->mSweepMask);
1140 auto kernel = [&](LeafT& leaf,
size_t ) {
1142 SdfAccT acc(mSdfGrid->tree());
1145 SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1146 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1147 if (
math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1149 voxelIter.setValue(
false);
1153 leafManager.foreach( kernel );
1156 mParent->computeSweepMaskLeafOrigins();
1158 #ifdef BENCHMARK_FAST_SWEEPING
1169 template <
typename SdfGr
idT,
typename ExtValueT>
1177 template<
typename HashOp>
1180 #ifdef BENCHMARK_FAST_SWEEPING
1185 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1188 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1189 LeafManagerT leafManager(maskTree);
1196 constexpr
int maskOffset = LeafT::DIM * 3;
1197 constexpr
int maskRange = maskOffset * 2;
1200 std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1201 auto kernel1 = [&](
const LeafT& leaf,
size_t leafIdx) {
1202 const size_t leafOffset = leafIdx * maskRange;
1203 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1204 const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1205 leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1208 leafManager.foreach( kernel1 );
1213 using ThreadLocalMap = std::unordered_map<int64_t, std::deque<size_t>>;
1214 tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1215 auto kernel2 = [&](
const LeafT& leaf,
size_t leafIdx) {
1216 ThreadLocalMap& map = pool.local();
1217 const Coord& origin = leaf.origin();
1218 const int64_t leafKey = hash(origin);
1219 const size_t leafOffset = leafIdx * maskRange;
1220 for (
int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1221 if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1222 const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1223 assert(voxelSliceKey >= 0);
1224 map[voxelSliceKey].emplace_back(leafIdx);
1228 leafManager.foreach( kernel2 );
1233 for (
auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1234 const ThreadLocalMap& map = *poolIt;
1235 for (
const auto& it : map) {
1236 for (
const size_t leafIdx : it.second) {
1237 mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1243 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1244 for (
const auto& it : mVoxelSliceMap) {
1245 mVoxelSliceKeys.push_back(it.first);
1249 auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1250 for (
size_t i = range.begin(); i < range.end(); i++) {
1251 const int64_t key = mVoxelSliceKeys[i];
1252 for (
auto& it : mVoxelSliceMap[key]) {
1253 it.second = std::make_unique<NodeMaskT>();
1257 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1264 auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1265 for (
size_t i = range.begin(); i < range.end(); i++) {
1266 const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1267 LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1268 for (LeafSlice& leafSlice : leafSliceArray) {
1269 const size_t leafIdx = leafSlice.first;
1270 NodeMaskPtrT& nodeMask = leafSlice.second;
1271 const LeafT& leaf = leafManager.leaf(leafIdx);
1272 const Coord& origin = leaf.origin();
1273 const int64_t leafKey = hash(origin);
1274 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1275 const Index voxelIdx = voxelIter.pos();
1276 const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1277 const int64_t key = leafKey + hash(ijk);
1278 if (key == voxelSliceKey) {
1279 nodeMask->setOn(voxelIdx);
1285 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1292 inline static Coord ijk(
const Coord &p,
int i) {
return p + FastSweeping::mOffset[i]; }
1297 inline operator bool()
const {
return v < SdfValueT(1000); }
1302 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() :
nullptr;
1304 const SdfValueT h =
static_cast<SdfValueT
>(mParent->mSdfGrid->voxelSize()[0]);
1305 const SdfValueT sqrt2h =
math::Sqrt(SdfValueT(2))*h;
1307 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1309 int64_t voxelSliceIndex(0);
1311 auto kernel = [&](
const tbb::blocked_range<size_t>& range) {
1312 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1314 SdfAccT acc1(mParent->mSdfGrid->tree());
1315 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ?
new ExtAccT(*tree2) :
nullptr);
1316 SdfValueT absV, sign, update, D;
1319 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1323 for (
size_t i = range.begin(); i < range.end(); ++i) {
1328 const LeafSlice& leafSlice = leafSliceArray[i];
1329 const size_t leafIdx = leafSlice.first;
1330 const NodeMaskPtrT& nodeMask = leafSlice.second;
1332 const Coord& origin = leafNodeOrigins[leafIdx];
1335 for (
auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1338 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1345 if (!(d1 || d2 || d3))
continue;
1350 SdfValueT &value =
const_cast<SdfValueT&
>(acc1.
getValue(ijk));
1353 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1359 if (d2 < d1) std::swap(d1, d2);
1360 if (d3 < d2) std::swap(d2, d3);
1361 if (d2 < d1) std::swap(d1, d2);
1367 if (update <= d2.
v) {
1368 if (update < absV) {
1369 value = sign * update;
1370 if (acc2) acc2->setValue(ijk, acc2->getValue(d1(ijk)));
1379 if (d2.
v <= sqrt2h + d1.
v) {
1381 update = SdfValueT(0.5) * (d1.
v + d2.
v + std::sqrt(D));
1382 if (update > d2.
v && update <= d3.
v) {
1383 if (update < absV) {
1384 value = sign * update;
1389 const SdfValueT w = SdfValueT(1)/(d1.
v+d2.
v);
1390 acc2->setValue(ijk, w*(d1.
v*acc2->getValue(d1(ijk)) +
1391 d2.
v*acc2->getValue(d2(ijk))));
1402 const SdfValueT d123 = d1.
v + d2.
v + d3.
v;
1403 D = d123*d123 - SdfValueT(3)*(d1.
v*d1.
v + d2.
v*d2.
v + d3.
v*d3.
v - h * h);
1404 if (D >= SdfValueT(0)) {
1405 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));
1407 if (update < absV) {
1408 value = sign * update;
1414 const SdfValueT w = SdfValueT(1)/(d1.
v+d2.
v+d3.
v);
1415 acc2->setValue(ijk, w*(d1.
v*acc2->getValue(d1(ijk)) +
1416 d2.
v*acc2->getValue(d2(ijk)) +
1417 d3.
v*acc2->getValue(d3(ijk))));
1425 #ifdef BENCHMARK_FAST_SWEEPING
1429 for (
size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1430 voxelSliceIndex = mVoxelSliceKeys[i];
1431 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1434 #ifdef BENCHMARK_FAST_SWEEPING
1435 timer.
restart(
"Backward sweeps");
1437 for (
size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1438 voxelSliceIndex = mVoxelSliceKeys[i-1];
1439 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1442 #ifdef BENCHMARK_FAST_SWEEPING
1448 using NodeMaskT =
typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1449 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1452 using LeafSlice = std::pair<size_t, NodeMaskPtrT>;
1453 using LeafSliceArray = std::deque<LeafSlice>;
1454 using VoxelSliceMap = std::map<int64_t, LeafSliceArray>;
1458 VoxelSliceMap mVoxelSliceMap;
1459 std::vector<int64_t> mVoxelSliceKeys;
1464 template<
typename Gr
idT>
1467 typename GridT::ValueType isoValue,
1471 if (fs.initSdf(fogGrid, isoValue,
false)) fs.sweep(nIter);
1472 return fs.sdfGrid();
1475 template<
typename Gr
idT>
1478 typename GridT::ValueType isoValue,
1482 if (fs.initSdf(sdfGrid, isoValue,
true)) fs.sweep(nIter);
1483 return fs.sdfGrid();
1486 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1487 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1490 const ExtValueT& background,
1491 typename FogGridT::ValueType isoValue,
1495 if (fs.initExt(fogGrid, op, background, isoValue,
false)) fs.sweep(nIter);
1496 return fs.extGrid();
1499 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1500 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1503 const ExtValueT &background,
1504 typename SdfGridT::ValueType isoValue,
1508 if (fs.initExt(sdfGrid, op, background, isoValue,
true)) fs.sweep(nIter);
1509 return fs.extGrid();
1512 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1513 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1516 const ExtValueT &background,
1517 typename FogGridT::ValueType isoValue,
1521 if (fs.initExt(fogGrid, op, background, isoValue,
false)) fs.sweep(nIter);
1522 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1525 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1526 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1529 const ExtValueT &background,
1530 typename SdfGridT::ValueType isoValue,
1534 if (fs.initExt(sdfGrid, op, background, isoValue,
true)) fs.sweep(nIter);
1535 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1538 template<
typename Gr
idT>
1546 if (fs.initDilate(sdfGrid, dilation, nn)) fs.sweep(nIter);
1547 return fs.sdfGrid();
1550 template<
typename Gr
idT,
typename MaskTreeT>
1554 bool ignoreActiveTiles,
1558 if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1559 return fs.sdfGrid();
1566 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED