62 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
63 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
65 #include <tbb/parallel_reduce.h>
66 #include <tbb/blocked_range.h>
80 #include <type_traits>
93 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
94 inline void particlesToSdf(
const ParticleListT&, GridT&, InterrupterT* =
nullptr);
100 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
101 inline void particlesToSdf(
const ParticleListT&, GridT&,
Real radius, InterrupterT* =
nullptr);
110 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
117 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
118 inline void particlesToMask(
const ParticleListT&, GridT&, InterrupterT* =
nullptr);
124 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
125 inline void particlesToMask(
const ParticleListT&, GridT&,
Real radius, InterrupterT* =
nullptr);
134 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
141 namespace p2ls_internal {
146 template<
typename VisibleT,
typename BlindT>
class BlindData;
150 template<
typename SdfGridT,
151 typename AttributeT =
void,
237 template<
typename ParticleListT>
246 template<
typename ParticleListT>
264 template<
typename ParticleListT>
268 using BlindType = p2ls_internal::BlindData<SdfType, AttType>;
272 template<
typename ParticleListT,
typename Gr
idT>
struct Raster;
275 typename AttGridType::Ptr mAttGrid;
276 BlindGridType* mBlindGrid;
277 InterrupterT* mInterrupter;
278 Real mDx, mHalfWidth;
280 size_t mMinCount, mMaxCount;
285 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
290 mInterrupter(interrupter),
291 mDx(grid.voxelSize()[0]),
292 mHalfWidth(grid.background()/mDx),
299 if (!mSdfGrid->hasUniformVoxels()) {
300 OPENVDB_THROW(RuntimeError,
"ParticlesToLevelSet only supports uniform voxels!");
303 mBlindGrid =
new BlindGridType(BlindType(grid.background()));
304 mBlindGrid->setTransform(mSdfGrid->transform().copy());
308 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
309 template<
typename ParticleListT>
314 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
315 r.rasterizeSpheres();
317 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
318 r.rasterizeSpheres();
322 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
323 template<
typename ParticleListT>
328 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
329 r.rasterizeSpheres(radius/mDx);
331 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
332 r.rasterizeSpheres(radius/mDx);
336 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
337 template<
typename ParticleListT>
342 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
343 r.rasterizeTrails(delta);
345 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
346 r.rasterizeTrails(delta);
351 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
370 using AttTreeT =
typename AttGridType::TreeType;
371 using AttLeafT =
typename AttTreeT::LeafNodeType;
372 using BlindTreeT =
typename BlindGridType::TreeType;
373 using BlindLeafIterT =
typename BlindTreeT::LeafCIter;
374 using BlindLeafT =
typename BlindTreeT::LeafNodeType;
375 using SdfTreeT =
typename SdfGridType::TreeType;
376 using SdfLeafT =
typename SdfTreeT::LeafNodeType;
379 const BlindTreeT& blindTree = mBlindGrid->tree();
382 typename AttTreeT::Ptr attTree(
new AttTreeT(
383 blindTree, blindTree.background().blind(), openvdb::TopologyCopy()));
385 mAttGrid =
typename AttGridType::Ptr(
new AttGridType(attTree));
386 mAttGrid->setTransform(mBlindGrid->transform().copy());
388 typename SdfTreeT::Ptr sdfTree;
392 sdfTree.reset(
new SdfTreeT(blindTree,
397 leafNodes.
foreach([&](AttLeafT& attLeaf,
size_t ) {
398 if (
const auto* blindLeaf = blindTree.probeConstLeaf(attLeaf.origin())) {
399 for (
auto iter = attLeaf.beginValueOn(); iter; ++iter) {
400 const auto pos = iter.pos();
401 attLeaf.setValueOnly(pos, blindLeaf->getValue(pos).blind());
406 const auto blindAcc = mBlindGrid->getConstAccessor();
407 auto iter = attTree->beginValueOn();
408 iter.setMaxDepth(AttTreeT::ValueOnIter::LEAF_DEPTH - 1);
409 for ( ; iter; ++iter) {
410 iter.modifyValue([&](
AttType&
v) { v = blindAcc.getValue(iter.getCoord()).blind(); });
415 sdfTree.reset(
new SdfTreeT(blindTree, blindTree.background().visible(),
TopologyCopy()));
416 for (BlindLeafIterT
n = blindTree.cbeginLeaf();
n; ++
n) {
417 const BlindLeafT& leaf = *
n;
418 const openvdb::Coord xyz = leaf.origin();
420 SdfLeafT* sdfLeaf = sdfTree->probeLeaf(xyz);
421 AttLeafT* attLeaf = attTree->probeLeaf(xyz);
423 typename BlindLeafT::ValueOnCIter
m=leaf.cbeginValueOn();
426 const BlindType&
v = leaf.getValue(k);
427 sdfLeaf->setValueOnly(k, v.visible());
428 attLeaf->setValueOnly(k, v.blind());
433 const BlindType&
v = *
m;
434 sdfLeaf->setValueOnly(k, v.visible());
435 attLeaf->setValueOnly(k, v.blind());
442 if (mSdfGrid->empty()) {
443 mSdfGrid->setTree(sdfTree);
446 mSdfGrid->tree().topologyUnion(*sdfTree);
460 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
461 template<
typename ParticleListT,
typename Gr
idT>
466 using SdfT =
typename ParticlesToLevelSetT::SdfType;
467 using AttT =
typename ParticlesToLevelSetT::AttType;
468 using ValueT =
typename GridT::ValueType;
469 using AccessorT =
typename GridT::Accessor;
470 using TreeT =
typename GridT::TreeType;
471 using LeafNodeT =
typename TreeT::LeafNodeType;
479 Raster(ParticlesToLevelSetT& parent, GridT* grid,
const ParticleListT& particles)
481 , mParticles(particles)
488 mPointPartitioner =
new PointPartitionerT;
489 mPointPartitioner->construct(particles, mGrid->transform());
494 : mParent(other.mParent)
495 , mParticles(other.mParticles)
496 , mGrid(new GridT(*other.mGrid, openvdb::
ShallowCopy()))
502 , mPointPartitioner(other.mPointPartitioner)
514 delete mPointPartitioner;
520 mMinCount = mMaxCount = 0;
521 if (mParent.mInterrupter) {
522 mParent.mInterrupter->start(
"Rasterizing particles to level set using spheres");
524 mTask = std::bind(&Raster::rasterSpheres, std::placeholders::_1, std::placeholders::_2);
526 if (mParent.mInterrupter) mParent.mInterrupter->end();
531 mMinCount = radius < mParent.mRmin ? mParticles.size() : 0;
532 mMaxCount = radius > mParent.mRmax ? mParticles.size() : 0;
533 if (mMinCount>0 || mMaxCount>0) {
534 mParent.mMinCount = mMinCount;
535 mParent.mMaxCount = mMaxCount;
537 if (mParent.mInterrupter) {
538 mParent.mInterrupter->start(
539 "Rasterizing particles to level set using const spheres");
541 mTask = std::bind(&Raster::rasterFixedSpheres,
542 std::placeholders::_1, std::placeholders::_2, radius);
544 if (mParent.mInterrupter) mParent.mInterrupter->end();
550 mMinCount = mMaxCount = 0;
551 if (mParent.mInterrupter) {
552 mParent.mInterrupter->start(
"Rasterizing particles to level set using trails");
554 mTask = std::bind(&Raster::rasterTrails,
555 std::placeholders::_1, std::placeholders::_2, delta);
557 if (mParent.mInterrupter) mParent.mInterrupter->end();
561 void operator()(
const tbb::blocked_range<size_t>&
r)
565 mParent.mMinCount = mMinCount;
566 mParent.mMaxCount = mMaxCount;
570 void join(Raster& other)
577 mGrid->topologyUnion(*other.mGrid);
583 mMinCount += other.mMinCount;
584 mMaxCount += other.mMaxCount;
589 Raster& operator=(
const Raster&) {
return *
this; }
592 bool ignoreParticle(
Real R)
594 if (R < mParent.mRmin) {
598 if (R > mParent.mRmax) {
607 void rasterSpheres(
const tbb::blocked_range<size_t>&
r)
609 AccessorT acc = mGrid->getAccessor();
611 const Real invDx = 1 / mParent.mDx;
617 for (
size_t n = r.begin(),
N = r.
end();
n <
N; ++
n) {
619 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(
n);
620 for ( ; run && iter; ++iter) {
622 mParticles.getPosRad(
id, pos, rad);
623 const Real R = invDx * rad;
624 if (this->ignoreParticle(R))
continue;
625 const Vec3R P = mMap.applyInverseMap(pos);
626 this->getAtt<DisableT>(
id, att);
627 run = this->makeSphere(P, R, att, acc);
635 void rasterFixedSpheres(
const tbb::blocked_range<size_t>& r,
Real R)
637 AccessorT acc = mGrid->getAccessor();
642 for (
size_t n = r.begin(), N = r.end();
n <
N; ++
n) {
644 for (
auto iter = mPointPartitioner->indices(
n); iter; ++iter) {
646 this->getAtt<DisableT>(
id, att);
647 mParticles.getPos(
id, pos);
648 const Vec3R P = mMap.applyInverseMap(pos);
649 this->makeSphere(P, R, att, acc);
657 void rasterTrails(
const tbb::blocked_range<size_t>& r,
Real delta)
659 AccessorT acc = mGrid->getAccessor();
664 const Vec3R origin = mMap.applyInverseMap(
Vec3R(0,0,0));
665 const Real Rmin = mParent.mRmin, invDx = 1 / mParent.mDx;
668 for (
size_t n = r.begin(), N = r.end();
n <
N; ++
n) {
670 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(
n);
671 for ( ; run && iter; ++iter) {
673 mParticles.getPosRadVel(
id, pos, rad, vel);
674 const Real R0 = invDx * rad;
675 if (this->ignoreParticle(R0))
continue;
676 this->getAtt<DisableT>(
id, att);
677 const Vec3R P0 = mMap.applyInverseMap(pos);
678 const Vec3R V = mMap.applyInverseMap(vel) - origin;
679 const Real speed = V.length(), invSpeed = 1.0 / speed;
680 const Vec3R Nrml = -V * invSpeed;
683 for (
size_t m = 0; run && d <= speed ; ++
m) {
684 run = this->makeSphere(P, R, att, acc);
685 P += 0.5 * delta * R * Nrml;
687 R = R0 - (R0 - Rmin) * d * invSpeed;
698 if (mParent.mGrainSize>0) {
699 tbb::parallel_reduce(
700 tbb::blocked_range<size_t>(0, bucketCount, mParent.mGrainSize), *
this);
702 (*this)(tbb::blocked_range<size_t>(0, bucketCount));
720 template <
bool IsMaskT = OutputIsMask>
722 makeSphere(
const Vec3R& P,
Real R,
const AttT& att, AccessorT& acc)
726 w = mParent.mHalfWidth,
734 const ValueT inside = -mGrid->background();
738 for (Coord
c = lo;
c.x() <= hi.x(); ++
c.x()) {
741 tbb::task::self().cancel_group_execution();
745 for (
c.y() = lo.y();
c.y() <= hi.y(); ++
c.y()) {
747 for (
c.z() = lo.z();
c.z() <= hi.z(); ++
c.z()) {
749 #if defined __INTEL_COMPILER
750 _Pragma(
"warning (push)")
751 _Pragma("warning (disable:186)")
753 if (x2y2z2 >= max2 || (!acc.probeValue(
c, v) && (v < ValueT(0))))
755 #if defined __INTEL_COMPILER
756 _Pragma(
"warning (pop)")
758 if (x2y2z2 <= min2) {
759 acc.setValueOff(
c, inside);
764 const ValueT d = Merge(static_cast<SdfT>(dx*(
math::Sqrt(x2y2z2)-R)), att);
765 if (d < v) acc.setValue(
c, d);
774 template <
bool IsMaskT = OutputIsMask>
776 makeSphere(
const Vec3R&
p,
Real r,
const AttT& att, AccessorT& acc)
790 const std::vector<CoordBBox> padding{
791 CoordBBox(outLo.x(), outLo.y(), outLo.z(), inLo.x()-1, outHi.y(), outHi.z()),
792 CoordBBox(inHi.x()+1, outLo.y(), outLo.z(), outHi.x(), outHi.y(), outHi.z()),
793 CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), inLo.y()-1, outHi.z()),
794 CoordBBox(outLo.x(), inHi.y()+1, outLo.z(), outHi.x(), outHi.y(), outHi.z()),
795 CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), outHi.y(), inLo.z()-1),
796 CoordBBox(outLo.x(), outLo.y(), inHi.z()+1, outHi.x(), outHi.y(), outHi.z()),
798 const ValueT onValue = Merge(SdfT(1), att);
802 acc.tree().sparseFill(CoordBBox(inLo, inHi), onValue);
805 for (
const auto& bbox: padding) {
807 tbb::task::self().cancel_group_execution();
810 const Coord &bmin = bbox.min(), &bmax = bbox.max();
813 for (c = bmin, cx = c.x(); c.x() <= bmax.x(); ++c.x(), cx += 1) {
815 for (c.y() = bmin.y(), cy = c.y(); c.y() <= bmax.y(); ++c.y(), cy += 1) {
817 for (c.z() = bmin.z(), cz = c.z(); c.z() <= bmax.z(); ++c.z(), cz += 1) {
819 if (x2y2z2 < rSquared) {
820 acc.setValue(c, onValue);
829 using FuncType =
typename std::function<void (Raster*, const tbb::blocked_range<size_t>&)>;
831 template<
typename DisableType>
833 getAtt(
size_t, AttT&)
const {}
835 template<
typename DisableType>
837 getAtt(
size_t n, AttT&
a)
const { mParticles.getAtt(n, a); }
841 Merge(T
s,
const AttT&)
const {
return s; }
845 Merge(T
s,
const AttT& a)
const {
return ValueT(s,a); }
847 ParticlesToLevelSetT& mParent;
848 const ParticleListT& mParticles;
850 const math::MapBase& mMap;
851 size_t mMinCount, mMaxCount;
854 PointPartitionerT* mPointPartitioner;
861 namespace p2ls_internal {
867 template<
typename VisibleT,
typename BlindT>
871 using type = VisibleT;
872 using VisibleType = VisibleT;
873 using BlindType = BlindT;
876 explicit BlindData(VisibleT v) : mVisible(v), mBlind(
zeroVal<BlindType>()) {}
877 BlindData(VisibleT v, BlindT
b) : mVisible(v), mBlind(b) {}
878 BlindData(
const BlindData&) =
default;
879 BlindData& operator=(
const BlindData&) =
default;
880 const VisibleT& visible()
const {
return mVisible; }
881 const BlindT& blind()
const {
return mBlind; }
883 bool operator==(
const BlindData& rhs)
const {
return mVisible == rhs.mVisible; }
885 bool operator< (
const BlindData& rhs)
const {
return mVisible < rhs.mVisible; }
886 bool operator> (
const BlindData& rhs)
const {
return mVisible > rhs.mVisible; }
887 BlindData
operator+(
const BlindData& rhs)
const {
return BlindData(mVisible + rhs.mVisible); }
888 BlindData
operator-(
const BlindData& rhs)
const {
return BlindData(mVisible - rhs.mVisible); }
889 BlindData
operator-()
const {
return BlindData(-mVisible, mBlind); }
898 template<
typename VisibleT,
typename BlindT>
899 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
901 ostr << rhs.visible();
907 template<
typename VisibleT,
typename BlindT>
908 inline BlindData<VisibleT, BlindT>
Abs(
const BlindData<VisibleT, BlindT>&
x)
910 return BlindData<VisibleT, BlindT>(
math::Abs(x.visible()), x.blind());
915 template<
typename VisibleT,
typename BlindT,
typename T>
916 inline BlindData<VisibleT, BlindT>
917 operator+(
const BlindData<VisibleT, BlindT>& x,
const T& rhs)
919 return BlindData<VisibleT, BlindT>(x.visible() +
static_cast<VisibleT
>(rhs), x.blind());
930 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
935 "particlesToSdf requires an SDF grid with floating-point values");
939 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
947 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
952 "particlesToSdf requires an SDF grid with floating-point values");
956 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
964 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
969 "particleTrailsToSdf requires an SDF grid with floating-point values");
973 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
981 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
986 "particlesToMask requires a boolean-valued grid");
992 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
997 "particlesToMask requires a boolean-valued grid");
1003 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
1008 "particleTrailsToMask requires a boolean-valued grid");
1018 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
size_t getMaxCount() const
Return the number of particles that were ignored because they were larger than the maximum radius...
vint4 max(const vint4 &a, const vint4 &b)
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
int Ceil(float x)
Return the ceiling of x.
#define OPENVDB_LOG_WARN(mesg)
Functions to efficiently perform various compositing operations on grids.
Tag dispatch class that distinguishes shallow copy constructors from deep copy constructors.
Type Pow2(Type x)
Return x2.
Real getRmax() const
Return the largest radius allowed in voxel units.
AttGridType::Ptr attributeGrid()
Return a pointer to the grid containing the optional user-defined attribute.
GLuint GLenum GLenum transform
Real getVoxelSize() const
Return the size of a voxel in world units.
void finalize(bool prune=false)
This method syncs up the level set and attribute grids and therefore needs to be called before any of...
GLboolean GLboolean GLboolean GLboolean a
Plane3< T > operator-(const Plane3< T > &plane)
#define OPENVDB_USE_VERSION_NAMESPACE
GA_API const UT_StringHolder P
Dummy NOOP interrupter class defining interface.
void setGrainSize(int grainSize)
Set the grain size used for threading.
typename SdfGridT::ValueType SdfType
ParticlesToLevelSet(SdfGridT &grid, InterrupterT *interrupt=nullptr)
Constructor using an existing boolean or narrow-band level set grid.
Spatially partitions points using a parallel radix-based sorting algorithm.
SYS_FORCE_INLINE const_iterator end() const
bool operator==(const BaseDimensions< T > &a, const BaseDimensions< Y > &b)
int getGrainSize() const
Return the grain size used for threading.
OIIO_FORCEINLINE vbool4 operator>(const vint4 &a, const vint4 &b)
void OIIO_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
GLint GLint GLint GLint GLint x
Coord Abs(const Coord &xyz)
void rasterizeSpheres(const ParticleListT &pa)
Rasterize each particle as a sphere with the particle's position and radius.
typename std::conditional< DisableT::value, size_t, AttributeT >::type AttType
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
static const bool OutputIsMask
void setRmax(Real Rmax)
Set the largest radius allowed in voxel units.
GLubyte GLubyte GLubyte GLubyte w
float Sqrt(float x)
Return the square root of a floating-point value.
Defined various multi-threaded utility functions for trees.
GLuint GLsizei GLsizei * length
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
bool operator<(const GU_TetrahedronFacet &a, const GU_TetrahedronFacet &b)
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
typename SdfGridT::template ValueConverter< AttType >::Type AttGridType
GLuint GLuint GLsizei GLenum type
size_t getMinCount() const
Return the number of particles that were ignored because they were smaller than the minimum radius...
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
bool ignoredParticles() const
Return true if any particles were ignored due to their size.
GLdouble GLdouble GLdouble b
void setRmin(Real Rmin)
Set the smallest radius allowed in voxel units.
void rasterizeTrails(const ParticleListT &pa, Real delta=1.0)
Rasterize each particle as a trail comprising the CSG union of spheres of decreasing radius...
GLdouble GLdouble GLdouble r
GLuint GLuint GLsizei count
GA_API const UT_StringHolder N
InterrupterT InterrupterType
int Floor(float x)
Return the floor of x.
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...
T zeroVal()
Return the value of type T that corresponds to zero.
typename std::is_void< AttributeT >::type DisableT
bool wasInterrupted(T *i, int percent=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
GLsizei const GLfloat * value
Real getRmin() const
Return the smallest radius allowed in voxel units.
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
#define OPENVDB_THROW(exception, message)
Quat< T > operator+(const Quat< T > &q1, const Quat< T > &q2)
#define OPENVDB_NO_FP_EQUALITY_WARNING_END