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.
testPointIndexIterator234 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 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 {
ValidPartitioningOpValidPartitioningOp348 ValidPartitioningOp(std::atomic<bool>& hasChanged,
349 const PointArrayT& points, const math::Transform& xform)
350 : mPoints(&points)
351 , mTransform(&xform)
352 , mHasChanged(&hasChanged)
353 {
354 }
355
356 template <typename LeafT>
operatorValidPartitioningOp357 void operator()(LeafT &leaf, size_t /*leafIndex*/) const
358 {
359 if ((*mHasChanged)) {
360 thread::cancelGroupExecution();
361 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 *begin = static_cast<IndexT*>(nullptr),
374 *end = static_cast<IndexT*>(nullptr);
375
376 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
377
378 if ((*mHasChanged)) break;
379
380 voxelCoord = iter.getCoord();
381 leaf.getIndices(iter.pos(), begin, end);
382
383 while (begin < end) {
384
385 mPoints->getPos(*begin, point);
386 if (voxelCoord != mTransform->worldToIndexCellCentered(point)) {
387 mHasChanged->store(true);
388 break;
389 }
390
391 ++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
PopulateLeafNodesOpPopulateLeafNodesOp409 PopulateLeafNodesOp(std::unique_ptr<LeafNodeT*[]>& leafNodes,
410 const Partitioner& partitioner)
411 : mLeafNodes(leafNodes.get())
412 , mPartitioner(&partitioner)
413 {
414 }
415
operatorPopulateLeafNodesOp416 void operator()(const tbb::blocked_range<size_t>& range) const {
417
418 using VoxelOffsetT = typename Partitioner::VoxelOffsetType;
419
420 size_t maxPointCount = 0;
421 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
422 maxPointCount = std::max(maxPointCount, mPartitioner->indices(n).size());
423 }
424
425 const IndexT voxelCount = LeafNodeT::SIZE;
426
427 // allocate histogram buffers
428 std::unique_ptr<VoxelOffsetT[]> offsets{new VoxelOffsetT[maxPointCount]};
429 std::unique_ptr<IndexT[]> histogram{new IndexT[voxelCount]};
430
431 VoxelOffsetT const * const voxelOffsets = mPartitioner->voxelOffsets().get();
432
433 for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
434
435 LeafNodeT* node = new LeafNodeT();
436 node->setOrigin(mPartitioner->origin(n));
437
438 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 for (IndexT i = 0; i < pointCount; ++i) {
445 offsets[i] = voxelOffsets[ indices[i] ];
446 }
447
448 // compute voxel-offset histogram
449 memset(&histogram[0], 0, voxelCount * sizeof(IndexT));
450 for (IndexT i = 0; i < pointCount; ++i) {
451 ++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 for (int i = 0; i < int(voxelCount); ++i) {
460 if (histogram[i] > 0) {
461 startOffset = count;
462 count += histogram[i];
463 histogram[i] = startOffset;
464 mask.setOn(i);
465 }
466 buffer.setValue(i, count);
467 }
468
469 // allocate point-index array
470 node->indices().resize(pointCount);
471 typename LeafNodeT::ValueType * const orderedIndices = node->indices().data();
472
473 // rank and permute
474 for (IndexT i = 0; i < pointCount; ++i) {
475 orderedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
476 }
477
478 mLeafNodes[n] = node;
479 }
480 }
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
constructPointTree(TreeType & tree,const math::Transform & xform,const PointArray & points)492 constructPointTree(TreeType& tree, const math::Transform& xform, const PointArray& points)
493 {
494 using LeafType = typename TreeType::LeafNodeType;
495
496 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 PointPartitioner<uint32_t, LeafType::LOG2DIM> partitioner;
504 partitioner.construct(points, xform, /*voxelOrder=*/false, /*recordVoxelOffsets=*/true);
505
506 if (!partitioner.usingCellCenteredTransform()) {
507 OPENVDB_THROW(LookupError, "The PointIndexGrid requires a "
508 "cell-centered transform.");
509 }
510
511 leafNodeCount = partitioner.size();
512 leafNodes.reset(new LeafType*[leafNodeCount]);
513
514 const tbb::blocked_range<size_t> range(0, leafNodeCount);
515 tbb::parallel_for(range, PopulateLeafNodesOp<LeafType>(leafNodes, partitioner));
516 }
517
518 tree::ValueAccessor<TreeType> acc(tree);
519 for (size_t n = 0; n < leafNodeCount; ++n) {
520 acc.addLeaf(leafNodes[n]);
521 }
522 }
523
524
525 ////////////////////////////////////////
526
527
528 template<typename T>
529 inline void
dequeToArray(const std::deque<T> & d,std::unique_ptr<T[]> & a,size_t & size)530 dequeToArray(const std::deque<T>& d, std::unique_ptr<T[]>& a, size_t& size)
531 {
532 size = d.size();
533 a.reset(new T[size]);
534 typename std::deque<T>::const_iterator it = d.begin(), itEnd = d.end();
535 T* item = a.get();
536 for ( ; it != itEnd; ++it, ++item) *item = *it;
537 }
538
539
540 inline void
constructExclusiveRegions(std::vector<CoordBBox> & regions,const CoordBBox & bbox,const CoordBBox & ibox)541 constructExclusiveRegions(std::vector<CoordBBox>& regions,
542 const CoordBBox& bbox, const CoordBBox& ibox)
543 {
544 regions.clear();
545 regions.reserve(6);
546 Coord cmin = ibox.min();
547 Coord cmax = ibox.max();
548
549 // left-face bbox
550 regions.push_back(bbox);
551 regions.back().max().z() = cmin.z();
552
553 // right-face bbox
554 regions.push_back(bbox);
555 regions.back().min().z() = cmax.z();
556
557 --cmax.z(); // accounting for cell centered bucketing.
558 ++cmin.z();
559
560 // front-face bbox
561 regions.push_back(bbox);
562 CoordBBox* lastRegion = ®ions.back();
563 lastRegion->min().z() = cmin.z();
564 lastRegion->max().z() = cmax.z();
565 lastRegion->max().x() = cmin.x();
566
567 // back-face bbox
568 regions.push_back(*lastRegion);
569 lastRegion = ®ions.back();
570 lastRegion->min().x() = cmax.x();
571 lastRegion->max().x() = bbox.max().x();
572
573 --cmax.x();
574 ++cmin.x();
575
576 // bottom-face bbox
577 regions.push_back(*lastRegion);
578 lastRegion = ®ions.back();
579 lastRegion->min().x() = cmin.x();
580 lastRegion->max().x() = cmax.x();
581 lastRegion->max().y() = cmin.y();
582
583 // top-face bbox
584 regions.push_back(*lastRegion);
585 lastRegion = ®ions.back();
586 lastRegion->min().y() = cmax.y();
587 lastRegion->max().y() = bbox.max().y();
588 }
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
BBoxFilterBBoxFilter600 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 , mMap(*xform.baseMap())
607 {
608 }
609
610 template <typename LeafNodeType>
filterLeafNodeBBoxFilter611 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
filterVoxelBBoxFilter623 void filterVoxel(const Coord&, const IndexT* begin, const IndexT* end)
624 {
625 PosType vec;
626
627 for (; begin < end; ++begin) {
628 mPoints.getPos(*begin, vec);
629
630 if (mRegion.isInside(mMap.applyInverseMap(vec))) {
631 mIndices.push_back(*begin);
632 }
633 }
634 }
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
RadialRangeFilterRadialRangeFilter654 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 , mWSRadiusSqr(ScalarType(radius * xform.voxelSize()[0]))
666 , mPoints(points)
667 , mSubvoxelAccuracy(subvoxelAccuracy)
668 {
669 const ScalarType voxelRadius = ScalarType(std::sqrt(3.0) * 0.5);
670 mVoxelDist1 = voxelRadius + ScalarType(radius);
671 mVoxelDist1 *= mVoxelDist1;
672
673 if (radius > voxelRadius) {
674 mVoxelDist2 = ScalarType(radius) - voxelRadius;
675 mVoxelDist2 *= mVoxelDist2;
676 }
677
678 const ScalarType leafNodeRadius = ScalarType(leafNodeDim * std::sqrt(3.0) * 0.5);
679 mLeafNodeDist1 = leafNodeRadius + ScalarType(radius);
680 mLeafNodeDist1 *= mLeafNodeDist1;
681
682 if (radius > leafNodeRadius) {
683 mLeafNodeDist2 = ScalarType(radius) - leafNodeRadius;
684 mLeafNodeDist2 *= mLeafNodeDist2;
685 }
686
687 mWSRadiusSqr *= mWSRadiusSqr;
688 }
689
690 template <typename LeafNodeType>
filterLeafNodeRadialRangeFilter691 void filterLeafNode(const LeafNodeType& leaf)
692 {
693 {
694 const Coord& ijk = leaf.origin();
695 PosType vec;
696 vec[0] = ScalarType(ijk[0]);
697 vec[1] = ScalarType(ijk[1]);
698 vec[2] = ScalarType(ijk[2]);
699 vec += ScalarType(LeafNodeType::DIM - 1) * 0.5;
700 vec -= mCenter;
701
702 const ScalarType dist = vec.lengthSqr();
703 if (dist > mLeafNodeDist1) return;
704
705 if (mLeafNodeDist2 > 0.0 && dist < mLeafNodeDist2) {
706 const IndexT* begin = &leaf.indices().front();
707 mRanges.push_back(Range(begin, begin + leaf.indices().size()));
708 return;
709 }
710 }
711
712 typename LeafNodeType::ValueOnCIter iter;
713 const IndexT
714 *begin = static_cast<IndexT*>(nullptr),
715 *end = static_cast<IndexT*>(nullptr);
716 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
717 leaf.getIndices(iter.pos(), begin, end);
718 filterVoxel(iter.getCoord(), begin, end);
719 }
720 }
721
filterVoxelRadialRangeFilter722 void filterVoxel(const Coord& ijk, const IndexT* begin, const IndexT* end)
723 {
724 PosType vec;
725
726 {
727 vec[0] = mCenter[0] - ScalarType(ijk[0]);
728 vec[1] = mCenter[1] - ScalarType(ijk[1]);
729 vec[2] = mCenter[2] - ScalarType(ijk[2]);
730
731 const ScalarType dist = vec.lengthSqr();
732 if (dist > mVoxelDist1) return;
733
734 if (!mSubvoxelAccuracy || (mVoxelDist2 > 0.0 && dist < mVoxelDist2)) {
735 if (!mRanges.empty() && mRanges.back().second == begin) {
736 mRanges.back().second = end;
737 } else {
738 mRanges.push_back(Range(begin, end));
739 }
740 return;
741 }
742 }
743
744
745 while (begin < end) {
746 mPoints.getPos(*begin, vec);
747 vec = mWSCenter - vec;
748
749 if (vec.lengthSqr() < mWSRadiusSqr) {
750 mIndices.push_back(*begin);
751 }
752 ++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
filteredPointIndexSearchVoxels(RangeFilterType & filter,const LeafNodeType & leaf,const Coord & min,const Coord & max)771 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 for (ijk[0] = min[0]; ijk[0] <= max[0]; ++ijk[0]) {
782 xPos = (ijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
783 for (ijk[1] = min[1]; ijk[1] <= max[1]; ++ijk[1]) {
784 yPos = xPos + ((ijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
785 for (ijk[2] = min[2]; ijk[2] <= max[2]; ++ijk[2]) {
786 pos = yPos + (ijk[2] & (LeafNodeType::DIM - 1u));
787
788 beginOffset = (pos == 0 ? PointIndexT(0) : leaf.getValue(pos - 1));
789 endOffset = leaf.getValue(pos);
790
791 if (endOffset > beginOffset) {
792 filter.filterVoxel(ijk, dataPtr + beginOffset, dataPtr + endOffset);
793 }
794 }
795 }
796 }
797 }
798
799
800 template<typename RangeFilterType, typename ConstAccessor>
801 inline void
filteredPointIndexSearch(RangeFilterType & filter,ConstAccessor & acc,const CoordBBox & bbox)802 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 for (ijk[0] = leafMin[0]; ijk[0] <= leafMax[0]; ijk[0] += LeafNodeType::DIM) {
810 for (ijk[1] = leafMin[1]; ijk[1] <= leafMax[1]; ijk[1] += LeafNodeType::DIM) {
811 for (ijk[2] = leafMin[2]; ijk[2] <= leafMax[2]; ijk[2] += LeafNodeType::DIM) {
812
813 if (const LeafNodeType* leaf = acc.probeConstLeaf(ijk)) {
814 ijkMax = ijk;
815 ijkMax.offset(LeafNodeType::DIM - 1);
816
817 // intersect leaf bbox with search region.
818 ijkA = Coord::maxComponent(bbox.min(), ijk);
819 ijkB = Coord::minComponent(bbox.max(), ijkMax);
820
821 if (ijkA != ijk || ijkB != ijkMax) {
822 filteredPointIndexSearchVoxels(filter, *leaf, ijkA, ijkB);
823 } else { // leaf bbox is inside the search region
824 filter.filterLeafNode(*leaf);
825 }
826 }
827 }
828 }
829 }
830 }
831
832
833 ////////////////////////////////////////
834
835
836 template<typename RangeDeque, typename LeafNodeType>
837 inline void
pointIndexSearchVoxels(RangeDeque & rangeList,const LeafNodeType & leaf,const Coord & min,const Coord & max)838 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 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 previousOffset(static_cast<IntT>(leaf.indices().size() + 1u));
849 Coord ijk(0);
850
851 for (ijk[0] = min[0]; ijk[0] <= max[0]; ++ijk[0]) {
852 xPos = (ijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
853
854 for (ijk[1] = min[1]; ijk[1] <= max[1]; ++ijk[1]) {
855 pos = xPos + ((ijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
856 pos += (min[2] & (LeafNodeType::DIM - 1u));
857
858 beginOffset = (pos == 0 ? PointIndexT(0) : leaf.getValue(pos - 1));
859 endOffset = leaf.getValue(pos+zStride);
860
861 if (endOffset > beginOffset) {
862
863 if (beginOffset == previousOffset) {
864 rangeList.back().second = dataPtr + endOffset;
865 } else {
866 rangeList.push_back(Range(dataPtr + beginOffset, dataPtr + endOffset));
867 }
868
869 previousOffset = endOffset;
870 }
871 }
872 }
873 }
874
875
876 template<typename RangeDeque, typename ConstAccessor>
877 inline void
pointIndexSearch(RangeDeque & rangeList,ConstAccessor & acc,const CoordBBox & bbox)878 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 for (ijk[0] = leafMin[0]; ijk[0] <= leafMax[0]; ijk[0] += LeafNodeType::DIM) {
889 for (ijk[1] = leafMin[1]; ijk[1] <= leafMax[1]; ijk[1] += LeafNodeType::DIM) {
890 for (ijk[2] = leafMin[2]; ijk[2] <= leafMax[2]; ijk[2] += LeafNodeType::DIM) {
891
892 if (const LeafNodeType* leaf = acc.probeConstLeaf(ijk)) {
893 ijkMax = ijk;
894 ijkMax.offset(LeafNodeType::DIM - 1);
895
896 // intersect leaf bbox with search region.
897 ijkA = Coord::maxComponent(bbox.min(), ijk);
898 ijkB = Coord::minComponent(bbox.max(), ijkMax);
899
900 if (ijkA != ijk || ijkB != ijkMax) {
901 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 rangeList.push_back(Range(begin, (begin + leaf->indices().size())));
906 }
907 }
908 }
909 }
910 }
911 }
912
913
914 } // namespace point_index_grid_internal
915
916 /// @endcond
917
918 // PointIndexIterator implementation
919
920 template<typename TreeType>
921 inline
PointIndexIterator()922 PointIndexIterator<TreeType>::PointIndexIterator()
923 : mRange(static_cast<ValueType*>(nullptr), static_cast<ValueType*>(nullptr))
924 , mRangeList()
925 , mIter(mRangeList.begin())
926 , mIndexArray()
927 , mIndexArraySize(0)
928 {
929 }
930
931
932 template<typename TreeType>
933 inline
PointIndexIterator(const PointIndexIterator & rhs)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
PointIndexIterator(const Coord & ijk,ConstAccessor & acc)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
PointIndexIterator(const CoordBBox & bbox,ConstAccessor & acc)987 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 , mIndexArraySize(0)
993 {
994 point_index_grid_internal::pointIndexSearch(mRangeList, acc, bbox);
995
996 if (!mRangeList.empty()) {
997 mIter = mRangeList.begin();
998 mRange = mRangeList.front();
999 }
1000 }
1001
1002
1003 template<typename TreeType>
1004 inline void
reset()1005 PointIndexIterator<TreeType>::reset()
1006 {
1007 mIter = mRangeList.begin();
1008 if (!mRangeList.empty()) {
1009 mRange = mRangeList.front();
1010 } else if (mIndexArray) {
1011 mRange.first = mIndexArray.get();
1012 mRange.second = mRange.first + mIndexArraySize;
1013 } else {
1014 mRange.first = static_cast<ValueType*>(nullptr);
1015 mRange.second = static_cast<ValueType*>(nullptr);
1016 }
1017 }
1018
1019
1020 template<typename TreeType>
1021 inline void
increment()1022 PointIndexIterator<TreeType>::increment()
1023 {
1024 ++mRange.first;
1025 if (mRange.first >= mRange.second && mIter != mRangeList.end()) {
1026 ++mIter;
1027 if (mIter != mRangeList.end()) {
1028 mRange = *mIter;
1029 } else if (mIndexArray) {
1030 mRange.first = mIndexArray.get();
1031 mRange.second = mRange.first + mIndexArraySize;
1032 }
1033 }
1034 }
1035
1036
1037 template<typename TreeType>
1038 inline bool
next()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
size()1049 PointIndexIterator<TreeType>::size() const
1050 {
1051 size_t count = 0;
1052 typename RangeDeque::const_iterator it = mRangeList.begin();
1053
1054 for ( ; it != mRangeList.end(); ++it) {
1055 count += it->second - it->first;
1056 }
1057
1058 return count + mIndexArraySize;
1059 }
1060
1061
1062 template<typename TreeType>
1063 inline void
clear()1064 PointIndexIterator<TreeType>::clear()
1065 {
1066 mRange.first = static_cast<ValueType*>(nullptr);
1067 mRange.second = static_cast<ValueType*>(nullptr);
1068 mRangeList.clear();
1069 mIter = mRangeList.end();
1070 mIndexArray.reset();
1071 mIndexArraySize = 0;
1072 }
1073
1074
1075 template<typename TreeType>
1076 inline void
searchAndUpdate(const Coord & ijk,ConstAccessor & acc)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
searchAndUpdate(const CoordBBox & bbox,ConstAccessor & acc)1090 PointIndexIterator<TreeType>::searchAndUpdate(const CoordBBox& bbox, ConstAccessor& acc)
1091 {
1092 this->clear();
1093 point_index_grid_internal::pointIndexSearch(mRangeList, acc, bbox);
1094
1095 if (!mRangeList.empty()) {
1096 mIter = mRangeList.begin();
1097 mRange = mRangeList.front();
1098 }
1099 }
1100
1101
1102 template<typename TreeType>
1103 template<typename PointArray>
1104 inline void
searchAndUpdate(const BBoxd & bbox,ConstAccessor & acc,const PointArray & points,const math::Transform & xform)1105 PointIndexIterator<TreeType>::searchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
1106 const PointArray& points, const math::Transform& xform)
1107 {
1108 this->clear();
1109
1110 std::vector<CoordBBox> searchRegions;
1111 CoordBBox region(Coord::round(bbox.min()), Coord::round(bbox.max()));
1112
1113 const Coord dim = region.dim();
1114 const int minExtent = std::min(dim[0], std::min(dim[1], dim[2]));
1115
1116 if (minExtent > 2) {
1117 // collect indices that don't need to be tested
1118 CoordBBox ibox = region;
1119 ibox.expand(-1);
1120
1121 point_index_grid_internal::pointIndexSearch(mRangeList, acc, ibox);
1122
1123 // define regions for the filtered search
1124 ibox.expand(1);
1125 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 filter(mRangeList, filteredIndices, bbox, points, xform);
1134
1135 for (size_t n = 0, N = searchRegions.size(); n < N; ++n) {
1136 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[n]);
1137 }
1138
1139 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1140
1141 this->reset();
1142 }
1143
1144
1145 template<typename TreeType>
1146 template<typename PointArray>
1147 inline void
searchAndUpdate(const Vec3d & center,double radius,ConstAccessor & acc,const PointArray & points,const math::Transform & xform,bool subvoxelAccuracy)1148 PointIndexIterator<TreeType>::searchAndUpdate(const Vec3d& center, double radius,
1149 ConstAccessor& acc, const PointArray& points, const math::Transform& xform,
1150 bool subvoxelAccuracy)
1151 {
1152 this->clear();
1153 std::vector<CoordBBox> searchRegions;
1154
1155 // bounding box
1156 CoordBBox bbox(
1157 Coord::round(Vec3d(center[0] - radius, center[1] - radius, center[2] - radius)),
1158 Coord::round(Vec3d(center[0] + radius, center[1] + radius, center[2] + radius)));
1159 bbox.expand(1);
1160
1161 const double iRadius = radius * double(1.0 / std::sqrt(3.0));
1162 if (iRadius > 2.0) {
1163 // inscribed box
1164 CoordBBox ibox(
1165 Coord::round(Vec3d(center[0] - iRadius, center[1] - iRadius, center[2] - iRadius)),
1166 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 point_index_grid_internal::pointIndexSearch(mRangeList, acc, ibox);
1171
1172 ibox.expand(1);
1173 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 FilterT filter(mRangeList, filteredIndices,
1185 center, radius, points, xform, leafNodeDim, subvoxelAccuracy);
1186
1187 for (size_t n = 0, N = searchRegions.size(); n < N; ++n) {
1188 point_index_grid_internal::filteredPointIndexSearch(filter, acc, searchRegions[n]);
1189 }
1190
1191 point_index_grid_internal::dequeToArray(filteredIndices, mIndexArray, mIndexArraySize);
1192
1193 this->reset();
1194 }
1195
1196
1197 template<typename TreeType>
1198 template<typename PointArray>
1199 inline void
worldSpaceSearchAndUpdate(const BBoxd & bbox,ConstAccessor & acc,const PointArray & points,const math::Transform & xform)1200 PointIndexIterator<TreeType>::worldSpaceSearchAndUpdate(const BBoxd& bbox, ConstAccessor& acc,
1201 const PointArray& points, const math::Transform& xform)
1202 {
1203 this->searchAndUpdate(
1204 BBoxd(xform.worldToIndex(bbox.min()), xform.worldToIndex(bbox.max())), acc, points, xform);
1205 }
1206
1207
1208 template<typename TreeType>
1209 template<typename PointArray>
1210 inline void
worldSpaceSearchAndUpdate(const Vec3d & center,double radius,ConstAccessor & acc,const PointArray & points,const math::Transform & xform,bool subvoxelAccuracy)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
PointIndexFilter(const PointArray & points,const TreeType & tree,const math::Transform & xform)1226 PointIndexFilter<PointArray, TreeType>::PointIndexFilter(
1227 const PointArray& points, const TreeType& tree, const math::Transform& xform)
1228 : mPoints(&points), mAcc(tree), mXform(xform), mInvVoxelSize(1.0/xform.voxelSize()[0])
1229 {
1230 }
1231
1232
1233 template<typename PointArray, typename TreeType>
1234 inline
PointIndexFilter(const PointIndexFilter & rhs)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
searchAndApply(const PosType & center,ScalarType radius,FilterType & op)1247 PointIndexFilter<PointArray, TreeType>::searchAndApply(
1248 const PosType& center, ScalarType radius, FilterType& op)
1249 {
1250 if (radius * mInvVoxelSize < ScalarType(8.0)) {
1251 mIter.searchAndUpdate(openvdb::CoordBBox(
1252 mXform.worldToIndexCellCentered(center - radius),
1253 mXform.worldToIndexCellCentered(center + radius)), mAcc);
1254 } else {
1255 mIter.worldSpaceSearchAndUpdate(
1256 center, radius, mAcc, *mPoints, mXform, /*subvoxelAccuracy=*/false);
1257 }
1258
1259 const ScalarType radiusSqr = radius * radius;
1260 ScalarType distSqr = 0.0;
1261 PosType pos;
1262 for (; mIter; ++mIter) {
1263 mPoints->getPos(*mIter, pos);
1264 pos -= center;
1265 distSqr = pos.lengthSqr();
1266
1267 if (distSqr < radiusSqr) {
1268 op(distSqr, *mIter);
1269 }
1270 }
1271 }
1272
1273
1274 ////////////////////////////////////////
1275
1276
1277 template<typename GridT, typename PointArrayT>
1278 inline typename GridT::Ptr
createPointIndexGrid(const PointArrayT & points,const math::Transform & xform)1279 createPointIndexGrid(const PointArrayT& points, const math::Transform& xform)
1280 {
1281 typename GridT::Ptr grid = GridT::create(typename GridT::ValueType(0));
1282 grid->setTransform(xform.copy());
1283
1284 if (points.size() > 0) {
1285 point_index_grid_internal::constructPointTree(
1286 grid->tree(), grid->transform(), points);
1287 }
1288
1289 return grid;
1290 }
1291
1292
1293 template<typename GridT, typename PointArrayT>
1294 inline typename GridT::Ptr
createPointIndexGrid(const PointArrayT & points,double voxelSize)1295 createPointIndexGrid(const PointArrayT& points, double voxelSize)
1296 {
1297 math::Transform::Ptr xform = math::Transform::createLinearTransform(voxelSize);
1298 return createPointIndexGrid<GridT>(points, *xform);
1299 }
1300
1301
1302 template<typename PointArrayT, typename GridT>
1303 inline bool
isValidPartition(const PointArrayT & points,const GridT & grid)1304 isValidPartition(const PointArrayT& points, const GridT& grid)
1305 {
1306 tree::LeafManager<const typename GridT::TreeType> leafs(grid.tree());
1307
1308 size_t pointCount = 0;
1309 for (size_t n = 0, N = leafs.leafCount(); n < N; ++n) {
1310 pointCount += leafs.leaf(n).indices().size();
1311 }
1312
1313 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 leafs.foreach(op);
1324
1325 return !bool(changed);
1326 }
1327
1328
1329 template<typename GridT, typename PointArrayT>
1330 inline typename GridT::ConstPtr
getValidPointIndexGrid(const PointArrayT & points,const typename GridT::ConstPtr & grid)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
getValidPointIndexGrid(const PointArrayT & points,const typename GridT::Ptr & grid)1343 getValidPointIndexGrid(const PointArrayT& points, const typename GridT::Ptr& grid)
1344 {
1345 if (isValidPartition(points, *grid)) {
1346 return grid;
1347 }
1348
1349 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
indicesPointIndexLeafNode1366 IndexArray& indices() { return mIndices; }
indicesPointIndexLeafNode1367 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
PointIndexLeafNodePointIndexLeafNode1398 PointIndexLeafNode() : BaseLeaf(), mIndices() {}
1399
1400 explicit
1401 PointIndexLeafNode(const Coord& coords, const T& value = zeroVal<T>(), bool active = false)
BaseLeafPointIndexLeafNode1402 : BaseLeaf(coords, value, active)
1403 , mIndices()
1404 {
1405 }
1406
1407 PointIndexLeafNode(PartialCreate, const Coord& coords,
1408 const T& value = zeroVal<T>(), bool active = false)
BaseLeafPointIndexLeafNode1409 : BaseLeaf(PartialCreate(), coords, value, active)
1410 , mIndices()
1411 {
1412 }
1413
1414 /// Deep copy constructor
PointIndexLeafNodePointIndexLeafNode1415 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>
hasSameTopologyPointIndexLeafNode1420 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
mergePointIndexLeafNode1429 template<MergePolicy Policy> void merge(const PointIndexLeafNode& rhs) {
1430 BaseLeaf::merge<Policy>(rhs);
1431 }
mergePointIndexLeafNode1432 template<MergePolicy Policy> void merge(const ValueType& tileValue, bool tileActive) {
1433 BaseLeaf::template merge<Policy>(tileValue, tileActive);
1434 }
1435
1436 template<MergePolicy Policy>
mergePointIndexLeafNode1437 void merge(const PointIndexLeafNode& other,
1438 const ValueType& /*bg*/, const ValueType& /*otherBG*/)
1439 {
1440 BaseLeaf::template merge<Policy>(other);
1441 }
1442
addLeafPointIndexLeafNode1443 void addLeaf(PointIndexLeafNode*) {}
1444 template<typename AccessorT>
addLeafAndCachePointIndexLeafNode1445 void addLeafAndCache(PointIndexLeafNode*, AccessorT&) {}
1446
1447 //@{
1448 /// @brief Return a pointer to this node.
touchLeafPointIndexLeafNode1449 PointIndexLeafNode* touchLeaf(const Coord&) { return this; }
1450 template<typename AccessorT>
touchLeafAndCachePointIndexLeafNode1451 PointIndexLeafNode* touchLeafAndCache(const Coord&, AccessorT&) { return this; }
1452
1453 template<typename NodeT, typename AccessorT>
probeNodeAndCachePointIndexLeafNode1454 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 }
probeLeafPointIndexLeafNode1461 PointIndexLeafNode* probeLeaf(const Coord&) { return this; }
1462 template<typename AccessorT>
probeLeafAndCachePointIndexLeafNode1463 PointIndexLeafNode* probeLeafAndCache(const Coord&, AccessorT&) { return this; }
1464 //@}
1465
1466 //@{
1467 /// @brief Return a @const pointer to this node.
probeConstLeafPointIndexLeafNode1468 const PointIndexLeafNode* probeConstLeaf(const Coord&) const { return this; }
1469 template<typename AccessorT>
probeConstLeafAndCachePointIndexLeafNode1470 const PointIndexLeafNode* probeConstLeafAndCache(const Coord&, AccessorT&) const {return this;}
1471 template<typename AccessorT>
probeLeafAndCachePointIndexLeafNode1472 const PointIndexLeafNode* probeLeafAndCache(const Coord&, AccessorT&) const { return this; }
probeLeafPointIndexLeafNode1473 const PointIndexLeafNode* probeLeaf(const Coord&) const { return this; }
1474 template<typename NodeT, typename AccessorT>
probeConstNodeAndCachePointIndexLeafNode1475 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
1494
1495 ////////////////////////////////////////
1496
1497 // Disable all write methods to avoid unintentional changes
1498 // to the point-array offsets.
1499
assertNonmodifiablePointIndexLeafNode1500 void assertNonmodifiable() {
1501 assert(false && "Cannot modify voxel values in a PointIndexTree.");
1502 }
1503
setActiveStatePointIndexLeafNode1504 void setActiveState(const Coord&, bool) { assertNonmodifiable(); }
setActiveStatePointIndexLeafNode1505 void setActiveState(Index, bool) { assertNonmodifiable(); }
1506
setValueOnlyPointIndexLeafNode1507 void setValueOnly(const Coord&, const ValueType&) { assertNonmodifiable(); }
setValueOnlyPointIndexLeafNode1508 void setValueOnly(Index, const ValueType&) { assertNonmodifiable(); }
1509
setValueOffPointIndexLeafNode1510 void setValueOff(const Coord&) { assertNonmodifiable(); }
setValueOffPointIndexLeafNode1511 void setValueOff(Index) { assertNonmodifiable(); }
1512
setValueOffPointIndexLeafNode1513 void setValueOff(const Coord&, const ValueType&) { assertNonmodifiable(); }
setValueOffPointIndexLeafNode1514 void setValueOff(Index, const ValueType&) { assertNonmodifiable(); }
1515
setValueOnPointIndexLeafNode1516 void setValueOn(const Coord&) { assertNonmodifiable(); }
setValueOnPointIndexLeafNode1517 void setValueOn(Index) { assertNonmodifiable(); }
1518
setValueOnPointIndexLeafNode1519 void setValueOn(const Coord&, const ValueType&) { assertNonmodifiable(); }
setValueOnPointIndexLeafNode1520 void setValueOn(Index, const ValueType&) { assertNonmodifiable(); }
1521
setValuePointIndexLeafNode1522 void setValue(const Coord&, const ValueType&) { assertNonmodifiable(); }
1523
setValuesOnPointIndexLeafNode1524 void setValuesOn() { assertNonmodifiable(); }
setValuesOffPointIndexLeafNode1525 void setValuesOff() { assertNonmodifiable(); }
1526
1527 template<typename ModifyOp>
modifyValuePointIndexLeafNode1528 void modifyValue(Index, const ModifyOp&) { assertNonmodifiable(); }
1529
1530 template<typename ModifyOp>
modifyValuePointIndexLeafNode1531 void modifyValue(const Coord&, const ModifyOp&) { assertNonmodifiable(); }
1532
1533 template<typename ModifyOp>
modifyValueAndActiveStatePointIndexLeafNode1534 void modifyValueAndActiveState(const Coord&, const ModifyOp&) { assertNonmodifiable(); }
1535
clipPointIndexLeafNode1536 void clip(const CoordBBox&, const ValueType&) { assertNonmodifiable(); }
1537
fillPointIndexLeafNode1538 void fill(const CoordBBox&, const ValueType&, bool) { assertNonmodifiable(); }
fillPointIndexLeafNode1539 void fill(const ValueType&) {}
fillPointIndexLeafNode1540 void fill(const ValueType&, bool) { assertNonmodifiable(); }
1541
1542 template<typename AccessorT>
setValueOnlyAndCachePointIndexLeafNode1543 void setValueOnlyAndCache(const Coord&, const ValueType&, AccessorT&) {assertNonmodifiable();}
1544
1545 template<typename ModifyOp, typename AccessorT>
modifyValueAndActiveStateAndCachePointIndexLeafNode1546 void modifyValueAndActiveStateAndCache(const Coord&, const ModifyOp&, AccessorT&) {
1547 assertNonmodifiable();
1548 }
1549
1550 template<typename AccessorT>
setValueOffAndCachePointIndexLeafNode1551 void setValueOffAndCache(const Coord&, const ValueType&, AccessorT&) { assertNonmodifiable(); }
1552
1553 template<typename AccessorT>
setActiveStateAndCachePointIndexLeafNode1554 void setActiveStateAndCache(const Coord&, bool, AccessorT&) { assertNonmodifiable(); }
1555
resetBackgroundPointIndexLeafNode1556 void resetBackground(const ValueType&, const ValueType&) { assertNonmodifiable(); }
1557
signedFloodFillPointIndexLeafNode1558 void signedFloodFill(const ValueType&) { assertNonmodifiable(); }
signedFloodFillPointIndexLeafNode1559 void signedFloodFill(const ValueType&, const ValueType&) { assertNonmodifiable(); }
1560
negatePointIndexLeafNode1561 void negate() { assertNonmodifiable(); }
1562
1563 protected:
1564 using ValueOn = typename BaseLeaf::ValueOn;
1565 using ValueOff = typename BaseLeaf::ValueOff;
1566 using ValueAll = typename BaseLeaf::ValueAll;
1567 using ChildOn = typename BaseLeaf::ChildOn;
1568 using ChildOff = typename BaseLeaf::ChildOff;
1569 using ChildAll = typename BaseLeaf::ChildAll;
1570
1571 using MaskOnIterator = typename NodeMaskType::OnIterator;
1572 using MaskOffIterator = typename NodeMaskType::OffIterator;
1573 using MaskDenseIterator = typename NodeMaskType::DenseIterator;
1574
1575 // During topology-only construction, access is needed
1576 // to protected/private members of other template instances.
1577 template<typename, Index> friend struct PointIndexLeafNode;
1578
1579 friend class tree::IteratorBase<MaskOnIterator, PointIndexLeafNode>;
1580 friend class tree::IteratorBase<MaskOffIterator, PointIndexLeafNode>;
1581 friend class tree::IteratorBase<MaskDenseIterator, PointIndexLeafNode>;
1582
1583 public:
1584 using ValueOnIter = typename BaseLeaf::template ValueIter<
1585 MaskOnIterator, PointIndexLeafNode, const ValueType, ValueOn>;
1586 using ValueOnCIter = typename BaseLeaf::template ValueIter<
1587 MaskOnIterator, const PointIndexLeafNode, const ValueType, ValueOn>;
1588 using ValueOffIter = typename BaseLeaf::template ValueIter<
1589 MaskOffIterator, PointIndexLeafNode, const ValueType, ValueOff>;
1590 using ValueOffCIter = typename BaseLeaf::template ValueIter<
1591 MaskOffIterator,const PointIndexLeafNode,const ValueType, ValueOff>;
1592 using ValueAllIter = typename BaseLeaf::template ValueIter<
1593 MaskDenseIterator, PointIndexLeafNode, const ValueType, ValueAll>;
1594 using ValueAllCIter = typename BaseLeaf::template ValueIter<
1595 MaskDenseIterator,const PointIndexLeafNode,const ValueType, ValueAll>;
1596 using ChildOnIter = typename BaseLeaf::template ChildIter<
1597 MaskOnIterator, PointIndexLeafNode, ChildOn>;
1598 using ChildOnCIter = typename BaseLeaf::template ChildIter<
1599 MaskOnIterator, const PointIndexLeafNode, ChildOn>;
1600 using ChildOffIter = typename BaseLeaf::template ChildIter<
1601 MaskOffIterator, PointIndexLeafNode, ChildOff>;
1602 using ChildOffCIter = typename BaseLeaf::template ChildIter<
1603 MaskOffIterator, const PointIndexLeafNode, ChildOff>;
1604 using ChildAllIter = typename BaseLeaf::template DenseIter<
1605 PointIndexLeafNode, ValueType, ChildAll>;
1606 using ChildAllCIter = typename BaseLeaf::template DenseIter<
1607 const PointIndexLeafNode, const ValueType, ChildAll>;
1608
1609 #define VMASK_ this->getValueMask()
cbeginValueOnPointIndexLeafNode1610 ValueOnCIter cbeginValueOn() const { return ValueOnCIter(VMASK_.beginOn(), this); }
beginValueOnPointIndexLeafNode1611 ValueOnCIter beginValueOn() const { return ValueOnCIter(VMASK_.beginOn(), this); }
beginValueOnPointIndexLeafNode1612 ValueOnIter beginValueOn() { return ValueOnIter(VMASK_.beginOn(), this); }
cbeginValueOffPointIndexLeafNode1613 ValueOffCIter cbeginValueOff() const { return ValueOffCIter(VMASK_.beginOff(), this); }
beginValueOffPointIndexLeafNode1614 ValueOffCIter beginValueOff() const { return ValueOffCIter(VMASK_.beginOff(), this); }
beginValueOffPointIndexLeafNode1615 ValueOffIter beginValueOff() { return ValueOffIter(VMASK_.beginOff(), this); }
cbeginValueAllPointIndexLeafNode1616 ValueAllCIter cbeginValueAll() const { return ValueAllCIter(VMASK_.beginDense(), this); }
beginValueAllPointIndexLeafNode1617 ValueAllCIter beginValueAll() const { return ValueAllCIter(VMASK_.beginDense(), this); }
beginValueAllPointIndexLeafNode1618 ValueAllIter beginValueAll() { return ValueAllIter(VMASK_.beginDense(), this); }
1619
cendValueOnPointIndexLeafNode1620 ValueOnCIter cendValueOn() const { return ValueOnCIter(VMASK_.endOn(), this); }
endValueOnPointIndexLeafNode1621 ValueOnCIter endValueOn() const { return ValueOnCIter(VMASK_.endOn(), this); }
endValueOnPointIndexLeafNode1622 ValueOnIter endValueOn() { return ValueOnIter(VMASK_.endOn(), this); }
cendValueOffPointIndexLeafNode1623 ValueOffCIter cendValueOff() const { return ValueOffCIter(VMASK_.endOff(), this); }
endValueOffPointIndexLeafNode1624 ValueOffCIter endValueOff() const { return ValueOffCIter(VMASK_.endOff(), this); }
endValueOffPointIndexLeafNode1625 ValueOffIter endValueOff() { return ValueOffIter(VMASK_.endOff(), this); }
cendValueAllPointIndexLeafNode1626 ValueAllCIter cendValueAll() const { return ValueAllCIter(VMASK_.endDense(), this); }
endValueAllPointIndexLeafNode1627 ValueAllCIter endValueAll() const { return ValueAllCIter(VMASK_.endDense(), this); }
endValueAllPointIndexLeafNode1628 ValueAllIter endValueAll() { return ValueAllIter(VMASK_.endDense(), this); }
1629
cbeginChildOnPointIndexLeafNode1630 ChildOnCIter cbeginChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
beginChildOnPointIndexLeafNode1631 ChildOnCIter beginChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
beginChildOnPointIndexLeafNode1632 ChildOnIter beginChildOn() { return ChildOnIter(VMASK_.endOn(), this); }
cbeginChildOffPointIndexLeafNode1633 ChildOffCIter cbeginChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
beginChildOffPointIndexLeafNode1634 ChildOffCIter beginChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
beginChildOffPointIndexLeafNode1635 ChildOffIter beginChildOff() { return ChildOffIter(VMASK_.endOff(), this); }
cbeginChildAllPointIndexLeafNode1636 ChildAllCIter cbeginChildAll() const { return ChildAllCIter(VMASK_.beginDense(), this); }
beginChildAllPointIndexLeafNode1637 ChildAllCIter beginChildAll() const { return ChildAllCIter(VMASK_.beginDense(), this); }
beginChildAllPointIndexLeafNode1638 ChildAllIter beginChildAll() { return ChildAllIter(VMASK_.beginDense(), this); }
1639
cendChildOnPointIndexLeafNode1640 ChildOnCIter cendChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
endChildOnPointIndexLeafNode1641 ChildOnCIter endChildOn() const { return ChildOnCIter(VMASK_.endOn(), this); }
endChildOnPointIndexLeafNode1642 ChildOnIter endChildOn() { return ChildOnIter(VMASK_.endOn(), this); }
cendChildOffPointIndexLeafNode1643 ChildOffCIter cendChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
endChildOffPointIndexLeafNode1644 ChildOffCIter endChildOff() const { return ChildOffCIter(VMASK_.endOff(), this); }
endChildOffPointIndexLeafNode1645 ChildOffIter endChildOff() { return ChildOffIter(VMASK_.endOff(), this); }
cendChildAllPointIndexLeafNode1646 ChildAllCIter cendChildAll() const { return ChildAllCIter(VMASK_.endDense(), this); }
endChildAllPointIndexLeafNode1647 ChildAllCIter endChildAll() const { return ChildAllCIter(VMASK_.endDense(), this); }
endChildAllPointIndexLeafNode1648 ChildAllIter endChildAll() { return ChildAllIter(VMASK_.endDense(), this); }
1649 #undef VMASK_
1650 }; // struct PointIndexLeafNode
1651
1652
1653 template<typename T, Index Log2Dim>
1654 inline bool
getIndices(const Coord & ijk,const ValueType * & begin,const ValueType * & end)1655 PointIndexLeafNode<T, Log2Dim>::getIndices(const Coord& ijk,
1656 const ValueType*& begin, const ValueType*& end) const
1657 {
1658 return getIndices(LeafNodeType::coordToOffset(ijk), begin, end);
1659 }
1660
1661
1662 template<typename T, Index Log2Dim>
1663 inline bool
getIndices(Index offset,const ValueType * & begin,const ValueType * & end)1664 PointIndexLeafNode<T, Log2Dim>::getIndices(Index offset,
1665 const ValueType*& begin, const ValueType*& end) const
1666 {
1667 if (this->isValueMaskOn(offset)) {
1668 const ValueType* dataPtr = &mIndices.front();
1669 begin = dataPtr + (offset == 0 ? ValueType(0) : this->buffer()[offset - 1]);
1670 end = dataPtr + this->buffer()[offset];
1671 return true;
1672 }
1673 return false;
1674 }
1675
1676
1677 template<typename T, Index Log2Dim>
1678 inline void
setOffsetOn(Index offset,const ValueType & val)1679 PointIndexLeafNode<T, Log2Dim>::setOffsetOn(Index offset, const ValueType& val)
1680 {
1681 this->buffer().setValue(offset, val);
1682 this->setValueMaskOn(offset);
1683 }
1684
1685
1686 template<typename T, Index Log2Dim>
1687 inline void
setOffsetOnly(Index offset,const ValueType & val)1688 PointIndexLeafNode<T, Log2Dim>::setOffsetOnly(Index offset, const ValueType& val)
1689 {
1690 this->buffer().setValue(offset, val);
1691 }
1692
1693
1694 template<typename T, Index Log2Dim>
1695 inline bool
isEmpty(const CoordBBox & bbox)1696 PointIndexLeafNode<T, Log2Dim>::isEmpty(const CoordBBox& bbox) const
1697 {
1698 Index xPos, pos, zStride = Index(bbox.max()[2] - bbox.min()[2]);
1699 Coord ijk;
1700
1701 for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
1702 xPos = (ijk[0] & (DIM - 1u)) << (2 * LOG2DIM);
1703
1704 for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
1705 pos = xPos + ((ijk[1] & (DIM - 1u)) << LOG2DIM);
1706 pos += (bbox.min()[2] & (DIM - 1u));
1707
1708 if (this->buffer()[pos+zStride] > (pos == 0 ? T(0) : this->buffer()[pos - 1])) {
1709 return false;
1710 }
1711 }
1712 }
1713
1714 return true;
1715 }
1716
1717
1718 template<typename T, Index Log2Dim>
1719 inline void
readBuffers(std::istream & is,bool fromHalf)1720 PointIndexLeafNode<T, Log2Dim>::readBuffers(std::istream& is, bool fromHalf)
1721 {
1722 BaseLeaf::readBuffers(is, fromHalf);
1723
1724 Index64 numIndices = Index64(0);
1725 is.read(reinterpret_cast<char*>(&numIndices), sizeof(Index64));
1726
1727 mIndices.resize(size_t(numIndices));
1728 is.read(reinterpret_cast<char*>(mIndices.data()), numIndices * sizeof(T));
1729 }
1730
1731
1732 template<typename T, Index Log2Dim>
1733 inline void
readBuffers(std::istream & is,const CoordBBox & bbox,bool fromHalf)1734 PointIndexLeafNode<T, Log2Dim>::readBuffers(std::istream& is, const CoordBBox& bbox, bool fromHalf)
1735 {
1736 // Read and clip voxel values.
1737 BaseLeaf::readBuffers(is, bbox, fromHalf);
1738
1739 Index64 numIndices = Index64(0);
1740 is.read(reinterpret_cast<char*>(&numIndices), sizeof(Index64));
1741
1742 const Index64 numBytes = numIndices * sizeof(T);
1743
1744 if (bbox.hasOverlap(this->getNodeBoundingBox())) {
1745 mIndices.resize(size_t(numIndices));
1746 is.read(reinterpret_cast<char*>(mIndices.data()), numBytes);
1747
1748 /// @todo If any voxels were deactivated as a result of clipping in the call to
1749 /// BaseLeaf::readBuffers(), the point index list will need to be regenerated.
1750 } else {
1751 // Read and discard voxel values.
1752 std::unique_ptr<char[]> buf{new char[numBytes]};
1753 is.read(buf.get(), numBytes);
1754 }
1755
1756 // Reserved for future use
1757 Index64 auxDataBytes = Index64(0);
1758 is.read(reinterpret_cast<char*>(&auxDataBytes), sizeof(Index64));
1759 if (auxDataBytes > 0) {
1760 // For now, read and discard any auxiliary data.
1761 std::unique_ptr<char[]> auxData{new char[auxDataBytes]};
1762 is.read(auxData.get(), auxDataBytes);
1763 }
1764 }
1765
1766
1767 template<typename T, Index Log2Dim>
1768 inline void
writeBuffers(std::ostream & os,bool toHalf)1769 PointIndexLeafNode<T, Log2Dim>::writeBuffers(std::ostream& os, bool toHalf) const
1770 {
1771 BaseLeaf::writeBuffers(os, toHalf);
1772
1773 Index64 numIndices = Index64(mIndices.size());
1774 os.write(reinterpret_cast<const char*>(&numIndices), sizeof(Index64));
1775 os.write(reinterpret_cast<const char*>(mIndices.data()), numIndices * sizeof(T));
1776
1777 // Reserved for future use
1778 const Index64 auxDataBytes = Index64(0);
1779 os.write(reinterpret_cast<const char*>(&auxDataBytes), sizeof(Index64));
1780 }
1781
1782
1783 template<typename T, Index Log2Dim>
1784 inline Index64
memUsage()1785 PointIndexLeafNode<T, Log2Dim>::memUsage() const
1786 {
1787 return BaseLeaf::memUsage() + Index64((sizeof(T)*mIndices.capacity()) + sizeof(mIndices));
1788 }
1789
1790 } // namespace tools
1791
1792
1793 ////////////////////////////////////////
1794
1795
1796 namespace tree {
1797
1798 /// Helper metafunction used to implement LeafNode::SameConfiguration
1799 /// (which, as an inner class, can't be independently specialized)
1800 template<Index Dim1, typename T2>
1801 struct SameLeafConfig<Dim1, openvdb::tools::PointIndexLeafNode<T2, Dim1> >
1802 {
1803 static const bool value = true;
1804 };
1805
1806 } // namespace tree
1807 } // namespace OPENVDB_VERSION_NAME
1808 } // namespace openvdb
1809
1810 #endif // OPENVDB_TOOLS_POINT_INDEX_GRID_HAS_BEEN_INCLUDED
1811