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