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