45 #ifndef OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED 46 #define OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED 49 #include <openvdb/Types.h> 50 #include <openvdb/math/Transform.h> 58 #include <boost/integer.hpp> 59 #include <boost/scoped_array.hpp> 60 #include <boost/shared_ptr.hpp> 61 #include <boost/math/special_functions/fpclassify.hpp> 63 #include <tbb/blocked_range.h> 64 #include <tbb/parallel_for.h> 65 #include <tbb/task_scheduler_init.h> 100 template<
typename Po
intIndexType = u
int32_t, Index BucketLog2Dim = 3>
104 enum { LOG2DIM = BucketLog2Dim };
106 typedef boost::shared_ptr<PointPartitioner>
Ptr;
107 typedef boost::shared_ptr<const PointPartitioner>
ConstPtr;
125 template<
typename Po
intArray>
127 bool voxelOrder =
false,
bool recordVoxelOffsets =
false);
136 template<
typename Po
intArray>
138 bool voxelOrder =
false,
bool recordVoxelOffsets =
false);
142 size_t size()
const {
return mPageCount; }
145 bool empty()
const {
return mPageCount == 0; }
158 return CoordBBox::createCube(mPageCoordinates[n], (1u << BucketLog2Dim));
162 const Coord&
origin(
size_t n)
const {
return mPageCoordinates[n]; }
173 boost::scoped_array<IndexType> mPointIndices;
174 VoxelOffsetArray mVoxelOffsets;
176 boost::scoped_array<IndexType> mPageOffsets;
177 boost::scoped_array<Coord> mPageCoordinates;
178 IndexType mPageCount;
185 template<
typename Po
intIndexType, Index BucketLog2Dim>
192 : mBegin(begin), mEnd(end), mItem(begin) {}
198 size_t size()
const {
return mEnd - mBegin; }
201 IndexType&
operator*() { assert(mItem != NULL);
return *mItem; }
202 const IndexType&
operator*()
const { assert(mItem != NULL);
return *mItem; }
205 operator bool()
const {
return mItem < mEnd; }
206 bool test()
const {
return mItem < mEnd; }
212 bool next() { this->operator++();
return this->test(); }
220 IndexType *
const mBegin, *
const mEnd;
231 namespace point_partitioner_internal {
234 template<
typename Po
intIndexType>
238 const PointIndexType* bucketCounters,
const PointIndexType* bucketOffsets)
239 : mPointOrder(pointOrder)
240 , mBucketCounters(bucketCounters)
241 , mBucketOffsets(bucketOffsets)
245 void operator()(
const tbb::blocked_range<size_t>& range)
const {
246 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
247 mPointOrder[n] += mBucketCounters[mBucketOffsets[n]];
257 template<
typename Po
intIndexType>
261 const PointIndexType* pointOrder,
const PointIndexType* indices)
262 : mOrderedIndexArray(orderedIndexArray)
263 , mPointOrder(pointOrder)
268 void operator()(
const tbb::blocked_range<size_t>& range)
const {
269 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
270 mOrderedIndexArray[mPointOrder[n]] = mIndices[n];
280 template<
typename Po
intIndexType, Index BucketLog2Dim>
287 VoxelOrderOp(IndexArray& indices,
const IndexArray& pages,
const VoxelOffsetArray& offsets)
288 : mIndices(indices.get())
289 , mPages(pages.get())
290 , mVoxelOffsets(offsets.get())
294 void operator()(
const tbb::blocked_range<size_t>& range)
const {
296 PointIndexType pointCount = 0;
297 for (
size_t n(range.begin()), N(range.end()); n != N; ++n) {
298 pointCount =
std::max(pointCount, (mPages[n + 1] - mPages[n]));
301 const PointIndexType voxelCount = 1 << (3 * BucketLog2Dim);
304 boost::scoped_array<VoxelOffsetType> offsets(
new VoxelOffsetType[pointCount]);
305 boost::scoped_array<PointIndexType> sortedIndices(
new PointIndexType[pointCount]);
306 boost::scoped_array<PointIndexType>
histogram(
new PointIndexType[voxelCount]);
308 for (
size_t n(range.begin()), N(range.end()); n != N; ++n) {
310 PointIndexType *
const indices = mIndices + mPages[n];
311 pointCount = mPages[n + 1] - mPages[n];
314 for (PointIndexType i = 0; i < pointCount; ++i) {
315 offsets[i] = mVoxelOffsets[ indices[i] ];
319 memset(&histogram[0], 0, voxelCount *
sizeof(PointIndexType));
322 for (PointIndexType i = 0; i < pointCount; ++i) {
323 ++histogram[ offsets[i] ];
326 PointIndexType count = 0, startOffset;
327 for (
int i = 0; i < int(voxelCount); ++i) {
328 if (histogram[i] > 0) {
330 count += histogram[i];
331 histogram[i] = startOffset;
336 for (PointIndexType i = 0; i < pointCount; ++i) {
337 sortedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
340 memcpy(&indices[0], &sortedIndices[0],
sizeof(PointIndexType) * pointCount);
350 template<
typename Po
intArray,
typename Po
intIndexType>
357 const IndexArray& indices,
const IndexArray& pages,
359 : mCoordinates(coordinates.get())
360 , mIndices(indices.get())
361 , mPages(pages.get())
368 void operator()(
const tbb::blocked_range<size_t>& range)
const {
370 typedef typename PointArray::value_type PointType;
372 const int mask = ~((1 << mLog2Dim) - 1);
376 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
378 mPoints->getPos(mIndices[mPages[n]], pos);
380 if (boost::math::isfinite(pos[0]) &&
381 boost::math::isfinite(pos[1]) &&
382 boost::math::isfinite(pos[2])) {
384 ijk = mXForm.worldToIndexCellCentered(pos);
390 mCoordinates[n] = ijk;
410 typedef boost::shared_ptr<Array>
Ptr;
412 Array(
size_t size) : mSize(size), mData(new T[size]) { }
414 size_t size()
const {
return mSize; }
416 T*
data() {
return mData.get(); }
417 const T*
data()
const {
return mData.get(); }
419 void clear() { mSize = 0; mData.reset(); }
423 boost::scoped_array<T> mData;
427 template<
typename Po
intIndexType>
434 : mIndexLists(&indexLists[0]), mSegments(segments)
438 void operator()(
const tbb::blocked_range<size_t>& range)
const {
439 for (
size_t n(range.begin()), N(range.end()); n != N; ++n) {
440 PointIndexType* indices = mIndexLists[n];
441 SegmentPtr&
segment = mSegments[n];
443 tbb::parallel_for(tbb::blocked_range<size_t>(0, segment->size()),
444 CopyData(indices, segment->data()));
454 CopyData(PointIndexType* lhs,
const PointIndexType* rhs) : mLhs(lhs), mRhs(rhs) { }
456 void operator()(
const tbb::blocked_range<size_t>& range)
const {
457 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
462 PointIndexType *
const mLhs;
463 PointIndexType
const *
const mRhs;
466 PointIndexType *
const *
const mIndexLists;
467 SegmentPtr *
const mSegments;
471 template<
typename Po
intIndexType>
477 typedef std::pair<PointIndexType, PointIndexType>
IndexPair;
484 SegmentPtr* indexSegments,
485 SegmentPtr* offsetSegments,
489 , mIndexSegments(indexSegments)
490 , mOffsetSegments(offsetSegments)
492 , mNumSegments(numSegments)
496 void operator()(
const tbb::blocked_range<size_t>& range)
const {
498 std::vector<IndexPairListPtr*> data;
499 std::vector<PointIndexType> arrayOffsets;
501 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
503 const Coord& ijk = mCoords[n];
504 size_t numIndices = 0;
508 for (
size_t i = 0, I = mNumSegments; i < I; ++i) {
510 IndexPairListMap& idxMap = *mBins[i];
511 typename IndexPairListMap::iterator iter = idxMap.find(ijk);
513 if (iter != idxMap.end() && iter->second) {
514 IndexPairListPtr& idxListPtr = iter->second;
516 data.push_back(&idxListPtr);
517 numIndices += idxListPtr->size();
521 if (data.empty() || numIndices == 0)
continue;
523 SegmentPtr& indexSegment = mIndexSegments[n];
524 SegmentPtr& offsetSegment = mOffsetSegments[n];
526 indexSegment.reset(
new Segment(numIndices));
527 offsetSegment.reset(
new Segment(numIndices));
529 arrayOffsets.clear();
530 arrayOffsets.reserve(data.size());
532 for (
size_t i = 0, count = 0, I = data.size(); i < I; ++i) {
533 arrayOffsets.push_back(PointIndexType(count));
534 count += (*data[i])->size();
537 tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size()),
538 CopyData(&data[0], &arrayOffsets[0], indexSegment->data(), offsetSegment->data()));
546 CopyData(IndexPairListPtr** indexLists,
547 const PointIndexType* arrayOffsets,
548 PointIndexType* indices,
549 PointIndexType* offsets)
550 : mIndexLists(indexLists)
551 , mArrayOffsets(arrayOffsets)
557 void operator()(
const tbb::blocked_range<size_t>& range)
const {
559 typedef typename IndexPairList::const_iterator CIter;
561 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
563 const PointIndexType arrayOffset = mArrayOffsets[n];
564 PointIndexType* indexPtr = &mIndices[arrayOffset];
565 PointIndexType* offsetPtr = &mOffsets[arrayOffset];
567 IndexPairListPtr& list = *mIndexLists[n];
569 for (CIter it = list->begin(), end = list->end(); it != end; ++it) {
570 const IndexPair& data = *it;
571 *indexPtr++ = data.first;
572 *offsetPtr++ = data.second;
579 IndexPairListPtr *
const *
const mIndexLists;
580 PointIndexType
const *
const mArrayOffsets;
581 PointIndexType *
const mIndices;
582 PointIndexType *
const mOffsets;
585 IndexPairListMapPtr *
const mBins;
586 SegmentPtr *
const mIndexSegments;
587 SegmentPtr *
const mOffsetSegments;
588 Coord
const *
const mCoords;
589 size_t const mNumSegments;
593 template<
typename Po
intArray,
typename Po
intIndexType,
typename VoxelOffsetType>
599 typedef std::pair<PointIndexType, PointIndexType>
IndexPair;
607 VoxelOffsetType* voxelOffsets,
614 , mVoxelOffsets(voxelOffsets)
616 , mBinLog2Dim(binLog2Dim)
617 , mBucketLog2Dim(bucketLog2Dim)
618 , mNumSegments(numSegments)
622 void operator()(
const tbb::blocked_range<size_t>& range)
const {
624 const Index log2dim = mBucketLog2Dim;
625 const Index log2dim2 = 2 * log2dim;
626 const Index bucketMask = (1u << log2dim) - 1u;
628 const Index binLog2dim = mBinLog2Dim;
629 const Index binLog2dim2 = 2 * binLog2dim;
631 const Index binMask = (1u << (log2dim + binLog2dim)) - 1u;
632 const Index invBinMask = ~binMask;
634 IndexPairList * idxList = NULL;
635 Coord ijk(0, 0, 0), loc(0, 0, 0), binCoord(0, 0, 0), lastBinCoord(1, 2, 3);
638 PointIndexType bucketOffset = 0;
639 VoxelOffsetType voxelOffset = 0;
641 const size_t numPoints = mPoints->size();
642 const size_t segmentSize = numPoints / mNumSegments;
644 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
646 IndexPairListMapPtr& dataPtr = mData[n];
647 if (!dataPtr) dataPtr.reset(
new IndexPairListMap());
648 IndexPairListMap& idxMap = *dataPtr;
650 const bool isLastSegment = (n + 1) >= mNumSegments;
652 const size_t start = n * segmentSize;
653 const size_t end = isLastSegment ? numPoints : (start + segmentSize);
655 for (
size_t i = start; i != end; ++i) {
657 mPoints->getPos(i, pos);
659 if (boost::math::isfinite(pos[0]) &&
660 boost::math::isfinite(pos[1]) &&
661 boost::math::isfinite(pos[2])) {
663 ijk = mXForm.worldToIndexCellCentered(pos);
666 loc[0] = ijk[0] & bucketMask;
667 loc[1] = ijk[1] & bucketMask;
668 loc[2] = ijk[2] & bucketMask;
669 voxelOffset = VoxelOffsetType((loc[0] << log2dim2) + (loc[1] << log2dim) + loc[2]);
672 binCoord[0] = ijk[0] & invBinMask;
673 binCoord[1] = ijk[1] & invBinMask;
674 binCoord[2] = ijk[2] & invBinMask;
684 bucketOffset = PointIndexType((ijk[0] << binLog2dim2) + (ijk[1] << binLog2dim) + ijk[2]);
686 if (lastBinCoord != binCoord) {
687 lastBinCoord = binCoord;
688 IndexPairListPtr& idxListPtr = idxMap[lastBinCoord];
689 if (!idxListPtr) idxListPtr.reset(
new IndexPairList());
690 idxList = idxListPtr.get();
693 idxList->push_back(IndexPair(PointIndexType(i), bucketOffset));
694 if (mVoxelOffsets) mVoxelOffsets[i] = voxelOffset;
710 template<
typename Po
intIndexType>
717 IndexArray* pageOffsetArrays,
Index binVolume)
718 : mIndexSegments(indexSegments)
719 , mOffsetSegments(offestSegments)
720 , mPageOffsetArrays(pageOffsetArrays)
721 , mBinVolume(binVolume)
725 void operator()(
const tbb::blocked_range<size_t>& range)
const {
727 const size_t bucketCountersSize = size_t(mBinVolume);
728 IndexArray bucketCounters(
new PointIndexType[bucketCountersSize]);
730 size_t maxSegmentSize = 0;
731 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
732 maxSegmentSize =
std::max(maxSegmentSize, mIndexSegments[n]->size());
735 IndexArray bucketIndices(
new PointIndexType[maxSegmentSize]);
738 for (
size_t n = range.begin(), N = range.end(); n != N; ++n) {
740 memset(bucketCounters.get(), 0,
sizeof(PointIndexType) * bucketCountersSize);
742 const size_t segmentSize = mOffsetSegments[n]->size();
743 PointIndexType* offsets = mOffsetSegments[n]->data();
747 for (
size_t i = 0; i < segmentSize; ++i) {
748 bucketIndices[i] = bucketCounters[offsets[i]]++;
751 PointIndexType nonemptyBucketCount = 0;
752 for (
size_t i = 0; i < bucketCountersSize; ++i) {
753 nonemptyBucketCount +=
static_cast<PointIndexType
>(bucketCounters[i] != 0);
757 IndexArray& pageOffsets = mPageOffsetArrays[n];
758 pageOffsets.reset(
new PointIndexType[nonemptyBucketCount + 1]);
759 pageOffsets[0] = nonemptyBucketCount + 1;
762 PointIndexType count = 0, idx = 1;
763 for (
size_t i = 0; i < bucketCountersSize; ++i) {
764 if (bucketCounters[i] != 0) {
765 pageOffsets[idx] = bucketCounters[i];
766 bucketCounters[i] = count;
767 count += pageOffsets[idx];
772 PointIndexType* indices = mIndexSegments[n]->data();
773 const tbb::blocked_range<size_t> segmentRange(0, segmentSize);
778 bucketIndices.get(), bucketCounters.get(), offsets));
781 offsets, bucketIndices.get(), indices));
783 mIndexSegments[n]->clear();
798 template<
typename Po
intIndexType,
typename VoxelOffsetType,
typename Po
intArray>
804 size_t& segmentCount,
805 const Index binLog2Dim,
806 const Index bucketLog2Dim,
807 VoxelOffsetType* voxelOffsets = NULL)
809 typedef std::pair<PointIndexType, PointIndexType> IndexPair;
810 typedef std::deque<IndexPair> IndexPairList;
811 typedef boost::shared_ptr<IndexPairList> IndexPairListPtr;
812 typedef std::map<Coord, IndexPairListPtr> IndexPairListMap;
813 typedef boost::shared_ptr<IndexPairListMap> IndexPairListMapPtr;
815 size_t numTasks = 1, numThreads = size_t(tbb::task_scheduler_init::default_num_threads());
816 if (points.size() > (numThreads * 2)) numTasks = numThreads * 2;
817 else if (points.size() > numThreads) numTasks = numThreads;
819 boost::scoped_array<IndexPairListMapPtr> bins(
new IndexPairListMapPtr[numTasks]);
823 tbb::parallel_for(tbb::blocked_range<size_t>(0, numTasks),
824 BinOp(bins.get(), points, voxelOffsets, xform, binLog2Dim, bucketLog2Dim, numTasks));
826 std::set<Coord> uniqueCoords;
828 for (
size_t i = 0; i < numTasks; ++i) {
829 IndexPairListMap& idxMap = *bins[i];
830 for (
typename IndexPairListMap::iterator it = idxMap.begin(); it != idxMap.end(); ++it) {
831 uniqueCoords.insert(it->first);
835 std::vector<Coord> coords(uniqueCoords.begin(), uniqueCoords.end());
836 uniqueCoords.clear();
838 segmentCount = coords.size();
842 indexSegments.reset(
new SegmentPtr[segmentCount]);
843 offsetSegments.reset(
new SegmentPtr[segmentCount]);
847 tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
848 MergeOp(bins.get(), indexSegments.get(), offsetSegments.get(), &coords[0], numTasks));
852 template<
typename Po
intIndexType,
typename VoxelOffsetType,
typename Po
intArray>
856 const Index bucketLog2Dim,
857 boost::scoped_array<PointIndexType>& pointIndices,
858 boost::scoped_array<PointIndexType>& pageOffsets,
859 PointIndexType& pageCount,
860 boost::scoped_array<VoxelOffsetType>& voxelOffsets,
861 bool recordVoxelOffsets)
863 if (recordVoxelOffsets) voxelOffsets.reset(
new VoxelOffsetType[points.size()]);
864 else voxelOffsets.reset();
866 const Index binLog2Dim = 5u;
872 size_t numSegments = 0;
874 boost::scoped_array<typename Array<PointIndexType>::Ptr> indexSegments;
875 boost::scoped_array<typename Array<PointIndexType>::Ptr> offestSegments;
877 binAndSegment<PointIndexType, VoxelOffsetType, PointArray>(points, xform,
878 indexSegments, offestSegments, numSegments, binLog2Dim, bucketLog2Dim, voxelOffsets.get());
880 const tbb::blocked_range<size_t> segmentRange(0, numSegments);
882 typedef boost::scoped_array<PointIndexType> IndexArray;
883 boost::scoped_array<IndexArray> pageOffsetArrays(
new IndexArray[numSegments]);
885 const Index binVolume = 1u << (3u * binLog2Dim);
888 (indexSegments.get(), offestSegments.get(), pageOffsetArrays.get(), binVolume));
890 indexSegments.reset();
893 for (
size_t n = 0; n < numSegments; ++n) {
894 pageCount += pageOffsetArrays[n][0] - 1;
897 pageOffsets.reset(
new PointIndexType[pageCount + 1]);
899 PointIndexType count = 0;
900 for (
size_t n = 0, idx = 0; n < numSegments; ++n) {
902 PointIndexType* offsets = pageOffsetArrays[n].get();
903 size_t size = size_t(offsets[0]);
905 for (
size_t i = 1; i < size; ++i) {
906 pageOffsets[idx++] = count;
911 pageOffsets[pageCount] = count;
913 pointIndices.reset(
new PointIndexType[points.size()]);
915 std::vector<PointIndexType*> indexArray;
916 indexArray.reserve(numSegments);
918 PointIndexType* index = pointIndices.get();
919 for (
size_t n = 0; n < numSegments; ++n) {
920 indexArray.push_back(index);
921 index += offestSegments[n]->size();
934 template<
typename Po
intIndexType, Index BucketLog2Dim>
936 : mPointIndices(NULL)
937 , mVoxelOffsets(NULL)
939 , mPageCoordinates(NULL)
945 template<
typename Po
intIndexType, Index BucketLog2Dim>
950 mPointIndices.reset();
951 mVoxelOffsets.reset();
952 mPageOffsets.reset();
953 mPageCoordinates.reset();
957 template<
typename Po
intIndexType, Index BucketLog2Dim>
961 const IndexType tmpLhsPageCount = mPageCount;
962 mPageCount = rhs.mPageCount;
963 rhs.mPageCount = tmpLhsPageCount;
965 mPointIndices.swap(rhs.mPointIndices);
966 mVoxelOffsets.swap(rhs.mVoxelOffsets);
967 mPageOffsets.swap(rhs.mPageOffsets);
968 mPageCoordinates.swap(rhs.mPageCoordinates);
972 template<
typename Po
intIndexType, Index BucketLog2Dim>
976 assert(
bool(mPointIndices) &&
bool(mPageCount));
978 mPointIndices.get() + mPageOffsets[n],
979 mPointIndices.get() + mPageOffsets[n + 1]);
983 template<
typename Po
intIndexType, Index BucketLog2Dim>
984 template<
typename Po
intArray>
987 const math::Transform& xform,
bool voxelOrder,
bool recordVoxelOffsets)
990 mPointIndices, mPageOffsets, mPageCount, mVoxelOffsets, (voxelOrder || recordVoxelOffsets));
992 const tbb::blocked_range<size_t> pageRange(0, mPageCount);
993 mPageCoordinates.reset(
new Coord[mPageCount]);
995 tbb::parallel_for(pageRange,
997 (mPageCoordinates, mPointIndices, mPageOffsets, points, xform, BucketLog2Dim));
999 if (mVoxelOffsets && voxelOrder) {
1001 IndexType, BucketLog2Dim>(mPointIndices, mPageOffsets, mVoxelOffsets));
1004 if (mVoxelOffsets && !recordVoxelOffsets) {
1005 mVoxelOffsets.reset();
1010 template<
typename Po
intIndexType, Index BucketLog2Dim>
1011 template<
typename Po
intArray>
1014 bool voxelOrder,
bool recordVoxelOffsets)
1017 ret->construct(points, xform, voxelOrder, recordVoxelOffsets);
1030 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
Index32 Index
Definition: Types.h:58
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:450
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Definition: Exceptions.h:39
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71