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>
55#include <unordered_map>
103template<
typename Gr
idT>
105fogToSdf(
const GridT &fogGrid,
106 typename GridT::ValueType isoValue,
136template<
typename Gr
idT>
138sdfToSdf(
const GridT &sdfGrid,
139 typename GridT::ValueType isoValue = 0,
192template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
193typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
194fogToExt(
const FogGridT &fogGrid,
196 const ExtValueT& background,
197 typename FogGridT::ValueType isoValue,
199 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
200 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid =
nullptr);
250template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
251typename 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);
313template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
314std::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);
376template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
377std::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);
414template<
typename Gr
idT>
440template<
typename Gr
idT,
typename MaskTreeT>
444 bool ignoreActiveTiles =
false,
459template<
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>
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();
692 struct PruneMinMaxFltKernel;
693 struct SweepingKernel;
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;
712template <
typename SdfGr
idT,
typename ExtValueT>
713const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
717template <
typename SdfGr
idT,
typename ExtValueT>
719 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(
FastSweepingDomain::
SWEEP_ALL), mIsInputSdf(true)
723template <
typename SdfGr
idT,
typename ExtValueT>
729 if (mExtGridInput) mExtGridInput.reset();
730 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
735template <
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;
764template <
typename SdfGr
idT,
typename ExtValueT>
768 mSdfGrid = fogGrid.deepCopy();
769 mIsInputSdf = isInputSdf;
771 kernel.
run(isoValue);
772 return this->isValid();
775template <
typename SdfGr
idT,
typename ExtValueT>
776template <
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;
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();
801template <
typename SdfGr
idT,
typename ExtValueT>
805 mSdfGrid = sdfGrid.deepCopy();
806 mSweepDirection = mode;
808 kernel.
run(dilate, nn);
809 return this->isValid();
812template <
typename SdfGr
idT,
typename ExtValueT>
813template<
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);
841 return this->isValid();
844template <
typename SdfGr
idT,
typename ExtValueT>
852 " a non-null reference extension grid input.");
854 if (this->boundaryVoxelCount() == 0) {
856 }
else if (this->sweepingVoxelCount() == 0) {
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);
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
916template <
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;
936 const bool vEqFltMin = v == -std::numeric_limits<SdfValueT>::max();
937 const bool vEqFltMax = v == std::numeric_limits<SdfValueT>::max();
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;
960template <
typename SdfGr
idT,
typename ExtValueT>
965 void operator()(
typename SdfTreeT::RootNodeType& node,
size_t = 1)
const {
966 for (
auto iter = node.beginValueAll(); iter; ++iter) {
967 if (*iter == -std::numeric_limits<SdfValueT>::max()) {
970 if (*iter == std::numeric_limits<SdfValueT>::max()) {
977 template<
typename NodeT>
980 for (
auto iter = node.beginValueAll(); iter; ++iter) {
981 if (*iter == -std::numeric_limits<SdfValueT>::max()) {
984 if (*iter == std::numeric_limits<SdfValueT>::max()) {
991 void operator()(
typename SdfTreeT::LeafNodeType& leaf,
size_t = 1)
const
993 for (
auto iter = leaf.beginValueOn(); iter; ++iter) {
994 if (*iter == -std::numeric_limits<SdfValueT>::max()) {
997 if (*iter == std::numeric_limits<SdfValueT>::max()) {
1008template <
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");
1031 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
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 ) {
1060 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
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
1118template <
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);
1156 tree.root().setBackground(std::numeric_limits<SdfValueT>::max(),
false);
1157 if (hasActiveTiles) tree.voxelizeActiveTiles();
1161 mParent->computeSweepMaskLeafOrigins();
1168 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
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();
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;) {
1192 SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1194 if (mask.test(i++)) {
1198 if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1208 template<
typename RootOrInternalNodeT>
1211 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1212 for (
auto it = node.cbeginValueAll(); it; ++it) {
1213 SdfValueT& v =
const_cast<SdfValueT&
>(*it);
1214 v = v > isoValue ? above : -above;
1227template <
typename SdfGr
idT,
typename ExtValueT>
1228template <
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
1256 util::CpuTimer timer(
"Initialize voxels");
1259 mParent->mSweepMask.clear();
1260 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1264 OpPoolT opPool(opPrototype);
1268 tbb::parallel_for(mgr.leafRange(32), *
this);
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);
1278 tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(),
false);
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());
1314 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1315 math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1316 const math::Transform& xform = mExtGrid->transform();
1317 typename OpPoolT::reference op = mOpPool->local();
1318 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
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);
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>();
1349 SdfValueT minD = std::numeric_limits<SdfValueT>::max();
1350 for (
int n=0, i=0; i<6;) {
1351 SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
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);
1375 sdf[voxelIter.pos()] = isAbove ? h /
math::Sqrt(sum1) : -h / math::
Sqrt(sum1);
1383 template<
typename RootOrInternalNodeT>
1384 void operator()(
const RootOrInternalNodeT& node)
const
1386 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
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;
1402template <
typename SdfGr
idT,
typename ExtValueT>
1403template <
typename MaskTreeT>
1404struct 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();
1419 static const SdfValueT
Unknown = std::numeric_limits<SdfValueT>::max();
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 ) {
1446 static const SdfValueT
Unknown = std::numeric_limits<SdfValueT>::max();
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
1474template <
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) {}
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());
1676 d1 = std::min(
NN(acc1, ijk, 0),
NN(acc1, ijk, 1));
1677 d2 = std::min(
NN(acc1, ijk, 2),
NN(acc1, ijk, 3));
1678 d3 = std::min(
NN(acc1, ijk, 4),
NN(acc1, ijk, 5));
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) {
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;
1835template<
typename Gr
idT>
1838 typename GridT::ValueType isoValue,
1842 if (fs.initSdf(fogGrid, isoValue,
false)) fs.
sweep(nIter);
1843 return fs.sdfGrid();
1846template<
typename Gr
idT>
1849 typename GridT::ValueType isoValue,
1853 if (fs.initSdf(sdfGrid, isoValue,
true)) fs.
sweep(nIter);
1854 return fs.sdfGrid();
1857template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1858typename 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);
1870 return fs.extGrid();
1873template<
typename SdfGr
idT,
typename OpT,
typename ExtValueT>
1874typename 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);
1886 return fs.extGrid();
1889template<
typename FogGr
idT,
typename ExtOpT,
typename ExtValueT>
1890std::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);
1902 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1905template<
typename SdfGr
idT,
typename ExtOpT,
typename ExtValueT>
1906std::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);
1918 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1921template<
typename Gr
idT>
1930 if (fs.initDilate(sdfGrid, dilation, nn, mode)) fs.
sweep(nIter);
1931 return fs.sdfGrid();
1934template<
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)
ValueT value
Definition: GridBuilder.h:1290
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean,...
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:411
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:571
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:908
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Definition: Exceptions.h:63
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:98
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
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:48
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Definition: Stencils.h:1233
Templated class to compute the minimum and maximum values.
Definition: Stats.h:31
Definition: LeafManager.h:102
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node's buffer with the nth corresponding auxiliary buffer, where n = bufferIdx.
Definition: LeafManager.h:359
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays,...
Definition: NodeManager.h:531
void foreachTopDown(const NodeOp &op, bool threaded=true, size_t grainSize=1)
Definition: NodeManager.h:631
Definition: ValueAccessor.h:191
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
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:235
Simple timer for basic profiling.
Definition: CpuTimer.h:67
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
double restart()
Re-start timer.
Definition: CpuTimer.h:150
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:13
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 isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance.
Definition: Math.h:349
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:931
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
Index32 Index
Definition: Types.h:54
@ GRID_LEVEL_SET
Definition: Types.h:416
math::Vec3< Real > Vec3R
Definition: Types.h:72
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
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
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157