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
99 template<
typename Gr
idT>
102 typename GridT::ValueType isoValue,
132 template<
typename Gr
idT>
135 typename GridT::ValueType isoValue = 0,
188 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
189 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
192 const ExtValueT& background,
193 typename FogGridT::ValueType isoValue,
196 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
246 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
247 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
250 const ExtValueT &background,
251 typename SdfGridT::ValueType isoValue = 0,
254 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
309 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
310 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
313 const ExtValueT &background,
314 typename FogGridT::ValueType isoValue,
317 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
372 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
373 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
376 const ExtValueT &background,
377 typename SdfGridT::ValueType isoValue = 0,
380 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
405 template<
typename Gr
idT>
431 template<
typename Gr
idT,
typename MaskTreeT>
435 bool ignoreActiveTiles =
false,
450 template<
typename SdfGr
idT,
typename ExtValueT =
typename SdfGr
idT::ValueType>
453 static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
454 "FastSweeping requires SdfGridT to have floating-point values");
456 using SdfValueT =
typename SdfGridT::ValueType;
457 using SdfTreeT =
typename SdfGridT::TreeType;
462 using ExtGridT =
typename SdfGridT::template ValueConverter<ExtValueT>::Type;
463 using ExtTreeT =
typename ExtGridT::TreeType;
467 using SweepMaskTreeT =
typename SdfTreeT::template ValueConverter<ValueMask>::Type;
490 typename SdfGridT::Ptr
sdfGrid() {
return mSdfGrid; }
498 typename ExtGridT::Ptr
extGrid() {
return mExtGrid; }
528 bool initSdf(
const SdfGridT &sdfGrid, SdfValueT isoValue,
bool isInputSdf);
576 template <
typename ExtOpT>
579 const ExtValueT &background,
583 const typename ExtGridT::ConstPtr extGrid =
nullptr);
606 bool initDilate(
const SdfGridT &sdfGrid,
629 template<
typename MaskTreeT>
630 bool initMask(
const SdfGridT &sdfGrid,
const Grid<MaskTreeT> &mask,
bool ignoreActiveTiles =
false);
644 void sweep(
int nIter = 1,
645 bool finalize =
true);
657 bool isValid()
const {
return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
673 void computeSweepMaskLeafOrigins();
683 struct SweepingKernel;
686 static const Coord mOffset[6];
689 typename SdfGridT::Ptr mSdfGrid;
690 typename ExtGridT::Ptr mExtGrid;
691 typename ExtGridT::Ptr mExtGridInput;
692 SweepMaskTreeT mSweepMask;
693 std::vector<Coord> mSweepMaskLeafOrigins;
694 size_t mSweepingVoxelCount, mBoundaryVoxelCount;
702 template <
typename SdfGr
idT,
typename ExtValueT>
703 const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
707 template <
typename SdfGr
idT,
typename ExtValueT>
709 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(
FastSweepingDomain::
SWEEP_ALL), mIsInputSdf(true)
713 template <
typename SdfGr
idT,
typename ExtValueT>
719 if (mExtGridInput) mExtGridInput.reset();
720 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
725 template <
typename SdfGr
idT,
typename ExtValueT>
731 mSweepMask.voxelizeActiveTiles();
734 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
735 LeafManagerT leafManager(mSweepMask);
737 mSweepMaskLeafOrigins.resize(leafManager.leafCount());
738 std::atomic<size_t> sweepingVoxelCount{0};
739 auto kernel = [&](
const LeafT& leaf,
size_t leafIdx) {
740 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
741 sweepingVoxelCount += leaf.onVoxelCount();
743 leafManager.foreach(kernel,
true, 1024);
745 mBoundaryVoxelCount = 0;
746 mSweepingVoxelCount = sweepingVoxelCount;
748 const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
749 assert( totalCount >= mSweepingVoxelCount );
750 mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
754 template <
typename SdfGr
idT,
typename ExtValueT>
758 mSdfGrid = fogGrid.deepCopy();
759 mIsInputSdf = isInputSdf;
761 kernel.
run(isoValue);
762 return this->isValid();
765 template <
typename SdfGr
idT,
typename ExtValueT>
766 template <
typename OpT>
772 if (extGrid->transform() != fogGrid.transform())
773 OPENVDB_THROW(
RuntimeError,
"FastSweeping::initExt extension grid input should have the same transform as Fog/SDF grid!");
777 mSdfGrid = fogGrid.deepCopy();
778 mExtGrid = createGrid<ExtGridT>( background );
779 mSweepDirection = mode;
780 mIsInputSdf = isInputSdf;
782 mExtGridInput = extGrid->deepCopy();
784 mExtGrid->topologyUnion( *mSdfGrid );
785 InitExt<OpT> kernel(*
this);
786 kernel.run(isoValue, op);
787 return this->isValid();
791 template <
typename SdfGr
idT,
typename ExtValueT>
795 mSdfGrid = sdfGrid.deepCopy();
796 mSweepDirection = mode;
798 kernel.
run(dilate, nn);
799 return this->isValid();
802 template <
typename SdfGr
idT,
typename ExtValueT>
803 template<
typename MaskTreeT>
807 mSdfGrid = sdfGrid.deepCopy();
809 if (mSdfGrid->transform() != mask.
transform()) {
814 using T =
typename MaskTreeT::template ValueConverter<bool>::Type;
816 tmp->
tree().voxelizeActiveTiles();
817 MaskKernel<T> kernel(*
this);
818 kernel.run(tmp->
tree());
820 if (ignoreActiveTiles || !mask.
tree().hasActiveTiles()) {
821 MaskKernel<MaskTreeT> kernel(*
this);
822 kernel.run(mask.
tree());
824 using T =
typename MaskTreeT::template ValueConverter<ValueMask>::Type;
826 tmp.voxelizeActiveTiles(
true);
827 MaskKernel<T> kernel(*
this);
831 return this->isValid();
834 template <
typename SdfGr
idT,
typename ExtValueT>
842 " a non-null reference extension grid input.");
844 if (this->boundaryVoxelCount() == 0) {
846 }
else if (this->sweepingVoxelCount() == 0) {
851 std::deque<SweepingKernel> kernels;
852 for (
int i = 0; i < 4; i++) kernels.emplace_back(*
this);
855 #ifdef BENCHMARK_FAST_SWEEPING
860 tbb::task_group tasks;
861 tasks.run([&] { kernels[0].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]+a[2]; }); });
862 tasks.run([&] { kernels[1].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]-a[2]; }); });
863 tasks.run([&] { kernels[2].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]+a[2]; }); });
864 tasks.run([&] { kernels[3].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]-a[2]; }); });
867 #ifdef BENCHMARK_FAST_SWEEPING
873 for (
int i = 0; i < nIter; ++i) {
878 #ifdef BENCHMARK_FAST_SWEEPING
882 auto e = kernel.
run(*mSdfGrid);
884 #ifdef BENCHMARK_FAST_SWEEPING
885 std::cerr <<
"Min = " << e.min() <<
" Max = " << e.max() << std::endl;
886 timer.
restart(
"Changing asymmetric background value");
890 #ifdef BENCHMARK_FAST_SWEEPING
899 template <
typename SdfGr
idT,
typename ExtValueT>
910 tbb::parallel_reduce(mgr.
leafRange(), *
this);
916 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
917 for (
auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
918 const SdfValueT v = *voxelIter;
919 if (v < mMin) mMin = v;
920 if (v > mMax) mMax = v;
927 if (other.
mMin < mMin) mMin = other.
mMin;
928 if (other.
mMax > mMax) mMax = other.
mMax;
937 template <
typename SdfGr
idT,
typename ExtValueT>
943 mBackground(parent.mSdfGrid->background())
945 mSdfGridInput = mParent->mSdfGrid->deepCopy();
952 #ifdef BENCHMARK_FAST_SWEEPING
957 #ifdef BENCHMARK_FAST_SWEEPING
958 timer.
restart(
"Changing background value");
963 #ifdef BENCHMARK_FAST_SWEEPING
964 timer.
restart(
"Dilating and updating mgr (parallel)");
974 #ifdef BENCHMARK_FAST_SWEEPING
975 timer.
restart(
"Initializing grid and sweep mask");
978 mParent->mSweepMask.clear();
979 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
982 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
986 LeafManagerT leafManager(mParent->mSdfGrid->tree());
988 auto kernel = [&](LeafT& leaf,
size_t ) {
990 const SdfValueT background = mBackground;
991 auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
992 SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
994 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
995 const SdfValueT value = *voxelIter;
996 SdfValueT inputValue;
997 const Coord ijk = voxelIter.getCoord();
1000 maskLeaf->setValueOff(voxelIter.pos());
1004 voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1007 if (value > 0) voxelIter.setValue(Unknown);
1009 maskLeaf->setValueOff(voxelIter.pos());
1010 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1011 if ( !isInputOn ) voxelIter.setValueOff();
1012 else voxelIter.setValue(inputValue);
1016 if (value < 0) voxelIter.setValue(-Unknown);
1018 maskLeaf->setValueOff(voxelIter.pos());
1019 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1020 if ( !isInputOn ) voxelIter.setValueOff();
1021 else voxelIter.setValue(inputValue);
1029 leafManager.foreach( kernel );
1032 mParent->computeSweepMaskLeafOrigins();
1034 #ifdef BENCHMARK_FAST_SWEEPING
1047 template <
typename SdfGr
idT,
typename ExtValueT>
1052 mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
1058 mIsoValue = isoValue;
1059 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1060 SdfTreeT &tree = mSdfGrid->tree();
1061 const bool hasActiveTiles = tree.hasActiveTiles();
1063 if (mParent->mIsInputSdf && hasActiveTiles) {
1067 #ifdef BENCHMARK_FAST_SWEEPING
1070 mParent->mSweepMask.clear();
1071 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1075 tbb::parallel_for(mgr.
leafRange(32), *
this);
1079 #ifdef BENCHMARK_FAST_SWEEPING
1080 timer.
restart(
"Initialize tiles - new");
1084 mgr.foreachBottomUp(*
this);
1086 if (hasActiveTiles) tree.voxelizeActiveTiles();
1090 mParent->computeSweepMaskLeafOrigins();
1098 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1099 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1100 SdfValueT* sdf = leafIter.buffer(1).data();
1101 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1102 const SdfValueT value = *voxelIter;
1103 const bool isAbove = value > isoValue;
1104 if (!voxelIter.isValueOn()) {
1105 sdf[voxelIter.pos()] = isAbove ? above : -above;
1107 const Coord ijk = voxelIter.getCoord();
1108 stencil.
moveTo(ijk, value);
1111 sdf[voxelIter.pos()] = isAbove ? above : -above;
1115 const SdfValueT delta = value - isoValue;
1117 sdf[voxelIter.pos()] = 0;
1120 for (
int i=0; i<6;) {
1123 if (mask.test(i++)) {
1137 template<
typename RootOrInternalNodeT>
1141 for (
auto it = node.cbeginValueAll(); it; ++it) {
1142 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1143 v = v > isoValue ? above : -above;
1156 template <
typename SdfGr
idT,
typename ExtValueT>
1157 template <
typename OpT>
1161 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1163 mOpPool(
nullptr), mSdfGrid(parent.mSdfGrid.get()),
1164 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1165 InitExt(
const InitExt&) =
default;
1166 InitExt&
operator=(
const InitExt&) =
delete;
1167 void run(SdfValueT isoValue,
const OpT &opPrototype)
1169 static_assert(std::is_convertible<decltype(opPrototype(
Vec3d(0))),ExtValueT>::value,
"Invalid return type of functor");
1174 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1175 mIsoValue = isoValue;
1176 auto &tree1 = mSdfGrid->tree();
1177 auto &tree2 = mExtGrid->tree();
1178 const bool hasActiveTiles = tree1.hasActiveTiles();
1180 if (mParent->mIsInputSdf && hasActiveTiles) {
1184 #ifdef BENCHMARK_FAST_SWEEPING
1188 mParent->mSweepMask.clear();
1189 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1193 OpPoolT opPool(opPrototype);
1197 tbb::parallel_for(mgr.leafRange(32), *
this);
1198 mgr.swapLeafBuffer(1);
1201 #ifdef BENCHMARK_FAST_SWEEPING
1202 timer.restart(
"Initialize tiles");
1205 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
1206 mgr.foreachBottomUp(*
this);
1208 if (hasActiveTiles) {
1209 #ifdef BENCHMARK_FAST_SWEEPING
1210 timer.restart(
"Voxelizing active tiles");
1212 tree1.voxelizeActiveTiles();
1213 tree2.voxelizeActiveTiles();
1219 mParent->computeSweepMaskLeafOrigins();
1221 #ifdef BENCHMARK_FAST_SWEEPING
1227 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1228 void sumHelper(ExtT& sum2, ExtT ext,
bool update,
const SdfT& )
const {
if (update) sum2 = ext; }
1231 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1232 void sumHelper(ExtT& sum2, ExtT ext,
bool ,
const SdfT& d2)
const { sum2 +=
static_cast<ExtValueT
>(d2 * ext); }
1234 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1235 ExtT extValHelper(ExtT extSum,
const SdfT& )
const {
return extSum; }
1237 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1238 ExtT extValHelper(ExtT extSum,
const SdfT& sdfSum)
const {
return ExtT((SdfT(1) / sdfSum) * extSum); }
1240 void operator()(
const LeafRange& r)
const
1242 ExtAccT acc(mExtGrid->tree());
1243 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1244 math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1245 const math::Transform& xform = mExtGrid->transform();
1246 typename OpPoolT::reference op = mOpPool->local();
1248 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1249 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1250 SdfValueT *sdf = leafIter.buffer(1).data();
1251 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();
1252 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1253 const SdfValueT value = *voxelIter;
1254 const bool isAbove = value > isoValue;
1255 if (!voxelIter.isValueOn()) {
1256 sdf[voxelIter.pos()] = isAbove ? above : -above;
1258 const Coord ijk = voxelIter.getCoord();
1259 stencil.moveTo(ijk, value);
1260 const auto mask = stencil.intersectionMask( isoValue );
1262 sdf[voxelIter.pos()] = isAbove ? above : -above;
1266 sweepMaskAcc.setValueOff(ijk);
1267 const SdfValueT delta = value - isoValue;
1269 sdf[voxelIter.pos()] = 0;
1270 ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1273 ExtValueT sum2 = zeroVal<ExtValueT>();
1279 for (
int n=0, i=0; i<6;) {
1281 if (mask.test(i++)) {
1282 d =
math::Abs(delta/(value-stencil.getValue(i)));
1285 if (mask.test(i++)) {
1286 d2 =
math::Abs(delta/(value-stencil.getValue(i)));
1295 const Vec3R xyz(
static_cast<SdfValueT
>(ijk[0])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][0]),
1296 static_cast<SdfValueT
>(ijk[1])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][1]),
1297 static_cast<SdfValueT
>(ijk[2])+d*
static_cast<SdfValueT
>(FastSweeping::mOffset[n][2]));
1299 sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1300 if (d < minD) minD = d;
1303 ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1304 sdf[voxelIter.pos()] = isAbove ? h /
math::Sqrt(sum1) : -h / math::
Sqrt(sum1);
1312 template<
typename RootOrInternalNodeT>
1313 void operator()(
const RootOrInternalNodeT& node)
const
1316 for (
auto it = node.cbeginValueAll(); it; ++it) {
1317 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1318 v = v > isoValue ? above : -above;
1326 SdfValueT mIsoValue;
1327 SdfValueT mAboveSign;
1331 template <
typename SdfGr
idT,
typename ExtValueT>
1332 template <
typename MaskTreeT>
1333 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1335 using LeafRange =
typename tree::LeafManager<const MaskTreeT>::LeafRange;
1337 mSdfGrid(parent.mSdfGrid.get()) {}
1338 MaskKernel(
const MaskKernel &parent) =
default;
1339 MaskKernel&
operator=(
const MaskKernel&) =
delete;
1341 void run(
const MaskTreeT &mask)
1343 #ifdef BENCHMARK_FAST_SWEEPING
1344 util::CpuTimer timer;
1346 auto &lsTree = mSdfGrid->tree();
1350 #ifdef BENCHMARK_FAST_SWEEPING
1351 timer.restart(
"Changing background value");
1355 #ifdef BENCHMARK_FAST_SWEEPING
1356 timer.restart(
"Union with mask");
1358 lsTree.topologyUnion(mask);
1361 tree::LeafManager<const MaskTreeT> mgr(mask);
1363 #ifdef BENCHMARK_FAST_SWEEPING
1364 timer.restart(
"Initializing grid and sweep mask");
1367 mParent->mSweepMask.clear();
1368 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1370 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1371 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1372 LeafManagerT leafManager(mParent->mSweepMask);
1374 auto kernel = [&](LeafT& leaf,
size_t ) {
1376 SdfAccT acc(mSdfGrid->tree());
1379 SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1380 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1381 if (
math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1383 voxelIter.setValue(
false);
1387 leafManager.foreach( kernel );
1390 mParent->computeSweepMaskLeafOrigins();
1392 #ifdef BENCHMARK_FAST_SWEEPING
1403 template <
typename SdfGr
idT,
typename ExtValueT>
1411 template<
typename HashOp>
1414 #ifdef BENCHMARK_FAST_SWEEPING
1419 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1422 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1423 LeafManagerT leafManager(maskTree);
1430 constexpr
int maskOffset = LeafT::DIM * 3;
1431 constexpr
int maskRange = maskOffset * 2;
1434 std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1435 auto kernel1 = [&](
const LeafT& leaf,
size_t leafIdx) {
1436 const size_t leafOffset = leafIdx * maskRange;
1437 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1438 const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1439 leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1442 leafManager.foreach( kernel1 );
1447 using ThreadLocalMap = std::unordered_map<int64_t, std::deque<size_t>>;
1448 tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1449 auto kernel2 = [&](
const LeafT& leaf,
size_t leafIdx) {
1450 ThreadLocalMap& map = pool.local();
1451 const Coord& origin = leaf.origin();
1452 const int64_t leafKey = hash(origin);
1453 const size_t leafOffset = leafIdx * maskRange;
1454 for (
int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1455 if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1456 const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1457 map[voxelSliceKey].emplace_back(leafIdx);
1461 leafManager.foreach( kernel2 );
1466 for (
auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1467 const ThreadLocalMap& map = *poolIt;
1468 for (
const auto& it : map) {
1469 for (
const size_t leafIdx : it.second) {
1470 mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1476 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1477 for (
const auto& it : mVoxelSliceMap) {
1478 mVoxelSliceKeys.push_back(it.first);
1482 auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1483 for (
size_t i = range.begin(); i < range.end(); i++) {
1484 const int64_t key = mVoxelSliceKeys[i];
1485 for (
auto& it : mVoxelSliceMap[key]) {
1486 it.second = std::make_unique<NodeMaskT>();
1490 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1497 auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1498 for (
size_t i = range.begin(); i < range.end(); i++) {
1499 const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1500 LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1501 for (LeafSlice& leafSlice : leafSliceArray) {
1502 const size_t leafIdx = leafSlice.first;
1503 NodeMaskPtrT& nodeMask = leafSlice.second;
1504 const LeafT& leaf = leafManager.leaf(leafIdx);
1505 const Coord& origin = leaf.origin();
1506 const int64_t leafKey = hash(origin);
1507 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1508 const Index voxelIdx = voxelIter.pos();
1509 const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1510 const int64_t key = leafKey + hash(ijk);
1511 if (key == voxelSliceKey) {
1512 nodeMask->setOn(voxelIdx);
1518 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1525 inline static Coord ijk(
const Coord &p,
int i) {
return p + FastSweeping::mOffset[i]; }
1530 inline operator bool()
const {
return v < SdfValueT(1000); }
1534 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1535 ExtT
twoNghbr(
const NN& d1,
const NN& d2,
const SdfT& ,
const ExtT& v1,
const ExtT& v2)
const {
return d1.
v < d2.
v ? v1 : v2; }
1538 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1539 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)); }
1542 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1543 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& ,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1550 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1551 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& w,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1552 return ExtT(w*(d1.
v*v1 + d2.
v*v2 + d3.
v*v3));
1557 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() :
nullptr;
1558 typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() :
nullptr;
1560 const SdfValueT h =
static_cast<SdfValueT
>(mParent->mSdfGrid->voxelSize()[0]);
1561 const SdfValueT sqrt2h =
math::Sqrt(SdfValueT(2))*h;
1563 const bool isInputSdf = mParent->mIsInputSdf;
1570 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1572 int64_t voxelSliceIndex(0);
1574 auto kernel = [&](
const tbb::blocked_range<size_t>& range) {
1575 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1577 SdfAccT acc1(mParent->mSdfGrid->tree());
1578 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ?
new ExtAccT(*tree2) :
nullptr);
1579 auto acc3 = std::unique_ptr<ExtAccT>(tree3 ?
new ExtAccT(*tree3) :
nullptr);
1580 SdfValueT absV, sign, update, D;
1583 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1587 for (
size_t i = range.begin(); i < range.end(); ++i) {
1592 const LeafSlice& leafSlice = leafSliceArray[i];
1593 const size_t leafIdx = leafSlice.first;
1594 const NodeMaskPtrT& nodeMask = leafSlice.second;
1596 const Coord& origin = leafNodeOrigins[leafIdx];
1599 for (
auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1602 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1609 if (!(d1 || d2 || d3))
continue;
1614 SdfValueT &value =
const_cast<SdfValueT&
>(acc1.
getValue(ijk));
1617 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1623 if (d2 < d1) std::swap(d1, d2);
1624 if (d3 < d2) std::swap(d2, d3);
1625 if (d2 < d1) std::swap(d1, d2);
1631 if (update <= d2.
v) {
1632 if (update < absV) {
1633 value = sign * update;
1636 ExtValueT updateExt = acc2->getValue(d1(ijk));
1638 if (
isInputSdf) updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1639 else updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1642 if (
isInputSdf) updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1643 else updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1645 acc2->setValue(ijk, updateExt);
1655 if (d2.
v <= sqrt2h + d1.
v) {
1657 update = SdfValueT(0.5) * (d1.
v + d2.
v + std::sqrt(D));
1658 if (update > d2.
v && update <= d3.
v) {
1659 if (update < absV) {
1660 value = sign * update;
1665 const SdfValueT w = SdfValueT(1)/(d1.
v+d2.
v);
1666 const ExtValueT v1 = acc2->getValue(d1(ijk));
1667 const ExtValueT v2 = acc2->getValue(d2(ijk));
1668 const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1670 ExtValueT updateExt = extVal;
1672 if (
isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1673 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1676 if (
isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1677 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1679 acc2->setValue(ijk, updateExt);
1690 const SdfValueT d123 = d1.
v + d2.
v + d3.
v;
1691 D = d123*d123 - SdfValueT(3)*(d1.
v*d1.
v + d2.
v*d2.
v + d3.
v*d3.
v - h * h);
1692 if (D >= SdfValueT(0)) {
1693 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));
1695 if (update < absV) {
1696 value = sign * update;
1702 const SdfValueT w = SdfValueT(1)/(d1.
v+d2.
v+d3.
v);
1703 const ExtValueT v1 = acc2->getValue(d1(ijk));
1704 const ExtValueT v2 = acc2->getValue(d2(ijk));
1705 const ExtValueT v3 = acc2->getValue(d3(ijk));
1706 const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1708 ExtValueT updateExt = extVal;
1710 if (
isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1711 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1714 if (
isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1715 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1717 acc2->setValue(ijk, updateExt);
1725 #ifdef BENCHMARK_FAST_SWEEPING
1729 for (
size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1730 voxelSliceIndex = mVoxelSliceKeys[i];
1731 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1734 #ifdef BENCHMARK_FAST_SWEEPING
1735 timer.
restart(
"Backward sweeps");
1737 for (
size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1738 voxelSliceIndex = mVoxelSliceKeys[i-1];
1739 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1742 #ifdef BENCHMARK_FAST_SWEEPING
1748 using NodeMaskT =
typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1749 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1752 using LeafSlice = std::pair<size_t, NodeMaskPtrT>;
1753 using LeafSliceArray = std::deque<LeafSlice>;
1754 using VoxelSliceMap = std::map<int64_t, LeafSliceArray>;
1758 VoxelSliceMap mVoxelSliceMap;
1759 std::vector<int64_t> mVoxelSliceKeys;
1764 template<
typename Gr
idT>
1767 typename GridT::ValueType isoValue,
1771 if (fs.initSdf(fogGrid, isoValue,
false)) fs.
sweep(nIter);
1772 return fs.sdfGrid();
1775 template<
typename Gr
idT>
1778 typename GridT::ValueType isoValue,
1782 if (fs.initSdf(sdfGrid, isoValue,
true)) fs.
sweep(nIter);
1783 return fs.sdfGrid();
1786 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1787 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1790 const ExtValueT& background,
1791 typename FogGridT::ValueType isoValue,
1794 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1797 if (fs.initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1798 fs.
sweep(nIter,
true);
1799 return fs.extGrid();
1802 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1803 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1806 const ExtValueT &background,
1807 typename SdfGridT::ValueType isoValue,
1810 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1813 if (fs.initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1814 fs.
sweep(nIter,
true);
1815 return fs.extGrid();
1818 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1819 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1822 const ExtValueT &background,
1823 typename FogGridT::ValueType isoValue,
1826 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1829 if (fs.initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1830 fs.
sweep(nIter,
true);
1831 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1834 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1835 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1838 const ExtValueT &background,
1839 typename SdfGridT::ValueType isoValue,
1842 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1845 if (fs.initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1846 fs.
sweep(nIter,
true);
1847 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1850 template<
typename Gr
idT>
1859 if (fs.initDilate(sdfGrid, dilation, nn, mode)) fs.
sweep(nIter);
1860 return fs.sdfGrid();
1863 template<
typename Gr
idT,
typename MaskTreeT>
1867 bool ignoreActiveTiles,
1871 if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.
sweep(nIter);
1872 return fs.sdfGrid();
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean,...
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:415
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:577
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:917
SharedPtr< Grid > Ptr
Definition: Grid.h:579
Definition: Exceptions.h:63
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:564
std::bitset< 6 > intersectionMask(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true a bit-mask where the 6 bits indicates if the center of the stencil intersects the iso-con...
Definition: Stencils.h:188
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Definition: Stencils.h:1232
Templated class to compute the minimum and maximum values.
Definition: Stats.h:31
Definition: LeafManager.h:102
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node's buffer with the nth corresponding auxiliary buffer, where n = bufferIdx.
Definition: LeafManager.h:359
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays,...
Definition: NodeManager.h:531
Definition: ValueAccessor.h:183
void setValueOff(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as inactive.
Definition: ValueAccessor.h:266
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:219
Simple timer for basic profiling.
Definition: CpuTimer.h:67
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
double restart()
Re-start timer.
Definition: CpuTimer.h:150
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
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
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Vec3< double > Vec3d
Definition: Vec3.h:668
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:938
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
Index32 Index
Definition: Types.h:54
@ GRID_LEVEL_SET
Definition: Types.h:337
math::Vec3< Real > Vec3R
Definition: Types.h:72
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:180