HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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_scheduler_init.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 
225 namespace point_partitioner_internal {
226 
227 
228 template<typename PointIndexType>
230 {
231  ComputePointOrderOp(PointIndexType* pointOrder,
232  const PointIndexType* bucketCounters, const PointIndexType* bucketOffsets)
233  : mPointOrder(pointOrder)
234  , mBucketCounters(bucketCounters)
235  , mBucketOffsets(bucketOffsets)
236  {
237  }
238 
239  void operator()(const tbb::blocked_range<size_t>& range) const {
240  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
242  }
243  }
244 
245  PointIndexType * const mPointOrder;
246  PointIndexType const * const mBucketCounters;
247  PointIndexType const * const mBucketOffsets;
248 }; // struct ComputePointOrderOp
249 
250 
251 template<typename PointIndexType>
253 {
254  CreateOrderedPointIndexArrayOp(PointIndexType* orderedIndexArray,
255  const PointIndexType* pointOrder, const PointIndexType* indices)
256  : mOrderedIndexArray(orderedIndexArray)
257  , mPointOrder(pointOrder)
258  , mIndices(indices)
259  {
260  }
261 
262  void operator()(const tbb::blocked_range<size_t>& range) const {
263  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
265  }
266  }
267 
268  PointIndexType * const mOrderedIndexArray;
269  PointIndexType const * const mPointOrder;
270  PointIndexType const * const mIndices;
271 }; // struct CreateOrderedPointIndexArrayOp
272 
273 
274 template<typename PointIndexType, Index BucketLog2Dim>
276 {
277  static constexpr Index bits = 1 + (3 * BucketLog2Dim);
278  // signed, so if bits is exactly 16, int32 is required
279  using VoxelOffsetType = typename std::conditional<(bits < 16),
280  int16_t, typename std::conditional<(bits < 32), int32_t, int64_t>::type>::type;
281 
282  using VoxelOffsetArray = std::unique_ptr<VoxelOffsetType[]>;
283  using IndexArray = std::unique_ptr<PointIndexType[]>;
284 
286  : mIndices(indices.get())
287  , mPages(pages.get())
288  , mVoxelOffsets(offsets.get())
289  {
290  }
291 
292  void operator()(const tbb::blocked_range<size_t>& range) const {
293 
294  PointIndexType pointCount = 0;
295  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
296  pointCount = std::max(pointCount, (mPages[n + 1] - mPages[n]));
297  }
298 
299  const PointIndexType voxelCount = 1 << (3 * BucketLog2Dim);
300 
301  // allocate histogram buffers
302  std::unique_ptr<VoxelOffsetType[]> offsets(new VoxelOffsetType[pointCount]);
303  std::unique_ptr<PointIndexType[]> sortedIndices(new PointIndexType[pointCount]);
304  std::unique_ptr<PointIndexType[]> histogram(new PointIndexType[voxelCount]);
305 
306  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
307 
308  PointIndexType * const indices = mIndices + mPages[n];
309  pointCount = mPages[n + 1] - mPages[n];
310 
311  // local copy of voxel offsets.
312  for (PointIndexType i = 0; i < pointCount; ++i) {
313  offsets[i] = mVoxelOffsets[ indices[i] ];
314  }
315 
316  // reset histogram
317  memset(&histogram[0], 0, voxelCount * sizeof(PointIndexType));
318 
319  // compute histogram
320  for (PointIndexType i = 0; i < pointCount; ++i) {
321  ++histogram[ offsets[i] ];
322  }
323 
324  PointIndexType count = 0, startOffset;
325  for (int i = 0; i < int(voxelCount); ++i) {
326  if (histogram[i] > 0) {
327  startOffset = count;
328  count += histogram[i];
329  histogram[i] = startOffset;
330  }
331  }
332 
333  // sort indices based on voxel offset
334  for (PointIndexType i = 0; i < pointCount; ++i) {
335  sortedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
336  }
337 
338  memcpy(&indices[0], &sortedIndices[0], sizeof(PointIndexType) * pointCount);
339  }
340  }
341 
342  PointIndexType * const mIndices;
343  PointIndexType const * const mPages;
345 }; // struct VoxelOrderOp
346 
347 
348 ////////////////////////////////////////
349 
350 
351 template<typename T>
352 struct Array
353 {
354  using Ptr = std::unique_ptr<Array>;
355 
356  Array(size_t size) : mSize(size), mData(new T[size]) { }
357 
358  size_t size() const { return mSize; }
359 
360  T* data() { return mData.get(); }
361  const T* data() const { return mData.get(); }
362 
363  void clear() { mSize = 0; mData.reset(); }
364 
365 private:
366  size_t mSize;
367  std::unique_ptr<T[]> mData;
368 }; // struct Array
369 
370 
371 template<typename PointIndexType>
373 {
375 
376  MoveSegmentDataOp(std::vector<PointIndexType*>& indexLists, SegmentPtr* segments)
377  : mIndexLists(&indexLists[0]), mSegments(segments)
378  {
379  }
380 
381  void operator()(const tbb::blocked_range<size_t>& range) const {
382  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
383  PointIndexType* indices = mIndexLists[n];
384  SegmentPtr& segment = mSegments[n];
385 
386  tbb::parallel_for(tbb::blocked_range<size_t>(0, segment->size()),
387  CopyData(indices, segment->data()));
388 
389  segment.reset(); // clear data
390  }
391  }
392 
393 private:
394 
395  struct CopyData
396  {
397  CopyData(PointIndexType* lhs, const PointIndexType* rhs) : mLhs(lhs), mRhs(rhs) { }
398 
399  void operator()(const tbb::blocked_range<size_t>& range) const {
400  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
401  mLhs[n] = mRhs[n];
402  }
403  }
404 
405  PointIndexType * const mLhs;
406  PointIndexType const * const mRhs;
407  };
408 
409  PointIndexType * const * const mIndexLists;
410  SegmentPtr * const mSegments;
411 }; // struct MoveSegmentDataOp
412 
413 
414 template<typename PointIndexType>
416 {
418  using SegmentPtr = typename Segment::Ptr;
419 
420  using IndexPair = std::pair<PointIndexType, PointIndexType>;
421  using IndexPairList = std::deque<IndexPair>;
422  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
423  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
424  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
425 
427  SegmentPtr* indexSegments,
428  SegmentPtr* offsetSegments,
429  Coord* coords,
430  size_t numSegments)
431  : mBins(bins)
432  , mIndexSegments(indexSegments)
433  , mOffsetSegments(offsetSegments)
434  , mCoords(coords)
435  , mNumSegments(numSegments)
436  {
437  }
438 
439  void operator()(const tbb::blocked_range<size_t>& range) const {
440 
441  std::vector<IndexPairListPtr*> data;
442  std::vector<PointIndexType> arrayOffsets;
443 
444  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
445 
446  const Coord& ijk = mCoords[n];
447  size_t numIndices = 0;
448 
449  data.clear();
450 
451  for (size_t i = 0, I = mNumSegments; i < I; ++i) {
452 
453  IndexPairListMap& idxMap = *mBins[i];
454  typename IndexPairListMap::iterator iter = idxMap.find(ijk);
455 
456  if (iter != idxMap.end() && iter->second) {
457  IndexPairListPtr& idxListPtr = iter->second;
458 
459  data.push_back(&idxListPtr);
460  numIndices += idxListPtr->size();
461  }
462  }
463 
464  if (data.empty() || numIndices == 0) continue;
465 
466  SegmentPtr& indexSegment = mIndexSegments[n];
467  SegmentPtr& offsetSegment = mOffsetSegments[n];
468 
469  indexSegment.reset(new Segment(numIndices));
470  offsetSegment.reset(new Segment(numIndices));
471 
472  arrayOffsets.clear();
473  arrayOffsets.reserve(data.size());
474 
475  for (size_t i = 0, count = 0, I = data.size(); i < I; ++i) {
476  arrayOffsets.push_back(PointIndexType(count));
477  count += (*data[i])->size();
478  }
479 
480  tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size()),
481  CopyData(&data[0], &arrayOffsets[0], indexSegment->data(), offsetSegment->data()));
482  }
483  }
484 
485 private:
486 
487  struct CopyData
488  {
489  CopyData(IndexPairListPtr** indexLists,
490  const PointIndexType* arrayOffsets,
491  PointIndexType* indices,
492  PointIndexType* offsets)
493  : mIndexLists(indexLists)
494  , mArrayOffsets(arrayOffsets)
495  , mIndices(indices)
496  , mOffsets(offsets)
497  {
498  }
499 
500  void operator()(const tbb::blocked_range<size_t>& range) const {
501 
502  using CIter = typename IndexPairList::const_iterator;
503 
504  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
505 
506  const PointIndexType arrayOffset = mArrayOffsets[n];
507  PointIndexType* indexPtr = &mIndices[arrayOffset];
508  PointIndexType* offsetPtr = &mOffsets[arrayOffset];
509 
510  IndexPairListPtr& list = *mIndexLists[n];
511 
512  for (CIter it = list->begin(), end = list->end(); it != end; ++it) {
513  const IndexPair& data = *it;
514  *indexPtr++ = data.first;
515  *offsetPtr++ = data.second;
516  }
517 
518  list.reset(); // clear data
519  }
520  }
521 
522  IndexPairListPtr * const * const mIndexLists;
523  PointIndexType const * const mArrayOffsets;
524  PointIndexType * const mIndices;
525  PointIndexType * const mOffsets;
526  }; // struct CopyData
527 
528  IndexPairListMapPtr * const mBins;
529  SegmentPtr * const mIndexSegments;
530  SegmentPtr * const mOffsetSegments;
531  Coord const * const mCoords;
532  size_t const mNumSegments;
533 }; // struct MergeBinsOp
534 
535 
536 template<typename PointArray, typename PointIndexType, typename VoxelOffsetType>
538 {
539  using PosType = typename PointArray::PosType;
540  using IndexPair = std::pair<PointIndexType, PointIndexType>;
541  using IndexPairList = std::deque<IndexPair>;
542  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
543  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
544  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
545 
547  const PointArray& points,
548  VoxelOffsetType* voxelOffsets,
549  const math::Transform& m,
550  Index binLog2Dim,
551  Index bucketLog2Dim,
552  size_t numSegments,
553  bool cellCenteredTransform)
554  : mData(data)
555  , mPoints(&points)
556  , mVoxelOffsets(voxelOffsets)
557  , mXForm(m)
558  , mBinLog2Dim(binLog2Dim)
559  , mBucketLog2Dim(bucketLog2Dim)
560  , mNumSegments(numSegments)
561  , mCellCenteredTransform(cellCenteredTransform)
562  {
563  }
564 
565  void operator()(const tbb::blocked_range<size_t>& range) const {
566 
567  const Index log2dim = mBucketLog2Dim;
568  const Index log2dim2 = 2 * log2dim;
569  const Index bucketMask = (1u << log2dim) - 1u;
570 
571  const Index binLog2dim = mBinLog2Dim;
572  const Index binLog2dim2 = 2 * binLog2dim;
573 
574  const Index binMask = (1u << (log2dim + binLog2dim)) - 1u;
575  const Index invBinMask = ~binMask;
576 
577  IndexPairList * idxList = nullptr;
578  Coord ijk(0, 0, 0), loc(0, 0, 0), binCoord(0, 0, 0), lastBinCoord(1, 2, 3);
579  PosType pos;
580 
581  PointIndexType bucketOffset = 0;
582  VoxelOffsetType voxelOffset = 0;
583 
584  const bool cellCentered = mCellCenteredTransform;
585 
586  const size_t numPoints = mPoints->size();
587  const size_t segmentSize = numPoints / mNumSegments;
588 
589  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
590 
591  IndexPairListMapPtr& dataPtr = mData[n];
592  if (!dataPtr) dataPtr.reset(new IndexPairListMap());
593  IndexPairListMap& idxMap = *dataPtr;
594 
595  const bool isLastSegment = (n + 1) >= mNumSegments;
596 
597  const size_t start = n * segmentSize;
598  const size_t end = isLastSegment ? numPoints : (start + segmentSize);
599 
600  for (size_t i = start; i != end; ++i) {
601 
602  mPoints->getPos(i, pos);
603 
604  if (std::isfinite(pos[0]) && std::isfinite(pos[1]) && std::isfinite(pos[2])) {
605  ijk = cellCentered ? mXForm.worldToIndexCellCentered(pos) :
607 
608  if (mVoxelOffsets) {
609  loc[0] = ijk[0] & bucketMask;
610  loc[1] = ijk[1] & bucketMask;
611  loc[2] = ijk[2] & bucketMask;
612  voxelOffset = VoxelOffsetType(
613  (loc[0] << log2dim2) + (loc[1] << log2dim) + loc[2]);
614  }
615 
616  binCoord[0] = ijk[0] & invBinMask;
617  binCoord[1] = ijk[1] & invBinMask;
618  binCoord[2] = ijk[2] & invBinMask;
619 
620  ijk[0] &= binMask;
621  ijk[1] &= binMask;
622  ijk[2] &= binMask;
623 
624  ijk[0] >>= log2dim;
625  ijk[1] >>= log2dim;
626  ijk[2] >>= log2dim;
627 
628  bucketOffset = PointIndexType(
629  (ijk[0] << binLog2dim2) + (ijk[1] << binLog2dim) + ijk[2]);
630 
631  if (lastBinCoord != binCoord) {
632  lastBinCoord = binCoord;
633  IndexPairListPtr& idxListPtr = idxMap[lastBinCoord];
634  if (!idxListPtr) idxListPtr.reset(new IndexPairList());
635  idxList = idxListPtr.get();
636  }
637 
638  idxList->push_back(IndexPair(PointIndexType(i), bucketOffset));
639  if (mVoxelOffsets) mVoxelOffsets[i] = voxelOffset;
640  }
641  }
642  }
643  }
644 
646  PointArray const * const mPoints;
647  VoxelOffsetType * const mVoxelOffsets;
651  size_t const mNumSegments;
653 }; // struct BinPointIndicesOp
654 
655 
656 template<typename PointIndexType>
658 {
659  using IndexArray = std::unique_ptr<PointIndexType[]>;
661 
662  OrderSegmentsOp(SegmentPtr* indexSegments, SegmentPtr* offsetSegments,
663  IndexArray* pageOffsetArrays, IndexArray* pageIndexArrays, Index binVolume)
664  : mIndexSegments(indexSegments)
665  , mOffsetSegments(offsetSegments)
666  , mPageOffsetArrays(pageOffsetArrays)
667  , mPageIndexArrays(pageIndexArrays)
668  , mBinVolume(binVolume)
669  {
670  }
671 
672  void operator()(const tbb::blocked_range<size_t>& range) const {
673 
674  const size_t bucketCountersSize = size_t(mBinVolume);
675  IndexArray bucketCounters(new PointIndexType[bucketCountersSize]);
676 
677  size_t maxSegmentSize = 0;
678  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
679  maxSegmentSize = std::max(maxSegmentSize, mIndexSegments[n]->size());
680  }
681 
682  IndexArray bucketIndices(new PointIndexType[maxSegmentSize]);
683 
684  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
685 
686  memset(bucketCounters.get(), 0, sizeof(PointIndexType) * bucketCountersSize);
687 
688  const size_t segmentSize = mOffsetSegments[n]->size();
689  PointIndexType* offsets = mOffsetSegments[n]->data();
690 
691  // Count the number of points per bucket and assign a local bucket index
692  // to each point.
693  for (size_t i = 0; i < segmentSize; ++i) {
694  bucketIndices[i] = bucketCounters[offsets[i]]++;
695  }
696 
697  PointIndexType nonemptyBucketCount = 0;
698  for (size_t i = 0; i < bucketCountersSize; ++i) {
699  nonemptyBucketCount += static_cast<PointIndexType>(bucketCounters[i] != 0);
700  }
701 
702 
703  IndexArray& pageOffsets = mPageOffsetArrays[n];
704  pageOffsets.reset(new PointIndexType[nonemptyBucketCount + 1]);
705  pageOffsets[0] = nonemptyBucketCount + 1; // stores array size in first element
706 
707  IndexArray& pageIndices = mPageIndexArrays[n];
708  pageIndices.reset(new PointIndexType[nonemptyBucketCount]);
709 
710  // Compute bucket counter prefix sum
711  PointIndexType count = 0, idx = 0;
712  for (size_t i = 0; i < bucketCountersSize; ++i) {
713  if (bucketCounters[i] != 0) {
714  pageIndices[idx] = static_cast<PointIndexType>(i);
715  pageOffsets[idx+1] = bucketCounters[i];
716  bucketCounters[i] = count;
717  count += pageOffsets[idx+1];
718  ++idx;
719  }
720  }
721 
722  PointIndexType* indices = mIndexSegments[n]->data();
723  const tbb::blocked_range<size_t> segmentRange(0, segmentSize);
724 
725  // Compute final point order by incrementing the local bucket point index
726  // with the prefix sum offset.
728  bucketIndices.get(), bucketCounters.get(), offsets));
729 
731  offsets, bucketIndices.get(), indices));
732 
733  mIndexSegments[n]->clear(); // clear data
734  }
735  }
736 
742 }; // struct OrderSegmentsOp
743 
744 
745 ////////////////////////////////////////
746 
747 
748 /// @brief Segment points using one level of least significant digit radix bins.
749 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
750 inline void binAndSegment(
751  const PointArray& points,
752  const math::Transform& xform,
753  std::unique_ptr<typename Array<PointIndexType>::Ptr[]>& indexSegments,
754  std::unique_ptr<typename Array<PointIndexType>::Ptr[]>& offsetSegments,
755  std::vector<Coord>& coords,
756  const Index binLog2Dim,
757  const Index bucketLog2Dim,
758  VoxelOffsetType* voxelOffsets = nullptr,
759  bool cellCenteredTransform = true)
760 {
761  using IndexPair = std::pair<PointIndexType, PointIndexType>;
762  using IndexPairList = std::deque<IndexPair>;
763  using IndexPairListPtr = std::shared_ptr<IndexPairList>;
764  using IndexPairListMap = std::map<Coord, IndexPairListPtr>;
765  using IndexPairListMapPtr = std::shared_ptr<IndexPairListMap>;
766 
767  size_t numTasks = 1, numThreads = size_t(tbb::task_scheduler_init::default_num_threads());
768  if (points.size() > (numThreads * 2)) numTasks = numThreads * 2;
769  else if (points.size() > numThreads) numTasks = numThreads;
770 
771  std::unique_ptr<IndexPairListMapPtr[]> bins(new IndexPairListMapPtr[numTasks]);
772 
774 
775  tbb::parallel_for(tbb::blocked_range<size_t>(0, numTasks),
776  BinOp(bins.get(), points, voxelOffsets, xform, binLog2Dim, bucketLog2Dim,
777  numTasks, cellCenteredTransform));
778 
779  std::set<Coord> uniqueCoords;
780 
781  for (size_t i = 0; i < numTasks; ++i) {
782  IndexPairListMap& idxMap = *bins[i];
783  for (typename IndexPairListMap::iterator it = idxMap.begin(); it != idxMap.end(); ++it) {
784  uniqueCoords.insert(it->first);
785  }
786  }
787 
788  coords.assign(uniqueCoords.begin(), uniqueCoords.end());
789  uniqueCoords.clear();
790 
791  size_t segmentCount = coords.size();
792 
793  using SegmentPtr = typename Array<PointIndexType>::Ptr;
794 
795  indexSegments.reset(new SegmentPtr[segmentCount]);
796  offsetSegments.reset(new SegmentPtr[segmentCount]);
797 
798  using MergeOp = MergeBinsOp<PointIndexType>;
799 
800  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
801  MergeOp(bins.get(), indexSegments.get(), offsetSegments.get(), &coords[0], numTasks));
802 }
803 
804 
805 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
806 inline void partition(
807  const PointArray& points,
808  const math::Transform& xform,
809  const Index bucketLog2Dim,
810  std::unique_ptr<PointIndexType[]>& pointIndices,
811  std::unique_ptr<PointIndexType[]>& pageOffsets,
812  std::unique_ptr<Coord[]>& pageCoordinates,
813  PointIndexType& pageCount,
814  std::unique_ptr<VoxelOffsetType[]>& voxelOffsets,
815  bool recordVoxelOffsets,
816  bool cellCenteredTransform)
817 {
818  using SegmentPtr = typename Array<PointIndexType>::Ptr;
819 
820  if (recordVoxelOffsets) voxelOffsets.reset(new VoxelOffsetType[points.size()]);
821  else voxelOffsets.reset();
822 
823  const Index binLog2Dim = 5u;
824  // note: Bins span a (2^(binLog2Dim + bucketLog2Dim))^3 voxel region,
825  // i.e. bucketLog2Dim = 3 and binLog2Dim = 5 corresponds to a
826  // (2^8)^3 = 256^3 voxel region.
827 
828 
829  std::vector<Coord> segmentCoords;
830 
831  std::unique_ptr<SegmentPtr[]> indexSegments;
832  std::unique_ptr<SegmentPtr[]> offsetSegments;
833 
834  binAndSegment<PointIndexType, VoxelOffsetType, PointArray>(points, xform,
835  indexSegments, offsetSegments, segmentCoords, binLog2Dim, bucketLog2Dim,
836  voxelOffsets.get(), cellCenteredTransform);
837 
838  size_t numSegments = segmentCoords.size();
839 
840  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
841 
842  using IndexArray = std::unique_ptr<PointIndexType[]>;
843  std::unique_ptr<IndexArray[]> pageOffsetArrays(new IndexArray[numSegments]);
844  std::unique_ptr<IndexArray[]> pageIndexArrays(new IndexArray[numSegments]);
845 
846  const Index binVolume = 1u << (3u * binLog2Dim);
847 
849  (indexSegments.get(), offsetSegments.get(),
850  pageOffsetArrays.get(), pageIndexArrays.get(), binVolume));
851 
852  indexSegments.reset();
853 
854  std::vector<Index> segmentOffsets;
855  segmentOffsets.reserve(numSegments);
856 
857  pageCount = 0;
858  for (size_t n = 0; n < numSegments; ++n) {
859  segmentOffsets.push_back(pageCount);
860  pageCount += pageOffsetArrays[n][0] - 1;
861  }
862 
863  pageOffsets.reset(new PointIndexType[pageCount + 1]);
864 
865  PointIndexType count = 0;
866  for (size_t n = 0, idx = 0; n < numSegments; ++n) {
867 
868  PointIndexType* offsets = pageOffsetArrays[n].get();
869  size_t size = size_t(offsets[0]);
870 
871  for (size_t i = 1; i < size; ++i) {
872  pageOffsets[idx++] = count;
873  count += offsets[i];
874  }
875  }
876 
877  pageOffsets[pageCount] = count;
878 
879  pointIndices.reset(new PointIndexType[points.size()]);
880 
881  std::vector<PointIndexType*> indexArray;
882  indexArray.reserve(numSegments);
883 
884  PointIndexType* index = pointIndices.get();
885  for (size_t n = 0; n < numSegments; ++n) {
886  indexArray.push_back(index);
887  index += offsetSegments[n]->size();
888  }
889 
890  // compute leaf node origin for each page
891 
892  pageCoordinates.reset(new Coord[pageCount]);
893 
894  tbb::parallel_for(segmentRange,
895  [&](tbb::blocked_range<size_t>& range)
896  {
897  for (size_t n = range.begin(); n < range.end(); n++)
898  {
899  Index segmentOffset = segmentOffsets[n];
900  PointIndexType* indices = pageIndexArrays[n].get();
901 
902  const Coord& segmentCoord = segmentCoords[n];
903 
904  // segment size stored in the first value of the offset array
905  const size_t segmentSize = pageOffsetArrays[n][0] - 1;
906  tbb::blocked_range<size_t> copyRange(0, segmentSize);
907  tbb::parallel_for(copyRange,
908  [&](tbb::blocked_range<size_t>& r)
909  {
910  for (size_t i = r.begin(); i < r.end(); i++)
911  {
912  Index pageIndex = indices[i];
913  Coord& ijk = pageCoordinates[segmentOffset+i];
914 
915  ijk[0] = pageIndex >> (2 * binLog2Dim);
916  Index pageIndexModulo = pageIndex - (ijk[0] << (2 * binLog2Dim));
917  ijk[1] = pageIndexModulo >> binLog2Dim;
918  ijk[2] = pageIndexModulo - (ijk[1] << binLog2Dim);
919 
920  ijk = (ijk << bucketLog2Dim) + segmentCoord;
921  }
922  }
923  );
924  }
925  }
926  );
927 
928  // move segment data
929 
930  tbb::parallel_for(segmentRange,
931  MoveSegmentDataOp<PointIndexType>(indexArray, offsetSegments.get()));
932 }
933 
934 
935 } // namespace point_partitioner_internal
936 
937 
938 ////////////////////////////////////////
939 
940 
941 template<typename PointIndexType, Index BucketLog2Dim>
943  : mPointIndices(nullptr)
944  , mVoxelOffsets(nullptr)
945  , mPageOffsets(nullptr)
946  , mPageCoordinates(nullptr)
947  , mPageCount(0)
948  , mUsingCellCenteredTransform(true)
949 {
950 }
951 
952 
953 template<typename PointIndexType, Index BucketLog2Dim>
954 inline void
956 {
957  mPageCount = 0;
958  mUsingCellCenteredTransform = true;
959  mPointIndices.reset();
960  mVoxelOffsets.reset();
961  mPageOffsets.reset();
962  mPageCoordinates.reset();
963 }
964 
965 
966 template<typename PointIndexType, Index BucketLog2Dim>
967 inline void
969 {
970  const IndexType tmpLhsPageCount = mPageCount;
971  mPageCount = rhs.mPageCount;
972  rhs.mPageCount = tmpLhsPageCount;
973 
974  mPointIndices.swap(rhs.mPointIndices);
975  mVoxelOffsets.swap(rhs.mVoxelOffsets);
976  mPageOffsets.swap(rhs.mPageOffsets);
977  mPageCoordinates.swap(rhs.mPageCoordinates);
978 
979  bool lhsCellCenteredTransform = mUsingCellCenteredTransform;
980  mUsingCellCenteredTransform = rhs.mUsingCellCenteredTransform;
981  rhs.mUsingCellCenteredTransform = lhsCellCenteredTransform;
982 }
983 
984 
985 template<typename PointIndexType, Index BucketLog2Dim>
988 {
989  assert(bool(mPointIndices) && bool(mPageCount));
990  return IndexIterator(
991  mPointIndices.get() + mPageOffsets[n],
992  mPointIndices.get() + mPageOffsets[n + 1]);
993 }
994 
995 
996 template<typename PointIndexType, Index BucketLog2Dim>
997 template<typename PointArray>
998 inline void
1000  const PointArray& points,
1001  const math::Transform& xform,
1002  bool voxelOrder,
1003  bool recordVoxelOffsets,
1004  bool cellCenteredTransform)
1005 {
1006  mUsingCellCenteredTransform = cellCenteredTransform;
1007 
1008  point_partitioner_internal::partition(points, xform, BucketLog2Dim,
1009  mPointIndices, mPageOffsets, mPageCoordinates, mPageCount, mVoxelOffsets,
1010  (voxelOrder || recordVoxelOffsets), cellCenteredTransform);
1011 
1012  const tbb::blocked_range<size_t> pageRange(0, mPageCount);
1013 
1014  if (mVoxelOffsets && voxelOrder) {
1016  IndexType, BucketLog2Dim>(mPointIndices, mPageOffsets, mVoxelOffsets));
1017  }
1018 
1019  if (mVoxelOffsets && !recordVoxelOffsets) {
1020  mVoxelOffsets.reset();
1021  }
1022 }
1023 
1024 
1025 template<typename PointIndexType, Index BucketLog2Dim>
1026 template<typename PointArray>
1029  const PointArray& points,
1030  const math::Transform& xform,
1031  bool voxelOrder,
1032  bool recordVoxelOffsets,
1033  bool cellCenteredTransform)
1034 {
1035  Ptr ret(new PointPartitioner());
1036  ret->construct(points, xform, voxelOrder, recordVoxelOffsets, cellCenteredTransform);
1037  return ret;
1038 }
1039 
1040 
1041 ////////////////////////////////////////
1042 
1043 
1044 } // namespace tools
1045 } // namespace OPENVDB_VERSION_NAME
1046 } // namespace openvdb
1047 
1048 
1049 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
vint4 max(const vint4 &a, const vint4 &b)
Definition: simd.h:4703
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:153
GLenum GLint * range
Definition: glew.h:3500
ComputePointOrderOp(PointIndexType *pointOrder, const PointIndexType *bucketCounters, const PointIndexType *bucketOffsets)
cvex test(vector P=0;int unbound=3;export float s=0;export vector Cf=0;)
Definition: test.vfl:11
GLsizeiptr size
Definition: glew.h:1681
size_t size() const
Returns the number of buckets.
void operator()(const tbb::blocked_range< size_t > &range) const
std::unique_ptr< VoxelOffsetType[]> VoxelOffsetArray
FMT_CONSTEXPR auto begin(const C &c) -> decltype(c.begin())
Definition: format.h:251
GLuint index
Definition: glew.h:1814
GLuint GLsizei const GLuint const GLintptr * offsets
Definition: glew.h:4117
IndexType & operator*()
Returns the item to which this iterator is currently pointing.
GLenum GLuint coords
Definition: glew.h:7906
bool usingCellCenteredTransform() const
Returns true if this point partitioning was constructed using a cell-centered transform.
OrderSegmentsOp(SegmentPtr *indexSegments, SegmentPtr *offsetSegments, IndexArray *pageOffsetArrays, IndexArray *pageIndexArrays, Index binVolume)
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:166
bool empty() const
true if the container size is 0, false otherwise.
const GLdouble * m
Definition: glew.h:9124
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.
void binAndSegment(const PointArray &points, const math::Transform &xform, std::unique_ptr< typename Array< PointIndexType >::Ptr[]> &indexSegments, std::unique_ptr< typename Array< PointIndexType >::Ptr[]> &offsetSegments, std::vector< Coord > &coords, const Index binLog2Dim, const Index bucketLog2Dim, VoxelOffsetType *voxelOffsets=nullptr, bool cellCenteredTransform=true)
Segment points using one level of least significant digit radix bins.
void operator()(const tbb::blocked_range< size_t > &range) const
GLsizei GLsizei numSegments
Definition: glew.h:13600
VoxelOrderOp(IndexArray &indices, const IndexArray &pages, const VoxelOffsetArray &offsets)
std::shared_ptr< T > SharedPtr
Definition: Types.h:91
SYS_FORCE_INLINE const_iterator end() const
bool operator==(const BaseDimensions< T > &a, const BaseDimensions< Y > &b)
Definition: Dimensions.h:137
CoordBBox getBBox(size_t n) const
Returns the coordinate-aligned bounding box for bucket n.
void partition(const PointArray &points, const math::Transform &xform, const Index bucketLog2Dim, std::unique_ptr< PointIndexType[]> &pointIndices, std::unique_ptr< PointIndexType[]> &pageOffsets, std::unique_ptr< Coord[]> &pageCoordinates, PointIndexType &pageCount, std::unique_ptr< VoxelOffsetType[]> &voxelOffsets, bool recordVoxelOffsets, bool cellCenteredTransform)
MoveSegmentDataOp(std::vector< PointIndexType * > &indexLists, SegmentPtr *segments)
size_t size() const
Number of point indices in the iterator range.
IndexIterator(IndexType *begin=nullptr, IndexType *end=nullptr)
IndexIterator indices(size_t n) const
Returns the point indices for bucket n.
GLint GLenum GLsizei GLint GLsizei const void * data
Definition: glew.h:1379
GLuint GLuint GLsizei GLenum const void * indices
Definition: glew.h:1253
GLuint GLuint end
Definition: glew.h:1253
void clear()
Removes all data and frees up memory.
GLsizei n
Definition: glew.h:4040
typename std::conditional<(bits< 16), int16_t, typename std::conditional<(bits< 32), int32_t, int64_t >::type >::type VoxelOffsetType
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
const Coord & origin(size_t n) const
Returns the origin coordinate for bucket n.
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: PointCount.h:88
Coord worldToIndexNodeCentered(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:112
GLuint GLuint GLsizei GLenum type
Definition: glew.h:1253
GLuint start
Definition: glew.h:1253
bool operator==(const IndexIterator &other) const
Equality operators.
typename std::conditional<(bits< 16), int16_t, typename std::conditional<(bits< 32), int32_t, int64_t >::type >::type VoxelOffsetType
const VoxelOffsetArray & voxelOffsets() const
Returns a list of LeafNode voxel offsets for the points.
std::pair< exint, exint > IndexPair
Definition: USD_Utils.h:222
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glew.h:3446
GLdouble GLdouble GLdouble r
Definition: glew.h:1406
GLuint GLuint GLsizei count
Definition: glew.h:1253
Coord worldToIndexCellCentered(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:111
GA_API const UT_StringHolder N
CreateOrderedPointIndexArrayOp(PointIndexType *orderedIndexArray, const PointIndexType *pointOrder, const PointIndexType *indices)
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.
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:341
void swap(PointPartitioner &)
Exchanges the content of the container by another.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:112
void operator()(const tbb::blocked_range< size_t > &range) const
MergeBinsOp(IndexPairListMapPtr *bins, SegmentPtr *indexSegments, SegmentPtr *offsetSegments, Coord *coords, size_t numSegments)
BinPointIndicesOp(IndexPairListMapPtr *data, const PointArray &points, VoxelOffsetType *voxelOffsets, const math::Transform &m, Index binLog2Dim, Index bucketLog2Dim, size_t numSegments, bool cellCenteredTransform)
type
Definition: core.h:528