HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ParticlesToLevelSet.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2017 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 //
29 ///////////////////////////////////////////////////////////////////////////
30 
31 /// @author Ken Museth
32 ///
33 /// @file ParticlesToLevelSet.h
34 ///
35 /// @brief This tool converts particles (with position, radius and velocity)
36 /// into a signed distance field encoded as a narrow band level set.
37 /// Optionally, arbitrary attributes on the particles can be transferred
38 /// resulting in an additional attribute grid with the same topology
39 /// as the level set grid.
40 ///
41 /// @note This fast particle to level set converter is always intended
42 /// to be combined with some kind of surface post processing,
43 /// i.e. tools::Filter. Without such post processing the generated
44 /// surface is typically too noisy and blobby. However it serves as a
45 /// great and fast starting point for subsequent level set surface
46 /// processing and convolution.
47 ///
48 /// The @c ParticleListT template argument below refers to any class
49 /// with the following interface (see unittest/TestParticlesToLevelSet.cc
50 /// and SOP_DW_OpenVDBParticleVoxelizer for practical examples):
51 /// @code
52 ///
53 /// class ParticleList {
54 /// ...
55 /// public:
56 /// typedef openvdb::Vec3R PosType;
57 ///
58 /// // Return the total number of particles in list.
59 /// // Always required!
60 /// size_t size() const;
61 ///
62 /// // Get the world space position of the nth particle.
63 /// // Required by ParticledToLevelSet::rasterizeSphere(*this,radius).
64 /// void getPos(size_t n, Vec3R& xyz) const;
65 ///
66 /// // Get the world space position and radius of the nth particle.
67 /// // Required by ParticledToLevelSet::rasterizeSphere(*this).
68 /// void getPosRad(size_t n, Vec3R& xyz, Real& rad) const;
69 ///
70 /// // Get the world space position, radius and velocity of the nth particle.
71 /// // Required by ParticledToLevelSet::rasterizeSphere(*this,radius).
72 /// void getPosRadVel(size_t n, Vec3R& xyz, Real& rad, Vec3R& vel) const;
73 ///
74 /// // Get the attribute of the nth particle. AttributeType is user-defined!
75 /// // Only required if attribute transfer is enabled in ParticlesToLevelSet.
76 /// void getAtt(size_t n, AttributeType& att) const;
77 /// };
78 /// @endcode
79 ///
80 /// @note See unittest/TestParticlesToLevelSet.cc for an example.
81 ///
82 /// The @c InterruptT template argument below refers to any class
83 /// with the following interface:
84 /// @code
85 /// class Interrupter {
86 /// ...
87 /// public:
88 /// void start(const char* name = nullptr) // called when computations begin
89 /// void end() // called when computations end
90 /// bool wasInterrupted(int percent=-1) // return true to break computation
91 /// };
92 /// @endcode
93 ///
94 /// @note If no template argument is provided for this InterruptT
95 /// the util::NullInterrupter is used which implies that all
96 /// interrupter calls are no-ops (i.e. incurs no computational overhead).
97 
98 #ifndef OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
99 #define OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
100 
101 #include <tbb/parallel_reduce.h>
102 #include <tbb/blocked_range.h>
103 #include <hboost/bind.hpp>
104 #include <hboost/function.hpp>
105 #include <hboost/type_traits/is_floating_point.hpp>
106 #include <hboost/utility/enable_if.hpp>
107 #include <hboost/mpl/if.hpp>
108 #include <openvdb/Types.h>
109 #include <openvdb/Grid.h>
110 #include <openvdb/math/Math.h>
111 #include <openvdb/math/Transform.h>
113 #include "Composite.h" // for csgUnion()
114 #include "PointPartitioner.h"
115 #include "Prune.h"
116 #include "SignedFloodFill.h"
117 
118 namespace openvdb {
120 namespace OPENVDB_VERSION_NAME {
121 namespace tools {
122 
123 namespace p2ls_internal {
124 // This is a simple type that combines a distance value and a particle
125 // attribute. It's required for attribute transfer which is performed
126 // in the ParticlesToLevelSet::Raster member class defined below.
127 template<typename VisibleT, typename BlindT> class BlindData;
128 }// namespace p2ls_internal
129 
130 
131 template<typename SdfGridT,
132  typename AttributeT = void,
133  typename InterrupterT = util::NullInterrupter>
135 {
136 public:
138  typedef InterrupterT InterrupterType;
139 
140  typedef SdfGridT SdfGridType;
141  typedef typename SdfGridT::ValueType SdfType;
142 
144  typedef typename SdfGridT::template ValueConverter<AttType>::Type AttGridType;
145 
147 
148  /// @brief Constructor using an exiting signed distance,
149  /// i.e. narrow band level set, grid.
150  ///
151  /// @param grid Level set grid in which particles are rasterized
152  /// @param interrupt Callback to interrupt a long-running process
153  ///
154  /// @note The input grid is assumed to be a valid level set and if
155  /// it already contains voxels (with SDF values) particles are unioned
156  /// onto the existing level set surface. However, if attribute transfer
157  /// is enabled, i.e. AttributeT != void, attributes are only
158  /// generated for voxels that overlap with particles, not the existing
159  /// voxels in the input grid (for which no attributes exist!).
160  ///
161  /// @details The width in voxel units of the generated narrow band level set is
162  /// given by 2*background/dx, where background is the background value
163  /// stored in the grid, and dx is the voxel size derived from the
164  /// transform also stored in the grid. Also note that -background
165  /// corresponds to the constant value inside the generated narrow
166  /// band level sets. Finally the default NullInterrupter should
167  /// compile out interruption checks during optimization, thus
168  /// incurring no run-time overhead.
169  explicit ParticlesToLevelSet(SdfGridT& grid, InterrupterT* interrupt = nullptr);
170 
171  /// Destructor
172  ~ParticlesToLevelSet() { delete mBlindGrid; }
173 
174  /// @brief This methods syncs up the level set and attribute grids
175  /// and therefore needs to be called before any of these grids are
176  /// used and after the last call to any of the rasterizer methods.
177  ///
178  /// @note Avoid calling this method more than once and only after
179  /// all the particles have been rasterized. It has no effect or
180  /// overhead if attribute transfer is disabled, i.e. AttributeT =
181  /// void and prune is false.
182  void finalize(bool prune = false);
183 
184  /// @brief Return a shared pointer to the grid containing the
185  /// (optional) attribute.
186  ///
187  /// @warning If attribute transfer was disabled, i.e. AttributeT =
188  /// void, or finalize() was not called the pointer is nullptr!
189  typename AttGridType::Ptr attributeGrid() { return mAttGrid; }
190 
191  /// @brief Return the size of a voxel in world units
192  Real getVoxelSize() const { return mDx; }
193 
194  /// @brief Return the half-width of the narrow band in voxel units
195  Real getHalfWidth() const { return mHalfWidth; }
196 
197  /// @brief Return the smallest radius allowed in voxel units
198  Real getRmin() const { return mRmin; }
199  /// @brief Return the largest radius allowed in voxel units
200  Real getRmax() const { return mRmax; }
201 
202  /// @brief Return true if any particles were ignored due to their size
203  bool ignoredParticles() const { return mMinCount>0 || mMaxCount>0; }
204  /// @brief Return number of small particles that were ignore due to Rmin
205  size_t getMinCount() const { return mMinCount; }
206  /// @brief Return number of large particles that were ignore due to Rmax
207  size_t getMaxCount() const { return mMaxCount; }
208 
209  /// @brief set the smallest radius allowed in voxel units
210  void setRmin(Real Rmin) { mRmin = math::Max(Real(0),Rmin); }
211  /// @brief set the largest radius allowed in voxel units
212  void setRmax(Real Rmax) { mRmax = math::Max(mRmin,Rmax); }
213 
214  /// @brief Returns the grain-size used for multi-threading
215  int getGrainSize() const { return mGrainSize; }
216  /// @brief Set the grain-size used for multi-threading.
217  /// @note A grainsize of 0 or less disables multi-threading!
218  void setGrainSize(int grainSize) { mGrainSize = grainSize; }
219 
220  /// @brief Rasterize a sphere per particle derived from their
221  /// position and radius. All spheres are CSG unioned.
222  ///
223  /// @param pa Particles with position and radius.
224  template <typename ParticleListT>
225  void rasterizeSpheres(const ParticleListT& pa);
226 
227  /// @brief Rasterize a sphere per particle derived from their
228  /// position and constant radius. All spheres are CSG unioned.
229  ///
230  /// @param pa Particles with position.
231  /// @param radius Constant particle radius in world units.
232  template <typename ParticleListT>
233  void rasterizeSpheres(const ParticleListT& pa, Real radius);
234 
235  /// @brief Rasterize a trail per particle derived from their
236  /// position, radius and velocity. Each trail is generated
237  /// as CSG unions of sphere instances with decreasing radius.
238  ///
239  /// @param pa particles with position, radius and velocity.
240  /// @param delta controls distance between sphere instances
241  /// (default=1). Be careful not to use too small values since this
242  /// can lead to excessive computation per trail (which the
243  /// interrupter can't stop).
244  ///
245  /// @note The direction of a trail is inverse to the direction of
246  /// the velocity vector, and the length is given by |V|. The radius
247  /// at the head of the trail is given by the radius of the particle
248  /// and the radius at the tail of the trail is Rmin voxel units which
249  /// has a default value of 1.5 corresponding to the Nyquist
250  /// frequency!
251  template <typename ParticleListT>
252  void rasterizeTrails(const ParticleListT& pa, Real delta=1.0);
253 
254 private:
256  typedef typename SdfGridT::template ValueConverter<BlindType>::Type BlindGridType;
257 
258  /// Class with multi-threaded implementation of particle rasterization
259  template<typename ParticleListT, typename GridT> struct Raster;
260 
261  SdfGridType* mSdfGrid;
262  typename AttGridType::Ptr mAttGrid;
263  BlindGridType* mBlindGrid;
264  InterrupterT* mInterrupter;
265  Real mDx, mHalfWidth;
266  Real mRmin, mRmax;//ignore particles outside this range of radii in voxel
267  size_t mMinCount, mMaxCount;//counters for ignored particles!
268  int mGrainSize;
269 
270 };//end of ParticlesToLevelSet class
271 
272 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
274 ParticlesToLevelSet(SdfGridT& grid, InterrupterT* interrupter) :
275  mSdfGrid(&grid),
276  mBlindGrid(nullptr),
277  mInterrupter(interrupter),
278  mDx(grid.voxelSize()[0]),
279  mHalfWidth(grid.background()/mDx),
280  mRmin(1.5),// corresponds to the Nyquist grid sampling frequency
281  mRmax(100.0),// corresponds to a huge particle (probably too large!)
282  mMinCount(0),
283  mMaxCount(0),
284  mGrainSize(1)
285 {
286  if (!mSdfGrid->hasUniformVoxels() ) {
287  OPENVDB_THROW(RuntimeError,
288  "ParticlesToLevelSet only supports uniform voxels!");
289  }
290  if (mSdfGrid->getGridClass() != GRID_LEVEL_SET) {
291  OPENVDB_THROW(RuntimeError,
292  "ParticlesToLevelSet only supports level sets!"
293  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
294  }
295 
296  if (!DisableT::value) {
297  mBlindGrid = new BlindGridType(BlindType(grid.background()));
298  mBlindGrid->setTransform(mSdfGrid->transform().copy());
299  }
300 }
301 
302 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
303 template <typename ParticleListT>
305 rasterizeSpheres(const ParticleListT& pa)
306 {
307  if (DisableT::value) {
308  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
309  r.rasterizeSpheres();
310  } else {
311  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
312  r.rasterizeSpheres();
313  }
314 }
315 
316 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
317 template <typename ParticleListT>
319 rasterizeSpheres(const ParticleListT& pa, Real radius)
320 {
321  if (DisableT::value) {
322  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
323  r.rasterizeSpheres(radius/mDx);
324  } else {
325  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
326  r.rasterizeSpheres(radius/mDx);
327  }
328 }
329 
330 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
331 template <typename ParticleListT>
333 rasterizeTrails(const ParticleListT& pa, Real delta)
334 {
335  if (DisableT::value) {
336  Raster<ParticleListT, SdfGridT> r(*this, mSdfGrid, pa);
337  r.rasterizeTrails(delta);
338  } else {
339  Raster<ParticleListT, BlindGridType> r(*this, mBlindGrid, pa);
340  r.rasterizeTrails(delta);
341  }
342 }
343 
344 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
345 inline void
347 {
348  if (mBlindGrid == nullptr) {
349  if (prune) tools::pruneLevelSet(mSdfGrid->tree());
350  return;
351  } else {
352  if (prune) tools::prune(mBlindGrid->tree());
353  }
354 
355  typedef typename SdfGridType::TreeType SdfTreeT;
356  typedef typename AttGridType::TreeType AttTreeT;
357  typedef typename BlindGridType::TreeType BlindTreeT;
358  // Use topology copy constructors since output grids have the same topology as mBlindDataGrid
359  const BlindTreeT& tree = mBlindGrid->tree();
360 
361  // New level set tree
362  typename SdfTreeT::Ptr sdfTree(new SdfTreeT(
363  tree, tree.background().visible(), openvdb::TopologyCopy()));
364 
365  // Note this overwrites any existing attribute grids!
366  typename AttTreeT::Ptr attTree(new AttTreeT(
367  tree, tree.background().blind(), openvdb::TopologyCopy()));
368  mAttGrid = typename AttGridType::Ptr(new AttGridType(attTree));
369  mAttGrid->setTransform(mBlindGrid->transform().copy());
370 
371  // Extract the level set and IDs from mBlindDataGrid. We will
372  // explore the fact that by design active values always live
373  // at the leaf node level, i.e. level sets have no active tiles!
374  typedef typename BlindTreeT::LeafCIter LeafIterT;
375  typedef typename BlindTreeT::LeafNodeType LeafT;
376  typedef typename SdfTreeT::LeafNodeType SdfLeafT;
377  typedef typename AttTreeT::LeafNodeType AttLeafT;
378  for (LeafIterT n = tree.cbeginLeaf(); n; ++n) {
379  const LeafT& leaf = *n;
380  const openvdb::Coord xyz = leaf.origin();
381  // Get leafnodes that were allocated during topology construction!
382  SdfLeafT* sdfLeaf = sdfTree->probeLeaf(xyz);
383  AttLeafT* attLeaf = attTree->probeLeaf(xyz);
384  // Use linear offset (vs coordinate) access for better performance!
385  typename LeafT::ValueOnCIter m=leaf.cbeginValueOn();
386  if (!m) {//no active values in leaf node so copy everything
387  for (openvdb::Index k = 0; k!=LeafT::SIZE; ++k) {
388  const BlindType& v = leaf.getValue(k);
389  sdfLeaf->setValueOnly(k, v.visible());
390  attLeaf->setValueOnly(k, v.blind());
391  }
392  } else {//only copy active values (using flood fill for the inactive values)
393  for(; m; ++m) {
394  const openvdb::Index k = m.pos();
395  const BlindType& v = *m;
396  sdfLeaf->setValueOnly(k, v.visible());
397  attLeaf->setValueOnly(k, v.blind());
398  }
399  }
400  }
401 
402  tools::signedFloodFill(*sdfTree);//required since we only transferred active voxels!
403 
404  if (mSdfGrid->empty()) {
405  mSdfGrid->setTree(sdfTree);
406  } else {
407  tools::csgUnion(mSdfGrid->tree(), *sdfTree, /*prune=*/true);
408  }
409 }
410 
411 ///////////////////////////////////////////////////////////
412 
413 template<typename SdfGridT, typename AttributeT, typename InterrupterT>
414 template<typename ParticleListT, typename GridT>
415 struct ParticlesToLevelSet<SdfGridT, AttributeT, InterrupterT>::Raster
416 {
418  typedef ParticlesToLevelSet<SdfGridT, AttributeT, InterrupterT> ParticlesToLevelSetT;
419  typedef typename ParticlesToLevelSetT::SdfType SdfT;//type of signed distance values
420  typedef typename ParticlesToLevelSetT::AttType AttT;//type of particle attribute
421  typedef typename GridT::ValueType ValueT;
422  typedef typename GridT::Accessor AccessorT;
423  typedef typename GridT::TreeType TreeT;
424  typedef typename TreeT::LeafNodeType LeafNodeT;
425  typedef PointPartitioner<Index32, LeafNodeT::LOG2DIM> PointPartitionerT;
426 
427 
428  /// @brief Main constructor
429  Raster(ParticlesToLevelSetT& parent, GridT* grid, const ParticleListT& particles)
430  : mParent(parent)
431  , mParticles(particles)
432  , mGrid(grid)
433  , mMap(*(mGrid->transform().baseMap()))
434  , mMinCount(0)
435  , mMaxCount(0)
436  , mIsCopy(false)
437  {
438  mPointPartitioner = new PointPartitionerT();
439  mPointPartitioner->construct(particles, mGrid->transform());
440  }
441 
442  /// @brief Copy constructor called by tbb threads
443  Raster(Raster& other, tbb::split)
444  : mParent(other.mParent)
445  , mParticles(other.mParticles)
446  , mGrid(new GridT(*other.mGrid, openvdb::ShallowCopy()))
447  , mMap(other.mMap)
448  , mMinCount(0)
449  , mMaxCount(0)
450  , mTask(other.mTask)
451  , mIsCopy(true)
452  , mPointPartitioner(other.mPointPartitioner)
453  {
454  mGrid->newTree();
455  }
456 
457  virtual ~Raster() {
458 
459  // Copies construct temporary grids that have to be deleted
460  // but the original has ownership of the bucket array
461  if (mIsCopy) {
462  delete mGrid;
463  } else {
464  delete mPointPartitioner;
465  }
466  }
467 
468  /// @brief Rasterize a sphere per particle derived from their
469  /// position and radius. All spheres are CSG unioned.
470  void rasterizeSpheres()
471  {
472  mMinCount = mMaxCount = 0;
473  if (mParent.mInterrupter) {
474  mParent.mInterrupter->start("Rasterizing particles to level set using spheres");
475  }
476  mTask = hboost::bind(&Raster::rasterSpheres, hboost::placeholders::_1, hboost::placeholders::_2);
477  this->cook();
478  if (mParent.mInterrupter) mParent.mInterrupter->end();
479  }
480  /// @brief Rasterize a sphere per particle derived from their
481  /// position and constant radius. All spheres are CSG unioned.
482  /// @param radius constant radius of all particles in voxel units.
483  void rasterizeSpheres(Real radius)
484  {
485  mMinCount = radius < mParent.mRmin ? mParticles.size() : 0;
486  mMaxCount = radius > mParent.mRmax ? mParticles.size() : 0;
487  if (mMinCount>0 || mMaxCount>0) {//skipping all particles!
488  mParent.mMinCount = mMinCount;
489  mParent.mMaxCount = mMaxCount;
490  } else {
491  if (mParent.mInterrupter) {
492  mParent.mInterrupter->start(
493  "Rasterizing particles to level set using const spheres");
494  }
495  mTask = hboost::bind(&Raster::rasterFixedSpheres, hboost::placeholders::_1, hboost::placeholders::_2, SdfT(radius));
496  this->cook();
497  if (mParent.mInterrupter) mParent.mInterrupter->end();
498  }
499  }
500  /// @brief Rasterize a trail per particle derived from their
501  /// position, radius and velocity. Each trail is generated
502  /// as CSG unions of sphere instances with decreasing radius.
503  ///
504  /// @param delta controls distance between sphere instances
505  /// (default=1). Be careful not to use too small values since this
506  /// can lead to excessive computation per trail (which the
507  /// interrupter can't stop).
508  ///
509  /// @note The direction of a trail is inverse to the direction of
510  /// the velocity vector, and the length is given by |V|. The radius
511  /// at the head of the trail is given by the radius of the particle
512  /// and the radius at the tail of the trail is Rmin voxel units which
513  /// has a default value of 1.5 corresponding to the Nyquist frequency!
514  void rasterizeTrails(Real delta=1.0)
515  {
516  mMinCount = mMaxCount = 0;
517  if (mParent.mInterrupter) {
518  mParent.mInterrupter->start("Rasterizing particles to level set using trails");
519  }
520  mTask = hboost::bind(&Raster::rasterTrails, hboost::placeholders::_1, hboost::placeholders::_2, SdfT(delta));
521  this->cook();
522  if (mParent.mInterrupter) mParent.mInterrupter->end();
523  }
524 
525  /// @brief Kicks off the optionally multithreaded computation
526  void operator()(const tbb::blocked_range<size_t>& r)
527  {
528  assert(mTask);
529  mTask(this, r);
530  mParent.mMinCount = mMinCount;
531  mParent.mMaxCount = mMaxCount;
532  }
533 
534  /// @brief Reguired by tbb::parallel_reduce
535  void join(Raster& other)
536  {
537  tools::csgUnion(*mGrid, *other.mGrid, /*prune=*/true);
538  mMinCount += other.mMinCount;
539  mMaxCount += other.mMaxCount;
540  }
541 private:
542  /// Disallow assignment since some of the members are references
543  Raster& operator=(const Raster&) { return *this; }
544 
545  /// @return true if the particle is too small or too large
546  bool ignoreParticle(SdfT R)
547  {
548  if (R < mParent.mRmin) {// below the cutoff radius
549  ++mMinCount;
550  return true;
551  }
552  if (R > mParent.mRmax) {// above the cutoff radius
553  ++mMaxCount;
554  return true;
555  }
556  return false;
557  }
558  /// @brief Reguired by tbb::parallel_reduce to multithreaded
559  /// rasterization of particles as spheres with variable radius
560  ///
561  /// @param r tbb's default range referring to the list of particles
562  void rasterSpheres(const tbb::blocked_range<size_t>& r)
563  {
564  AccessorT acc = mGrid->getAccessor(); // local accessor
565  bool run = true;
566  const SdfT invDx = SdfT(1/mParent.mDx);
567  AttT att;
568  Vec3R pos;
569  Real rad;
570 
571  // Loop over buckets
572  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
573  // Loop over particles in bucket n.
574  typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
575  for ( ; run && iter; ++iter) {
576  const Index32& id = *iter;
577  mParticles.getPosRad(id, pos, rad);
578  const SdfT R = SdfT(invDx * rad);// in voxel units
579  if (this->ignoreParticle(R)) continue;
580  const Vec3R P = mMap.applyInverseMap(pos);
581  this->getAtt<DisableT>(id, att);
582  run = this->makeSphere(P, R, att, acc);
583  }//end loop over particles
584  }//end loop over buckets
585  }
586 
587  /// @brief Reguired by tbb::parallel_reduce to multithreaded
588  /// rasterization of particles as spheres with a fixed radius
589  ///
590  /// @param r tbb's default range referring to the list of particles
591  void rasterFixedSpheres(const tbb::blocked_range<size_t>& r, SdfT R)
592  {
593  const SdfT
594  dx = static_cast<SdfT>(mParent.mDx),
595  w = static_cast<SdfT>(mParent.mHalfWidth); // in voxel units
596  AccessorT acc = mGrid->getAccessor(); // local accessor
597  const ValueT inside = -mGrid->background();
598  const SdfT max = R + w;// maximum distance in voxel units
599  const SdfT max2 = math::Pow2(max);//square of maximum distance in voxel units
600  const SdfT min2 = math::Pow2(math::Max(SdfT(0), R - w));//square of minimum distance
601  ValueT v;
602  size_t count = 0;
603  AttT att;
604  Vec3R pos;
605 
606  // Loop over buckets
607  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
608  // Loop over particles in bucket n.
609  typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
610  for ( ; iter; ++iter) {
611  const Index32& id = *iter;
612  this->getAtt<DisableT>(id, att);
613  mParticles.getPos(id, pos);
614  const Vec3R P = mMap.applyInverseMap(pos);
615  const Coord a(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max));
616  const Coord b(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
617  for (Coord c = a; c.x() <= b.x(); ++c.x()) {
618  //only check interrupter every 32'th scan in x
619  if (!(count++ & ((1<<5)-1)) && util::wasInterrupted(mParent.mInterrupter)) {
620  tbb::task::self().cancel_group_execution();
621  return;
622  }
623  SdfT x2 = static_cast<SdfT>(math::Pow2(c.x() - P[0]));
624  for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
625  SdfT x2y2 = static_cast<SdfT>(x2 + math::Pow2(c.y() - P[1]));
626  for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
627  SdfT x2y2z2 = static_cast<SdfT>(
628  x2y2 + math::Pow2(c.z()- P[2])); // square distance from c to P
629  if (x2y2z2 >= max2 || (!acc.probeValue(c,v) && v<ValueT(0)))
630  continue;//outside narrow band of particle or inside existing level set
631  if (x2y2z2 <= min2) {//inside narrow band of the particle.
632  acc.setValueOff(c, inside);
633  continue;
634  }
635  // convert signed distance from voxel units to world units
636  const ValueT d=Merge(dx*(math::Sqrt(x2y2z2) - R), att);
637  if (d < v) acc.setValue(c, d);//CSG union
638  }//end loop over z
639  }//end loop over y
640  }//end loop over x
641  }//end loop over particles
642  }// end loop over buckts
643  }
644 
645  /// @brief Reguired by tbb::parallel_reduce to multithreaded
646  /// rasterization of particles as spheres with velocity blurring
647  ///
648  /// @param r tbb's default range referring to the list of particles
649  void rasterTrails(const tbb::blocked_range<size_t>& r, SdfT delta)
650  {
651  AccessorT acc = mGrid->getAccessor(); // local accessor
652  bool run = true;
653  AttT att;
654  Vec3R pos, vel;
655  Real rad;
656  const Vec3R origin = mMap.applyInverseMap(Vec3R(0,0,0));
657  const SdfT Rmin = SdfT(mParent.mRmin), invDx = SdfT(1/mParent.mDx);
658 
659  // Loop over buckets
660  for (size_t n = r.begin(), N = r.end(); n < N; ++n) {
661  // Loop over particles in bucket n.
662  typename PointPartitionerT::IndexIterator iter = mPointPartitioner->indices(n);
663  for ( ; run && iter; ++iter) {
664  const Index32& id = *iter;
665  mParticles.getPosRadVel(id, pos, rad, vel);
666  const SdfT R0 = SdfT(invDx*rad);
667  if (this->ignoreParticle(R0)) continue;
668  this->getAtt<DisableT>(id, att);
669  const Vec3R P0 = mMap.applyInverseMap(pos);
670  const Vec3R V = mMap.applyInverseMap(vel) - origin;//exclude translation
671  const SdfT speed = SdfT(V.length()), inv_speed = SdfT(1.0/speed);
672  const Vec3R Nrml = -V*inv_speed;// inverse normalized direction
673  Vec3R P = P0;// local position of instance
674  SdfT R = R0, d=0;// local radius and length of trail
675  for (size_t m=0; run && d <= speed ; ++m) {
676  run = this->makeSphere(P, R, att, acc);
677  P += 0.5*delta*R*Nrml;// adaptive offset along inverse velocity direction
678  d = SdfT((P-P0).length());// current length of trail
679  R = R0-(R0-Rmin)*d*inv_speed;// R = R0 -> mRmin(e.g. 1.5)
680  }//end loop over sphere instances
681  }//end loop over particles
682  }//end loop over buckets
683  }
684 
685  void cook()
686  {
687  // parallelize over the point buckets
688  const Index32 bucketCount = Index32(mPointPartitioner->size());
689 
690  if (mParent.mGrainSize>0) {
691  tbb::parallel_reduce(
692  tbb::blocked_range<size_t>(0, bucketCount, mParent.mGrainSize), *this);
693  } else {
694  (*this)(tbb::blocked_range<size_t>(0, bucketCount));
695  }
696  }
697 
698  /// @brief Rasterize sphere at position P and radius R into a
699  /// narrow-band level set with half-width, mHalfWidth.
700  /// @return false if it was interrupted
701  ///
702  /// @param P coordinates of the particle position in voxel units
703  /// @param R radius of particle in voxel units
704  /// @param id
705  /// @param accessor grid accessor with a private copy of the grid
706  ///
707  /// @note For best performance all computations are performed in
708  /// voxel-space with the important exception of the final level set
709  /// value that is converted to world units (e.g. the grid stores
710  /// the closest Euclidean signed distances measured in world
711  /// units). Also note we use the convention of positive distances
712  /// outside the surface and negative distances inside the surface.
713  bool makeSphere(const Vec3R &P, SdfT R, const AttT& att, AccessorT& acc)
714  {
715  const ValueT inside = -mGrid->background();
716  const SdfT dx = SdfT(mParent.mDx), w = SdfT(mParent.mHalfWidth);
717  const SdfT max = R + w;// maximum distance in voxel units
718  const Coord a(math::Floor(P[0]-max),math::Floor(P[1]-max),math::Floor(P[2]-max));
719  const Coord b(math::Ceil( P[0]+max),math::Ceil( P[1]+max),math::Ceil( P[2]+max));
720  const SdfT max2 = math::Pow2(max);//square of maximum distance in voxel units
721  const SdfT min2 = math::Pow2(math::Max(SdfT(0), R - w));//square of minimum distance
722  ValueT v;
723  size_t count = 0;
724  for ( Coord c = a; c.x() <= b.x(); ++c.x() ) {
725  //only check interrupter every 32'th scan in x
726  if (!(count++ & ((1<<5)-1)) && util::wasInterrupted(mParent.mInterrupter)) {
727  tbb::task::self().cancel_group_execution();
728  return false;
729  }
730  SdfT x2 = SdfT(math::Pow2(c.x() - P[0]));
731  for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
732  SdfT x2y2 = SdfT(x2 + math::Pow2(c.y() - P[1]));
733  for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
734  SdfT x2y2z2 = SdfT(x2y2 + math::Pow2(c.z()-P[2]));//square distance from c to P
735  if (x2y2z2 >= max2 || (!acc.probeValue(c,v) && v<ValueT(0)))
736  continue;//outside narrow band of the particle or inside existing level set
737  if (x2y2z2 <= min2) {//inside narrow band of the particle.
738  acc.setValueOff(c, inside);
739  continue;
740  }
741  // convert signed distance from voxel units to world units
742  //const ValueT d=dx*(math::Sqrt(x2y2z2) - R);
743  const ValueT d=Merge(dx*(math::Sqrt(x2y2z2) - R), att);
744  if (d < v) acc.setValue(c, d);//CSG union
745  }//end loop over z
746  }//end loop over y
747  }//end loop over x
748  return true;
749  }
750  typedef typename hboost::function<void (Raster*, const tbb::blocked_range<size_t>&)> FuncType;
751 
752  template <typename DisableType>
754  getAtt(size_t, AttT&) const {;}
755 
756  template <typename DisableType>
758  getAtt(size_t n, AttT& a) const { mParticles.getAtt(n, a); }
759 
760  template <typename T>
761  typename hboost::enable_if<hboost::is_same<T,ValueT>, ValueT>::type
762  Merge(T s, const AttT&) const { return s; }
763 
764  template <typename T>
765  typename hboost::disable_if<hboost::is_same<T,ValueT>, ValueT>::type
766  Merge(T s, const AttT& a) const { return ValueT(s,a); }
767 
768  ParticlesToLevelSetT& mParent;
769  const ParticleListT& mParticles;//list of particles
770  GridT* mGrid;
771  const math::MapBase& mMap;
772  size_t mMinCount, mMaxCount;//counters for ignored particles!
773  FuncType mTask;
774  const bool mIsCopy;
775  PointPartitionerT* mPointPartitioner;
776 };//end of Raster struct
777 
778 
779 ///////////////////// YOU CAN SAFELY IGNORE THIS SECTION /////////////////////
780 
781 namespace p2ls_internal {
782 
783 // This is a simple type that combines a distance value and a particle
784 // attribute. It's required for attribute transfer which is defined in the
785 // Raster class above.
786 template<typename VisibleT, typename BlindT>
787 class BlindData
788 {
789 public:
790  typedef VisibleT type;
791  typedef VisibleT VisibleType;
792  typedef BlindT BlindType;
793 
795  explicit BlindData(VisibleT v) : mVisible(v), mBlind(zeroVal<BlindType>()) {}
796  BlindData(VisibleT v, BlindT b) : mVisible(v), mBlind(b) {}
797  BlindData(const BlindData&) = default;
798  BlindData& operator=(const BlindData&) = default;
799  const VisibleT& visible() const { return mVisible; }
800  const BlindT& blind() const { return mBlind; }
802  bool operator==(const BlindData& rhs) const { return mVisible == rhs.mVisible; }
804  bool operator< (const BlindData& rhs) const { return mVisible < rhs.mVisible; }
805  bool operator> (const BlindData& rhs) const { return mVisible > rhs.mVisible; }
806  BlindData operator+(const BlindData& rhs) const { return BlindData(mVisible + rhs.mVisible); }
807  BlindData operator+(const VisibleT& rhs) const { return BlindData(mVisible + rhs); }
808  BlindData operator-(const BlindData& rhs) const { return BlindData(mVisible - rhs.mVisible); }
809  BlindData operator-() const { return BlindData(-mVisible, mBlind); }
810 
811 protected:
812  VisibleT mVisible;
813  BlindT mBlind;
814 };
815 
816 // Required by several of the tree nodes
817 template<typename VisibleT, typename BlindT>
818 inline std::ostream& operator<<(std::ostream& ostr, const BlindData<VisibleT, BlindT>& rhs)
819 {
820  ostr << rhs.visible();
821  return ostr;
822 }
823 
824 // Required by math::Abs
825 template<typename VisibleT, typename BlindT>
827 {
828  return BlindData<VisibleT, BlindT>(math::Abs(x.visible()), x.blind());
829 }
830 
831 } // namespace p2ls_internal
832 
833 //////////////////////////////////////////////////////////////////////////////
834 
835 } // namespace tools
836 } // namespace OPENVDB_VERSION_NAME
837 } // namespace openvdb
838 
839 #endif // OPENVDB_TOOLS_PARTICLES_TO_LEVELSET_HAS_BEEN_INCLUDED
840 
841 // Copyright (c) 2012-2017 DreamWorks Animation LLC
842 // All rights reserved. This software is distributed under the
843 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
size_t getMaxCount() const
Return number of large particles that were ignore due to Rmax.
Real getHalfWidth() const
Return the half-width of the narrow band in voxel units.
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:822
Functions to efficiently perform various compositing operations on grids.
math::Vec3< Real > Vec3R
Definition: Types.h:75
Tag dispatch class that distinguishes shallow copy constructors from deep copy constructors.
Definition: Types.h:500
Type Pow2(Type x)
Return .
Definition: Math.h:514
Real getRmax() const
Return the largest radius allowed in voxel units.
AttGridType::Ptr attributeGrid()
Return a shared pointer to the grid containing the (optional) attribute.
Real getVoxelSize() const
Return the size of a voxel in world units.
const GLdouble * v
Definition: glcorearb.h:836
void finalize(bool prune=false)
This methods syncs up the level set and attribute grids and therefore needs to be called before any o...
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
png_infop png_color_16p * background
Definition: png.h:2326
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
GA_API const UT_StringHolder P
Dummy NOOP interrupter class defining interface.
void setGrainSize(int grainSize)
Set the grain-size used for multi-threading.
ParticlesToLevelSet(SdfGridT &grid, InterrupterT *interrupt=nullptr)
Constructor using an exiting signed distance, i.e. narrow band level set, grid.
Spatially partitions points using a parallel radix-based sorting algorithm.
const hboost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
GLuint id
Definition: glcorearb.h:654
SYS_FORCE_INLINE const_iterator end() const
int getGrainSize() const
Returns the grain-size used for multi-threading.
GLdouble n
Definition: glcorearb.h:2007
BlindData< VisibleT, BlindT > Abs(const BlindData< VisibleT, BlindT > &x)
Coord Abs(const Coord &xyz)
Definition: Coord.h:254
#define OPENVDB_VERSION_NAME
Definition: version.h:43
void rasterizeSpheres(const ParticleListT &pa)
Rasterize a sphere per particle derived from their position and radius. All spheres are CSG unioned...
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
void setRmax(Real Rmax)
set the largest radius allowed in voxel units
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
Defined various multi-threaded utility functions for trees.
hboost::mpl::if_< DisableT, size_t, AttributeT >::type AttType
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:82
OPENVDB_NO_FP_EQUALITY_WARNING_END bool operator<(const BlindData &rhs) const
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
GA_API const UT_StringHolder transform
SdfGridT::template ValueConverter< AttType >::Type AttGridType
size_t getMinCount() const
Return number of small particles that were ignore due to Rmin.
GLint GLsizei count
Definition: glcorearb.h:404
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
bool ignoredParticles() const
Return true if any particles were ignored due to their size.
HBOOST_STATIC_ASSERT(hboost::is_floating_point< SdfType >::value)
GLsizei const GLfloat * value
Definition: glcorearb.h:823
OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN bool operator==(const BlindData &rhs) const
void setRmin(Real Rmin)
set the smallest radius allowed in voxel units
void rasterizeTrails(const ParticleListT &pa, Real delta=1.0)
Rasterize a trail per particle derived from their position, radius and velocity. Each trail is genera...
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...
OPENVDB_STATIC_SPECIALIZATION 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:1022
GLint GLenum GLint x
Definition: glcorearb.h:408
GA_API const UT_StringHolder N
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
#define SIZE
Definition: simple.C:40
int Floor(float x)
Return the floor of x.
Definition: Math.h:814
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:856
GLboolean r
Definition: glcorearb.h:1221
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
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:374
T zeroVal()
Return the value of type T that corresponds to zero.
Definition: Math.h:94
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:429
bool wasInterrupted(T *i, int percent=-1)
Real getRmin() const
Return the smallest radius allowed in voxel units.
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:101
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:794
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:83