18 #ifndef OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
19 #define OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
24 #include <tbb/blocked_range.h>
25 #include <tbb/parallel_for.h>
26 #include <tbb/task_scheduler_init.h>
75 template<
typename Po
intIndexType = u
int32_t, Index BucketLog2Dim = 3>
79 enum { LOG2DIM = BucketLog2Dim };
86 static constexpr
Index bits = 1 + (3 * BucketLog2Dim);
89 int16_t,
typename std::conditional<(bits < 32), int32_t, int64_t>::type>::type;
108 template<
typename Po
intArray>
110 bool voxelOrder =
false,
bool recordVoxelOffsets =
false,
111 bool cellCenteredTransform =
true);
123 template<
typename Po
intArray>
125 bool voxelOrder =
false,
bool recordVoxelOffsets =
false,
126 bool cellCenteredTransform =
true);
130 size_t size()
const {
return mPageCount; }
133 bool empty()
const {
return mPageCount == 0; }
142 IndexIterator indices(
size_t n)
const;
150 const Coord&
origin(
size_t n)
const {
return mPageCoordinates[n]; }
166 std::unique_ptr<IndexType[]> mPointIndices;
167 VoxelOffsetArray mVoxelOffsets;
169 std::unique_ptr<IndexType[]> mPageOffsets;
170 std::unique_ptr<Coord[]> mPageCoordinates;
171 IndexType mPageCount;
172 bool mUsingCellCenteredTransform;
179 template<
typename Po
intIndexType, Index BucketLog2Dim>
186 : mBegin(begin), mEnd(end), mItem(begin) {}
192 size_t size()
const {
return mEnd - mBegin; }
199 operator bool()
const {
return mItem < mEnd; }
200 bool test()
const {
return mItem < mEnd; }
206 bool next() { this->operator++();
return this->test(); }
225 namespace point_partitioner_internal {
228 template<
typename Po
intIndexType>
232 const PointIndexType* bucketCounters,
const PointIndexType* bucketOffsets)
233 : mPointOrder(pointOrder)
234 , mBucketCounters(bucketCounters)
235 , mBucketOffsets(bucketOffsets)
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]];
251 template<
typename Po
intIndexType>
255 const PointIndexType* pointOrder,
const PointIndexType* indices)
256 : mOrderedIndexArray(orderedIndexArray)
257 , mPointOrder(pointOrder)
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];
274 template<
typename Po
intIndexType, Index BucketLog2Dim>
277 static constexpr
Index bits = 1 + (3 * BucketLog2Dim);
280 int16_t,
typename std::conditional<(bits < 32), int32_t, int64_t>::type>::type;
286 : mIndices(indices.get())
287 , mPages(pages.get())
288 , mVoxelOffsets(offsets.get())
292 void operator()(
const tbb::blocked_range<size_t>& range)
const {
295 for (
size_t n(range.begin()), N(range.end()); n != N; ++n) {
299 const PointIndexType voxelCount = 1 << (3 * BucketLog2Dim);
303 std::unique_ptr<PointIndexType[]> sortedIndices(
new PointIndexType[
pointCount]);
304 std::unique_ptr<PointIndexType[]>
histogram(
new PointIndexType[voxelCount]);
306 for (
size_t n(range.begin()), N(range.end()); n != N; ++n) {
308 PointIndexType *
const indices = mIndices + mPages[n];
312 for (PointIndexType i = 0; i <
pointCount; ++i) {
313 offsets[i] = mVoxelOffsets[ indices[i] ];
317 memset(&
histogram[0], 0, voxelCount *
sizeof(PointIndexType));
320 for (PointIndexType i = 0; i <
pointCount; ++i) {
324 PointIndexType count = 0, startOffset;
325 for (
int i = 0; i < int(voxelCount); ++i) {
334 for (PointIndexType i = 0; i <
pointCount; ++i) {
335 sortedIndices[
histogram[ offsets[i] ]++ ] = indices[i];
338 memcpy(&indices[0], &sortedIndices[0],
sizeof(PointIndexType) *
pointCount);
354 using Ptr = std::unique_ptr<Array>;
356 Array(
size_t size) : mSize(size), mData(new T[size]) { }
358 size_t size()
const {
return mSize; }
360 T*
data() {
return mData.get(); }
361 const T*
data()
const {
return mData.get(); }
363 void clear() { mSize = 0; mData.reset(); }
367 std::unique_ptr<T[]> mData;
371 template<
typename Po
intIndexType>
377 : mIndexLists(&indexLists[0]), mSegments(segments)
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];
386 tbb::parallel_for(tbb::blocked_range<size_t>(0, segment->size()),
387 CopyData(indices, segment->data()));
397 CopyData(PointIndexType* lhs,
const PointIndexType* rhs) : mLhs(lhs), mRhs(rhs) { }
399 void operator()(
const tbb::blocked_range<size_t>& range)
const {
400 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
405 PointIndexType *
const mLhs;
406 PointIndexType
const *
const mRhs;
409 PointIndexType *
const *
const mIndexLists;
410 SegmentPtr *
const mSegments;
414 template<
typename Po
intIndexType>
420 using IndexPair = std::pair<PointIndexType, PointIndexType>;
432 , mIndexSegments(indexSegments)
433 , mOffsetSegments(offsetSegments)
435 , mNumSegments(numSegments)
439 void operator()(
const tbb::blocked_range<size_t>& range)
const {
441 std::vector<IndexPairListPtr*> data;
442 std::vector<PointIndexType> arrayOffsets;
444 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
446 const Coord& ijk = mCoords[n];
447 size_t numIndices = 0;
451 for (
size_t i = 0, I = mNumSegments; i < I; ++i) {
454 typename IndexPairListMap::iterator iter = idxMap.find(ijk);
456 if (iter != idxMap.end() && iter->second) {
459 data.push_back(&idxListPtr);
460 numIndices += idxListPtr->size();
464 if (data.empty() || numIndices == 0)
continue;
467 SegmentPtr& offsetSegment = mOffsetSegments[n];
472 arrayOffsets.clear();
473 arrayOffsets.reserve(data.size());
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();
480 tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size()),
481 CopyData(&data[0], &arrayOffsets[0], indexSegment->data(), offsetSegment->data()));
489 CopyData(IndexPairListPtr** indexLists,
490 const PointIndexType* arrayOffsets,
491 PointIndexType* indices,
492 PointIndexType* offsets)
493 : mIndexLists(indexLists)
494 , mArrayOffsets(arrayOffsets)
500 void operator()(
const tbb::blocked_range<size_t>& range)
const {
502 using CIter =
typename IndexPairList::const_iterator;
504 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
506 const PointIndexType arrayOffset = mArrayOffsets[n];
507 PointIndexType* indexPtr = &mIndices[arrayOffset];
508 PointIndexType* offsetPtr = &mOffsets[arrayOffset];
510 IndexPairListPtr& list = *mIndexLists[n];
512 for (CIter it = list->begin(), end = list->end(); it != end; ++it) {
514 *indexPtr++ = data.first;
515 *offsetPtr++ = data.second;
522 IndexPairListPtr *
const *
const mIndexLists;
523 PointIndexType
const *
const mArrayOffsets;
524 PointIndexType *
const mIndices;
525 PointIndexType *
const mOffsets;
528 IndexPairListMapPtr *
const mBins;
529 SegmentPtr *
const mIndexSegments;
530 SegmentPtr *
const mOffsetSegments;
531 Coord
const *
const mCoords;
532 size_t const mNumSegments;
536 template<
typename Po
intArray,
typename Po
intIndexType,
typename VoxelOffsetType>
540 using IndexPair = std::pair<PointIndexType, PointIndexType>;
548 VoxelOffsetType* voxelOffsets,
553 bool cellCenteredTransform)
556 , mVoxelOffsets(voxelOffsets)
558 , mBinLog2Dim(binLog2Dim)
559 , mBucketLog2Dim(bucketLog2Dim)
560 , mNumSegments(numSegments)
561 , mCellCenteredTransform(cellCenteredTransform)
565 void operator()(
const tbb::blocked_range<size_t>& range)
const {
567 const Index log2dim = mBucketLog2Dim;
568 const Index log2dim2 = 2 * log2dim;
569 const Index bucketMask = (1u << log2dim) - 1u;
571 const Index binLog2dim = mBinLog2Dim;
572 const Index binLog2dim2 = 2 * binLog2dim;
574 const Index binMask = (1u << (log2dim + binLog2dim)) - 1u;
575 const Index invBinMask = ~binMask;
578 Coord ijk(0, 0, 0), loc(0, 0, 0), binCoord(0, 0, 0), lastBinCoord(1, 2, 3);
581 PointIndexType bucketOffset = 0;
582 VoxelOffsetType voxelOffset = 0;
584 const bool cellCentered = mCellCenteredTransform;
586 const size_t numPoints = mPoints->size();
587 const size_t segmentSize = numPoints / mNumSegments;
589 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
595 const bool isLastSegment = (n + 1) >= mNumSegments;
597 const size_t start = n * segmentSize;
598 const size_t end = isLastSegment ? numPoints : (start + segmentSize);
600 for (
size_t i = start; i != end; ++i) {
602 mPoints->getPos(i, pos);
604 if (std::isfinite(pos[0]) && std::isfinite(pos[1]) && std::isfinite(pos[2])) {
605 ijk = cellCentered ? mXForm.worldToIndexCellCentered(pos) :
606 mXForm.worldToIndexNodeCentered(pos);
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]);
616 binCoord[0] = ijk[0] & invBinMask;
617 binCoord[1] = ijk[1] & invBinMask;
618 binCoord[2] = ijk[2] & invBinMask;
628 bucketOffset = PointIndexType(
629 (ijk[0] << binLog2dim2) + (ijk[1] << binLog2dim) + ijk[2]);
631 if (lastBinCoord != binCoord) {
632 lastBinCoord = binCoord;
635 idxList = idxListPtr.get();
638 idxList->push_back(
IndexPair(PointIndexType(i), bucketOffset));
639 if (mVoxelOffsets) mVoxelOffsets[i] = voxelOffset;
656 template<
typename Po
intIndexType>
664 : mIndexSegments(indexSegments)
665 , mOffsetSegments(offsetSegments)
666 , mPageOffsetArrays(pageOffsetArrays)
667 , mPageIndexArrays(pageIndexArrays)
668 , mBinVolume(binVolume)
672 void operator()(
const tbb::blocked_range<size_t>& range)
const {
674 const size_t bucketCountersSize = size_t(mBinVolume);
675 IndexArray bucketCounters(
new PointIndexType[bucketCountersSize]);
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());
682 IndexArray bucketIndices(
new PointIndexType[maxSegmentSize]);
684 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
686 memset(bucketCounters.get(), 0,
sizeof(PointIndexType) * bucketCountersSize);
688 const size_t segmentSize = mOffsetSegments[n]->size();
689 PointIndexType* offsets = mOffsetSegments[n]->data();
693 for (
size_t i = 0; i < segmentSize; ++i) {
694 bucketIndices[i] = bucketCounters[offsets[i]]++;
697 PointIndexType nonemptyBucketCount = 0;
698 for (
size_t i = 0; i < bucketCountersSize; ++i) {
699 nonemptyBucketCount +=
static_cast<PointIndexType
>(bucketCounters[i] != 0);
703 IndexArray& pageOffsets = mPageOffsetArrays[n];
704 pageOffsets.reset(
new PointIndexType[nonemptyBucketCount + 1]);
705 pageOffsets[0] = nonemptyBucketCount + 1;
707 IndexArray& pageIndices = mPageIndexArrays[n];
708 pageIndices.reset(
new PointIndexType[nonemptyBucketCount]);
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];
722 PointIndexType* indices = mIndexSegments[n]->data();
723 const tbb::blocked_range<size_t> segmentRange(0, segmentSize);
728 bucketIndices.get(), bucketCounters.get(), offsets));
731 offsets, bucketIndices.get(), indices));
733 mIndexSegments[n]->clear();
749 template<
typename Po
intIndexType,
typename VoxelOffsetType,
typename Po
intArray>
755 std::vector<Coord>& coords,
756 const Index binLog2Dim,
757 const Index bucketLog2Dim,
758 VoxelOffsetType* voxelOffsets =
nullptr,
759 bool cellCenteredTransform =
true)
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>;
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;
771 std::unique_ptr<IndexPairListMapPtr[]> bins(
new IndexPairListMapPtr[numTasks]);
775 tbb::parallel_for(tbb::blocked_range<size_t>(0, numTasks),
776 BinOp(bins.get(), points, voxelOffsets, xform, binLog2Dim, bucketLog2Dim,
777 numTasks, cellCenteredTransform));
779 std::set<Coord> uniqueCoords;
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);
788 coords.assign(uniqueCoords.begin(), uniqueCoords.end());
789 uniqueCoords.clear();
791 size_t segmentCount = coords.size();
795 indexSegments.reset(
new SegmentPtr[segmentCount]);
796 offsetSegments.reset(
new SegmentPtr[segmentCount]);
800 tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
801 MergeOp(bins.get(), indexSegments.get(), offsetSegments.get(), &coords[0], numTasks));
805 template<
typename Po
intIndexType,
typename VoxelOffsetType,
typename Po
intArray>
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)
820 if (recordVoxelOffsets) voxelOffsets.reset(
new VoxelOffsetType[points.size()]);
821 else voxelOffsets.reset();
823 const Index binLog2Dim = 5u;
829 std::vector<Coord> segmentCoords;
831 std::unique_ptr<SegmentPtr[]> indexSegments;
832 std::unique_ptr<SegmentPtr[]> offsetSegments;
834 binAndSegment<PointIndexType, VoxelOffsetType, PointArray>(points, xform,
835 indexSegments, offsetSegments, segmentCoords, binLog2Dim, bucketLog2Dim,
836 voxelOffsets.get(), cellCenteredTransform);
838 size_t numSegments = segmentCoords.size();
840 const tbb::blocked_range<size_t> segmentRange(0, numSegments);
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]);
846 const Index binVolume = 1u << (3u * binLog2Dim);
849 (indexSegments.get(), offsetSegments.get(),
850 pageOffsetArrays.get(), pageIndexArrays.get(), binVolume));
852 indexSegments.reset();
854 std::vector<Index> segmentOffsets;
855 segmentOffsets.reserve(numSegments);
858 for (
size_t n = 0; n < numSegments; ++n) {
859 segmentOffsets.push_back(pageCount);
860 pageCount += pageOffsetArrays[n][0] - 1;
863 pageOffsets.reset(
new PointIndexType[pageCount + 1]);
865 PointIndexType count = 0;
866 for (
size_t n = 0, idx = 0; n < numSegments; ++n) {
868 PointIndexType* offsets = pageOffsetArrays[n].get();
869 size_t size = size_t(offsets[0]);
871 for (
size_t i = 1; i < size; ++i) {
872 pageOffsets[idx++] = count;
877 pageOffsets[pageCount] = count;
879 pointIndices.reset(
new PointIndexType[points.size()]);
881 std::vector<PointIndexType*> indexArray;
882 indexArray.reserve(numSegments);
884 PointIndexType* index = pointIndices.get();
885 for (
size_t n = 0; n < numSegments; ++n) {
886 indexArray.push_back(index);
887 index += offsetSegments[n]->size();
892 pageCoordinates.reset(
new Coord[pageCount]);
894 tbb::parallel_for(segmentRange,
895 [&](tbb::blocked_range<size_t>& range)
897 for (
size_t n = range.begin(); n < range.end(); n++)
899 Index segmentOffset = segmentOffsets[n];
900 PointIndexType* indices = pageIndexArrays[n].get();
902 const Coord& segmentCoord = segmentCoords[n];
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)
910 for (size_t i = r.begin(); i < r.end(); i++)
912 Index pageIndex = indices[i];
913 Coord& ijk = pageCoordinates[segmentOffset+i];
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);
920 ijk = (ijk << bucketLog2Dim) + segmentCoord;
930 tbb::parallel_for(segmentRange,
931 MoveSegmentDataOp<PointIndexType>(indexArray, offsetSegments.get()));
941 template<
typename Po
intIndexType, Index BucketLog2Dim>
943 : mPointIndices(nullptr)
944 , mVoxelOffsets(nullptr)
945 , mPageOffsets(nullptr)
946 , mPageCoordinates(nullptr)
948 , mUsingCellCenteredTransform(true)
953 template<
typename Po
intIndexType, Index BucketLog2Dim>
958 mUsingCellCenteredTransform =
true;
959 mPointIndices.reset();
960 mVoxelOffsets.reset();
961 mPageOffsets.reset();
962 mPageCoordinates.reset();
966 template<
typename Po
intIndexType, Index BucketLog2Dim>
970 const IndexType tmpLhsPageCount = mPageCount;
971 mPageCount = rhs.mPageCount;
972 rhs.mPageCount = tmpLhsPageCount;
974 mPointIndices.swap(rhs.mPointIndices);
975 mVoxelOffsets.swap(rhs.mVoxelOffsets);
976 mPageOffsets.swap(rhs.mPageOffsets);
977 mPageCoordinates.swap(rhs.mPageCoordinates);
979 bool lhsCellCenteredTransform = mUsingCellCenteredTransform;
980 mUsingCellCenteredTransform = rhs.mUsingCellCenteredTransform;
981 rhs.mUsingCellCenteredTransform = lhsCellCenteredTransform;
985 template<
typename Po
intIndexType, Index BucketLog2Dim>
989 assert(
bool(mPointIndices) &&
bool(mPageCount));
991 mPointIndices.get() + mPageOffsets[n],
992 mPointIndices.get() + mPageOffsets[n + 1]);
996 template<
typename Po
intIndexType, Index BucketLog2Dim>
997 template<
typename Po
intArray>
1003 bool recordVoxelOffsets,
1004 bool cellCenteredTransform)
1006 mUsingCellCenteredTransform = cellCenteredTransform;
1009 mPointIndices, mPageOffsets, mPageCoordinates, mPageCount, mVoxelOffsets,
1010 (voxelOrder || recordVoxelOffsets), cellCenteredTransform);
1012 const tbb::blocked_range<size_t> pageRange(0, mPageCount);
1014 if (mVoxelOffsets && voxelOrder) {
1016 IndexType, BucketLog2Dim>(mPointIndices, mPageOffsets, mVoxelOffsets));
1019 if (mVoxelOffsets && !recordVoxelOffsets) {
1020 mVoxelOffsets.reset();
1025 template<
typename Po
intIndexType, Index BucketLog2Dim>
1026 template<
typename Po
intArray>
1032 bool recordVoxelOffsets,
1033 bool cellCenteredTransform)
1036 ret->construct(points, xform, voxelOrder, recordVoxelOffsets, cellCenteredTransform);
1049 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED