OpenVDB  7.0.0
ParticleAtlas.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
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 
55 
56 
82 template<typename PointIndexGridType = PointIndexGrid>
84 {
87 
88  using PointIndexGridPtr = typename PointIndexGridType::Ptr;
89  using IndexType = typename PointIndexGridType::ValueType;
90 
91  struct Iterator;
92 
94 
95  ParticleAtlas() : mIndexGridArray(), mMinRadiusArray(), mMaxRadiusArray() {}
96 
102  template<typename ParticleArrayType>
103  void construct(const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels = 50);
104 
110  template<typename ParticleArrayType>
111  static Ptr create(const ParticleArrayType& particles,
112  double minVoxelSize, size_t maxLevels = 50);
113 
115  size_t levels() const { return mIndexGridArray.size(); }
117  bool empty() const { return mIndexGridArray.empty(); }
118 
120  double minRadius(size_t n) const { return mMinRadiusArray[n]; }
122  double maxRadius(size_t n) const { return mMaxRadiusArray[n]; }
123 
125  PointIndexGridType& pointIndexGrid(size_t n) { return *mIndexGridArray[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 
143 
144 
150 template<typename PointIndexGridType>
151 struct ParticleAtlas<PointIndexGridType>::Iterator
152 {
153  using TreeType = typename PointIndexGridType::TreeType;
155  using ConstAccessorPtr = std::unique_ptr<ConstAccessor>;
156 
158 
160  explicit Iterator(const ParticleAtlas& atlas);
161 
167  template<typename ParticleArrayType>
168  void worldSpaceSearchAndUpdate(const Vec3d& center, double radius,
169  const ParticleArrayType& particles);
170 
175  template<typename ParticleArrayType>
176  void worldSpaceSearchAndUpdate(const BBoxd& bbox, const ParticleArrayType& particles);
177 
179  size_t levels() const { return mAtlas->levels(); }
180 
183  void updateFromLevel(size_t level);
184 
186  void reset();
187 
189  const IndexType& operator*() const { return *mRange.first; }
190 
193  bool test() const { return mRange.first < mRange.second || mIter != mRangeList.end(); }
194  operator bool() const { return this->test(); }
196 
198  void increment();
199 
201  void operator++() { this->increment(); }
202 
205  bool next();
206 
208  size_t size() const;
209 
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 
239 
240 // Internal operators and implementation details
241 
242 
243 namespace particle_atlas_internal {
244 
245 
246 template<typename ParticleArrayT>
248 {
249  using PosType = typename ParticleArrayT::PosType;
250  using ScalarType = typename PosType::value_type;
251 
252  ComputeExtremas(const ParticleArrayT& particles)
253  : particleArray(&particles)
254  , minRadius(std::numeric_limits<ScalarType>::max())
255  , maxRadius(-std::numeric_limits<ScalarType>::max())
256  {
257  }
258 
260  : particleArray(rhs.particleArray)
261  , minRadius(std::numeric_limits<ScalarType>::max())
262  , maxRadius(-std::numeric_limits<ScalarType>::max())
263  {
264  }
265 
266  void operator()(const tbb::blocked_range<size_t>& range) {
267 
268  ScalarType radius, tmpMin = minRadius, tmpMax = maxRadius;
269 
270  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
271  particleArray->getRadius(n, radius);
272  tmpMin = std::min(radius, tmpMin);
273  tmpMax = std::max(radius, tmpMax);
274  }
275 
276  minRadius = std::min(minRadius, tmpMin);
277  maxRadius = std::max(maxRadius, tmpMax);
278  }
279 
280  void join(const ComputeExtremas& rhs) {
281  minRadius = std::min(minRadius, rhs.minRadius);
282  maxRadius = std::max(maxRadius, rhs.maxRadius);
283  }
284 
285  ParticleArrayT const * const particleArray;
286  ScalarType minRadius, maxRadius;
287 }; // struct ComputeExtremas
288 
289 
290 template<typename ParticleArrayT, typename PointIndex>
292 {
295  using ParticleArray = ParticleArrayT;
296 
297  using PosType = typename ParticleArray::PosType;
298  using ScalarType = typename PosType::value_type;
299 
300  SplittableParticleArray(const ParticleArrayT& particles)
301  : mIndexMap(), mParticleArray(&particles), mSize(particles.size())
302  {
303  updateExtremas();
304  }
305 
306  SplittableParticleArray(const ParticleArrayT& particles, double minR, double maxR)
307  : mIndexMap(), mParticleArray(&particles), mSize(particles.size())
308  {
309  mMinRadius = ScalarType(minR);
310  mMaxRadius = ScalarType(maxR);
311  }
312 
313  const ParticleArrayT& particleArray() const { return *mParticleArray; }
314 
315  size_t size() const { return mSize; }
316 
317  void getPos(size_t n, PosType& xyz) const
318  { return mParticleArray->getPos(getGlobalIndex(n), xyz); }
319  void getRadius(size_t n, ScalarType& radius) const
320  { return mParticleArray->getRadius(getGlobalIndex(n), radius); }
321 
322  ScalarType minRadius() const { return mMinRadius; }
323  ScalarType maxRadius() const { return mMaxRadius; }
324 
325  size_t getGlobalIndex(size_t n) const { return mIndexMap ? size_t(mIndexMap[n]) : n; }
326 
329  Ptr split(ScalarType maxRadiusLimit) {
330 
331  if (mMaxRadius < maxRadiusLimit) return Ptr();
332 
333  std::unique_ptr<bool[]> mask{new bool[mSize]};
334 
335  tbb::parallel_for(tbb::blocked_range<size_t>(0, mSize),
336  MaskParticles(*this, mask, maxRadiusLimit));
337 
338  Ptr output(new SplittableParticleArray(*this, mask));
339  if (output->size() == 0) return Ptr();
340 
341  size_t newSize = 0;
342  for (size_t n = 0, N = mSize; n < N; ++n) {
343  newSize += size_t(!mask[n]);
344  }
345 
346  std::unique_ptr<PointIndex[]> newIndexMap{new PointIndex[newSize]};
347 
348  setIndexMap(newIndexMap, mask, false);
349 
350  mSize = newSize;
351  mIndexMap.swap(newIndexMap);
352  updateExtremas();
353 
354  return output;
355  }
356 
357 
358 private:
359  // Disallow copying
362 
363  // Masked copy constructor
365  const std::unique_ptr<bool[]>& mask)
366  : mIndexMap(), mParticleArray(&other.particleArray()), mSize(0)
367  {
368  for (size_t n = 0, N = other.size(); n < N; ++n) {
369  mSize += size_t(mask[n]);
370  }
371 
372  if (mSize != 0) {
373  mIndexMap.reset(new PointIndex[mSize]);
374  other.setIndexMap(mIndexMap, mask, true);
375  }
376 
377  updateExtremas();
378  }
379 
380  struct MaskParticles {
381  MaskParticles(const SplittableParticleArray& particles,
382  const std::unique_ptr<bool[]>& mask, ScalarType radius)
383  : particleArray(&particles)
384  , particleMask(mask.get())
385  , radiusLimit(radius)
386  {
387  }
388 
389  void operator()(const tbb::blocked_range<size_t>& range) const {
390  const ScalarType maxRadius = radiusLimit;
391  ScalarType radius;
392  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
393  particleArray->getRadius(n, radius);
394  particleMask[n] = !(radius < maxRadius);
395  }
396  }
397 
398  SplittableParticleArray const * const particleArray;
399  bool * const particleMask;
400  ScalarType const radiusLimit;
401  }; // struct MaskParticles
402 
403  inline void updateExtremas() {
405  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mSize), op);
406  mMinRadius = op.minRadius;
407  mMaxRadius = op.maxRadius;
408  }
409 
410  void setIndexMap(std::unique_ptr<PointIndex[]>& newIndexMap,
411  const std::unique_ptr<bool[]>& mask, bool maskValue) const
412  {
413  if (mIndexMap.get() != nullptr) {
414  const PointIndex* indices = mIndexMap.get();
415  for (size_t idx = 0, n = 0, N = mSize; n < N; ++n) {
416  if (mask[n] == maskValue) newIndexMap[idx++] = indices[n];
417  }
418  } else {
419  for (size_t idx = 0, n = 0, N = mSize; n < N; ++n) {
420  if (mask[n] == maskValue) {
421  newIndexMap[idx++] = PointIndex(static_cast<typename PointIndex::IntType>(n));
422  }
423  }
424  }
425  }
426 
427 
429 
430  std::unique_ptr<PointIndex[]> mIndexMap;
431  ParticleArrayT const * const mParticleArray;
432  size_t mSize;
433  ScalarType mMinRadius, mMaxRadius;
434 }; // struct SplittableParticleArray
435 
436 
437 template<typename ParticleArrayType, typename PointIndexLeafNodeType>
438 struct RemapIndices {
439 
440  RemapIndices(const ParticleArrayType& particles, std::vector<PointIndexLeafNodeType*>& nodes)
441  : mParticles(&particles)
442  , mNodes(nodes.empty() ? nullptr : &nodes.front())
443  {
444  }
445 
446  void operator()(const tbb::blocked_range<size_t>& range) const
447  {
448  using PointIndexType = typename PointIndexLeafNodeType::ValueType;
449  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
450 
451  PointIndexLeafNodeType& node = *mNodes[n];
452  const size_t numIndices = node.indices().size();
453 
454  if (numIndices > 0) {
455  PointIndexType* begin = &node.indices().front();
456  const PointIndexType* end = begin + numIndices;
457 
458  while (begin < end) {
459  *begin = PointIndexType(static_cast<typename PointIndexType::IntType>(
460  mParticles->getGlobalIndex(*begin)));
461  ++begin;
462  }
463  }
464  }
465  }
466 
467  ParticleArrayType const * const mParticles;
468  PointIndexLeafNodeType * const * const mNodes;
469 }; // struct RemapIndices
470 
471 
472 template<typename ParticleArrayType, typename IndexT>
474 {
475  using PosType = typename ParticleArrayType::PosType;
476  using ScalarType = typename PosType::value_type;
477 
478  using Range = std::pair<const IndexT*, const IndexT*>;
479  using RangeDeque = std::deque<Range>;
480  using IndexDeque = std::deque<IndexT>;
481 
482  RadialRangeFilter(RangeDeque& ranges, IndexDeque& indices, const PosType& xyz,
483  ScalarType radius, const ParticleArrayType& particles, bool hasUniformRadius = false)
484  : mRanges(ranges)
485  , mIndices(indices)
486  , mCenter(xyz)
487  , mRadius(radius)
488  , mParticles(&particles)
489  , mHasUniformRadius(hasUniformRadius)
490  {
491  if (mHasUniformRadius) {
492  ScalarType uniformRadius;
493  mParticles->getRadius(0, uniformRadius);
494  mRadius = mRadius + uniformRadius;
495  mRadius *= mRadius;
496  }
497  }
498 
499  template <typename LeafNodeType>
500  void filterLeafNode(const LeafNodeType& leaf)
501  {
502  const size_t numIndices = leaf.indices().size();
503  if (numIndices > 0) {
504  const IndexT* begin = &leaf.indices().front();
505  filterVoxel(leaf.origin(), begin, begin + numIndices);
506  }
507  }
508 
509  void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
510  {
511  PosType pos;
512 
513  if (mHasUniformRadius) {
514 
515  const ScalarType searchRadiusSqr = mRadius;
516 
517  while (begin < end) {
518  mParticles->getPos(size_t(*begin), pos);
519  const ScalarType distSqr = (mCenter - pos).lengthSqr();
520  if (distSqr < searchRadiusSqr) {
521  mIndices.push_back(*begin);
522  }
523  ++begin;
524  }
525  } else {
526  while (begin < end) {
527  const size_t idx = size_t(*begin);
528  mParticles->getPos(idx, pos);
529 
530  ScalarType radius;
531  mParticles->getRadius(idx, radius);
532 
533  ScalarType searchRadiusSqr = mRadius + radius;
534  searchRadiusSqr *= searchRadiusSqr;
535 
536  const ScalarType distSqr = (mCenter - pos).lengthSqr();
537 
538  if (distSqr < searchRadiusSqr) {
539  mIndices.push_back(*begin);
540  }
541 
542  ++begin;
543  }
544  }
545  }
546 
547 private:
549  RadialRangeFilter& operator=(const RadialRangeFilter&);
550 
551  RangeDeque& mRanges;
552  IndexDeque& mIndices;
553  PosType const mCenter;
554  ScalarType mRadius;
555  ParticleArrayType const * const mParticles;
556  bool const mHasUniformRadius;
557 }; // struct RadialRangeFilter
558 
559 
560 template<typename ParticleArrayType, typename IndexT>
562 {
563  using PosType = typename ParticleArrayType::PosType;
564  using ScalarType = typename PosType::value_type;
565 
566  using Range = std::pair<const IndexT*, const IndexT*>;
567  using RangeDeque = std::deque<Range>;
568  using IndexDeque = std::deque<IndexT>;
569 
570  BBoxFilter(RangeDeque& ranges, IndexDeque& indices,
571  const BBoxd& bbox, const ParticleArrayType& particles, bool hasUniformRadius = false)
572  : mRanges(ranges)
573  , mIndices(indices)
574  , mBBox(PosType(bbox.min()), PosType(bbox.max()))
575  , mCenter(mBBox.getCenter())
576  , mParticles(&particles)
577  , mHasUniformRadius(hasUniformRadius)
578  , mUniformRadiusSqr(ScalarType(0.0))
579  {
580  if (mHasUniformRadius) {
581  mParticles->getRadius(0, mUniformRadiusSqr);
582  mUniformRadiusSqr *= mUniformRadiusSqr;
583  }
584  }
585 
586  template <typename LeafNodeType>
587  void filterLeafNode(const LeafNodeType& leaf)
588  {
589  const size_t numIndices = leaf.indices().size();
590  if (numIndices > 0) {
591  const IndexT* begin = &leaf.indices().front();
592  filterVoxel(leaf.origin(), begin, begin + numIndices);
593  }
594  }
595 
596  void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
597  {
598  PosType pos;
599 
600  if (mHasUniformRadius) {
601  const ScalarType radiusSqr = mUniformRadiusSqr;
602 
603  while (begin < end) {
604 
605  mParticles->getPos(size_t(*begin), pos);
606  if (mBBox.isInside(pos)) {
607  mIndices.push_back(*begin++);
608  continue;
609  }
610 
611  const ScalarType distSqr = pointToBBoxDistSqr(pos);
612  if (!(distSqr > radiusSqr)) {
613  mIndices.push_back(*begin);
614  }
615 
616  ++begin;
617  }
618 
619  } else {
620  while (begin < end) {
621 
622  const size_t idx = size_t(*begin);
623  mParticles->getPos(idx, pos);
624  if (mBBox.isInside(pos)) {
625  mIndices.push_back(*begin++);
626  continue;
627  }
628 
629  ScalarType radius;
630  mParticles->getRadius(idx, radius);
631  const ScalarType distSqr = pointToBBoxDistSqr(pos);
632  if (!(distSqr > (radius * radius))) {
633  mIndices.push_back(*begin);
634  }
635 
636  ++begin;
637  }
638  }
639  }
640 
641 private:
642  BBoxFilter(const BBoxFilter&);
643  BBoxFilter& operator=(const BBoxFilter&);
644 
645  ScalarType pointToBBoxDistSqr(const PosType& pos) const
646  {
647  ScalarType distSqr = ScalarType(0.0);
648 
649  for (int i = 0; i < 3; ++i) {
650 
651  const ScalarType a = pos[i];
652 
653  ScalarType b = mBBox.min()[i];
654  if (a < b) {
655  ScalarType delta = b - a;
656  distSqr += delta * delta;
657  }
658 
659  b = mBBox.max()[i];
660  if (a > b) {
661  ScalarType delta = a - b;
662  distSqr += delta * delta;
663  }
664  }
665 
666  return distSqr;
667  }
668 
669  RangeDeque& mRanges;
670  IndexDeque& mIndices;
671  math::BBox<PosType> const mBBox;
672  PosType const mCenter;
673  ParticleArrayType const * const mParticles;
674  bool const mHasUniformRadius;
675  ScalarType mUniformRadiusSqr;
676 }; // struct BBoxFilter
677 
678 
679 } // namespace particle_atlas_internal
680 
681 
683 
684 
685 template<typename PointIndexGridType>
686 template<typename ParticleArrayType>
687 inline void
689  const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels)
690 {
691  using SplittableParticleArray =
693  using SplittableParticleArrayPtr = typename SplittableParticleArray::Ptr;
694  using ScalarType = typename ParticleArrayType::ScalarType;
695 
697 
699  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, particles.size()), extremas);
700  const double firstMin = extremas.minRadius;
701  const double firstMax = extremas.maxRadius;
702  const double firstVoxelSize = std::max(minVoxelSize, firstMin);
703 
704  if (!(firstMax < (firstVoxelSize * double(2.0))) && maxLevels > 1) {
705 
706  std::vector<SplittableParticleArrayPtr> levels;
707  levels.push_back(SplittableParticleArrayPtr(
708  new SplittableParticleArray(particles, firstMin, firstMax)));
709 
710  std::vector<double> voxelSizeArray;
711  voxelSizeArray.push_back(firstVoxelSize);
712 
713  for (size_t n = 0; n < maxLevels; ++n) {
714 
715  const double maxParticleRadius = double(levels.back()->maxRadius());
716  const double particleRadiusLimit = voxelSizeArray.back() * double(2.0);
717  if (maxParticleRadius < particleRadiusLimit) break;
718 
719  SplittableParticleArrayPtr newLevel =
720  levels.back()->split(ScalarType(particleRadiusLimit));
721  if (!newLevel) break;
722 
723  levels.push_back(newLevel);
724  voxelSizeArray.push_back(double(newLevel->minRadius()));
725  }
726 
727  size_t numPoints = 0;
728 
729  using PointIndexTreeType = typename PointIndexGridType::TreeType;
730  using PointIndexLeafNodeType = typename PointIndexTreeType::LeafNodeType;
731 
732  std::vector<PointIndexLeafNodeType*> nodes;
733 
734  for (size_t n = 0, N = levels.size(); n < N; ++n) {
735 
736  const SplittableParticleArray& particleArray = *levels[n];
737 
738  numPoints += particleArray.size();
739 
740  mMinRadiusArray.push_back(double(particleArray.minRadius()));
741  mMaxRadiusArray.push_back(double(particleArray.maxRadius()));
742 
743  PointIndexGridPtr grid =
744  createPointIndexGrid<PointIndexGridType>(particleArray, voxelSizeArray[n]);
745 
746  nodes.clear();
747  grid->tree().getNodes(nodes);
748 
749  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
750  particle_atlas_internal::RemapIndices<SplittableParticleArray,
751  PointIndexLeafNodeType>(particleArray, nodes));
752 
753  mIndexGridArray.push_back(grid);
754  }
755 
756  } else {
757  mMinRadiusArray.push_back(firstMin);
758  mMaxRadiusArray.push_back(firstMax);
759  mIndexGridArray.push_back(
760  createPointIndexGrid<PointIndexGridType>(particles, firstVoxelSize));
761  }
762 }
763 
764 
765 template<typename PointIndexGridType>
766 template<typename ParticleArrayType>
769  const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels)
770 {
771  Ptr ret(new ParticleAtlas());
772  ret->construct(particles, minVoxelSize, maxLevels);
773  return ret;
774 }
775 
776 
778 
779 // ParticleAtlas::Iterator implementation
780 
781 template<typename PointIndexGridType>
782 inline
784  : mAtlas(&atlas)
785  , mAccessorList()
786  , mRange(static_cast<IndexType*>(nullptr), static_cast<IndexType*>(nullptr))
787  , mRangeList()
788  , mIter(mRangeList.begin())
789  , mIndexArray()
790  , mIndexArraySize(0)
791  , mAccessorListSize(atlas.levels())
792 {
793  if (mAccessorListSize > 0) {
794  mAccessorList.reset(new ConstAccessorPtr[mAccessorListSize]);
795  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
796  mAccessorList[n].reset(new ConstAccessor(atlas.pointIndexGrid(n).tree()));
797  }
798  }
799 }
800 
801 
802 template<typename PointIndexGridType>
803 inline void
805 {
806  mIter = mRangeList.begin();
807  if (!mRangeList.empty()) {
808  mRange = mRangeList.front();
809  } else if (mIndexArray) {
810  mRange.first = mIndexArray.get();
811  mRange.second = mRange.first + mIndexArraySize;
812  } else {
813  mRange.first = static_cast<IndexType*>(nullptr);
814  mRange.second = static_cast<IndexType*>(nullptr);
815  }
816 }
817 
818 
819 template<typename PointIndexGridType>
820 inline void
822 {
823  ++mRange.first;
824  if (mRange.first >= mRange.second && mIter != mRangeList.end()) {
825  ++mIter;
826  if (mIter != mRangeList.end()) {
827  mRange = *mIter;
828  } else if (mIndexArray) {
829  mRange.first = mIndexArray.get();
830  mRange.second = mRange.first + mIndexArraySize;
831  }
832  }
833 }
834 
835 
836 template<typename PointIndexGridType>
837 inline bool
839 {
840  if (!this->test()) return false;
841  this->increment();
842  return this->test();
843 }
844 
845 
846 template<typename PointIndexGridType>
847 inline size_t
849 {
850  size_t count = 0;
851  typename RangeDeque::const_iterator it =
852  mRangeList.begin(), end = mRangeList.end();
853 
854  for ( ; it != end; ++it) {
855  count += it->second - it->first;
856  }
857 
858  return count + mIndexArraySize;
859 }
860 
861 
862 template<typename PointIndexGridType>
863 inline void
865 {
866  mRange.first = static_cast<IndexType*>(nullptr);
867  mRange.second = static_cast<IndexType*>(nullptr);
868  mRangeList.clear();
869  mIter = mRangeList.end();
870  mIndexArray.reset();
871  mIndexArraySize = 0;
872 }
873 
874 
875 template<typename PointIndexGridType>
876 inline void
878 {
879  using TreeT = typename PointIndexGridType::TreeType;
880  using LeafNodeType = typename TreeType::LeafNodeType;
881 
882  this->clear();
883 
884  if (mAccessorListSize > 0) {
885  const size_t levelIdx = std::min(mAccessorListSize - 1, level);
886 
887  const TreeT& tree = mAtlas->pointIndexGrid(levelIdx).tree();
888 
889  std::vector<const LeafNodeType*> nodes;
890  tree.getNodes(nodes);
891 
892  for (size_t n = 0, N = nodes.size(); n < N; ++n) {
893 
894  const LeafNodeType& node = *nodes[n];
895  const size_t numIndices = node.indices().size();
896 
897  if (numIndices > 0) {
898  const IndexType* begin = &node.indices().front();
899  mRangeList.push_back(Range(begin, (begin + numIndices)));
900  }
901  }
902  }
903 
904  this->reset();
905 }
906 
907 
908 template<typename PointIndexGridType>
909 template<typename ParticleArrayType>
910 inline void
912  const Vec3d& center, double radius, const ParticleArrayType& particles)
913 {
914  using PosType = typename ParticleArrayType::PosType;
915  using ScalarType = typename ParticleArrayType::ScalarType;
916 
918 
919  this->clear();
920 
921  std::deque<IndexType> filteredIndices;
922  std::vector<CoordBBox> searchRegions;
923 
924  const double iRadius = radius * double(1.0 / std::sqrt(3.0));
925 
926  const Vec3d ibMin(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius);
927  const Vec3d ibMax(center[0] + iRadius, center[1] + iRadius, center[2] + iRadius);
928 
929  const Vec3d bMin(center[0] - radius, center[1] - radius, center[2] - radius);
930  const Vec3d bMax(center[0] + radius, center[1] + radius, center[2] + radius);
931 
932  const PosType pos = PosType(center);
933  const ScalarType dist = ScalarType(radius);
934 
935  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
936 
937  const double maxRadius = mAtlas->maxRadius(n);
938  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
939 
940  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
941 
942  ConstAccessor& acc = *mAccessorList[n];
943 
944  openvdb::CoordBBox inscribedRegion(
945  xform.worldToIndexCellCentered(ibMin),
946  xform.worldToIndexCellCentered(ibMax));
947 
948  inscribedRegion.expand(-1); // erode by one voxel
949 
950  // collect indices that don't need to be tested
951  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
952 
953  searchRegions.clear();
954 
955  const openvdb::CoordBBox region(
956  xform.worldToIndexCellCentered(bMin - maxRadius),
957  xform.worldToIndexCellCentered(bMax + maxRadius));
958 
959  inscribedRegion.expand(1);
961  searchRegions, region, inscribedRegion);
962 
964  FilterType filter(mRangeList, filteredIndices, pos, dist, particles, uniformRadius);
965 
966  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
967  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
968  }
969  }
970 
971  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
972 
973  this->reset();
974 }
975 
976 
977 template<typename PointIndexGridType>
978 template<typename ParticleArrayType>
979 inline void
981  const BBoxd& bbox, const ParticleArrayType& particles)
982 {
983  this->clear();
984 
985  std::deque<IndexType> filteredIndices;
986  std::vector<CoordBBox> searchRegions;
987 
988  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
989 
990  const double maxRadius = mAtlas->maxRadius(n);
991  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
992  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
993 
994  ConstAccessor& acc = *mAccessorList[n];
995 
996  openvdb::CoordBBox inscribedRegion(
997  xform.worldToIndexCellCentered(bbox.min()),
998  xform.worldToIndexCellCentered(bbox.max()));
999 
1000  inscribedRegion.expand(-1); // erode by one voxel
1001 
1002  // collect indices that don't need to be tested
1003  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
1004 
1005  searchRegions.clear();
1006 
1007  const openvdb::CoordBBox region(
1008  xform.worldToIndexCellCentered(bbox.min() - maxRadius),
1009  xform.worldToIndexCellCentered(bbox.max() + maxRadius));
1010 
1011  inscribedRegion.expand(1);
1013  searchRegions, region, inscribedRegion);
1014 
1016  FilterType filter(mRangeList, filteredIndices, bbox, particles, uniformRadius);
1017 
1018  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
1019  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
1020  }
1021  }
1022 
1023  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1024 
1025  this->reset();
1026 }
1027 
1028 
1029 } // namespace tools
1030 } // namespace OPENVDB_VERSION_NAME
1031 } // namespace openvdb
1032 
1033 #endif // OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
BBoxFilter(RangeDeque &ranges, IndexDeque &indices, const BBoxd &bbox, const ParticleArrayType &particles, bool hasUniformRadius=false)
Definition: ParticleAtlas.h:570
PointIndexLeafNodeType *const *const mNodes
Definition: ParticleAtlas.h:468
std::deque< Range > RangeDeque
Definition: ParticleAtlas.h:479
tree::ValueAccessor< const TreeType > ConstAccessor
Definition: ParticleAtlas.h:154
std::pair< const IndexT *, const IndexT * > Range
Definition: ParticleAtlas.h:566
void filterLeafNode(const LeafNodeType &leaf)
Definition: ParticleAtlas.h:587
Integer wrapper, required to distinguish PointIndexGrid and PointDataGrid from Int32Grid and Int64Gri...
Definition: Types.h:133
typename PointIndexGridType::Ptr PointIndexGridPtr
Definition: ParticleAtlas.h:88
void join(const ComputeExtremas &rhs)
Definition: ParticleAtlas.h:280
const ParticleArrayT & particleArray() const
Definition: ParticleAtlas.h:313
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
typename ParticleArrayType::PosType PosType
Definition: ParticleAtlas.h:475
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:351
bool next()
Advance iterator to next item.
Definition: ParticleAtlas.h:838
size_t levels() const
Returns the number of resolution levels.
Definition: ParticleAtlas.h:115
Vec3< double > Vec3d
Definition: Vec3.h:662
Space-partitioning acceleration structure for points. Partitions the points into voxels to accelerate...
std::deque< Range > RangeDeque
Definition: ParticleAtlas.h:567
typename ParticleArray::PosType PosType
Definition: ParticleAtlas.h:297
std::shared_ptr< T > SharedPtr
Definition: Types.h:91
Provides accelerated range and nearest-neighbor searches for particles that are partitioned using the...
Definition: ParticleAtlas.h:151
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:877
Ptr split(ScalarType maxRadiusLimit)
Definition: ParticleAtlas.h:329
typename PointIndexGridType::TreeType TreeType
Definition: ParticleAtlas.h:153
Definition: Coord.h:587
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:62
std::pair< const IndexT *, const IndexT * > Range
Definition: ParticleAtlas.h:478
void getPos(size_t n, PosType &xyz) const
Definition: ParticleAtlas.h:317
size_t size() const
Return the number of point indices in the iterator range.
Definition: ParticleAtlas.h:848
SharedPtr< SplittableParticleArray > Ptr
Definition: ParticleAtlas.h:293
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:476
SharedPtr< const SplittableParticleArray > ConstPtr
Definition: ParticleAtlas.h:294
void increment()
Advance iterator to next item.
Definition: ParticleAtlas.h:821
bool test() const
Return true if this iterator is not yet exhausted.
Definition: ParticleAtlas.h:193
PointIndexGridType & pointIndexGrid(size_t n)
Returns the PointIndexGrid that represents the given level n.
Definition: ParticleAtlas.h:125
SharedPtr< ParticleAtlas > Ptr
Definition: ParticleAtlas.h:85
size_t getGlobalIndex(size_t n) const
Definition: ParticleAtlas.h:325
Selectively extract and filter point data using a custom filter operator.
void dequeToArray(const std::deque< T > &d, std::unique_ptr< T[]> &a, size_t &size)
Definition: PointIndexGrid.h:528
size_t levels() const
Returns the total number of resolution levels.
Definition: ParticleAtlas.h:179
std::unique_ptr< ConstAccessor > ConstAccessorPtr
Definition: ParticleAtlas.h:155
ComputeExtremas(const ParticleArrayT &particles)
Definition: ParticleAtlas.h:252
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
ParticleArrayT const *const particleArray
Definition: ParticleAtlas.h:285
ScalarType minRadius
Definition: ParticleAtlas.h:286
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:911
bool operator!=(const Iterator &p) const
Definition: ParticleAtlas.h:212
void filteredPointIndexSearch(RangeFilterType &filter, ConstAccessor &acc, const CoordBBox &bbox)
Definition: PointIndexGrid.h:800
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:298
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:564
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
Definition: ParticleAtlas.h:596
const PointIndexGridType & pointIndexGrid(size_t n) const
Returns the PointIndexGrid that represents the given level n.
Definition: ParticleAtlas.h:127
Definition: Exceptions.h:13
std::deque< IndexT > IndexDeque
Definition: ParticleAtlas.h:480
ComputeExtremas(ComputeExtremas &rhs, tbb::split)
Definition: ParticleAtlas.h:259
typename ParticleArrayType::PosType PosType
Definition: ParticleAtlas.h:563
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:64
typename PointIndexGridType::ValueType IndexType
Definition: ParticleAtlas.h:89
std::deque< IndexT > IndexDeque
Definition: ParticleAtlas.h:568
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: ParticleAtlas.h:446
void operator()(const tbb::blocked_range< size_t > &range)
Definition: ParticleAtlas.h:266
Partition particles and performs range and nearest-neighbor searches.
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:250
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
Definition: ParticleAtlas.h:509
SplittableParticleArray(const ParticleArrayT &particles, double minR, double maxR)
Definition: ParticleAtlas.h:306
RadialRangeFilter(RangeDeque &ranges, IndexDeque &indices, const PosType &xyz, ScalarType radius, const ParticleArrayType &particles, bool hasUniformRadius=false)
Definition: ParticleAtlas.h:482
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:471
ScalarType maxRadius
Definition: ParticleAtlas.h:286
SplittableParticleArray(const ParticleArrayT &particles)
Definition: ParticleAtlas.h:300
ScalarType minRadius() const
Definition: ParticleAtlas.h:322
ParticleAtlas()
Definition: ParticleAtlas.h:95
SharedPtr< const ParticleAtlas > ConstPtr
Definition: ParticleAtlas.h:86
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:102
void reset()
Reset the iterator to point to the first item.
Definition: ParticleAtlas.h:804
void constructExclusiveRegions(std::vector< CoordBBox > &regions, const CoordBBox &bbox, const CoordBBox &ibox)
Definition: PointIndexGrid.h:539
void pointIndexSearch(RangeDeque &rangeList, ConstAccessor &acc, const CoordBBox &bbox)
Definition: PointIndexGrid.h:876
bool operator==(const Iterator &p) const
Return true if both iterators point to the same element.
Definition: ParticleAtlas.h:211
typename ParticleArrayT::PosType PosType
Definition: ParticleAtlas.h:249
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
void filterLeafNode(const LeafNodeType &leaf)
Definition: ParticleAtlas.h:500
double minRadius(size_t n) const
Returns minimum particle radius for level n.
Definition: ParticleAtlas.h:120
ScalarType maxRadius() const
Definition: ParticleAtlas.h:323
size_t size() const
Definition: ParticleAtlas.h:315
ParticleArrayType const *const mParticles
Definition: ParticleAtlas.h:467
RemapIndices(const ParticleArrayType &particles, std::vector< PointIndexLeafNodeType * > &nodes)
Definition: ParticleAtlas.h:440
void getRadius(size_t n, ScalarType &radius) const
Definition: ParticleAtlas.h:319
const IndexType & operator*() const
Return a const reference to the item to which this iterator is pointing.
Definition: ParticleAtlas.h:189
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106
Definition: ParticleAtlas.h:83