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>
333 const Grid<MaskTreeT> &mask,
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 std::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;
944 template <
typename SdfGr
idT,
typename ExtValueT>
945 template <
typename OpT>
949 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
951 mOpPool(
nullptr), mSdfGrid(parent.mSdfGrid.get()),
952 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
953 InitExt(
const InitExt&) =
default;
954 InitExt&
operator=(
const InitExt&) =
delete;
955 void run(SdfValueT isoValue,
const OpT &opPrototype,
bool isInputSdf)
957 static_assert(std::is_convertible<decltype(opPrototype(
Vec3d(0))),ExtValueT>::value,
"Invalid return type of functor");
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();
968 if (isInputSdf && hasActiveTiles) {
972 #ifdef BENCHMARK_FAST_SWEEPING
973 util::CpuTimer timer(
"Initialize voxels");
976 mParent->mSweepMask.clear();
977 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
981 OpPoolT opPool(opPrototype);
984 tree::LeafManager<SdfTreeT> mgr(tree1, 1);
985 tbb::parallel_for(mgr.leafRange(32), *
this);
986 mgr.swapLeafBuffer(1);
989 #ifdef BENCHMARK_FAST_SWEEPING
990 timer.restart(
"Initialize tiles");
993 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
994 mgr.foreachBottomUp(*
this);
996 if (hasActiveTiles) {
997 #ifdef BENCHMARK_FAST_SWEEPING
998 timer.restart(
"Voxelizing active tiles");
1000 tree1.voxelizeActiveTiles();
1001 tree2.voxelizeActiveTiles();
1007 mParent->computeSweepMaskLeafOrigins();
1009 #ifdef BENCHMARK_FAST_SWEEPING
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& )
const {
if (update) sum2 = ext; }
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 ,
const SdfT& d2)
const { sum2 +=
static_cast<ExtValueT
>(d2 * ext); }
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& )
const {
return extSum; }
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); }
1028 void operator()(
const LeafRange& r)
const
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();
1036 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1037 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1038 SdfValueT *sdf = leafIter.buffer(1).data();
1039 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();
1040 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1041 const SdfValueT value = *voxelIter;
1042 const bool isAbove = value > isoValue;
1043 if (!voxelIter.isValueOn()) {
1044 sdf[voxelIter.pos()] = isAbove ? above : -above;
1046 const Coord ijk = voxelIter.getCoord();
1047 stencil.moveTo(ijk, value);
1048 const auto mask = stencil.intersectionMask( isoValue );
1050 sdf[voxelIter.pos()] = isAbove ? above : -above;
1054 sweepMaskAcc.setValueOff(ijk);
1055 const SdfValueT delta = value - isoValue;
1057 sdf[voxelIter.pos()] = 0;
1058 ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1061 ExtValueT sum2 = zeroVal<ExtValueT>();
1067 for (
int n=0, i=0; i<6;) {
1069 if (mask.test(i++)) {
1070 d =
math::Abs(delta/(value-stencil.getValue(i)));
1073 if (mask.test(i++)) {
1074 d2 =
math::Abs(delta/(value-stencil.getValue(i)));
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]));
1087 sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1088 if (d < minD) minD = d;
1091 ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1092 sdf[voxelIter.pos()] = isAbove ? h /
math::Sqrt(sum1) : -h / math::
Sqrt(sum1);
1100 template<
typename RootOrInternalNodeT>
1101 void operator()(
const RootOrInternalNodeT& node)
const
1104 for (
auto it = node.cbeginValueAll(); it; ++it) {
1105 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1106 v = v > isoValue ? above : -above;
1114 SdfValueT mIsoValue;
1115 SdfValueT mAboveSign;
1119 template <
typename SdfGr
idT,
typename ExtValueT>
1120 template <
typename MaskTreeT>
1121 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1123 using LeafRange =
typename tree::LeafManager<const MaskTreeT>::LeafRange;
1125 mSdfGrid(parent.mSdfGrid.get()) {}
1126 MaskKernel(
const MaskKernel &parent) =
default;
1127 MaskKernel&
operator=(
const MaskKernel&) =
delete;
1129 void run(
const MaskTreeT &mask)
1131 #ifdef BENCHMARK_FAST_SWEEPING
1132 util::CpuTimer timer;
1134 auto &lsTree = mSdfGrid->tree();
1138 #ifdef BENCHMARK_FAST_SWEEPING
1139 timer.restart(
"Changing background value");
1143 #ifdef BENCHMARK_FAST_SWEEPING
1144 timer.restart(
"Union with mask");
1146 lsTree.topologyUnion(mask);
1149 tree::LeafManager<const MaskTreeT> mgr(mask);
1151 #ifdef BENCHMARK_FAST_SWEEPING
1152 timer.restart(
"Initializing grid and sweep mask");
1155 mParent->mSweepMask.clear();
1156 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1158 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1159 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1160 LeafManagerT leafManager(mParent->mSweepMask);
1162 auto kernel = [&](LeafT& leaf,
size_t ) {
1164 SdfAccT acc(mSdfGrid->tree());
1167 SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1168 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1169 if (
math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1171 voxelIter.setValue(
false);
1175 leafManager.foreach( kernel );
1178 mParent->computeSweepMaskLeafOrigins();
1180 #ifdef BENCHMARK_FAST_SWEEPING
1191 template <
typename SdfGr
idT,
typename ExtValueT>
1199 template<
typename HashOp>
1202 #ifdef BENCHMARK_FAST_SWEEPING
1207 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1210 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1211 LeafManagerT leafManager(maskTree);
1218 constexpr
int maskOffset = LeafT::DIM * 3;
1219 constexpr
int maskRange = maskOffset * 2;
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);
1230 leafManager.foreach( kernel1 );
1235 using ThreadLocalMap = std::unordered_map<int64_t, 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);
1249 leafManager.foreach( kernel2 );
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());
1264 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1265 for (
const auto& it : mVoxelSliceMap) {
1266 mVoxelSliceKeys.push_back(it.first);
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>();
1278 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
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);
1306 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1313 inline static Coord ijk(
const Coord &p,
int i) {
return p + FastSweeping::mOffset[i]; }
1318 inline operator bool()
const {
return v < SdfValueT(1000); }
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& ,
const ExtT& v1,
const ExtT& v2)
const {
return d1.
v < d2.
v ? v1 : v2; }
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)); }
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& ,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
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));
1345 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() :
nullptr;
1347 const SdfValueT h =
static_cast<SdfValueT
>(mParent->mSdfGrid->voxelSize()[0]);
1348 const SdfValueT sqrt2h =
math::Sqrt(SdfValueT(2))*h;
1350 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1352 int64_t voxelSliceIndex(0);
1354 auto kernel = [&](
const tbb::blocked_range<size_t>& range) {
1355 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1357 SdfAccT acc1(mParent->mSdfGrid->tree());
1358 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ?
new ExtAccT(*tree2) :
nullptr);
1359 SdfValueT absV, sign, update, D;
1362 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1366 for (
size_t i = range.begin(); i < range.end(); ++i) {
1371 const LeafSlice& leafSlice = leafSliceArray[i];
1372 const size_t leafIdx = leafSlice.first;
1373 const NodeMaskPtrT& nodeMask = leafSlice.second;
1375 const Coord& origin = leafNodeOrigins[leafIdx];
1378 for (
auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1381 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1388 if (!(d1 || d2 || d3))
continue;
1393 SdfValueT &value =
const_cast<SdfValueT&
>(acc1.
getValue(ijk));
1396 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1402 if (d2 < d1) std::swap(d1, d2);
1403 if (d3 < d2) std::swap(d2, d3);
1404 if (d2 < d1) std::swap(d1, d2);
1410 if (update <= d2.
v) {
1411 if (update < absV) {
1412 value = sign * update;
1413 if (acc2) acc2->setValue(ijk, acc2->getValue(d1(ijk)));
1422 if (d2.
v <= sqrt2h + d1.
v) {
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;
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);
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)) {
1450 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));
1452 if (update < absV) {
1453 value = sign * update;
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);
1472 #ifdef BENCHMARK_FAST_SWEEPING
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);
1481 #ifdef BENCHMARK_FAST_SWEEPING
1482 timer.
restart(
"Backward sweeps");
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);
1489 #ifdef BENCHMARK_FAST_SWEEPING
1495 using NodeMaskT =
typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1496 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1499 using LeafSlice = std::pair<size_t, NodeMaskPtrT>;
1500 using LeafSliceArray = std::deque<LeafSlice>;
1501 using VoxelSliceMap = std::map<int64_t, LeafSliceArray>;
1505 VoxelSliceMap mVoxelSliceMap;
1506 std::vector<int64_t> mVoxelSliceKeys;
1511 template<
typename Gr
idT>
1514 typename GridT::ValueType isoValue,
1518 if (fs.initSdf(fogGrid, isoValue,
false)) fs.sweep(nIter);
1519 return fs.sdfGrid();
1522 template<
typename Gr
idT>
1525 typename GridT::ValueType isoValue,
1529 if (fs.initSdf(sdfGrid, isoValue,
true)) fs.sweep(nIter);
1530 return fs.sdfGrid();
1533 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1534 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1537 const ExtValueT& background,
1538 typename FogGridT::ValueType isoValue,
1542 if (fs.initExt(fogGrid, op, background, isoValue,
false)) fs.sweep(nIter);
1543 return fs.extGrid();
1546 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1547 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1550 const ExtValueT &background,
1551 typename SdfGridT::ValueType isoValue,
1555 if (fs.initExt(sdfGrid, op, background, isoValue,
true)) fs.sweep(nIter);
1556 return fs.extGrid();
1559 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1560 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1563 const ExtValueT &background,
1564 typename FogGridT::ValueType isoValue,
1568 if (fs.initExt(fogGrid, op, background, isoValue,
false)) fs.sweep(nIter);
1569 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1572 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1573 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1576 const ExtValueT &background,
1577 typename SdfGridT::ValueType isoValue,
1581 if (fs.initExt(sdfGrid, op, background, isoValue,
true)) fs.sweep(nIter);
1582 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1585 template<
typename Gr
idT>
1593 if (fs.initDilate(sdfGrid, dilation, nn)) fs.sweep(nIter);
1594 return fs.sdfGrid();
1597 template<
typename Gr
idT,
typename MaskTreeT>
1601 bool ignoreActiveTiles,
1605 if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1606 return fs.sdfGrid();
1613 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED