OpenVDB  1.1.0
ParticlesToLevelSet.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
81 
82 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
83 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
84 
85 #include <tbb/parallel_reduce.h>
86 #include <tbb/blocked_range.h>
87 #include <boost/bind.hpp>
88 #include <boost/function.hpp>
89 #include <openvdb/util/Util.h>
90 #include <openvdb/Types.h>
91 #include <openvdb/Grid.h>
92 #include <openvdb/math/Math.h>
93 #include <openvdb/math/Transform.h>
94 #include <openvdb/util/NullInterrupter.h>
95 #include "Composite.h" // for csgUnion()
96 
97 
98 namespace openvdb {
100 namespace OPENVDB_VERSION_NAME {
101 namespace tools {
102 namespace local {
107 template <typename T>
108 struct DualTrait
109 {
110  static T merge(T dist, Index32) { return dist; }
111  static T split(T dist) { return dist; }
112 };
113 }// namespace local
114 
115 template<typename GridT,
116  typename ParticleListT,
117  typename InterruptT=util::NullInterrupter,
118  typename RealT = typename GridT::ValueType>
120 {
121 public:
136  ParticlesToLevelSet(GridT& grid, InterruptT* interrupt = NULL) :
137  mGrid(&grid),
138  mPa(NULL),
139  mDx(grid.transform().voxelSize()[0]),
140  mHalfWidth(local::DualTrait<ValueT>::split(grid.background()) / mDx),
141  mRmin(1.5),// corresponds to the Nyquist grid sampling frequency
142  mRmax(100.0f),// corresponds to a huge particle (probably too large!)
143  mGrainSize(1),
144  mInterrupter(interrupt),
145  mMinCount(0),
146  mMaxCount(0),
147  mIsSlave(false)
148  {
149  if ( !mGrid->hasUniformVoxels() ) {
151  "The transform must have uniform scale for ParticlesToLevelSet to function!");
152  }
153  if (mGrid->getGridClass() != GRID_LEVEL_SET) {
155  "ParticlesToLevelSet only supports level sets!"
156  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
157  }
159  //mGrid->newTree();
160  }
163  mGrid(new GridT(*other.mGrid, openvdb::ShallowCopy())),
164  mPa(other.mPa),
165  mDx(other.mDx),
166  mHalfWidth(other.mHalfWidth),
167  mRmin(other.mRmin),
168  mRmax(other.mRmax),
169  mGrainSize(other.mGrainSize),
170  mTask(other.mTask),
171  mInterrupter(other.mInterrupter),
172  mMinCount(0),
173  mMaxCount(0),
174  mIsSlave(true)
175  {
176  mGrid->newTree();
177  }
178  virtual ~ParticlesToLevelSet() { if (mIsSlave) delete mGrid; }
179 
181  RealT getHalfWidth() const { return mHalfWidth; }
182 
184  RealT getVoxelSize() const { return mDx; }
185 
187  RealT getRmin() const { return mRmin; }
189  RealT getRmax() const { return mRmax; }
190 
192  bool ignoredParticles() const { return mMinCount>0 || mMaxCount>0; }
194  size_t getMinCount() const { return mMinCount; }
196  size_t getMaxCount() const { return mMaxCount; }
197 
199  void setRmin(RealT Rmin) { mRmin = math::Max(RealT(0),Rmin); }
201  void setRmax(RealT Rmax) { mRmax = math::Max(mRmin,Rmax); }
202 
204  int getGrainSize() const { return mGrainSize; }
207  void setGrainSize(int grainSize) { mGrainSize = grainSize; }
208 
212  void rasterizeSpheres(const ParticleListT& pa)
213  {
214  mPa = &pa;
215  if (mInterrupter) mInterrupter->start("Rasterizing particles to level set using spheres");
216  mTask = boost::bind(&ParticlesToLevelSet::rasterSpheres, _1, _2);
217  this->cook();
218  if (mInterrupter) mInterrupter->end();
219  }
220 
236  void rasterizeTrails(const ParticleListT& pa, Real delta=1.0)
237  {
238  mPa = &pa;
239  if (mInterrupter) mInterrupter->start("Rasterizing particles to level set using trails");
240  mTask = boost::bind(&ParticlesToLevelSet::rasterTrails, _1, _2, RealT(delta));
241  this->cook();
242  if (mInterrupter) mInterrupter->end();
243  }
244 
245  // ========> DO NOT CALL ANY OF THE PUBLIC METHODS BELOW! <=============
246 
250  void operator()(const tbb::blocked_range<size_t>& r)
251  {
252  if (!mTask) {
253  OPENVDB_THROW(ValueError, "mTask is undefined - don't call operator() directly!");
254  }
255  mTask(this, r);
256  }
257 
261  void join(ParticlesToLevelSet& other)
262  {
263  tools::csgUnion(*mGrid, *other.mGrid, /*prune=*/true);
264  mMinCount += other.mMinCount;
265  mMaxCount += other.mMaxCount;
266  }
267 
268 private:
269  typedef typename GridT::ValueType ValueT;
270  typedef typename GridT::Accessor Accessor;
271  typedef typename boost::function<void (ParticlesToLevelSet*,
272  const tbb::blocked_range<size_t>&)> FuncType;
273  GridT* mGrid;
274  const ParticleListT* mPa;//list of particles
275  const RealT mDx;//size of voxel in world units
276  const RealT mHalfWidth;//half width of narrow band LS in voxels
277  RealT mRmin;//ignore particles smaller than this radius in voxels
278  RealT mRmax;//ignore particles larger than this radius in voxels
279  int mGrainSize;
280  FuncType mTask;
281  InterruptT* mInterrupter;
282  size_t mMinCount, mMaxCount;//counters for ignored particles!
283  const bool mIsSlave;
284 
285  void cook()
286  {
287  if (mGrainSize>0) {
288  tbb::parallel_reduce(tbb::blocked_range<size_t>(0,mPa->size(),mGrainSize), *this);
289  } else {
290  (*this)(tbb::blocked_range<size_t>(0, mPa->size()));
291  }
292  }
294  inline bool ignoreParticle(RealT R)
295  {
296  if (R<mRmin) {// below the cutoff radius
297  ++mMinCount;
298  return true;
299  }
300  if (R>mRmax) {// above the cutoff radius
301  ++mMaxCount;
302  return true;
303  }
304  return false;
305  }
321  inline bool rasterSphere(const Vec3R &P, RealT R, Index32 id, Accessor& accessor)
322  {
323  const ValueT inside = -mGrid->background();
324  const RealT dx = mDx;
325  const RealT max = R + mHalfWidth;// maximum distance in voxel units
326  const Coord a(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max));
327  const Coord b(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
328  const RealT max2 = math::Pow2(max);//square of maximum distance in voxel units
329  const RealT min2 = math::Pow2(math::Max(RealT(0), R - mHalfWidth));//square of minimum distance
330  ValueT v;
331  size_t count = 0;
332  for ( Coord c = a; c.x() <= b.x(); ++c.x() ) {
333  //only check interrupter every 32'th scan in x
334  if (!(count++ & (1<<5)-1) && util::wasInterrupted(mInterrupter)) {
335  tbb::task::self().cancel_group_execution();
336  return false;
337  }
338  RealT x2 = math::Pow2( c.x() - P[0] );
339  for ( c.y() = a.y(); c.y() <= b.y(); ++c.y() ) {
340  RealT x2y2 = x2 + math::Pow2( c.y() - P[1] );
341  for ( c.z() = a.z(); c.z() <= b.z(); ++c.z() ) {
342  RealT x2y2z2 = x2y2 + math::Pow2( c.z() - P[2] );//square distance from c to P
343  if ( x2y2z2 >= max2 || (!accessor.probeValue(c,v) && v<ValueT(0)) )
344  continue;//outside narrow band of the particle or inside existing ls
345  if ( x2y2z2 <= min2 ) {//inside narrowband of the particle.
346  accessor.setValueOff(c, inside);
347  continue;
348  }
349  // distance in world units
350  const ValueT d = local::DualTrait<ValueT>::merge(dx*(math::Sqrt(x2y2z2) - R), id);
351  if (d < v) accessor.setValue(c, d);//CSG union
352  }//end loop over z
353  }//end loop over y
354  }//end loop over x
355  return true;
356  }
357 
361  void rasterSpheres(const tbb::blocked_range<size_t> &r)
362  {
363  Accessor accessor = mGrid->getAccessor(); // local accessor
364  const RealT inv_dx = RealT(1)/mDx;
365  bool run = true;
366  for (Index32 id = r.begin(), e=r.end(); run && id != e; ++id) {
367  const RealT R = inv_dx*mPa->radius(id);// in voxel units
368  if (this->ignoreParticle(R)) continue;
369  const Vec3R P = mGrid->transform().worldToIndex(mPa->pos(id));
370  run = this->rasterSphere(P, R, id, accessor);
371  }//end loop over particles
372  }
373 
384  void rasterTrails(const tbb::blocked_range<size_t> &r,
385  RealT delta)//scale distance between instances (eg 1.0)
386  {
387  Accessor accessor = mGrid->getAccessor(); // local accessor
388  const RealT inv_dx = RealT(1)/mDx, Rmin = mRmin;
389  bool run = true;
390  for (Index32 id = r.begin(), e=r.end(); run && id != e; ++id) {
391  const RealT R0 = inv_dx*mPa->radius(id);
392  if (this->ignoreParticle(R0)) continue;
393  const Vec3R P0 = mGrid->transform().worldToIndex(mPa->pos(id)),
394  V = inv_dx*mPa->vel(id);
395  const RealT speed = V.length(), inv_speed=1.0/speed;
396  const Vec3R N = -V*inv_speed;// inverse normalized direction
397  Vec3R P = P0;// local position of instance
398  RealT R = R0, d=0;// local radius and length of trail
399  for (size_t m=0; run && d < speed ; ++m) {
400  run = this->rasterSphere(P, R, id, accessor);
401  P += 0.5*delta*R*N;// adaptive offset along inverse velocity direction
402  d = (P-P0).length();// current length of trail
403  R = R0-(R0-Rmin)*d*inv_speed;// R = R0 -> mRmin(e.g. 1.5)
404  }//loop over sphere instances
405  }//end loop over particles
406  }
407 };//end of ParticlesToLevelSet class
408 
410 namespace local {
411 // This is a simple type that combines a distance value and a particle
412 // id. It's required for attribute transfer which is defined in the
413 // ParticlesToLevelSetAndId class below.
414 template <typename RealT>
415 class Dual
416 {
417 public:
418  explicit Dual() : mId(util::INVALID_IDX) {}
419  explicit Dual(RealT d) : mDist(d), mId(util::INVALID_IDX) {}
420  explicit Dual(RealT d, Index32 id) : mDist(d), mId(id) {}
421  Dual& operator=(const Dual& rhs) { mDist = rhs.mDist; mId = rhs.mId; return *this;}
422  bool isIdValid() const { return mId != util::INVALID_IDX; }
423  Index32 id() const { return mId; }
424  RealT dist() const { return mDist; }
425  bool operator!=(const Dual& rhs) const { return mDist != rhs.mDist; }
426  bool operator==(const Dual& rhs) const { return mDist == rhs.mDist; }
427  bool operator< (const Dual& rhs) const { return mDist < rhs.mDist; };
428  bool operator<=(const Dual& rhs) const { return mDist <= rhs.mDist; };
429  bool operator> (const Dual& rhs) const { return mDist > rhs.mDist; };
430  Dual operator+ (const Dual& rhs) const { return Dual(mDist+rhs.mDist); };
431  Dual operator+ (const RealT& rhs) const { return Dual(mDist+rhs); };
432  Dual operator- (const Dual& rhs) const { return Dual(mDist-rhs.mDist); };
433  Dual operator-() const { return Dual(-mDist); }
434 protected:
435  RealT mDist;
437 };
438 // Required by several of the tree nodes
439 template <typename RealT>
440 inline std::ostream& operator<<(std::ostream& ostr, const Dual<RealT>& rhs)
441 {
442  ostr << rhs.dist();
443  return ostr;
444 }
445 // Required by math::Abs
446 template <typename RealT>
447 inline Dual<RealT> Abs(const Dual<RealT>& x)
448 {
449  return Dual<RealT>(math::Abs(x.dist()),x.id());
450 }
451 // Specialization of trait class used to merge and split a distance and a particle id
452 template <typename T>
453 struct DualTrait<Dual<T> >
454 {
455  static Dual<T> merge(T dist, Index32 id) { return Dual<T>(dist, id); }
456  static T split(Dual<T> dual) { return dual.dist(); }
457 };
458 }// local namespace
460 
467 template<typename LevelSetGridT,
468  typename ParticleListT,
469  typename InterruptT = util::NullInterrupter>
471 {
472 public:
473  typedef typename LevelSetGridT::ValueType RealT;
474  typedef typename local::Dual<RealT> DualT;
475  typedef typename LevelSetGridT::TreeType RealTreeT;
476  typedef typename RealTreeT::template ValueConverter<DualT>::Type DualTreeT;
480 
481  ParticlesToLevelSetAndId(LevelSetGridT& ls, InterruptT* interrupter = NULL) :
482  mRealGrid(ls),
483  mDualGrid(DualT(ls.background()))
484  {
485  mDualGrid.setGridClass(ls.getGridClass());
486  mDualGrid.setTransform(ls.transformPtr());
487  mRaster = new RasterT(mDualGrid, interrupter);
488  }
489 
490  virtual ~ParticlesToLevelSetAndId() { delete mRaster; }
491 
493  RealT getHalfWidth() const { return mRaster->getHalfWidth(); }
494 
496  RealT getVoxelSize() const { return mRaster->getVoxelSize(); }
497 
499  bool ignoredParticles() const { return mRaster->ignoredParticles(); }
501  size_t getMinCount() const { return mRaster->getMinCount(); }
503  size_t getMaxCount() const { return mRaster->getMaxCount(); }
504 
506  RealT getRmin() const { return mRaster->getRmin(); }
508  RealT getRmax() const { return mRaster->getRmax(); }
509 
511  void setRmin(RealT Rmin) { mRaster->setRmin(Rmin); }
513  void setRmax(RealT Rmax) { mRaster->setRmax(Rmax); }
514 
516  int getGrainSize() const { return mRaster->getGrainSize(); }
519  void setGrainSize(int grainSize) { mRaster->setGrainSize(grainSize); }
520 
525  typename IndxGridT::Ptr rasterizeSpheres(const ParticleListT& pa)
526  {
527  mRaster->rasterizeSpheres(pa);
528  return this->extract();
529  }
545  typename IndxGridT::Ptr rasterizeTrails(const ParticleListT& pa, Real delta=1.0)
546  {
547  mRaster->rasterizeTrails(pa, delta);
548  return this->extract();
549  }
550 
551 private:
552 
556  ParticlesToLevelSetAndId& operator=(const ParticlesToLevelSetAndId& rhs)
557  {
558  return *this;
559  }
561  typename IndxGridT::Ptr extract()
562  {
563  // Use topology copy constructors since output grids have the
564  // same topology as mDualGrid
565  const DualTreeT& dualTree = mDualGrid.tree();
566  typename IndxTreeT::Ptr indxTree(new IndxTreeT(dualTree,util::INVALID_IDX,TopologyCopy()));
567  typename IndxGridT::Ptr indxGrid = typename IndxGridT::Ptr(new IndxGridT(indxTree));
568  indxGrid->setTransform(mDualGrid.transformPtr());
569  typename RealTreeT::Ptr realTree(new RealTreeT(dualTree,mRealGrid.background(),TopologyCopy()));
570  mRealGrid.setTree(realTree);
571 
572  // Extract the level set and IDs from mDualGrid. We will
573  // explore the fact that by design active values always live
574  // at the leaf node level, i.e. no active tiles exist in level sets
575  typedef typename DualGridT::TreeType::LeafCIter LeafIterT;
576  typedef typename DualGridT::TreeType::LeafNodeType LeafT;
577  typedef typename LevelSetGridT::TreeType::LeafNodeType RealLeafT;
578  typedef typename IndxGridT::TreeType::LeafNodeType IndxLeafT;
579  RealTreeT& realTreeRef = *realTree;
580  IndxTreeT& indxTreeRef = *indxTree;
581  for (LeafIterT n = mDualGrid.tree().cbeginLeaf(); n; ++n) {
582  const LeafT& leaf = *n;
583  const Coord xyz = leaf.getOrigin();
584  // Get leafnodes that were allocated during topology contruction!
585  RealLeafT& i = *realTreeRef.probeLeaf(xyz);
586  IndxLeafT& j = *indxTreeRef.probeLeaf(xyz);
587  for (typename LeafT::ValueOnCIter m=leaf.cbeginValueOn(); m; ++m) {
588  // Use linear offset (vs coordinate) access for better performance!
589  const Index k = m.pos();
590  const DualT& v = *m;
591  i.setValueOnly(k, v.dist());
592  j.setValueOnly(k, v.id());
593  }
594  }
595  mRealGrid.signedFloodFill();//required since we only transferred active voxels!
596  return indxGrid;
597  }
598 
599  typedef ParticlesToLevelSet<DualGridT, ParticleListT, InterruptT, RealT> RasterT;
600  LevelSetGridT& mRealGrid;// input level set grid
601  DualGridT mDualGrid;// grid encoding both the level set and the point id
602  RasterT* mRaster;
603 };//end of ParticlesToLevelSetAndId
604 
605 } // namespace tools
606 } // namespace OPENVDB_VERSION_NAME
607 } // namespace openvdb
608 
609 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
610 
611 // Copyright (c) 2012-2013 DreamWorks Animation LLC
612 // All rights reserved. This software is distributed under the
613 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )