HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ParticleAtlas.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 /// @file ParticleAtlas.h
32 ///
33 /// @brief Space-partitioning acceleration structure for particles, points with
34 /// radius. Partitions particle indices into voxels to accelerate range
35 /// and nearest neighbor searches.
36 ///
37 /// @note   This acceleration structure only stores integer offsets into an external particle
38 /// data structure that conforms to the ParticleArray interface. 
39 ///
40 /// @details Constructs and maintains a sequence of @c PointIndexGrids each of progressively
41 /// lower resolution. Particles are uniquely assigned to a particular resolution
42 /// level based on their radius. This strategy has proven efficient for accelerating
43 /// spatial queries on particle data sets with varying radii.
44 ///
45 /// @details The data structure automatically detects and adapts to particle data sets with
46 /// uniform radii. The construction is simplified and spatial queries pre-cache the
47 /// uniform particle radius to avoid redundant access calls to the
48 /// ParticleArray::getRadius method.
49 ///
50 /// @author Mihai Alden
51 
52 #ifndef OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
53 #define OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
54 
55 #include "PointIndexGrid.h"
56 
57 #include <openvdb/Grid.h>
58 #include <openvdb/Types.h>
59 #include <openvdb/math/Transform.h>
60 #include <openvdb/tree/Tree.h>
61 #include <openvdb/tree/LeafNode.h>
62 
63 #include <tbb/blocked_range.h>
64 #include <tbb/parallel_for.h>
65 #include <tbb/parallel_reduce.h>
66 #include <algorithm> // for std::min(), std::max()
67 #include <cmath> // for std::sqrt()
68 #include <deque>
69 #include <limits>
70 #include <memory>
71 #include <utility> // for std::pair
72 #include <vector>
73 
74 
75 namespace openvdb {
77 namespace OPENVDB_VERSION_NAME {
78 namespace tools {
79 
80 
81 ////////////////////////////////////////
82 
83 
84 /// @brief Partition particles and performs range and nearest-neighbor searches.
85 ///
86 /// @interface ParticleArray
87 /// Expected interface for the ParticleArray container:
88 /// @code
89 /// template<typename VectorType>
90 /// struct ParticleArray
91 /// {
92 /// // The type used to represent world-space positions
93 /// using PosType = VectorType;
94 /// using ScalarType = typename PosType::value_type;
95 ///
96 /// // Return the number of particles in the array
97 /// size_t size() const;
98 ///
99 /// // Return the world-space position for the nth particle.
100 /// void getPos(size_t n, PosType& xyz) const;
101 ///
102 /// // Return the world-space radius for the nth particle.
103 /// void getRadius(size_t n, ScalarType& radius) const;
104 /// };
105 /// @endcode
106 ///
107 /// @details Constructs a collection of @c PointIndexGrids of different resolutions
108 /// to accelerate spatial searches for particles with varying radius.
109 template<typename PointIndexGridType = PointIndexGrid>
111 {
114 
115  using PointIndexGridPtr = typename PointIndexGridType::Ptr;
116  using IndexType = typename PointIndexGridType::ValueType;
117 
118  struct Iterator;
119 
120  //////////
121 
122  ParticleAtlas() : mIndexGridArray(), mMinRadiusArray(), mMaxRadiusArray() {}
123 
124  /// @brief Partitions particle indices
125  ///
126  /// @param particles container conforming to the ParticleArray interface
127  /// @param minVoxelSize minimum voxel size limit
128  /// @param maxLevels maximum number of resolution levels
129  template<typename ParticleArrayType>
130  void construct(const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels = 50);
131 
132  /// @brief Create a new @c ParticleAtlas from the given @a particles.
133  ///
134  /// @param particles container conforming to the ParticleArray interface
135  /// @param minVoxelSize minimum voxel size limit
136  /// @param maxLevels maximum number of resolution levels
137  template<typename ParticleArrayType>
138  static Ptr create(const ParticleArrayType& particles,
139  double minVoxelSize, size_t maxLevels = 50);
140 
141  /// @brief Returns the number of resolution levels.
142  size_t levels() const { return mIndexGridArray.size(); }
143  /// @brief true if the container size is 0, false otherwise.
144  bool empty() const { return mIndexGridArray.empty(); }
145 
146  /// @brief Returns minimum particle radius for level @a n.
147  double minRadius(size_t n) const { return mMinRadiusArray[n]; }
148  /// @brief Returns maximum particle radius for level @a n.
149  double maxRadius(size_t n) const { return mMaxRadiusArray[n]; }
150 
151  /// @brief Returns the @c PointIndexGrid that represents the given level @a n.
152  PointIndexGridType& pointIndexGrid(size_t n) { return *mIndexGridArray[n]; }
153  /// @brief Returns the @c PointIndexGrid that represents the given level @a n.
154  const PointIndexGridType& pointIndexGrid(size_t n) const { return *mIndexGridArray[n]; }
155 
156 private:
157  // Disallow copying
159  ParticleAtlas& operator=(const ParticleAtlas&);
160 
161  std::vector<PointIndexGridPtr> mIndexGridArray;
162  std::vector<double> mMinRadiusArray, mMaxRadiusArray;
163 }; // struct ParticleAtlas
164 
165 
167 
168 
169 ////////////////////////////////////////
170 
171 
172 /// @brief Provides accelerated range and nearest-neighbor searches for
173 /// particles that are partitioned using the ParticleAtlas.
174 ///
175 /// @note Prefer to construct the iterator object once and reuse it
176 /// for subsequent queries.
177 template<typename PointIndexGridType>
178 struct ParticleAtlas<PointIndexGridType>::Iterator
179 {
180  using TreeType = typename PointIndexGridType::TreeType;
182  using ConstAccessorPtr = std::unique_ptr<ConstAccessor>;
183 
184  /////
185 
186  /// @brief Construct an iterator from the given @a atlas.
187  explicit Iterator(const ParticleAtlas& atlas);
188 
189  /// @brief Clear the iterator and update it with the result of the given
190  /// world-space radial query.
191  /// @param center world-space center
192  /// @param radius world-space search radius
193  /// @param particles container conforming to the ParticleArray interface
194  template<typename ParticleArrayType>
195  void worldSpaceSearchAndUpdate(const Vec3d& center, double radius,
196  const ParticleArrayType& particles);
197 
198  /// @brief Clear the iterator and update it with the result of the given
199  /// world-space radial query.
200  /// @param bbox world-space bounding box
201  /// @param particles container conforming to the ParticleArray interface
202  template<typename ParticleArrayType>
203  void worldSpaceSearchAndUpdate(const BBoxd& bbox, const ParticleArrayType& particles);
204 
205  /// @brief Returns the total number of resolution levels.
206  size_t levels() const { return mAtlas->levels(); }
207 
208  /// @brief Clear the iterator and update it with all particles that reside
209  /// at the given resolution @a level.
210  void updateFromLevel(size_t level);
211 
212  /// Reset the iterator to point to the first item.
213  void reset();
214 
215  /// Return a const reference to the item to which this iterator is pointing.
216  const IndexType& operator*() const { return *mRange.first; }
217 
218  /// @{
219  /// @brief Return @c true if this iterator is not yet exhausted.
220  bool test() const { return mRange.first < mRange.second || mIter != mRangeList.end(); }
221  operator bool() const { return this->test(); }
222  /// @}
223 
224  /// Advance iterator to next item.
225  void increment();
226 
227  /// Advance iterator to next item.
228  void operator++() { this->increment(); }
229 
230  /// @brief Advance iterator to next item.
231  /// @return @c true if this iterator is not yet exhausted.
232  bool next();
233 
234  /// Return the number of point indices in the iterator range.
235  size_t size() const;
236 
237  /// Return @c true if both iterators point to the same element.
238  bool operator==(const Iterator& p) const { return mRange.first == p.mRange.first; }
239  bool operator!=(const Iterator& p) const { return !this->operator==(p); }
240 
241 private:
242  Iterator(const Iterator& rhs);
243  Iterator& operator=(const Iterator& rhs);
244 
245  void clear();
246 
247  using Range = std::pair<const IndexType*, const IndexType*>;
248  using RangeDeque = std::deque<Range>;
249  using RangeDequeCIter = typename RangeDeque::const_iterator;
250  using IndexArray = std::unique_ptr<IndexType[]>;
251 
252  ParticleAtlas const * const mAtlas;
253  std::unique_ptr<ConstAccessorPtr[]> mAccessorList;
254 
255  // Primary index collection
256  Range mRange;
257  RangeDeque mRangeList;
258  RangeDequeCIter mIter;
259  // Secondary index collection
260  IndexArray mIndexArray;
261  size_t mIndexArraySize, mAccessorListSize;
262 }; // struct ParticleAtlas::Iterator
263 
264 
265 ////////////////////////////////////////
266 
267 // Internal operators and implementation details
268 
269 
270 namespace particle_atlas_internal {
271 
272 
273 template<typename ParticleArrayT>
275 {
276  using PosType = typename ParticleArrayT::PosType;
277  using ScalarType = typename PosType::value_type;
278 
279  ComputeExtremas(const ParticleArrayT& particles)
280  : particleArray(&particles)
281  , minRadius(std::numeric_limits<ScalarType>::max())
282  , maxRadius(-std::numeric_limits<ScalarType>::max())
283  {
284  }
285 
288  , minRadius(std::numeric_limits<ScalarType>::max())
289  , maxRadius(-std::numeric_limits<ScalarType>::max())
290  {
291  }
292 
293  void operator()(const tbb::blocked_range<size_t>& range) {
294 
295  ScalarType radius, tmpMin = minRadius, tmpMax = maxRadius;
296 
297  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
298  particleArray->getRadius(n, radius);
299  tmpMin = std::min(radius, tmpMin);
300  tmpMax = std::max(radius, tmpMax);
301  }
302 
303  minRadius = std::min(minRadius, tmpMin);
304  maxRadius = std::max(maxRadius, tmpMax);
305  }
306 
307  void join(const ComputeExtremas& rhs) {
310  }
311 
312  ParticleArrayT const * const particleArray;
314 }; // struct ComputeExtremas
315 
316 
317 template<typename ParticleArrayT, typename PointIndex>
319 {
322  using ParticleArray = ParticleArrayT;
323 
324  using PosType = typename ParticleArray::PosType;
325  using ScalarType = typename PosType::value_type;
326 
327  SplittableParticleArray(const ParticleArrayT& particles)
328  : mIndexMap(), mParticleArray(&particles), mSize(particles.size())
329  {
330  updateExtremas();
331  }
332 
333  SplittableParticleArray(const ParticleArrayT& particles, double minR, double maxR)
334  : mIndexMap(), mParticleArray(&particles), mSize(particles.size())
335  {
336  mMinRadius = ScalarType(minR);
337  mMaxRadius = ScalarType(maxR);
338  }
339 
340  const ParticleArrayT& particleArray() const { return *mParticleArray; }
341 
342  size_t size() const { return mSize; }
343 
344  void getPos(size_t n, PosType& xyz) const
345  { return mParticleArray->getPos(getGlobalIndex(n), xyz); }
346  void getRadius(size_t n, ScalarType& radius) const
347  { return mParticleArray->getRadius(getGlobalIndex(n), radius); }
348 
349  ScalarType minRadius() const { return mMinRadius; }
350  ScalarType maxRadius() const { return mMaxRadius; }
351 
352  size_t getGlobalIndex(size_t n) const { return mIndexMap ? size_t(mIndexMap[n]) : n; }
353 
354  /// Move all particle indices that have a radius larger or equal to @a maxRadiusLimit
355  /// into a separate container.
356  Ptr split(ScalarType maxRadiusLimit) {
357 
358  if (mMaxRadius < maxRadiusLimit) return Ptr();
359 
360  std::unique_ptr<bool[]> mask{new bool[mSize]};
361 
362  tbb::parallel_for(tbb::blocked_range<size_t>(0, mSize),
363  MaskParticles(*this, mask, maxRadiusLimit));
364 
365  Ptr output(new SplittableParticleArray(*this, mask));
366  if (output->size() == 0) return Ptr();
367 
368  size_t newSize = 0;
369  for (size_t n = 0, N = mSize; n < N; ++n) {
370  newSize += size_t(!mask[n]);
371  }
372 
373  std::unique_ptr<PointIndex[]> newIndexMap{new PointIndex[newSize]};
374 
375  setIndexMap(newIndexMap, mask, false);
376 
377  mSize = newSize;
378  mIndexMap.swap(newIndexMap);
379  updateExtremas();
380 
381  return output;
382  }
383 
384 
385 private:
386  // Disallow copying
389 
390  // Masked copy constructor
392  const std::unique_ptr<bool[]>& mask)
393  : mIndexMap(), mParticleArray(&other.particleArray()), mSize(0)
394  {
395  for (size_t n = 0, N = other.size(); n < N; ++n) {
396  mSize += size_t(mask[n]);
397  }
398 
399  if (mSize != 0) {
400  mIndexMap.reset(new PointIndex[mSize]);
401  other.setIndexMap(mIndexMap, mask, true);
402  }
403 
404  updateExtremas();
405  }
406 
407  struct MaskParticles {
408  MaskParticles(const SplittableParticleArray& particles,
409  const std::unique_ptr<bool[]>& mask, ScalarType radius)
410  : particleArray(&particles)
411  , particleMask(mask.get())
412  , radiusLimit(radius)
413  {
414  }
415 
416  void operator()(const tbb::blocked_range<size_t>& range) const {
417  const ScalarType maxRadius = radiusLimit;
418  ScalarType radius;
419  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
420  particleArray->getRadius(n, radius);
421  particleMask[n] = !(radius < maxRadius);
422  }
423  }
424 
426  bool * const particleMask;
427  ScalarType const radiusLimit;
428  }; // struct MaskParticles
429 
430  inline void updateExtremas() {
431  ComputeExtremas<SplittableParticleArray> op(*this);
432  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mSize), op);
433  mMinRadius = op.minRadius;
434  mMaxRadius = op.maxRadius;
435  }
436 
437  void setIndexMap(std::unique_ptr<PointIndex[]>& newIndexMap,
438  const std::unique_ptr<bool[]>& mask, bool maskValue) const
439  {
440  if (mIndexMap.get() != nullptr) {
441  const PointIndex* indices = mIndexMap.get();
442  for (size_t idx = 0, n = 0, N = mSize; n < N; ++n) {
443  if (mask[n] == maskValue) newIndexMap[idx++] = indices[n];
444  }
445  } else {
446  for (size_t idx = 0, n = 0, N = mSize; n < N; ++n) {
447  if (mask[n] == maskValue) {
448  newIndexMap[idx++] = PointIndex(static_cast<typename PointIndex::IntType>(n));
449  }
450  }
451  }
452  }
453 
454 
455  //////////
456 
457  std::unique_ptr<PointIndex[]> mIndexMap;
458  ParticleArrayT const * const mParticleArray;
459  size_t mSize;
460  ScalarType mMinRadius, mMaxRadius;
461 }; // struct SplittableParticleArray
462 
463 
464 template<typename ParticleArrayType, typename PointIndexLeafNodeType>
465 struct RemapIndices {
466 
467  RemapIndices(const ParticleArrayType& particles, std::vector<PointIndexLeafNodeType*>& nodes)
468  : mParticles(&particles)
469  , mNodes(nodes.empty() ? nullptr : &nodes.front())
470  {
471  }
472 
473  void operator()(const tbb::blocked_range<size_t>& range) const
474  {
475  using PointIndexType = typename PointIndexLeafNodeType::ValueType;
476  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
477 
478  PointIndexLeafNodeType& node = *mNodes[n];
479  const size_t numIndices = node.indices().size();
480 
481  if (numIndices > 0) {
482  PointIndexType* begin = &node.indices().front();
483  const PointIndexType* end = begin + numIndices;
484 
485  while (begin < end) {
486  *begin = PointIndexType(static_cast<typename PointIndexType::IntType>(
487  mParticles->getGlobalIndex(*begin)));
488  ++begin;
489  }
490  }
491  }
492  }
493 
494  ParticleArrayType const * const mParticles;
495  PointIndexLeafNodeType * const * const mNodes;
496 }; // struct RemapIndices
497 
498 
499 template<typename ParticleArrayType, typename IndexT>
501 {
502  using PosType = typename ParticleArrayType::PosType;
503  using ScalarType = typename PosType::value_type;
504 
505  using Range = std::pair<const IndexT*, const IndexT*>;
506  using RangeDeque = std::deque<Range>;
507  using IndexDeque = std::deque<IndexT>;
508 
509  RadialRangeFilter(RangeDeque& ranges, IndexDeque& indices, const PosType& xyz,
510  ScalarType radius, const ParticleArrayType& particles, bool hasUniformRadius = false)
511  : mRanges(ranges)
512  , mIndices(indices)
513  , mCenter(xyz)
514  , mRadius(radius)
515  , mParticles(&particles)
516  , mHasUniformRadius(hasUniformRadius)
517  {
518  if (mHasUniformRadius) {
519  ScalarType uniformRadius;
520  mParticles->getRadius(0, uniformRadius);
521  mRadius = mRadius + uniformRadius;
522  mRadius *= mRadius;
523  }
524  }
525 
526  template <typename LeafNodeType>
527  void filterLeafNode(const LeafNodeType& leaf)
528  {
529  const size_t numIndices = leaf.indices().size();
530  if (numIndices > 0) {
531  const IndexT* begin = &leaf.indices().front();
532  filterVoxel(leaf.origin(), begin, begin + numIndices);
533  }
534  }
535 
536  void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
537  {
538  PosType pos;
539 
540  if (mHasUniformRadius) {
541 
542  const ScalarType searchRadiusSqr = mRadius;
543 
544  while (begin < end) {
545  mParticles->getPos(size_t(*begin), pos);
546  const ScalarType distSqr = (mCenter - pos).lengthSqr();
547  if (distSqr < searchRadiusSqr) {
548  mIndices.push_back(*begin);
549  }
550  ++begin;
551  }
552  } else {
553  while (begin < end) {
554  const size_t idx = size_t(*begin);
555  mParticles->getPos(idx, pos);
556 
557  ScalarType radius;
558  mParticles->getRadius(idx, radius);
559 
560  ScalarType searchRadiusSqr = mRadius + radius;
561  searchRadiusSqr *= searchRadiusSqr;
562 
563  const ScalarType distSqr = (mCenter - pos).lengthSqr();
564 
565  if (distSqr < searchRadiusSqr) {
566  mIndices.push_back(*begin);
567  }
568 
569  ++begin;
570  }
571  }
572  }
573 
574 private:
576  RadialRangeFilter& operator=(const RadialRangeFilter&);
577 
578  RangeDeque& mRanges;
579  IndexDeque& mIndices;
580  PosType const mCenter;
581  ScalarType mRadius;
582  ParticleArrayType const * const mParticles;
583  bool const mHasUniformRadius;
584 }; // struct RadialRangeFilter
585 
586 
587 template<typename ParticleArrayType, typename IndexT>
589 {
590  using PosType = typename ParticleArrayType::PosType;
591  using ScalarType = typename PosType::value_type;
592 
593  using Range = std::pair<const IndexT*, const IndexT*>;
594  using RangeDeque = std::deque<Range>;
595  using IndexDeque = std::deque<IndexT>;
596 
597  BBoxFilter(RangeDeque& ranges, IndexDeque& indices,
598  const BBoxd& bbox, const ParticleArrayType& particles, bool hasUniformRadius = false)
599  : mRanges(ranges)
600  , mIndices(indices)
601  , mBBox(PosType(bbox.min()), PosType(bbox.max()))
602  , mCenter(mBBox.getCenter())
603  , mParticles(&particles)
604  , mHasUniformRadius(hasUniformRadius)
605  , mUniformRadiusSqr(ScalarType(0.0))
606  {
607  if (mHasUniformRadius) {
608  mParticles->getRadius(0, mUniformRadiusSqr);
609  mUniformRadiusSqr *= mUniformRadiusSqr;
610  }
611  }
612 
613  template <typename LeafNodeType>
614  void filterLeafNode(const LeafNodeType& leaf)
615  {
616  const size_t numIndices = leaf.indices().size();
617  if (numIndices > 0) {
618  const IndexT* begin = &leaf.indices().front();
619  filterVoxel(leaf.origin(), begin, begin + numIndices);
620  }
621  }
622 
623  void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
624  {
625  PosType pos;
626 
627  if (mHasUniformRadius) {
628  const ScalarType radiusSqr = mUniformRadiusSqr;
629 
630  while (begin < end) {
631 
632  mParticles->getPos(size_t(*begin), pos);
633  if (mBBox.isInside(pos)) {
634  mIndices.push_back(*begin++);
635  continue;
636  }
637 
638  const ScalarType distSqr = pointToBBoxDistSqr(pos);
639  if (!(distSqr > radiusSqr)) {
640  mIndices.push_back(*begin);
641  }
642 
643  ++begin;
644  }
645 
646  } else {
647  while (begin < end) {
648 
649  const size_t idx = size_t(*begin);
650  mParticles->getPos(idx, pos);
651  if (mBBox.isInside(pos)) {
652  mIndices.push_back(*begin++);
653  continue;
654  }
655 
656  ScalarType radius;
657  mParticles->getRadius(idx, radius);
658  const ScalarType distSqr = pointToBBoxDistSqr(pos);
659  if (!(distSqr > (radius * radius))) {
660  mIndices.push_back(*begin);
661  }
662 
663  ++begin;
664  }
665  }
666  }
667 
668 private:
669  BBoxFilter(const BBoxFilter&);
670  BBoxFilter& operator=(const BBoxFilter&);
671 
672  ScalarType pointToBBoxDistSqr(const PosType& pos) const
673  {
674  ScalarType distSqr = ScalarType(0.0);
675 
676  for (int i = 0; i < 3; ++i) {
677 
678  const ScalarType a = pos[i];
679 
680  ScalarType b = mBBox.min()[i];
681  if (a < b) {
682  ScalarType delta = b - a;
683  distSqr += delta * delta;
684  }
685 
686  b = mBBox.max()[i];
687  if (a > b) {
688  ScalarType delta = a - b;
689  distSqr += delta * delta;
690  }
691  }
692 
693  return distSqr;
694  }
695 
696  RangeDeque& mRanges;
697  IndexDeque& mIndices;
698  math::BBox<PosType> const mBBox;
699  PosType const mCenter;
700  ParticleArrayType const * const mParticles;
701  bool const mHasUniformRadius;
702  ScalarType mUniformRadiusSqr;
703 }; // struct BBoxFilter
704 
705 
706 } // namespace particle_atlas_internal
707 
708 
709 ////////////////////////////////////////
710 
711 
712 template<typename PointIndexGridType>
713 template<typename ParticleArrayType>
714 inline void
716  const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels)
717 {
718  using SplittableParticleArray =
720  using SplittableParticleArrayPtr = typename SplittableParticleArray::Ptr;
721  using ScalarType = typename ParticleArrayType::ScalarType;
722 
723  /////
724 
726  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, particles.size()), extremas);
727  const double firstMin = extremas.minRadius;
728  const double firstMax = extremas.maxRadius;
729  const double firstVoxelSize = std::max(minVoxelSize, firstMin);
730 
731  if (!(firstMax < (firstVoxelSize * double(2.0))) && maxLevels > 1) {
732 
733  std::vector<SplittableParticleArrayPtr> levels;
734  levels.push_back(SplittableParticleArrayPtr(
735  new SplittableParticleArray(particles, firstMin, firstMax)));
736 
737  std::vector<double> voxelSizeArray;
738  voxelSizeArray.push_back(firstVoxelSize);
739 
740  for (size_t n = 0; n < maxLevels; ++n) {
741 
742  const double maxParticleRadius = double(levels.back()->maxRadius());
743  const double particleRadiusLimit = voxelSizeArray.back() * double(2.0);
744  if (maxParticleRadius < particleRadiusLimit) break;
745 
746  SplittableParticleArrayPtr newLevel =
747  levels.back()->split(ScalarType(particleRadiusLimit));
748  if (!newLevel) break;
749 
750  levels.push_back(newLevel);
751  voxelSizeArray.push_back(double(newLevel->minRadius()));
752  }
753 
754  size_t numPoints = 0;
755 
756  using PointIndexTreeType = typename PointIndexGridType::TreeType;
757  using PointIndexLeafNodeType = typename PointIndexTreeType::LeafNodeType;
758 
759  std::vector<PointIndexLeafNodeType*> nodes;
760 
761  for (size_t n = 0, N = levels.size(); n < N; ++n) {
762 
763  const SplittableParticleArray& particleArray = *levels[n];
764 
765  numPoints += particleArray.size();
766 
767  mMinRadiusArray.push_back(double(particleArray.minRadius()));
768  mMaxRadiusArray.push_back(double(particleArray.maxRadius()));
769 
770  PointIndexGridPtr grid =
771  createPointIndexGrid<PointIndexGridType>(particleArray, voxelSizeArray[n]);
772 
773  nodes.clear();
774  grid->tree().getNodes(nodes);
775 
776  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
777  particle_atlas_internal::RemapIndices<SplittableParticleArray,
778  PointIndexLeafNodeType>(particleArray, nodes));
779 
780  mIndexGridArray.push_back(grid);
781  }
782 
783  } else {
784  mMinRadiusArray.push_back(firstMin);
785  mMaxRadiusArray.push_back(firstMax);
786  mIndexGridArray.push_back(
787  createPointIndexGrid<PointIndexGridType>(particles, firstVoxelSize));
788  }
789 }
790 
791 
792 template<typename PointIndexGridType>
793 template<typename ParticleArrayType>
796  const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels)
797 {
798  Ptr ret(new ParticleAtlas());
799  ret->construct(particles, minVoxelSize, maxLevels);
800  return ret;
801 }
802 
803 
804 ////////////////////////////////////////
805 
806 // ParticleAtlas::Iterator implementation
807 
808 template<typename PointIndexGridType>
809 inline
811  : mAtlas(&atlas)
812  , mAccessorList()
813  , mRange(static_cast<IndexType*>(nullptr), static_cast<IndexType*>(nullptr))
814  , mRangeList()
815  , mIter(mRangeList.begin())
816  , mIndexArray()
817  , mIndexArraySize(0)
818  , mAccessorListSize(atlas.levels())
819 {
820  if (mAccessorListSize > 0) {
821  mAccessorList.reset(new ConstAccessorPtr[mAccessorListSize]);
822  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
823  mAccessorList[n].reset(new ConstAccessor(atlas.pointIndexGrid(n).tree()));
824  }
825  }
826 }
827 
828 
829 template<typename PointIndexGridType>
830 inline void
832 {
833  mIter = mRangeList.begin();
834  if (!mRangeList.empty()) {
835  mRange = mRangeList.front();
836  } else if (mIndexArray) {
837  mRange.first = mIndexArray.get();
838  mRange.second = mRange.first + mIndexArraySize;
839  } else {
840  mRange.first = static_cast<IndexType*>(nullptr);
841  mRange.second = static_cast<IndexType*>(nullptr);
842  }
843 }
844 
845 
846 template<typename PointIndexGridType>
847 inline void
849 {
850  ++mRange.first;
851  if (mRange.first >= mRange.second && mIter != mRangeList.end()) {
852  ++mIter;
853  if (mIter != mRangeList.end()) {
854  mRange = *mIter;
855  } else if (mIndexArray) {
856  mRange.first = mIndexArray.get();
857  mRange.second = mRange.first + mIndexArraySize;
858  }
859  }
860 }
861 
862 
863 template<typename PointIndexGridType>
864 inline bool
866 {
867  if (!this->test()) return false;
868  this->increment();
869  return this->test();
870 }
871 
872 
873 template<typename PointIndexGridType>
874 inline size_t
876 {
877  size_t count = 0;
878  typename RangeDeque::const_iterator it =
879  mRangeList.begin(), end = mRangeList.end();
880 
881  for ( ; it != end; ++it) {
882  count += it->second - it->first;
883  }
884 
885  return count + mIndexArraySize;
886 }
887 
888 
889 template<typename PointIndexGridType>
890 inline void
892 {
893  mRange.first = static_cast<IndexType*>(nullptr);
894  mRange.second = static_cast<IndexType*>(nullptr);
895  mRangeList.clear();
896  mIter = mRangeList.end();
897  mIndexArray.reset();
898  mIndexArraySize = 0;
899 }
900 
901 
902 template<typename PointIndexGridType>
903 inline void
905 {
906  using TreeType = typename PointIndexGridType::TreeType;
907  using LeafNodeType = typename TreeType::LeafNodeType;
908 
909  this->clear();
910 
911  if (mAccessorListSize > 0) {
912  const size_t levelIdx = std::min(mAccessorListSize - 1, level);
913 
914  const TreeType& tree = mAtlas->pointIndexGrid(levelIdx).tree();
915 
916 
917  std::vector<const LeafNodeType*> nodes;
918  tree.getNodes(nodes);
919 
920  for (size_t n = 0, N = nodes.size(); n < N; ++n) {
921 
922  const LeafNodeType& node = *nodes[n];
923  const size_t numIndices = node.indices().size();
924 
925  if (numIndices > 0) {
926  const IndexType* begin = &node.indices().front();
927  mRangeList.push_back(Range(begin, (begin + numIndices)));
928  }
929  }
930  }
931 
932  this->reset();
933 }
934 
935 
936 template<typename PointIndexGridType>
937 template<typename ParticleArrayType>
938 inline void
940  const Vec3d& center, double radius, const ParticleArrayType& particles)
941 {
942  using PosType = typename ParticleArrayType::PosType;
943  using ScalarType = typename ParticleArrayType::ScalarType;
944 
945  /////
946 
947  this->clear();
948 
949  std::deque<IndexType> filteredIndices;
950  std::vector<CoordBBox> searchRegions;
951 
952  const double iRadius = radius * double(1.0 / std::sqrt(3.0));
953 
954  const Vec3d ibMin(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius);
955  const Vec3d ibMax(center[0] + iRadius, center[1] + iRadius, center[2] + iRadius);
956 
957  const Vec3d bMin(center[0] - radius, center[1] - radius, center[2] - radius);
958  const Vec3d bMax(center[0] + radius, center[1] + radius, center[2] + radius);
959 
960  const PosType pos = PosType(center);
961  const ScalarType dist = ScalarType(radius);
962 
963  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
964 
965  const double maxRadius = mAtlas->maxRadius(n);
966  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
967 
968  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
969 
970  ConstAccessor& acc = *mAccessorList[n];
971 
972  openvdb::CoordBBox inscribedRegion(
973  xform.worldToIndexCellCentered(ibMin),
974  xform.worldToIndexCellCentered(ibMax));
975 
976  inscribedRegion.expand(-1); // erode by one voxel
977 
978  // collect indices that don't need to be tested
979  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
980 
981  searchRegions.clear();
982 
983  const openvdb::CoordBBox region(
984  xform.worldToIndexCellCentered(bMin - maxRadius),
985  xform.worldToIndexCellCentered(bMax + maxRadius));
986 
987  inscribedRegion.expand(1);
989  searchRegions, region, inscribedRegion);
990 
992  FilterType filter(mRangeList, filteredIndices, pos, dist, particles, uniformRadius);
993 
994  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
995  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
996  }
997  }
998 
999  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1000 
1001  this->reset();
1002 }
1003 
1004 
1005 template<typename PointIndexGridType>
1006 template<typename ParticleArrayType>
1007 inline void
1009  const BBoxd& bbox, const ParticleArrayType& particles)
1010 {
1011  this->clear();
1012 
1013  std::deque<IndexType> filteredIndices;
1014  std::vector<CoordBBox> searchRegions;
1015 
1016  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
1017 
1018  const double maxRadius = mAtlas->maxRadius(n);
1019  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
1020  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
1021 
1022  ConstAccessor& acc = *mAccessorList[n];
1023 
1024  openvdb::CoordBBox inscribedRegion(
1025  xform.worldToIndexCellCentered(bbox.min()),
1026  xform.worldToIndexCellCentered(bbox.max()));
1027 
1028  inscribedRegion.expand(-1); // erode by one voxel
1029 
1030  // collect indices that don't need to be tested
1031  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
1032 
1033  searchRegions.clear();
1034 
1035  const openvdb::CoordBBox region(
1036  xform.worldToIndexCellCentered(bbox.min() - maxRadius),
1037  xform.worldToIndexCellCentered(bbox.max() + maxRadius));
1038 
1039  inscribedRegion.expand(1);
1041  searchRegions, region, inscribedRegion);
1042 
1044  FilterType filter(mRangeList, filteredIndices, bbox, particles, uniformRadius);
1045 
1046  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
1047  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
1048  }
1049  }
1050 
1051  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1052 
1053  this->reset();
1054 }
1055 
1056 
1057 } // namespace tools
1058 } // namespace OPENVDB_VERSION_NAME
1059 } // namespace openvdb
1060 
1061 #endif // OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
1062 
1063 // Copyright (c) 2012-2018 DreamWorks Animation LLC
1064 // All rights reserved. This software is distributed under the
1065 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
GA_API const UT_StringHolder dist
cvex test(vector P=0;int unbound=3;export float s=0;export vector Cf=0;)
Definition: test.vfl:11
GLsizei GLenum const void * indices
Definition: glcorearb.h:405
GLenum GLint * range
Definition: glcorearb.h:1924
SplittableParticleArray(const ParticleArrayT &particles, double minR, double maxR)
Provides accelerated range and nearest-neighbor searches for particles that are partitioned using the...
BBoxFilter(RangeDeque &ranges, IndexDeque &indices, const BBoxd &bbox, const ParticleArrayType &particles, bool hasUniformRadius=false)
void dequeToArray(const std::deque< T > &d, std::unique_ptr< T[]> &a, size_t &size)
RemapIndices(const ParticleArrayType &particles, std::vector< PointIndexLeafNodeType * > &nodes)
RadialRangeFilter(RangeDeque &ranges, IndexDeque &indices, const PosType &xyz, ScalarType radius, const ParticleArrayType &particles, bool hasUniformRadius=false)
size_t size() const
Return the number of point indices in the iterator range.
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
GLint level
Definition: glcorearb.h:107
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
GLint GLuint mask
Definition: glcorearb.h:123
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
Selectively extract and filter point data using a custom filter operator.
png_uint_32 i
Definition: png.h:2877
uint64 value_type
Definition: GA_PrimCompat.h:29
GLsizeiptr size
Definition: glcorearb.h:663
PointIndexGridType & pointIndexGrid(size_t n)
Returns the PointIndexGrid that represents the given level n.
void updateFromLevel(size_t level)
Clear the iterator and update it with all particles that reside at the given resolution level...
std::shared_ptr< T > SharedPtr
Definition: Types.h:139
SYS_FORCE_INLINE const_iterator end() const
bool operator==(const BaseDimensions< T > &a, const BaseDimensions< Y > &b)
Definition: Dimensions.h:137
void operator()(const tbb::blocked_range< size_t > &range) const
double minRadius(size_t n) const
Returns minimum particle radius for level n.
GLdouble n
Definition: glcorearb.h:2007
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:89
const IndexType & operator*() const
Return a const reference to the item to which this iterator is pointing.
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
Iterator(const ParticleAtlas &atlas)
Construct an iterator from the given atlas.
size_t levels() const
Returns the number of resolution levels.
GLuint GLuint end
Definition: glcorearb.h:474
const PointIndexGridType & pointIndexGrid(size_t n) const
Returns the PointIndexGrid that represents the given level n.
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
void construct(const ParticleArrayType &particles, double minVoxelSize, size_t maxLevels=50)
Partitions particle indices.
static Ptr create(const ParticleArrayType &particles, double minVoxelSize, size_t maxLevels=50)
Create a new ParticleAtlas from the given particles.
bool empty() const
true if the container size is 0, false otherwise.
size_t levels() const
Returns the total number of resolution levels.
void pointIndexSearch(RangeDeque &rangeList, ConstAccessor &acc, const CoordBBox &bbox)
void constructExclusiveRegions(std::vector< CoordBBox > &regions, const CoordBBox &bbox, const CoordBBox &ibox)
Partition particles and performs range and nearest-neighbor searches.
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
GLsizei levels
Definition: glcorearb.h:2223
GLint GLsizei count
Definition: glcorearb.h:404
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:91
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:358
void worldSpaceSearchAndUpdate(const Vec3d &center, double radius, const ParticleArrayType &particles)
Clear the iterator and update it with the result of the given world-space radial query.
void filteredPointIndexSearch(RangeFilterType &filter, ConstAccessor &acc, const CoordBBox &bbox)
Space-partitioning acceleration structure for points. Partitions the points into voxels to accelerate...
void reset()
Reset the iterator to point to the first item.
GA_API const UT_StringHolder N
Integer wrapper, required to distinguish PointIndexGrid and PointDataGrid from Int32Grid and Int64Gri...
Definition: Types.h:183
bool test() const
Return true if this iterator is not yet exhausted.
typename PointIndexGridType::ValueType IndexType
double maxRadius(size_t n) const
Returns maximum particle radius for level n.
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
bool operator==(const Iterator &p) const
Return true if both iterators point to the same element.
bool isInside(const Vec3T &xyz) const
Return true if the given point is inside this bounding box.
Definition: BBox.h:294
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition: glcorearb.h:1296
typename PointIndexGridType::Ptr PointIndexGridPtr