OpenVDB  8.1.1
ParticlesToLevelSet.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
61 
62 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
63 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
64 
65 #include <tbb/parallel_reduce.h>
66 #include <tbb/blocked_range.h>
67 #include <openvdb/Types.h>
68 #include <openvdb/Grid.h>
69 #include <openvdb/math/Math.h>
70 #include <openvdb/math/Transform.h>
72 #include <openvdb/util/logging.h>
74 #include "Composite.h" // for csgUnion()
75 #include "PointPartitioner.h"
76 #include "Prune.h"
77 #include "SignedFloodFill.h"
78 #include <functional>
79 #include <iostream>
80 #include <type_traits>
81 #include <vector>
82 
83 
84 namespace openvdb {
86 namespace OPENVDB_VERSION_NAME {
87 namespace tools {
88 
93 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
94 inline void particlesToSdf(const ParticleListT&, GridT&, InterrupterT* = nullptr);
95 
100 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
101 inline void particlesToSdf(const ParticleListT&, GridT&, Real radius, InterrupterT* = nullptr);
102 
110 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
111 inline void particleTrailsToSdf(const ParticleListT&, GridT&, Real delta=1, InterrupterT* =nullptr);
112 
117 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
118 inline void particlesToMask(const ParticleListT&, GridT&, InterrupterT* = nullptr);
119 
124 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
125 inline void particlesToMask(const ParticleListT&, GridT&, Real radius, InterrupterT* = nullptr);
126 
134 template<typename GridT, typename ParticleListT, typename InterrupterT = util::NullInterrupter>
135 inline void particleTrailsToMask(const ParticleListT&, GridT&,Real delta=1,InterrupterT* =nullptr);
136 
137 
139 
140 
141 namespace p2ls_internal {
142 // This is a simple type that combines a distance value and a particle
143 // attribute. It's required for attribute transfer which is performed
144 // in the ParticlesToLevelSet::Raster member class defined below.
146 template<typename VisibleT, typename BlindT> class BlindData;
147 }
148 
149 
150 template<typename SdfGridT,
151  typename AttributeT = void,
152  typename InterrupterT = util::NullInterrupter>
154 {
155 public:
156  using DisableT = typename std::is_void<AttributeT>::type;
157  using InterrupterType = InterrupterT;
158 
159  using SdfGridType = SdfGridT;
160  using SdfType = typename SdfGridT::ValueType;
161 
162  using AttType = typename std::conditional<DisableT::value, size_t, AttributeT>::type;
163  using AttGridType = typename SdfGridT::template ValueConverter<AttType>::Type;
164 
165  static const bool OutputIsMask = std::is_same<SdfType, bool>::value;
166 
185  explicit ParticlesToLevelSet(SdfGridT& grid, InterrupterT* interrupt = nullptr);
186 
187  ~ParticlesToLevelSet() { delete mBlindGrid; }
188 
197  void finalize(bool prune = false);
198 
202  typename AttGridType::Ptr attributeGrid() { return mAttGrid; }
203 
205  Real getVoxelSize() const { return mDx; }
206 
208  Real getHalfWidth() const { return mHalfWidth; }
209 
211  Real getRmin() const { return mRmin; }
213  void setRmin(Real Rmin) { mRmin = math::Max(Real(0),Rmin); }
214 
216  Real getRmax() const { return mRmax; }
218  void setRmax(Real Rmax) { mRmax = math::Max(mRmin,Rmax); }
219 
221  bool ignoredParticles() const { return mMinCount>0 || mMaxCount>0; }
224  size_t getMinCount() const { return mMinCount; }
227  size_t getMaxCount() const { return mMaxCount; }
228 
230  int getGrainSize() const { return mGrainSize; }
233  void setGrainSize(int grainSize) { mGrainSize = grainSize; }
234 
237  template<typename ParticleListT>
238  void rasterizeSpheres(const ParticleListT& pa);
239 
246  template<typename ParticleListT>
247  void rasterizeSpheres(const ParticleListT& pa, Real radius);
248 
264  template<typename ParticleListT>
265  void rasterizeTrails(const ParticleListT& pa, Real delta=1.0);
266 
267 private:
268  using BlindType = p2ls_internal::BlindData<SdfType, AttType>;
269  using BlindGridType = typename SdfGridT::template ValueConverter<BlindType>::Type;
270 
272  template<typename ParticleListT, typename GridT> struct Raster;
273 
274  SdfGridType* mSdfGrid;
275  typename AttGridType::Ptr mAttGrid;
276  BlindGridType* mBlindGrid;
277  InterrupterT* mInterrupter;
278  Real mDx, mHalfWidth;
279  Real mRmin, mRmax; // ignore particles outside this range of radii in voxel
280  size_t mMinCount, mMaxCount; // counters for ignored particles
281  int mGrainSize;
282 }; // class ParticlesToLevelSet
283 
284 
285 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
287 ParticlesToLevelSet(SdfGridT& grid, InterrupterT* interrupter) :
288  mSdfGrid(&grid),
289  mBlindGrid(nullptr),
290  mInterrupter(interrupter),
291  mDx(grid.voxelSize()[0]),
292  mHalfWidth(grid.background()/mDx),
293  mRmin(1.5),// corresponds to the Nyquist grid sampling frequency
294  mRmax(100.0),// corresponds to a huge particle (probably too large!)
295  mMinCount(0),
296  mMaxCount(0),
297  mGrainSize(1)
298 {
299  if (!mSdfGrid->hasUniformVoxels()) {
300  OPENVDB_THROW(RuntimeError, "ParticlesToLevelSet only supports uniform voxels!");
301  }
302  if (!DisableT::value) {
303  mBlindGrid = new BlindGridType(BlindType(grid.background()));
304  mBlindGrid->setTransform(mSdfGrid->transform().copy());
305  }
306 }
307 
308 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
309 template<typename ParticleListT>
311 rasterizeSpheres(const ParticleListT& pa)
312 {
313  if (DisableT::value) {
314  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
315  r.rasterizeSpheres();
316  } else {
317  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
318  r.rasterizeSpheres();
319  }
320 }
321 
322 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
323 template<typename ParticleListT>
325 rasterizeSpheres(const ParticleListT& pa, Real radius)
326 {
327  if (DisableT::value) {
328  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
329  r.rasterizeSpheres(radius/mDx);
330  } else {
331  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
332  r.rasterizeSpheres(radius/mDx);
333  }
334 }
335 
336 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
337 template<typename ParticleListT>
339 rasterizeTrails(const ParticleListT& pa, Real delta)
340 {
341  if (DisableT::value) {
342  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
343  r.rasterizeTrails(delta);
344  } else {
345  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
346  r.rasterizeTrails(delta);
347  }
348 }
349 
350 
351 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
352 inline void
354 {
356 
357  if (!mBlindGrid) {
358  if (prune) {
359  if (OutputIsMask) {
360  tools::prune(mSdfGrid->tree());
361  } else {
362  tools::pruneLevelSet(mSdfGrid->tree());
363  }
364  }
365  return;
366  }
367 
368  if (prune) tools::prune(mBlindGrid->tree());
369 
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;
377 
378  // Use topology copy constructors since output grids have the same topology as mBlindDataGrid
379  const BlindTreeT& blindTree = mBlindGrid->tree();
380 
381  // Create the output attribute grid.
382  typename AttTreeT::Ptr attTree(new AttTreeT(
383  blindTree, blindTree.background().blind(), openvdb::TopologyCopy()));
384  // Note this overwrites any existing attribute grids!
385  mAttGrid = typename AttGridType::Ptr(new AttGridType(attTree));
386  mAttGrid->setTransform(mBlindGrid->transform().copy());
387 
388  typename SdfTreeT::Ptr sdfTree; // the output mask or level set tree
389 
390  // Extract the attribute grid and the mask or level set grid from mBlindDataGrid.
391  if (OutputIsMask) {
392  sdfTree.reset(new SdfTreeT(blindTree,
393  /*off=*/SdfType(0), /*on=*/SdfType(1), TopologyCopy()));
394 
395  // Copy leaf voxels in parallel.
396  tree::LeafManager<AttTreeT> leafNodes(*attTree);
397  leafNodes.foreach([&](AttLeafT& attLeaf, size_t /*leafIndex*/) {
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());
402  }
403  }
404  });
405  // Copy tiles serially.
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(); });
411  }
412  } else {
413  // Here we exploit the fact that by design level sets have no active tiles.
414  // Only leaf voxels can be active.
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();
419  // Get leafnodes that were allocated during topology construction!
420  SdfLeafT* sdfLeaf = sdfTree->probeLeaf(xyz);
421  AttLeafT* attLeaf = attTree->probeLeaf(xyz);
422  // Use linear offset (vs coordinate) access for better performance!
423  typename BlindLeafT::ValueOnCIter m=leaf.cbeginValueOn();
424  if (!m) {//no active values in leaf node so copy everything
425  for (openvdb::Index k = 0; k!=BlindLeafT::SIZE; ++k) {
426  const BlindType& v = leaf.getValue(k);
427  sdfLeaf->setValueOnly(k, v.visible());
428  attLeaf->setValueOnly(k, v.blind());
429  }
430  } else {//only copy active values (using flood fill for the inactive values)
431  for(; m; ++m) {
432  const openvdb::Index k = m.pos();
433  const BlindType& v = *m;
434  sdfLeaf->setValueOnly(k, v.visible());
435  attLeaf->setValueOnly(k, v.blind());
436  }
437  }
438  }
439  tools::signedFloodFill(*sdfTree);//required since we only transferred active voxels!
440  }
441 
442  if (mSdfGrid->empty()) {
443  mSdfGrid->setTree(sdfTree);
444  } else {
445  if (OutputIsMask) {
446  mSdfGrid->tree().topologyUnion(*sdfTree);
447  tools::prune(mSdfGrid->tree());
448  } else {
449  tools::csgUnion(mSdfGrid->tree(), *sdfTree, /*prune=*/true);
450  }
451  }
452 
454 }
455 
456 
458 
459 
460 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
461 template<typename ParticleListT, typename GridT>
462 struct ParticlesToLevelSet<SdfGridT, AttributeT, InterrupterT>::Raster
463 {
464  using DisableT = typename std::is_void<AttributeT>::type;
465  using ParticlesToLevelSetT = ParticlesToLevelSet<SdfGridT, AttributeT, InterrupterT>;
466  using SdfT = typename ParticlesToLevelSetT::SdfType; // type of signed distance values
467  using AttT = typename ParticlesToLevelSetT::AttType; // type of particle attribute
468  using ValueT = typename GridT::ValueType;
469  using AccessorT = typename GridT::Accessor;
470  using TreeT = typename GridT::TreeType;
471  using LeafNodeT = typename TreeT::LeafNodeType;
472  using PointPartitionerT = PointPartitioner<Index32, LeafNodeT::LOG2DIM>;
473 
474  static const bool
475  OutputIsMask = std::is_same<SdfT, bool>::value,
476  DoAttrXfer = !DisableT::value;
477 
479  Raster(ParticlesToLevelSetT& parent, GridT* grid, const ParticleListT& particles)
480  : mParent(parent)
481  , mParticles(particles)
482  , mGrid(grid)
483  , mMap(*(mGrid->transform().baseMap()))
484  , mMinCount(0)
485  , mMaxCount(0)
486  , mIsCopy(false)
487  {
488  mPointPartitioner = new PointPartitionerT;
489  mPointPartitioner->construct(particles, mGrid->transform());
490  }
491 
493  Raster(Raster& other, tbb::split)
494  : mParent(other.mParent)
495  , mParticles(other.mParticles)
496  , mGrid(new GridT(*other.mGrid, openvdb::ShallowCopy()))
497  , mMap(other.mMap)
498  , mMinCount(0)
499  , mMaxCount(0)
500  , mTask(other.mTask)
501  , mIsCopy(true)
502  , mPointPartitioner(other.mPointPartitioner)
503  {
504  mGrid->newTree();
505  }
506 
507  virtual ~Raster()
508  {
509  // Copy-constructed Rasters own temporary grids that have to be deleted,
510  // while the original has ownership of the bucket array.
511  if (mIsCopy) {
512  delete mGrid;
513  } else {
514  delete mPointPartitioner;
515  }
516  }
517 
518  void rasterizeSpheres()
519  {
520  mMinCount = mMaxCount = 0;
521  if (mParent.mInterrupter) {
522  mParent.mInterrupter->start("Rasterizing particles to level set using spheres");
523  }
524  mTask = std::bind(&Raster::rasterSpheres, std::placeholders::_1, std::placeholders::_2);
525  this->cook();
526  if (mParent.mInterrupter) mParent.mInterrupter->end();
527  }
528 
529  void rasterizeSpheres(Real radius)
530  {
531  mMinCount = radius < mParent.mRmin ? mParticles.size() : 0;
532  mMaxCount = radius > mParent.mRmax ? mParticles.size() : 0;
533  if (mMinCount>0 || mMaxCount>0) {//skipping all particles!
534  mParent.mMinCount = mMinCount;
535  mParent.mMaxCount = mMaxCount;
536  } else {
537  if (mParent.mInterrupter) {
538  mParent.mInterrupter->start(
539  "Rasterizing particles to level set using const spheres");
540  }
541  mTask = std::bind(&Raster::rasterFixedSpheres,
542  std::placeholders::_1, std::placeholders::_2, radius);
543  this->cook();
544  if (mParent.mInterrupter) mParent.mInterrupter->end();
545  }
546  }
547 
548  void rasterizeTrails(Real delta=1.0)
549  {
550  mMinCount = mMaxCount = 0;
551  if (mParent.mInterrupter) {
552  mParent.mInterrupter->start("Rasterizing particles to level set using trails");
553  }
554  mTask = std::bind(&Raster::rasterTrails,
555  std::placeholders::_1, std::placeholders::_2, delta);
556  this->cook();
557  if (mParent.mInterrupter) mParent.mInterrupter->end();
558  }
559 
561  void operator()(const tbb::blocked_range<size_t>& r)
562  {
563  assert(mTask);
564  mTask(this, r);
565  mParent.mMinCount = mMinCount;
566  mParent.mMaxCount = mMaxCount;
567  }
568 
570  void join(Raster& other)
571  {
573  if (OutputIsMask) {
574  if (DoAttrXfer) {
575  tools::compMax(*mGrid, *other.mGrid);
576  } else {
577  mGrid->topologyUnion(*other.mGrid);
578  }
579  } else {
580  tools::csgUnion(*mGrid, *other.mGrid, /*prune=*/true);
581  }
583  mMinCount += other.mMinCount;
584  mMaxCount += other.mMaxCount;
585  }
586 
587 private:
589  Raster& operator=(const Raster&) { return *this; }
590 
592  bool ignoreParticle(Real R)
593  {
594  if (R < mParent.mRmin) {// below the cutoff radius
595  ++mMinCount;
596  return true;
597  }
598  if (R > mParent.mRmax) {// above the cutoff radius
599  ++mMaxCount;
600  return true;
601  }
602  return false;
603  }
604 
607  void rasterSpheres(const tbb::blocked_range<size_t>& r)
608  {
609  AccessorT acc = mGrid->getAccessor(); // local accessor
610  bool run = true;
611  const Real invDx = 1 / mParent.mDx;
612  AttT att;
613  Vec3R pos;
614  Real rad;
615 
616  // Loop over buckets
617  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
618  // Loop over particles in bucket n.
619  typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
620  for ( ; run && iter; ++iter) {
621  const Index32& id = *iter;
622  mParticles.getPosRad(id, pos, rad);
623  const Real R = invDx * rad;// in voxel units
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);
628  }//end loop over particles
629  }//end loop over buckets
630  }
631 
635  void rasterFixedSpheres(const tbb::blocked_range<size_t>& r, Real R)
636  {
637  AccessorT acc = mGrid->getAccessor(); // local accessor
638  AttT att;
639  Vec3R pos;
640 
641  // Loop over buckets
642  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
643  // Loop over particles in bucket n.
644  for (auto iter = mPointPartitioner->indices(n); iter; ++iter) {
645  const Index32& id = *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);
650  }
651  }
652  }
653 
657  void rasterTrails(const tbb::blocked_range<size_t>& r, Real delta)
658  {
659  AccessorT acc = mGrid->getAccessor(); // local accessor
660  bool run = true;
661  AttT att;
662  Vec3R pos, vel;
663  Real rad;
664  const Vec3R origin = mMap.applyInverseMap(Vec3R(0,0,0));
665  const Real Rmin = mParent.mRmin, invDx = 1 / mParent.mDx;
666 
667  // Loop over buckets
668  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
669  // Loop over particles in bucket n.
670  typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
671  for ( ; run && iter; ++iter) {
672  const Index32& id = *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; // exclude translation
679  const Real speed = V.length(), invSpeed = 1.0 / speed;
680  const Vec3R Nrml = -V * invSpeed; // inverse normalized direction
681  Vec3R P = P0; // local position of instance
682  Real R = R0, d = 0; // local radius and length of trail
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; // adaptive offset along inverse velocity direction
686  d = (P - P0).length(); // current length of trail
687  R = R0 - (R0 - Rmin) * d * invSpeed; // R = R0 -> mRmin(e.g. 1.5)
688  }//end loop over sphere instances
689  }//end loop over particles
690  }//end loop over buckets
691  }
692 
693  void cook()
694  {
695  // parallelize over the point buckets
696  const Index32 bucketCount = Index32(mPointPartitioner->size());
697 
698  if (mParent.mGrainSize>0) {
699  tbb::parallel_reduce(
700  tbb::blocked_range<size_t>(0, bucketCount, mParent.mGrainSize), *this);
701  } else {
702  (*this)(tbb::blocked_range<size_t>(0, bucketCount));
703  }
704  }
705 
720  template <bool IsMaskT = OutputIsMask>
721  typename std::enable_if<!IsMaskT, bool>::type
722  makeSphere(const Vec3R& P, Real R, const AttT& att, AccessorT& acc)
723  {
724  const Real
725  dx = mParent.mDx,
726  w = mParent.mHalfWidth,
727  max = R + w, // maximum distance in voxel units
728  max2 = math::Pow2(max), // square of maximum distance in voxel units
729  min2 = math::Pow2(math::Max(Real(0), R - w)); // square of minimum distance
730  // Bounding box of the sphere
731  const Coord
732  lo(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max)),
733  hi(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
734  const ValueT inside = -mGrid->background();
735 
736  ValueT v;
737  size_t count = 0;
738  for (Coord c = lo; c.x() <= hi.x(); ++c.x()) {
739  //only check interrupter every 32'th scan in x
740  if (!(count++ & ((1<<5)-1)) && util::wasInterrupted(mParent.mInterrupter)) {
741  tbb::task::self().cancel_group_execution();
742  return false;
743  }
744  const Real x2 = math::Pow2(c.x() - P[0]);
745  for (c.y() = lo.y(); c.y() <= hi.y(); ++c.y()) {
746  const Real x2y2 = x2 + math::Pow2(c.y() - P[1]);
747  for (c.z() = lo.z(); c.z() <= hi.z(); ++c.z()) {
748  const Real x2y2z2 = x2y2 + math::Pow2(c.z()-P[2]); // squared dist from c to P
749 #if defined __INTEL_COMPILER
750  _Pragma("warning (push)")
751  _Pragma("warning (disable:186)") // "pointless comparison of unsigned integer with zero"
752 #endif
753  if (x2y2z2 >= max2 || (!acc.probeValue(c, v) && (v < ValueT(0))))
754  continue;//outside narrow band of the particle or inside existing level set
755 #if defined __INTEL_COMPILER
756  _Pragma("warning (pop)")
757 #endif
758  if (x2y2z2 <= min2) {//inside narrow band of the particle.
759  acc.setValueOff(c, inside);
760  continue;
761  }
762  // convert signed distance from voxel units to world units
763  //const ValueT d=dx*(math::Sqrt(x2y2z2) - R);
764  const ValueT d = Merge(static_cast<SdfT>(dx*(math::Sqrt(x2y2z2)-R)), att);
765  if (d < v) acc.setValue(c, d);//CSG union
766  }//end loop over z
767  }//end loop over y
768  }//end loop over x
769  return true;
770  }
771 
774  template <bool IsMaskT = OutputIsMask>
775  typename std::enable_if<IsMaskT, bool>::type
776  makeSphere(const Vec3R& p, Real r, const AttT& att, AccessorT& acc)
777  {
778  const Real
779  rSquared = r * r, // sphere radius squared, in voxel units
780  inW = r / math::Sqrt(6.0); // half the width in voxel units of an inscribed cube
781  const Coord
782  // Bounding box of the sphere
783  outLo(math::Floor(p[0] - r), math::Floor(p[1] - r), math::Floor(p[2] - r)),
784  outHi(math::Ceil(p[0] + r), math::Ceil(p[1] + r), math::Ceil(p[2] + r)),
785  // Bounds of the inscribed cube
786  inLo(math::Ceil(p[0] - inW), math::Ceil(p[1] - inW), math::Ceil(p[2] - inW)),
787  inHi(math::Floor(p[0] + inW), math::Floor(p[1] + inW), math::Floor(p[2] + inW));
788  // Bounding boxes of regions comprising out - in
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()),
797  };
798  const ValueT onValue = Merge(SdfT(1), att);
799 
800  // Sparsely fill the inscribed cube.
802  acc.tree().sparseFill(CoordBBox(inLo, inHi), onValue);
803 
804  // Densely fill the remaining regions.
805  for (const auto& bbox: padding) {
806  if (util::wasInterrupted(mParent.mInterrupter)) {
807  tbb::task::self().cancel_group_execution();
808  return false;
809  }
810  const Coord &bmin = bbox.min(), &bmax = bbox.max();
811  Coord c;
812  Real cx, cy, cz;
813  for (c = bmin, cx = c.x(); c.x() <= bmax.x(); ++c.x(), cx += 1) {
814  const Real x2 = math::Pow2(cx - p[0]);
815  for (c.y() = bmin.y(), cy = c.y(); c.y() <= bmax.y(); ++c.y(), cy += 1) {
816  const Real x2y2 = x2 + math::Pow2(cy - p[1]);
817  for (c.z() = bmin.z(), cz = c.z(); c.z() <= bmax.z(); ++c.z(), cz += 1) {
818  const Real x2y2z2 = x2y2 + math::Pow2(cz - p[2]);
819  if (x2y2z2 < rSquared) {
820  acc.setValue(c, onValue);
821  }
822  }
823  }
824  }
825  }
826  return true;
827  }
828 
829  using FuncType = typename std::function<void (Raster*, const tbb::blocked_range<size_t>&)>;
830 
831  template<typename DisableType>
832  typename std::enable_if<DisableType::value>::type
833  getAtt(size_t, AttT&) const {}
834 
835  template<typename DisableType>
836  typename std::enable_if<!DisableType::value>::type
837  getAtt(size_t n, AttT& a) const { mParticles.getAtt(n, a); }
838 
839  template<typename T>
840  typename std::enable_if<std::is_same<T, ValueT>::value, ValueT>::type
841  Merge(T s, const AttT&) const { return s; }
842 
843  template<typename T>
844  typename std::enable_if<!std::is_same<T, ValueT>::value, ValueT>::type
845  Merge(T s, const AttT& a) const { return ValueT(s,a); }
846 
847  ParticlesToLevelSetT& mParent;
848  const ParticleListT& mParticles;//list of particles
849  GridT* mGrid;
850  const math::MapBase& mMap;
851  size_t mMinCount, mMaxCount;//counters for ignored particles!
852  FuncType mTask;
853  const bool mIsCopy;
854  PointPartitionerT* mPointPartitioner;
855 }; // struct ParticlesToLevelSet::Raster
856 
857 
859 
860 
861 namespace p2ls_internal {
862 
863 // This is a simple type that combines a distance value and a particle
864 // attribute. It's required for attribute transfer which is defined in the
865 // Raster class above.
867 template<typename VisibleT, typename BlindT>
868 class BlindData
869 {
870 public:
871  using type = VisibleT;
872  using VisibleType = VisibleT;
873  using BlindType = BlindT;
874 
875  BlindData() {}
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); }
890 
891 protected:
892  VisibleT mVisible;
893  BlindT mBlind;
894 };
895 
897 // Required by several of the tree nodes
898 template<typename VisibleT, typename BlindT>
899 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
900 {
901  ostr << rhs.visible();
902  return ostr;
903 }
904 
906 // Required by math::Abs
907 template<typename VisibleT, typename BlindT>
908 inline BlindData<VisibleT, BlindT> Abs(const BlindData<VisibleT, BlindT>& x)
909 {
910  return BlindData<VisibleT, BlindT>(math::Abs(x.visible()), x.blind());
911 }
912 
914 // Required to support the (zeroVal<BlindData>() + val) idiom.
915 template<typename VisibleT, typename BlindT, typename T>
916 inline BlindData<VisibleT, BlindT>
917 operator+(const BlindData<VisibleT, BlindT>& x, const T& rhs)
918 {
919  return BlindData<VisibleT, BlindT>(x.visible() + static_cast<VisibleT>(rhs), x.blind());
920 }
921 
922 } // namespace p2ls_internal
923 
924 
926 
927 
928 // The following are convenience functions for common use cases.
929 
930 template<typename GridT, typename ParticleListT, typename InterrupterT>
931 inline void
932 particlesToSdf(const ParticleListT& plist, GridT& grid, InterrupterT* interrupt)
933 {
934  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
935  "particlesToSdf requires an SDF grid with floating-point values");
936 
937  if (grid.getGridClass() != GRID_LEVEL_SET) {
938  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
939  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
940  }
941 
942  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
943  p2ls.rasterizeSpheres(plist);
944  tools::pruneLevelSet(grid.tree());
945 }
946 
947 template<typename GridT, typename ParticleListT, typename InterrupterT>
948 inline void
949 particlesToSdf(const ParticleListT& plist, GridT& grid, Real radius, InterrupterT* interrupt)
950 {
951  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
952  "particlesToSdf requires an SDF grid with floating-point values");
953 
954  if (grid.getGridClass() != GRID_LEVEL_SET) {
955  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
956  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
957  }
958 
959  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
960  p2ls.rasterizeSpheres(plist, radius);
961  tools::pruneLevelSet(grid.tree());
962 }
963 
964 template<typename GridT, typename ParticleListT, typename InterrupterT>
965 inline void
966 particleTrailsToSdf(const ParticleListT& plist, GridT& grid, Real delta, InterrupterT* interrupt)
967 {
968  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
969  "particleTrailsToSdf requires an SDF grid with floating-point values");
970 
971  if (grid.getGridClass() != GRID_LEVEL_SET) {
972  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
973  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
974  }
975 
976  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
977  p2ls.rasterizeTrails(plist, delta);
978  tools::pruneLevelSet(grid.tree());
979 }
980 
981 template<typename GridT, typename ParticleListT, typename InterrupterT>
982 inline void
983 particlesToMask(const ParticleListT& plist, GridT& grid, InterrupterT* interrupt)
984 {
985  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
986  "particlesToMask requires a boolean-valued grid");
987  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
988  p2ls.rasterizeSpheres(plist);
989  tools::prune(grid.tree());
990 }
991 
992 template<typename GridT, typename ParticleListT, typename InterrupterT>
993 inline void
994 particlesToMask(const ParticleListT& plist, GridT& grid, Real radius, InterrupterT* interrupt)
995 {
996  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
997  "particlesToMask requires a boolean-valued grid");
998  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
999  p2ls.rasterizeSpheres(plist, radius);
1000  tools::prune(grid.tree());
1001 }
1002 
1003 template<typename GridT, typename ParticleListT, typename InterrupterT>
1004 inline void
1005 particleTrailsToMask(const ParticleListT& plist, GridT& grid, Real delta, InterrupterT* interrupt)
1006 {
1007  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
1008  "particleTrailsToMask requires a boolean-valued grid");
1009  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
1010  p2ls.rasterizeTrails(plist, delta);
1011  tools::prune(grid.tree());
1012 }
1013 
1014 } // namespace tools
1015 } // namespace OPENVDB_VERSION_NAME
1016 } // namespace openvdb
1017 
1018 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
void rasterizeSpheres(const ParticleListT &pa)
Rasterize each particle as a sphere with the particle&#39;s position and radius.
Definition: ParticlesToLevelSet.h:311
size_t getMinCount() const
Return the number of particles that were ignored because they were smaller than the minimum radius...
Definition: ParticlesToLevelSet.h:224
void particlesToSdf(const ParticleListT &, GridT &, InterrupterT *=nullptr)
Populate a scalar, floating-point grid with CSG-unioned level set spheres described by the given part...
Definition: ParticlesToLevelSet.h:932
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
uint32_t Index32
Definition: openvdb/Types.h:48
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
void prune(TreeT &tree, typename TreeT::ValueType tolerance=zeroVal< typename TreeT::ValueType >(), bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with tiles any nodes whose values are all the same...
Definition: Prune.h:334
Functions to efficiently perform various compositing operations on grids.
void compMax(GridOrTreeT &a, GridOrTreeT &b)
Given grids A and B, compute max(a, b) per voxel (using sparse traversal). Store the result in the A ...
Definition: Composite.h:752
int Floor(float x)
Return the floor of x.
Definition: Math.h:851
InterrupterT InterrupterType
Definition: ParticlesToLevelSet.h:157
Real getVoxelSize() const
Return the size of a voxel in world units.
Definition: ParticlesToLevelSet.h:205
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:266
void particlesToMask(const ParticleListT &, GridT &, InterrupterT *=nullptr)
Activate a boolean grid wherever it intersects the spheres described by the given particle positions ...
Definition: ParticlesToLevelSet.h:983
void csgUnion(GridOrTreeT &a, GridOrTreeT &b, bool prune=true)
Given two level set grids, replace the A grid with the union of A and B.
Definition: Composite.h:885
void particleTrailsToSdf(const ParticleListT &, GridT &, Real delta=1, InterrupterT *=nullptr)
Populate a scalar, floating-point grid with CSG-unioned trails of level set spheres with decreasing r...
Definition: ParticlesToLevelSet.h:966
Real getRmin() const
Return the smallest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:211
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
Definition: ParticlesToLevelSet.h:208
typename std::is_void< AttributeT >::type DisableT
Definition: ParticlesToLevelSet.h:156
void pruneLevelSet(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing nodes whose values are all inactive with inactive ...
Definition: Prune.h:389
math::Vec3< Real > Vec3R
Definition: openvdb/Types.h:68
typename SdfGridT::ValueType SdfType
Definition: ParticlesToLevelSet.h:160
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 threading.
Definition: ParticlesToLevelSet.h:233
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Index32 Index
Definition: openvdb/Types.h:50
typename SdfGridT::template ValueConverter< AttType >::Type AttGridType
Definition: ParticlesToLevelSet.h:163
void particlesToSdf(const ParticleListT &, GridT &, Real radius, InterrupterT *=nullptr)
Populate a scalar, floating-point grid with fixed-size, CSG-unioned level set spheres described by th...
Definition: ParticlesToLevelSet.h:949
typename std::conditional< DisableT::value, size_t, AttributeT >::type AttType
Definition: ParticlesToLevelSet.h:162
void particlesToMask(const ParticleListT &, GridT &, Real radius, InterrupterT *=nullptr)
Activate a boolean grid wherever it intersects the fixed-size spheres described by the given particle...
Definition: ParticlesToLevelSet.h:994
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
void setRmax(Real Rmax)
Set the largest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:218
void setRmin(Real Rmin)
Set the smallest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:213
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:189
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Definition: openvdb/Exceptions.h:63
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:859
SdfGridT SdfGridType
Definition: ParticlesToLevelSet.h:159
Definition: openvdb/Exceptions.h:13
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
SIMD Intrinsic Headers.
Definition: Platform.h:115
void particleTrailsToMask(const ParticleListT &, GridT &, Real delta=1, InterrupterT *=nullptr)
Activate a boolean grid wherever it intersects trails of spheres with decreasing radius, where the starting position and radius and the direction of each trail is given by particle attributes.
Definition: ParticlesToLevelSet.h:1005
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:48
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
Real getRmax() const
Return the largest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:216
~ParticlesToLevelSet()
Definition: ParticlesToLevelSet.h:187
double Real
Definition: openvdb/Types.h:56
Vec3< typename promote< T, Coord::ValueType >::type > operator-(const Vec3< T > &v0, const Coord &v1)
Allow a Coord to be subtracted from a Vec3.
Definition: Coord.h:551
Abstract base class for maps.
Definition: Maps.h:134
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3705
Tag dispatch class that distinguishes shallow copy constructors from deep copy constructors.
Definition: openvdb/Types.h:557
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:47
bool ignoredParticles() const
Return true if any particles were ignored due to their size.
Definition: ParticlesToLevelSet.h:221
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Definition: openvdb/Types.h:333
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:116
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:477
Definition: ParticlesToLevelSet.h:153
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form &#39;someVar << "some text" << ...&#39;.
Definition: logging.h:253
int getGrainSize() const
Return the grain size used for threading.
Definition: ParticlesToLevelSet.h:230
void rasterizeTrails(const ParticleListT &pa, Real delta=1.0)
Rasterize each particle as a trail comprising the CSG union of spheres of decreasing radius...
Definition: ParticlesToLevelSet.h:339
T length() const
Length of the vector.
Definition: Vec3.h:204
size_t getMaxCount() const
Return the number of particles that were ignored because they were larger than the maximum radius...
Definition: ParticlesToLevelSet.h:227
void finalize(bool prune=false)
This method syncs up the level set and attribute grids and therefore needs to be called before any of...
Definition: ParticlesToLevelSet.h:353
AttGridType::Ptr attributeGrid()
Return a pointer to the grid containing the optional user-defined attribute.
Definition: ParticlesToLevelSet.h:202
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:598
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Definition: PointPartitioner.h:76
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.
Definition: LeafManager.h:483
static const bool OutputIsMask
Definition: ParticlesToLevelSet.h:165
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:201