98 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
99 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
101 #include <tbb/parallel_reduce.h>
102 #include <tbb/blocked_range.h>
103 #include <openvdb/Types.h>
104 #include <openvdb/Grid.h>
105 #include <openvdb/math/Math.h>
106 #include <openvdb/math/Transform.h>
107 #include <openvdb/util/NullInterrupter.h>
113 #include <type_traits>
121 namespace p2ls_internal {
125 template<
typename VisibleT,
typename BlindT>
class BlindData;
129 template<
typename SdfGridT,
130 typename AttributeT = void,
135 using DisableT =
typename std::is_void<AttributeT>::type;
141 using AttType =
typename std::conditional<DisableT::value, size_t, AttributeT>::type;
142 using AttGridType =
typename SdfGridT::template ValueConverter<AttType>::Type;
144 static_assert(std::is_floating_point<SdfType>::value,
145 "ParticlesToLevelSet requires an SDF grid with floating-point values");
181 void finalize(
bool prune =
false);
223 template <
typename ParticleListT>
224 void rasterizeSpheres(
const ParticleListT& pa);
231 template <
typename ParticleListT>
232 void rasterizeSpheres(
const ParticleListT& pa,
Real radius);
250 template <
typename ParticleListT>
251 void rasterizeTrails(
const ParticleListT& pa,
Real delta=1.0);
255 using BlindGridType =
typename SdfGridT::template ValueConverter<BlindType>::Type;
258 template<
typename ParticleListT,
typename Gr
idT>
struct Raster;
260 SdfGridType* mSdfGrid;
261 typename AttGridType::Ptr mAttGrid;
262 BlindGridType* mBlindGrid;
263 InterrupterT* mInterrupter;
264 Real mDx, mHalfWidth;
266 size_t mMinCount, mMaxCount;
271 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
272 inline ParticlesToLevelSet<SdfGridT, AttributeT, InterrupterT>::
273 ParticlesToLevelSet(SdfGridT& grid, InterrupterT* interrupter) :
276 mInterrupter(interrupter),
277 mDx(grid.voxelSize()[0]),
278 mHalfWidth(grid.background()/mDx),
285 if (!mSdfGrid->hasUniformVoxels() ) {
287 "ParticlesToLevelSet only supports uniform voxels!");
291 "ParticlesToLevelSet only supports level sets!"
292 "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
295 if (!DisableT::value) {
296 mBlindGrid =
new BlindGridType(
BlindType(grid.background()));
297 mBlindGrid->setTransform(mSdfGrid->transform().copy());
301 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
302 template <
typename ParticleListT>
306 if (DisableT::value) {
307 Raster<ParticleListT, SdfGridT> r(*
this, mSdfGrid, pa);
308 r.rasterizeSpheres();
310 Raster<ParticleListT, BlindGridType> r(*
this, mBlindGrid, pa);
311 r.rasterizeSpheres();
315 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
316 template <
typename ParticleListT>
320 if (DisableT::value) {
321 Raster<ParticleListT, SdfGridT> r(*
this, mSdfGrid, pa);
322 r.rasterizeSpheres(radius/mDx);
324 Raster<ParticleListT, BlindGridType> r(*
this, mBlindGrid, pa);
325 r.rasterizeSpheres(radius/mDx);
329 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
330 template <
typename ParticleListT>
334 if (DisableT::value) {
335 Raster<ParticleListT, SdfGridT> r(*
this, mSdfGrid, pa);
336 r.rasterizeTrails(delta);
338 Raster<ParticleListT, BlindGridType> r(*
this, mBlindGrid, pa);
339 r.rasterizeTrails(delta);
343 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
347 if (mBlindGrid ==
nullptr) {
354 using SdfTreeT =
typename SdfGridType::TreeType;
355 using AttTreeT =
typename AttGridType::TreeType;
356 using BlindTreeT =
typename BlindGridType::TreeType;
358 const BlindTreeT& tree = mBlindGrid->tree();
361 typename SdfTreeT::Ptr sdfTree(
new SdfTreeT(
365 typename AttTreeT::Ptr attTree(
new AttTreeT(
367 mAttGrid =
typename AttGridType::Ptr(
new AttGridType(attTree));
368 mAttGrid->setTransform(mBlindGrid->transform().copy());
373 using LeafIterT =
typename BlindTreeT::LeafCIter;
374 using LeafT =
typename BlindTreeT::LeafNodeType;
375 using SdfLeafT =
typename SdfTreeT::LeafNodeType;
376 using AttLeafT =
typename AttTreeT::LeafNodeType;
377 for (LeafIterT n = tree.cbeginLeaf(); n; ++n) {
378 const LeafT& leaf = *n;
381 SdfLeafT* sdfLeaf = sdfTree->probeLeaf(xyz);
382 AttLeafT* attLeaf = attTree->probeLeaf(xyz);
384 typename LeafT::ValueOnCIter m=leaf.cbeginValueOn();
388 sdfLeaf->setValueOnly(k, v.
visible());
389 attLeaf->setValueOnly(k, v.
blind());
395 sdfLeaf->setValueOnly(k, v.
visible());
396 attLeaf->setValueOnly(k, v.
blind());
403 if (mSdfGrid->empty()) {
404 mSdfGrid->setTree(sdfTree);
412 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
413 template<
typename ParticleListT,
typename Gr
idT>
416 using DisableT =
typename std::is_void<AttributeT>::type;
418 using SdfT =
typename ParticlesToLevelSetT::SdfType;
419 using AttT =
typename ParticlesToLevelSetT::AttType;
420 using ValueT =
typename GridT::ValueType;
421 using AccessorT =
typename GridT::Accessor;
422 using TreeT =
typename GridT::TreeType;
423 using LeafNodeT =
typename TreeT::LeafNodeType;
428 Raster(ParticlesToLevelSetT& parent, GridT* grid,
const ParticleListT& particles)
430 , mParticles(particles)
432 , mMap(*(mGrid->transform().baseMap()))
437 mPointPartitioner =
new PointPartitionerT();
438 mPointPartitioner->construct(particles, mGrid->transform());
442 Raster(Raster& other, tbb::split)
443 : mParent(other.mParent)
444 , mParticles(other.mParticles)
451 , mPointPartitioner(other.mPointPartitioner)
463 delete mPointPartitioner;
471 mMinCount = mMaxCount = 0;
472 if (mParent.mInterrupter) {
473 mParent.mInterrupter->start(
"Rasterizing particles to level set using spheres");
475 mTask = std::bind(&Raster::rasterSpheres, std::placeholders::_1, std::placeholders::_2);
477 if (mParent.mInterrupter) mParent.mInterrupter->end();
484 mMinCount = radius < mParent.mRmin ? mParticles.size() : 0;
485 mMaxCount = radius > mParent.mRmax ? mParticles.size() : 0;
486 if (mMinCount>0 || mMaxCount>0) {
487 mParent.mMinCount = mMinCount;
488 mParent.mMaxCount = mMaxCount;
490 if (mParent.mInterrupter) {
491 mParent.mInterrupter->start(
492 "Rasterizing particles to level set using const spheres");
494 mTask = std::bind(&Raster::rasterFixedSpheres,
495 std::placeholders::_1, std::placeholders::_2, SdfT(radius));
497 if (mParent.mInterrupter) mParent.mInterrupter->end();
516 mMinCount = mMaxCount = 0;
517 if (mParent.mInterrupter) {
518 mParent.mInterrupter->start(
"Rasterizing particles to level set using trails");
520 mTask = std::bind(&Raster::rasterTrails,
521 std::placeholders::_1, std::placeholders::_2, SdfT(delta));
523 if (mParent.mInterrupter) mParent.mInterrupter->end();
527 void operator()(
const tbb::blocked_range<size_t>& r)
531 mParent.mMinCount = mMinCount;
532 mParent.mMaxCount = mMaxCount;
536 void join(Raster& other)
539 mMinCount += other.mMinCount;
540 mMaxCount += other.mMaxCount;
544 Raster& operator=(
const Raster&) {
return *
this; }
547 bool ignoreParticle(SdfT R)
549 if (R < mParent.mRmin) {
553 if (R > mParent.mRmax) {
563 void rasterSpheres(
const tbb::blocked_range<size_t>& r)
565 AccessorT acc = mGrid->getAccessor();
567 const SdfT invDx = SdfT(1/mParent.mDx);
573 for (
size_t n = r.begin(), N = r.end(); n < N; ++n) {
575 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
576 for ( ; run && iter; ++iter) {
578 mParticles.getPosRad(
id, pos, rad);
579 const SdfT R = SdfT(invDx * rad);
580 if (this->ignoreParticle(R))
continue;
581 const Vec3R P = mMap.applyInverseMap(pos);
582 this->getAtt<DisableT>(
id, att);
583 run = this->makeSphere(P, R, att, acc);
592 void rasterFixedSpheres(
const tbb::blocked_range<size_t>& r, SdfT R)
595 dx = static_cast<SdfT>(mParent.mDx),
596 w = static_cast<SdfT>(mParent.mHalfWidth);
597 AccessorT acc = mGrid->getAccessor();
598 const ValueT inside = -mGrid->background();
599 const SdfT
max = R + w;
608 for (
size_t n = r.begin(), N = r.end(); n < N; ++n) {
610 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
611 for ( ; iter; ++iter) {
613 this->getAtt<DisableT>(
id, att);
614 mParticles.getPos(
id, pos);
615 const Vec3R P = mMap.applyInverseMap(pos);
618 for (Coord c = a; c.x() <= b.x(); ++c.x()) {
621 tbb::task::self().cancel_group_execution();
624 SdfT x2 = static_cast<SdfT>(
math::Pow2(c.x() - P[0]));
625 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
626 SdfT x2y2 = static_cast<SdfT>(x2 +
math::Pow2(c.y() - P[1]));
627 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
628 SdfT x2y2z2 = static_cast<SdfT>(
630 if (x2y2z2 >= max2 || (!acc.probeValue(c,v) && v<ValueT(0)))
632 if (x2y2z2 <= min2) {
633 acc.setValueOff(c, inside);
637 const ValueT d=Merge(dx*(
math::Sqrt(x2y2z2) - R), att);
638 if (d < v) acc.setValue(c, d);
650 void rasterTrails(
const tbb::blocked_range<size_t>& r, SdfT delta)
652 AccessorT acc = mGrid->getAccessor();
657 const Vec3R origin = mMap.applyInverseMap(
Vec3R(0,0,0));
658 const SdfT Rmin = SdfT(mParent.mRmin), invDx = SdfT(1/mParent.mDx);
661 for (
size_t n = r.begin(), N = r.end(); n < N; ++n) {
663 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
664 for ( ; run && iter; ++iter) {
666 mParticles.getPosRadVel(
id, pos, rad, vel);
667 const SdfT R0 = SdfT(invDx*rad);
668 if (this->ignoreParticle(R0))
continue;
669 this->getAtt<DisableT>(
id, att);
670 const Vec3R P0 = mMap.applyInverseMap(pos);
671 const Vec3R V = mMap.applyInverseMap(vel) - origin;
672 const SdfT speed = SdfT(V.length()), inv_speed = SdfT(1.0/speed);
673 const Vec3R Nrml = -V*inv_speed;
676 for (
size_t m=0; run && d <= speed ; ++m) {
677 run = this->makeSphere(P, R, att, acc);
678 P += 0.5*delta*R*Nrml;
679 d = SdfT((P-P0).length());
680 R = R0-(R0-Rmin)*d*inv_speed;
691 if (mParent.mGrainSize>0) {
692 tbb::parallel_reduce(
693 tbb::blocked_range<size_t>(0, bucketCount, mParent.mGrainSize), *
this);
695 (*this)(tbb::blocked_range<size_t>(0, bucketCount));
714 bool makeSphere(
const Vec3R &P, SdfT R,
const AttT& att, AccessorT& acc)
716 const ValueT inside = -mGrid->background();
717 const SdfT dx = SdfT(mParent.mDx), w = SdfT(mParent.mHalfWidth);
718 const SdfT
max = R + w;
725 for ( Coord c = a; c.x() <= b.x(); ++c.x() ) {
728 tbb::task::self().cancel_group_execution();
732 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
733 SdfT x2y2 = SdfT(x2 +
math::Pow2(c.y() - P[1]));
734 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
735 SdfT x2y2z2 = SdfT(x2y2 +
math::Pow2(c.z()-P[2]));
736 if (x2y2z2 >= max2 || (!acc.probeValue(c,v) && v<ValueT(0)))
738 if (x2y2z2 <= min2) {
739 acc.setValueOff(c, inside);
744 const ValueT d=Merge(dx*(
math::Sqrt(x2y2z2) - R), att);
745 if (d < v) acc.setValue(c, d);
751 using FuncType =
typename std::function<void (Raster*,
const tbb::blocked_range<size_t>&)>;
753 template<
typename DisableType>
754 typename std::enable_if<DisableType::value>::type
755 getAtt(
size_t, AttT&)
const {}
757 template<
typename DisableType>
758 typename std::enable_if<!DisableType::value>::type
759 getAtt(
size_t n, AttT& a)
const { mParticles.getAtt(n, a); }
762 typename std::enable_if<std::is_same<T, ValueT>::value, ValueT>::type
763 Merge(T s,
const AttT&)
const {
return s; }
766 typename std::enable_if<!std::is_same<T, ValueT>::value, ValueT>::type
767 Merge(T s,
const AttT& a)
const {
return ValueT(s,a); }
769 ParticlesToLevelSetT& mParent;
770 const ParticleListT& mParticles;
772 const math::MapBase& mMap;
773 size_t mMinCount, mMaxCount;
776 PointPartitionerT* mPointPartitioner;
782 namespace p2ls_internal {
787 template<
typename VisibleT,
typename BlindT>
818 template<
typename VisibleT,
typename BlindT>
821 ostr << rhs.visible();
826 template<
typename VisibleT,
typename BlindT>
840 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED