OpenVDB  0.104.0
ParticlesToLevelSet.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012 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 //
34 
35 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
37 
38 #include <tbb/parallel_reduce.h>
39 #include <tbb/blocked_range.h>
40 #include <boost/bind.hpp>
41 #include <boost/function.hpp>
42 #include <openvdb/util/Util.h>
43 #include <openvdb/Types.h>
44 #include <openvdb/Grid.h>
45 #include <openvdb/math/Math.h>
46 #include <openvdb/math/Transform.h>
47 #include <openvdb/util/NullInterrupter.h>
48 #include "Composite.h" // for csgUnion()
49 
50 
51 namespace openvdb {
53 namespace OPENVDB_VERSION_NAME {
54 namespace tools {
55 
103 
104 namespace local {
109 template <typename T>
110 struct DualTrait
111 {
112  static T merge(T dist, Index32) { return dist; }
113  static T split(T dist) { return dist; }
114 };
115 }// namespace local
116 
117 template<typename GridT,
118  typename ParticleListT,
119  typename InterruptT=util::NullInterrupter,
120  typename RealT = typename GridT::ValueType>
122 {
123 public:
138  ParticlesToLevelSet(GridT& grid, InterruptT* interrupt = NULL) :
139  mGrid(&grid),
140  mPa(NULL),
141  mDx(grid.transform().voxelSize()[0]),
142  mHalfWidth(local::DualTrait<ValueT>::split(grid.background()) / mDx),
143  mRmin(1.5),// corresponds to the Nyquist grid sampling frequency
144  mRmax(100.0f),// corresponds to a huge particle (probably too large!)
145  mGrainSize(1),
146  mInterrupter(interrupt),
147  mMinCount(0),
148  mMaxCount(0),
149  mIsSlave(false)
150  {
151  if ( !mGrid->hasUniformVoxels() ) {
153  "The transform must have uniform scale for ParticlesToLevelSet to function!");
154  }
155  if (mGrid->getGridClass() != GRID_LEVEL_SET) {
157  "ParticlesToLevelSet only supports level sets!"
158  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
159  }
161  //mGrid->newTree();
162  }
165  mGrid(new GridT(*other.mGrid, openvdb::ShallowCopy())),
166  mPa(other.mPa),
167  mDx(other.mDx),
168  mHalfWidth(other.mHalfWidth),
169  mRmin(other.mRmin),
170  mRmax(other.mRmax),
171  mGrainSize(other.mGrainSize),
172  mTask(other.mTask),
173  mInterrupter(other.mInterrupter),
174  mMinCount(0),
175  mMaxCount(0),
176  mIsSlave(true)
177  {
178  mGrid->newTree();
179  }
180  virtual ~ParticlesToLevelSet() { if (mIsSlave) delete mGrid; }
181 
183  RealT getHalfWidth() const { return mHalfWidth; }
184 
186  RealT getVoxelSize() const { return mDx; }
187 
189  RealT getRmin() const { return mRmin; }
191  RealT getRmax() const { return mRmax; }
192 
194  bool ignoredParticles() const { return mMinCount>0 || mMaxCount>0; }
196  size_t getMinCount() const { return mMinCount; }
198  size_t getMaxCount() const { return mMaxCount; }
199 
201  void setRmin(RealT Rmin) { mRmin = math::Max(RealT(0),Rmin); }
203  void setRmax(RealT Rmax) { mRmax = math::Max(mRmin,Rmax); }
204 
206  int getGrainSize() const { return mGrainSize; }
209  void setGrainSize(int grainSize) { mGrainSize = grainSize; }
210 
214  void rasterizeSpheres(const ParticleListT& pa)
215  {
216  mPa = &pa;
217  if (mInterrupter) mInterrupter->start("Rasterizing particles to level set using spheres");
218  mTask = boost::bind(&ParticlesToLevelSet::rasterSpheres, _1, _2);
219  this->cook();
220  if (mInterrupter) mInterrupter->end();
221  }
222 
238  void rasterizeTrails(const ParticleListT& pa, Real delta=1.0)
239  {
240  mPa = &pa;
241  if (mInterrupter) mInterrupter->start("Rasterizing particles to level set using trails");
242  mTask = boost::bind(&ParticlesToLevelSet::rasterTrails, _1, _2, RealT(delta));
243  this->cook();
244  if (mInterrupter) mInterrupter->end();
245  }
246 
247  // ========> DO NOT CALL ANY OF THE PUBLIC METHODS BELOW! <=============
248 
252  void operator()(const tbb::blocked_range<size_t>& r)
253  {
254  if (!mTask) {
255  OPENVDB_THROW(ValueError, "mTask is undefined - don't call operator() directly!");
256  }
257  mTask(this, r);
258  }
259 
263  void join(ParticlesToLevelSet& other)
264  {
265  tools::csgUnion(*mGrid, *other.mGrid, /*prune=*/true);
266  mMinCount += other.mMinCount;
267  mMaxCount += other.mMaxCount;
268  }
269 
270 private:
271  typedef typename GridT::ValueType ValueT;
272  typedef typename GridT::Accessor Accessor;
273  typedef typename boost::function<void (ParticlesToLevelSet*,
274  const tbb::blocked_range<size_t>&)> FuncType;
275  GridT* mGrid;
276  const ParticleListT* mPa;//list of particles
277  const RealT mDx;//size of voxel in world units
278  const RealT mHalfWidth;//half width of narrow band LS in voxels
279  RealT mRmin;//ignore particles smaller than this radius in voxels
280  RealT mRmax;//ignore particles larger than this radius in voxels
281  int mGrainSize;
282  FuncType mTask;
283  InterruptT* mInterrupter;
284  size_t mMinCount, mMaxCount;//counters for ignored particles!
285  const bool mIsSlave;
286 
287  void cook()
288  {
289  if (mGrainSize>0) {
290  tbb::parallel_reduce(tbb::blocked_range<size_t>(0,mPa->size(),mGrainSize), *this);
291  } else {
292  (*this)(tbb::blocked_range<size_t>(0, mPa->size()));
293  }
294  }
296  inline bool ignoreParticle(RealT R)
297  {
298  if (R<mRmin) {// below the cutoff radius
299  ++mMinCount;
300  return true;
301  }
302  if (R>mRmax) {// above the cutoff radius
303  ++mMaxCount;
304  return true;
305  }
306  return false;
307  }
323  inline bool rasterSphere(const Vec3R &P, RealT R, Index32 id, Accessor& accessor)
324  {
325  const ValueT inside = -mGrid->background();
326  const RealT dx = mDx;
327  const RealT max = R + mHalfWidth;// maximum distance in voxel units
328  const Coord a(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max));
329  const Coord b(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
330  const RealT max2 = math::Pow2(max);//square of maximum distance in voxel units
331  const RealT min2 = math::Pow2(math::Max(RealT(0), R - mHalfWidth));//square of minimum distance
332  ValueT v;
333  size_t count = 0;
334  for ( Coord c = a; c.x() <= b.x(); ++c.x() ) {
335  //only check interrupter every 32'th scan in x
336  if (!(count++ & (1<<5)-1) && util::wasInterrupted(mInterrupter)) {
337  tbb::task::self().cancel_group_execution();
338  return false;
339  }
340  RealT x2 = math::Pow2( c.x() - P[0] );
341  for ( c.y() = a.y(); c.y() <= b.y(); ++c.y() ) {
342  RealT x2y2 = x2 + math::Pow2( c.y() - P[1] );
343  for ( c.z() = a.z(); c.z() <= b.z(); ++c.z() ) {
344  RealT x2y2z2 = x2y2 + math::Pow2( c.z() - P[2] );//square distance from c to P
345  if ( x2y2z2 >= max2 || (!accessor.probeValue(c,v) && v<ValueT(0)) )
346  continue;//outside narrow band of the particle or inside existing ls
347  if ( x2y2z2 <= min2 ) {//inside narrowband of the particle.
348  accessor.setValueOff(c, inside);
349  continue;
350  }
351  // distance in world units
352  const ValueT d = local::DualTrait<ValueT>::merge(dx*(math::Sqrt(x2y2z2) - R), id);
353  if (d < v) accessor.setValue(c, d);//CSG union
354  }//end loop over z
355  }//end loop over y
356  }//end loop over x
357  return true;
358  }
359 
363  void rasterSpheres(const tbb::blocked_range<size_t> &r)
364  {
365  Accessor accessor = mGrid->getAccessor(); // local accessor
366  const RealT inv_dx = RealT(1)/mDx;
367  bool run = true;
368  for (Index32 id = r.begin(), e=r.end(); run && id != e; ++id) {
369  const RealT R = inv_dx*mPa->radius(id);// in voxel units
370  if (this->ignoreParticle(R)) continue;
371  const Vec3R P = mGrid->transform().worldToIndex(mPa->pos(id));
372  run = this->rasterSphere(P, R, id, accessor);
373  }//end loop over particles
374  }
375 
386  void rasterTrails(const tbb::blocked_range<size_t> &r,
387  RealT delta)//scale distance between instances (eg 1.0)
388  {
389  Accessor accessor = mGrid->getAccessor(); // local accessor
390  const RealT inv_dx = RealT(1)/mDx, Rmin = mRmin;
391  bool run = true;
392  for (Index32 id = r.begin(), e=r.end(); run && id != e; ++id) {
393  const RealT R0 = inv_dx*mPa->radius(id);
394  if (this->ignoreParticle(R0)) continue;
395  const Vec3R P0 = mGrid->transform().worldToIndex(mPa->pos(id)),
396  V = inv_dx*mPa->vel(id);
397  const RealT speed = V.length(), inv_speed=1.0/speed;
398  const Vec3R N = -V*inv_speed;// inverse normalized direction
399  Vec3R P = P0;// local position of instance
400  RealT R = R0, d=0;// local radius and length of trail
401  for (size_t m=0; run && d < speed ; ++m) {
402  run = this->rasterSphere(P, R, id, accessor);
403  P += 0.5*delta*R*N;// adaptive offset along inverse velocity direction
404  d = (P-P0).length();// current length of trail
405  R = R0-(R0-Rmin)*d*inv_speed;// R = R0 -> mRmin(e.g. 1.5)
406  }//loop over sphere instances
407  }//end loop over particles
408  }
409 };//end of ParticlesToLevelSet class
410 
412 namespace local {
413 // This is a simple type that combines a distance value and a particle
414 // id. It's required for attribute transfer which is defined in the
415 // ParticlesToLevelSetAndId class below.
416 template <typename RealT>
417 class Dual
418 {
419 public:
420  explicit Dual() : mId(util::INVALID_IDX) {}
421  explicit Dual(RealT d) : mDist(d), mId(util::INVALID_IDX) {}
422  explicit Dual(RealT d, Index32 id) : mDist(d), mId(id) {}
423  Dual& operator=(const Dual& rhs) { mDist = rhs.mDist; mId = rhs.mId; return *this;}
424  bool isIdValid() const { return mId != util::INVALID_IDX; }
425  Index32 id() const { return mId; }
426  RealT dist() const { return 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  bool operator<=(const Dual& rhs) const { return mDist <= rhs.mDist; };
431  bool operator> (const Dual& rhs) const { return mDist > rhs.mDist; };
432  Dual operator+ (const Dual& rhs) const { return Dual(mDist+rhs.mDist); };
433  Dual operator+ (const RealT& rhs) const { return Dual(mDist+rhs); };
434  Dual operator- (const Dual& rhs) const { return Dual(mDist-rhs.mDist); };
435  Dual operator-() const { return Dual(-mDist); }
436 protected:
437  RealT mDist;
439 };
440 // Required by several of the tree nodes
441 template <typename RealT>
442 inline std::ostream& operator<<(std::ostream& ostr, const Dual<RealT>& rhs)
443 {
444  ostr << rhs.dist();
445  return ostr;
446 }
447 // Required by math::Abs
448 template <typename RealT>
449 inline Dual<RealT> Abs(const Dual<RealT>& x)
450 {
451  return Dual<RealT>(math::Abs(x.dist()),x.id());
452 }
453 // Specialization of trait class used to merge and split a distance and a particle id
454 template <typename T>
455 struct DualTrait<Dual<T> >
456 {
457  static Dual<T> merge(T dist, Index32 id) { return Dual<T>(dist, id); }
458  static T split(Dual<T> dual) { return dual.dist(); }
459 };
460 }// local namespace
462 
469 template<typename LevelSetGridT,
470  typename ParticleListT,
471  typename InterruptT = util::NullInterrupter>
473 {
474 public:
475  typedef typename LevelSetGridT::ValueType RealT;
476  typedef typename local::Dual<RealT> DualT;
477  typedef typename LevelSetGridT::TreeType RealTreeT;
478  typedef typename RealTreeT::template ValueConverter<DualT>::Type DualTreeT;
482 
483  ParticlesToLevelSetAndId(LevelSetGridT& ls, InterruptT* interrupter = NULL) :
484  mRealGrid(ls),
485  mDualGrid(DualT(ls.background()))
486  {
487  mDualGrid.setGridClass(ls.getGridClass());
488  mDualGrid.setTransform(ls.transformPtr());
489  mRaster = new RasterT(mDualGrid, interrupter);
490  }
491 
492  virtual ~ParticlesToLevelSetAndId() { delete mRaster; }
493 
495  RealT getHalfWidth() const { return mRaster->getHalfWidth(); }
496 
498  RealT getVoxelSize() const { return mRaster->getVoxelSize(); }
499 
501  bool ignoredParticles() const { return mRaster->ignoredParticles(); }
503  size_t getMinCount() const { return mRaster->getMinCount(); }
505  size_t getMaxCount() const { return mRaster->getMaxCount(); }
506 
508  RealT getRmin() const { return mRaster->getRmin(); }
510  RealT getRmax() const { return mRaster->getRmax(); }
511 
513  void setRmin(RealT Rmin) { mRaster->setRmin(Rmin); }
515  void setRmax(RealT Rmax) { mRaster->setRmax(Rmax); }
516 
518  int getGrainSize() const { return mRaster->getGrainSize(); }
521  void setGrainSize(int grainSize) { mRaster->setGrainSize(grainSize); }
522 
527  typename IndxGridT::Ptr rasterizeSpheres(const ParticleListT& pa)
528  {
529  mRaster->rasterizeSpheres(pa);
530  return this->extract();
531  }
547  typename IndxGridT::Ptr rasterizeTrails(const ParticleListT& pa, Real delta=1.0)
548  {
549  mRaster->rasterizeTrails(pa, delta);
550  return this->extract();
551  }
552 
553 private:
554 
558  ParticlesToLevelSetAndId& operator=(const ParticlesToLevelSetAndId& rhs)
559  {
560  return *this;
561  }
563  typename IndxGridT::Ptr extract()
564  {
565  // Use topology copy constructors since output grids have the
566  // same topology as mDualGrid
567  const DualTreeT& dualTree = mDualGrid.tree();
568  typename IndxTreeT::Ptr indxTree(new IndxTreeT(dualTree,util::INVALID_IDX,TopologyCopy()));
569  typename IndxGridT::Ptr indxGrid = typename IndxGridT::Ptr(new IndxGridT(indxTree));
570  indxGrid->setTransform(mDualGrid.transformPtr());
571  typename RealTreeT::Ptr realTree(new RealTreeT(dualTree,mRealGrid.background(),TopologyCopy()));
572  mRealGrid.setTree(realTree);
573 
574  // Extract the level set and IDs from mDualGrid. We will
575  // explore the fact that by design active values always live
576  // at the leaf node level, i.e. no active tiles exist in level sets
577  typedef typename DualGridT::TreeType::LeafCIter LeafIterT;
578  typedef typename DualGridT::TreeType::LeafNodeType LeafT;
579  typedef typename LevelSetGridT::TreeType::LeafNodeType RealLeafT;
580  typedef typename IndxGridT::TreeType::LeafNodeType IndxLeafT;
581  RealTreeT& realTreeRef = *realTree;
582  IndxTreeT& indxTreeRef = *indxTree;
583  for (LeafIterT n = mDualGrid.tree().cbeginLeaf(); n; ++n) {
584  const LeafT& leaf = *n;
585  const Coord xyz = leaf.getOrigin();
586  // Get leafnodes that were allocated during topology contruction!
587  RealLeafT& i = *realTreeRef.probeLeaf(xyz);
588  IndxLeafT& j = *indxTreeRef.probeLeaf(xyz);
589  for (typename LeafT::ValueOnCIter m=leaf.cbeginValueOn(); m; ++m) {
590  // Use linear offset (vs coordinate) access for better performance!
591  const Index k = m.pos();
592  const DualT& v = *m;
593  i.setValueOnly(k, v.dist());
594  j.setValueOnly(k, v.id());
595  }
596  }
597  mRealGrid.signedFloodFill();//required since we only transferred active voxels!
598  return indxGrid;
599  }
600 
601  typedef ParticlesToLevelSet<DualGridT, ParticleListT, InterruptT, RealT> RasterT;
602  LevelSetGridT& mRealGrid;// input level set grid
603  DualGridT mDualGrid;// grid encoding both the level set and the point id
604  RasterT* mRaster;
605 };//end of ParticlesToLevelSetAndId
606 
607 } // namespace tools
608 } // namespace OPENVDB_VERSION_NAME
609 } // namespace openvdb
610 
611 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
612 
613 // Copyright (c) 2012 DreamWorks Animation LLC
614 // All rights reserved. This software is distributed under the
615 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )