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