OpenVDB  8.1.1
PointPartitioner.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
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 
44 
45 
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 
96 
98 
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 
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 
130  size_t size() const { return mPageCount; }
131 
133  bool empty() const { return mPageCount == 0; }
134 
136  void clear();
137 
139  void swap(PointPartitioner&);
140 
142  IndexIterator indices(size_t n) const;
143 
145  CoordBBox getBBox(size_t n) const {
146  return CoordBBox::createCube(mPageCoordinates[n], (1u << BucketLog2Dim));
147  }
148 
150  const Coord& origin(size_t n) const { return mPageCoordinates[n]; }
151 
154  const VoxelOffsetArray& voxelOffsets() const { return mVoxelOffsets; }
155 
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 
189  void reset() { mItem = mBegin; }
190 
192  size_t size() const { return mEnd - mBegin; }
193 
195  IndexType& operator*() { assert(mItem != nullptr); return *mItem; }
196  const IndexType& operator*() const { assert(mItem != nullptr); return *mItem; }
197 
199  operator bool() const { return mItem < mEnd; }
200  bool test() const { return mItem < mEnd; }
201 
203  IndexIterator& operator++() { assert(this->test()); ++mItem; return *this; }
204 
206  bool next() { this->operator++(); return this->test(); }
207  bool increment() { this->next(); return this->test(); }
208 
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 
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) {
241  mPointOrder[n] += mBucketCounters[mBucketOffsets[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) {
264  mOrderedIndexArray[mPointOrder[n]] = mIndices[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 
285  VoxelOrderOp(IndexArray& indices, const IndexArray& pages,const VoxelOffsetArray& offsets)
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 
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) :
606  mXForm.worldToIndexNodeCentered(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.
727  tbb::parallel_for(segmentRange, ComputePointOrderOp<PointIndexType>(
728  bucketIndices.get(), bucketCounters.get(), offsets));
729 
730  tbb::parallel_for(segmentRange, CreateOrderedPointIndexArrayOp<PointIndexType>(
731  offsets, bucketIndices.get(), indices));
732 
733  mIndexSegments[n]->clear(); // clear data
734  }
735  }
736 
742 }; // struct OrderSegmentsOp
743 
744 
746 
747 
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 
848  tbb::parallel_for(segmentRange, OrderSegmentsOp<PointIndexType>
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 
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) {
1015  tbb::parallel_for(pageRange, point_partitioner_internal::VoxelOrderOp<
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 
1042 
1043 
1044 } // namespace tools
1045 } // namespace OPENVDB_VERSION_NAME
1046 } // namespace openvdb
1047 
1048 
1049 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:292
IndexPairListMapPtr *const mData
Definition: PointPartitioner.h:645
std::unique_ptr< VoxelOffsetType[]> VoxelOffsetArray
Definition: PointPartitioner.h:91
Index const mBinVolume
Definition: PointPartitioner.h:741
IndexArray *const mPageIndexArrays
Definition: PointPartitioner.h:740
IndexIterator & operator++()
Advance to the next item.
Definition: PointPartitioner.h:203
SharedPtr< const PointPartitioner > ConstPtr
Definition: PointPartitioner.h:82
OrderSegmentsOp(SegmentPtr *indexSegments, SegmentPtr *offsetSegments, IndexArray *pageOffsetArrays, IndexArray *pageIndexArrays, Index binVolume)
Definition: PointPartitioner.h:662
PointIndexType const *const mBucketOffsets
Definition: PointPartitioner.h:247
typename std::conditional<(bits< 16), int16_t, typename std::conditional<(bits< 32), int32_t, int64_t >::type >::type VoxelOffsetType
Definition: PointPartitioner.h:280
PointIndexType const *const mBucketCounters
Definition: PointPartitioner.h:246
std::unique_ptr< PointIndexType[]> IndexArray
Definition: PointPartitioner.h:659
const VoxelOffsetArray & voxelOffsets() const
Returns a list of LeafNode voxel offsets for the points.
Definition: PointPartitioner.h:154
void clear()
Definition: PointPartitioner.h:363
IndexIterator(IndexType *begin=nullptr, IndexType *end=nullptr)
Definition: PointPartitioner.h:185
std::pair< PointIndexType, PointIndexType > IndexPair
Definition: PointPartitioner.h:420
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:381
std::shared_ptr< IndexPairListMap > IndexPairListMapPtr
Definition: PointPartitioner.h:544
const IndexType & operator*() const
Definition: PointPartitioner.h:196
PointIndexType *const mIndices
Definition: PointPartitioner.h:342
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:239
SegmentPtr *const mIndexSegments
Definition: PointPartitioner.h:737
Array(size_t size)
Definition: PointPartitioner.h:356
bool test() const
Definition: PointPartitioner.h:200
size_t size() const
Returns the number of buckets.
Definition: PointPartitioner.h:130
BinPointIndicesOp(IndexPairListMapPtr *data, const PointArray &points, VoxelOffsetType *voxelOffsets, const math::Transform &m, Index binLog2Dim, Index bucketLog2Dim, size_t numSegments, bool cellCenteredTransform)
Definition: PointPartitioner.h:546
PointIndexType const *const mPages
Definition: PointPartitioner.h:343
CoordBBox getBBox(size_t n) const
Returns the coordinate-aligned bounding box for bucket n.
Definition: PointPartitioner.h:145
std::pair< PointIndexType, PointIndexType > IndexPair
Definition: PointPartitioner.h:540
bool operator==(const IndexIterator &other) const
Equality operators.
Definition: PointPartitioner.h:210
std::shared_ptr< IndexPairList > IndexPairListPtr
Definition: PointPartitioner.h:422
void reset()
Rewind to first item.
Definition: PointPartitioner.h:189
Index32 Index
Definition: openvdb/Types.h:50
std::map< Coord, IndexPairListPtr > IndexPairListMap
Definition: PointPartitioner.h:423
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:672
std::shared_ptr< T > SharedPtr
Definition: openvdb/Types.h:110
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:999
std::unique_ptr< Array > Ptr
Definition: PointPartitioner.h:354
SegmentPtr *const mOffsetSegments
Definition: PointPartitioner.h:738
IndexType & operator*()
Returns the item to which this iterator is currently pointing.
Definition: PointPartitioner.h:195
PointIndexType const *const mPointOrder
Definition: PointPartitioner.h:269
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)
Definition: PointPartitioner.h:806
void clear()
Removes all data and frees up memory.
Definition: PointPartitioner.h:955
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
size_t size() const
Number of point indices in the iterator range.
Definition: PointPartitioner.h:192
void swap(PointPartitioner &)
Exchanges the content of the container by another.
Definition: PointPartitioner.h:968
bool usingCellCenteredTransform() const
Returns true if this point partitioning was constructed using a cell-centered transform.
Definition: PointPartitioner.h:159
std::deque< IndexPair > IndexPairList
Definition: PointPartitioner.h:421
Index const mBinLog2Dim
Definition: PointPartitioner.h:649
typename std::conditional<(bits< 16), int16_t, typename std::conditional<(bits< 32), int32_t, int64_t >::type >::type VoxelOffsetType
Definition: PointPartitioner.h:89
VoxelOffsetType const *const mVoxelOffsets
Definition: PointPartitioner.h:344
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:565
CreateOrderedPointIndexArrayOp(PointIndexType *orderedIndexArray, const PointIndexType *pointOrder, const PointIndexType *indices)
Definition: PointPartitioner.h:254
const Coord & origin(size_t n) const
Returns the origin coordinate for bucket n.
Definition: PointPartitioner.h:150
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:1028
Definition: openvdb/Exceptions.h:13
typename Array< PointIndexType >::Ptr SegmentPtr
Definition: PointPartitioner.h:660
IndexArray *const mPageOffsetArrays
Definition: PointPartitioner.h:739
typename Segment::Ptr SegmentPtr
Definition: PointPartitioner.h:418
VoxelOrderOp(IndexArray &indices, const IndexArray &pages, const VoxelOffsetArray &offsets)
Definition: PointPartitioner.h:285
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
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:439
bool operator!=(const IndexIterator &other) const
Definition: PointPartitioner.h:211
std::vector< Index > IndexArray
Definition: PointMove.h:161
ComputePointOrderOp(PointIndexType *pointOrder, const PointIndexType *bucketCounters, const PointIndexType *bucketOffsets)
Definition: PointPartitioner.h:231
bool const mCellCenteredTransform
Definition: PointPartitioner.h:652
typename Array< PointIndexType >::Ptr SegmentPtr
Definition: PointPartitioner.h:374
MergeBinsOp(IndexPairListMapPtr *bins, SegmentPtr *indexSegments, SegmentPtr *offsetSegments, Coord *coords, size_t numSegments)
Definition: PointPartitioner.h:426
MoveSegmentDataOp(std::vector< PointIndexType * > &indexLists, SegmentPtr *segments)
Definition: PointPartitioner.h:376
std::pair< Index, Index > IndexPair
Definition: PointMove.h:168
typename PointArray::PosType PosType
Definition: PointPartitioner.h:539
VoxelOffsetType *const mVoxelOffsets
Definition: PointPartitioner.h:647
bool increment()
Definition: PointPartitioner.h:207
SharedPtr< PointPartitioner > Ptr
Definition: PointPartitioner.h:81
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:477
std::map< Coord, IndexPairListPtr > IndexPairListMap
Definition: PointPartitioner.h:543
math::Transform const mXForm
Definition: PointPartitioner.h:648
PointIndexType *const mPointOrder
Definition: PointPartitioner.h:245
PointIndexType const *const mIndices
Definition: PointPartitioner.h:270
IndexIterator indices(size_t n) const
Returns the point indices for bucket n.
Definition: PointPartitioner.h:987
PointPartitioner()
Definition: PointPartitioner.h:942
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:262
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.
Definition: PointPartitioner.h:750
PointArray const *const mPoints
Definition: PointPartitioner.h:646
const T * data() const
Definition: PointPartitioner.h:361
size_t size() const
Definition: PointPartitioner.h:358
PointIndexType IndexType
Definition: PointPartitioner.h:84
std::unique_ptr< PointIndexType[]> IndexArray
Definition: PointPartitioner.h:283
PointIndexType *const mOrderedIndexArray
Definition: PointPartitioner.h:268
std::shared_ptr< IndexPairListMap > IndexPairListMapPtr
Definition: PointPartitioner.h:424
size_t const mNumSegments
Definition: PointPartitioner.h:651
std::deque< IndexPair > IndexPairList
Definition: PointPartitioner.h:541
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
std::shared_ptr< IndexPairList > IndexPairListPtr
Definition: PointPartitioner.h:542
T * data()
Definition: PointPartitioner.h:360
PointIndexType IndexType
Definition: PointPartitioner.h:183
Definition: PointPartitioner.h:76
bool next()
Advance to the next item.
Definition: PointPartitioner.h:206
Index const mBucketLog2Dim
Definition: PointPartitioner.h:650
bool empty() const
true if the container size is 0, false otherwise.
Definition: PointPartitioner.h:133
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
Definition: Transform.h:39
std::unique_ptr< VoxelOffsetType[]> VoxelOffsetArray
Definition: PointPartitioner.h:282
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