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