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