27 #ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED 28 #define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED 43 #ifdef BENCHMARK_FAST_SWEEPING 47 #include <tbb/parallel_for.h> 48 #include <tbb/enumerable_thread_specific.h> 49 #include <tbb/task_group.h> 51 #include <type_traits> 55 #include <unordered_map> 103 template<
typename Gr
idT>
106 typename GridT::ValueType isoValue,
136 template<
typename Gr
idT>
139 typename GridT::ValueType isoValue = 0,
192 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
193 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
196 const ExtValueT& background,
197 typename FogGridT::ValueType isoValue,
200 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
250 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
251 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
254 const ExtValueT &background,
255 typename SdfGridT::ValueType isoValue = 0,
258 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
313 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
314 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
317 const ExtValueT &background,
318 typename FogGridT::ValueType isoValue,
321 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
376 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
377 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
380 const ExtValueT &background,
381 typename SdfGridT::ValueType isoValue = 0,
384 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
414 template<
typename Gr
idT>
440 template<
typename Gr
idT,
typename MaskTreeT>
444 bool ignoreActiveTiles =
false,
459 template<
typename SdfGr
idT,
typename ExtValueT =
typename SdfGr
idT::ValueType>
463 "FastSweeping requires SdfGridT to have floating-point values");
465 using SdfValueT =
typename SdfGridT::ValueType;
466 using SdfTreeT =
typename SdfGridT::TreeType;
471 using ExtGridT =
typename SdfGridT::template ValueConverter<ExtValueT>::Type;
472 using ExtTreeT =
typename ExtGridT::TreeType;
476 using SweepMaskTreeT =
typename SdfTreeT::template ValueConverter<ValueMask>::Type;
499 typename SdfGridT::Ptr
sdfGrid() {
return mSdfGrid; }
507 typename ExtGridT::Ptr
extGrid() {
return mExtGrid; }
537 bool initSdf(
const SdfGridT &sdfGrid, SdfValueT isoValue,
bool isInputSdf);
585 template <
typename ExtOpT>
586 bool initExt(
const SdfGridT &sdfGrid,
588 const ExtValueT &background,
592 const typename ExtGridT::ConstPtr extGrid =
nullptr);
615 bool initDilate(
const SdfGridT &sdfGrid,
638 template<
typename MaskTreeT>
639 bool initMask(
const SdfGridT &sdfGrid,
const Grid<MaskTreeT> &mask,
bool ignoreActiveTiles =
false);
653 void sweep(
int nIter = 1,
654 bool finalize =
true);
666 bool isValid()
const {
return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
682 void computeSweepMaskLeafOrigins();
696 static const Coord mOffset[6];
699 typename SdfGridT::Ptr mSdfGrid;
700 typename ExtGridT::Ptr mExtGrid;
701 typename ExtGridT::Ptr mExtGridInput;
702 SweepMaskTreeT mSweepMask;
703 std::vector<Coord> mSweepMaskLeafOrigins;
704 size_t mSweepingVoxelCount, mBoundaryVoxelCount;
712 template <
typename SdfGr
idT,
typename ExtValueT>
717 template <
typename SdfGr
idT,
typename ExtValueT>
719 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(
FastSweepingDomain::
SWEEP_ALL), mIsInputSdf(true)
723 template <
typename SdfGr
idT,
typename ExtValueT>
729 if (mExtGridInput) mExtGridInput.reset();
730 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
735 template <
typename SdfGr
idT,
typename ExtValueT>
741 mSweepMask.voxelizeActiveTiles();
744 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
745 LeafManagerT leafManager(mSweepMask);
747 mSweepMaskLeafOrigins.resize(leafManager.leafCount());
749 auto kernel = [&](
const LeafT& leaf,
size_t leafIdx) {
750 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
753 leafManager.foreach(kernel,
true, 1024);
755 mBoundaryVoxelCount = 0;
758 const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
759 assert( totalCount >= mSweepingVoxelCount );
760 mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
764 template <
typename SdfGr
idT,
typename ExtValueT>
768 mSdfGrid = fogGrid.deepCopy();
771 kernel.
run(isoValue);
775 template <
typename SdfGr
idT,
typename ExtValueT>
776 template <
typename OpT>
782 if (extGrid->transform() != fogGrid.transform())
783 OPENVDB_THROW(
RuntimeError,
"FastSweeping::initExt extension grid input should have the same transform as Fog/SDF grid!");
787 mSdfGrid = fogGrid.deepCopy();
788 mExtGrid = createGrid<ExtGridT>( background );
789 mSweepDirection = mode;
792 mExtGridInput = extGrid->deepCopy();
794 mExtGrid->topologyUnion( *mSdfGrid );
795 InitExt<OpT> kernel(*
this);
796 kernel.run(isoValue, op);
801 template <
typename SdfGr
idT,
typename ExtValueT>
805 mSdfGrid = sdfGrid.deepCopy();
806 mSweepDirection = mode;
808 kernel.
run(dilate, nn);
812 template <
typename SdfGr
idT,
typename ExtValueT>
813 template<
typename MaskTreeT>
817 mSdfGrid = sdfGrid.deepCopy();
819 if (mSdfGrid->transform() != mask.
transform()) {
824 using T =
typename MaskTreeT::template ValueConverter<bool>::Type;
826 tmp->
tree().voxelizeActiveTiles();
827 MaskKernel<T> kernel(*
this);
828 kernel.run(tmp->
tree());
830 if (ignoreActiveTiles || !mask.
tree().hasActiveTiles()) {
831 MaskKernel<MaskTreeT> kernel(*
this);
832 kernel.run(mask.
tree());
834 using T =
typename MaskTreeT::template ValueConverter<ValueMask>::Type;
836 tmp.voxelizeActiveTiles(
true);
837 MaskKernel<T> kernel(*
this);
844 template <
typename SdfGr
idT,
typename ExtValueT>
852 " a non-null reference extension grid input.");
861 std::deque<SweepingKernel> kernels;
862 for (
int i = 0; i < 4; i++) kernels.emplace_back(*
this);
865 #ifdef BENCHMARK_FAST_SWEEPING 870 tbb::task_group tasks;
871 tasks.run([&] { kernels[0].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]+a[2]; }); });
872 tasks.run([&] { kernels[1].computeVoxelSlices([](
const Coord &a){
return a[0]+a[1]-a[2]; }); });
873 tasks.run([&] { kernels[2].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]+a[2]; }); });
874 tasks.run([&] { kernels[3].computeVoxelSlices([](
const Coord &a){
return a[0]-a[1]-a[2]; }); });
877 #ifdef BENCHMARK_FAST_SWEEPING 883 for (
int i = 0; i < nIter; ++i) {
888 #ifdef BENCHMARK_FAST_SWEEPING 892 auto e = kernel.
run(*mSdfGrid);
897 nodeManager.foreachTopDown(op,
true , 1 );
899 #ifdef BENCHMARK_FAST_SWEEPING 900 std::cerr <<
"Min = " << e.min() <<
" Max = " << e.max() << std::endl;
901 timer.
restart(
"Changing asymmetric background value");
905 #ifdef BENCHMARK_FAST_SWEEPING 916 template <
typename SdfGr
idT,
typename ExtValueT>
921 MinMaxKernel() : mMin(
std::numeric_limits<SdfValueT>::
max()), mMax(-mMin), mFltMinExists(false), mFltMaxExists(true) {}
922 MinMaxKernel(
MinMaxKernel& other, tbb::split) : mMin(other.mMin), mMax(other.mMax), mFltMinExists(other.mFltMinExists), mFltMaxExists(other.mFltMaxExists) {}
927 tbb::parallel_reduce(mgr.leafRange(), *
this);
933 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
934 for (
auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
935 const SdfValueT v = *voxelIter;
938 if (v < mMin && !vEqFltMin) mMin = v;
939 if (v > mMax && !vEqFltMax) mMax = v;
940 if (vEqFltMin) mFltMinExists =
true;
941 if (vEqFltMax) mFltMaxExists =
true;
948 if (other.
mMin < mMin) mMin = other.
mMin;
949 if (other.
mMax > mMax) mMax = other.
mMax;
960 template <
typename SdfGr
idT,
typename ExtValueT>
965 void operator()(
typename SdfTreeT::RootNodeType& node,
size_t = 1)
const {
966 for (
auto iter = node.beginValueAll(); iter; ++iter) {
977 template<
typename NodeT>
980 for (
auto iter = node.beginValueAll(); iter; ++iter) {
991 void operator()(
typename SdfTreeT::LeafNodeType& leaf,
size_t = 1)
const 993 for (
auto iter = leaf.beginValueOn(); iter; ++iter) {
1008 template <
typename SdfGr
idT,
typename ExtValueT>
1014 mBackground(parent.mSdfGrid->background())
1016 mSdfGridInput = mParent->mSdfGrid->deepCopy();
1023 #ifdef BENCHMARK_FAST_SWEEPING 1028 #ifdef BENCHMARK_FAST_SWEEPING 1029 timer.
restart(
"Changing background value");
1034 #ifdef BENCHMARK_FAST_SWEEPING 1035 timer.
restart(
"Dilating and updating mgr (parallel)");
1039 const int delta = 5;
1045 #ifdef BENCHMARK_FAST_SWEEPING 1046 timer.
restart(
"Initializing grid and sweep mask");
1049 mParent->mSweepMask.clear();
1050 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1053 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1057 LeafManagerT leafManager(mParent->mSdfGrid->tree());
1059 auto kernel = [&](LeafT& leaf,
size_t ) {
1061 const SdfValueT background = mBackground;
1062 auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
1063 SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
1065 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1066 const SdfValueT
value = *voxelIter;
1067 SdfValueT inputValue;
1068 const Coord ijk = voxelIter.getCoord();
1071 maskLeaf->setValueOff(voxelIter.pos());
1075 voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1078 if (value > 0) voxelIter.setValue(Unknown);
1080 maskLeaf->setValueOff(voxelIter.pos());
1081 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1082 if ( !isInputOn ) voxelIter.setValueOff();
1083 else voxelIter.setValue(inputValue);
1087 if (value < 0) voxelIter.setValue(-Unknown);
1089 maskLeaf->setValueOff(voxelIter.pos());
1090 bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1091 if ( !isInputOn ) voxelIter.setValueOff();
1092 else voxelIter.setValue(inputValue);
1100 leafManager.foreach( kernel );
1103 mParent->computeSweepMaskLeafOrigins();
1105 #ifdef BENCHMARK_FAST_SWEEPING 1118 template <
typename SdfGr
idT,
typename ExtValueT>
1123 mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
1129 mIsoValue = isoValue;
1130 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1131 SdfTreeT &tree = mSdfGrid->tree();
1132 const bool hasActiveTiles = tree.hasActiveTiles();
1134 if (mParent->mIsInputSdf && hasActiveTiles) {
1138 #ifdef BENCHMARK_FAST_SWEEPING 1141 mParent->mSweepMask.clear();
1142 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1146 tbb::parallel_for(mgr.
leafRange(32), *
this);
1150 #ifdef BENCHMARK_FAST_SWEEPING 1151 timer.
restart(
"Initialize tiles - new");
1155 mgr.foreachBottomUp(*
this);
1157 if (hasActiveTiles) tree.voxelizeActiveTiles();
1161 mParent->computeSweepMaskLeafOrigins();
1169 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1170 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1171 SdfValueT* sdf = leafIter.buffer(1).data();
1172 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1173 const SdfValueT
value = *voxelIter;
1174 const bool isAbove = value > isoValue;
1175 if (!voxelIter.isValueOn()) {
1176 sdf[voxelIter.pos()] = isAbove ? above : -above;
1178 const Coord ijk = voxelIter.getCoord();
1179 stencil.
moveTo(ijk, value);
1182 sdf[voxelIter.pos()] = isAbove ? above : -above;
1186 const SdfValueT delta = value - isoValue;
1188 sdf[voxelIter.pos()] = 0;
1191 for (
int i=0; i<6;) {
1194 if (mask.test(i++)) {
1198 if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1208 template<
typename RootOrInternalNodeT>
1212 for (
auto it = node.cbeginValueAll(); it; ++it) {
1213 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1214 v = v > isoValue ? above : -above;
1227 template <
typename SdfGr
idT,
typename ExtValueT>
1228 template <
typename OpT>
1232 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1234 mOpPool(
nullptr), mSdfGrid(parent.mSdfGrid.get()),
1235 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1236 InitExt(
const InitExt&) =
default;
1237 InitExt&
operator=(
const InitExt&) =
delete;
1238 void run(SdfValueT isoValue,
const OpT &opPrototype)
1240 static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::
value,
"Invalid return type of functor");
1245 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1246 mIsoValue = isoValue;
1247 auto &tree1 = mSdfGrid->tree();
1248 auto &tree2 = mExtGrid->tree();
1249 const bool hasActiveTiles = tree1.hasActiveTiles();
1251 if (mParent->mIsInputSdf && hasActiveTiles) {
1255 #ifdef BENCHMARK_FAST_SWEEPING 1259 mParent->mSweepMask.clear();
1260 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1264 OpPoolT opPool(opPrototype);
1268 tbb::parallel_for(mgr.
leafRange(32), *
this);
1272 #ifdef BENCHMARK_FAST_SWEEPING 1273 timer.
restart(
"Initialize tiles");
1277 mgr.foreachBottomUp(*
this);
1279 if (hasActiveTiles) {
1280 #ifdef BENCHMARK_FAST_SWEEPING 1281 timer.
restart(
"Voxelizing active tiles");
1283 tree1.voxelizeActiveTiles();
1284 tree2.voxelizeActiveTiles();
1290 mParent->computeSweepMaskLeafOrigins();
1292 #ifdef BENCHMARK_FAST_SWEEPING 1299 void sumHelper(ExtT& sum2, ExtT ext,
bool update,
const SdfT& )
const {
if (update) sum2 = ext; }
1303 void sumHelper(ExtT& sum2, ExtT ext,
bool ,
const SdfT& d2)
const { sum2 +=
static_cast<ExtValueT
>(d2 * ext); }
1306 ExtT extValHelper(ExtT extSum,
const SdfT& )
const {
return extSum; }
1309 ExtT extValHelper(ExtT extSum,
const SdfT& sdfSum)
const {
return ExtT((SdfT(1) / sdfSum) * extSum); }
1311 void operator()(
const LeafRange& r)
const 1313 ExtAccT acc(mExtGrid->tree());
1317 typename OpPoolT::reference op = mOpPool->local();
1319 const SdfValueT h = mAboveSign*
static_cast<SdfValueT
>(mSdfGrid->voxelSize()[0]);
1320 for (
auto leafIter = r.begin(); leafIter; ++leafIter) {
1321 SdfValueT *sdf = leafIter.buffer(1).data();
1322 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();
1323 for (
auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1324 const SdfValueT
value = *voxelIter;
1325 const bool isAbove = value > isoValue;
1326 if (!voxelIter.isValueOn()) {
1327 sdf[voxelIter.pos()] = isAbove ? above : -above;
1329 const Coord ijk = voxelIter.getCoord();
1330 stencil.
moveTo(ijk, value);
1333 sdf[voxelIter.pos()] = isAbove ? above : -above;
1337 sweepMaskAcc.setValueOff(ijk);
1338 const SdfValueT delta = value - isoValue;
1340 sdf[voxelIter.pos()] = 0;
1341 ext[voxelIter.pos()] = ExtValueT(op(xform.
indexToWorld(ijk)));
1344 ExtValueT sum2 = zeroVal<ExtValueT>();
1350 for (
int n=0, i=0; i<6;) {
1352 if (mask.test(i++)) {
1356 if (mask.test(i++)) {
1363 if (d < std::numeric_limits<SdfValueT>::max()) {
1366 const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1367 static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1368 static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1370 sumHelper(sum2, ExtValueT(op(xform.
indexToWorld(xyz))), d < minD, d2);
1371 if (d < minD) minD = d;
1374 ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1383 template<
typename RootOrInternalNodeT>
1384 void operator()(
const RootOrInternalNodeT& node)
const 1387 for (
auto it = node.cbeginValueAll(); it; ++it) {
1388 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1389 v = v > isoValue ? above : -above;
1397 SdfValueT mIsoValue;
1398 SdfValueT mAboveSign;
1402 template <
typename SdfGr
idT,
typename ExtValueT>
1403 template <
typename MaskTreeT>
1408 mSdfGrid(parent.mSdfGrid.get()) {}
1409 MaskKernel(
const MaskKernel &parent) =
default;
1410 MaskKernel&
operator=(
const MaskKernel&) =
delete;
1412 void run(
const MaskTreeT &mask)
1414 #ifdef BENCHMARK_FAST_SWEEPING 1417 auto &lsTree = mSdfGrid->tree();
1421 #ifdef BENCHMARK_FAST_SWEEPING 1422 timer.
restart(
"Changing background value");
1426 #ifdef BENCHMARK_FAST_SWEEPING 1427 timer.
restart(
"Union with mask");
1429 lsTree.topologyUnion(mask);
1434 #ifdef BENCHMARK_FAST_SWEEPING 1435 timer.
restart(
"Initializing grid and sweep mask");
1438 mParent->mSweepMask.clear();
1439 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1442 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1443 LeafManagerT leafManager(mParent->mSweepMask);
1445 auto kernel = [&](LeafT& leaf,
size_t ) {
1447 SdfAccT acc(mSdfGrid->tree());
1450 SdfValueT *data = acc.
probeLeaf(leaf.origin())->buffer().data();
1451 for (
auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1452 if (
math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1454 voxelIter.setValue(
false);
1458 leafManager.foreach( kernel );
1461 mParent->computeSweepMaskLeafOrigins();
1463 #ifdef BENCHMARK_FAST_SWEEPING 1474 template <
typename SdfGr
idT,
typename ExtValueT>
1482 template<
typename HashOp>
1485 #ifdef BENCHMARK_FAST_SWEEPING 1490 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1493 using LeafT =
typename SweepMaskTreeT::LeafNodeType;
1494 LeafManagerT leafManager(maskTree);
1501 constexpr
int maskOffset = LeafT::DIM * 3;
1502 constexpr
int maskRange = maskOffset * 2;
1505 std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1506 auto kernel1 = [&](
const LeafT& leaf,
size_t leafIdx) {
1507 const size_t leafOffset = leafIdx * maskRange;
1508 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1509 const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1510 leafSliceMasks[leafOffset +
hash(ijk) + maskOffset] = uint8_t(1);
1513 leafManager.foreach( kernel1 );
1518 using ThreadLocalMap = std::unordered_map<int64_t, std::deque<size_t>>;
1519 tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1520 auto kernel2 = [&](
const LeafT& leaf,
size_t leafIdx) {
1521 ThreadLocalMap& map = pool.local();
1522 const Coord& origin = leaf.origin();
1523 const int64_t leafKey =
hash(origin);
1524 const size_t leafOffset = leafIdx * maskRange;
1525 for (
int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1526 if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1527 const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1528 map[voxelSliceKey].emplace_back(leafIdx);
1532 leafManager.foreach( kernel2 );
1537 for (
auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1538 const ThreadLocalMap& map = *poolIt;
1539 for (
const auto& it : map) {
1540 for (
const size_t leafIdx : it.second) {
1541 mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1547 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1548 for (
const auto& it : mVoxelSliceMap) {
1549 mVoxelSliceKeys.push_back(it.first);
1553 auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1554 for (
size_t i = range.begin(); i < range.end(); i++) {
1555 const int64_t key = mVoxelSliceKeys[i];
1556 for (
auto& it : mVoxelSliceMap[key]) {
1557 it.second = std::make_unique<NodeMaskT>();
1561 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1568 auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1569 for (
size_t i = range.begin(); i < range.end(); i++) {
1570 const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1571 LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1572 for (LeafSlice& leafSlice : leafSliceArray) {
1573 const size_t leafIdx = leafSlice.first;
1574 NodeMaskPtrT& nodeMask = leafSlice.second;
1575 const LeafT& leaf = leafManager.leaf(leafIdx);
1576 const Coord& origin = leaf.origin();
1577 const int64_t leafKey =
hash(origin);
1578 for (
auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1579 const Index voxelIdx = voxelIter.pos();
1580 const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1581 const int64_t key = leafKey +
hash(ijk);
1582 if (key == voxelSliceKey) {
1583 nodeMask->setOn(voxelIdx);
1589 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1596 inline static Coord
ijk(
const Coord &p,
int i) {
return p + FastSweeping::mOffset[i]; }
1598 NN(
const SdfAccT &a,
const Coord &p,
int i) : v(math::
Abs(a.getValue(ijk(p,i)))), n(i) {}
1599 inline Coord
operator()(
const Coord &p)
const {
return ijk(p, n); }
1601 inline operator bool()
const {
return v < SdfValueT(1000); }
1606 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 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 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& ,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1622 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& w,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1623 return ExtT(w*(d1.
v*v1 + d2.
v*v2 + d3.
v*v3));
1628 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() :
nullptr;
1629 typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() :
nullptr;
1631 const SdfValueT h =
static_cast<SdfValueT
>(mParent->mSdfGrid->voxelSize()[0]);
1632 const SdfValueT sqrt2h =
math::Sqrt(SdfValueT(2))*h;
1634 const bool isInputSdf = mParent->mIsInputSdf;
1641 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1643 int64_t voxelSliceIndex(0);
1645 auto kernel = [&](
const tbb::blocked_range<size_t>& range) {
1646 using LeafT =
typename SdfGridT::TreeType::LeafNodeType;
1648 SdfAccT acc1(mParent->mSdfGrid->tree());
1649 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ?
new ExtAccT(*tree2) :
nullptr);
1650 auto acc3 = std::unique_ptr<ExtAccT>(tree3 ?
new ExtAccT(*tree3) :
nullptr);
1651 SdfValueT absV, sign, update, D;
1654 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1658 for (
size_t i = range.begin(); i < range.end(); ++i) {
1663 const LeafSlice& leafSlice = leafSliceArray[i];
1664 const size_t leafIdx = leafSlice.first;
1665 const NodeMaskPtrT& nodeMask = leafSlice.second;
1667 const Coord& origin = leafNodeOrigins[leafIdx];
1670 for (
auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1673 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1680 if (!(d1 || d2 || d3))
continue;
1685 SdfValueT &
value =
const_cast<SdfValueT&
>(acc1.getValue(ijk));
1688 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1694 if (d2 < d1) std::swap(d1, d2);
1695 if (d3 < d2) std::swap(d2, d3);
1696 if (d2 < d1) std::swap(d1, d2);
1702 if (update <= d2.
v) {
1703 if (update < absV) {
1704 value = sign * update;
1707 ExtValueT updateExt = acc2->getValue(d1(ijk));
1709 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1710 else updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1713 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1714 else updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1716 acc2->setValue(ijk, updateExt);
1726 if (d2.
v <= sqrt2h + d1.v) {
1727 D = SdfValueT(2) * h * h -
math::Pow2(d1.v - d2.
v);
1728 update = SdfValueT(0.5) * (d1.v + d2.
v + std::sqrt(D));
1729 if (update > d2.
v && update <= d3.
v) {
1730 if (update < absV) {
1731 value = sign * update;
1736 const SdfValueT w = SdfValueT(1)/(d1.v+d2.
v);
1737 const ExtValueT v1 = acc2->getValue(d1(ijk));
1738 const ExtValueT v2 = acc2->getValue(d2(ijk));
1739 const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1741 ExtValueT updateExt = extVal;
1743 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1744 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1747 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1748 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1750 acc2->setValue(ijk, updateExt);
1761 const SdfValueT d123 = d1.v + d2.
v + d3.
v;
1762 D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.
v*d2.
v + d3.
v*d3.
v - h * h);
1763 if (D >= SdfValueT(0)) {
1764 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));
1766 if (update < absV) {
1767 value = sign * update;
1773 const SdfValueT w = SdfValueT(1)/(d1.v+d2.
v+d3.
v);
1774 const ExtValueT v1 = acc2->getValue(d1(ijk));
1775 const ExtValueT v2 = acc2->getValue(d2(ijk));
1776 const ExtValueT v3 = acc2->getValue(d3(ijk));
1777 const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1779 ExtValueT updateExt = extVal;
1781 if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1782 else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1785 if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1786 else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1788 acc2->setValue(ijk, updateExt);
1796 #ifdef BENCHMARK_FAST_SWEEPING 1800 for (
size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1801 voxelSliceIndex = mVoxelSliceKeys[i];
1802 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1805 #ifdef BENCHMARK_FAST_SWEEPING 1806 timer.
restart(
"Backward sweeps");
1808 for (
size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1809 voxelSliceIndex = mVoxelSliceKeys[i-1];
1810 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1813 #ifdef BENCHMARK_FAST_SWEEPING 1819 using NodeMaskT =
typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1820 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1823 using LeafSlice = std::pair<size_t, NodeMaskPtrT>;
1824 using LeafSliceArray = std::deque<LeafSlice>;
1825 using VoxelSliceMap = std::map<int64_t, LeafSliceArray>;
1829 VoxelSliceMap mVoxelSliceMap;
1830 std::vector<int64_t> mVoxelSliceKeys;
1835 template<
typename Gr
idT>
1838 typename GridT::ValueType isoValue,
1842 if (fs.
initSdf(fogGrid, isoValue,
false)) fs.
sweep(nIter);
1846 template<
typename Gr
idT>
1849 typename GridT::ValueType isoValue,
1853 if (fs.
initSdf(sdfGrid, isoValue,
true)) fs.
sweep(nIter);
1857 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1858 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1861 const ExtValueT& background,
1862 typename FogGridT::ValueType isoValue,
1865 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1868 if (fs.
initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1869 fs.
sweep(nIter,
true);
1873 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1874 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1877 const ExtValueT &background,
1878 typename SdfGridT::ValueType isoValue,
1881 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1884 if (fs.
initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1885 fs.
sweep(nIter,
true);
1889 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1890 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1893 const ExtValueT &background,
1894 typename FogGridT::ValueType isoValue,
1897 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1900 if (fs.
initExt(fogGrid, op, background, isoValue,
false, mode, extGrid))
1901 fs.
sweep(nIter,
true);
1905 template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1906 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1909 const ExtValueT &background,
1910 typename SdfGridT::ValueType isoValue,
1913 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr
extGrid)
1916 if (fs.
initExt(sdfGrid, op, background, isoValue,
true, mode, extGrid))
1917 fs.
sweep(nIter,
true);
1921 template<
typename Gr
idT>
1934 template<
typename Gr
idT,
typename MaskTreeT>
1938 bool ignoreActiveTiles,
1942 if (fs.
initMask(sdfGrid, mask, ignoreActiveTiles)) fs.
sweep(nIter);
1952 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 1954 #ifdef OPENVDB_INSTANTIATE_FASTSWEEPING 1958 #define _FUNCTION(TreeT) \ 1959 Grid<TreeT>::Ptr fogToSdf(const Grid<TreeT>&, TreeT::ValueType, int) 1963 #define _FUNCTION(TreeT) \ 1964 Grid<TreeT>::Ptr sdfToSdf(const Grid<TreeT>&, TreeT::ValueType, int) 1968 #define _FUNCTION(TreeT) \ 1969 Grid<TreeT>::Ptr dilateSdf(const Grid<TreeT>&, int, NearestNeighbors, int, FastSweepingDomain) 1973 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 1980 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED To facilitate threading over the nodes of a tree, cache node pointers in linear arrays, one for each level of the tree.
Definition: NodeManager.h:30
Index32 Index
Definition: Types.h:54
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Definition: Stencils.h:1232
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Templated class to compute the minimum and maximum values.
Definition: Stats.h:30
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
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:282
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
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
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.
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
Definition: ValueAccessor.h:395
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
Definition: Exceptions.h:63
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:908
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:48
Implementation of morphological dilation and erosion.
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:931
Definition: Exceptions.h:13
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
ValueT value
Definition: GridBuilder.h:1290
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:13
Definition: LeafManager.h:101
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:411
Definition: ValueAccessor.h:190
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
double restart()
Re-start timer.
Definition: CpuTimer.h:150
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:98
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
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
Simple timer for basic profiling.
Definition: CpuTimer.h:66
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
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:189