27 #ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED 28 #define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED 44 #ifdef BENCHMARK_FAST_SWEEPING 48 #include <tbb/parallel_for.h> 49 #include <tbb/enumerable_thread_specific.h> 50 #include <tbb/task_group.h> 52 #include <type_traits> 56 #include <unordered_map> 104 template<
typename Gr
idT>
107 typename GridT::ValueType isoValue,
137 template<
typename Gr
idT>
140 typename GridT::ValueType isoValue = 0,
193 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
194 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
197 const ExtValueT& background,
198 typename FogGridT::ValueType isoValue,
201 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
251 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
252 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
255 const ExtValueT &background,
256 typename SdfGridT::ValueType isoValue = 0,
259 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
314 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
315 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
318 const ExtValueT &background,
319 typename FogGridT::ValueType isoValue,
322 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
377 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
378 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
381 const ExtValueT &background,
382 typename SdfGridT::ValueType isoValue = 0,
385 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
415 template<
typename Gr
idT>
441 template<
typename Gr
idT,
typename MaskTreeT>
445 bool ignoreActiveTiles =
false,
460 template<
typename SdfGr
idT,
typename ExtValueT =
typename SdfGr
idT::ValueType>
463 static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
464 "FastSweeping requires SdfGridT to have floating-point values");
466 using SdfValueT =
typename SdfGridT::ValueType;
467 using SdfTreeT =
typename SdfGridT::TreeType;
472 using ExtGridT =
typename SdfGridT::template ValueConverter<ExtValueT>::Type;
473 using ExtTreeT =
typename ExtGridT::TreeType;
477 using SweepMaskTreeT =
typename SdfTreeT::template ValueConverter<ValueMask>::Type;
500 typename SdfGridT::Ptr
sdfGrid() {
return mSdfGrid; }
508 typename ExtGridT::Ptr
extGrid() {
return mExtGrid; }
538 bool initSdf(
const SdfGridT &sdfGrid, SdfValueT isoValue,
bool isInputSdf);
586 template <
typename ExtOpT>
587 bool initExt(
const SdfGridT &sdfGrid,
589 const ExtValueT &background,
593 const typename ExtGridT::ConstPtr extGrid =
nullptr);
616 bool initDilate(
const SdfGridT &sdfGrid,
639 template<
typename MaskTreeT>
640 bool initMask(
const SdfGridT &sdfGrid,
const Grid<MaskTreeT> &mask,
bool ignoreActiveTiles =
false);
654 void sweep(
int nIter = 1,
655 bool finalize =
true);
667 bool isValid()
const {
return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
683 void computeSweepMaskLeafOrigins();
697 static const Coord mOffset[6];
700 typename SdfGridT::Ptr mSdfGrid;
701 typename ExtGridT::Ptr mExtGrid;
702 typename ExtGridT::Ptr mExtGridInput;
703 SweepMaskTreeT mSweepMask;
704 std::vector<Coord> mSweepMaskLeafOrigins;
705 size_t mSweepingVoxelCount, mBoundaryVoxelCount;
713 template <
typename SdfGr
idT,
typename ExtValueT>
718 template <
typename SdfGr
idT,
typename ExtValueT>
720 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(
FastSweepingDomain::
SWEEP_ALL), mIsInputSdf(true)
724 template <
typename SdfGr
idT,
typename ExtValueT>
730 if (mExtGridInput) mExtGridInput.reset();
731 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
736 template <
typename SdfGr
idT,
typename ExtValueT>
742 mSweepMask.voxelizeActiveTiles();
745 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
746 LeafManagerT leafManager(mSweepMask);
748 mSweepMaskLeafOrigins.resize(leafManager.leafCount());
750 auto kernel = [&](
const LeafT& leaf,
size_t leafIdx) {
751 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
754 leafManager.foreach(kernel,
true, 1024);
756 mBoundaryVoxelCount = 0;
759 const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
761 mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
765 template <
typename SdfGr
idT,
typename ExtValueT>
769 mSdfGrid = fogGrid.deepCopy();
772 kernel.
run(isoValue);
776 template <
typename SdfGr
idT,
typename ExtValueT>
777 template <
typename OpT>
783 if (extGrid->transform() != fogGrid.transform())
784 OPENVDB_THROW(
RuntimeError,
"FastSweeping::initExt extension grid input should have the same transform as Fog/SDF grid!");
788 mSdfGrid = fogGrid.deepCopy();
789 mExtGrid = createGrid<ExtGridT>( background );
790 mSweepDirection = mode;
793 mExtGridInput = extGrid->deepCopy();
795 mExtGrid->topologyUnion( *mSdfGrid );
796 InitExt<OpT> kernel(*
this);
797 kernel.run(isoValue, op);
802 template <
typename SdfGr
idT,
typename ExtValueT>
806 mSdfGrid = sdfGrid.deepCopy();
807 mSweepDirection = mode;
809 kernel.
run(dilate, nn);
813 template <
typename SdfGr
idT,
typename ExtValueT>
814 template<
typename MaskTreeT>
818 mSdfGrid = sdfGrid.deepCopy();
820 if (mSdfGrid->transform() != mask.
transform()) {
825 using T =
typename MaskTreeT::template ValueConverter<bool>::Type;
827 tmp->
tree().voxelizeActiveTiles();
828 MaskKernel<T> kernel(*
this);
829 kernel.run(tmp->
tree());
831 if (ignoreActiveTiles || !mask.
tree().hasActiveTiles()) {
832 MaskKernel<MaskTreeT> kernel(*
this);
833 kernel.run(mask.
tree());
835 using T =
typename MaskTreeT::template ValueConverter<ValueMask>::Type;
837 tmp.voxelizeActiveTiles(
true);
838 MaskKernel<T> kernel(*
this);
845 template <
typename SdfGr
idT,
typename ExtValueT>
853 " a non-null reference extension grid input.");
862 std::deque<SweepingKernel> kernels;
863 for (
int i = 0; i < 4; i++) kernels.emplace_back(*
this);
866 #ifdef BENCHMARK_FAST_SWEEPING 871 tbb::task_group tasks;
872 tasks.run([&] { kernels[0].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]+a[2]; }); });
873 tasks.run([&] { kernels[1].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]-a[2]; }); });
874 tasks.run([&] { kernels[2].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]+a[2]; }); });
875 tasks.run([&] { kernels[3].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]-a[2]; }); });
878 #ifdef BENCHMARK_FAST_SWEEPING 884 for (
int i = 0; i < nIter; ++i) {
889 #ifdef BENCHMARK_FAST_SWEEPING 893 auto e = kernel.
run(*mSdfGrid);
898 nodeManager.foreachTopDown(op,
true , 1 );
900 #ifdef BENCHMARK_FAST_SWEEPING 901 std::cerr <<
"Min = " << e.min() <<
" Max = " << e.max() << std::endl;
902 timer.
restart(
"Changing asymmetric background value");
906 #ifdef BENCHMARK_FAST_SWEEPING 917 template <
typename SdfGr
idT,
typename ExtValueT>
922 MinMaxKernel() : mMin(
std::numeric_limits<SdfValueT>::
max()), mMax(-mMin), mFltMinExists(false), mFltMaxExists(true) {}
928 tbb::parallel_reduce(mgr.leafRange(), *
this);
934 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
935 for (
auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
936 const SdfValueT v = *voxelIter;
939 if (v < mMin && !vEqFltMin) mMin = v;
940 if (v > mMax && !vEqFltMax) mMax = v;
941 if (vEqFltMin) mFltMinExists =
true;
942 if (vEqFltMax) mFltMaxExists =
true;
949 if (other.
mMin < mMin) mMin = other.
mMin;
950 if (other.
mMax > mMax) mMax = other.
mMax;
961 template <
typename SdfGr
idT,
typename ExtValueT>
966 void operator()(
typename SdfTreeT::RootNodeType& node,
size_t = 1)
const {
967 for (
auto iter = node.beginValueAll(); iter; ++iter) {
978 template<
typename NodeT>
981 for (
auto iter = node.beginValueAll(); iter; ++iter) {
992 void operator()(
typename SdfTreeT::LeafNodeType& leaf,
size_t = 1)
const 994 for (
auto iter = leaf.beginValueOn(); iter; ++iter) {
1009 template <
typename SdfGr
idT,
typename ExtValueT>
1015 mBackground(parent.mSdfGrid->background())
1017 mSdfGridInput = mParent->mSdfGrid->deepCopy();
1024 #ifdef BENCHMARK_FAST_SWEEPING 1029 #ifdef BENCHMARK_FAST_SWEEPING 1030 timer.
restart(
"Changing background value");
1035 #ifdef BENCHMARK_FAST_SWEEPING 1036 timer.
restart(
"Dilating and updating mgr (parallel)");
1040 const int delta = 5;
1046 #ifdef BENCHMARK_FAST_SWEEPING 1047 timer.
restart(
"Initializing grid and sweep mask");
1050 mParent->mSweepMask.clear();
1051 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1054 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1058 LeafManagerT leafManager(mParent->mSdfGrid->tree());
1060 auto kernel = [&](LeafT& leaf,
size_t ) {
1062 const SdfValueT background = mBackground;
1063 auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
1064 SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
1066 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1067 const SdfValueT value = *voxelIter;
1068 SdfValueT inputValue;
1069 const Coord ijk = voxelIter.getCoord();
1072 maskLeaf->setValueOff(voxelIter.pos());
1076 voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1079 if (value > 0) voxelIter.setValue(Unknown);
1081 maskLeaf->setValueOff(voxelIter.pos());
1082 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1083 if ( !isInputOn ) voxelIter.setValueOff();
1084 else voxelIter.setValue(inputValue);
1088 if (value < 0) voxelIter.setValue(-Unknown);
1090 maskLeaf->setValueOff(voxelIter.pos());
1091 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1092 if ( !isInputOn ) voxelIter.setValueOff();
1093 else voxelIter.setValue(inputValue);
1101 leafManager.foreach( kernel );
1104 mParent->computeSweepMaskLeafOrigins();
1106 #ifdef BENCHMARK_FAST_SWEEPING 1119 template <
typename SdfGr
idT,
typename ExtValueT>
1124 mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
1130 mIsoValue = isoValue;
1131 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1132 SdfTreeT &tree = mSdfGrid->tree();
1133 const bool hasActiveTiles = tree.hasActiveTiles();
1135 if (mParent->mIsInputSdf && hasActiveTiles) {
1139 #ifdef BENCHMARK_FAST_SWEEPING 1142 mParent->mSweepMask.clear();
1143 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1147 tbb::parallel_for(mgr.
leafRange(32), *
this);
1151 #ifdef BENCHMARK_FAST_SWEEPING 1152 timer.
restart(
"Initialize tiles - new");
1156 mgr.foreachBottomUp(*
this);
1158 if (hasActiveTiles) tree.voxelizeActiveTiles();
1162 mParent->computeSweepMaskLeafOrigins();
1170 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1171 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1172 SdfValueT* sdf = leafIter.buffer(1).data();
1173 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1174 const SdfValueT value = *voxelIter;
1175 const bool isAbove = value > isoValue;
1176 if (!voxelIter.isValueOn()) {
1177 sdf[voxelIter.pos()] = isAbove ? above : -above;
1179 const Coord ijk = voxelIter.getCoord();
1180 stencil.
moveTo(ijk, value);
1183 sdf[voxelIter.pos()] = isAbove ? above : -above;
1187 const SdfValueT delta = value - isoValue;
1189 sdf[voxelIter.pos()] = 0;
1192 for (
int i=0; i<6;) {
1195 if (mask.test(i++)) {
1199 if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1209 template<
typename RootOrInternalNodeT>
1213 for (
auto it = node.cbeginValueAll(); it; ++it) {
1214 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1215 v = v > isoValue ? above : -above;
1228 template <
typename SdfGr
idT,
typename ExtValueT>
1229 template <
typename OpT>
1233 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1235 mOpPool(
nullptr), mSdfGrid(parent.mSdfGrid.get()),
1236 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1237 InitExt(
const InitExt&) =
default;
1238 InitExt&
operator=(
const InitExt&) =
delete;
1239 void run(SdfValueT isoValue,
const OpT &opPrototype)
1241 static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::value,
"Invalid return type of functor");
1246 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1247 mIsoValue = isoValue;
1248 auto &tree1 = mSdfGrid->tree();
1249 auto &tree2 = mExtGrid->tree();
1250 const bool hasActiveTiles = tree1.hasActiveTiles();
1252 if (mParent->mIsInputSdf && hasActiveTiles) {
1256 #ifdef BENCHMARK_FAST_SWEEPING 1260 mParent->mSweepMask.clear();
1261 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1265 OpPoolT opPool(opPrototype);
1269 tbb::parallel_for(mgr.
leafRange(32), *
this);
1273 #ifdef BENCHMARK_FAST_SWEEPING 1274 timer.
restart(
"Initialize tiles");
1278 mgr.foreachBottomUp(*
this);
1280 if (hasActiveTiles) {
1281 #ifdef BENCHMARK_FAST_SWEEPING 1282 timer.
restart(
"Voxelizing active tiles");
1284 tree1.voxelizeActiveTiles();
1285 tree2.voxelizeActiveTiles();
1291 mParent->computeSweepMaskLeafOrigins();
1293 #ifdef BENCHMARK_FAST_SWEEPING 1299 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1300 void sumHelper(ExtT& sum2, ExtT ext,
bool update,
const SdfT& )
const {
if (update) sum2 = ext; }
1303 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1304 void sumHelper(ExtT& sum2, ExtT ext,
bool ,
const SdfT& d2)
const { sum2 +=
static_cast<ExtValueT
>(d2 * ext); }
1306 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1307 ExtT extValHelper(ExtT extSum,
const SdfT& )
const {
return extSum; }
1309 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1310 ExtT extValHelper(ExtT extSum,
const SdfT& sdfSum)
const {
return ExtT((SdfT(1) / sdfSum) * extSum); }
1312 void operator()(
const LeafRange& r)
const 1314 ExtAccT acc(mExtGrid->tree());
1318 typename OpPoolT::reference op = mOpPool->local();
1320 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1321 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1322 SdfValueT *sdf = leafIter.buffer(1).data();
1323 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();
1324 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1325 const SdfValueT value = *voxelIter;
1326 const bool isAbove = value > isoValue;
1327 if (!voxelIter.isValueOn()) {
1328 sdf[voxelIter.pos()] = isAbove ? above : -above;
1330 const Coord ijk = voxelIter.getCoord();
1331 stencil.
moveTo(ijk, value);
1334 sdf[voxelIter.pos()] = isAbove ? above : -above;
1338 sweepMaskAcc.setValueOff(ijk);
1339 const SdfValueT delta = value - isoValue;
1341 sdf[voxelIter.pos()] = 0;
1345 ExtValueT sum2 = zeroVal<ExtValueT>();
1351 for (
int n=0, i=0; i<6;) {
1353 if (mask.test(i++)) {
1357 if (mask.test(i++)) {
1364 if (d < std::numeric_limits<SdfValueT>::max()) {
1367 const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1368 static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1369 static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1371 sumHelper(sum2, ExtValueT(
op(xform.
indexToWorld(xyz))), d < minD, d2);
1372 if (d < minD) minD = d;
1375 ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1384 template<
typename RootOrInternalNodeT>
1385 void operator()(
const RootOrInternalNodeT& node)
const 1388 for (
auto it = node.cbeginValueAll(); it; ++it) {
1389 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1390 v = v > isoValue ? above : -above;
1398 SdfValueT mIsoValue;
1399 SdfValueT mAboveSign;
1403 template <
typename SdfGr
idT,
typename ExtValueT>
1404 template <
typename MaskTreeT>
1409 mSdfGrid(parent.mSdfGrid.get()) {}
1410 MaskKernel(
const MaskKernel &parent) =
default;
1411 MaskKernel&
operator=(
const MaskKernel&) =
delete;
1413 void run(
const MaskTreeT &mask)
1415 #ifdef BENCHMARK_FAST_SWEEPING 1418 auto &lsTree = mSdfGrid->tree();
1422 #ifdef BENCHMARK_FAST_SWEEPING 1423 timer.
restart(
"Changing background value");
1427 #ifdef BENCHMARK_FAST_SWEEPING 1428 timer.
restart(
"Union with mask");
1430 lsTree.topologyUnion(mask);
1435 #ifdef BENCHMARK_FAST_SWEEPING 1436 timer.
restart(
"Initializing grid and sweep mask");
1439 mParent->mSweepMask.clear();
1440 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1443 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1444 LeafManagerT leafManager(mParent->mSweepMask);
1446 auto kernel = [&](LeafT& leaf,
size_t ) {
1448 SdfAccT acc(mSdfGrid->tree());
1451 SdfValueT *data = acc.
probeLeaf(leaf.origin())->buffer().data();
1452 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1453 if (
math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1455 voxelIter.setValue(
false);
1459 leafManager.foreach( kernel );
1462 mParent->computeSweepMaskLeafOrigins();
1464 #ifdef BENCHMARK_FAST_SWEEPING 1475 template <
typename SdfGr
idT,
typename ExtValueT>
1483 template<
typename HashOp>
1486 #ifdef BENCHMARK_FAST_SWEEPING 1491 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1494 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1495 LeafManagerT leafManager(maskTree);
1502 constexpr
int maskOffset = LeafT::DIM * 3;
1503 constexpr
int maskRange = maskOffset * 2;
1506 std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1507 auto kernel1 = [&](
const LeafT& leaf,
size_t leafIdx) {
1508 const size_t leafOffset = leafIdx * maskRange;
1509 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1510 const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1511 leafSliceMasks[leafOffset +
hash(ijk) + maskOffset] = uint8_t(1);
1514 leafManager.foreach( kernel1 );
1519 using ThreadLocalMap = std::unordered_map<int64_t, std::deque<size_t>>;
1520 tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1521 auto kernel2 = [&](
const LeafT& leaf,
size_t leafIdx) {
1522 ThreadLocalMap& map = pool.local();
1523 const Coord& origin = leaf.origin();
1524 const int64_t leafKey =
hash(origin);
1525 const size_t leafOffset = leafIdx * maskRange;
1526 for (
int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1527 if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1528 const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1529 map[voxelSliceKey].emplace_back(leafIdx);
1533 leafManager.foreach( kernel2 );
1538 for (
auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1539 const ThreadLocalMap& map = *poolIt;
1540 for (
const auto& it : map) {
1541 for (
const size_t leafIdx : it.second) {
1542 mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1548 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1549 for (
const auto& it : mVoxelSliceMap) {
1550 mVoxelSliceKeys.push_back(it.first);
1554 auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1555 for (
size_t i = range.begin(); i < range.end(); i++) {
1556 const int64_t key = mVoxelSliceKeys[i];
1557 for (
auto& it : mVoxelSliceMap[key]) {
1558 it.second = std::make_unique<NodeMaskT>();
1562 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1569 auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1570 for (
size_t i = range.begin(); i < range.end(); i++) {
1571 const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1572 LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1573 for (LeafSlice& leafSlice : leafSliceArray) {
1574 const size_t leafIdx = leafSlice.first;
1575 NodeMaskPtrT& nodeMask = leafSlice.second;
1576 const LeafT& leaf = leafManager.leaf(leafIdx);
1577 const Coord& origin = leaf.origin();
1578 const int64_t leafKey =
hash(origin);
1579 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1580 const Index voxelIdx = voxelIter.pos();
1581 const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1582 const int64_t key = leafKey +
hash(ijk);
1583 if (key == voxelSliceKey) {
1584 nodeMask->setOn(voxelIdx);
1590 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1597 inline static Coord
ijk(
const Coord &p,
int i) {
return p + FastSweeping::mOffset[i]; }
1599 NN(
const SdfAccT &a,
const Coord &p,
int i) : v(math::
Abs(a.getValue(ijk(p,i)))), n(i) {}
1600 inline Coord
operator()(
const Coord &p)
const {
return ijk(p, n); }
1602 inline operator bool()
const {
return v < SdfValueT(1000); }
1606 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1607 ExtT
twoNghbr(
const NN& d1,
const NN& d2,
const SdfT& ,
const ExtT& v1,
const ExtT& v2)
const {
return d1.
v < d2.
v ? v1 : v2; }
1610 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1611 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)); }
1614 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1615 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& ,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1622 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1623 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& w,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1624 return ExtT(w*(d1.
v*v1 + d2.
v*v2 + d3.
v*v3));
1629 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() :
nullptr;
1630 typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() :
nullptr;
1632 const SdfValueT h =
static_cast<SdfValueT
>(mParent->mSdfGrid->voxelSize()[0]);
1633 const SdfValueT sqrt2h =
math::Sqrt(SdfValueT(2))*h;
1635 const bool isInputSdf = mParent->mIsInputSdf;
1642 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1644 int64_t voxelSliceIndex(0);
1646 auto kernel = [&](
const tbb::blocked_range<size_t>& range) {
1647 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1649 SdfAccT acc1(mParent->mSdfGrid->tree());
1650 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ?
new ExtAccT(*tree2) :
nullptr);
1651 auto acc3 = std::unique_ptr<ExtAccT>(tree3 ?
new ExtAccT(*tree3) :
nullptr);
1652 SdfValueT absV, sign, update, D;
1655 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1659 for (
size_t i = range.begin(); i < range.end(); ++i) {
1664 const LeafSlice& leafSlice = leafSliceArray[i];
1665 const size_t leafIdx = leafSlice.first;
1666 const NodeMaskPtrT& nodeMask = leafSlice.second;
1668 const Coord& origin = leafNodeOrigins[leafIdx];
1671 for (
auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1674 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1681 if (!(d1 || d2 || d3))
continue;
1686 SdfValueT &value =
const_cast<SdfValueT&
>(acc1.getValue(ijk));
1689 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1695 if (d2 < d1) std::swap(d1, d2);
1696 if (d3 < d2) std::swap(d2, d3);
1697 if (d2 < d1) std::swap(d1, d2);
1703 if (update <= d2.
v) {
1704 if (update < absV) {
1705 value = sign * update;
1708 ExtValueT updateExt = acc2->getValue(d1(ijk));
1710 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1711 else updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1714 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1715 else updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1717 acc2->setValue(ijk, updateExt);
1727 if (d2.
v <= sqrt2h + d1.v) {
1728 D = SdfValueT(2) * h * h -
math::Pow2(d1.v - d2.
v);
1729 update = SdfValueT(0.5) * (d1.v + d2.
v + std::sqrt(D));
1730 if (update > d2.
v && update <= d3.
v) {
1731 if (update < absV) {
1732 value = sign * update;
1737 const SdfValueT w = SdfValueT(1)/(d1.v+d2.
v);
1738 const ExtValueT v1 = acc2->getValue(d1(ijk));
1739 const ExtValueT v2 = acc2->getValue(d2(ijk));
1740 const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1742 ExtValueT updateExt = extVal;
1744 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1745 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1748 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1749 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1751 acc2->setValue(ijk, updateExt);
1762 const SdfValueT d123 = d1.v + d2.
v + d3.
v;
1763 D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.
v*d2.
v + d3.
v*d3.
v - h * h);
1764 if (D >= SdfValueT(0)) {
1765 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));
1767 if (update < absV) {
1768 value = sign * update;
1774 const SdfValueT w = SdfValueT(1)/(d1.v+d2.
v+d3.
v);
1775 const ExtValueT v1 = acc2->getValue(d1(ijk));
1776 const ExtValueT v2 = acc2->getValue(d2(ijk));
1777 const ExtValueT v3 = acc2->getValue(d3(ijk));
1778 const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1780 ExtValueT updateExt = extVal;
1782 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1783 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1786 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1787 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1789 acc2->setValue(ijk, updateExt);
1797 #ifdef BENCHMARK_FAST_SWEEPING 1801 for (
size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1802 voxelSliceIndex = mVoxelSliceKeys[i];
1803 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1806 #ifdef BENCHMARK_FAST_SWEEPING 1807 timer.
restart(
"Backward sweeps");
1809 for (
size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1810 voxelSliceIndex = mVoxelSliceKeys[i-1];
1811 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1814 #ifdef BENCHMARK_FAST_SWEEPING 1820 using NodeMaskT =
typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1821 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1824 using LeafSlice = std::pair<size_t, NodeMaskPtrT>;
1825 using LeafSliceArray = std::deque<LeafSlice>;
1826 using VoxelSliceMap = std::map<int64_t, LeafSliceArray>;
1830 VoxelSliceMap mVoxelSliceMap;
1831 std::vector<int64_t> mVoxelSliceKeys;
1836 template<
typename Gr
idT>
1839 typename GridT::ValueType isoValue,
1843 if (fs.
initSdf(fogGrid, isoValue,
false)) fs.
sweep(nIter);
1847 template<
typename Gr
idT>
1850 typename GridT::ValueType isoValue,
1854 if (fs.
initSdf(sdfGrid, isoValue,
true)) fs.
sweep(nIter);
1858 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1859 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1862 const ExtValueT& background,
1863 typename FogGridT::ValueType isoValue,
1866 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1869 if (fs.
initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1870 fs.
sweep(nIter,
true);
1874 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1875 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1878 const ExtValueT &background,
1879 typename SdfGridT::ValueType isoValue,
1882 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1885 if (fs.
initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1886 fs.
sweep(nIter,
true);
1890 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1891 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1894 const ExtValueT &background,
1895 typename FogGridT::ValueType isoValue,
1898 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1901 if (fs.
initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1902 fs.
sweep(nIter,
true);
1906 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1907 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1910 const ExtValueT &background,
1911 typename SdfGridT::ValueType isoValue,
1914 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1917 if (fs.
initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1918 fs.
sweep(nIter,
true);
1922 template<
typename Gr
idT>
1935 template<
typename Gr
idT,
typename MaskTreeT>
1939 bool ignoreActiveTiles,
1943 if (fs.
initMask(sdfGrid, mask, ignoreActiveTiles)) fs.
sweep(nIter);
1953 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 1955 #ifdef OPENVDB_INSTANTIATE_FASTSWEEPING 1959 #define _FUNCTION(TreeT) \ 1960 Grid<TreeT>::Ptr fogToSdf(const Grid<TreeT>&, TreeT::ValueType, int) 1964 #define _FUNCTION(TreeT) \ 1965 Grid<TreeT>::Ptr sdfToSdf(const Grid<TreeT>&, TreeT::ValueType, int) 1969 #define _FUNCTION(TreeT) \ 1970 Grid<TreeT>::Ptr dilateSdf(const Grid<TreeT>&, int, NearestNeighbors, int, FastSweepingDomain) 1974 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 1981 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
Simple timer for basic profiling.
Definition: CpuTimer.h:66
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:360
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:31
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
void setValueOff(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as inactive.
Definition: ValueAccessor.h:607
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Coord Abs(const Coord &xyz)
Definition: Coord.h:518
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Definition: Stencils.h:1231
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
Templated class to compute the minimum and maximum values.
Definition: Stats.h:31
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
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
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Index32 Index
Definition: Types.h:54
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:908
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:349
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:346
Implementation of morphological dilation and erosion.
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Definition: Exceptions.h:13
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:14
OPENVDB_AX_API void run(const char *ax, openvdb::GridBase &grid, const AttributeBindings &bindings={})
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
Definition: LeafManager.h:102
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:162
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:411
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
Definition: ValueAccessor.h:836
double restart()
Re-start timer.
Definition: CpuTimer.h:150
NodeManager produces linear arrays of all tree nodes allowing for efficient threading and bottom-up p...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:930
Definition: Exceptions.h:63