31 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 32 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 34 #include <openvdb/tree/ValueAccessor.h> 35 #include <openvdb/tree/LeafManager.h> 36 #include <openvdb/tools/Morphology.h> 38 #include <openvdb/tools/PointScatter.h> 39 #include <openvdb/tools/LevelSetUtil.h> 40 #include <openvdb/tools/VolumeToMesh.h> 42 #include <boost/scoped_array.hpp> 43 #include <boost/scoped_ptr.hpp> 44 #include <tbb/blocked_range.h> 45 #include <tbb/parallel_for.h> 46 #include <tbb/parallel_reduce.h> 89 template<
typename Gr
idT,
typename InterrupterT>
93 std::vector<openvdb::Vec4s>& spheres,
95 bool overlapping =
false,
96 float minRadius = 1.0,
99 int instanceCount = 10000,
100 InterrupterT* interrupter = NULL);
105 template<
typename Gr
idT>
109 std::vector<openvdb::Vec4s>& spheres,
111 bool overlapping =
false,
112 float minRadius = 1.0,
114 float isovalue = 0.0,
115 int instanceCount = 10000)
117 fillWithSpheres<GridT, util::NullInterrupter>(grid, spheres,
118 maxSphereCount, overlapping, minRadius, maxRadius, isovalue, instanceCount);
128 template<
typename Gr
idT>
132 typedef typename GridT::TreeType
TreeT;
133 typedef typename TreeT::template ValueConverter<bool>::Type
BoolTreeT;
134 typedef typename TreeT::template ValueConverter<Index32>::Type
Index32TreeT;
135 typedef typename TreeT::template ValueConverter<Int16>::Type
Int16TreeT;
152 template<
typename InterrupterT>
153 void initialize(
const GridT& grid,
float isovalue = 0.0, InterrupterT* interrupter = NULL);
158 void initialize(
const GridT& grid,
float isovalue = 0.0);
168 bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
178 bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
183 const Index32TreeT&
indexTree()
const {
return *mIdxTreePt; }
184 const Int16TreeT&
signTree()
const {
return *mSignTreePt; }
188 typedef typename Index32TreeT::LeafNodeType Index32LeafT;
189 typedef std::pair<size_t, size_t> IndexRange;
192 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
193 std::vector<IndexRange> mLeafRanges;
194 std::vector<const Index32LeafT*> mLeafNodes;
196 size_t mPointListSize, mMaxNodeLeafs;
198 typename Index32TreeT::Ptr mIdxTreePt;
199 typename Int16TreeT::Ptr mSignTreePt;
201 bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
224 mPoints.push_back(pos);
227 std::vector<Vec3R>& mPoints;
231 template<
typename Index32LeafT>
236 LeafBS(std::vector<Vec4R>& leafBoundingSpheres,
237 const std::vector<const Index32LeafT*>& leafNodes,
241 void run(
bool threaded =
true);
244 void operator()(
const tbb::blocked_range<size_t>&)
const;
247 std::vector<Vec4R>& mLeafBoundingSpheres;
248 const std::vector<const Index32LeafT*>& mLeafNodes;
253 template<
typename Index32LeafT>
255 std::vector<Vec4R>& leafBoundingSpheres,
256 const std::vector<const Index32LeafT*>& leafNodes,
259 : mLeafBoundingSpheres(leafBoundingSpheres)
260 , mLeafNodes(leafNodes)
261 , mTransform(transform)
262 , mSurfacePointList(surfacePointList)
266 template<
typename Index32LeafT>
271 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
273 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
277 template<
typename Index32LeafT>
281 typename Index32LeafT::ValueOnCIter iter;
284 for (
size_t n = range.begin(); n != range.end(); ++n) {
291 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
292 avg += mSurfacePointList[iter.getValue()];
296 if (count > 1) avg *= float(1.0 /
double(count));
300 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
301 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
302 if (tmpDist > maxDist) maxDist = tmpDist;
305 Vec4R& sphere = mLeafBoundingSpheres[n];
310 sphere[3] = maxDist * 2.0;
320 NodeBS(std::vector<Vec4R>& nodeBoundingSpheres,
321 const std::vector<IndexRange>& leafRanges,
322 const std::vector<Vec4R>& leafBoundingSpheres);
324 inline void run(
bool threaded =
true);
326 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
329 std::vector<Vec4R>& mNodeBoundingSpheres;
330 const std::vector<IndexRange>& mLeafRanges;
331 const std::vector<Vec4R>& mLeafBoundingSpheres;
336 const std::vector<IndexRange>& leafRanges,
337 const std::vector<Vec4R>& leafBoundingSpheres)
338 : mNodeBoundingSpheres(nodeBoundingSpheres)
339 , mLeafRanges(leafRanges)
340 , mLeafBoundingSpheres(leafBoundingSpheres)
348 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
350 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
359 for (
size_t n = range.begin(); n != range.end(); ++n) {
365 int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
367 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
368 avg[0] += mLeafBoundingSpheres[i][0];
369 avg[1] += mLeafBoundingSpheres[i][1];
370 avg[2] += mLeafBoundingSpheres[i][2];
373 if (count > 1) avg *= float(1.0 /
double(count));
376 double maxDist = 0.0;
378 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
379 pos[0] = mLeafBoundingSpheres[i][0];
380 pos[1] = mLeafBoundingSpheres[i][1];
381 pos[2] = mLeafBoundingSpheres[i][2];
383 double tmpDist = (pos - avg).lengthSqr() + mLeafBoundingSpheres[i][3];
384 if (tmpDist > maxDist) maxDist = tmpDist;
387 Vec4R& sphere = mNodeBoundingSpheres[n];
392 sphere[3] = maxDist * 2.0;
401 template<
typename Index32LeafT>
408 std::vector<Vec3R>& instancePoints,
409 std::vector<float>& instanceDistances,
411 const std::vector<const Index32LeafT*>& leafNodes,
412 const std::vector<IndexRange>& leafRanges,
413 const std::vector<Vec4R>& leafBoundingSpheres,
414 const std::vector<Vec4R>& nodeBoundingSpheres,
416 bool transformPoints =
false);
419 void run(
bool threaded =
true);
422 void operator()(
const tbb::blocked_range<size_t>&)
const;
426 void evalLeaf(
size_t index,
const Index32LeafT& leaf)
const;
427 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
430 std::vector<Vec3R>& mInstancePoints;
431 std::vector<float>& mInstanceDistances;
435 const std::vector<const Index32LeafT*>& mLeafNodes;
436 const std::vector<IndexRange>& mLeafRanges;
437 const std::vector<Vec4R>& mLeafBoundingSpheres;
438 const std::vector<Vec4R>& mNodeBoundingSpheres;
440 std::vector<float> mLeafDistances, mNodeDistances;
442 const bool mTransformPoints;
443 size_t mClosestPointIndex;
447 template<
typename Index32LeafT>
449 std::vector<Vec3R>& instancePoints,
450 std::vector<float>& instanceDistances,
452 const std::vector<const Index32LeafT*>& leafNodes,
453 const std::vector<IndexRange>& leafRanges,
454 const std::vector<Vec4R>& leafBoundingSpheres,
455 const std::vector<Vec4R>& nodeBoundingSpheres,
457 bool transformPoints)
458 : mInstancePoints(instancePoints)
459 , mInstanceDistances(instanceDistances)
460 , mSurfacePointList(surfacePointList)
461 , mLeafNodes(leafNodes)
462 , mLeafRanges(leafRanges)
463 , mLeafBoundingSpheres(leafBoundingSpheres)
464 , mNodeBoundingSpheres(nodeBoundingSpheres)
465 , mLeafDistances(maxNodeLeafs, 0.0)
466 , mNodeDistances(leafRanges.size(), 0.0)
467 , mTransformPoints(transformPoints)
468 , mClosestPointIndex(0)
473 template<
typename Index32LeafT>
478 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
480 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
484 template<
typename Index32LeafT>
488 typename Index32LeafT::ValueOnCIter iter;
489 const Vec3s center = mInstancePoints[index];
490 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
492 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
494 const Vec3s& point = mSurfacePointList[iter.getValue()];
495 float tmpDist = (point - center).lengthSqr();
497 if (tmpDist < mInstanceDistances[index]) {
498 mInstanceDistances[index] = tmpDist;
499 closestPointIndex = iter.getValue();
505 template<
typename Index32LeafT>
509 const Vec3R& pos = mInstancePoints[pointIndex];
510 float minDist = mInstanceDistances[pointIndex];
511 size_t minDistIdx = 0;
513 bool updatedDist =
false;
515 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
516 i < mLeafRanges[nodeIndex].second; ++i, ++n)
518 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
520 center[0] = mLeafBoundingSpheres[i][0];
521 center[1] = mLeafBoundingSpheres[i][1];
522 center[2] = mLeafBoundingSpheres[i][2];
524 distToLeaf = float((pos - center).lengthSqr() - mLeafBoundingSpheres[i][3]);
526 if (distToLeaf < minDist) {
527 minDist = distToLeaf;
533 if (!updatedDist)
return;
535 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
537 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
538 i < mLeafRanges[nodeIndex].second; ++i, ++n)
540 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
541 evalLeaf(pointIndex, *mLeafNodes[i]);
547 template<
typename Index32LeafT>
552 for (
size_t n = range.begin(); n != range.end(); ++n) {
554 const Vec3R& pos = mInstancePoints[n];
555 float minDist = mInstanceDistances[n];
556 size_t minDistIdx = 0;
558 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
559 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
561 center[0] = mNodeBoundingSpheres[i][0];
562 center[1] = mNodeBoundingSpheres[i][1];
563 center[2] = mNodeBoundingSpheres[i][2];
565 distToNode = float((pos - center).lengthSqr() - mNodeBoundingSpheres[i][3]);
567 if (distToNode < minDist) {
568 minDist = distToNode;
573 evalNode(n, minDistIdx);
575 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
576 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
581 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
583 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
593 const std::vector<Vec3R>& points,
594 std::vector<float>& distances,
595 std::vector<unsigned char>& mask,
599 int index()
const {
return mIndex; }
601 inline void run(
bool threaded =
true);
605 inline void operator()(
const tbb::blocked_range<size_t>& range);
608 if (rhs.mRadius > mRadius) {
609 mRadius = rhs.mRadius;
616 const Vec4s& mSphere;
617 const std::vector<Vec3R>& mPoints;
619 std::vector<float>& mDistances;
620 std::vector<unsigned char>& mMask;
630 const std::vector<Vec3R>& points,
631 std::vector<float>& distances,
632 std::vector<unsigned char>& mask,
636 , mDistances(distances)
638 , mOverlapping(overlapping)
646 : mSphere(rhs.mSphere)
647 , mPoints(rhs.mPoints)
648 , mDistances(rhs.mDistances)
650 , mOverlapping(rhs.mOverlapping)
651 , mRadius(rhs.mRadius)
660 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
662 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
670 for (
size_t n = range.begin(); n != range.end(); ++n) {
671 if (mMask[n])
continue;
673 pos.x() = float(mPoints[n].x()) - mSphere[0];
674 pos.y() = float(mPoints[n].y()) - mSphere[1];
675 pos.z() = float(mPoints[n].z()) - mSphere[2];
677 float dist = pos.length();
679 if (dist < mSphere[3]) {
685 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
688 if (mDistances[n] > mRadius) {
689 mRadius = mDistances[n];
702 template<
typename Gr
idT,
typename InterrupterT>
706 std::vector<openvdb::Vec4s>& spheres,
713 InterrupterT* interrupter)
716 spheres.reserve(maxSphereCount);
718 const bool addNBPoints = grid.activeVoxelCount() < 10000;
719 int instances =
std::max(instanceCount, maxSphereCount);
721 typedef typename GridT::TreeType TreeT;
722 typedef typename GridT::ValueType ValueT;
724 typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
725 typedef typename TreeT::template ValueConverter<Index32>::Type Index32TreeT;
726 typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
728 typedef boost::mt11213b RandGen;
731 const TreeT& tree = grid.tree();
734 std::vector<Vec3R> instancePoints;
744 interiorMaskPtr->
tree().topologyUnion(tree);
747 if (interrupter && interrupter->wasInterrupted())
return;
751 instancePoints.reserve(instances);
755 ptnAcc,
Index64(addNBPoints ? (instances / 2) : instances), mtRand, interrupter);
757 scatter(*interiorMaskPtr);
760 if (interrupter && interrupter->wasInterrupted())
return;
762 std::vector<float> instanceRadius;
768 if (instancePoints.size() < size_t(instances)) {
769 const Int16TreeT& signTree = csp.
signTree();
770 typename Int16TreeT::LeafNodeType::ValueOnCIter it;
771 typename Int16TreeT::LeafCIter leafIt = signTree.cbeginLeaf();
773 for (; leafIt; ++leafIt) {
774 for (it = leafIt->cbeginValueOn(); it; ++it) {
776 const int flags = int(it.getValue());
778 instancePoints.push_back(transform.
indexToWorld(it.getCoord()));
781 if (instancePoints.size() == size_t(instances))
break;
783 if (instancePoints.size() == size_t(instances))
break;
788 if (interrupter && interrupter->wasInterrupted())
return;
790 if (!csp.
search(instancePoints, instanceRadius))
return;
792 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
793 float largestRadius = 0.0;
794 int largestRadiusIdx = 0;
796 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
797 if (instanceRadius[n] > largestRadius) {
798 largestRadius = instanceRadius[n];
799 largestRadiusIdx = int(n);
805 minRadius = float(minRadius * transform.
voxelSize()[0]);
806 maxRadius = float(maxRadius * transform.
voxelSize()[0]);
808 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
810 if (interrupter && interrupter->wasInterrupted())
return;
812 largestRadius =
std::min(maxRadius, largestRadius);
814 if (s != 0 && largestRadius < minRadius)
break;
816 sphere[0] = float(instancePoints[largestRadiusIdx].x());
817 sphere[1] = float(instancePoints[largestRadiusIdx].y());
818 sphere[2] = float(instancePoints[largestRadiusIdx].z());
819 sphere[3] = largestRadius;
821 spheres.push_back(sphere);
822 instanceMask[largestRadiusIdx] = 1;
825 sphere, instancePoints, instanceRadius, instanceMask, overlapping);
828 largestRadius = op.
radius();
829 largestRadiusIdx = op.
index();
836 template<
typename Gr
idT>
838 : mIsInitialized(false)
839 , mLeafBoundingSpheres(0)
840 , mNodeBoundingSpheres(0)
843 , mSurfacePointList()
851 template<
typename Gr
idT>
855 initialize<GridT, util::NullInterrupter>(grid, isovalue, NULL);
859 template<
typename Gr
idT>
860 template<
typename InterrupterT>
863 const GridT& grid,
float isovalue, InterrupterT* interrupter)
865 mIsInitialized =
false;
869 typedef typename GridT::ValueType ValueT;
877 mIdxTreePt.reset(
new Index32TreeT(boost::integer_traits<Index32>::const_max));
883 if (interrupter && interrupter->wasInterrupted())
return;
887 typedef typename Int16TreeT::LeafNodeType Int16LeafNodeType;
888 typedef typename Index32TreeT::LeafNodeType Index32LeafNodeType;
890 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
891 mSignTreePt->getNodes(signFlagsLeafNodes);
893 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
895 boost::scoped_array<Index32> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
897 tbb::parallel_for(auxiliaryLeafNodeRange,
899 (signFlagsLeafNodes, leafNodeOffsets));
903 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
904 const Index32 tmp = leafNodeOffsets[n];
905 leafNodeOffsets[n] = pointCount;
909 mPointListSize = size_t(pointCount);
910 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
914 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
915 mIdxTreePt->getNodes(pointIndexLeafNodes);
917 tbb::parallel_for(auxiliaryLeafNodeRange,
919 signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
922 if (interrupter && interrupter->wasInterrupted())
return;
925 CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
930 dim[0] = std::abs(dim[0]);
931 dim[1] = std::abs(dim[1]);
932 dim[2] = std::abs(dim[2]);
935 mMaxRadiusSqr *= 0.51f;
936 mMaxRadiusSqr *= mMaxRadiusSqr;
939 Index32LeafManagerT idxLeafs(*mIdxTreePt);
941 typedef typename Index32TreeT::RootNodeType Index32RootNodeT;
942 typedef typename Index32RootNodeT::NodeChainType Index32NodeChainT;
943 BOOST_STATIC_ASSERT(boost::mpl::size<Index32NodeChainT>::value > 1);
944 typedef typename boost::mpl::at<Index32NodeChainT, boost::mpl::int_<1> >::type Index32InternalNodeT;
946 typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
947 nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
948 nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
950 std::vector<const Index32InternalNodeT*> internalNodes;
952 const Index32InternalNodeT* node = NULL;
955 if (node) internalNodes.push_back(node);
958 std::vector<IndexRange>().swap(mLeafRanges);
959 mLeafRanges.resize(internalNodes.size());
961 std::vector<const Index32LeafT*>().swap(mLeafNodes);
962 mLeafNodes.reserve(idxLeafs.leafCount());
964 typename Index32InternalNodeT::ChildOnCIter leafIt;
966 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
968 mLeafRanges[n].first = mLeafNodes.size();
970 size_t leafCount = 0;
971 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
972 mLeafNodes.push_back(&(*leafIt));
976 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
978 mLeafRanges[n].second = mLeafNodes.size();
981 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
982 mLeafBoundingSpheres.resize(mLeafNodes.size());
985 mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
989 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
990 mNodeBoundingSpheres.resize(internalNodes.size());
992 internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
994 mIsInitialized =
true;
998 template<
typename Gr
idT>
1001 std::vector<float>& distances,
bool transformPoints)
1003 if (!mIsInitialized)
return false;
1006 distances.resize(points.size(), mMaxRadiusSqr);
1009 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1010 mMaxNodeLeafs, transformPoints);
1018 template<
typename Gr
idT>
1022 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
1026 template<
typename Gr
idT>
1029 std::vector<float>& distances)
1031 return search(points, distances,
true);
1039 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
Vec3< float > Vec3s
Definition: Vec3.h:650
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:115
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:761
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:338
uint64_t Index64
Definition: Types.h:57
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3< double > Vec3d
Definition: Vec3.h:651
Vec4< float > Vec4s
Definition: Vec4.h:583
Definition: Exceptions.h:39
uint32_t Index32
Definition: Types.h:56
OPENVDB_IMPORT void initialize()
Global registration of basic types.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1039
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54