OpenVDB  11.0.0
PointPartitioner.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 /// @file PointPartitioner.h
5 ///
6 /// @brief Spatially partitions points using a parallel radix-based
7 /// sorting algorithm.
8 ///
9 /// @details Performs a stable deterministic sort; partitioning the same
10 /// point sequence will produce the same result each time.
11 /// @details The algorithm is unbounded meaning that points may be
12 /// distributed anywhere in index space.
13 /// @details The actual points are never stored in the tool, only
14 /// offsets into an external array.
15 ///
16 /// @author Mihai Alden
17 
18 #ifndef OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
19 #define OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
20 
21 #include <openvdb/Types.h>
22 #include <openvdb/math/Transform.h>
23 
24 #include <tbb/blocked_range.h>
25 #include <tbb/parallel_for.h>
26 #include <tbb/task_arena.h>
27 
28 #include <algorithm>
29 #include <cmath> // for std::isfinite()
30 #include <deque>
31 #include <map>
32 #include <set>
33 #include <utility> // std::pair
34 #include <vector>
35 
36 
37 namespace openvdb {
39 namespace OPENVDB_VERSION_NAME {
40 namespace tools {
41 
42 
43 ////////////////////////////////////////
44 
45 
46 /// @brief Partitions points into @c BucketLog2Dim aligned buckets
47 /// using a parallel radix-based sorting algorithm.
48 ///
49 /// @interface PointArray
50 /// Expected interface for the PointArray container:
51 /// @code
52 /// template<typename VectorType>
53 /// struct PointArray
54 /// {
55 /// // The type used to represent world-space point positions
56 /// using PosType = VectorType;
57 ///
58 /// // Return the number of points in the array
59 /// size_t size() const;
60 ///
61 /// // Return the world-space position of the nth point in the array.
62 /// void getPos(size_t n, PosType& xyz) const;
63 /// };
64 /// @endcode
65 ///
66 /// @details Performs a stable deterministic sort; partitioning the same
67 /// point sequence will produce the same result each time.
68 /// @details The algorithm is unbounded meaning that points may be
69 /// distributed anywhere in index space.
70 /// @details The actual points are never stored in the tool, only
71 /// offsets into an external array.
72 /// @details @c BucketLog2Dim defines the bucket coordinate dimensions,
73 /// i.e. BucketLog2Dim = 3 corresponds to a bucket that spans
74 /// a (2^3)^3 = 8^3 voxel region.
75 template<typename PointIndexType = uint32_t, Index BucketLog2Dim = 3>
77 {
78 public:
79  enum { LOG2DIM = BucketLog2Dim };
80 
83 
84  using IndexType = PointIndexType;
85 
86  static constexpr Index bits = 1 + (3 * BucketLog2Dim);
87  // signed, so if bits is exactly 16, int32 is required
88  using VoxelOffsetType = typename std::conditional<(bits < 16),
89  int16_t, typename std::conditional<(bits < 32), int32_t, int64_t>::type>::type;
90 
91  using VoxelOffsetArray = std::unique_ptr<VoxelOffsetType[]>;
92 
93  class IndexIterator;
94 
95  //////////
96 
98 
99  /// @brief Partitions point indices into @c BucketLog2Dim aligned buckets.
100  ///
101  /// @param points list of world space points.
102  /// @param xform world to index space transform.
103  /// @param voxelOrder sort point indices by local voxel offsets.
104  /// @param recordVoxelOffsets construct local voxel offsets
105  /// @param cellCenteredTransform toggle the cell-centered interpretation that imagines world
106  /// space as divided into discrete cells (e.g., cubes) centered
107  /// on the image of the index-space lattice points.
108  template<typename PointArray>
109  void construct(const PointArray& points, const math::Transform& xform,
110  bool voxelOrder = false, bool recordVoxelOffsets = false,
111  bool cellCenteredTransform = true);
112 
113 
114  /// @brief Partitions point indices into @c BucketLog2Dim aligned buckets.
115  ///
116  /// @param points list of world space points.
117  /// @param xform world to index space transform.
118  /// @param voxelOrder sort point indices by local voxel offsets.
119  /// @param recordVoxelOffsets construct local voxel offsets
120  /// @param cellCenteredTransform toggle the cell-centered interpretation that imagines world
121  /// space as divided into discrete cells (e.g., cubes) centered
122  /// on the image of the index-space lattice points.
123  template<typename PointArray>
124  static Ptr create(const PointArray& points, const math::Transform& xform,
125  bool voxelOrder = false, bool recordVoxelOffsets = false,
126  bool cellCenteredTransform = true);
127 
128 
129  /// @brief Returns the number of buckets.
130  size_t size() const { return mPageCount; }
131 
132  /// @brief true if the container size is 0, false otherwise.
133  bool empty() const { return mPageCount == 0; }
134 
135  /// @brief Removes all data and frees up memory.
136  void clear();
137 
138  /// @brief Exchanges the content of the container by another.
139  void swap(PointPartitioner&);
140 
141  /// @brief Returns the point indices for bucket @a n
142  IndexIterator indices(size_t n) const;
143 
144  /// @brief Returns the coordinate-aligned bounding box for bucket @a n
145  CoordBBox getBBox(size_t n) const {
146  return CoordBBox::createCube(mPageCoordinates[n], (1u << BucketLog2Dim));
147  }
148 
149  /// @brief Returns the origin coordinate for bucket @a n
150  const Coord& origin(size_t n) const { return mPageCoordinates[n]; }
151 
152  /// @brief Returns a list of @c LeafNode voxel offsets for the points.
153  /// @note The list is optionally constructed.
154  const VoxelOffsetArray& voxelOffsets() const { return mVoxelOffsets; }
155 
156  /// @brief Returns @c true if this point partitioning was constructed
157  /// using a cell-centered transform.
158  /// @note Cell-centered interpretation is the default behavior.
159  bool usingCellCenteredTransform() const { return mUsingCellCenteredTransform; }
160 
161 private:
162  // Disallow copying
164  PointPartitioner& operator=(const PointPartitioner&);
165 
166  std::unique_ptr<IndexType[]> mPointIndices;
167  VoxelOffsetArray mVoxelOffsets;
168 
169  std::unique_ptr<IndexType[]> mPageOffsets;
170  std::unique_ptr<Coord[]> mPageCoordinates;
171  IndexType mPageCount;
172  bool mUsingCellCenteredTransform;
173 }; // class PointPartitioner
174 
175 
177 
178 
179 template<typename PointIndexType, Index BucketLog2Dim>
180 class PointPartitioner<PointIndexType, BucketLog2Dim>::IndexIterator
181 {
182 public:
183  using IndexType = PointIndexType;
184 
185  IndexIterator(IndexType* begin = nullptr, IndexType* end = nullptr)
186  : mBegin(begin), mEnd(end), mItem(begin) {}
187 
188  /// @brief Rewind to first item.
189  void reset() { mItem = mBegin; }
190 
191  /// @brief Number of point indices in the iterator range.
192  size_t size() const { return mEnd - mBegin; }
193 
194  /// @brief Returns the item to which this iterator is currently pointing.
195  IndexType& operator*() { assert(mItem != nullptr); return *mItem; }
196  const IndexType& operator*() const { assert(mItem != nullptr); return *mItem; }
197 
198  /// @brief Return @c true if this iterator is not yet exhausted.
199  operator bool() const { return mItem < mEnd; }
200  bool test() const { return mItem < mEnd; }
201 
202  /// @brief Advance to the next item.
203  IndexIterator& operator++() { assert(this->test()); ++mItem; return *this; }
204 
205  /// @brief Advance to the next item.
206  bool next() { this->operator++(); return this->test(); }
207  bool increment() { this->next(); return this->test(); }
208 
209  /// @brief Equality operators
210  bool operator==(const IndexIterator& other) const { return mItem == other.mItem; }
211  bool operator!=(const IndexIterator& other) const { return !this->operator==(other); }
212 
213 private:
214  IndexType * const mBegin, * const mEnd;
215  IndexType * mItem;
216 }; // class PointPartitioner::IndexIterator
217 
218 
219 ////////////////////////////////////////
220 ////////////////////////////////////////
221 
222 // Implementation details
223 
224 /// @cond OPENVDB_DOCS_INTERNAL
225 
226 namespace point_partitioner_internal {
227 
228 
229 template<typename PointIndexType>
230 struct ComputePointOrderOp
231 {
232  ComputePointOrderOp(PointIndexType* pointOrder,
233  const PointIndexType* bucketCounters, const PointIndexType* bucketOffsets)
234  : mPointOrder(pointOrder)
235  , mBucketCounters(bucketCounters)
236  , mBucketOffsets(bucketOffsets)
237  {
238  }
239 
240  void operator()(const tbb::blocked_range<size_t>& range) const {
241  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
242  mPointOrder[n] += mBucketCounters[mBucketOffsets[n]];
243  }
244  }
245 
246  PointIndexType * const mPointOrder;
247  PointIndexType const * const mBucketCounters;
248  PointIndexType const * const mBucketOffsets;
249 }; // struct ComputePointOrderOp
250 
251 
252 template<typename PointIndexType>
253 struct CreateOrderedPointIndexArrayOp
254 {
255  CreateOrderedPointIndexArrayOp(PointIndexType* orderedIndexArray,
256  const PointIndexType* pointOrder, const PointIndexType* indices)
257  : mOrderedIndexArray(orderedIndexArray)
258  , mPointOrder(pointOrder)
259  , mIndices(indices)
260  {
261  }
262 
263  void operator()(const tbb::blocked_range<size_t>& range) const {
264  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
265  mOrderedIndexArray[mPointOrder[n]] = mIndices[n];
266  }
267  }
268 
269  PointIndexType * const mOrderedIndexArray;
270  PointIndexType const * const mPointOrder;
271  PointIndexType const * const mIndices;
272 }; // struct CreateOrderedPointIndexArrayOp
273 
274 
275 template<typename PointIndexType, Index BucketLog2Dim>
276 struct VoxelOrderOp
277 {
278  static constexpr Index bits = 1 + (3 * BucketLog2Dim);
279  // signed, so if bits is exactly 16, int32 is required
280  using VoxelOffsetType = typename std::conditional<(bits < 16),
281  int16_t, typename std::conditional<(bits < 32), int32_t, int64_t>::type>::type;
282 
283  using VoxelOffsetArray = std::unique_ptr<VoxelOffsetType[]>;
284  using IndexArray = std::unique_ptr<PointIndexType[]>;
285 
286  VoxelOrderOp(IndexArray& indices, const IndexArray& pages,const VoxelOffsetArray& offsets)
287  : mIndices(indices.get())
288  , mPages(pages.get())
289  , mVoxelOffsets(offsets.get())
290  {
291  }
292 
293  void operator()(const tbb::blocked_range<size_t>& range) const {
294 
295  PointIndexType pointCount = 0;
296  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
297  pointCount = std::max(pointCount, (mPages[n + 1] - mPages[n]));
298  }
299 
300  const PointIndexType voxelCount = 1 << (3 * BucketLog2Dim);
301 
302  // allocate histogram buffers
303  std::unique_ptr<VoxelOffsetType[]> offsets(new VoxelOffsetType[pointCount]);
304  std::unique_ptr<PointIndexType[]> sortedIndices(new PointIndexType[pointCount]);
305  std::unique_ptr<PointIndexType[]> histogram(new PointIndexType[voxelCount]);
306 
307  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
308 
309  PointIndexType * const indices = mIndices + mPages[n];
310  pointCount = mPages[n + 1] - mPages[n];
311 
312  // local copy of voxel offsets.
313  for (PointIndexType i = 0; i < pointCount; ++i) {
314  offsets[i] = mVoxelOffsets[ indices[i] ];
315  }
316 
317  // reset histogram
318  memset(&histogram[0], 0, voxelCount * sizeof(PointIndexType));
319 
320  // compute histogram
321  for (PointIndexType i = 0; i < pointCount; ++i) {
322  ++histogram[ offsets[i] ];
323  }
324 
325  PointIndexType count = 0, startOffset;
326  for (int i = 0; i < int(voxelCount); ++i) {
327  if (histogram[i] > 0) {
328  startOffset = count;
329  count += histogram[i];
330  histogram[i] = startOffset;
331  }
332  }
333 
334  // sort indices based on voxel offset
335  for (PointIndexType i = 0; i < pointCount; ++i) {
336  sortedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
337  }
338 
339  memcpy(&indices[0], &sortedIndices[0], sizeof(PointIndexType) * pointCount);
340  }
341  }
342 
343  PointIndexType * const mIndices;
344  PointIndexType const * const mPages;
345  VoxelOffsetType const * const mVoxelOffsets;
346 }; // struct VoxelOrderOp
347 
348 
349 ////////////////////////////////////////
350 
351 
352 template<typename T>
353 struct Array
354 {
355  using Ptr = std::unique_ptr<Array>;
356 
357  Array(size_t size) : mSize(size), mData(new T[size]) { }
358 
359  size_t size() const { return mSize; }
360 
361  T* data() { return mData.get(); }
362  const T* data() const { return mData.get(); }
363 
364  void clear() { mSize = 0; mData.reset(); }
365 
366 private:
367  size_t mSize;
368  std::unique_ptr<T[]> mData;
369 }; // struct Array
370 
371 
372 template<typename PointIndexType>
373 struct MoveSegmentDataOp
374 {
375  using SegmentPtr = typename Array<PointIndexType>::Ptr;
376 
377  MoveSegmentDataOp(std::vector<PointIndexType*>& indexLists, SegmentPtr* segments)
378  : mIndexLists(&indexLists[0]), mSegments(segments)
379  {
380  }
381 
382  void operator()(const tbb::blocked_range<size_t>& range) const {
383  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
384  PointIndexType* indices = mIndexLists[n];
385  SegmentPtr& segment = mSegments[n];
386 
387  tbb::parallel_for(tbb::blocked_range<size_t>(0, segment->size()),
388  CopyData(indices, segment->data()));
389 
390  segment.reset(); // clear data
391  }
392  }
393 
394 private:
395 
396  struct CopyData
397  {
398  CopyData(PointIndexType* lhs, const PointIndexType* rhs) : mLhs(lhs), mRhs(rhs) { }
399 
400  void operator()(const tbb::blocked_range<size_t>& range) const {
401  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
402  mLhs[n] = mRhs[n];
403  }
404  }
405 
406  PointIndexType * const mLhs;
407  PointIndexType const * const mRhs;
408  };
409 
410  PointIndexType * const * const mIndexLists;
411  SegmentPtr * const mSegments;
412 }; // struct MoveSegmentDataOp
413 
414 
415 template<typename PointIndexType>
416 struct MergeBinsOp
417 {
418  using Segment = Array<PointIndexType>;
419  using SegmentPtr = typename Segment::Ptr;
420 
421  using IndexPair = std::pair<PointIndexType, PointIndexType>;
422  using IndexPairList = std::deque<IndexPair>;
423  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
424  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
425  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
426 
427  MergeBinsOp(IndexPairListMapPtr* bins,
428  SegmentPtr* indexSegments,
429  SegmentPtr* offsetSegments,
430  Coord* coords,
431  size_t numSegments)
432  : mBins(bins)
433  , mIndexSegments(indexSegments)
434  , mOffsetSegments(offsetSegments)
435  , mCoords(coords)
436  , mNumSegments(numSegments)
437  {
438  }
439 
440  void operator()(const tbb::blocked_range<size_t>& range) const {
441 
442  std::vector<IndexPairListPtr*> data;
443  std::vector<PointIndexType> arrayOffsets;
444 
445  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
446 
447  const Coord& ijk = mCoords[n];
448  size_t numIndices = 0;
449 
450  data.clear();
451 
452  for (size_t i = 0, I = mNumSegments; i < I; ++i) {
453 
454  IndexPairListMap& idxMap = *mBins[i];
455  typename IndexPairListMap::iterator iter = idxMap.find(ijk);
456 
457  if (iter != idxMap.end() && iter->second) {
458  IndexPairListPtr& idxListPtr = iter->second;
459 
460  data.push_back(&idxListPtr);
461  numIndices += idxListPtr->size();
462  }
463  }
464 
465  if (data.empty() || numIndices == 0) continue;
466 
467  SegmentPtr& indexSegment = mIndexSegments[n];
468  SegmentPtr& offsetSegment = mOffsetSegments[n];
469 
470  indexSegment.reset(new Segment(numIndices));
471  offsetSegment.reset(new Segment(numIndices));
472 
473  arrayOffsets.clear();
474  arrayOffsets.reserve(data.size());
475 
476  for (size_t i = 0, count = 0, I = data.size(); i < I; ++i) {
477  arrayOffsets.push_back(PointIndexType(count));
478  count += (*data[i])->size();
479  }
480 
481  tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size()),
482  CopyData(&data[0], &arrayOffsets[0], indexSegment->data(), offsetSegment->data()));
483  }
484  }
485 
486 private:
487 
488  struct CopyData
489  {
490  CopyData(IndexPairListPtr** indexLists,
491  const PointIndexType* arrayOffsets,
492  PointIndexType* indices,
493  PointIndexType* offsets)
494  : mIndexLists(indexLists)
495  , mArrayOffsets(arrayOffsets)
496  , mIndices(indices)
497  , mOffsets(offsets)
498  {
499  }
500 
501  void operator()(const tbb::blocked_range<size_t>& range) const {
502 
503  using CIter = typename IndexPairList::const_iterator;
504 
505  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
506 
507  const PointIndexType arrayOffset = mArrayOffsets[n];
508  PointIndexType* indexPtr = &mIndices[arrayOffset];
509  PointIndexType* offsetPtr = &mOffsets[arrayOffset];
510 
511  IndexPairListPtr& list = *mIndexLists[n];
512 
513  for (CIter it = list->begin(), end = list->end(); it != end; ++it) {
514  const IndexPair& data = *it;
515  *indexPtr++ = data.first;
516  *offsetPtr++ = data.second;
517  }
518 
519  list.reset(); // clear data
520  }
521  }
522 
523  IndexPairListPtr * const * const mIndexLists;
524  PointIndexType const * const mArrayOffsets;
525  PointIndexType * const mIndices;
526  PointIndexType * const mOffsets;
527  }; // struct CopyData
528 
529  IndexPairListMapPtr * const mBins;
530  SegmentPtr * const mIndexSegments;
531  SegmentPtr * const mOffsetSegments;
532  Coord const * const mCoords;
533  size_t const mNumSegments;
534 }; // struct MergeBinsOp
535 
536 
537 template<typename PointArray, typename PointIndexType, typename VoxelOffsetType>
538 struct BinPointIndicesOp
539 {
540  using PosType = typename PointArray::PosType;
541  using IndexPair = std::pair<PointIndexType, PointIndexType>;
542  using IndexPairList = std::deque<IndexPair>;
543  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
544  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
545  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
546 
547  BinPointIndicesOp(IndexPairListMapPtr* data,
548  const PointArray& points,
549  VoxelOffsetType* voxelOffsets,
550  const math::Transform& m,
551  Index binLog2Dim,
552  Index bucketLog2Dim,
553  size_t numSegments,
554  bool cellCenteredTransform)
555  : mData(data)
556  , mPoints(&points)
557  , mVoxelOffsets(voxelOffsets)
558  , mXForm(m)
559  , mBinLog2Dim(binLog2Dim)
560  , mBucketLog2Dim(bucketLog2Dim)
561  , mNumSegments(numSegments)
562  , mCellCenteredTransform(cellCenteredTransform)
563  {
564  }
565 
566  void operator()(const tbb::blocked_range<size_t>& range) const {
567 
568  const Index log2dim = mBucketLog2Dim;
569  const Index log2dim2 = 2 * log2dim;
570  const Index bucketMask = (1u << log2dim) - 1u;
571 
572  const Index binLog2dim = mBinLog2Dim;
573  const Index binLog2dim2 = 2 * binLog2dim;
574 
575  const Index binMask = (1u << (log2dim + binLog2dim)) - 1u;
576  const Index invBinMask = ~binMask;
577 
578  IndexPairList * idxList = nullptr;
579  Coord ijk(0, 0, 0), loc(0, 0, 0), binCoord(0, 0, 0), lastBinCoord(1, 2, 3);
580  PosType pos;
581 
582  PointIndexType bucketOffset = 0;
583  VoxelOffsetType voxelOffset = 0;
584 
585  const bool cellCentered = mCellCenteredTransform;
586 
587  const size_t numPoints = mPoints->size();
588  const size_t segmentSize = numPoints / mNumSegments;
589 
590  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
591 
592  IndexPairListMapPtr& dataPtr = mData[n];
593  if (!dataPtr) dataPtr.reset(new IndexPairListMap());
594  IndexPairListMap& idxMap = *dataPtr;
595 
596  const bool isLastSegment = (n + 1) >= mNumSegments;
597 
598  const size_t start = n * segmentSize;
599  const size_t end = isLastSegment ? numPoints : (start + segmentSize);
600 
601  for (size_t i = start; i != end; ++i) {
602 
603  mPoints->getPos(i, pos);
604 
605  if (std::isfinite(pos[0]) && std::isfinite(pos[1]) && std::isfinite(pos[2])) {
606  ijk = cellCentered ? mXForm.worldToIndexCellCentered(pos) :
607  mXForm.worldToIndexNodeCentered(pos);
608 
609  if (mVoxelOffsets) {
610  loc[0] = ijk[0] & bucketMask;
611  loc[1] = ijk[1] & bucketMask;
612  loc[2] = ijk[2] & bucketMask;
613  voxelOffset = VoxelOffsetType(
614  (loc[0] << log2dim2) + (loc[1] << log2dim) + loc[2]);
615  }
616 
617  binCoord[0] = ijk[0] & invBinMask;
618  binCoord[1] = ijk[1] & invBinMask;
619  binCoord[2] = ijk[2] & invBinMask;
620 
621  ijk[0] &= binMask;
622  ijk[1] &= binMask;
623  ijk[2] &= binMask;
624 
625  ijk[0] >>= log2dim;
626  ijk[1] >>= log2dim;
627  ijk[2] >>= log2dim;
628 
629  bucketOffset = PointIndexType(
630  (ijk[0] << binLog2dim2) + (ijk[1] << binLog2dim) + ijk[2]);
631 
632  if (lastBinCoord != binCoord) {
633  lastBinCoord = binCoord;
634  IndexPairListPtr& idxListPtr = idxMap[lastBinCoord];
635  if (!idxListPtr) idxListPtr.reset(new IndexPairList());
636  idxList = idxListPtr.get();
637  }
638 
639  idxList->push_back(IndexPair(PointIndexType(i), bucketOffset));
640  if (mVoxelOffsets) mVoxelOffsets[i] = voxelOffset;
641  }
642  }
643  }
644  }
645 
646  IndexPairListMapPtr * const mData;
647  PointArray const * const mPoints;
648  VoxelOffsetType * const mVoxelOffsets;
649  math::Transform const mXForm;
650  Index const mBinLog2Dim;
651  Index const mBucketLog2Dim;
652  size_t const mNumSegments;
653  bool const mCellCenteredTransform;
654 }; // struct BinPointIndicesOp
655 
656 
657 template<typename PointIndexType>
658 struct OrderSegmentsOp
659 {
660  using IndexArray = std::unique_ptr<PointIndexType[]>;
661  using SegmentPtr = typename Array<PointIndexType>::Ptr;
662 
663  OrderSegmentsOp(SegmentPtr* indexSegments, SegmentPtr* offsetSegments,
664  IndexArray* pageOffsetArrays, IndexArray* pageIndexArrays, Index binVolume)
665  : mIndexSegments(indexSegments)
666  , mOffsetSegments(offsetSegments)
667  , mPageOffsetArrays(pageOffsetArrays)
668  , mPageIndexArrays(pageIndexArrays)
669  , mBinVolume(binVolume)
670  {
671  }
672 
673  void operator()(const tbb::blocked_range<size_t>& range) const {
674 
675  const size_t bucketCountersSize = size_t(mBinVolume);
676  IndexArray bucketCounters(new PointIndexType[bucketCountersSize]);
677 
678  size_t maxSegmentSize = 0;
679  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
680  maxSegmentSize = std::max(maxSegmentSize, mIndexSegments[n]->size());
681  }
682 
683  IndexArray bucketIndices(new PointIndexType[maxSegmentSize]);
684 
685  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
686 
687  memset(bucketCounters.get(), 0, sizeof(PointIndexType) * bucketCountersSize);
688 
689  const size_t segmentSize = mOffsetSegments[n]->size();
690  PointIndexType* offsets = mOffsetSegments[n]->data();
691 
692  // Count the number of points per bucket and assign a local bucket index
693  // to each point.
694  for (size_t i = 0; i < segmentSize; ++i) {
695  bucketIndices[i] = bucketCounters[offsets[i]]++;
696  }
697 
698  PointIndexType nonemptyBucketCount = 0;
699  for (size_t i = 0; i < bucketCountersSize; ++i) {
700  nonemptyBucketCount += static_cast<PointIndexType>(bucketCounters[i] != 0);
701  }
702 
703 
704  IndexArray& pageOffsets = mPageOffsetArrays[n];
705  pageOffsets.reset(new PointIndexType[nonemptyBucketCount + 1]);
706  pageOffsets[0] = nonemptyBucketCount + 1; // stores array size in first element
707 
708  IndexArray& pageIndices = mPageIndexArrays[n];
709  pageIndices.reset(new PointIndexType[nonemptyBucketCount]);
710 
711  // Compute bucket counter prefix sum
712  PointIndexType count = 0, idx = 0;
713  for (size_t i = 0; i < bucketCountersSize; ++i) {
714  if (bucketCounters[i] != 0) {
715  pageIndices[idx] = static_cast<PointIndexType>(i);
716  pageOffsets[idx+1] = bucketCounters[i];
717  bucketCounters[i] = count;
718  count += pageOffsets[idx+1];
719  ++idx;
720  }
721  }
722 
723  PointIndexType* indices = mIndexSegments[n]->data();
724  const tbb::blocked_range<size_t> segmentRange(0, segmentSize);
725 
726  // Compute final point order by incrementing the local bucket point index
727  // with the prefix sum offset.
728  tbb::parallel_for(segmentRange, ComputePointOrderOp<PointIndexType>(
729  bucketIndices.get(), bucketCounters.get(), offsets));
730 
731  tbb::parallel_for(segmentRange, CreateOrderedPointIndexArrayOp<PointIndexType>(
732  offsets, bucketIndices.get(), indices));
733 
734  mIndexSegments[n]->clear(); // clear data
735  }
736  }
737 
738  SegmentPtr * const mIndexSegments;
739  SegmentPtr * const mOffsetSegments;
740  IndexArray * const mPageOffsetArrays;
741  IndexArray * const mPageIndexArrays;
742  Index const mBinVolume;
743 }; // struct OrderSegmentsOp
744 
745 
746 ////////////////////////////////////////
747 
748 
749 /// @brief Segment points using one level of least significant digit radix bins.
750 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
751 inline void binAndSegment(
752  const PointArray& points,
753  const math::Transform& xform,
754  std::unique_ptr<typename Array<PointIndexType>::Ptr[]>& indexSegments,
755  std::unique_ptr<typename Array<PointIndexType>::Ptr[]>& offsetSegments,
756  std::vector<Coord>& coords,
757  const Index binLog2Dim,
758  const Index bucketLog2Dim,
759  VoxelOffsetType* voxelOffsets = nullptr,
760  bool cellCenteredTransform = true)
761 {
762  using IndexPair = std::pair<PointIndexType, PointIndexType>;
763  using IndexPairList = std::deque<IndexPair>;
764  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
765  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
766  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
767 
768  size_t numTasks = 1, numThreads = size_t(tbb::this_task_arena::max_concurrency());
769  if (points.size() > (numThreads * 2)) numTasks = numThreads * 2;
770  else if (points.size() > numThreads) numTasks = numThreads;
771 
772  std::unique_ptr<IndexPairListMapPtr[]> bins(new IndexPairListMapPtr[numTasks]);
773 
774  using BinOp = BinPointIndicesOp<PointArray, PointIndexType, VoxelOffsetType>;
775 
776  tbb::parallel_for(tbb::blocked_range<size_t>(0, numTasks),
777  BinOp(bins.get(), points, voxelOffsets, xform, binLog2Dim, bucketLog2Dim,
778  numTasks, cellCenteredTransform));
779 
780  std::set<Coord> uniqueCoords;
781 
782  for (size_t i = 0; i < numTasks; ++i) {
783  IndexPairListMap& idxMap = *bins[i];
784  for (typename IndexPairListMap::iterator it = idxMap.begin(); it != idxMap.end(); ++it) {
785  uniqueCoords.insert(it->first);
786  }
787  }
788 
789  coords.assign(uniqueCoords.begin(), uniqueCoords.end());
790  uniqueCoords.clear();
791 
792  size_t segmentCount = coords.size();
793 
794  using SegmentPtr = typename Array<PointIndexType>::Ptr;
795 
796  indexSegments.reset(new SegmentPtr[segmentCount]);
797  offsetSegments.reset(new SegmentPtr[segmentCount]);
798 
799  using MergeOp = MergeBinsOp<PointIndexType>;
800 
801  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
802  MergeOp(bins.get(), indexSegments.get(), offsetSegments.get(), &coords[0], numTasks));
803 }
804 
805 
806 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
807 inline void partition(
808  const PointArray& points,
809  const math::Transform& xform,
810  const Index bucketLog2Dim,
811  std::unique_ptr<PointIndexType[]>& pointIndices,
812  std::unique_ptr<PointIndexType[]>& pageOffsets,
813  std::unique_ptr<Coord[]>& pageCoordinates,
814  PointIndexType& pageCount,
815  std::unique_ptr<VoxelOffsetType[]>& voxelOffsets,
816  bool recordVoxelOffsets,
817  bool cellCenteredTransform)
818 {
819  using SegmentPtr = typename Array<PointIndexType>::Ptr;
820 
821  if (recordVoxelOffsets) voxelOffsets.reset(new VoxelOffsetType[points.size()]);
822  else voxelOffsets.reset();
823 
824  const Index binLog2Dim = 5u;
825  // note: Bins span a (2^(binLog2Dim + bucketLog2Dim))^3 voxel region,
826  // i.e. bucketLog2Dim = 3 and binLog2Dim = 5 corresponds to a
827  // (2^8)^3 = 256^3 voxel region.
828 
829 
830  std::vector<Coord> segmentCoords;
831 
832  std::unique_ptr<SegmentPtr[]> indexSegments;
833  std::unique_ptr<SegmentPtr[]> offsetSegments;
834 
835  binAndSegment<PointIndexType, VoxelOffsetType, PointArray>(points, xform,
836  indexSegments, offsetSegments, segmentCoords, binLog2Dim, bucketLog2Dim,
837  voxelOffsets.get(), cellCenteredTransform);
838 
839  size_t numSegments = segmentCoords.size();
840 
841  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
842 
843  using IndexArray = std::unique_ptr<PointIndexType[]>;
844  std::unique_ptr<IndexArray[]> pageOffsetArrays(new IndexArray[numSegments]);
845  std::unique_ptr<IndexArray[]> pageIndexArrays(new IndexArray[numSegments]);
846 
847  const Index binVolume = 1u << (3u * binLog2Dim);
848 
849  tbb::parallel_for(segmentRange, OrderSegmentsOp<PointIndexType>
850  (indexSegments.get(), offsetSegments.get(),
851  pageOffsetArrays.get(), pageIndexArrays.get(), binVolume));
852 
853  indexSegments.reset();
854 
855  std::vector<Index> segmentOffsets;
856  segmentOffsets.reserve(numSegments);
857 
858  pageCount = 0;
859  for (size_t n = 0; n < numSegments; ++n) {
860  segmentOffsets.push_back(pageCount);
861  pageCount += pageOffsetArrays[n][0] - 1;
862  }
863 
864  pageOffsets.reset(new PointIndexType[pageCount + 1]);
865 
866  PointIndexType count = 0;
867  for (size_t n = 0, idx = 0; n < numSegments; ++n) {
868 
869  PointIndexType* offsets = pageOffsetArrays[n].get();
870  size_t size = size_t(offsets[0]);
871 
872  for (size_t i = 1; i < size; ++i) {
873  pageOffsets[idx++] = count;
874  count += offsets[i];
875  }
876  }
877 
878  pageOffsets[pageCount] = count;
879 
880  pointIndices.reset(new PointIndexType[points.size()]);
881 
882  std::vector<PointIndexType*> indexArray;
883  indexArray.reserve(numSegments);
884 
885  PointIndexType* index = pointIndices.get();
886  for (size_t n = 0; n < numSegments; ++n) {
887  indexArray.push_back(index);
888  index += offsetSegments[n]->size();
889  }
890 
891  // compute leaf node origin for each page
892 
893  pageCoordinates.reset(new Coord[pageCount]);
894 
895  tbb::parallel_for(segmentRange,
896  [&](tbb::blocked_range<size_t>& range)
897  {
898  for (size_t n = range.begin(); n < range.end(); n++)
899  {
900  Index segmentOffset = segmentOffsets[n];
901  PointIndexType* indices = pageIndexArrays[n].get();
902 
903  const Coord& segmentCoord = segmentCoords[n];
904 
905  // segment size stored in the first value of the offset array
906  const size_t segmentSize = pageOffsetArrays[n][0] - 1;
907  tbb::blocked_range<size_t> copyRange(0, segmentSize);
908  tbb::parallel_for(copyRange,
909  [&](tbb::blocked_range<size_t>& r)
910  {
911  for (size_t i = r.begin(); i < r.end(); i++)
912  {
913  Index pageIndex = indices[i];
914  Coord& ijk = pageCoordinates[segmentOffset+i];
915 
916  ijk[0] = pageIndex >> (2 * binLog2Dim);
917  Index pageIndexModulo = pageIndex - (ijk[0] << (2 * binLog2Dim));
918  ijk[1] = pageIndexModulo >> binLog2Dim;
919  ijk[2] = pageIndexModulo - (ijk[1] << binLog2Dim);
920 
921  ijk = (ijk << bucketLog2Dim) + segmentCoord;
922  }
923  }
924  );
925  }
926  }
927  );
928 
929  // move segment data
930 
931  tbb::parallel_for(segmentRange,
932  MoveSegmentDataOp<PointIndexType>(indexArray, offsetSegments.get()));
933 }
934 
935 
936 } // namespace point_partitioner_internal
937 
938 /// @endcond
939 
940 ////////////////////////////////////////
941 
942 
943 template<typename PointIndexType, Index BucketLog2Dim>
945  : mPointIndices(nullptr)
946  , mVoxelOffsets(nullptr)
947  , mPageOffsets(nullptr)
948  , mPageCoordinates(nullptr)
949  , mPageCount(0)
950  , mUsingCellCenteredTransform(true)
951 {
952 }
953 
954 
955 template<typename PointIndexType, Index BucketLog2Dim>
956 inline void
958 {
959  mPageCount = 0;
960  mUsingCellCenteredTransform = true;
961  mPointIndices.reset();
962  mVoxelOffsets.reset();
963  mPageOffsets.reset();
964  mPageCoordinates.reset();
965 }
966 
967 
968 template<typename PointIndexType, Index BucketLog2Dim>
969 inline void
971 {
972  const IndexType tmpLhsPageCount = mPageCount;
973  mPageCount = rhs.mPageCount;
974  rhs.mPageCount = tmpLhsPageCount;
975 
976  mPointIndices.swap(rhs.mPointIndices);
977  mVoxelOffsets.swap(rhs.mVoxelOffsets);
978  mPageOffsets.swap(rhs.mPageOffsets);
979  mPageCoordinates.swap(rhs.mPageCoordinates);
980 
981  bool lhsCellCenteredTransform = mUsingCellCenteredTransform;
982  mUsingCellCenteredTransform = rhs.mUsingCellCenteredTransform;
983  rhs.mUsingCellCenteredTransform = lhsCellCenteredTransform;
984 }
985 
986 
987 template<typename PointIndexType, Index BucketLog2Dim>
990 {
991  assert(bool(mPointIndices) && bool(mPageCount));
992  return IndexIterator(
993  mPointIndices.get() + mPageOffsets[n],
994  mPointIndices.get() + mPageOffsets[n + 1]);
995 }
996 
997 
998 template<typename PointIndexType, Index BucketLog2Dim>
999 template<typename PointArray>
1000 inline void
1002  const PointArray& points,
1003  const math::Transform& xform,
1004  bool voxelOrder,
1005  bool recordVoxelOffsets,
1006  bool cellCenteredTransform)
1007 {
1008  mUsingCellCenteredTransform = cellCenteredTransform;
1009 
1010  point_partitioner_internal::partition(points, xform, BucketLog2Dim,
1011  mPointIndices, mPageOffsets, mPageCoordinates, mPageCount, mVoxelOffsets,
1012  (voxelOrder || recordVoxelOffsets), cellCenteredTransform);
1013 
1014  const tbb::blocked_range<size_t> pageRange(0, mPageCount);
1015 
1016  if (mVoxelOffsets && voxelOrder) {
1017  tbb::parallel_for(pageRange, point_partitioner_internal::VoxelOrderOp<
1018  IndexType, BucketLog2Dim>(mPointIndices, mPageOffsets, mVoxelOffsets));
1019  }
1020 
1021  if (mVoxelOffsets && !recordVoxelOffsets) {
1022  mVoxelOffsets.reset();
1023  }
1024 }
1025 
1026 
1027 template<typename PointIndexType, Index BucketLog2Dim>
1028 template<typename PointArray>
1031  const PointArray& points,
1032  const math::Transform& xform,
1033  bool voxelOrder,
1034  bool recordVoxelOffsets,
1035  bool cellCenteredTransform)
1036 {
1037  Ptr ret(new PointPartitioner());
1038  ret->construct(points, xform, voxelOrder, recordVoxelOffsets, cellCenteredTransform);
1039  return ret;
1040 }
1041 
1042 
1043 ////////////////////////////////////////
1044 
1045 
1046 } // namespace tools
1047 } // namespace OPENVDB_VERSION_NAME
1048 } // namespace openvdb
1049 
1050 
1051 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
size_t size() const
Returns the number of buckets.
Definition: PointPartitioner.h:130
const VoxelOffsetArray & voxelOffsets() const
Returns a list of LeafNode voxel offsets for the points.
Definition: PointPartitioner.h:154
math::Histogram histogram(const IterT &iter, double minVal, double maxVal, size_t numBins=10, bool threaded=true)
Iterate over a scalar grid and compute a histogram of the values of the voxels that are visited...
Definition: Statistics.h:343
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Definition: PointCountImpl.h:18
std::pair< Index, Index > IndexPair
Definition: PointMoveImpl.h:95
std::shared_ptr< T > SharedPtr
Definition: Types.h:114
IndexIterator(IndexType *begin=nullptr, IndexType *end=nullptr)
Definition: PointPartitioner.h:185
std::unique_ptr< VoxelOffsetType[]> VoxelOffsetArray
Definition: PointPartitioner.h:91
bool next()
Advance to the next item.
Definition: PointPartitioner.h:206
bool increment()
Definition: PointPartitioner.h:207
BBox< Coord > CoordBBox
Definition: NanoVDB.h:2535
SharedPtr< PointPartitioner > Ptr
Definition: PointPartitioner.h:81
PointIndexType IndexType
Definition: PointPartitioner.h:84
IndexIterator & operator++()
Advance to the next item.
Definition: PointPartitioner.h:203
bool empty() const
true if the container size is 0, false otherwise.
Definition: PointPartitioner.h:133
bool operator==(const IndexIterator &other) const
Equality operators.
Definition: PointPartitioner.h:210
const Coord & origin(size_t n) const
Returns the origin coordinate for bucket n.
Definition: PointPartitioner.h:150
void construct(const PointArray &points, const math::Transform &xform, bool voxelOrder=false, bool recordVoxelOffsets=false, bool cellCenteredTransform=true)
Partitions point indices into BucketLog2Dim aligned buckets.
Definition: PointPartitioner.h:1001
Definition: PointPartitioner.h:76
PointIndexType IndexType
Definition: PointPartitioner.h:183
Definition: Exceptions.h:13
Definition: Transform.h:39
PointPartitioner()
Definition: PointPartitioner.h:944
CoordBBox getBBox(size_t n) const
Returns the coordinate-aligned bounding box for bucket n.
Definition: PointPartitioner.h:145
bool usingCellCenteredTransform() const
Returns true if this point partitioning was constructed using a cell-centered transform.
Definition: PointPartitioner.h:159
bool test() const
Definition: PointPartitioner.h:200
Index32 Index
Definition: Types.h:54
IndexIterator indices(size_t n) const
Returns the point indices for bucket n.
Definition: PointPartitioner.h:989
static Ptr create(const PointArray &points, const math::Transform &xform, bool voxelOrder=false, bool recordVoxelOffsets=false, bool cellCenteredTransform=true)
Partitions point indices into BucketLog2Dim aligned buckets.
Definition: PointPartitioner.h:1030
void clear()
Removes all data and frees up memory.
Definition: PointPartitioner.h:957
bool operator!=(const IndexIterator &other) const
Definition: PointPartitioner.h:211
typename std::conditional<(bits< 16), int16_t, typename std::conditional<(bits< 32), int32_t, int64_t >::type >::type VoxelOffsetType
Definition: PointPartitioner.h:89
const IndexType & operator*() const
Definition: PointPartitioner.h:196
SharedPtr< const PointPartitioner > ConstPtr
Definition: PointPartitioner.h:82
void swap(PointPartitioner &)
Exchanges the content of the container by another.
Definition: PointPartitioner.h:970
void reset()
Rewind to first item.
Definition: PointPartitioner.h:189
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
IndexType & operator*()
Returns the item to which this iterator is currently pointing.
Definition: PointPartitioner.h:195
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:473
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
size_t size() const
Number of point indices in the iterator range.
Definition: PointPartitioner.h:192