OpenVDB  6.1.0
ParticleAtlas.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2019 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 //
30 
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 
82 
83 
109 template<typename PointIndexGridType = PointIndexGrid>
111 {
114 
115  using PointIndexGridPtr = typename PointIndexGridType::Ptr;
116  using IndexType = typename PointIndexGridType::ValueType;
117 
118  struct Iterator;
119 
121 
122  ParticleAtlas() : mIndexGridArray(), mMinRadiusArray(), mMaxRadiusArray() {}
123 
129  template<typename ParticleArrayType>
130  void construct(const ParticleArrayType& particles, double minVoxelSize, size_t maxLevels = 50);
131 
137  template<typename ParticleArrayType>
138  static Ptr create(const ParticleArrayType& particles,
139  double minVoxelSize, size_t maxLevels = 50);
140 
142  size_t levels() const { return mIndexGridArray.size(); }
144  bool empty() const { return mIndexGridArray.empty(); }
145 
147  double minRadius(size_t n) const { return mMinRadiusArray[n]; }
149  double maxRadius(size_t n) const { return mMaxRadiusArray[n]; }
150 
152  PointIndexGridType& pointIndexGrid(size_t n) { return *mIndexGridArray[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 
170 
171 
177 template<typename PointIndexGridType>
178 struct ParticleAtlas<PointIndexGridType>::Iterator
179 {
180  using TreeType = typename PointIndexGridType::TreeType;
182  using ConstAccessorPtr = std::unique_ptr<ConstAccessor>;
183 
185 
187  explicit Iterator(const ParticleAtlas& atlas);
188 
194  template<typename ParticleArrayType>
195  void worldSpaceSearchAndUpdate(const Vec3d& center, double radius,
196  const ParticleArrayType& particles);
197 
202  template<typename ParticleArrayType>
203  void worldSpaceSearchAndUpdate(const BBoxd& bbox, const ParticleArrayType& particles);
204 
206  size_t levels() const { return mAtlas->levels(); }
207 
210  void updateFromLevel(size_t level);
211 
213  void reset();
214 
216  const IndexType& operator*() const { return *mRange.first; }
217 
220  bool test() const { return mRange.first < mRange.second || mIter != mRangeList.end(); }
221  operator bool() const { return this->test(); }
223 
225  void increment();
226 
228  void operator++() { this->increment(); }
229 
232  bool next();
233 
235  size_t size() const;
236 
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 
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 
287  : particleArray(rhs.particleArray)
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) {
308  minRadius = std::min(minRadius, rhs.minRadius);
309  maxRadius = std::max(maxRadius, rhs.maxRadius);
310  }
311 
312  ParticleArrayT const * const particleArray;
313  ScalarType minRadius, maxRadius;
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 
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 
425  SplittableParticleArray const * const particleArray;
426  bool * const particleMask;
427  ScalarType const radiusLimit;
428  }; // struct MaskParticles
429 
430  inline void updateExtremas() {
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 
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 
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 
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 
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 TreeT = 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 TreeT& tree = mAtlas->pointIndexGrid(levelIdx).tree();
915 
916  std::vector<const LeafNodeType*> nodes;
917  tree.getNodes(nodes);
918 
919  for (size_t n = 0, N = nodes.size(); n < N; ++n) {
920 
921  const LeafNodeType& node = *nodes[n];
922  const size_t numIndices = node.indices().size();
923 
924  if (numIndices > 0) {
925  const IndexType* begin = &node.indices().front();
926  mRangeList.push_back(Range(begin, (begin + numIndices)));
927  }
928  }
929  }
930 
931  this->reset();
932 }
933 
934 
935 template<typename PointIndexGridType>
936 template<typename ParticleArrayType>
937 inline void
939  const Vec3d& center, double radius, const ParticleArrayType& particles)
940 {
941  using PosType = typename ParticleArrayType::PosType;
942  using ScalarType = typename ParticleArrayType::ScalarType;
943 
945 
946  this->clear();
947 
948  std::deque<IndexType> filteredIndices;
949  std::vector<CoordBBox> searchRegions;
950 
951  const double iRadius = radius * double(1.0 / std::sqrt(3.0));
952 
953  const Vec3d ibMin(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius);
954  const Vec3d ibMax(center[0] + iRadius, center[1] + iRadius, center[2] + iRadius);
955 
956  const Vec3d bMin(center[0] - radius, center[1] - radius, center[2] - radius);
957  const Vec3d bMax(center[0] + radius, center[1] + radius, center[2] + radius);
958 
959  const PosType pos = PosType(center);
960  const ScalarType dist = ScalarType(radius);
961 
962  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
963 
964  const double maxRadius = mAtlas->maxRadius(n);
965  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
966 
967  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
968 
969  ConstAccessor& acc = *mAccessorList[n];
970 
971  openvdb::CoordBBox inscribedRegion(
972  xform.worldToIndexCellCentered(ibMin),
973  xform.worldToIndexCellCentered(ibMax));
974 
975  inscribedRegion.expand(-1); // erode by one voxel
976 
977  // collect indices that don't need to be tested
978  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
979 
980  searchRegions.clear();
981 
982  const openvdb::CoordBBox region(
983  xform.worldToIndexCellCentered(bMin - maxRadius),
984  xform.worldToIndexCellCentered(bMax + maxRadius));
985 
986  inscribedRegion.expand(1);
988  searchRegions, region, inscribedRegion);
989 
991  FilterType filter(mRangeList, filteredIndices, pos, dist, particles, uniformRadius);
992 
993  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
994  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
995  }
996  }
997 
998  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
999 
1000  this->reset();
1001 }
1002 
1003 
1004 template<typename PointIndexGridType>
1005 template<typename ParticleArrayType>
1006 inline void
1008  const BBoxd& bbox, const ParticleArrayType& particles)
1009 {
1010  this->clear();
1011 
1012  std::deque<IndexType> filteredIndices;
1013  std::vector<CoordBBox> searchRegions;
1014 
1015  for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
1016 
1017  const double maxRadius = mAtlas->maxRadius(n);
1018  const bool uniformRadius = math::isApproxEqual(mAtlas->minRadius(n), maxRadius);
1019  const openvdb::math::Transform& xform = mAtlas->pointIndexGrid(n).transform();
1020 
1021  ConstAccessor& acc = *mAccessorList[n];
1022 
1023  openvdb::CoordBBox inscribedRegion(
1024  xform.worldToIndexCellCentered(bbox.min()),
1025  xform.worldToIndexCellCentered(bbox.max()));
1026 
1027  inscribedRegion.expand(-1); // erode by one voxel
1028 
1029  // collect indices that don't need to be tested
1030  point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
1031 
1032  searchRegions.clear();
1033 
1034  const openvdb::CoordBBox region(
1035  xform.worldToIndexCellCentered(bbox.min() - maxRadius),
1036  xform.worldToIndexCellCentered(bbox.max() + maxRadius));
1037 
1038  inscribedRegion.expand(1);
1040  searchRegions, region, inscribedRegion);
1041 
1043  FilterType filter(mRangeList, filteredIndices, bbox, particles, uniformRadius);
1044 
1045  for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
1046  point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
1047  }
1048  }
1049 
1050  point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1051 
1052  this->reset();
1053 }
1054 
1055 
1056 } // namespace tools
1057 } // namespace OPENVDB_VERSION_NAME
1058 } // namespace openvdb
1059 
1060 #endif // OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
1061 
1062 // Copyright (c) 2012-2019 DreamWorks Animation LLC
1063 // All rights reserved. This software is distributed under the
1064 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
typename PointIndexGridType::TreeType TreeType
Definition: ParticleAtlas.h:180
ParticleAtlas()
Definition: ParticleAtlas.h:122
void filterLeafNode(const LeafNodeType &leaf)
Definition: ParticleAtlas.h:614
std::deque< IndexT > IndexDeque
Definition: ParticleAtlas.h:595
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: ParticleAtlas.h:473
size_t levels() const
Returns the number of resolution levels.
Definition: ParticleAtlas.h:142
bool test() const
Return true if this iterator is not yet exhausted.
Definition: ParticleAtlas.h:220
PointIndexGridType & pointIndexGrid(size_t n)
Returns the PointIndexGrid that represents the given level n.
Definition: ParticleAtlas.h:152
typename PointIndexGridType::Ptr PointIndexGridPtr
Definition: ParticleAtlas.h:115
ComputeExtremas(ComputeExtremas &rhs, tbb::split)
Definition: ParticleAtlas.h:286
Space-partitioning acceleration structure for points. Partitions the points into voxels to accelerate...
void updateFromLevel(size_t level)
Clear the iterator and update it with all particles that reside at the given resolution level...
Definition: ParticleAtlas.h:904
SharedPtr< const ParticleAtlas > ConstPtr
Definition: ParticleAtlas.h:113
const ParticleArrayT & particleArray() const
Definition: ParticleAtlas.h:340
void filterLeafNode(const LeafNodeType &leaf)
Definition: ParticleAtlas.h:527
Integer wrapper, required to distinguish PointIndexGrid and PointDataGrid from Int32Grid and Int64Gri...
Definition: Types.h:183
double minRadius(size_t n) const
Returns minimum particle radius for level n.
Definition: ParticleAtlas.h:147
Definition: Coord.h:608
Ptr split(ScalarType maxRadiusLimit)
Definition: ParticleAtlas.h:356
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:277
size_t getGlobalIndex(size_t n) const
Definition: ParticleAtlas.h:352
typename ParticleArrayT::PosType PosType
Definition: ParticleAtlas.h:276
RadialRangeFilter(RangeDeque &ranges, IndexDeque &indices, const PosType &xyz, ScalarType radius, const ParticleArrayType &particles, bool hasUniformRadius=false)
Definition: ParticleAtlas.h:509
void constructExclusiveRegions(std::vector< CoordBBox > &regions, const CoordBBox &bbox, const CoordBBox &ibox)
Definition: PointIndexGrid.h:566
ScalarType maxRadius
Definition: ParticleAtlas.h:313
void getPos(size_t n, PosType &xyz) const
Definition: ParticleAtlas.h:344
std::shared_ptr< T > SharedPtr
Definition: Types.h:139
void getRadius(size_t n, ScalarType &radius) const
Definition: ParticleAtlas.h:346
size_t size() const
Definition: ParticleAtlas.h:342
void operator++()
Advance iterator to next item.
Definition: ParticleAtlas.h:228
const PointIndexGridType & pointIndexGrid(size_t n) const
Returns the PointIndexGrid that represents the given level n.
Definition: ParticleAtlas.h:154
RemapIndices(const ParticleArrayType &particles, std::vector< PointIndexLeafNodeType * > &nodes)
Definition: ParticleAtlas.h:467
bool next()
Advance iterator to next item.
Definition: ParticleAtlas.h:865
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
Definition: ParticleAtlas.h:623
Selectively extract and filter point data using a custom filter operator.
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
ParticleArrayType const *const mParticles
Definition: ParticleAtlas.h:494
std::pair< const IndexT *, const IndexT * > Range
Definition: ParticleAtlas.h:593
Vec3< double > Vec3d
Definition: Vec3.h:689
std::pair< const IndexT *, const IndexT * > Range
Definition: ParticleAtlas.h:505
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
SplittableParticleArray(const ParticleArrayT &particles, double minR, double maxR)
Definition: ParticleAtlas.h:333
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:125
typename PointIndexGridType::ValueType IndexType
Definition: ParticleAtlas.h:116
void filteredPointIndexSearch(RangeFilterType &filter, ConstAccessor &acc, const CoordBBox &bbox)
Definition: PointIndexGrid.h:827
SharedPtr< const SplittableParticleArray > ConstPtr
Definition: ParticleAtlas.h:321
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:325
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:938
Provides accelerated range and nearest-neighbor searches for particles that are partitioned using the...
Definition: ParticleAtlas.h:178
std::deque< Range > RangeDeque
Definition: ParticleAtlas.h:594
ScalarType minRadius() const
Definition: ParticleAtlas.h:349
std::deque< IndexT > IndexDeque
Definition: ParticleAtlas.h:507
Definition: Exceptions.h:40
tree::ValueAccessor< const TreeType > ConstAccessor
Definition: ParticleAtlas.h:181
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:503
void increment()
Advance iterator to next item.
Definition: ParticleAtlas.h:848
void dequeToArray(const std::deque< T > &d, std::unique_ptr< T[]> &a, size_t &size)
Definition: PointIndexGrid.h:555
size_t size() const
Return the number of point indices in the iterator range.
Definition: ParticleAtlas.h:875
bool operator!=(const Iterator &p) const
Definition: ParticleAtlas.h:239
typename PosType::value_type ScalarType
Definition: ParticleAtlas.h:591
Partition particles and performs range and nearest-neighbor searches.
SharedPtr< ParticleAtlas > Ptr
Definition: ParticleAtlas.h:112
ScalarType maxRadius() const
Definition: ParticleAtlas.h:350
PointIndexLeafNodeType *const *const mNodes
Definition: ParticleAtlas.h:495
SharedPtr< SplittableParticleArray > Ptr
Definition: ParticleAtlas.h:320
SplittableParticleArray(const ParticleArrayT &particles)
Definition: ParticleAtlas.h:327
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:91
void filterVoxel(const Coord &, const IndexT *begin, const IndexT *end)
Definition: ParticleAtlas.h:536
double maxRadius(size_t n) const
Returns maximum particle radius for level n.
Definition: ParticleAtlas.h:149
size_t levels() const
Returns the total number of resolution levels.
Definition: ParticleAtlas.h:206
bool operator==(const Iterator &p) const
Return true if both iterators point to the same element.
Definition: ParticleAtlas.h:238
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
void pointIndexSearch(RangeDeque &rangeList, ConstAccessor &acc, const CoordBBox &bbox)
Definition: PointIndexGrid.h:903
void operator()(const tbb::blocked_range< size_t > &range)
Definition: ParticleAtlas.h:293
std::unique_ptr< ConstAccessor > ConstAccessorPtr
Definition: ParticleAtlas.h:182
BBoxFilter(RangeDeque &ranges, IndexDeque &indices, const BBoxd &bbox, const ParticleArrayType &particles, bool hasUniformRadius=false)
Definition: ParticleAtlas.h:597
const IndexType & operator*() const
Return a const reference to the item to which this iterator is pointing.
Definition: ParticleAtlas.h:216
typename ParticleArrayType::PosType PosType
Definition: ParticleAtlas.h:590
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:89
ComputeExtremas(const ParticleArrayT &particles)
Definition: ParticleAtlas.h:279
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:177
ScalarType minRadius
Definition: ParticleAtlas.h:313
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:498
std::deque< Range > RangeDeque
Definition: ParticleAtlas.h:506
ParticleArrayT const *const particleArray
Definition: ParticleAtlas.h:312
void join(const ComputeExtremas &rhs)
Definition: ParticleAtlas.h:307
typename ParticleArray::PosType PosType
Definition: ParticleAtlas.h:324
bool empty() const
true if the container size is 0, false otherwise.
Definition: ParticleAtlas.h:144
Definition: ParticleAtlas.h:110
typename ParticleArrayType::PosType PosType
Definition: ParticleAtlas.h:502
void reset()
Reset the iterator to point to the first item.
Definition: ParticleAtlas.h:831