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>
462 static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
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) {}
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 1298 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1299 void sumHelper(ExtT& sum2, ExtT ext,
bool update,
const SdfT& )
const {
if (update) sum2 = ext; }
1302 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
1303 void sumHelper(ExtT& sum2, ExtT ext,
bool ,
const SdfT& d2)
const { sum2 +=
static_cast<ExtValueT
>(d2 * ext); }
1305 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1306 ExtT extValHelper(ExtT extSum,
const SdfT& )
const {
return extSum; }
1308 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
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); }
1605 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
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; }
1609 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
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)); }
1613 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value,
int>::type = 0>
1614 ExtT
threeNghbr(
const NN& d1,
const NN& d2,
const NN& d3,
const SdfT& ,
const ExtT& v1,
const ExtT& v2,
const ExtT& v3)
const {
1621 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value,
int>::type = 0>
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 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
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
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
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
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
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Definition: Stencils.h:1231
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
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.
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
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
Definition: LeafManager.h:101
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:411
Coord Abs(const Coord &xyz)
Definition: Coord.h:517
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:908
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
Simple timer for basic profiling.
Definition: CpuTimer.h:66
SharedPtr< Grid > Ptr
Definition: Grid.h:573
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:14
Index32 Index
Definition: Types.h:54
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.
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
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
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
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
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
double restart()
Re-start timer.
Definition: CpuTimer.h:150
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:63
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212