62 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
63 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
79 #include <tbb/parallel_reduce.h>
80 #include <tbb/blocked_range.h>
84 #include <type_traits>
97 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
98 inline void particlesToSdf(
const ParticleListT&, GridT&, InterrupterT* =
nullptr);
104 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
105 inline void particlesToSdf(
const ParticleListT&, GridT&,
Real radius, InterrupterT* =
nullptr);
114 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
121 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
122 inline void particlesToMask(
const ParticleListT&, GridT&, InterrupterT* =
nullptr);
128 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
129 inline void particlesToMask(
const ParticleListT&, GridT&,
Real radius, InterrupterT* =
nullptr);
138 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT = util::NullInterrupter>
146 namespace p2ls_internal {
151 template<
typename VisibleT,
typename BlindT>
class BlindData;
156 template<
typename SdfGridT,
157 typename AttributeT =
void,
158 typename InterrupterT = util::NullInterrupter>
243 template<
typename ParticleListT>
252 template<
typename ParticleListT>
270 template<
typename ParticleListT>
274 using BlindType = p2ls_internal::BlindData<SdfType, AttType>;
278 template<
typename ParticleListT,
typename Gr
idT>
struct Raster;
281 typename AttGridType::Ptr mAttGrid;
282 BlindGridType* mBlindGrid;
283 InterrupterT* mInterrupter;
284 Real mDx, mHalfWidth;
286 size_t mMinCount, mMaxCount;
291 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
296 mInterrupter(interrupter),
297 mDx(grid.voxelSize()[0]),
298 mHalfWidth(grid.background()/mDx),
305 if (!mSdfGrid->hasUniformVoxels()) {
306 OPENVDB_THROW(RuntimeError,
"ParticlesToLevelSet only supports uniform voxels!");
309 mBlindGrid =
new BlindGridType(BlindType(grid.background()));
310 mBlindGrid->setTransform(mSdfGrid->transform().copy());
314 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
315 template<
typename ParticleListT>
320 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
321 r.rasterizeSpheres();
323 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
324 r.rasterizeSpheres();
328 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
329 template<
typename ParticleListT>
334 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
335 r.rasterizeSpheres(radius/mDx);
337 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
338 r.rasterizeSpheres(radius/mDx);
342 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
343 template<
typename ParticleListT>
348 Raster<ParticleListT, SdfGridT>
r(*
this, mSdfGrid, pa);
349 r.rasterizeTrails(delta);
351 Raster<ParticleListT, BlindGridType>
r(*
this, mBlindGrid, pa);
352 r.rasterizeTrails(delta);
357 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
376 using AttTreeT =
typename AttGridType::TreeType;
377 using AttLeafT =
typename AttTreeT::LeafNodeType;
378 using BlindTreeT =
typename BlindGridType::TreeType;
379 using BlindLeafIterT =
typename BlindTreeT::LeafCIter;
380 using BlindLeafT =
typename BlindTreeT::LeafNodeType;
381 using SdfTreeT =
typename SdfGridType::TreeType;
382 using SdfLeafT =
typename SdfTreeT::LeafNodeType;
385 const BlindTreeT& blindTree = mBlindGrid->tree();
388 typename AttTreeT::Ptr attTree(
new AttTreeT(
389 blindTree, blindTree.background().blind(), openvdb::TopologyCopy()));
391 mAttGrid =
typename AttGridType::Ptr(
new AttGridType(attTree));
392 mAttGrid->setTransform(mBlindGrid->transform().copy());
394 typename SdfTreeT::Ptr sdfTree;
398 sdfTree.reset(
new SdfTreeT(blindTree,
403 leafNodes.
foreach([&](AttLeafT& attLeaf,
size_t ) {
404 if (
const auto* blindLeaf = blindTree.probeConstLeaf(attLeaf.origin())) {
405 for (
auto iter = attLeaf.beginValueOn(); iter; ++iter) {
406 const auto pos = iter.pos();
407 attLeaf.setValueOnly(pos, blindLeaf->getValue(pos).blind());
412 const auto blindAcc = mBlindGrid->getConstAccessor();
413 auto iter = attTree->beginValueOn();
414 iter.setMaxDepth(AttTreeT::ValueOnIter::LEAF_DEPTH - 1);
415 for ( ; iter; ++iter) {
416 iter.modifyValue([&](
AttType&
v) { v = blindAcc.getValue(iter.getCoord()).blind(); });
421 sdfTree.reset(
new SdfTreeT(blindTree, blindTree.background().visible(),
TopologyCopy()));
422 for (BlindLeafIterT
n = blindTree.cbeginLeaf();
n; ++
n) {
423 const BlindLeafT& leaf = *
n;
424 const openvdb::Coord xyz = leaf.origin();
426 SdfLeafT* sdfLeaf = sdfTree->probeLeaf(xyz);
427 AttLeafT* attLeaf = attTree->probeLeaf(xyz);
429 typename BlindLeafT::ValueOnCIter m=leaf.cbeginValueOn();
432 const BlindType&
v = leaf.getValue(k);
433 sdfLeaf->setValueOnly(k, v.visible());
434 attLeaf->setValueOnly(k, v.blind());
438 const openvdb::Index k = m.pos();
439 const BlindType&
v = *m;
440 sdfLeaf->setValueOnly(k, v.visible());
441 attLeaf->setValueOnly(k, v.blind());
448 if (mSdfGrid->empty()) {
449 mSdfGrid->setTree(sdfTree);
452 mSdfGrid->tree().topologyUnion(*sdfTree);
466 template<
typename SdfGr
idT,
typename AttributeT,
typename InterrupterT>
467 template<
typename ParticleListT,
typename Gr
idT>
472 using SdfT =
typename ParticlesToLevelSetT::SdfType;
473 using AttT =
typename ParticlesToLevelSetT::AttType;
475 using AccessorT =
typename GridT::Accessor;
476 using TreeT =
typename GridT::TreeType;
477 using LeafNodeT =
typename TreeT::LeafNodeType;
485 Raster(ParticlesToLevelSetT& parent, GridT* grid,
const ParticleListT& particles)
487 , mParticles(particles)
494 mPointPartitioner =
new PointPartitionerT;
495 mPointPartitioner->construct(particles, mGrid->transform());
500 : mParent(other.mParent)
501 , mParticles(other.mParticles)
502 , mGrid(new GridT(*other.mGrid, openvdb::
ShallowCopy()))
508 , mPointPartitioner(other.mPointPartitioner)
520 delete mPointPartitioner;
526 mMinCount = mMaxCount = 0;
527 if (mParent.mInterrupter) {
528 mParent.mInterrupter->start(
"Rasterizing particles to level set using spheres");
530 mTask = std::bind(&Raster::rasterSpheres, std::placeholders::_1, std::placeholders::_2);
532 if (mParent.mInterrupter) mParent.mInterrupter->end();
537 mMinCount = radius < mParent.mRmin ? mParticles.size() : 0;
538 mMaxCount = radius > mParent.mRmax ? mParticles.size() : 0;
539 if (mMinCount>0 || mMaxCount>0) {
540 mParent.mMinCount = mMinCount;
541 mParent.mMaxCount = mMaxCount;
543 if (mParent.mInterrupter) {
544 mParent.mInterrupter->start(
545 "Rasterizing particles to level set using const spheres");
547 mTask = std::bind(&Raster::rasterFixedSpheres,
548 std::placeholders::_1, std::placeholders::_2, radius);
550 if (mParent.mInterrupter) mParent.mInterrupter->end();
556 mMinCount = mMaxCount = 0;
557 if (mParent.mInterrupter) {
558 mParent.mInterrupter->start(
"Rasterizing particles to level set using trails");
560 mTask = std::bind(&Raster::rasterTrails,
561 std::placeholders::_1, std::placeholders::_2, delta);
563 if (mParent.mInterrupter) mParent.mInterrupter->end();
567 void operator()(
const tbb::blocked_range<size_t>&
r)
571 mParent.mMinCount = mMinCount;
572 mParent.mMaxCount = mMaxCount;
576 void join(Raster& other)
583 mGrid->topologyUnion(*other.mGrid);
589 mMinCount += other.mMinCount;
590 mMaxCount += other.mMaxCount;
595 Raster&
operator=(
const Raster&) {
return *
this; }
598 bool ignoreParticle(
Real R)
600 if (R < mParent.mRmin) {
604 if (R > mParent.mRmax) {
613 void rasterSpheres(
const tbb::blocked_range<size_t>&
r)
615 AccessorT acc = mGrid->getAccessor();
617 const Real invDx = 1 / mParent.mDx;
623 for (
size_t n = r.begin(),
N = r.
end();
n <
N; ++
n) {
625 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(
n);
626 for ( ; run && iter; ++iter) {
628 mParticles.getPosRad(
id, pos, rad);
629 const Real R = invDx * rad;
630 if (this->ignoreParticle(R))
continue;
631 const Vec3R P = mMap.applyInverseMap(pos);
632 this->getAtt<DisableT>(
id, att);
633 run = this->makeSphere(P, R, att, acc);
641 void rasterFixedSpheres(
const tbb::blocked_range<size_t>& r,
Real R)
643 AccessorT acc = mGrid->getAccessor();
648 for (
size_t n = r.begin(), N = r.end();
n <
N; ++
n) {
650 for (
auto iter = mPointPartitioner->indices(
n); iter; ++iter) {
652 this->getAtt<DisableT>(
id, att);
653 mParticles.getPos(
id, pos);
654 const Vec3R P = mMap.applyInverseMap(pos);
655 this->makeSphere(P, R, att, acc);
663 void rasterTrails(
const tbb::blocked_range<size_t>& r,
Real delta)
665 AccessorT acc = mGrid->getAccessor();
670 const Vec3R origin = mMap.applyInverseMap(
Vec3R(0,0,0));
671 const Real Rmin = mParent.mRmin, invDx = 1 / mParent.mDx;
674 for (
size_t n = r.begin(), N = r.end();
n <
N; ++
n) {
676 typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(
n);
677 for ( ; run && iter; ++iter) {
679 mParticles.getPosRadVel(
id, pos, rad, vel);
680 const Real R0 = invDx * rad;
681 if (this->ignoreParticle(R0))
continue;
682 this->getAtt<DisableT>(
id, att);
683 const Vec3R P0 = mMap.applyInverseMap(pos);
684 const Vec3R V = mMap.applyInverseMap(vel) - origin;
685 const Real speed = V.length(), invSpeed = 1.0 / speed;
686 const Vec3R Nrml = -V * invSpeed;
689 while (run && d <= speed) {
690 run = this->makeSphere(P, R, att, acc);
691 P += 0.5 * delta * R * Nrml;
693 R = R0 - (R0 - Rmin) * d * invSpeed;
704 if (mParent.mGrainSize>0) {
705 tbb::parallel_reduce(
706 tbb::blocked_range<size_t>(0, bucketCount, mParent.mGrainSize), *
this);
708 (*this)(tbb::blocked_range<size_t>(0, bucketCount));
726 template <
bool IsMaskT = OutputIsMask>
728 makeSphere(
const Vec3R& P,
Real R,
const AttT& att, AccessorT& acc)
732 w = mParent.mHalfWidth,
740 const ValueT inside = -mGrid->background();
744 for (Coord
c = lo;
c.x() <= hi.x(); ++
c.x()) {
747 thread::cancelGroupExecution();
751 for (
c.y() = lo.y();
c.y() <= hi.y(); ++
c.y()) {
753 for (
c.z() = lo.z();
c.z() <= hi.z(); ++
c.z()) {
755 #if defined __INTEL_COMPILER
756 _Pragma(
"warning (push)")
757 _Pragma("warning (disable:186)")
759 if (x2y2z2 >= max2 || (!acc.probeValue(
c, v) && (v < ValueT(0))))
761 #if defined __INTEL_COMPILER
762 _Pragma(
"warning (pop)")
764 if (x2y2z2 <= min2) {
765 acc.setValueOff(
c, inside);
770 const ValueT d = Merge(static_cast<SdfT>(dx*(
math::Sqrt(x2y2z2)-R)), att);
771 if (d < v) acc.setValue(
c, d);
780 template <
bool IsMaskT = OutputIsMask>
782 makeSphere(
const Vec3R& p,
Real r,
const AttT& att, AccessorT& acc)
796 const std::vector<CoordBBox> padding{
797 CoordBBox(outLo.x(), outLo.y(), outLo.z(), inLo.x()-1, outHi.y(), outHi.z()),
798 CoordBBox(inHi.x()+1, outLo.y(), outLo.z(), outHi.x(), outHi.y(), outHi.z()),
799 CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), inLo.y()-1, outHi.z()),
800 CoordBBox(outLo.x(), inHi.y()+1, outLo.z(), outHi.x(), outHi.y(), outHi.z()),
801 CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), outHi.y(), inLo.z()-1),
802 CoordBBox(outLo.x(), outLo.y(), inHi.z()+1, outHi.x(), outHi.y(), outHi.z()),
804 const ValueT onValue = Merge(SdfT(1), att);
808 acc.tree().sparseFill(
CoordBBox(inLo, inHi), onValue);
811 for (
const auto& bbox: padding) {
813 thread::cancelGroupExecution();
816 const Coord &bmin = bbox.min(), &bmax = bbox.max();
819 for (c = bmin, cx = c.x(); c.x() <= bmax.x(); ++c.x(), cx += 1) {
821 for (c.y() = bmin.y(), cy = c.y(); c.y() <= bmax.y(); ++c.y(), cy += 1) {
823 for (c.z() = bmin.z(), cz = c.z(); c.z() <= bmax.z(); ++c.z(), cz += 1) {
825 if (x2y2z2 < rSquared) {
826 acc.setValue(c, onValue);
835 using FuncType =
typename std::function<void (Raster*, const tbb::blocked_range<size_t>&)>;
837 template<
typename DisableType>
839 getAtt(
size_t, AttT&)
const {}
841 template<
typename DisableType>
843 getAtt(
size_t n, AttT&
a)
const { mParticles.getAtt(n, a); }
847 Merge(T
s,
const AttT&)
const {
return s; }
851 Merge(T
s,
const AttT& a)
const {
return ValueT(s,a); }
853 ParticlesToLevelSetT& mParent;
854 const ParticleListT& mParticles;
856 const math::MapBase& mMap;
857 size_t mMinCount, mMaxCount;
860 PointPartitionerT* mPointPartitioner;
868 namespace p2ls_internal {
874 template<
typename VisibleT,
typename BlindT>
878 using type = VisibleT;
879 using VisibleType = VisibleT;
880 using BlindType = BlindT;
883 explicit BlindData(VisibleT v) : mVisible(v), mBlind(
zeroVal<BlindType>()) {}
884 BlindData(VisibleT v, BlindT
b) : mVisible(v), mBlind(b) {}
885 BlindData(
const BlindData&) =
default;
886 BlindData&
operator=(
const BlindData&) =
default;
887 const VisibleT& visible()
const {
return mVisible; }
888 const BlindT& blind()
const {
return mBlind; }
890 bool operator==(
const BlindData& rhs)
const {
return mVisible == rhs.mVisible; }
892 bool operator< (
const BlindData& rhs)
const {
return mVisible < rhs.mVisible; }
893 bool operator> (
const BlindData& rhs)
const {
return mVisible > rhs.mVisible; }
894 BlindData
operator+(
const BlindData& rhs)
const {
return BlindData(mVisible + rhs.mVisible); }
895 BlindData
operator-(
const BlindData& rhs)
const {
return BlindData(mVisible - rhs.mVisible); }
896 BlindData
operator-()
const {
return BlindData(-mVisible, mBlind); }
905 template<
typename VisibleT,
typename BlindT>
906 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
908 ostr << rhs.visible();
914 template<
typename VisibleT,
typename BlindT>
915 inline BlindData<VisibleT, BlindT>
Abs(
const BlindData<VisibleT, BlindT>&
x)
917 return BlindData<VisibleT, BlindT>(
math::Abs(x.visible()), x.blind());
922 template<
typename VisibleT,
typename BlindT,
typename T>
923 inline BlindData<VisibleT, BlindT>
924 operator+(
const BlindData<VisibleT, BlindT>& x,
const T& rhs)
926 return BlindData<VisibleT, BlindT>(x.visible() +
static_cast<VisibleT
>(rhs), x.blind());
938 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
943 "particlesToSdf requires an SDF grid with floating-point values");
947 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
955 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
960 "particlesToSdf requires an SDF grid with floating-point values");
964 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
972 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
977 "particleTrailsToSdf requires an SDF grid with floating-point values");
981 " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
989 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
994 "particlesToMask requires a boolean-valued grid");
1000 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
1005 "particlesToMask requires a boolean-valued grid");
1011 template<
typename Gr
idT,
typename ParticleListT,
typename InterrupterT>
1016 "particleTrailsToMask requires a boolean-valued grid");
1026 #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...
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
int Ceil(float x)
Return the ceiling of x.
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.
Real getVoxelSize() const
Return the size of a voxel in world units.
GLsizei const GLfloat * value
void finalize(bool prune=false)
This method syncs up the level set and attribute grids and therefore needs to be called before any of...
IMATH_HOSTDEVICE constexpr Plane3< T > operator-(const Plane3< T > &plane) IMATH_NOEXCEPT
Reflect the pla.
#define OPENVDB_LOG_WARN(mesg)
GLboolean GLboolean GLboolean GLboolean a
GLuint GLsizei GLsizei * length
#define OPENVDB_USE_VERSION_NAMESPACE
GA_API const UT_StringHolder P
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)
constexpr T zeroVal()
Return the value of type T that corresponds to zero.
int getGrainSize() const
Return the grain size used for threading.
OIIO_FORCEINLINE vbool4 operator>(const vint4 &a, const vint4 &b)
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.
float Sqrt(float x)
Return the square root of a floating-point value.
Defined various multi-threaded utility functions for trees.
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
GLboolean GLboolean GLboolean b
GA_API const UT_StringHolder transform
IMATH_HOSTDEVICE constexpr Quat< T > operator+(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
Quaterion addition.
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.
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...
LeafData & operator=(const LeafData &)=delete
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GA_API const UT_StringHolder N
InterrupterT InterrupterType
int Floor(float x)
Return the floor of x.
GLubyte GLubyte GLubyte GLubyte w
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
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.
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)
#define OPENVDB_NO_FP_EQUALITY_WARNING_END