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