OpenVDB  3.2.0
PointPartitioner.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2016 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
44 
45 #ifndef OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
46 #define OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
47 
48 
49 #include <openvdb/Types.h>
50 #include <openvdb/math/Transform.h>
51 
52 #include <deque>
53 #include <map>
54 #include <set>
55 #include <utility> // std::pair
56 #include <vector>
57 
58 #include <boost/integer.hpp> // boost::int_t<N>::least
59 #include <boost/scoped_array.hpp>
60 #include <boost/shared_ptr.hpp>
61 #include <boost/math/special_functions/fpclassify.hpp> // boost::math::isfinite
62 
63 #include <tbb/blocked_range.h>
64 #include <tbb/parallel_for.h>
65 #include <tbb/task_scheduler_init.h>
66 
67 
68 namespace openvdb {
70 namespace OPENVDB_VERSION_NAME {
71 namespace tools {
72 
73 
75 
76 
106 template<typename PointIndexType = uint32_t, Index BucketLog2Dim = 3>
108 {
109 public:
110  enum { LOG2DIM = BucketLog2Dim };
111 
112  typedef boost::shared_ptr<PointPartitioner> Ptr;
113  typedef boost::shared_ptr<const PointPartitioner> ConstPtr;
114 
115  typedef PointIndexType IndexType;
116  typedef typename boost::int_t<1 + (3 * BucketLog2Dim)>::least VoxelOffsetType;
117  typedef boost::scoped_array<VoxelOffsetType> VoxelOffsetArray;
118 
119  class IndexIterator;
120 
122 
124 
134  template<typename PointArray>
135  void construct(const PointArray& points, const math::Transform& xform,
136  bool voxelOrder = false, bool recordVoxelOffsets = false, bool cellCenteredTransform = true);
137 
138 
148  template<typename PointArray>
149  static Ptr create(const PointArray& points, const math::Transform& xform,
150  bool voxelOrder = false, bool recordVoxelOffsets = false, bool cellCenteredTransform = true);
151 
152 
154  size_t size() const { return mPageCount; }
155 
157  bool empty() const { return mPageCount == 0; }
158 
160  void clear();
161 
163  void swap(PointPartitioner&);
164 
166  IndexIterator indices(size_t n) const;
167 
169  CoordBBox getBBox(size_t n) const {
170  return CoordBBox::createCube(mPageCoordinates[n], (1u << BucketLog2Dim));
171  }
172 
174  const Coord& origin(size_t n) const { return mPageCoordinates[n]; }
175 
178  const VoxelOffsetArray& voxelOffsets() const { return mVoxelOffsets; }
179 
183  bool usingCellCenteredTransform() const { return mUsingCellCenteredTransform; }
184 
185 private:
186  // Disallow copying
188  PointPartitioner& operator=(const PointPartitioner&);
189 
190  boost::scoped_array<IndexType> mPointIndices;
191  VoxelOffsetArray mVoxelOffsets;
192 
193  boost::scoped_array<IndexType> mPageOffsets;
194  boost::scoped_array<Coord> mPageCoordinates;
195  IndexType mPageCount;
196  bool mUsingCellCenteredTransform;
197 }; // class PointPartitioner
198 
199 
201 
202 
203 template<typename PointIndexType, Index BucketLog2Dim>
204 class PointPartitioner<PointIndexType, BucketLog2Dim>::IndexIterator
205 {
206 public:
207  typedef PointIndexType IndexType;
208 
209  IndexIterator(IndexType* begin = NULL, IndexType* end = NULL)
210  : mBegin(begin), mEnd(end), mItem(begin) {}
211 
213  void reset() { mItem = mBegin; }
214 
216  size_t size() const { return mEnd - mBegin; }
217 
219  IndexType& operator*() { assert(mItem != NULL); return *mItem; }
220  const IndexType& operator*() const { assert(mItem != NULL); return *mItem; }
221 
223  operator bool() const { return mItem < mEnd; }
224  bool test() const { return mItem < mEnd; }
225 
227  IndexIterator& operator++() { assert(this->test()); ++mItem; return *this; }
228 
230  bool next() { this->operator++(); return this->test(); }
231  bool increment() { this->next(); return this->test(); }
232 
234  bool operator==(const IndexIterator& other) const { return mItem == other.mItem; }
235  bool operator!=(const IndexIterator& other) const { return !this->operator==(other); }
236 
237 private:
238  IndexType * const mBegin, * const mEnd;
239  IndexType * mItem;
240 }; // class PointPartitioner::IndexIterator
241 
242 
245 
246 // Implementation details
247 
248 
249 namespace point_partitioner_internal {
250 
251 
252 template<typename PointIndexType>
254 {
255  ComputePointOrderOp(PointIndexType* pointOrder,
256  const PointIndexType* bucketCounters, const PointIndexType* bucketOffsets)
257  : mPointOrder(pointOrder)
258  , mBucketCounters(bucketCounters)
259  , mBucketOffsets(bucketOffsets)
260  {
261  }
262 
263  void operator()(const tbb::blocked_range<size_t>& range) const {
264  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
265  mPointOrder[n] += mBucketCounters[mBucketOffsets[n]];
266  }
267  }
268 
269  PointIndexType * const mPointOrder;
270  PointIndexType const * const mBucketCounters;
271  PointIndexType const * const mBucketOffsets;
272 }; // struct ComputePointOrderOp
273 
274 
275 template<typename PointIndexType>
277 {
278  CreateOrderedPointIndexArrayOp(PointIndexType* orderedIndexArray,
279  const PointIndexType* pointOrder, const PointIndexType* indices)
280  : mOrderedIndexArray(orderedIndexArray)
281  , mPointOrder(pointOrder)
282  , mIndices(indices)
283  {
284  }
285 
286  void operator()(const tbb::blocked_range<size_t>& range) const {
287  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
288  mOrderedIndexArray[mPointOrder[n]] = mIndices[n];
289  }
290  }
291 
292  PointIndexType * const mOrderedIndexArray;
293  PointIndexType const * const mPointOrder;
294  PointIndexType const * const mIndices;
295 }; // struct CreateOrderedPointIndexArrayOp
296 
297 
298 template<typename PointIndexType, Index BucketLog2Dim>
300 {
301  typedef typename boost::int_t<1 + (3 * BucketLog2Dim)>::least VoxelOffsetType;
302  typedef boost::scoped_array<VoxelOffsetType> VoxelOffsetArray;
303  typedef boost::scoped_array<PointIndexType> IndexArray;
304 
305  VoxelOrderOp(IndexArray& indices, const IndexArray& pages,const VoxelOffsetArray& offsets)
306  : mIndices(indices.get())
307  , mPages(pages.get())
308  , mVoxelOffsets(offsets.get())
309  {
310  }
311 
312  void operator()(const tbb::blocked_range<size_t>& range) const {
313 
314  PointIndexType pointCount = 0;
315  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
316  pointCount = std::max(pointCount, (mPages[n + 1] - mPages[n]));
317  }
318 
319  const PointIndexType voxelCount = 1 << (3 * BucketLog2Dim);
320 
321  // allocate histogram buffers
322  boost::scoped_array<VoxelOffsetType> offsets(new VoxelOffsetType[pointCount]);
323  boost::scoped_array<PointIndexType> sortedIndices(new PointIndexType[pointCount]);
324  boost::scoped_array<PointIndexType> histogram(new PointIndexType[voxelCount]);
325 
326  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
327 
328  PointIndexType * const indices = mIndices + mPages[n];
329  pointCount = mPages[n + 1] - mPages[n];
330 
331  // local copy of voxel offsets.
332  for (PointIndexType i = 0; i < pointCount; ++i) {
333  offsets[i] = mVoxelOffsets[ indices[i] ];
334  }
335 
336  // reset histogram
337  memset(&histogram[0], 0, voxelCount * sizeof(PointIndexType));
338 
339  // compute histogram
340  for (PointIndexType i = 0; i < pointCount; ++i) {
341  ++histogram[ offsets[i] ];
342  }
343 
344  PointIndexType count = 0, startOffset;
345  for (int i = 0; i < int(voxelCount); ++i) {
346  if (histogram[i] > 0) {
347  startOffset = count;
348  count += histogram[i];
349  histogram[i] = startOffset;
350  }
351  }
352 
353  // sort indices based on voxel offset
354  for (PointIndexType i = 0; i < pointCount; ++i) {
355  sortedIndices[ histogram[ offsets[i] ]++ ] = indices[i];
356  }
357 
358  memcpy(&indices[0], &sortedIndices[0], sizeof(PointIndexType) * pointCount);
359  }
360  }
361 
362  PointIndexType * const mIndices;
363  PointIndexType const * const mPages;
364  VoxelOffsetType const * const mVoxelOffsets;
365 }; // struct VoxelOrderOp
366 
367 
368 template<typename PointArray, typename PointIndexType>
370 {
371  typedef boost::scoped_array<PointIndexType> IndexArray;
372  typedef boost::scoped_array<Coord> CoordArray;
373 
374  LeafNodeOriginOp(CoordArray& coordinates,
375  const IndexArray& indices, const IndexArray& pages,
376  const PointArray& points, const math::Transform& m, int log2dim, bool cellCenteredTransform)
377  : mCoordinates(coordinates.get())
378  , mIndices(indices.get())
379  , mPages(pages.get())
380  , mPoints(&points)
381  , mXForm(m)
382  , mLog2Dim(log2dim)
383  , mCellCenteredTransform(cellCenteredTransform)
384  {
385  }
386 
387  void operator()(const tbb::blocked_range<size_t>& range) const {
388 
389  typedef typename PointArray::PosType PosType;
390 
391  const bool cellCentered = mCellCenteredTransform;
392  const int mask = ~((1 << mLog2Dim) - 1);
393  Coord ijk;
394  PosType pos;
395 
396  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
397 
398  mPoints->getPos(mIndices[mPages[n]], pos);
399 
400  if (boost::math::isfinite(pos[0]) &&
401  boost::math::isfinite(pos[1]) &&
402  boost::math::isfinite(pos[2])) {
403 
404  ijk = cellCentered ? mXForm.worldToIndexCellCentered(pos) :
405  mXForm.worldToIndexNodeCentered(pos);
406 
407  ijk[0] &= mask;
408  ijk[1] &= mask;
409  ijk[2] &= mask;
410 
411  mCoordinates[n] = ijk;
412  }
413  }
414  }
415 
416  Coord * const mCoordinates;
417  PointIndexType const * const mIndices;
418  PointIndexType const * const mPages;
419  PointArray const * const mPoints;
421  int const mLog2Dim;
423 }; // struct LeafNodeOriginOp
424 
425 
427 
428 
429 template<typename T>
430 struct Array
431 {
432  typedef boost::shared_ptr<Array> Ptr;
433 
434  Array(size_t size) : mSize(size), mData(new T[size]) { }
435 
436  size_t size() const { return mSize; }
437 
438  T* data() { return mData.get(); }
439  const T* data() const { return mData.get(); }
440 
441  void clear() { mSize = 0; mData.reset(); }
442 
443 private:
444  size_t mSize;
445  boost::scoped_array<T> mData;
446 }; // struct Array
447 
448 
449 template<typename PointIndexType>
451 {
453  typedef typename Segment::Ptr SegmentPtr;
454 
455  MoveSegmentDataOp(std::vector<PointIndexType*>& indexLists, SegmentPtr* segments)
456  : mIndexLists(&indexLists[0]), mSegments(segments)
457  {
458  }
459 
460  void operator()(const tbb::blocked_range<size_t>& range) const {
461  for (size_t n(range.begin()), N(range.end()); n != N; ++n) {
462  PointIndexType* indices = mIndexLists[n];
463  SegmentPtr& segment = mSegments[n];
464 
465  tbb::parallel_for(tbb::blocked_range<size_t>(0, segment->size()),
466  CopyData(indices, segment->data()));
467 
468  segment.reset(); // clear data
469  }
470  }
471 
472 private:
473 
474  struct CopyData
475  {
476  CopyData(PointIndexType* lhs, const PointIndexType* rhs) : mLhs(lhs), mRhs(rhs) { }
477 
478  void operator()(const tbb::blocked_range<size_t>& range) const {
479  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
480  mLhs[n] = mRhs[n];
481  }
482  }
483 
484  PointIndexType * const mLhs;
485  PointIndexType const * const mRhs;
486  };
487 
488  PointIndexType * const * const mIndexLists;
489  SegmentPtr * const mSegments;
490 }; // struct MoveSegmentDataOp
491 
492 
493 template<typename PointIndexType>
495 {
497  typedef typename Segment::Ptr SegmentPtr;
498 
499  typedef std::pair<PointIndexType, PointIndexType> IndexPair;
500  typedef std::deque<IndexPair> IndexPairList;
501  typedef boost::shared_ptr<IndexPairList> IndexPairListPtr;
502  typedef std::map<Coord, IndexPairListPtr> IndexPairListMap;
503  typedef boost::shared_ptr<IndexPairListMap> IndexPairListMapPtr;
504 
505  MergeBinsOp(IndexPairListMapPtr* bins,
506  SegmentPtr* indexSegments,
507  SegmentPtr* offsetSegments,
508  Coord* coords,
509  size_t numSegments)
510  : mBins(bins)
511  , mIndexSegments(indexSegments)
512  , mOffsetSegments(offsetSegments)
513  , mCoords(coords)
514  , mNumSegments(numSegments)
515  {
516  }
517 
518  void operator()(const tbb::blocked_range<size_t>& range) const {
519 
520  std::vector<IndexPairListPtr*> data;
521  std::vector<PointIndexType> arrayOffsets;
522 
523  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
524 
525  const Coord& ijk = mCoords[n];
526  size_t numIndices = 0;
527 
528  data.clear();
529 
530  for (size_t i = 0, I = mNumSegments; i < I; ++i) {
531 
532  IndexPairListMap& idxMap = *mBins[i];
533  typename IndexPairListMap::iterator iter = idxMap.find(ijk);
534 
535  if (iter != idxMap.end() && iter->second) {
536  IndexPairListPtr& idxListPtr = iter->second;
537 
538  data.push_back(&idxListPtr);
539  numIndices += idxListPtr->size();
540  }
541  }
542 
543  if (data.empty() || numIndices == 0) continue;
544 
545  SegmentPtr& indexSegment = mIndexSegments[n];
546  SegmentPtr& offsetSegment = mOffsetSegments[n];
547 
548  indexSegment.reset(new Segment(numIndices));
549  offsetSegment.reset(new Segment(numIndices));
550 
551  arrayOffsets.clear();
552  arrayOffsets.reserve(data.size());
553 
554  for (size_t i = 0, count = 0, I = data.size(); i < I; ++i) {
555  arrayOffsets.push_back(PointIndexType(count));
556  count += (*data[i])->size();
557  }
558 
559  tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size()),
560  CopyData(&data[0], &arrayOffsets[0], indexSegment->data(), offsetSegment->data()));
561  }
562  }
563 
564 private:
565 
566  struct CopyData
567  {
568  CopyData(IndexPairListPtr** indexLists,
569  const PointIndexType* arrayOffsets,
570  PointIndexType* indices,
571  PointIndexType* offsets)
572  : mIndexLists(indexLists)
573  , mArrayOffsets(arrayOffsets)
574  , mIndices(indices)
575  , mOffsets(offsets)
576  {
577  }
578 
579  void operator()(const tbb::blocked_range<size_t>& range) const {
580 
581  typedef typename IndexPairList::const_iterator CIter;
582 
583  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
584 
585  const PointIndexType arrayOffset = mArrayOffsets[n];
586  PointIndexType* indexPtr = &mIndices[arrayOffset];
587  PointIndexType* offsetPtr = &mOffsets[arrayOffset];
588 
589  IndexPairListPtr& list = *mIndexLists[n];
590 
591  for (CIter it = list->begin(), end = list->end(); it != end; ++it) {
592  const IndexPair& data = *it;
593  *indexPtr++ = data.first;
594  *offsetPtr++ = data.second;
595  }
596 
597  list.reset(); // clear data
598  }
599  }
600 
601  IndexPairListPtr * const * const mIndexLists;
602  PointIndexType const * const mArrayOffsets;
603  PointIndexType * const mIndices;
604  PointIndexType * const mOffsets;
605  }; // struct CopyData
606 
607  IndexPairListMapPtr * const mBins;
608  SegmentPtr * const mIndexSegments;
609  SegmentPtr * const mOffsetSegments;
610  Coord const * const mCoords;
611  size_t const mNumSegments;
612 }; // struct MergeBinsOp
613 
614 
615 template<typename PointArray, typename PointIndexType, typename VoxelOffsetType>
617 {
618  typedef typename PointArray::PosType PosType;
619  typedef std::pair<PointIndexType, PointIndexType> IndexPair;
620  typedef std::deque<IndexPair> IndexPairList;
621  typedef boost::shared_ptr<IndexPairList> IndexPairListPtr;
622  typedef std::map<Coord, IndexPairListPtr> IndexPairListMap;
623  typedef boost::shared_ptr<IndexPairListMap> IndexPairListMapPtr;
624 
625  BinPointIndicesOp(IndexPairListMapPtr* data,
626  const PointArray& points,
627  VoxelOffsetType* voxelOffsets,
628  const math::Transform& m,
629  Index binLog2Dim,
630  Index bucketLog2Dim,
631  size_t numSegments,
632  bool cellCenteredTransform)
633  : mData(data)
634  , mPoints(&points)
635  , mVoxelOffsets(voxelOffsets)
636  , mXForm(m)
637  , mBinLog2Dim(binLog2Dim)
638  , mBucketLog2Dim(bucketLog2Dim)
639  , mNumSegments(numSegments)
640  , mCellCenteredTransform(cellCenteredTransform)
641  {
642  }
643 
644  void operator()(const tbb::blocked_range<size_t>& range) const {
645 
646  const Index log2dim = mBucketLog2Dim;
647  const Index log2dim2 = 2 * log2dim;
648  const Index bucketMask = (1u << log2dim) - 1u;
649 
650  const Index binLog2dim = mBinLog2Dim;
651  const Index binLog2dim2 = 2 * binLog2dim;
652 
653  const Index binMask = (1u << (log2dim + binLog2dim)) - 1u;
654  const Index invBinMask = ~binMask;
655 
656  IndexPairList * idxList = NULL;
657  Coord ijk(0, 0, 0), loc(0, 0, 0), binCoord(0, 0, 0), lastBinCoord(1, 2, 3);
658  PosType pos;
659 
660  PointIndexType bucketOffset = 0;
661  VoxelOffsetType voxelOffset = 0;
662 
663  const bool cellCentered = mCellCenteredTransform;
664 
665  const size_t numPoints = mPoints->size();
666  const size_t segmentSize = numPoints / mNumSegments;
667 
668  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
669 
670  IndexPairListMapPtr& dataPtr = mData[n];
671  if (!dataPtr) dataPtr.reset(new IndexPairListMap());
672  IndexPairListMap& idxMap = *dataPtr;
673 
674  const bool isLastSegment = (n + 1) >= mNumSegments;
675 
676  const size_t start = n * segmentSize;
677  const size_t end = isLastSegment ? numPoints : (start + segmentSize);
678 
679  for (size_t i = start; i != end; ++i) {
680 
681  mPoints->getPos(i, pos);
682 
683  if (boost::math::isfinite(pos[0]) &&
684  boost::math::isfinite(pos[1]) &&
685  boost::math::isfinite(pos[2])) {
686 
687  ijk = cellCentered ? mXForm.worldToIndexCellCentered(pos) :
688  mXForm.worldToIndexNodeCentered(pos);
689 
690  if (mVoxelOffsets) {
691  loc[0] = ijk[0] & bucketMask;
692  loc[1] = ijk[1] & bucketMask;
693  loc[2] = ijk[2] & bucketMask;
694  voxelOffset = VoxelOffsetType((loc[0] << log2dim2) + (loc[1] << log2dim) + loc[2]);
695  }
696 
697  binCoord[0] = ijk[0] & invBinMask;
698  binCoord[1] = ijk[1] & invBinMask;
699  binCoord[2] = ijk[2] & invBinMask;
700 
701  ijk[0] &= binMask;
702  ijk[1] &= binMask;
703  ijk[2] &= binMask;
704 
705  ijk[0] >>= log2dim;
706  ijk[1] >>= log2dim;
707  ijk[2] >>= log2dim;
708 
709  bucketOffset = PointIndexType((ijk[0] << binLog2dim2) + (ijk[1] << binLog2dim) + ijk[2]);
710 
711  if (lastBinCoord != binCoord) {
712  lastBinCoord = binCoord;
713  IndexPairListPtr& idxListPtr = idxMap[lastBinCoord];
714  if (!idxListPtr) idxListPtr.reset(new IndexPairList());
715  idxList = idxListPtr.get();
716  }
717 
718  idxList->push_back(IndexPair(PointIndexType(i), bucketOffset));
719  if (mVoxelOffsets) mVoxelOffsets[i] = voxelOffset;
720  }
721  }
722  }
723  }
724 
725  IndexPairListMapPtr * const mData;
726  PointArray const * const mPoints;
727  VoxelOffsetType * const mVoxelOffsets;
731  size_t const mNumSegments;
733 }; // struct BinPointIndicesOp
734 
735 
736 template<typename PointIndexType>
738 {
739  typedef boost::scoped_array<PointIndexType> IndexArray;
741 
742  OrderSegmentsOp(SegmentPtr* indexSegments, SegmentPtr* offestSegments,
743  IndexArray* pageOffsetArrays, Index binVolume)
744  : mIndexSegments(indexSegments)
745  , mOffsetSegments(offestSegments)
746  , mPageOffsetArrays(pageOffsetArrays)
747  , mBinVolume(binVolume)
748  {
749  }
750 
751  void operator()(const tbb::blocked_range<size_t>& range) const {
752 
753  const size_t bucketCountersSize = size_t(mBinVolume);
754  IndexArray bucketCounters(new PointIndexType[bucketCountersSize]);
755 
756  size_t maxSegmentSize = 0;
757  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
758  maxSegmentSize = std::max(maxSegmentSize, mIndexSegments[n]->size());
759  }
760 
761  IndexArray bucketIndices(new PointIndexType[maxSegmentSize]);
762 
763 
764  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
765 
766  memset(bucketCounters.get(), 0, sizeof(PointIndexType) * bucketCountersSize);
767 
768  const size_t segmentSize = mOffsetSegments[n]->size();
769  PointIndexType* offsets = mOffsetSegments[n]->data();
770 
771  // Count the number of points per bucket and assign a local bucket index
772  // to each point.
773  for (size_t i = 0; i < segmentSize; ++i) {
774  bucketIndices[i] = bucketCounters[offsets[i]]++;
775  }
776 
777  PointIndexType nonemptyBucketCount = 0;
778  for (size_t i = 0; i < bucketCountersSize; ++i) {
779  nonemptyBucketCount += static_cast<PointIndexType>(bucketCounters[i] != 0);
780  }
781 
782 
783  IndexArray& pageOffsets = mPageOffsetArrays[n];
784  pageOffsets.reset(new PointIndexType[nonemptyBucketCount + 1]);
785  pageOffsets[0] = nonemptyBucketCount + 1; // stores array size in first element
786 
787  // Compute bucket counter prefix sum
788  PointIndexType count = 0, idx = 1;
789  for (size_t i = 0; i < bucketCountersSize; ++i) {
790  if (bucketCounters[i] != 0) {
791  pageOffsets[idx] = bucketCounters[i];
792  bucketCounters[i] = count;
793  count += pageOffsets[idx];
794  ++idx;
795  }
796  }
797 
798  PointIndexType* indices = mIndexSegments[n]->data();
799  const tbb::blocked_range<size_t> segmentRange(0, segmentSize);
800 
801  // Compute final point order by incrementing the local bucket point index
802  // with the prefix sum offset.
803  tbb::parallel_for(segmentRange, ComputePointOrderOp<PointIndexType>(
804  bucketIndices.get(), bucketCounters.get(), offsets));
805 
806  tbb::parallel_for(segmentRange, CreateOrderedPointIndexArrayOp<PointIndexType>(
807  offsets, bucketIndices.get(), indices));
808 
809  mIndexSegments[n]->clear(); // clear data
810  }
811  }
812 
813  SegmentPtr * const mIndexSegments;
814  SegmentPtr * const mOffsetSegments;
815  IndexArray * const mPageOffsetArrays;
817 }; // struct OrderSegmentsOp
818 
819 
821 
822 
824 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
825 inline void binAndSegment(
826  const PointArray& points,
827  const math::Transform& xform,
828  boost::scoped_array<typename Array<PointIndexType>::Ptr>& indexSegments,
829  boost::scoped_array<typename Array<PointIndexType>::Ptr>& offsetSegments,
830  size_t& segmentCount,
831  const Index binLog2Dim,
832  const Index bucketLog2Dim,
833  VoxelOffsetType* voxelOffsets = NULL,
834  bool cellCenteredTransform = true)
835 {
836  typedef std::pair<PointIndexType, PointIndexType> IndexPair;
837  typedef std::deque<IndexPair> IndexPairList;
838  typedef boost::shared_ptr<IndexPairList> IndexPairListPtr;
839  typedef std::map<Coord, IndexPairListPtr> IndexPairListMap;
840  typedef boost::shared_ptr<IndexPairListMap> IndexPairListMapPtr;
841 
842  size_t numTasks = 1, numThreads = size_t(tbb::task_scheduler_init::default_num_threads());
843  if (points.size() > (numThreads * 2)) numTasks = numThreads * 2;
844  else if (points.size() > numThreads) numTasks = numThreads;
845 
846  boost::scoped_array<IndexPairListMapPtr> bins(new IndexPairListMapPtr[numTasks]);
847 
849 
850  tbb::parallel_for(tbb::blocked_range<size_t>(0, numTasks),
851  BinOp(bins.get(), points, voxelOffsets, xform, binLog2Dim, bucketLog2Dim,
852  numTasks, cellCenteredTransform));
853 
854  std::set<Coord> uniqueCoords;
855 
856  for (size_t i = 0; i < numTasks; ++i) {
857  IndexPairListMap& idxMap = *bins[i];
858  for (typename IndexPairListMap::iterator it = idxMap.begin(); it != idxMap.end(); ++it) {
859  uniqueCoords.insert(it->first);
860  }
861  }
862 
863  std::vector<Coord> coords(uniqueCoords.begin(), uniqueCoords.end());
864  uniqueCoords.clear();
865 
866  segmentCount = coords.size();
867 
868  typedef typename Array<PointIndexType>::Ptr SegmentPtr;
869 
870  indexSegments.reset(new SegmentPtr[segmentCount]);
871  offsetSegments.reset(new SegmentPtr[segmentCount]);
872 
873  typedef MergeBinsOp<PointIndexType> MergeOp;
874 
875  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
876  MergeOp(bins.get(), indexSegments.get(), offsetSegments.get(), &coords[0], numTasks));
877 }
878 
879 
880 template<typename PointIndexType, typename VoxelOffsetType, typename PointArray>
881 inline void partition(
882  const PointArray& points,
883  const math::Transform& xform,
884  const Index bucketLog2Dim,
885  boost::scoped_array<PointIndexType>& pointIndices,
886  boost::scoped_array<PointIndexType>& pageOffsets,
887  PointIndexType& pageCount,
888  boost::scoped_array<VoxelOffsetType>& voxelOffsets,
889  bool recordVoxelOffsets,
890  bool cellCenteredTransform)
891 {
892  if (recordVoxelOffsets) voxelOffsets.reset(new VoxelOffsetType[points.size()]);
893  else voxelOffsets.reset();
894 
895  const Index binLog2Dim = 5u;
896  // note: Bins span a (2^(binLog2Dim + bucketLog2Dim))^3 voxel region,
897  // i.e. bucketLog2Dim = 3 and binLog2Dim = 5 corresponds to a
898  // (2^8)^3 = 256^3 voxel region.
899 
900 
901  size_t numSegments = 0;
902 
903  boost::scoped_array<typename Array<PointIndexType>::Ptr> indexSegments;
904  boost::scoped_array<typename Array<PointIndexType>::Ptr> offestSegments;
905 
906  binAndSegment<PointIndexType, VoxelOffsetType, PointArray>(points, xform,
907  indexSegments, offestSegments, numSegments, binLog2Dim, bucketLog2Dim,
908  voxelOffsets.get(), cellCenteredTransform);
909 
910  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
911 
912  typedef boost::scoped_array<PointIndexType> IndexArray;
913  boost::scoped_array<IndexArray> pageOffsetArrays(new IndexArray[numSegments]);
914 
915  const Index binVolume = 1u << (3u * binLog2Dim);
916 
917  tbb::parallel_for(segmentRange, OrderSegmentsOp<PointIndexType>
918  (indexSegments.get(), offestSegments.get(), pageOffsetArrays.get(), binVolume));
919 
920  indexSegments.reset();
921 
922  pageCount = 0;
923  for (size_t n = 0; n < numSegments; ++n) {
924  pageCount += pageOffsetArrays[n][0] - 1;
925  }
926 
927  pageOffsets.reset(new PointIndexType[pageCount + 1]);
928 
929  PointIndexType count = 0;
930  for (size_t n = 0, idx = 0; n < numSegments; ++n) {
931 
932  PointIndexType* offsets = pageOffsetArrays[n].get();
933  size_t size = size_t(offsets[0]);
934 
935  for (size_t i = 1; i < size; ++i) {
936  pageOffsets[idx++] = count;
937  count += offsets[i];
938  }
939  }
940 
941  pageOffsets[pageCount] = count;
942 
943  pointIndices.reset(new PointIndexType[points.size()]);
944 
945  std::vector<PointIndexType*> indexArray;
946  indexArray.reserve(numSegments);
947 
948  PointIndexType* index = pointIndices.get();
949  for (size_t n = 0; n < numSegments; ++n) {
950  indexArray.push_back(index);
951  index += offestSegments[n]->size();
952  }
953 
954  tbb::parallel_for(segmentRange, MoveSegmentDataOp<PointIndexType>(indexArray, offestSegments.get()));
955 }
956 
957 
958 } // namespace point_partitioner_internal
959 
960 
962 
963 
964 template<typename PointIndexType, Index BucketLog2Dim>
966  : mPointIndices(NULL)
967  , mVoxelOffsets(NULL)
968  , mPageOffsets(NULL)
969  , mPageCoordinates(NULL)
970  , mPageCount(0)
971  , mUsingCellCenteredTransform(true)
972 {
973 }
974 
975 
976 template<typename PointIndexType, Index BucketLog2Dim>
977 inline void
979 {
980  mPageCount = 0;
981  mUsingCellCenteredTransform = true;
982  mPointIndices.reset();
983  mVoxelOffsets.reset();
984  mPageOffsets.reset();
985  mPageCoordinates.reset();
986 }
987 
988 
989 template<typename PointIndexType, Index BucketLog2Dim>
990 inline void
992 {
993  const IndexType tmpLhsPageCount = mPageCount;
994  mPageCount = rhs.mPageCount;
995  rhs.mPageCount = tmpLhsPageCount;
996 
997  mPointIndices.swap(rhs.mPointIndices);
998  mVoxelOffsets.swap(rhs.mVoxelOffsets);
999  mPageOffsets.swap(rhs.mPageOffsets);
1000  mPageCoordinates.swap(rhs.mPageCoordinates);
1001 
1002  bool lhsCellCenteredTransform = mUsingCellCenteredTransform;
1003  mUsingCellCenteredTransform = rhs.mUsingCellCenteredTransform;
1004  rhs.mUsingCellCenteredTransform = lhsCellCenteredTransform;
1005 }
1006 
1007 
1008 template<typename PointIndexType, Index BucketLog2Dim>
1011 {
1012  assert(bool(mPointIndices) && bool(mPageCount));
1013  return IndexIterator(
1014  mPointIndices.get() + mPageOffsets[n],
1015  mPointIndices.get() + mPageOffsets[n + 1]);
1016 }
1017 
1018 
1019 template<typename PointIndexType, Index BucketLog2Dim>
1020 template<typename PointArray>
1021 inline void
1023  const math::Transform& xform, bool voxelOrder, bool recordVoxelOffsets, bool cellCenteredTransform)
1024 {
1025  mUsingCellCenteredTransform = cellCenteredTransform;
1026 
1027  point_partitioner_internal::partition(points, xform, BucketLog2Dim,
1028  mPointIndices, mPageOffsets, mPageCount, mVoxelOffsets,
1029  (voxelOrder || recordVoxelOffsets), cellCenteredTransform);
1030 
1031  const tbb::blocked_range<size_t> pageRange(0, mPageCount);
1032  mPageCoordinates.reset(new Coord[mPageCount]);
1033 
1034  tbb::parallel_for(pageRange,
1036  (mPageCoordinates, mPointIndices, mPageOffsets, points, xform,
1037  BucketLog2Dim, cellCenteredTransform));
1038 
1039  if (mVoxelOffsets && voxelOrder) {
1040  tbb::parallel_for(pageRange, point_partitioner_internal::VoxelOrderOp<
1041  IndexType, BucketLog2Dim>(mPointIndices, mPageOffsets, mVoxelOffsets));
1042  }
1043 
1044  if (mVoxelOffsets && !recordVoxelOffsets) {
1045  mVoxelOffsets.reset();
1046  }
1047 }
1048 
1049 
1050 template<typename PointIndexType, Index BucketLog2Dim>
1051 template<typename PointArray>
1054  bool voxelOrder, bool recordVoxelOffsets, bool cellCenteredTransform)
1055 {
1056  Ptr ret(new PointPartitioner());
1057  ret->construct(points, xform, voxelOrder, recordVoxelOffsets, cellCenteredTransform);
1058  return ret;
1059 }
1060 
1061 
1063 
1064 
1065 } // namespace tools
1066 } // namespace OPENVDB_VERSION_NAME
1067 } // namespace openvdb
1068 
1069 
1070 #endif // OPENVDB_TOOLS_POINT_PARTITIONER_HAS_BEEN_INCLUDED
1071 
1072 // Copyright (c) 2012-2016 DreamWorks Animation LLC
1073 // All rights reserved. This software is distributed under the
1074 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
Array< PointIndexType > Segment
Definition: PointPartitioner.h:496
size_t const mNumSegments
Definition: PointPartitioner.h:731
bool operator!=(const IndexIterator &other) const
Definition: PointPartitioner.h:235
PointIndexType const *const mIndices
Definition: PointPartitioner.h:294
std::deque< IndexPair > IndexPairList
Definition: PointPartitioner.h:500
Segment::Ptr SegmentPtr
Definition: PointPartitioner.h:453
VoxelOffsetType *const mVoxelOffsets
Definition: PointPartitioner.h:727
void clear()
Removes all data and frees up memory.
Definition: PointPartitioner.h:978
bool const mCellCenteredTransform
Definition: PointPartitioner.h:732
boost::int_t< 1+(3 *BucketLog2Dim)>::least VoxelOffsetType
Definition: PointPartitioner.h:301
Index const mBinVolume
Definition: PointPartitioner.h:816
boost::scoped_array< PointIndexType > IndexArray
Definition: PointPartitioner.h:739
void binAndSegment(const PointArray &points, const math::Transform &xform, boost::scoped_array< typename Array< PointIndexType >::Ptr > &indexSegments, boost::scoped_array< typename Array< PointIndexType >::Ptr > &offsetSegments, size_t &segmentCount, const Index binLog2Dim, const Index bucketLog2Dim, VoxelOffsetType *voxelOffsets=NULL, bool cellCenteredTransform=true)
Segment points using one level of least significant digit radix bins.
Definition: PointPartitioner.h:825
bool empty() const
true if the container size is 0, false otherwise.
Definition: PointPartitioner.h:157
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:460
IndexType & operator*()
Returns the item to which this iterator is currently pointing.
Definition: PointPartitioner.h:219
PointIndexType const *const mPointOrder
Definition: PointPartitioner.h:293
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:312
const T * data() const
Definition: PointPartitioner.h:439
bool const mCellCenteredTransform
Definition: PointPartitioner.h:422
boost::scoped_array< VoxelOffsetType > VoxelOffsetArray
Definition: PointPartitioner.h:302
OrderSegmentsOp(SegmentPtr *indexSegments, SegmentPtr *offestSegments, IndexArray *pageOffsetArrays, Index binVolume)
Definition: PointPartitioner.h:742
bool test() const
Definition: PointPartitioner.h:224
std::deque< IndexPair > IndexPairList
Definition: PointPartitioner.h:620
Array(size_t size)
Definition: PointPartitioner.h:434
PointIndexType IndexType
Definition: PointPartitioner.h:207
size_t size() const
Returns the number of buckets.
Definition: PointPartitioner.h:154
BinPointIndicesOp(IndexPairListMapPtr *data, const PointArray &points, VoxelOffsetType *voxelOffsets, const math::Transform &m, Index binLog2Dim, Index bucketLog2Dim, size_t numSegments, bool cellCenteredTransform)
Definition: PointPartitioner.h:625
Segment::Ptr SegmentPtr
Definition: PointPartitioner.h:497
const Coord & origin(size_t n) const
Returns the origin coordinate for bucket n.
Definition: PointPartitioner.h:174
PointIndexType const *const mBucketCounters
Definition: PointPartitioner.h:270
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:751
void clear()
Definition: PointPartitioner.h:441
PointIndexType *const mPointOrder
Definition: PointPartitioner.h:269
CoordBBox getBBox(size_t n) const
Returns the coordinate-aligned bounding box for bucket n.
Definition: PointPartitioner.h:169
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
Definition: PointPartitioner.h:107
PointIndexType const *const mBucketOffsets
Definition: PointPartitioner.h:271
bool usingCellCenteredTransform() const
Returns true if this point partitioning was constructed using a cell-centered transform.
Definition: PointPartitioner.h:183
PointIndexType const *const mPages
Definition: PointPartitioner.h:363
Coord *const mCoordinates
Definition: PointPartitioner.h:416
T * data()
Definition: PointPartitioner.h:438
VoxelOffsetType const *const mVoxelOffsets
Definition: PointPartitioner.h:364
Index32 Index
Definition: Types.h:58
boost::scoped_array< Coord > CoordArray
Definition: PointPartitioner.h:372
boost::shared_ptr< IndexPairListMap > IndexPairListMapPtr
Definition: PointPartitioner.h:503
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:387
math::Transform const mXForm
Definition: PointPartitioner.h:728
IndexIterator & operator++()
Advance to the next item.
Definition: PointPartitioner.h:227
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:1053
#define OPENVDB_VERSION_NAME
Definition: version.h:43
SegmentPtr *const mOffsetSegments
Definition: PointPartitioner.h:814
boost::scoped_array< PointIndexType > IndexArray
Definition: PointPartitioner.h:303
CreateOrderedPointIndexArrayOp(PointIndexType *orderedIndexArray, const PointIndexType *pointOrder, const PointIndexType *indices)
Definition: PointPartitioner.h:278
Index const mBinLog2Dim
Definition: PointPartitioner.h:729
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:518
VoxelOrderOp(IndexArray &indices, const IndexArray &pages, const VoxelOffsetArray &offsets)
Definition: PointPartitioner.h:305
void reset()
Rewind to first item.
Definition: PointPartitioner.h:213
boost::scoped_array< PointIndexType > IndexArray
Definition: PointPartitioner.h:371
PointPartitioner< uint32_t, 3 > UInt32PointPartitioner
Definition: PointPartitioner.h:200
boost::shared_ptr< IndexPairListMap > IndexPairListMapPtr
Definition: PointPartitioner.h:623
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
LeafNodeOriginOp(CoordArray &coordinates, const IndexArray &indices, const IndexArray &pages, const PointArray &points, const math::Transform &m, int log2dim, bool cellCenteredTransform)
Definition: PointPartitioner.h:374
std::map< Coord, IndexPairListPtr > IndexPairListMap
Definition: PointPartitioner.h:622
std::pair< PointIndexType, PointIndexType > IndexPair
Definition: PointPartitioner.h:499
Definition: Exceptions.h:39
boost::shared_ptr< Array > Ptr
Definition: PointPartitioner.h:432
bool operator==(const IndexIterator &other) const
Equality operators.
Definition: PointPartitioner.h:234
const VoxelOffsetArray & voxelOffsets() const
Returns a list of LeafNode voxel offsets for the points.
Definition: PointPartitioner.h:178
IndexIterator(IndexType *begin=NULL, IndexType *end=NULL)
Definition: PointPartitioner.h:209
boost::shared_ptr< IndexPairList > IndexPairListPtr
Definition: PointPartitioner.h:621
MergeBinsOp(IndexPairListMapPtr *bins, SegmentPtr *indexSegments, SegmentPtr *offsetSegments, Coord *coords, size_t numSegments)
Definition: PointPartitioner.h:505
PointIndexType const *const mPages
Definition: PointPartitioner.h:418
boost::shared_ptr< IndexPairList > IndexPairListPtr
Definition: PointPartitioner.h:501
PointArray::PosType PosType
Definition: PointPartitioner.h:618
PointIndexType IndexType
Definition: PointPartitioner.h:115
boost::shared_ptr< const PointPartitioner > ConstPtr
Definition: PointPartitioner.h:113
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:644
math::Transform const mXForm
Definition: PointPartitioner.h:420
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:1022
int const mLog2Dim
Definition: PointPartitioner.h:421
Array< PointIndexType > Segment
Definition: PointPartitioner.h:452
void partition(const PointArray &points, const math::Transform &xform, const Index bucketLog2Dim, boost::scoped_array< PointIndexType > &pointIndices, boost::scoped_array< PointIndexType > &pageOffsets, PointIndexType &pageCount, boost::scoped_array< VoxelOffsetType > &voxelOffsets, bool recordVoxelOffsets, bool cellCenteredTransform)
Definition: PointPartitioner.h:881
size_t size() const
Definition: PointPartitioner.h:436
Array< PointIndexType >::Ptr SegmentPtr
Definition: PointPartitioner.h:740
PointArray const *const mPoints
Definition: PointPartitioner.h:419
bool operator==(const Vec3< T0 > &v0, const Vec3< T1 > &v1)
Equality operator, does exact floating point comparisons.
Definition: Vec3.h:450
const IndexType & operator*() const
Definition: PointPartitioner.h:220
void swap(PointPartitioner &)
Exchanges the content of the container by another.
Definition: PointPartitioner.h:991
MoveSegmentDataOp(std::vector< PointIndexType * > &indexLists, SegmentPtr *segments)
Definition: PointPartitioner.h:455
IndexPairListMapPtr *const mData
Definition: PointPartitioner.h:725
PointIndexType const *const mIndices
Definition: PointPartitioner.h:417
PointIndexType *const mOrderedIndexArray
Definition: PointPartitioner.h:292
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:263
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: PointPartitioner.h:286
IndexArray *const mPageOffsetArrays
Definition: PointPartitioner.h:815
PointIndexType *const mIndices
Definition: PointPartitioner.h:362
boost::shared_ptr< PointPartitioner > Ptr
Definition: PointPartitioner.h:112
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:368
std::map< Coord, IndexPairListPtr > IndexPairListMap
Definition: PointPartitioner.h:502
IndexIterator indices(size_t n) const
Returns the point indices for bucket n.
Definition: PointPartitioner.h:1010
bool increment()
Definition: PointPartitioner.h:231
boost::int_t< 1+(3 *BucketLog2Dim)>::least VoxelOffsetType
Definition: PointPartitioner.h:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
std::pair< PointIndexType, PointIndexType > IndexPair
Definition: PointPartitioner.h:619
ComputePointOrderOp(PointIndexType *pointOrder, const PointIndexType *bucketCounters, const PointIndexType *bucketOffsets)
Definition: PointPartitioner.h:255
PointPartitioner()
Definition: PointPartitioner.h:965
size_t size() const
Number of point indices in the iterator range.
Definition: PointPartitioner.h:216
PointArray const *const mPoints
Definition: PointPartitioner.h:726
bool next()
Advance to the next item.
Definition: PointPartitioner.h:230
boost::scoped_array< VoxelOffsetType > VoxelOffsetArray
Definition: PointPartitioner.h:117
SegmentPtr *const mIndexSegments
Definition: PointPartitioner.h:813
Index const mBucketLog2Dim
Definition: PointPartitioner.h:730