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>
136 template<
typename Gr
idT>
192 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
193 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
196 const ExtValueT& background,
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,
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,
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,
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");
466 using SdfTreeT =
typename SdfGridT::TreeType;
472 using ExtTreeT =
typename ExtGridT::TreeType;
499 typename SdfGridT::Ptr
sdfGrid() {
return mSdfGrid; }
507 typename ExtGridT::Ptr
extGrid() {
return mExtGrid; }
585 template <
typename ExtOpT>
588 const ExtValueT &background,
592 const typename ExtGridT::ConstPtr
extGrid =
nullptr);
638 template<
typename MaskTreeT>
653 void sweep(
int nIter = 1,
654 bool finalize =
true);
666 bool isValid()
const {
return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
682 void computeSweepMaskLeafOrigins();
692 struct PruneMinMaxFltKernel;
693 struct SweepingKernel;
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());
748 std::atomic<size_t> sweepingVoxelCount{0};
749 auto kernel = [&](
const LeafT& leaf,
size_t leafIdx) {
750 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
751 sweepingVoxelCount += leaf.onVoxelCount();
753 leafManager.foreach(kernel,
true, 1024);
755 mBoundaryVoxelCount = 0;
756 mSweepingVoxelCount = sweepingVoxelCount;
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();
769 mIsInputSdf = isInputSdf;
771 kernel.
run(isoValue);
772 return this->isValid();
775 template <
typename SdfGr
idT,
typename ExtValueT>
776 template <
typename OpT>
781 OPENVDB_THROW(RuntimeError,
"FastSweeping::initExt Calling initExt with mode != SWEEP_ALL requires an extension grid!");
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;
790 mIsInputSdf = isInputSdf;
792 mExtGridInput = extGrid->deepCopy();
794 mExtGrid->topologyUnion( *mSdfGrid );
795 InitExt<OpT> kernel(*
this);
796 kernel.run(isoValue, op);
797 return this->isValid();
801 template <
typename SdfGr
idT,
typename ExtValueT>
805 mSdfGrid = sdfGrid.deepCopy();
806 mSweepDirection =
mode;
808 kernel.
run(dilate, nn);
809 return this->isValid();
812 template <
typename SdfGr
idT,
typename ExtValueT>
813 template<
typename MaskTreeT>
817 mSdfGrid = sdfGrid.deepCopy();
819 if (mSdfGrid->transform() != mask.
transform()) {
820 OPENVDB_THROW(RuntimeError,
"FastSweeping: Mask not aligned with the grid!");
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());
836 tmp.voxelizeActiveTiles(
true);
837 MaskKernel<T> kernel(*
this);
841 return this->isValid();
844 template <
typename SdfGr
idT,
typename ExtValueT>
848 OPENVDB_THROW(RuntimeError,
"FastSweeping::sweep called before initialization!");
851 OPENVDB_THROW(RuntimeError,
"FastSweeping: Trying to extend a field in one direction needs"
852 " a non-null reference extension grid input.");
854 if (this->boundaryVoxelCount() == 0) {
855 OPENVDB_THROW(RuntimeError,
"FastSweeping: No boundary voxels found!");
856 }
else if (this->sweepingVoxelCount() == 0) {
857 OPENVDB_THROW(RuntimeError,
"FastSweeping: No computing voxels found!");
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) {}
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) {
1135 OPENVDB_THROW(RuntimeError,
"FastSweeping: A SDF should not have active tiles!");
1138 #ifdef BENCHMARK_FAST_SWEEPING
1141 mParent->mSweepMask.clear();
1142 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
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");
1242 OPENVDB_THROW(RuntimeError,
"FastSweeping::InitExt expected an extension grid!");
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) {
1252 OPENVDB_THROW(RuntimeError,
"FastSweeping: A SDF should not have active tiles!");
1255 #ifdef BENCHMARK_FAST_SWEEPING
1256 util::CpuTimer timer(
"Initialize voxels");
1259 mParent->mSweepMask.clear();
1260 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1264 OpPoolT opPool(opPrototype);
1267 tree::LeafManager<SdfTreeT> mgr(tree1, 1);
1269 mgr.swapLeafBuffer(1);
1272 #ifdef BENCHMARK_FAST_SWEEPING
1273 timer.restart(
"Initialize tiles");
1276 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
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());
1314 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1315 math::GradStencil<SdfGridT, false>
stencil(*mSdfGrid);
1316 const math::Transform& xform = mExtGrid->transform();
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();
1331 const auto mask =
stencil.intersectionMask( isoValue );
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++)) {
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);
1375 sdf[voxelIter.pos()] = isAbove ? h /
math::Sqrt(sum1) : -h / math::
Sqrt(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>
1404 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1406 using LeafRange =
typename tree::LeafManager<const MaskTreeT>::LeafRange;
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
1415 util::CpuTimer timer;
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);
1432 tree::LeafManager<const MaskTreeT> mgr(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());
1441 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
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]; }
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 {
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);
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];
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];
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>
1842 if (fs.initSdf(fogGrid, isoValue,
false)) fs.sweep(nIter);
1843 return fs.sdfGrid();
1846 template<
typename Gr
idT>
1853 if (fs.initSdf(sdfGrid, isoValue,
true)) fs.sweep(nIter);
1854 return fs.sdfGrid();
1857 template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1858 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1861 const ExtValueT& background,
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);
1870 return fs.extGrid();
1873 template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1874 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1877 const ExtValueT &background,
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);
1886 return fs.extGrid();
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,
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);
1902 return std::make_pair(fs.sdfGrid(), fs.extGrid());
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,
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);
1918 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1921 template<
typename Gr
idT>
1930 if (fs.initDilate(sdfGrid, dilation, nn, mode)) fs.sweep(nIter);
1931 return fs.sdfGrid();
1934 template<
typename Gr
idT,
typename MaskTreeT>
1938 bool ignoreActiveTiles,
1942 if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1943 return fs.sdfGrid();
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
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Type Pow2(Type x)
Return x2.
void swap(UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &a, UT::ArraySet< Key, MULTI, MAX_LOAD_FACTOR_256, Clearer, Hash, KeyEqual > &b)
GLsizei const GLfloat * value
vfloat4 sqrt(const vfloat4 &a)
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
GLboolean GLboolean GLboolean GLboolean a
#define OPENVDB_USE_VERSION_NAMESPACE
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
GLfloat GLfloat GLfloat v2
NodeManager produces linear arrays of all tree nodes allowing for efficient threading and bottom-up p...
GLfloat GLfloat GLfloat GLfloat v3
__hostdev__ float getValue(uint32_t i) const
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays, one for each level of the tree.
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Coord Abs(const Coord &xyz)
ImageBuf OIIO_API dilate(const ImageBuf &src, int width=3, int height=-1, ROI roi={}, int nthreads=0)
Simple timer for basic profiling.
double restart()
Re-start timer.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
float Sqrt(float x)
Return the square root of a floating-point value.
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
IMATH_HOSTDEVICE constexpr int sign(T a) IMATH_NOEXCEPT
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Templated class to compute the minimum and maximum values.
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node's buffer with the nth corresponding auxiliary buffer, where n = bufferIdx...
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Container class that associates a tree with a transform and metadata.
GLfloat GLfloat GLfloat GLfloat h
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
double stop() const
Returns and prints time in milliseconds since construction or start was called.
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
GLubyte GLubyte GLubyte GLubyte w
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
GLint GLfloat GLint stencil
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
void setValueOff(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as inactive.
that also have some descendant prim *whose name begins with which in turn has a child named baz where *the predicate and *a name There is also one special expression reference
__hostdev__ float Sqrt(float x)
Return the square root of a floating-point value.
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...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
#define OPENVDB_THROW(exception, message)
**Note that the tasks the is the thread number *for the pool