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;
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>
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;
773 size_t mMinCount, mMaxCount;
776 PointPartitionerT* mPointPartitioner;
782 namespace p2ls_internal {
787 template<
typename VisibleT,
typename BlindT>
797 BlindData(VisibleT v, BlindT b) : mVisible(v), mBlind(b) {}
800 const VisibleT&
visible()
const {
return mVisible; }
801 const BlindT&
blind()
const {
return mBlind; }
818 template<
typename VisibleT,
typename BlindT>
819 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
821 ostr << rhs.visible();
826 template<
typename VisibleT,
typename BlindT>
840 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
Int32 x() const
Definition: Coord.h:157
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:219
double Real
Definition: Types.h:67
typename std::is_void< AttributeT >::type DisableT
Definition: ParticlesToLevelSet.h:135
int getGrainSize() const
Returns the grain-size used for multi-threading.
Definition: ParticlesToLevelSet.h:214
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
Real getRmin() const
Return the smallest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:197
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:549
void rasterizeSpheres(const ParticleListT &pa)
Rasterize a sphere per particle derived from their position and radius. All spheres are CSG unioned...
Definition: ParticlesToLevelSet.h:304
void finalize(bool prune=false)
This methods syncs up the level set and attribute grids and therefore needs to be called before any o...
Definition: ParticlesToLevelSet.h:345
Functions to efficiently perform various compositing operations on grids.
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
size_t getMinCount() const
Return number of small particles that were ignore due to Rmin.
Definition: ParticlesToLevelSet.h:204
T length() const
Length of the vector.
Definition: Vec3.h:225
typename SdfGridT::ValueType SdfType
Definition: ParticlesToLevelSet.h:139
Real getVoxelSize() const
Return the size of a voxel in world units.
Definition: ParticlesToLevelSet.h:191
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
Definition: ParticlesToLevelSet.h:194
SdfGridT SdfGridType
Definition: ParticlesToLevelSet.h:138
Spatially partitions points using a parallel radix-based sorting algorithm.
Defined various multi-threaded utility functions for trees.
void setGrainSize(int grainSize)
Set the grain-size used for multi-threading.
Definition: ParticlesToLevelSet.h:217
void setRmin(Real Rmin)
set the smallest radius allowed in voxel units
Definition: ParticlesToLevelSet.h:209
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
math::Vec3< Real > Vec3R
Definition: Types.h:79
Int32 z() const
Definition: Coord.h:159
Tag dispatch class that distinguishes shallow copy constructors from deep copy constructors.
Definition: Types.h:515
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
size_t getMaxCount() const
Return number of large particles that were ignore due to Rmax.
Definition: ParticlesToLevelSet.h:206
Definition: Exceptions.h:40
Definition: Exceptions.h:90
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:75
typename SdfGridT::template ValueConverter< AttType >::Type AttGridType
Definition: ParticlesToLevelSet.h:142
Index32 Index
Definition: Types.h:61
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:715
Definition: ParticlesToLevelSet.h:132
bool ignoredParticles() const
Return true if any particles were ignored due to their size.
Definition: ParticlesToLevelSet.h:202
Abstract base class for maps.
Definition: Maps.h:161
Int32 y() const
Definition: Coord.h:158
T zeroVal()
Return the value of type T that corresponds to zero.
Definition: Math.h:86
~ParticlesToLevelSet()
Destructor.
Definition: ParticlesToLevelSet.h:171
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:74
Type Pow2(Type x)
Return x2.
Definition: Math.h:502
void setRmax(Real Rmax)
set the largest radius allowed in voxel units
Definition: ParticlesToLevelSet.h:211
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:207
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
int Floor(float x)
Return the floor of x.
Definition: Math.h:802
InterrupterT InterrupterType
Definition: ParticlesToLevelSet.h:136
void rasterizeTrails(const ParticleListT &pa, Real delta=1.0)
Rasterize a trail per particle derived from their position, radius and velocity. Each trail is genera...
Definition: ParticlesToLevelSet.h:332
typename std::conditional< DisableT::value, size_t, AttributeT >::type AttType
Definition: ParticlesToLevelSet.h:141
AttGridType::Ptr attributeGrid()
Return a shared pointer to the grid containing the (optional) attribute.
Definition: ParticlesToLevelSet.h:188
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:518
Real getRmax() const
Return the largest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:199
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
uint32_t Index32
Definition: Types.h:59
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:810