GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/PointIndexGrid.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 315 369 85.4%
Functions: 55 73 75.3%
Branches: 218 320 68.1%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @file PointIndexGrid.h
5 ///
6 /// @brief Space-partitioning acceleration structure for points. Partitions
7 /// the points into voxels to accelerate range and nearest neighbor
8 /// searches.
9 ///
10 /// @note Leaf nodes store a single point-index array and the voxels are only
11 /// integer offsets into that array. The actual points are never stored
12 /// in the acceleration structure, only offsets into an external array.
13 ///
14 /// @author Mihai Alden
15
16 #ifndef OPENVDB_TOOLS_POINT_INDEX_GRID_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_POINT_INDEX_GRID_HAS_BEEN_INCLUDED
18
19 #include "openvdb/thread/Threading.h"
20 #include "PointPartitioner.h"
21
22 #include <openvdb/version.h>
23 #include <openvdb/Exceptions.h>
24 #include <openvdb/Grid.h>
25 #include <openvdb/Types.h>
26 #include <openvdb/math/Transform.h>
27 #include <openvdb/tree/LeafManager.h>
28 #include <openvdb/tree/LeafNode.h>
29 #include <openvdb/tree/Tree.h>
30
31 #include <tbb/blocked_range.h>
32 #include <tbb/parallel_for.h>
33 #include <atomic>
34 #include <algorithm> // for std::min(), std::max()
35 #include <cmath> // for std::sqrt()
36 #include <deque>
37 #include <iostream>
38 #include <type_traits> // for std::is_same
39 #include <utility> // for std::pair
40 #include <vector>
41
42
43 namespace openvdb {
44 OPENVDB_USE_VERSION_NAMESPACE
45 namespace OPENVDB_VERSION_NAME {
46
47 namespace tree {
48 template<Index, typename> struct SameLeafConfig; // forward declaration
49 }
50
51 namespace tools {
52
53 template<typename T, Index Log2Dim> struct PointIndexLeafNode; // forward declaration
54
55 /// Point index tree configured to match the default OpenVDB tree configuration
56 using PointIndexTree = tree::Tree<tree::RootNode<tree::InternalNode<tree::InternalNode
57 <PointIndexLeafNode<PointIndex32, 3>, 4>, 5>>>;
58
59 /// Point index grid
60 using PointIndexGrid = Grid<PointIndexTree>;
61
62
63 ////////////////////////////////////////
64
65
66 /// @interface PointArray
67 /// Expected interface for the PointArray container:
68 /// @code
69 /// template<typename VectorType>
70 /// struct PointArray
71 /// {
72 /// // The type used to represent world-space point positions
73 /// using PosType = VectorType;
74 ///
75 /// // Return the number of points in the array
76 /// size_t size() const;
77 ///
78 /// // Return the world-space position of the nth point in the array.
79 /// void getPos(size_t n, PosType& xyz) const;
80 /// };
81 /// @endcode
82
83
84 ////////////////////////////////////////
85
86
87 /// @brief Partition points into a point index grid to accelerate range and
88 /// nearest-neighbor searches.
89 ///
90 /// @param points world-space point array conforming to the PointArray interface
91 /// @param voxelSize voxel size in world units
92 template<typename GridT, typename PointArrayT>
93 inline typename GridT::Ptr
94 createPointIndexGrid(const PointArrayT& points, double voxelSize);
95
96
97 /// @brief Partition points into a point index grid to accelerate range and
98 /// nearest-neighbor searches.
99 ///
100 /// @param points world-space point array conforming to the PointArray interface
101 /// @param xform world-to-index-space transform
102 template<typename GridT, typename PointArrayT>
103 inline typename GridT::Ptr
104 createPointIndexGrid(const PointArrayT& points, const math::Transform& xform);
105
106
107 /// @brief Return @c true if the given point index grid represents a valid partitioning
108 /// of the given point array.
109 ///
110 /// @param points world-space point array conforming to the PointArray interface
111 /// @param grid point index grid to validate
112 template<typename PointArrayT, typename GridT>
113 inline bool
114 isValidPartition(const PointArrayT& points, const GridT& grid);
115
116
117 /// Repartition the @a points if needed, otherwise return the input @a grid.
118 template<typename GridT, typename PointArrayT>
119 inline typename GridT::ConstPtr
120 getValidPointIndexGrid(const PointArrayT& points, const typename GridT::ConstPtr& grid);
121
122 /// Repartition the @a points if needed, otherwise return the input @a grid.
123 template<typename GridT, typename PointArrayT>
124 inline typename GridT::Ptr
125 getValidPointIndexGrid(const PointArrayT& points, const typename GridT::Ptr& grid);
126
127
128 ////////////////////////////////////////
129
130
131 /// Accelerated range and nearest-neighbor searches for point index grids
132 template<typename TreeType = PointIndexTree>
133 struct PointIndexIterator
134 {
135 using ConstAccessor = tree::ValueAccessor<const TreeType>;
136 using LeafNodeType = typename TreeType::LeafNodeType;
137 using ValueType = typename TreeType::ValueType;
138
139
140 PointIndexIterator();
141 PointIndexIterator(const PointIndexIterator& rhs);
142 PointIndexIterator& operator=(const PointIndexIterator& rhs);
143
144
145 /// @brief Construct an iterator over the indices of the points contained in voxel (i, j, k).
146 /// @param ijk the voxel containing the points over which to iterate
147 /// @param acc an accessor for the grid or tree that holds the point indices
148 PointIndexIterator(const Coord& ijk, ConstAccessor& acc);
149
150
151 /// @brief Construct an iterator over the indices of the points contained in
152 /// the given bounding box.
153 /// @param bbox the bounding box of the voxels containing the points over which to iterate
154 /// @param acc an accessor for the grid or tree that holds the point indices
155 /// @note The range of the @a bbox is inclusive. Thus, a bounding box with
156 /// min = max is not empty but rather encloses a single voxel.
157 PointIndexIterator(const CoordBBox& bbox, ConstAccessor& acc);
158
159
160 /// @brief Clear the iterator and update it with the result of the given voxel query.
161 /// @param ijk the voxel containing the points over which to iterate
162 /// @param acc an accessor for the grid or tree that holds the point indices
163 void searchAndUpdate(const Coord& ijk, ConstAccessor& acc);
164
165
166 /// @brief Clear the iterator and update it with the result of the given voxel region query.
167 /// @param bbox the bounding box of the voxels containing the points over which to iterate
168 /// @param acc an accessor for the grid or tree that holds the point indices
169 /// @note The range of the @a bbox is inclusive. Thus, a bounding box with
170 /// min = max is not empty but rather encloses a single voxel.
171 void searchAndUpdate(const CoordBBox& bbox, ConstAccessor& acc);
172
173
174 /// @brief Clear the iterator and update it with the result of the given
175 /// index-space bounding box query.
176 /// @param bbox index-space bounding box
177 /// @param acc an accessor for the grid or tree that holds the point indices
178 /// @param points world-space point array conforming to the PointArray interface
179 /// @param xform linear, uniform-scale transform (i.e., cubical voxels)
180 template<typename PointArray>
181 void searchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
182 const PointArray& points, const math::Transform& xform);
183
184
185 /// @brief Clear the iterator and update it with the result of the given
186 /// index-space radial query.
187 /// @param center index-space center
188 /// @param radius index-space radius
189 /// @param acc an accessor for the grid or tree that holds the point indices
190 /// @param points world-space point array conforming to the PointArray interface
191 /// @param xform linear, uniform-scale transform (i.e., cubical voxels)
192 /// @param subvoxelAccuracy if true, check individual points against the search region,
193 /// otherwise return all points that reside in voxels that are inside
194 /// or intersect the search region
195 template<typename PointArray>
196 void searchAndUpdate(const Vec3d& center, double radius, ConstAccessor& acc,
197 const PointArray& points, const math::Transform& xform, bool subvoxelAccuracy = true);
198
199
200 /// @brief Clear the iterator and update it with the result of the given
201 /// world-space bounding box query.
202 /// @param bbox world-space bounding box
203 /// @param acc an accessor for the grid or tree that holds the point indices
204 /// @param points world-space point array conforming to the PointArray interface
205 /// @param xform linear, uniform-scale transform (i.e., cubical voxels)
206 template<typename PointArray>
207 void worldSpaceSearchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
208 const PointArray& points, const math::Transform& xform);
209
210
211 /// @brief Clear the iterator and update it with the result of the given
212 /// world-space radial query.
213 /// @param center world-space center
214 /// @param radius world-space radius
215 /// @param acc an accessor for the grid or tree that holds the point indices
216 /// @param points world-space point array conforming to the PointArray interface
217 /// @param xform linear, uniform-scale transform (i.e., cubical voxels)
218 /// @param subvoxelAccuracy if true, check individual points against the search region,
219 /// otherwise return all points that reside in voxels that are inside
220 /// or intersect the search region
221 template<typename PointArray>
222 void worldSpaceSearchAndUpdate(const Vec3d& center, double radius, ConstAccessor& acc,
223 const PointArray& points, const math::Transform& xform, bool subvoxelAccuracy = true);
224
225
226 /// Reset the iterator to point to the first item.
227 void reset();
228
229 /// Return a const reference to the item to which this iterator is pointing.
230 const ValueType& operator*() const { return *mRange.first; }
231
232 /// @{
233 /// @brief Return @c true if this iterator is not yet exhausted.
234
19/36
✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 97694 times.
✓ Branch 2 taken 10000 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 124979 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 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 40000 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 40000 times.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 1 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 2212 times.
✓ Branch 34 taken 1 times.
✗ Branch 35 not taken.
 314893 bool test() const { return mRange.first < mRange.second || mIter != mRangeList.end(); }
235 operator bool() const { return this->test(); }
236 /// @}
237
238 /// Advance iterator to next item.
239 void increment();
240
241 /// Advance iterator to next item.
242 207191 void operator++() { this->increment(); }
243
244
245 /// @brief Advance iterator to next item.
246 /// @return @c true if this iterator is not yet exhausted.
247 bool next();
248
249 /// Return the number of point indices in the iterator range.
250 size_t size() const;
251
252 /// Return @c true if both iterators point to the same element.
253 bool operator==(const PointIndexIterator& p) const { return mRange.first == p.mRange.first; }
254 bool operator!=(const PointIndexIterator& p) const { return !this->operator==(p); }
255
256
257 private:
258 using Range = std::pair<const ValueType*, const ValueType*>;
259 using RangeDeque = std::deque<Range>;
260 using RangeDequeCIter = typename RangeDeque::const_iterator;
261 using IndexArray = std::unique_ptr<ValueType[]>;
262
263 void clear();
264
265 // Primary index collection
266 Range mRange;
267 RangeDeque mRangeList;
268 RangeDequeCIter mIter;
269 // Secondary index collection
270 IndexArray mIndexArray;
271 size_t mIndexArraySize;
272 }; // struct PointIndexIterator
273
274
275 /// @brief Selectively extract and filter point data using a custom filter operator.
276 ///
277 /// @par FilterType example:
278 /// @interface FilterType
279 /// @code
280 /// template<typename T>
281 /// struct WeightedAverageAccumulator {
282 /// using ValueType = T;
283 ///
284 /// WeightedAverageAccumulator(T const * const array, const T radius)
285 /// : mValues(array), mInvRadius(1.0/radius), mWeightSum(0.0), mValueSum(0.0) {}
286 ///
287 /// void reset() { mWeightSum = mValueSum = T(0.0); }
288 ///
289 /// // the following method is invoked by the PointIndexFilter
290 /// void operator()(const T distSqr, const size_t pointIndex) {
291 /// const T weight = T(1.0) - openvdb::math::Sqrt(distSqr) * mInvRadius;
292 /// mWeightSum += weight;
293 /// mValueSum += weight * mValues[pointIndex];
294 /// }
295 ///
296 /// T result() const { return mWeightSum > T(0.0) ? mValueSum / mWeightSum : T(0.0); }
297 ///
298 /// private:
299 /// T const * const mValues;
300 /// const T mInvRadius;
301 /// T mWeightSum, mValueSum;
302 /// }; // struct WeightedAverageAccumulator
303 /// @endcode
304 template<typename PointArray, typename TreeType = PointIndexTree>
305 struct PointIndexFilter
306 {
307 using PosType = typename PointArray::PosType;
308 using ScalarType = typename PosType::value_type;
309 using ConstAccessor = tree::ValueAccessor<const TreeType>;
310
311 /// @brief Constructor
312 /// @param points world-space point array conforming to the PointArray interface
313 /// @param tree a point index tree
314 /// @param xform linear, uniform-scale transform (i.e., cubical voxels)
315 PointIndexFilter(const PointArray& points, const TreeType& tree, const math::Transform& xform);
316
317 /// Thread safe copy constructor
318 PointIndexFilter(const PointIndexFilter& rhs);
319
320 /// @brief Perform a radial search query and apply the given filter
321 /// operator to the selected points.
322 /// @param center world-space center
323 /// @param radius world-space radius
324 /// @param op custom filter operator (see the FilterType example for interface details)
325 template<typename FilterType>
326 void searchAndApply(const PosType& center, ScalarType radius, FilterType& op);
327
328 private:
329 PointArray const * const mPoints;
330 ConstAccessor mAcc;
331 const math::Transform mXform;
332 const ScalarType mInvVoxelSize;
333 PointIndexIterator<TreeType> mIter;
334 }; // struct PointIndexFilter
335
336
337 ////////////////////////////////////////
338
339 // Internal operators and implementation details
340
341 /// @cond OPENVDB_DOCS_INTERNAL
342
343 namespace point_index_grid_internal {
344
345 template<typename PointArrayT>
346 struct ValidPartitioningOp
347 {
348 4 ValidPartitioningOp(std::atomic<bool>& hasChanged,
349 const PointArrayT& points, const math::Transform& xform)
350 : mPoints(&points)
351 , mTransform(&xform)
352 4 , mHasChanged(&hasChanged)
353 {
354 }
355
356 template <typename LeafT>
357 6876 void operator()(LeafT &leaf, size_t /*leafIndex*/) const
358 {
359
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 6820 times.
 6876 if ((*mHasChanged)) {
360 56 thread::cancelGroupExecution();
361 56 return;
362 }
363
364 using IndexArrayT = typename LeafT::IndexArray;
365 using IndexT = typename IndexArrayT::value_type;
366 using PosType = typename PointArrayT::PosType;
367
368 typename LeafT::ValueOnCIter iter;
369 Coord voxelCoord;
370 PosType point;
371
372 const IndexT
373 6820 *begin = static_cast<IndexT*>(nullptr),
374 6820 *end = static_cast<IndexT*>(nullptr);
375
376
2/2
✓ Branch 0 taken 88736 times.
✓ Branch 1 taken 6816 times.
 102372 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
377
378
2/2
✓ Branch 0 taken 88732 times.
✓ Branch 1 taken 4 times.
 88736 if ((*mHasChanged)) break;
379
380 88732 voxelCoord = iter.getCoord();
381 88732 leaf.getIndices(iter.pos(), begin, end);
382
383
2/2
✓ Branch 0 taken 101243 times.
✓ Branch 1 taken 88730 times.
 189973 while (begin < end) {
384
385 101243 mPoints->getPos(*begin, point);
386 101245 if (voxelCoord != mTransform->worldToIndexCellCentered(point)) {
387 2 mHasChanged->store(true);
388 break;
389 }
390
391 101241 ++begin;
392 }
393 }
394 }
395
396 private:
397 PointArrayT const * const mPoints;
398 math::Transform const * const mTransform;
399 std::atomic<bool> * const mHasChanged;
400 };
401
402
403 template<typename LeafNodeT>
404 struct PopulateLeafNodesOp
405 {
406 using IndexT = uint32_t;
407 using Partitioner = PointPartitioner<IndexT, LeafNodeT::LOG2DIM>;
408
409 719 PopulateLeafNodesOp(std::unique_ptr<LeafNodeT*[]>& leafNodes,
410 const Partitioner& partitioner)
411 : mLeafNodes(leafNodes.get())
412 719 , mPartitioner(&partitioner)
413 {
414 }
415
416 3620 void operator()(const tbb::blocked_range<size_t>& range) const {
417
418 using VoxelOffsetT = typename Partitioner::VoxelOffsetType;
419
420 3620 size_t maxPointCount = 0;
421
2/2
✓ Branch 0 taken 43399 times.
✓ Branch 1 taken 3620 times.
 47019 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
422
2/2
✓ Branch 1 taken 6951 times.
✓ Branch 2 taken 36448 times.
 50350 maxPointCount = std::max(maxPointCount, mPartitioner->indices(n).size());
423 }
424
425 const IndexT voxelCount = LeafNodeT::SIZE;
426
427 // allocate histogram buffers
428
1/2
✓ Branch 0 taken 3620 times.
✗ Branch 1 not taken.
 3620 std::unique_ptr<VoxelOffsetT[]> offsets{new VoxelOffsetT[maxPointCount]};
429
1/2
✓ Branch 1 taken 3620 times.
✗ Branch 2 not taken.
 3620 std::unique_ptr<IndexT[]> histogram{new IndexT[voxelCount]};
430
431 3620 VoxelOffsetT const * const voxelOffsets = mPartitioner->voxelOffsets().get();
432
433
2/2
✓ Branch 0 taken 43399 times.
✓ Branch 1 taken 3620 times.
 47019 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
434
435
1/2
✓ Branch 1 taken 43399 times.
✗ Branch 2 not taken.
 43399 LeafNodeT* node = new LeafNodeT();
436 43399 node->setOrigin(mPartitioner->origin(n));
437
438 43399 typename Partitioner::IndexIterator it = mPartitioner->indices(n);
439
440 const size_t pointCount = it.size();
441 IndexT const * const indices = &*it;
442
443 // local copy of voxel offsets.
444
2/2
✓ Branch 0 taken 2653121 times.
✓ Branch 1 taken 43399 times.
 2696520 for (IndexT i = 0; i < pointCount; ++i) {
445 2653121 offsets[i] = voxelOffsets[ indices[i] ];
446 }
447
448 // compute voxel-offset histogram
449 43399 memset(&histogram[0], 0, voxelCount * sizeof(IndexT));
450
2/2
✓ Branch 0 taken 2653121 times.
✓ Branch 1 taken 43399 times.
 2696520 for (IndexT i = 0; i < pointCount; ++i) {
451 2653121 ++histogram[ offsets[i] ];
452 }
453
454 typename LeafNodeT::NodeMaskType& mask = node->getValueMask();
455 typename LeafNodeT::Buffer& buffer = node->buffer();
456
457 // scan histogram (all-prefix-sums)
458 IndexT count = 0, startOffset;
459
2/2
✓ Branch 0 taken 22223872 times.
✓ Branch 1 taken 43399 times.
 22267271 for (int i = 0; i < int(voxelCount); ++i) {
460
2/2
✓ Branch 0 taken 624780 times.
✓ Branch 1 taken 21599092 times.
 22223872 if (histogram[i] > 0) {
461 startOffset = count;
462 624780 count += histogram[i];
463 624780 histogram[i] = startOffset;
464 624780 mask.setOn(i);
465 }
466
1/4
✓ Branch 1 taken 22223872 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 22223872 buffer.setValue(i, count);
467 }
468
469 // allocate point-index array
470 43399 node->indices().resize(pointCount);
471 typename LeafNodeT::ValueType * const orderedIndices = node->indices().data();
472
473 // rank and permute
474
2/2
✓ Branch 0 taken 2653121 times.
✓ Branch 1 taken 43399 times.
 2696520 for (IndexT i = 0; i < pointCount; ++i) {
475 2653121 orderedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
476 }
477
478 43399 mLeafNodes[n] = node;
479 }
480 3620 }
481
482 //////////
483
484 LeafNodeT* * const mLeafNodes;
485 Partitioner const * const mPartitioner;
486 };
487
488
489 /// Construct a @c PointIndexTree
490 template<typename TreeType, typename PointArray>
491 inline void
492 734 constructPointTree(TreeType& tree, const math::Transform& xform, const PointArray& points)
493 {
494 using LeafType = typename TreeType::LeafNodeType;
495
496
1/2
✓ Branch 1 taken 719 times.
✗ Branch 2 not taken.
 734 std::unique_ptr<LeafType*[]> leafNodes;
497 size_t leafNodeCount = 0;
498
499 {
500 // Important: Do not disable the cell-centered transform in the PointPartitioner.
501 // This interpretation is assumed in the PointIndexGrid and all related
502 // search algorithms.
503 734 PointPartitioner<uint32_t, LeafType::LOG2DIM> partitioner;
504
1/2
✓ Branch 1 taken 719 times.
✗ Branch 2 not taken.
 734 partitioner.construct(points, xform, /*voxelOrder=*/false, /*recordVoxelOffsets=*/true);
505
506
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 719 times.
 734 if (!partitioner.usingCellCenteredTransform()) {
507 ✗ OPENVDB_THROW(LookupError, "The PointIndexGrid requires a "
508 "cell-centered transform.");
509 }
510
511 leafNodeCount = partitioner.size();
512
1/2
✓ Branch 1 taken 719 times.
✗ Branch 2 not taken.
 734 leafNodes.reset(new LeafType*[leafNodeCount]);
513
514 const tbb::blocked_range<size_t> range(0, leafNodeCount);
515
1/2
✓ Branch 1 taken 719 times.
✗ Branch 2 not taken.
 734 tbb::parallel_for(range, PopulateLeafNodesOp<LeafType>(leafNodes, partitioner));
516 }
517
518 tree::ValueAccessor<TreeType> acc(tree);
519
2/2
✓ Branch 0 taken 43399 times.
✓ Branch 1 taken 719 times.
 75361 for (size_t n = 0; n < leafNodeCount; ++n) {
520
1/2
✓ Branch 1 taken 43399 times.
✗ Branch 2 not taken.
 74627 acc.addLeaf(leafNodes[n]);
521 }
522 734 }
523
524
525 ////////////////////////////////////////
526
527
528 template<typename T>
529 inline void
530
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 dequeToArray(const std::deque<T>& d, std::unique_ptr<T[]>& a, size_t& size)
531 {
532 6 size = d.size();
533
3/4
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 63442 times.
✓ Branch 5 taken 6 times.
 63448 a.reset(new T[size]);
534 typename std::deque<T>::const_iterator it = d.begin(), itEnd = d.end();
535 T* item = a.get();
536
4/4
✓ Branch 0 taken 63442 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 492 times.
✓ Branch 3 taken 62950 times.
 63448 for ( ; it != itEnd; ++it, ++item) *item = *it;
537 6 }
538
539
540 inline void
541
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 constructExclusiveRegions(std::vector<CoordBBox>& regions,
542 const CoordBBox& bbox, const CoordBBox& ibox)
543 {
544 regions.clear();
545 8 regions.reserve(6);
546 8 Coord cmin = ibox.min();
547 8 Coord cmax = ibox.max();
548
549 // left-face bbox
550 8 regions.push_back(bbox);
551 8 regions.back().max().z() = cmin.z();
552
553 // right-face bbox
554 8 regions.push_back(bbox);
555 8 regions.back().min().z() = cmax.z();
556
557 8 --cmax.z(); // accounting for cell centered bucketing.
558 8 ++cmin.z();
559
560 // front-face bbox
561 8 regions.push_back(bbox);
562 CoordBBox* lastRegion = &regions.back();
563 8 lastRegion->min().z() = cmin.z();
564 8 lastRegion->max().z() = cmax.z();
565 8 lastRegion->max().x() = cmin.x();
566
567 // back-face bbox
568 8 regions.push_back(*lastRegion);
569 lastRegion = &regions.back();
570 8 lastRegion->min().x() = cmax.x();
571 8 lastRegion->max().x() = bbox.max().x();
572
573 8 --cmax.x();
574 8 ++cmin.x();
575
576 // bottom-face bbox
577 8 regions.push_back(*lastRegion);
578 lastRegion = &regions.back();
579 8 lastRegion->min().x() = cmin.x();
580 8 lastRegion->max().x() = cmax.x();
581 8 lastRegion->max().y() = cmin.y();
582
583 // top-face bbox
584 8 regions.push_back(*lastRegion);
585 lastRegion = &regions.back();
586 8 lastRegion->min().y() = cmax.y();
587 8 lastRegion->max().y() = bbox.max().y();
588 8 }
589
590
591 template<typename PointArray, typename IndexT>
592 struct BBoxFilter
593 {
594 using PosType = typename PointArray::PosType;
595 using ScalarType = typename PosType::value_type;
596 using Range = std::pair<const IndexT*, const IndexT*>;
597 using RangeDeque = std::deque<Range>;
598 using IndexDeque = std::deque<IndexT>;
599
600 2 BBoxFilter(RangeDeque& ranges, IndexDeque& indices, const BBoxd& bbox,
601 const PointArray& points, const math::Transform& xform)
602 : mRanges(ranges)
603 , mIndices(indices)
604 , mRegion(bbox)
605 , mPoints(points)
606
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 , mMap(*xform.baseMap())
607 {
608 2 }
609
610 template <typename LeafNodeType>
611 ✗ void filterLeafNode(const LeafNodeType& leaf)
612 {
613 typename LeafNodeType::ValueOnCIter iter;
614 const IndexT
615 ✗ *begin = static_cast<IndexT*>(nullptr),
616 ✗ *end = static_cast<IndexT*>(nullptr);
617 ✗ for (iter = leaf.cbeginValueOn(); iter; ++iter) {
618 ✗ leaf.getIndices(iter.pos(), begin, end);
619 ✗ filterVoxel(iter.getCoord(), begin, end);
620 }
621 }
622
623 3278 void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
624 {
625 PosType vec;
626
627
2/2
✓ Branch 0 taken 40320 times.
✓ Branch 1 taken 3278 times.
 43598 for (; begin < end; ++begin) {
628 40320 mPoints.getPos(*begin, vec);
629
630
2/2
✓ Branch 1 taken 20584 times.
✓ Branch 2 taken 19736 times.
 40320 if (mRegion.isInside(mMap.applyInverseMap(vec))) {
631 20584 mIndices.push_back(*begin);
632 }
633 }
634 3278 }
635
636 private:
637 RangeDeque& mRanges;
638 IndexDeque& mIndices;
639 const BBoxd mRegion;
640 const PointArray& mPoints;
641 const math::MapBase& mMap;
642 };
643
644
645 template<typename PointArray, typename IndexT>
646 struct RadialRangeFilter
647 {
648 using PosType = typename PointArray::PosType;
649 using ScalarType = typename PosType::value_type;
650 using Range = std::pair<const IndexT*, const IndexT*>;
651 using RangeDeque = std::deque<Range>;
652 using IndexDeque = std::deque<IndexT>;
653
654 2 RadialRangeFilter(RangeDeque& ranges, IndexDeque& indices, const Vec3d& xyz, double radius,
655 const PointArray& points, const math::Transform& xform,
656 const double leafNodeDim, const bool subvoxelAccuracy)
657 : mRanges(ranges)
658 , mIndices(indices)
659 , mCenter(xyz)
660 , mWSCenter(xform.indexToWorld(xyz))
661 , mVoxelDist1(ScalarType(0.0))
662 , mVoxelDist2(ScalarType(0.0))
663 , mLeafNodeDist1(ScalarType(0.0))
664 , mLeafNodeDist2(ScalarType(0.0))
665 2 , mWSRadiusSqr(ScalarType(radius * xform.voxelSize()[0]))
666 , mPoints(points)
667 2 , mSubvoxelAccuracy(subvoxelAccuracy)
668 {
669 const ScalarType voxelRadius = ScalarType(std::sqrt(3.0) * 0.5);
670 2 mVoxelDist1 = voxelRadius + ScalarType(radius);
671 2 mVoxelDist1 *= mVoxelDist1;
672
673
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (radius > voxelRadius) {
674 2 mVoxelDist2 = ScalarType(radius) - voxelRadius;
675 2 mVoxelDist2 *= mVoxelDist2;
676 }
677
678 2 const ScalarType leafNodeRadius = ScalarType(leafNodeDim * std::sqrt(3.0) * 0.5);
679 2 mLeafNodeDist1 = leafNodeRadius + ScalarType(radius);
680 2 mLeafNodeDist1 *= mLeafNodeDist1;
681
682
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (radius > leafNodeRadius) {
683 2 mLeafNodeDist2 = ScalarType(radius) - leafNodeRadius;
684 2 mLeafNodeDist2 *= mLeafNodeDist2;
685 }
686
687 2 mWSRadiusSqr *= mWSRadiusSqr;
688 2 }
689
690 template <typename LeafNodeType>
691
2/2
✓ Branch 0 taken 2211 times.
✓ Branch 1 taken 84 times.
 2295 void filterLeafNode(const LeafNodeType& leaf)
692 {
693 {
694 const Coord& ijk = leaf.origin();
695 PosType vec;
696 2295 vec[0] = ScalarType(ijk[0]);
697 2295 vec[1] = ScalarType(ijk[1]);
698
2/2
✓ Branch 0 taken 2211 times.
✓ Branch 1 taken 84 times.
 2295 vec[2] = ScalarType(ijk[2]);
699 vec += ScalarType(LeafNodeType::DIM - 1) * 0.5;
700 vec -= mCenter;
701
702 const ScalarType dist = vec.lengthSqr();
703
2/2
✓ Branch 0 taken 2211 times.
✓ Branch 1 taken 84 times.
 2295 if (dist > mLeafNodeDist1) return;
704
705
3/4
✓ Branch 0 taken 2211 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 78 times.
✓ Branch 3 taken 2133 times.
 2211 if (mLeafNodeDist2 > 0.0 && dist < mLeafNodeDist2) {
706 const IndexT* begin = &leaf.indices().front();
707 78 mRanges.push_back(Range(begin, begin + leaf.indices().size()));
708 return;
709 }
710 }
711
712 typename LeafNodeType::ValueOnCIter iter;
713 const IndexT
714 2133 *begin = static_cast<IndexT*>(nullptr),
715 2133 *end = static_cast<IndexT*>(nullptr);
716
2/2
✓ Branch 0 taken 26608 times.
✓ Branch 1 taken 2133 times.
 30874 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
717 26608 leaf.getIndices(iter.pos(), begin, end);
718 26608 filterVoxel(iter.getCoord(), begin, end);
719 }
720 }
721
722
2/2
✓ Branch 0 taken 35888 times.
✓ Branch 1 taken 904 times.
 36792 void filterVoxel(const Coord& ijk, const IndexT* begin, const IndexT* end)
723 {
724 PosType vec;
725
726 {
727 36792 vec[0] = mCenter[0] - ScalarType(ijk[0]);
728 36792 vec[1] = mCenter[1] - ScalarType(ijk[1]);
729
2/2
✓ Branch 0 taken 35888 times.
✓ Branch 1 taken 904 times.
 36792 vec[2] = mCenter[2] - ScalarType(ijk[2]);
730
731 const ScalarType dist = vec.lengthSqr();
732
2/2
✓ Branch 0 taken 35888 times.
✓ Branch 1 taken 904 times.
 37502 if (dist > mVoxelDist1) return;
733
734
4/6
✓ Branch 0 taken 35888 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 35888 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 710 times.
✓ Branch 5 taken 35178 times.
 35888 if (!mSubvoxelAccuracy || (mVoxelDist2 > 0.0 && dist < mVoxelDist2)) {
735
3/4
✓ Branch 0 taken 710 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 505 times.
✓ Branch 3 taken 205 times.
 1420 if (!mRanges.empty() && mRanges.back().second == begin) {
736 505 mRanges.back().second = end;
737 } else {
738 205 mRanges.push_back(Range(begin, end));
739 }
740 710 return;
741 }
742 }
743
744
745
2/2
✓ Branch 0 taken 40175 times.
✓ Branch 1 taken 35178 times.
 75353 while (begin < end) {
746
2/2
✓ Branch 0 taken 40084 times.
✓ Branch 1 taken 91 times.
 40175 mPoints.getPos(*begin, vec);
747
2/2
✓ Branch 0 taken 40084 times.
✓ Branch 1 taken 91 times.
 40175 vec = mWSCenter - vec;
748
749
2/2
✓ Branch 0 taken 40084 times.
✓ Branch 1 taken 91 times.
 40175 if (vec.lengthSqr() < mWSRadiusSqr) {
750 40084 mIndices.push_back(*begin);
751 }
752 40175 ++begin;
753 }
754 }
755
756 private:
757 RangeDeque& mRanges;
758 IndexDeque& mIndices;
759 const PosType mCenter, mWSCenter;
760 ScalarType mVoxelDist1, mVoxelDist2, mLeafNodeDist1, mLeafNodeDist2, mWSRadiusSqr;
761 const PointArray& mPoints;
762 const bool mSubvoxelAccuracy;
763 }; // struct RadialRangeFilter
764
765
766 ////////////////////////////////////////
767
768
769 template<typename RangeFilterType, typename LeafNodeType>
770 inline void
771 3894 filteredPointIndexSearchVoxels(RangeFilterType& filter,
772 const LeafNodeType& leaf, const Coord& min, const Coord& max)
773 {
774 using PointIndexT = typename LeafNodeType::ValueType;
775 Index xPos(0), yPos(0), pos(0);
776 Coord ijk(0);
777
778 const PointIndexT* dataPtr = &leaf.indices().front();
779 PointIndexT beginOffset, endOffset;
780
781
2/2
✓ Branch 0 taken 12514 times.
✓ Branch 1 taken 1947 times.
 28922 for (ijk[0] = min[0]; ijk[0] <= max[0]; ++ijk[0]) {
782 25028 xPos = (ijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
783
2/2
✓ Branch 0 taken 76021 times.
✓ Branch 1 taken 12514 times.
 177070 for (ijk[1] = min[1]; ijk[1] <= max[1]; ++ijk[1]) {
784 152042 yPos = xPos + ((ijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
785
2/2
✓ Branch 0 taken 372498 times.
✓ Branch 1 taken 76021 times.
 897038 for (ijk[2] = min[2]; ijk[2] <= max[2]; ++ijk[2]) {
786 744996 pos = yPos + (ijk[2] & (LeafNodeType::DIM - 1u));
787
788
2/2
✓ Branch 0 taken 731 times.
✓ Branch 1 taken 371767 times.
 744996 beginOffset = (pos == 0 ? PointIndexT(0) : leaf.getValue(pos - 1));
789
2/2
✓ Branch 1 taken 14791 times.
✓ Branch 2 taken 357707 times.
 744996 endOffset = leaf.getValue(pos);
790
791
2/2
✓ Branch 0 taken 14791 times.
✓ Branch 1 taken 357707 times.
 744996 if (endOffset > beginOffset) {
792 29582 filter.filterVoxel(ijk, dataPtr + beginOffset, dataPtr + endOffset);
793 }
794 }
795 }
796 }
797 3894 }
798
799
800 template<typename RangeFilterType, typename ConstAccessor>
801 inline void
802 96 filteredPointIndexSearch(RangeFilterType& filter, ConstAccessor& acc, const CoordBBox& bbox)
803 {
804 using LeafNodeType = typename ConstAccessor::TreeType::LeafNodeType;
805 Coord ijk(0), ijkMax(0), ijkA(0), ijkB(0);
806 const Coord leafMin = bbox.min() & ~(LeafNodeType::DIM - 1);
807 const Coord leafMax = bbox.max() & ~(LeafNodeType::DIM - 1);
808
809
2/2
✓ Branch 0 taken 504 times.
✓ Branch 1 taken 48 times.
 1104 for (ijk[0] = leafMin[0]; ijk[0] <= leafMax[0]; ijk[0] += LeafNodeType::DIM) {
810
2/2
✓ Branch 0 taken 6576 times.
✓ Branch 1 taken 504 times.
 14160 for (ijk[1] = leafMin[1]; ijk[1] <= leafMax[1]; ijk[1] += LeafNodeType::DIM) {
811
2/2
✓ Branch 0 taken 43956 times.
✓ Branch 1 taken 6576 times.
 101064 for (ijk[2] = leafMin[2]; ijk[2] <= leafMax[2]; ijk[2] += LeafNodeType::DIM) {
812
813
2/2
✓ Branch 1 taken 4242 times.
✓ Branch 2 taken 39714 times.
 87912 if (const LeafNodeType* leaf = acc.probeConstLeaf(ijk)) {
814 8484 ijkMax = ijk;
815 ijkMax.offset(LeafNodeType::DIM - 1);
816
817 // intersect leaf bbox with search region.
818 8484 ijkA = Coord::maxComponent(bbox.min(), ijk);
819 8484 ijkB = Coord::minComponent(bbox.max(), ijkMax);
820
821 if (ijkA != ijk || ijkB != ijkMax) {
822 3894 filteredPointIndexSearchVoxels(filter, *leaf, ijkA, ijkB);
823 } else { // leaf bbox is inside the search region
824 4590 filter.filterLeafNode(*leaf);
825 }
826 }
827 }
828 }
829 }
830 96 }
831
832
833 ////////////////////////////////////////
834
835
836 template<typename RangeDeque, typename LeafNodeType>
837 inline void
838 40175 pointIndexSearchVoxels(RangeDeque& rangeList,
839 const LeafNodeType& leaf, const Coord& min, const Coord& max)
840 {
841 using PointIndexT = typename LeafNodeType::ValueType;
842 using IntT = typename PointIndexT::IntType;
843 using Range = typename RangeDeque::value_type;
844
845 40175 Index xPos(0), pos(0), zStride = Index(max[2] - min[2]);
846 const PointIndexT* dataPtr = &leaf.indices().front();
847 PointIndexT beginOffset(0), endOffset(0),
848 40175 previousOffset(static_cast<IntT>(leaf.indices().size() + 1u));
849 Coord ijk(0);
850
851
2/2
✓ Branch 0 taken 169934 times.
✓ Branch 1 taken 40175 times.
 210109 for (ijk[0] = min[0]; ijk[0] <= max[0]; ++ijk[0]) {
852 169934 xPos = (ijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
853
854
2/2
✓ Branch 0 taken 725879 times.
✓ Branch 1 taken 169934 times.
 895813 for (ijk[1] = min[1]; ijk[1] <= max[1]; ++ijk[1]) {
855
2/2
✓ Branch 0 taken 5149 times.
✓ Branch 1 taken 720730 times.
 725879 pos = xPos + ((ijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
856 725879 pos += (min[2] & (LeafNodeType::DIM - 1u));
857
858
2/2
✓ Branch 0 taken 5149 times.
✓ Branch 1 taken 720730 times.
 725879 beginOffset = (pos == 0 ? PointIndexT(0) : leaf.getValue(pos - 1));
859
2/2
✓ Branch 1 taken 86629 times.
✓ Branch 2 taken 639250 times.
 725879 endOffset = leaf.getValue(pos+zStride);
860
861
2/2
✓ Branch 0 taken 86629 times.
✓ Branch 1 taken 639250 times.
 725879 if (endOffset > beginOffset) {
862
863
2/2
✓ Branch 0 taken 47832 times.
✓ Branch 1 taken 38797 times.
 86629 if (beginOffset == previousOffset) {
864
2/2
✓ Branch 0 taken 382 times.
✓ Branch 1 taken 47450 times.
 48214 rangeList.back().second = dataPtr + endOffset;
865 } else {
866 38797 rangeList.push_back(Range(dataPtr + beginOffset, dataPtr + endOffset));
867 }
868
869 86629 previousOffset = endOffset;
870 }
871 }
872 }
873 40175 }
874
875
876 template<typename RangeDeque, typename ConstAccessor>
877 inline void
878 10009 pointIndexSearch(RangeDeque& rangeList, ConstAccessor& acc, const CoordBBox& bbox)
879 {
880 using LeafNodeType = typename ConstAccessor::TreeType::LeafNodeType;
881 using PointIndexT = typename LeafNodeType::ValueType;
882 using Range = typename RangeDeque::value_type;
883
884 Coord ijk(0), ijkMax(0), ijkA(0), ijkB(0);
885 const Coord leafMin = bbox.min() & ~(LeafNodeType::DIM - 1);
886 const Coord leafMax = bbox.max() & ~(LeafNodeType::DIM - 1);
887
888
2/2
✓ Branch 0 taken 17719 times.
✓ Branch 1 taken 10009 times.
 27728 for (ijk[0] = leafMin[0]; ijk[0] <= leafMax[0]; ijk[0] += LeafNodeType::DIM) {
889
2/2
✓ Branch 0 taken 33129 times.
✓ Branch 1 taken 17719 times.
 50848 for (ijk[1] = leafMin[1]; ijk[1] <= leafMax[1]; ijk[1] += LeafNodeType::DIM) {
890
2/2
✓ Branch 0 taken 105488 times.
✓ Branch 1 taken 33129 times.
 138617 for (ijk[2] = leafMin[2]; ijk[2] <= leafMax[2]; ijk[2] += LeafNodeType::DIM) {
891
892
2/2
✓ Branch 1 taken 46121 times.
✓ Branch 2 taken 59367 times.
 105488 if (const LeafNodeType* leaf = acc.probeConstLeaf(ijk)) {
893 46121 ijkMax = ijk;
894 ijkMax.offset(LeafNodeType::DIM - 1);
895
896 // intersect leaf bbox with search region.
897 46121 ijkA = Coord::maxComponent(bbox.min(), ijk);
898 46121 ijkB = Coord::minComponent(bbox.max(), ijkMax);
899
900 if (ijkA != ijk || ijkB != ijkMax) {
901 40175 pointIndexSearchVoxels(rangeList, *leaf, ijkA, ijkB);
902 } else {
903 // leaf bbox is inside the search region, add all indices.
904 const PointIndexT* begin = &leaf->indices().front();
905 5946 rangeList.push_back(Range(begin, (begin + leaf->indices().size())));
906 }
907 }
908 }
909 }
910 }
911 10009 }
912
913
914 } // namespace point_index_grid_internal
915
916 /// @endcond
917
918 // PointIndexIterator implementation
919
920 template<typename TreeType>
921 inline
922 2 PointIndexIterator<TreeType>::PointIndexIterator()
923 : mRange(static_cast<ValueType*>(nullptr), static_cast<ValueType*>(nullptr))
924 , mRangeList()
925 , mIter(mRangeList.begin())
926 , mIndexArray()
927 2 , mIndexArraySize(0)
928 {
929 2 }
930
931
932 template<typename TreeType>
933 inline
934 PointIndexIterator<TreeType>::PointIndexIterator(const PointIndexIterator& rhs)
935 : mRange(rhs.mRange)
936 , mRangeList(rhs.mRangeList)
937 , mIter(mRangeList.begin())
938 , mIndexArray()
939 , mIndexArraySize(rhs.mIndexArraySize)
940 {
941 if (rhs.mIndexArray) {
942 mIndexArray.reset(new ValueType[mIndexArraySize]);
943 memcpy(mIndexArray.get(), rhs.mIndexArray.get(), mIndexArraySize * sizeof(ValueType));
944 }
945 }
946
947
948 template<typename TreeType>
949 inline PointIndexIterator<TreeType>&
950 PointIndexIterator<TreeType>::operator=(const PointIndexIterator& rhs)
951 {
952 if (&rhs != this) {
953 mRange = rhs.mRange;
954 mRangeList = rhs.mRangeList;
955 mIter = mRangeList.begin();
956 mIndexArray.reset();
957 mIndexArraySize = rhs.mIndexArraySize;
958
959 if (rhs.mIndexArray) {
960 mIndexArray.reset(new ValueType[mIndexArraySize]);
961 memcpy(mIndexArray.get(), rhs.mIndexArray.get(), mIndexArraySize * sizeof(ValueType));
962 }
963 }
964 return *this;
965 }
966
967
968 template<typename TreeType>
969 inline
970 PointIndexIterator<TreeType>::PointIndexIterator(const Coord& ijk, ConstAccessor& acc)
971 : mRange(static_cast<ValueType*>(nullptr), static_cast<ValueType*>(nullptr))
972 , mRangeList()
973 , mIter(mRangeList.begin())
974 , mIndexArray()
975 , mIndexArraySize(0)
976 {
977 const LeafNodeType* leaf = acc.probeConstLeaf(ijk);
978 if (leaf && leaf->getIndices(ijk, mRange.first, mRange.second)) {
979 mRangeList.push_back(mRange);
980 mIter = mRangeList.begin();
981 }
982 }
983
984
985 template<typename TreeType>
986 inline
987 1 PointIndexIterator<TreeType>::PointIndexIterator(const CoordBBox& bbox, ConstAccessor& acc)
988 : mRange(static_cast<ValueType*>(nullptr), static_cast<ValueType*>(nullptr))
989 , mRangeList()
990 , mIter(mRangeList.begin())
991 , mIndexArray()
992 1 , mIndexArraySize(0)
993 {
994
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 point_index_grid_internal::pointIndexSearch(mRangeList, acc, bbox);
995
996
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (!mRangeList.empty()) {
997 1 mIter = mRangeList.begin();
998 mRange = mRangeList.front();
999 }
1000 1 }
1001
1002
1003 template<typename TreeType>
1004 inline void
1005
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 PointIndexIterator<TreeType>::reset()
1006 {
1007
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 mIter = mRangeList.begin();
1008
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 if (!mRangeList.empty()) {
1009 mRange = mRangeList.front();
1010
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (mIndexArray) {
1011 1 mRange.first = mIndexArray.get();
1012 1 mRange.second = mRange.first + mIndexArraySize;
1013 } else {
1014 ✗ mRange.first = static_cast<ValueType*>(nullptr);
1015 ✗ mRange.second = static_cast<ValueType*>(nullptr);
1016 }
1017 4 }
1018
1019
1020 template<typename TreeType>
1021 inline void
1022 304885 PointIndexIterator<TreeType>::increment()
1023 {
1024 304885 ++mRange.first;
1025
4/4
✓ Branch 0 taken 40814 times.
✓ Branch 1 taken 264071 times.
✓ Branch 2 taken 40810 times.
✓ Branch 3 taken 4 times.
 304885 if (mRange.first >= mRange.second && mIter != mRangeList.end()) {
1026 ++mIter;
1027
2/2
✓ Branch 0 taken 30807 times.
✓ Branch 1 taken 10003 times.
 40810 if (mIter != mRangeList.end()) {
1028 mRange = *mIter;
1029
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 10000 times.
 10003 } else if (mIndexArray) {
1030 3 mRange.first = mIndexArray.get();
1031 3 mRange.second = mRange.first + mIndexArraySize;
1032 }
1033 }
1034 304885 }
1035
1036
1037 template<typename TreeType>
1038 inline bool
1039 PointIndexIterator<TreeType>::next()
1040 {
1041 if (!this->test()) return false;
1042 this->increment();
1043 return this->test();
1044 }
1045
1046
1047 template<typename TreeType>
1048 inline size_t
1049 9 PointIndexIterator<TreeType>::size() const
1050 {
1051 size_t count = 0;
1052 typename RangeDeque::const_iterator it = mRangeList.begin();
1053
1054
2/2
✓ Branch 0 taken 11770 times.
✓ Branch 1 taken 9 times.
 11779 for ( ; it != mRangeList.end(); ++it) {
1055
2/2
✓ Branch 0 taken 11405 times.
✓ Branch 1 taken 365 times.
 11770 count += it->second - it->first;
1056 }
1057
1058 9 return count + mIndexArraySize;
1059 }
1060
1061
1062 template<typename TreeType>
1063 inline void
1064 10004 PointIndexIterator<TreeType>::clear()
1065 {
1066 10004 mRange.first = static_cast<ValueType*>(nullptr);
1067 10004 mRange.second = static_cast<ValueType*>(nullptr);
1068 10004 mRangeList.clear();
1069
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 10002 times.
 10004 mIter = mRangeList.end();
1070 mIndexArray.reset();
1071 10004 mIndexArraySize = 0;
1072 10004 }
1073
1074
1075 template<typename TreeType>
1076 inline void
1077 PointIndexIterator<TreeType>::searchAndUpdate(const Coord& ijk, ConstAccessor& acc)
1078 {
1079 this->clear();
1080 const LeafNodeType* leaf = acc.probeConstLeaf(ijk);
1081 if (leaf && leaf->getIndices(ijk, mRange.first, mRange.second)) {
1082 mRangeList.push_back(mRange);
1083 mIter = mRangeList.begin();
1084 }
1085 }
1086
1087
1088 template<typename TreeType>
1089 inline void
1090 10000 PointIndexIterator<TreeType>::searchAndUpdate(const CoordBBox& bbox, ConstAccessor& acc)
1091 {
1092 10000 this->clear();
1093 10000 point_index_grid_internal::pointIndexSearch(mRangeList, acc, bbox);
1094
1095
1/2
✓ Branch 0 taken 10000 times.
✗ Branch 1 not taken.
 10000 if (!mRangeList.empty()) {
1096 10000 mIter = mRangeList.begin();
1097 mRange = mRangeList.front();
1098 }
1099 10000 }
1100
1101
1102 template<typename TreeType>
1103 template<typename PointArray>
1104 inline void
1105 2 PointIndexIterator<TreeType>::searchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
1106 const PointArray& points, const math::Transform& xform)
1107 {
1108 2 this->clear();
1109
1110 std::vector<CoordBBox> searchRegions;
1111 2 CoordBBox region(Coord::round(bbox.min()), Coord::round(bbox.max()));
1112
1113 2 const Coord dim = region.dim();
1114
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 2 const int minExtent = std::min(dim[0], std::min(dim[1], dim[2]));
1115
1116
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (minExtent > 2) {
1117 // collect indices that don't need to be tested
1118
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 CoordBBox ibox = region;
1119 ibox.expand(-1);
1120
1121
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::pointIndexSearch(mRangeList, acc, ibox);
1122
1123 // define regions for the filtered search
1124 ibox.expand(1);
1125
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::constructExclusiveRegions(searchRegions, region, ibox);
1126 } else {
1127 ✗ searchRegions.push_back(region);
1128 }
1129
1130 // filtered search
1131 std::deque<ValueType> filteredIndices;
1132 point_index_grid_internal::BBoxFilter<PointArray, ValueType>
1133 2 filter(mRangeList, filteredIndices, bbox, points, xform);
1134
1135
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 14 for (size_t n = 0, N = searchRegions.size(); n < N; ++n) {
1136
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[n]);
1137 }
1138
1139
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1140
1141 2 this->reset();
1142 2 }
1143
1144
1145 template<typename TreeType>
1146 template<typename PointArray>
1147 inline void
1148 2 PointIndexIterator<TreeType>::searchAndUpdate(const Vec3d& center, double radius,
1149 ConstAccessor& acc, const PointArray& points, const math::Transform& xform,
1150 bool subvoxelAccuracy)
1151 {
1152 2 this->clear();
1153 std::vector<CoordBBox> searchRegions;
1154
1155 // bounding box
1156 CoordBBox bbox(
1157 2 Coord::round(Vec3d(center[0] - radius, center[1] - radius, center[2] - radius)),
1158
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 Coord::round(Vec3d(center[0] + radius, center[1] + radius, center[2] + radius)));
1159 bbox.expand(1);
1160
1161 2 const double iRadius = radius * double(1.0 / std::sqrt(3.0));
1162
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (iRadius > 2.0) {
1163 // inscribed box
1164 CoordBBox ibox(
1165 2 Coord::round(Vec3d(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius)),
1166
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 Coord::round(Vec3d(center[0] + iRadius, center[1] + iRadius, center[2] + iRadius)));
1167 ibox.expand(-1);
1168
1169 // collect indices that don't need to be tested
1170
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::pointIndexSearch(mRangeList, acc, ibox);
1171
1172 ibox.expand(1);
1173
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::constructExclusiveRegions(searchRegions, bbox, ibox);
1174 } else {
1175 ✗ searchRegions.push_back(bbox);
1176 }
1177
1178 // filtered search
1179 std::deque<ValueType> filteredIndices;
1180 const double leafNodeDim = double(TreeType::LeafNodeType::DIM);
1181
1182 using FilterT = point_index_grid_internal::RadialRangeFilter<PointArray, ValueType>;
1183
1184
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FilterT filter(mRangeList, filteredIndices,
1185 center, radius, points, xform, leafNodeDim, subvoxelAccuracy);
1186
1187
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 14 for (size_t n = 0, N = searchRegions.size(); n < N; ++n) {
1188
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[n]);
1189 }
1190
1191
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1192
1193 2 this->reset();
1194 2 }
1195
1196
1197 template<typename TreeType>
1198 template<typename PointArray>
1199 inline void
1200 1 PointIndexIterator<TreeType>::worldSpaceSearchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
1201 const PointArray& points, const math::Transform& xform)
1202 {
1203 1 this->searchAndUpdate(
1204 1 BBoxd(xform.worldToIndex(bbox.min()), xform.worldToIndex(bbox.max())), acc, points, xform);
1205 1 }
1206
1207
1208 template<typename TreeType>
1209 template<typename PointArray>
1210 inline void
1211 ✗ PointIndexIterator<TreeType>::worldSpaceSearchAndUpdate(const Vec3d& center, double radius,
1212 ConstAccessor& acc, const PointArray& points, const math::Transform& xform,
1213 bool subvoxelAccuracy)
1214 {
1215 ✗ this->searchAndUpdate(xform.worldToIndex(center),
1216 ✗ (radius / xform.voxelSize()[0]), acc, points, xform, subvoxelAccuracy);
1217 }
1218
1219
1220 ////////////////////////////////////////
1221
1222 // PointIndexFilter implementation
1223
1224 template<typename PointArray, typename TreeType>
1225 inline
1226 1 PointIndexFilter<PointArray, TreeType>::PointIndexFilter(
1227 const PointArray& points, const TreeType& tree, const math::Transform& xform)
1228
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 : mPoints(&points), mAcc(tree), mXform(xform), mInvVoxelSize(1.0/xform.voxelSize()[0])
1229 {
1230 1 }
1231
1232
1233 template<typename PointArray, typename TreeType>
1234 inline
1235 PointIndexFilter<PointArray, TreeType>::PointIndexFilter(const PointIndexFilter& rhs)
1236 : mPoints(rhs.mPoints)
1237 , mAcc(rhs.mAcc.tree())
1238 , mXform(rhs.mXform)
1239 , mInvVoxelSize(rhs.mInvVoxelSize)
1240 {
1241 }
1242
1243
1244 template<typename PointArray, typename TreeType>
1245 template<typename FilterType>
1246 inline void
1247 10000 PointIndexFilter<PointArray, TreeType>::searchAndApply(
1248 const PosType& center, ScalarType radius, FilterType& op)
1249 {
1250
1/2
✓ Branch 0 taken 10000 times.
✗ Branch 1 not taken.
 10000 if (radius * mInvVoxelSize < ScalarType(8.0)) {
1251 10000 mIter.searchAndUpdate(openvdb::CoordBBox(
1252 10000 mXform.worldToIndexCellCentered(center - radius),
1253 20000 mXform.worldToIndexCellCentered(center + radius)), mAcc);
1254 } else {
1255 ✗ mIter.worldSpaceSearchAndUpdate(
1256 ✗ center, radius, mAcc, *mPoints, mXform, /*subvoxelAccuracy=*/false);
1257 }
1258
1259 10000 const ScalarType radiusSqr = radius * radius;
1260 ScalarType distSqr = 0.0;
1261 PosType pos;
1262 97694 for (; mIter; ++mIter) {
1263
2/2
✓ Branch 0 taken 56648 times.
✓ Branch 1 taken 41046 times.
 97694 mPoints->getPos(*mIter, pos);
1264 pos -= center;
1265 distSqr = pos.lengthSqr();
1266
1267
2/2
✓ Branch 0 taken 56648 times.
✓ Branch 1 taken 41046 times.
 97694 if (distSqr < radiusSqr) {
1268 56648 op(distSqr, *mIter);
1269 }
1270 }
1271 10000 }
1272
1273
1274 ////////////////////////////////////////
1275
1276
1277 template<typename GridT, typename PointArrayT>
1278 inline typename GridT::Ptr
1279 744 createPointIndexGrid(const PointArrayT& points, const math::Transform& xform)
1280 {
1281
1/2
✓ Branch 2 taken 729 times.
✗ Branch 3 not taken.
 744 typename GridT::Ptr grid = GridT::create(typename GridT::ValueType(0));
1282
4/8
✓ Branch 1 taken 729 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 729 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 719 times.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 1488 grid->setTransform(xform.copy());
1283
1284
2/2
✓ Branch 0 taken 719 times.
✓ Branch 1 taken 10 times.
 744 if (points.size() > 0) {
1285
1/2
✓ Branch 1 taken 719 times.
✗ Branch 2 not taken.
 734 point_index_grid_internal::constructPointTree(
1286 grid->tree(), grid->transform(), points);
1287 }
1288
1289 744 return grid;
1290 }
1291
1292
1293 template<typename GridT, typename PointArrayT>
1294 inline typename GridT::Ptr
1295 52 createPointIndexGrid(const PointArrayT& points, double voxelSize)
1296 {
1297 52 math::Transform::Ptr xform = math::Transform::createLinearTransform(voxelSize);
1298
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 104 return createPointIndexGrid<GridT>(points, *xform);
1299 }
1300
1301
1302 template<typename PointArrayT, typename GridT>
1303 inline bool
1304 4 isValidPartition(const PointArrayT& points, const GridT& grid)
1305 {
1306 8 tree::LeafManager<const typename GridT::TreeType> leafs(grid.tree());
1307
1308 size_t pointCount = 0;
1309
2/2
✓ Branch 0 taken 10802 times.
✓ Branch 1 taken 4 times.
 10806 for (size_t n = 0, N = leafs.leafCount(); n < N; ++n) {
1310 10802 pointCount += leafs.leaf(n).indices().size();
1311 }
1312
1313
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (points.size() != pointCount) {
1314 return false;
1315 }
1316
1317 std::atomic<bool> changed;
1318 changed = false;
1319
1320 point_index_grid_internal::ValidPartitioningOp<PointArrayT>
1321 op(changed, points, grid.transform());
1322
1323
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 leafs.foreach(op);
1324
1325 4 return !bool(changed);
1326 }
1327
1328
1329 template<typename GridT, typename PointArrayT>
1330 inline typename GridT::ConstPtr
1331 getValidPointIndexGrid(const PointArrayT& points, const typename GridT::ConstPtr& grid)
1332 {
1333 if (isValidPartition(points, *grid)) {
1334 return grid;
1335 }
1336
1337 return createPointIndexGrid<GridT>(points, grid->transform());
1338 }
1339
1340
1341 template<typename GridT, typename PointArrayT>
1342 inline typename GridT::Ptr
1343 1 getValidPointIndexGrid(const PointArrayT& points, const typename GridT::Ptr& grid)
1344 {
1345
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (isValidPartition(points, *grid)) {
1346 return grid;
1347 }
1348
1349 1 return createPointIndexGrid<GridT>(points, grid->transform());
1350 }
1351
1352
1353 ////////////////////////////////////////
1354
1355
1356 template<typename T, Index Log2Dim>
1357 struct PointIndexLeafNode : public tree::LeafNode<T, Log2Dim>
1358 {
1359 using LeafNodeType = PointIndexLeafNode<T, Log2Dim>;
1360 using Ptr = SharedPtr<PointIndexLeafNode>;
1361
1362 using ValueType = T;
1363 using IndexArray = std::vector<ValueType>;
1364
1365
1366
2/4
✓ Branch 1 taken 43399 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 43400 IndexArray& indices() { return mIndices; }
1367 const IndexArray& indices() const { return mIndices; }
1368
1369 bool getIndices(const Coord& ijk, const ValueType*& begin, const ValueType*& end) const;
1370 bool getIndices(Index offset, const ValueType*& begin, const ValueType*& end) const;
1371
1372 void setOffsetOn(Index offset, const ValueType& val);
1373 void setOffsetOnly(Index offset, const ValueType& val);
1374
1375 bool isEmpty(const CoordBBox& bbox) const;
1376
1377 private:
1378 IndexArray mIndices;
1379
1380 ////////////////////////////////////////
1381
1382 // The following methods had to be copied from the LeafNode class
1383 // to make the derived PointIndexLeafNode class compatible with the tree structure.
1384
1385 public:
1386 using BaseLeaf = tree::LeafNode<T, Log2Dim>;
1387 using NodeMaskType = util::NodeMask<Log2Dim>;
1388
1389 using BaseLeaf::LOG2DIM;
1390 using BaseLeaf::TOTAL;
1391 using BaseLeaf::DIM;
1392 using BaseLeaf::NUM_VALUES;
1393 using BaseLeaf::NUM_VOXELS;
1394 using BaseLeaf::SIZE;
1395 using BaseLeaf::LEVEL;
1396
1397 /// Default constructor
1398 43399 PointIndexLeafNode() : BaseLeaf(), mIndices() {}
1399
1400 explicit
1401 1 PointIndexLeafNode(const Coord& coords, const T& value = zeroVal<T>(), bool active = false)
1402 : BaseLeaf(coords, value, active)
1403
1/6
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 , mIndices()
1404 {
1405 }
1406
1407 ✗ PointIndexLeafNode(PartialCreate, const Coord& coords,
1408 const T& value = zeroVal<T>(), bool active = false)
1409 : BaseLeaf(PartialCreate(), coords, value, active)
1410 ✗ , mIndices()
1411 {
1412 }
1413
1414 /// Deep copy constructor
1415 ✗ PointIndexLeafNode(const PointIndexLeafNode& rhs) : BaseLeaf(rhs), mIndices(rhs.mIndices) {}
1416
1417 /// @brief Return @c true if the given node (which may have a different @c ValueType
1418 /// than this node) has the same active value topology as this node.
1419 template<typename OtherType, Index OtherLog2Dim>
1420 bool hasSameTopology(const PointIndexLeafNode<OtherType, OtherLog2Dim>* other) const {
1421 return BaseLeaf::hasSameTopology(other);
1422 }
1423
1424 /// Check for buffer, state and origin equivalence.
1425 bool operator==(const PointIndexLeafNode& other) const { return BaseLeaf::operator==(other); }
1426
1427 bool operator!=(const PointIndexLeafNode& other) const { return !(other == *this); }
1428
1429 template<MergePolicy Policy> void merge(const PointIndexLeafNode& rhs) {
1430 BaseLeaf::merge<Policy>(rhs);
1431 }
1432 template<MergePolicy Policy> void merge(const ValueType& tileValue, bool tileActive) {
1433 BaseLeaf::template merge<Policy>(tileValue, tileActive);
1434 }
1435
1436 template<MergePolicy Policy>
1437 void merge(const PointIndexLeafNode& other,
1438 const ValueType& /*bg*/, const ValueType& /*otherBG*/)
1439 {
1440 BaseLeaf::template merge<Policy>(other);
1441 }
1442
1443 void addLeaf(PointIndexLeafNode*) {}
1444 template<typename AccessorT>
1445 void addLeafAndCache(PointIndexLeafNode*, AccessorT&) {}
1446
1447 //@{
1448 /// @brief Return a pointer to this node.
1449 PointIndexLeafNode* touchLeaf(const Coord&) { return this; }
1450 template<typename AccessorT>
1451 PointIndexLeafNode* touchLeafAndCache(const Coord&, AccessorT&) { return this; }
1452
1453 template<typename NodeT, typename AccessorT>
1454 NodeT* probeNodeAndCache(const Coord&, AccessorT&)
1455 {
1456 OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
1457 if (!(std::is_same<NodeT, PointIndexLeafNode>::value)) return nullptr;
1458 return reinterpret_cast<NodeT*>(this);
1459 OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
1460 }
1461 PointIndexLeafNode* probeLeaf(const Coord&) { return this; }
1462 template<typename AccessorT>
1463 PointIndexLeafNode* probeLeafAndCache(const Coord&, AccessorT&) { return this; }
1464 //@}
1465
1466 //@{
1467 /// @brief Return a @const pointer to this node.
1468 const PointIndexLeafNode* probeConstLeaf(const Coord&) const { return this; }
1469 template<typename AccessorT>
1470 const PointIndexLeafNode* probeConstLeafAndCache(const Coord&, AccessorT&) const {return this;}
1471 template<typename AccessorT>
1472 const PointIndexLeafNode* probeLeafAndCache(const Coord&, AccessorT&) const { return this; }
1473 const PointIndexLeafNode* probeLeaf(const Coord&) const { return this; }
1474 template<typename NodeT, typename AccessorT>
1475 const NodeT* probeConstNodeAndCache(const Coord&, AccessorT&) const
1476 {
1477 OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
1478 if (!(std::is_same<NodeT, PointIndexLeafNode>::value)) return nullptr;
1479 return reinterpret_cast<const NodeT*>(this);
1480 OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
1481 }
1482 //@}
1483
1484
1485 // I/O methods
1486
1487 void readBuffers(std::istream& is, bool fromHalf = false);
1488 void readBuffers(std::istream& is, const CoordBBox&, bool fromHalf = false);
1489 void writeBuffers(std::ostream& os, bool toHalf = false) const;
1490
1491
1492 Index64 memUsage() const;
1493 Index64 memUsageIfLoaded() const;
1494
1495
1496 ////////////////////////////////////////
1497
1498 // Disable all write methods to avoid unintentional changes
1499 // to the point-array offsets.
1500
1501 ✗ void assertNonmodifiable() {
1502 ✗ assert(false && "Cannot modify voxel values in a PointIndexTree.");
1503 }
1504
1505 void setActiveState(const Coord&, bool) { assertNonmodifiable(); }
1506 void setActiveState(Index, bool) { assertNonmodifiable(); }
1507
1508 void setValueOnly(const Coord&, const ValueType&) { assertNonmodifiable(); }
1509 void setValueOnly(Index, const ValueType&) { assertNonmodifiable(); }
1510
1511 void setValueOff(const Coord&) { assertNonmodifiable(); }
1512 void setValueOff(Index) { assertNonmodifiable(); }
1513
1514 void setValueOff(const Coord&, const ValueType&) { assertNonmodifiable(); }
1515 void setValueOff(Index, const ValueType&) { assertNonmodifiable(); }
1516
1517 void setValueOn(const Coord&) { assertNonmodifiable(); }
1518 void setValueOn(Index) { assertNonmodifiable(); }
1519
1520 void setValueOn(const Coord&, const ValueType&) { assertNonmodifiable(); }
1521 void setValueOn(Index, const ValueType&) { assertNonmodifiable(); }
1522
1523 void setValue(const Coord&, const ValueType&) { assertNonmodifiable(); }
1524
1525 void setValuesOn() { assertNonmodifiable(); }
1526 void setValuesOff() { assertNonmodifiable(); }
1527
1528 template<typename ModifyOp>
1529 void modifyValue(Index, const ModifyOp&) { assertNonmodifiable(); }
1530
1531 template<typename ModifyOp>
1532 void modifyValue(const Coord&, const ModifyOp&) { assertNonmodifiable(); }
1533
1534 template<typename ModifyOp>
1535 void modifyValueAndActiveState(const Coord&, const ModifyOp&) { assertNonmodifiable(); }
1536
1537 ✗ void clip(const CoordBBox&, const ValueType&) { assertNonmodifiable(); }
1538
1539 ✗ void fill(const CoordBBox&, const ValueType&, bool) { assertNonmodifiable(); }
1540 void fill(const ValueType&) {}
1541 void fill(const ValueType&, bool) { assertNonmodifiable(); }
1542
1543 template<typename AccessorT>
1544 void setValueOnlyAndCache(const Coord&, const ValueType&, AccessorT&) {assertNonmodifiable();}
1545
1546 template<typename ModifyOp, typename AccessorT>
1547 void modifyValueAndActiveStateAndCache(const Coord&, const ModifyOp&, AccessorT&) {
1548 assertNonmodifiable();
1549 }
1550
1551 template<typename AccessorT>
1552 void setValueOffAndCache(const Coord&, const ValueType&, AccessorT&) { assertNonmodifiable(); }
1553
1554 template<typename AccessorT>
1555 void setActiveStateAndCache(const Coord&, bool, AccessorT&) { assertNonmodifiable(); }
1556
1557 void resetBackground(const ValueType&, const ValueType&) { assertNonmodifiable(); }
1558
1559 void signedFloodFill(const ValueType&) { assertNonmodifiable(); }
1560 void signedFloodFill(const ValueType&, const ValueType&) { assertNonmodifiable(); }
1561
1562 void negate() { assertNonmodifiable(); }
1563
1564 protected:
1565 using ValueOn = typename BaseLeaf::ValueOn;
1566 using ValueOff = typename BaseLeaf::ValueOff;
1567 using ValueAll = typename BaseLeaf::ValueAll;
1568 using ChildOn = typename BaseLeaf::ChildOn;
1569 using ChildOff = typename BaseLeaf::ChildOff;
1570 using ChildAll = typename BaseLeaf::ChildAll;
1571
1572 using MaskOnIterator = typename NodeMaskType::OnIterator;
1573 using MaskOffIterator = typename NodeMaskType::OffIterator;
1574 using MaskDenseIterator = typename NodeMaskType::DenseIterator;
1575
1576 // During topology-only construction, access is needed
1577 // to protected/private members of other template instances.
1578 template<typename, Index> friend struct PointIndexLeafNode;
1579
1580 friend class tree::IteratorBase<MaskOnIterator, PointIndexLeafNode>;
1581 friend class tree::IteratorBase<MaskOffIterator, PointIndexLeafNode>;
1582 friend class tree::IteratorBase<MaskDenseIterator, PointIndexLeafNode>;
1583
1584 public:
1585 using ValueOnIter = typename BaseLeaf::template ValueIter<
1586 MaskOnIterator, PointIndexLeafNode, const ValueType, ValueOn>;
1587 using ValueOnCIter = typename BaseLeaf::template ValueIter<
1588 MaskOnIterator, const PointIndexLeafNode, const ValueType, ValueOn>;
1589 using ValueOffIter = typename BaseLeaf::template ValueIter<
1590 MaskOffIterator, PointIndexLeafNode, const ValueType, ValueOff>;
1591 using ValueOffCIter = typename BaseLeaf::template ValueIter<
1592 MaskOffIterator,const PointIndexLeafNode,const ValueType, ValueOff>;
1593 using ValueAllIter = typename BaseLeaf::template ValueIter<
1594 MaskDenseIterator, PointIndexLeafNode, const ValueType, ValueAll>;
1595 using ValueAllCIter = typename BaseLeaf::template ValueIter<
1596 MaskDenseIterator,const PointIndexLeafNode,const ValueType, ValueAll>;
1597 using ChildOnIter = typename BaseLeaf::template ChildIter<
1598 MaskOnIterator, PointIndexLeafNode, ChildOn>;
1599 using ChildOnCIter = typename BaseLeaf::template ChildIter<
1600 MaskOnIterator, const PointIndexLeafNode, ChildOn>;
1601 using ChildOffIter = typename BaseLeaf::template ChildIter<
1602 MaskOffIterator, PointIndexLeafNode, ChildOff>;
1603 using ChildOffCIter = typename BaseLeaf::template ChildIter<
1604 MaskOffIterator, const PointIndexLeafNode, ChildOff>;
1605 using ChildAllIter = typename BaseLeaf::template DenseIter<
1606 PointIndexLeafNode, ValueType, ChildAll>;
1607 using ChildAllCIter = typename BaseLeaf::template DenseIter<
1608 const PointIndexLeafNode, const ValueType, ChildAll>;
1609
1610 #define VMASK_ this->getValueMask()
1611
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 39870 ValueOnCIter cbeginValueOn() const { return ValueOnCIter(VMASK_.beginOn(), this); }
1612 ValueOnCIter beginValueOn() const { return ValueOnCIter(VMASK_.beginOn(), this); }
1613 ValueOnIter beginValueOn() { return ValueOnIter(VMASK_.beginOn(), this); }
1614 ValueOffCIter cbeginValueOff() const { return ValueOffCIter(VMASK_.beginOff(), this); }
1615 ValueOffCIter beginValueOff() const { return ValueOffCIter(VMASK_.beginOff(), this); }
1616 ValueOffIter beginValueOff() { return ValueOffIter(VMASK_.beginOff(), this); }
1617 ValueAllCIter cbeginValueAll() const { return ValueAllCIter(VMASK_.beginDense(), this); }
1618 ValueAllCIter beginValueAll() const { return ValueAllCIter(VMASK_.beginDense(), this); }
1619 ValueAllIter beginValueAll() { return ValueAllIter(VMASK_.beginDense(), this); }
1620
1621 ValueOnCIter cendValueOn() const { return ValueOnCIter(VMASK_.endOn(), this); }
1622 ValueOnCIter endValueOn() const { return ValueOnCIter(VMASK_.endOn(), this); }
1623 ValueOnIter endValueOn() { return ValueOnIter(VMASK_.endOn(), this); }
1624 ValueOffCIter cendValueOff() const { return ValueOffCIter(VMASK_.endOff(), this); }
1625 ValueOffCIter endValueOff() const { return ValueOffCIter(VMASK_.endOff(), this); }
1626 ValueOffIter endValueOff() { return ValueOffIter(VMASK_.endOff(), this); }
1627 ValueAllCIter cendValueAll() const { return ValueAllCIter(VMASK_.endDense(), this); }
1628 ValueAllCIter endValueAll() const { return ValueAllCIter(VMASK_.endDense(), this); }
1629 ValueAllIter endValueAll() { return ValueAllIter(VMASK_.endDense(), this); }
1630
1631 2700 ChildOnCIter cbeginChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
1632 ChildOnCIter beginChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
1633 ✗ ChildOnIter beginChildOn() { return ChildOnIter(VMASK_.endOn(), this); }
1634 ChildOffCIter cbeginChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
1635 ChildOffCIter beginChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
1636 ChildOffIter beginChildOff() { return ChildOffIter(VMASK_.endOff(), this); }
1637 ChildAllCIter cbeginChildAll() const { return ChildAllCIter(VMASK_.beginDense(), this); }
1638 ChildAllCIter beginChildAll() const { return ChildAllCIter(VMASK_.beginDense(), this); }
1639 ChildAllIter beginChildAll() { return ChildAllIter(VMASK_.beginDense(), this); }
1640
1641 ChildOnCIter cendChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
1642 ChildOnCIter endChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
1643 ChildOnIter endChildOn() { return ChildOnIter(VMASK_.endOn(), this); }
1644 ChildOffCIter cendChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
1645 ChildOffCIter endChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
1646 ChildOffIter endChildOff() { return ChildOffIter(VMASK_.endOff(), this); }
1647 ChildAllCIter cendChildAll() const { return ChildAllCIter(VMASK_.endDense(), this); }
1648 ChildAllCIter endChildAll() const { return ChildAllCIter(VMASK_.endDense(), this); }
1649 ChildAllIter endChildAll() { return ChildAllIter(VMASK_.endDense(), this); }
1650 #undef VMASK_
1651 }; // struct PointIndexLeafNode
1652
1653
1654 template<typename T, Index Log2Dim>
1655 inline bool
1656 411255 PointIndexLeafNode<T, Log2Dim>::getIndices(const Coord& ijk,
1657 const ValueType*& begin, const ValueType*& end) const
1658 {
1659 411255 return getIndices(LeafNodeType::coordToOffset(ijk), begin, end);
1660 }
1661
1662
1663 template<typename T, Index Log2Dim>
1664 inline bool
1665 526595 PointIndexLeafNode<T, Log2Dim>::getIndices(Index offset,
1666 const ValueType*& begin, const ValueType*& end) const
1667 {
1668
1/2
✓ Branch 0 taken 526595 times.
✗ Branch 1 not taken.
 526595 if (this->isValueMaskOn(offset)) {
1669 const ValueType* dataPtr = &mIndices.front();
1670
2/2
✓ Branch 0 taken 20482 times.
✓ Branch 1 taken 506113 times.
 547077 begin = dataPtr + (offset == 0 ? ValueType(0) : this->buffer()[offset - 1]);
1671 526595 end = dataPtr + this->buffer()[offset];
1672 526595 return true;
1673 }
1674 return false;
1675 }
1676
1677
1678 template<typename T, Index Log2Dim>
1679 inline void
1680 PointIndexLeafNode<T, Log2Dim>::setOffsetOn(Index offset, const ValueType& val)
1681 {
1682 this->buffer().setValue(offset, val);
1683 this->setValueMaskOn(offset);
1684 }
1685
1686
1687 template<typename T, Index Log2Dim>
1688 inline void
1689 PointIndexLeafNode<T, Log2Dim>::setOffsetOnly(Index offset, const ValueType& val)
1690 {
1691 this->buffer().setValue(offset, val);
1692 }
1693
1694
1695 template<typename T, Index Log2Dim>
1696 inline bool
1697 PointIndexLeafNode<T, Log2Dim>::isEmpty(const CoordBBox& bbox) const
1698 {
1699 Index xPos, pos, zStride = Index(bbox.max()[2] - bbox.min()[2]);
1700 Coord ijk;
1701
1702 for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
1703 xPos = (ijk[0] & (DIM - 1u)) << (2 * LOG2DIM);
1704
1705 for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
1706 pos = xPos + ((ijk[1] & (DIM - 1u)) << LOG2DIM);
1707 pos += (bbox.min()[2] & (DIM - 1u));
1708
1709 if (this->buffer()[pos+zStride] > (pos == 0 ? T(0) : this->buffer()[pos - 1])) {
1710 return false;
1711 }
1712 }
1713 }
1714
1715 return true;
1716 }
1717
1718
1719 template<typename T, Index Log2Dim>
1720 inline void
1721 ✗ PointIndexLeafNode<T, Log2Dim>::readBuffers(std::istream& is, bool fromHalf)
1722 {
1723 ✗ BaseLeaf::readBuffers(is, fromHalf);
1724
1725 ✗ Index64 numIndices = Index64(0);
1726 ✗ is.read(reinterpret_cast<char*>(&numIndices), sizeof(Index64));
1727
1728 ✗ mIndices.resize(size_t(numIndices));
1729 ✗ is.read(reinterpret_cast<char*>(mIndices.data()), numIndices * sizeof(T));
1730 }
1731
1732
1733 template<typename T, Index Log2Dim>
1734 inline void
1735 ✗ PointIndexLeafNode<T, Log2Dim>::readBuffers(std::istream& is, const CoordBBox& bbox, bool fromHalf)
1736 {
1737 // Read and clip voxel values.
1738 ✗ BaseLeaf::readBuffers(is, bbox, fromHalf);
1739
1740 ✗ Index64 numIndices = Index64(0);
1741 ✗ is.read(reinterpret_cast<char*>(&numIndices), sizeof(Index64));
1742
1743 ✗ const Index64 numBytes = numIndices * sizeof(T);
1744
1745 ✗ if (bbox.hasOverlap(this->getNodeBoundingBox())) {
1746 ✗ mIndices.resize(size_t(numIndices));
1747 ✗ is.read(reinterpret_cast<char*>(mIndices.data()), numBytes);
1748
1749 /// @todo If any voxels were deactivated as a result of clipping in the call to
1750 /// BaseLeaf::readBuffers(), the point index list will need to be regenerated.
1751 } else {
1752 // Read and discard voxel values.
1753 ✗ std::unique_ptr<char[]> buf{new char[numBytes]};
1754 ✗ is.read(buf.get(), numBytes);
1755 }
1756
1757 // Reserved for future use
1758 ✗ Index64 auxDataBytes = Index64(0);
1759 ✗ is.read(reinterpret_cast<char*>(&auxDataBytes), sizeof(Index64));
1760 ✗ if (auxDataBytes > 0) {
1761 // For now, read and discard any auxiliary data.
1762 ✗ std::unique_ptr<char[]> auxData{new char[auxDataBytes]};
1763 ✗ is.read(auxData.get(), auxDataBytes);
1764 }
1765 }
1766
1767
1768 template<typename T, Index Log2Dim>
1769 inline void
1770 ✗ PointIndexLeafNode<T, Log2Dim>::writeBuffers(std::ostream& os, bool toHalf) const
1771 {
1772 ✗ BaseLeaf::writeBuffers(os, toHalf);
1773
1774 ✗ Index64 numIndices = Index64(mIndices.size());
1775 ✗ os.write(reinterpret_cast<const char*>(&numIndices), sizeof(Index64));
1776 ✗ os.write(reinterpret_cast<const char*>(mIndices.data()), numIndices * sizeof(T));
1777
1778 // Reserved for future use
1779 ✗ const Index64 auxDataBytes = Index64(0);
1780 ✗ os.write(reinterpret_cast<const char*>(&auxDataBytes), sizeof(Index64));
1781 }
1782
1783
1784 template<typename T, Index Log2Dim>
1785 inline Index64
1786 PointIndexLeafNode<T, Log2Dim>::memUsage() const
1787 {
1788 ✗ return BaseLeaf::memUsage() + Index64((sizeof(T)*mIndices.capacity()) + sizeof(mIndices));
1789 }
1790
1791 template<typename T, Index Log2Dim>
1792 inline Index64
1793 PointIndexLeafNode<T, Log2Dim>::memUsageIfLoaded() const
1794 {
1795 ✗ return BaseLeaf::memUsageIfLoaded() + Index64((sizeof(T)*mIndices.capacity()) + sizeof(mIndices));
1796 }
1797
1798 } // namespace tools
1799
1800
1801 ////////////////////////////////////////
1802
1803
1804 namespace tree {
1805
1806 /// Helper metafunction used to implement LeafNode::SameConfiguration
1807 /// (which, as an inner class, can't be independently specialized)
1808 template<Index Dim1, typename T2>
1809 struct SameLeafConfig<Dim1, openvdb::tools::PointIndexLeafNode<T2, Dim1> >
1810 {
1811 static const bool value = true;
1812 };
1813
1814 } // namespace tree
1815 } // namespace OPENVDB_VERSION_NAME
1816 } // namespace openvdb
1817
1818 #endif // OPENVDB_TOOLS_POINT_INDEX_GRID_HAS_BEEN_INCLUDED
1819