GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/ParticleAtlas.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 228 261 87.4%
Functions: 21 21 100.0%
Branches: 199 358 55.6%
Line Branch Exec Source
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
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
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 2 double minRadius(size_t n) const { return mMinRadiusArray[n]; }
121 /// @brief Returns maximum particle radius for level @a n.
122
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 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
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
16/32
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 20000 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 20000 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 1408 times.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 1 times.
✓ Branch 29 taken 20559 times.
✓ Branch 30 taken 1 times.
✗ Branch 31 not taken.
 61975 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 61967 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 3 ComputeExtremas(const ParticleArrayT& particles)
254 : particleArray(&particles)
255 , minRadius(std::numeric_limits<ScalarType>::max())
256 3 , maxRadius(-std::numeric_limits<ScalarType>::max())
257 {
258 }
259
260 11 ComputeExtremas(ComputeExtremas& rhs, tbb::split)
261 11 : particleArray(rhs.particleArray)
262 , minRadius(std::numeric_limits<ScalarType>::max())
263 11 , maxRadius(-std::numeric_limits<ScalarType>::max())
264 {
265 }
266
267 974 void operator()(const tbb::blocked_range<size_t>& range) {
268
269 974 ScalarType radius, tmpMin = minRadius, tmpMax = maxRadius;
270
271
2/2
✓ Branch 0 taken 80000 times.
✓ Branch 1 taken 487 times.
 160974 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
272
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40000 times.
 160000 particleArray->getRadius(n, radius);
273
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 80000 times.
 160000 tmpMin = std::min(radius, tmpMin);
274 160000 tmpMax = std::max(radius, tmpMax);
275 }
276
277
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 473 times.
 974 minRadius = std::min(minRadius, tmpMin);
278
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 473 times.
 974 maxRadius = std::max(maxRadius, tmpMax);
279 974 }
280
281 void join(const ComputeExtremas& rhs) {
282
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 11 minRadius = std::min(minRadius, rhs.minRadius);
283
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 2 times.
 11 maxRadius = std::max(maxRadius, rhs.maxRadius);
284 11 }
285
286 ParticleArrayT const * const particleArray;
287 ScalarType minRadius, maxRadius;
288 }; // struct ComputeExtremas
289
290
291 template<typename ParticleArrayT, typename PointIndex>
292 2 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 1 SplittableParticleArray(const ParticleArrayT& particles, double minR, double maxR)
308
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 : mIndexMap(), mParticleArray(&particles), mSize(particles.size())
309 {
310 1 mMinRadius = ScalarType(minR);
311
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mMaxRadius = ScalarType(maxR);
312 }
313
314 1 const ParticleArrayT& particleArray() const { return *mParticleArray; }
315
316
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 17 size_t size() const { return mSize; }
317
318 void getPos(size_t n, PosType& xyz) const
319
1/2
✓ Branch 0 taken 40000 times.
✗ Branch 1 not taken.
 40000 { return mParticleArray->getPos(getGlobalIndex(n), xyz); }
320 void getRadius(size_t n, ScalarType& radius) const
321
1/2
✓ Branch 0 taken 40000 times.
✗ Branch 1 not taken.
 40000 { return mParticleArray->getRadius(getGlobalIndex(n), radius); }
322
323 3 ScalarType minRadius() const { return mMinRadius; }
324 4 ScalarType maxRadius() const { return mMaxRadius; }
325
326
4/8
✓ Branch 0 taken 40000 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 40000 times.
✓ Branch 4 taken 40000 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 40000 times.
✗ Branch 7 not taken.
 160000 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 1 Ptr split(ScalarType maxRadiusLimit) {
331
332
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (mMaxRadius < maxRadiusLimit) return Ptr();
333
334 1 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
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 Ptr output(new SplittableParticleArray(*this, mask));
340
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (output->size() == 0) return Ptr();
341
342 size_t newSize = 0;
343
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 1 times.
 40001 for (size_t n = 0, N = mSize; n < N; ++n) {
344 40000 newSize += size_t(!mask[n]);
345 }
346
347
4/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20000 times.
✓ Branch 7 taken 1 times.
 20001 std::unique_ptr<PointIndex[]> newIndexMap{new PointIndex[newSize]};
348
349 1 setIndexMap(newIndexMap, mask, false);
350
351
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mSize = newSize;
352 mIndexMap.swap(newIndexMap);
353
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 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 1 SplittableParticleArray(const SplittableParticleArray& other,
366 const std::unique_ptr<bool[]>& mask)
367 1 : mIndexMap(), mParticleArray(&other.particleArray()), mSize(0)
368 {
369
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 1 times.
 40001 for (size_t n = 0, N = other.size(); n < N; ++n) {
370 40000 mSize += size_t(mask[n]);
371 }
372
373
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (mSize != 0) {
374
4/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20000 times.
✓ Branch 7 taken 1 times.
 20001 mIndexMap.reset(new PointIndex[mSize]);
375 1 other.setIndexMap(mIndexMap, mask, true);
376 }
377
378
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 updateExtremas();
379 1 }
380
381 struct MaskParticles {
382 1 MaskParticles(const SplittableParticleArray& particles,
383 const std::unique_ptr<bool[]>& mask, ScalarType radius)
384 : particleArray(&particles)
385 , particleMask(mask.get())
386
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 , radiusLimit(radius)
387 {
388 }
389
390 135 void operator()(const tbb::blocked_range<size_t>& range) const {
391 135 const ScalarType maxRadius = radiusLimit;
392 ScalarType radius;
393
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 135 times.
 40135 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
394
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40000 times.
 40000 particleArray->getRadius(n, radius);
395 40000 particleMask[n] = !(radius < maxRadius);
396 }
397 135 }
398
399 SplittableParticleArray const * const particleArray;
400 bool * const particleMask;
401 ScalarType const radiusLimit;
402 }; // struct MaskParticles
403
404 2 inline void updateExtremas() {
405 ComputeExtremas<SplittableParticleArray> op(*this);
406 2 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mSize), op);
407 2 mMinRadius = op.minRadius;
408 2 mMaxRadius = op.maxRadius;
409 2 }
410
411
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 void setIndexMap(std::unique_ptr<PointIndex[]>& newIndexMap,
412 const std::unique_ptr<bool[]>& mask, bool maskValue) const
413 {
414
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 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
2/2
✓ Branch 0 taken 80000 times.
✓ Branch 1 taken 2 times.
 80002 for (size_t idx = 0, n = 0, N = mSize; n < N; ++n) {
421
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 40000 times.
 80000 if (mask[n] == maskValue) {
422 40000 newIndexMap[idx++] = PointIndex(static_cast<typename PointIndex::IntType>(n));
423 }
424 }
425 }
426 2 }
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 2 RemapIndices(const ParticleArrayType& particles, std::vector<PointIndexLeafNodeType*>& nodes)
442 : mParticles(&particles)
443
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 4 , mNodes(nodes.empty() ? nullptr : &nodes.front())
444 {
445 }
446
447 309 void operator()(const tbb::blocked_range<size_t>& range) const
448 {
449 using PointIndexType = typename PointIndexLeafNodeType::ValueType;
450
2/2
✓ Branch 0 taken 1761 times.
✓ Branch 1 taken 309 times.
 2070 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
451
452
1/2
✓ Branch 0 taken 1761 times.
✗ Branch 1 not taken.
 1761 PointIndexLeafNodeType& node = *mNodes[n];
453 const size_t numIndices = node.indices().size();
454
455
1/2
✓ Branch 0 taken 1761 times.
✗ Branch 1 not taken.
 1761 if (numIndices > 0) {
456 PointIndexType* begin = &node.indices().front();
457 const PointIndexType* end = begin + numIndices;
458
459
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 1761 times.
 41761 while (begin < end) {
460
1/2
✓ Branch 0 taken 40000 times.
✗ Branch 1 not taken.
 80000 *begin = PointIndexType(static_cast<typename PointIndexType::IntType>(
461
1/2
✓ Branch 0 taken 40000 times.
✗ Branch 1 not taken.
 40000 mParticles->getGlobalIndex(*begin)));
462 40000 ++begin;
463 }
464 }
465 }
466 309 }
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 2 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 2 , mHasUniformRadius(hasUniformRadius)
491 {
492
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (mHasUniformRadius) {
493 ScalarType uniformRadius;
494 mParticles->getRadius(0, uniformRadius);
495 2 mRadius = mRadius + uniformRadius;
496 2 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 297 void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
511 {
512 PosType pos;
513
514
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 297 times.
 297 if (mHasUniformRadius) {
515
516 297 const ScalarType searchRadiusSqr = mRadius;
517
518
2/2
✓ Branch 0 taken 1040 times.
✓ Branch 1 taken 297 times.
 1337 while (begin < end) {
519
2/2
✓ Branch 0 taken 701 times.
✓ Branch 1 taken 339 times.
 1040 mParticles->getPos(size_t(*begin), pos);
520 1040 const ScalarType distSqr = (mCenter - pos).lengthSqr();
521
2/2
✓ Branch 0 taken 701 times.
✓ Branch 1 taken 339 times.
 1040 if (distSqr < searchRadiusSqr) {
522 701 mIndices.push_back(*begin);
523 }
524 1040 ++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 297 }
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 2 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
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 , mUniformRadiusSqr(ScalarType(0.0))
580 {
581
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (mHasUniformRadius) {
582 mParticles->getRadius(0, mUniformRadiusSqr);
583 2 mUniformRadiusSqr *= mUniformRadiusSqr;
584 }
585 2 }
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 1032 void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
598 {
599 PosType pos;
600
601
1/2
✓ Branch 0 taken 1032 times.
✗ Branch 1 not taken.
 1032 if (mHasUniformRadius) {
602 1032 const ScalarType radiusSqr = mUniformRadiusSqr;
603
604
2/2
✓ Branch 0 taken 2947 times.
✓ Branch 1 taken 1032 times.
 3979 while (begin < end) {
605
606
2/2
✓ Branch 0 taken 1514 times.
✓ Branch 1 taken 1433 times.
 2947 mParticles->getPos(size_t(*begin), pos);
607
2/2
✓ Branch 0 taken 1514 times.
✓ Branch 1 taken 1433 times.
 2947 if (mBBox.isInside(pos)) {
608 1514 mIndices.push_back(*begin++);
609 1514 continue;
610 }
611
612 1433 const ScalarType distSqr = pointToBBoxDistSqr(pos);
613
2/2
✓ Branch 0 taken 559 times.
✓ Branch 1 taken 874 times.
 1433 if (!(distSqr > radiusSqr)) {
614 559 mIndices.push_back(*begin);
615 }
616
617 1433 ++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 1032 }
641
642 private:
643 BBoxFilter(const BBoxFilter&);
644 BBoxFilter& operator=(const BBoxFilter&);
645
646 1433 ScalarType pointToBBoxDistSqr(const PosType& pos) const
647 {
648 ScalarType distSqr = ScalarType(0.0);
649
650
2/2
✓ Branch 0 taken 4299 times.
✓ Branch 1 taken 1433 times.
 5732 for (int i = 0; i < 3; ++i) {
651
652 const ScalarType a = pos[i];
653
654 ScalarType b = mBBox.min()[i];
655
2/2
✓ Branch 0 taken 1433 times.
✓ Branch 1 taken 2866 times.
 4299 if (a < b) {
656 1433 ScalarType delta = b - a;
657 1433 distSqr += delta * delta;
658 }
659
660 b = mBBox.max()[i];
661
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4298 times.
 4299 if (a > b) {
662 1 ScalarType delta = a - b;
663 1 distSqr += delta * delta;
664 }
665 }
666
667 1433 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 1 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 1 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, particles.size()), extremas);
702 1 const double firstMin = extremas.minRadius;
703
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 const double firstMax = extremas.maxRadius;
704 1 const double firstVoxelSize = std::max(minVoxelSize, firstMin);
705
706
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if (!(firstMax < (firstVoxelSize * double(2.0))) && maxLevels > 1) {
707
708 1 std::vector<SplittableParticleArrayPtr> levels;
709
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 levels.push_back(SplittableParticleArrayPtr(
710
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 new SplittableParticleArray(particles, firstMin, firstMax)));
711
712 std::vector<double> voxelSizeArray;
713
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 voxelSizeArray.push_back(firstVoxelSize);
714
715
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 for (size_t n = 0; n < maxLevels; ++n) {
716
717 const double maxParticleRadius = double(levels.back()->maxRadius());
718 2 const double particleRadiusLimit = voxelSizeArray.back() * double(2.0);
719
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (maxParticleRadius < particleRadiusLimit) break;
720
721
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 SplittableParticleArrayPtr newLevel =
722 levels.back()->split(ScalarType(particleRadiusLimit));
723
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (!newLevel) break;
724
725
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 levels.push_back(newLevel);
726
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 1 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
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 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
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mMinRadiusArray.push_back(double(particleArray.minRadius()));
743
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 mMaxRadiusArray.push_back(double(particleArray.maxRadius()));
744
745
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 PointIndexGridPtr grid =
746 createPointIndexGrid<PointIndexGridType>(particleArray, voxelSizeArray[n]);
747
748 nodes.clear();
749 grid->tree().getNodes(nodes);
750
751
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
752 particle_atlas_internal::RemapIndices<SplittableParticleArray,
753 PointIndexLeafNodeType>(particleArray, nodes));
754
755
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mIndexGridArray.push_back(grid);
756 1 }
757
758 } else {
759 ✗ mMinRadiusArray.push_back(firstMin);
760 ✗ mMaxRadiusArray.push_back(firstMax);
761 ✗ mIndexGridArray.push_back(
762 createPointIndexGrid<PointIndexGridType>(particles, firstVoxelSize));
763 }
764 1 }
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 1 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
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 , mAccessorListSize(atlas.levels())
794 {
795
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (mAccessorListSize > 0) {
796
4/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 1 times.
 3 mAccessorList.reset(new ConstAccessorPtr[mAccessorListSize]);
797
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
798
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 mAccessorList[n].reset(new ConstAccessor(atlas.pointIndexGrid(n).tree()));
799 }
800 }
801 1 }
802
803
804 template<typename PointIndexGridType>
805 inline void
806
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 ParticleAtlas<PointIndexGridType>::Iterator::reset()
807 {
808
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 mIter = mRangeList.begin();
809
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 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 4 }
819
820
821 template<typename PointIndexGridType>
822 inline void
823 61967 ParticleAtlas<PointIndexGridType>::Iterator::increment()
824 {
825 61967 ++mRange.first;
826
4/4
✓ Branch 0 taken 3279 times.
✓ Branch 1 taken 58688 times.
✓ Branch 2 taken 3277 times.
✓ Branch 3 taken 2 times.
 61967 if (mRange.first >= mRange.second && mIter != mRangeList.end()) {
827 ++mIter;
828
2/2
✓ Branch 0 taken 3273 times.
✓ Branch 1 taken 4 times.
 3277 if (mIter != mRangeList.end()) {
829 mRange = *mIter;
830
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 } else if (mIndexArray) {
831 2 mRange.first = mIndexArray.get();
832 2 mRange.second = mRange.first + mIndexArraySize;
833 }
834 }
835 61967 }
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 4 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
2/2
✓ Branch 0 taken 3277 times.
✓ Branch 1 taken 4 times.
 3281 for ( ; it != end; ++it) {
857
2/2
✓ Branch 0 taken 3177 times.
✓ Branch 1 taken 100 times.
 3277 count += it->second - it->first;
858 }
859
860 4 return count + mIndexArraySize;
861 }
862
863
864 template<typename PointIndexGridType>
865 inline void
866 4 ParticleAtlas<PointIndexGridType>::Iterator::clear()
867 {
868 4 mRange.first = static_cast<IndexType*>(nullptr);
869 4 mRange.second = static_cast<IndexType*>(nullptr);
870 4 mRangeList.clear();
871
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
 4 mIter = mRangeList.end();
872 mIndexArray.reset();
873 4 mIndexArraySize = 0;
874 4 }
875
876
877 template<typename PointIndexGridType>
878 inline void
879 2 ParticleAtlas<PointIndexGridType>::Iterator::updateFromLevel(size_t level)
880 {
881 using TreeT = typename PointIndexGridType::TreeType;
882 using LeafNodeType = typename TreeType::LeafNodeType;
883
884 2 this->clear();
885
886
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (mAccessorListSize > 0) {
887
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 const size_t levelIdx = std::min(mAccessorListSize - 1, level);
888
889
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const TreeT& tree = mAtlas->pointIndexGrid(levelIdx).tree();
890
891 std::vector<const LeafNodeType*> nodes;
892 tree.getNodes(nodes);
893
894
2/2
✓ Branch 0 taken 1761 times.
✓ Branch 1 taken 2 times.
 1763 for (size_t n = 0, N = nodes.size(); n < N; ++n) {
895
896
1/2
✓ Branch 0 taken 1761 times.
✗ Branch 1 not taken.
 1761 const LeafNodeType& node = *nodes[n];
897 const size_t numIndices = node.indices().size();
898
899
1/2
✓ Branch 0 taken 1761 times.
✗ Branch 1 not taken.
 1761 if (numIndices > 0) {
900 const IndexType* begin = &node.indices().front();
901
1/4
✓ Branch 1 taken 1761 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1761 mRangeList.push_back(Range(begin, (begin + numIndices)));
902 }
903 }
904 }
905
906 2 this->reset();
907 2 }
908
909
910 template<typename PointIndexGridType>
911 template<typename ParticleArrayType>
912 inline void
913 1 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 1 this->clear();
922
923 std::deque<IndexType> filteredIndices;
924 std::vector<CoordBBox> searchRegions;
925
926 1 const double iRadius = radius * double(1.0 / std::sqrt(3.0));
927
928 1 const Vec3d ibMin(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius);
929 1 const Vec3d ibMax(center[0] + iRadius, center[1] + iRadius, center[2] + iRadius);
930
931 1 const Vec3d bMin(center[0] - radius, center[1] - radius, center[2] - radius);
932 1 const Vec3d bMax(center[0] + radius, center[1] + radius, center[2] + radius);
933
934 1 const PosType pos = PosType(center);
935 const ScalarType dist = ScalarType(radius);
936
937
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
938
939
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 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
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xform.worldToIndexCellCentered(ibMin),
948
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 xform.worldToIndexCellCentered(ibMax));
949
950 inscribedRegion.expand(-1); // erode by one voxel
951
952 // collect indices that don't need to be tested
953
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
954
955 searchRegions.clear();
956
957 const openvdb::CoordBBox region(
958
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xform.worldToIndexCellCentered(bMin - maxRadius),
959
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 4 xform.worldToIndexCellCentered(bMax + maxRadius));
960
961 inscribedRegion.expand(1);
962
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 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
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 14 for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
969
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
970 }
971 }
972
973
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
974
975 1 this->reset();
976 1 }
977
978
979 template<typename PointIndexGridType>
980 template<typename ParticleArrayType>
981 inline void
982 1 ParticleAtlas<PointIndexGridType>::Iterator::worldSpaceSearchAndUpdate(
983 const BBoxd& bbox, const ParticleArrayType& particles)
984 {
985 1 this->clear();
986
987 std::deque<IndexType> filteredIndices;
988 std::vector<CoordBBox> searchRegions;
989
990
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (size_t n = 0, N = mAccessorListSize; n < N; ++n) {
991
992
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 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
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xform.worldToIndexCellCentered(bbox.min()),
1000
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 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
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::pointIndexSearch(mRangeList, acc, inscribedRegion);
1006
1007 searchRegions.clear();
1008
1009 const openvdb::CoordBBox region(
1010
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 xform.worldToIndexCellCentered(bbox.min() - maxRadius),
1011
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 4 xform.worldToIndexCellCentered(bbox.max() + maxRadius));
1012
1013 inscribedRegion.expand(1);
1014
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::constructExclusiveRegions(
1015 searchRegions, region, inscribedRegion);
1016
1017 using FilterType = particle_atlas_internal::BBoxFilter<ParticleArrayType, IndexType>;
1018 2 FilterType filter(mRangeList, filteredIndices, bbox, particles, uniformRadius);
1019
1020
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 14 for (size_t i = 0, I = searchRegions.size(); i < I; ++i) {
1021
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[i]);
1022 }
1023 }
1024
1025
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1026
1027 1 this->reset();
1028 1 }
1029
1030
1031 } // namespace tools
1032 } // namespace OPENVDB_VERSION_NAME
1033 } // namespace openvdb
1034
1035 #endif // OPENVDB_TOOLS_PARTICLE_ATLAS_HAS_BEEN_INCLUDED
1036