OpenVDB  7.0.0
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 
713  bool makeSphere(const Vec3R& P, Real R, const AttT& att, AccessorT& acc)
714  {
716  if (OutputIsMask) {
717  return makeSphereMask(P, R, att, acc);
718  } else {
719  return makeNarrowBandSphere(P, R, att, acc);
720  }
722  }
723 
738  bool makeNarrowBandSphere(const Vec3R& P, Real R, const AttT& att, AccessorT& acc)
739  {
740  const Real
741  dx = mParent.mDx,
742  w = mParent.mHalfWidth,
743  max = R + w, // maximum distance in voxel units
744  max2 = math::Pow2(max), // square of maximum distance in voxel units
745  min2 = math::Pow2(math::Max(Real(0), R - w)); // square of minimum distance
746  // Bounding box of the sphere
747  const Coord
748  lo(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max)),
749  hi(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
750  const ValueT inside = -mGrid->background();
751 
752  ValueT v;
753  size_t count = 0;
754  for (Coord c = lo; c.x() <= hi.x(); ++c.x()) {
755  //only check interrupter every 32'th scan in x
756  if (!(count++ & ((1<<5)-1)) && util::wasInterrupted(mParent.mInterrupter)) {
757  tbb::task::self().cancel_group_execution();
758  return false;
759  }
760  const Real x2 = math::Pow2(c.x() - P[0]);
761  for (c.y() = lo.y(); c.y() <= hi.y(); ++c.y()) {
762  const Real x2y2 = x2 + math::Pow2(c.y() - P[1]);
763  for (c.z() = lo.z(); c.z() <= hi.z(); ++c.z()) {
764  const Real x2y2z2 = x2y2 + math::Pow2(c.z()-P[2]); // squared dist from c to P
765 #if defined __INTEL_COMPILER
766  _Pragma("warning (push)")
767  _Pragma("warning (disable:186)") // "pointless comparison of unsigned integer with zero"
768 #endif
769  if (x2y2z2 >= max2 || (!acc.probeValue(c, v) && (v < ValueT(0))))
770  continue;//outside narrow band of the particle or inside existing level set
771 #if defined __INTEL_COMPILER
772  _Pragma("warning (pop)")
773 #endif
774  if (x2y2z2 <= min2) {//inside narrow band of the particle.
775  acc.setValueOff(c, inside);
776  continue;
777  }
778  // convert signed distance from voxel units to world units
779  //const ValueT d=dx*(math::Sqrt(x2y2z2) - R);
780  const ValueT d = Merge(static_cast<SdfT>(dx*(math::Sqrt(x2y2z2)-R)), att);
781  if (d < v) acc.setValue(c, d);//CSG union
782  }//end loop over z
783  }//end loop over y
784  }//end loop over x
785  return true;
786  }
787 
790  bool makeSphereMask(const Vec3R& p, Real r, const AttT& att, AccessorT& acc)
791  {
792  const Real
793  rSquared = r * r, // sphere radius squared, in voxel units
794  inW = r / math::Sqrt(6.0); // half the width in voxel units of an inscribed cube
795  const Coord
796  // Bounding box of the sphere
797  outLo(math::Floor(p[0] - r), math::Floor(p[1] - r), math::Floor(p[2] - r)),
798  outHi(math::Ceil(p[0] + r), math::Ceil(p[1] + r), math::Ceil(p[2] + r)),
799  // Bounds of the inscribed cube
800  inLo(math::Ceil(p[0] - inW), math::Ceil(p[1] - inW), math::Ceil(p[2] - inW)),
801  inHi(math::Floor(p[0] + inW), math::Floor(p[1] + inW), math::Floor(p[2] + inW));
802  // Bounding boxes of regions comprising out - in
804  const std::vector<CoordBBox> padding{
805  CoordBBox(outLo.x(), outLo.y(), outLo.z(), inLo.x()-1, outHi.y(), outHi.z()),
806  CoordBBox(inHi.x()+1, outLo.y(), outLo.z(), outHi.x(), outHi.y(), outHi.z()),
807  CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), inLo.y()-1, outHi.z()),
808  CoordBBox(outLo.x(), inHi.y()+1, outLo.z(), outHi.x(), outHi.y(), outHi.z()),
809  CoordBBox(outLo.x(), outLo.y(), outLo.z(), outHi.x(), outHi.y(), inLo.z()-1),
810  CoordBBox(outLo.x(), outLo.y(), inHi.z()+1, outHi.x(), outHi.y(), outHi.z()),
811  };
812  const ValueT onValue = Merge(SdfT(1), att);
813 
814  // Sparsely fill the inscribed cube.
816  acc.tree().sparseFill(CoordBBox(inLo, inHi), onValue);
817 
818  // Densely fill the remaining regions.
819  for (const auto& bbox: padding) {
820  if (util::wasInterrupted(mParent.mInterrupter)) {
821  tbb::task::self().cancel_group_execution();
822  return false;
823  }
824  const Coord &bmin = bbox.min(), &bmax = bbox.max();
825  Coord c;
826  Real cx, cy, cz;
827  for (c = bmin, cx = c.x(); c.x() <= bmax.x(); ++c.x(), cx += 1) {
828  const Real x2 = math::Pow2(cx - p[0]);
829  for (c.y() = bmin.y(), cy = c.y(); c.y() <= bmax.y(); ++c.y(), cy += 1) {
830  const Real x2y2 = x2 + math::Pow2(cy - p[1]);
831  for (c.z() = bmin.z(), cz = c.z(); c.z() <= bmax.z(); ++c.z(), cz += 1) {
832  const Real x2y2z2 = x2y2 + math::Pow2(cz - p[2]);
833  if (x2y2z2 < rSquared) {
834  acc.setValue(c, onValue);
835  }
836  }
837  }
838  }
839  }
840  return true;
841  }
842 
843  using FuncType = typename std::function<void (Raster*, const tbb::blocked_range<size_t>&)>;
844 
845  template<typename DisableType>
846  typename std::enable_if<DisableType::value>::type
847  getAtt(size_t, AttT&) const {}
848 
849  template<typename DisableType>
850  typename std::enable_if<!DisableType::value>::type
851  getAtt(size_t n, AttT& a) const { mParticles.getAtt(n, a); }
852 
853  template<typename T>
854  typename std::enable_if<std::is_same<T, ValueT>::value, ValueT>::type
855  Merge(T s, const AttT&) const { return s; }
856 
857  template<typename T>
858  typename std::enable_if<!std::is_same<T, ValueT>::value, ValueT>::type
859  Merge(T s, const AttT& a) const { return ValueT(s,a); }
860 
861  ParticlesToLevelSetT& mParent;
862  const ParticleListT& mParticles;//list of particles
863  GridT* mGrid;
864  const math::MapBase& mMap;
865  size_t mMinCount, mMaxCount;//counters for ignored particles!
866  FuncType mTask;
867  const bool mIsCopy;
868  PointPartitionerT* mPointPartitioner;
869 }; // struct ParticlesToLevelSet::Raster
870 
871 
873 
874 
875 namespace p2ls_internal {
876 
877 // This is a simple type that combines a distance value and a particle
878 // attribute. It's required for attribute transfer which is defined in the
879 // Raster class above.
881 template<typename VisibleT, typename BlindT>
882 class BlindData
883 {
884 public:
885  using type = VisibleT;
886  using VisibleType = VisibleT;
887  using BlindType = BlindT;
888 
889  BlindData() {}
890  explicit BlindData(VisibleT v) : mVisible(v), mBlind(zeroVal<BlindType>()) {}
891  BlindData(VisibleT v, BlindT b) : mVisible(v), mBlind(b) {}
892  BlindData(const BlindData&) = default;
893  BlindData& operator=(const BlindData&) = default;
894  const VisibleT& visible() const { return mVisible; }
895  const BlindT& blind() const { return mBlind; }
897  bool operator==(const BlindData& rhs) const { return mVisible == rhs.mVisible; }
899  bool operator< (const BlindData& rhs) const { return mVisible < rhs.mVisible; }
900  bool operator> (const BlindData& rhs) const { return mVisible > rhs.mVisible; }
901  BlindData operator+(const BlindData& rhs) const { return BlindData(mVisible + rhs.mVisible); }
902  BlindData operator-(const BlindData& rhs) const { return BlindData(mVisible - rhs.mVisible); }
903  BlindData operator-() const { return BlindData(-mVisible, mBlind); }
904 
905 protected:
906  VisibleT mVisible;
907  BlindT mBlind;
908 };
909 
911 // Required by several of the tree nodes
912 template<typename VisibleT, typename BlindT>
913 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
914 {
915  ostr << rhs.visible();
916  return ostr;
917 }
918 
920 // Required by math::Abs
921 template<typename VisibleT, typename BlindT>
922 inline BlindData<VisibleT, BlindT> Abs(const BlindData<VisibleT, BlindT>& x)
923 {
924  return BlindData<VisibleT, BlindT>(math::Abs(x.visible()), x.blind());
925 }
926 
928 // Required to support the (zeroVal<BlindData>() + val) idiom.
929 template<typename VisibleT, typename BlindT, typename T>
930 inline BlindData<VisibleT, BlindT>
931 operator+(const BlindData<VisibleT, BlindT>& x, const T& rhs)
932 {
933  return BlindData<VisibleT, BlindT>(x.visible() + static_cast<VisibleT>(rhs), x.blind());
934 }
935 
936 } // namespace p2ls_internal
937 
938 
940 
941 
942 // The following are convenience functions for common use cases.
943 
944 template<typename GridT, typename ParticleListT, typename InterrupterT>
945 inline void
946 particlesToSdf(const ParticleListT& plist, GridT& grid, InterrupterT* interrupt)
947 {
948  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
949  "particlesToSdf requires an SDF grid with floating-point values");
950 
951  if (grid.getGridClass() != GRID_LEVEL_SET) {
952  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
953  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
954  }
955 
956  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
957  p2ls.rasterizeSpheres(plist);
958  tools::pruneLevelSet(grid.tree());
959 }
960 
961 template<typename GridT, typename ParticleListT, typename InterrupterT>
962 inline void
963 particlesToSdf(const ParticleListT& plist, GridT& grid, Real radius, InterrupterT* interrupt)
964 {
965  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
966  "particlesToSdf requires an SDF grid with floating-point values");
967 
968  if (grid.getGridClass() != GRID_LEVEL_SET) {
969  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
970  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
971  }
972 
973  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
974  p2ls.rasterizeSpheres(plist, radius);
975  tools::pruneLevelSet(grid.tree());
976 }
977 
978 template<typename GridT, typename ParticleListT, typename InterrupterT>
979 inline void
980 particleTrailsToSdf(const ParticleListT& plist, GridT& grid, Real delta, InterrupterT* interrupt)
981 {
982  static_assert(std::is_floating_point<typename GridT::ValueType>::value,
983  "particleTrailsToSdf requires an SDF grid with floating-point values");
984 
985  if (grid.getGridClass() != GRID_LEVEL_SET) {
986  OPENVDB_LOG_WARN("particlesToSdf requires a level set grid;"
987  " try Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
988  }
989 
990  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
991  p2ls.rasterizeTrails(plist, delta);
992  tools::pruneLevelSet(grid.tree());
993 }
994 
995 template<typename GridT, typename ParticleListT, typename InterrupterT>
996 inline void
997 particlesToMask(const ParticleListT& plist, GridT& grid, InterrupterT* interrupt)
998 {
999  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
1000  "particlesToMask requires a boolean-valued grid");
1001  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
1002  p2ls.rasterizeSpheres(plist);
1003  tools::prune(grid.tree());
1004 }
1005 
1006 template<typename GridT, typename ParticleListT, typename InterrupterT>
1007 inline void
1008 particlesToMask(const ParticleListT& plist, GridT& grid, Real radius, InterrupterT* interrupt)
1009 {
1010  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
1011  "particlesToMask requires a boolean-valued grid");
1012  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
1013  p2ls.rasterizeSpheres(plist, radius);
1014  tools::prune(grid.tree());
1015 }
1016 
1017 template<typename GridT, typename ParticleListT, typename InterrupterT>
1018 inline void
1019 particleTrailsToMask(const ParticleListT& plist, GridT& grid, Real delta, InterrupterT* interrupt)
1020 {
1021  static_assert(std::is_same<bool, typename GridT::ValueType>::value,
1022  "particleTrailsToMask requires a boolean-valued grid");
1023  ParticlesToLevelSet<GridT> p2ls(grid, interrupt);
1024  p2ls.rasterizeTrails(plist, delta);
1025  tools::prune(grid.tree());
1026 }
1027 
1028 } // namespace tools
1029 } // namespace OPENVDB_VERSION_NAME
1030 } // namespace openvdb
1031 
1032 #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
Definition: PointPartitioner.h:79
void setGrainSize(int grainSize)
Set the grain size used for threading.
Definition: ParticlesToLevelSet.h:233
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
std::string operator+(const std::string &s, bool)
Needed to support the (zeroVal<ValueType>() + val) idiom when ValueType is std::string.
Definition: Math.h:69
uint32_t Index32
Definition: Types.h:29
Functions to efficiently perform various compositing operations on grids.
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:803
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
size_t getMinCount() const
Return the number of particles that were ignored because they were smaller than the minimum radius...
Definition: ParticlesToLevelSet.h:224
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3683
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
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:997
Abstract base class for maps.
Definition: Maps.h:134
typename SdfGridT::template ValueConverter< AttType >::Type AttGridType
Definition: ParticlesToLevelSet.h:163
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:1127
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
Definition: ParticlesToLevelSet.h:208
Type Pow2(Type x)
Return x2.
Definition: Math.h:495
Definition: ParticlesToLevelSet.h:153
T length() const
Length of the vector.
Definition: Vec3.h:198
typename std::conditional< DisableT::value, size_t, AttributeT >::type AttType
Definition: ParticlesToLevelSet.h:162
double Real
Definition: Types.h:37
AttGridType::Ptr attributeGrid()
Return a pointer to the grid containing the optional user-defined attribute.
Definition: ParticlesToLevelSet.h:202
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:542
Spatially partitions points using a parallel radix-based sorting algorithm.
Defined various multi-threaded utility functions for trees.
Real getRmax() const
Return the largest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:216
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:192
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
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:946
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
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
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:980
~ParticlesToLevelSet()
Definition: ParticlesToLevelSet.h:187
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Real getRmin() const
Return the smallest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:211
typename SdfGridT::ValueType SdfType
Definition: ParticlesToLevelSet.h:160
Definition: Exceptions.h:13
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
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
SIMD Intrinsic Headers.
Definition: Platform.h:114
int Floor(float x)
Return the floor of x.
Definition: Math.h:795
math::Vec3< Real > Vec3R
Definition: Types.h:49
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:48
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:82
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
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:1008
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:708
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:730
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:681
Definition: Exceptions.h:63
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:471
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:47
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:496
Index32 Index
Definition: Types.h:31
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:115
bool ignoredParticles() const
Return true if any particles were ignored due to their size.
Definition: ParticlesToLevelSet.h:221
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:963
static const bool OutputIsMask
Definition: ParticlesToLevelSet.h:165
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:180
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form &#39;someVar << "some text" << ...&#39;.
Definition: logging.h:253
void setRmin(Real Rmin)
Set the smallest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:213
SdfGridT SdfGridType
Definition: ParticlesToLevelSet.h:159
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:1019
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
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
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
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
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
Definition: Types.h:454
InterrupterT InterrupterType
Definition: ParticlesToLevelSet.h:157
size_t getMaxCount() const
Return the number of particles that were ignored because they were larger than the maximum radius...
Definition: ParticlesToLevelSet.h:227
typename std::is_void< AttributeT >::type DisableT
Definition: ParticlesToLevelSet.h:156
int getGrainSize() const
Return the grain size used for threading.
Definition: ParticlesToLevelSet.h:230
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106
Real getVoxelSize() const
Return the size of a voxel in world units.
Definition: ParticlesToLevelSet.h:205
Tag dispatch class that distinguishes shallow copy constructors from deep copy constructors.
Definition: Types.h:678
void setRmax(Real Rmax)
Set the largest radius allowed in voxel units.
Definition: ParticlesToLevelSet.h:218