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<int>::Type
IntTreeT;
134 typedef typename TreeT::template ValueConverter<Int16>::Type
Int16TreeT;
151 template<
typename InterrupterT>
152 void initialize(
const GridT& grid,
float isovalue = 0.0, InterrupterT* interrupter = NULL);
157 void initialize(
const GridT& grid,
float isovalue = 0.0);
167 bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
177 bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
182 const IntTreeT&
indexTree()
const {
return *mIdxTreePt; }
183 const Int16TreeT&
signTree()
const {
return *mSignTreePt; }
187 typedef typename IntTreeT::LeafNodeType IntLeafT;
188 typedef std::pair<size_t, size_t> IndexRange;
191 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
192 std::vector<IndexRange> mLeafRanges;
193 std::vector<const IntLeafT*> mLeafNodes;
195 size_t mPointListSize, mMaxNodeLeafs;
197 typename IntTreeT::Ptr mIdxTreePt;
198 typename Int16TreeT::Ptr mSignTreePt;
200 bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
223 mPoints.push_back(pos);
226 std::vector<Vec3R>& mPoints;
230 template<
typename IntLeafT>
235 LeafBS(std::vector<Vec4R>& leafBoundingSpheres,
236 const std::vector<const IntLeafT*>& leafNodes,
240 void run(
bool threaded =
true);
243 void operator()(
const tbb::blocked_range<size_t>&)
const;
246 std::vector<Vec4R>& mLeafBoundingSpheres;
247 const std::vector<const IntLeafT*>& mLeafNodes;
252 template<
typename IntLeafT>
254 std::vector<Vec4R>& leafBoundingSpheres,
255 const std::vector<const IntLeafT*>& leafNodes,
258 : mLeafBoundingSpheres(leafBoundingSpheres)
259 , mLeafNodes(leafNodes)
260 , mTransform(transform)
261 , mSurfacePointList(surfacePointList)
265 template<
typename IntLeafT>
270 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
272 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
276 template<
typename IntLeafT>
280 typename IntLeafT::ValueOnCIter iter;
283 for (
size_t n = range.begin(); n != range.end(); ++n) {
290 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
291 avg += mSurfacePointList[iter.getValue()];
295 if (count > 1) avg *= float(1.0 /
double(count));
299 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
300 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
301 if (tmpDist > maxDist) maxDist = tmpDist;
304 Vec4R& sphere = mLeafBoundingSpheres[n];
309 sphere[3] = maxDist * 2.0;
319 NodeBS(std::vector<Vec4R>& nodeBoundingSpheres,
320 const std::vector<IndexRange>& leafRanges,
321 const std::vector<Vec4R>& leafBoundingSpheres);
323 inline void run(
bool threaded =
true);
325 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
328 std::vector<Vec4R>& mNodeBoundingSpheres;
329 const std::vector<IndexRange>& mLeafRanges;
330 const std::vector<Vec4R>& mLeafBoundingSpheres;
335 const std::vector<IndexRange>& leafRanges,
336 const std::vector<Vec4R>& leafBoundingSpheres)
337 : mNodeBoundingSpheres(nodeBoundingSpheres)
338 , mLeafRanges(leafRanges)
339 , mLeafBoundingSpheres(leafBoundingSpheres)
347 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
349 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
358 for (
size_t n = range.begin(); n != range.end(); ++n) {
364 int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
366 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
367 avg[0] += mLeafBoundingSpheres[i][0];
368 avg[1] += mLeafBoundingSpheres[i][1];
369 avg[2] += mLeafBoundingSpheres[i][2];
372 if (count > 1) avg *= float(1.0 /
double(count));
375 double maxDist = 0.0;
377 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
378 pos[0] = mLeafBoundingSpheres[i][0];
379 pos[1] = mLeafBoundingSpheres[i][1];
380 pos[2] = mLeafBoundingSpheres[i][2];
382 double tmpDist = (pos - avg).lengthSqr() + mLeafBoundingSpheres[i][3];
383 if (tmpDist > maxDist) maxDist = tmpDist;
386 Vec4R& sphere = mNodeBoundingSpheres[n];
391 sphere[3] = maxDist * 2.0;
400 template<
typename IntLeafT>
407 std::vector<Vec3R>& instancePoints,
408 std::vector<float>& instanceDistances,
410 const std::vector<const IntLeafT*>& leafNodes,
411 const std::vector<IndexRange>& leafRanges,
412 const std::vector<Vec4R>& leafBoundingSpheres,
413 const std::vector<Vec4R>& nodeBoundingSpheres,
415 bool transformPoints =
false);
418 void run(
bool threaded =
true);
421 void operator()(
const tbb::blocked_range<size_t>&)
const;
425 void evalLeaf(
size_t index,
const IntLeafT& leaf)
const;
426 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
429 std::vector<Vec3R>& mInstancePoints;
430 std::vector<float>& mInstanceDistances;
434 const std::vector<const IntLeafT*>& mLeafNodes;
435 const std::vector<IndexRange>& mLeafRanges;
436 const std::vector<Vec4R>& mLeafBoundingSpheres;
437 const std::vector<Vec4R>& mNodeBoundingSpheres;
439 std::vector<float> mLeafDistances, mNodeDistances;
441 const bool mTransformPoints;
442 size_t mClosestPointIndex;
446 template<
typename IntLeafT>
448 std::vector<Vec3R>& instancePoints,
449 std::vector<float>& instanceDistances,
451 const std::vector<const IntLeafT*>& leafNodes,
452 const std::vector<IndexRange>& leafRanges,
453 const std::vector<Vec4R>& leafBoundingSpheres,
454 const std::vector<Vec4R>& nodeBoundingSpheres,
456 bool transformPoints)
457 : mInstancePoints(instancePoints)
458 , mInstanceDistances(instanceDistances)
459 , mSurfacePointList(surfacePointList)
460 , mLeafNodes(leafNodes)
461 , mLeafRanges(leafRanges)
462 , mLeafBoundingSpheres(leafBoundingSpheres)
463 , mNodeBoundingSpheres(nodeBoundingSpheres)
464 , mLeafDistances(maxNodeLeafs, 0.0)
465 , mNodeDistances(leafRanges.size(), 0.0)
466 , mTransformPoints(transformPoints)
467 , mClosestPointIndex(0)
472 template<
typename IntLeafT>
477 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
479 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
483 template<
typename IntLeafT>
487 typename IntLeafT::ValueOnCIter iter;
488 const Vec3s center = mInstancePoints[index];
489 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
491 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
493 const Vec3s& point = mSurfacePointList[iter.getValue()];
494 float tmpDist = (point - center).lengthSqr();
496 if (tmpDist < mInstanceDistances[index]) {
497 mInstanceDistances[index] = tmpDist;
498 closestPointIndex = iter.getValue();
504 template<
typename IntLeafT>
506 ClosestPointDist<IntLeafT>::evalNode(
size_t pointIndex,
size_t nodeIndex)
const
508 const Vec3R& pos = mInstancePoints[pointIndex];
509 float minDist = mInstanceDistances[pointIndex];
510 size_t minDistIdx = 0;
512 bool updatedDist =
false;
514 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
515 i < mLeafRanges[nodeIndex].second; ++i, ++n)
517 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
519 center[0] = mLeafBoundingSpheres[i][0];
520 center[1] = mLeafBoundingSpheres[i][1];
521 center[2] = mLeafBoundingSpheres[i][2];
523 distToLeaf = float((pos - center).lengthSqr() - mLeafBoundingSpheres[i][3]);
525 if (distToLeaf < minDist) {
526 minDist = distToLeaf;
532 if (!updatedDist)
return;
534 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
536 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
537 i < mLeafRanges[nodeIndex].second; ++i, ++n)
539 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
540 evalLeaf(pointIndex, *mLeafNodes[i]);
546 template<
typename IntLeafT>
551 for (
size_t n = range.begin(); n != range.end(); ++n) {
553 const Vec3R& pos = mInstancePoints[n];
554 float minDist = mInstanceDistances[n];
555 size_t minDistIdx = 0;
557 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
558 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
560 center[0] = mNodeBoundingSpheres[i][0];
561 center[1] = mNodeBoundingSpheres[i][1];
562 center[2] = mNodeBoundingSpheres[i][2];
564 distToNode = float((pos - center).lengthSqr() - mNodeBoundingSpheres[i][3]);
566 if (distToNode < minDist) {
567 minDist = distToNode;
572 evalNode(n, minDistIdx);
574 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
575 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
580 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
582 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
592 const std::vector<Vec3R>& points,
593 std::vector<float>& distances,
594 std::vector<unsigned char>& mask,
598 int index()
const {
return mIndex; };
600 inline void run(
bool threaded =
true);
604 inline void operator()(
const tbb::blocked_range<size_t>& range);
607 if (rhs.mRadius > mRadius) {
608 mRadius = rhs.mRadius;
615 const Vec4s& mSphere;
616 const std::vector<Vec3R>& mPoints;
618 std::vector<float>& mDistances;
619 std::vector<unsigned char>& mMask;
629 const std::vector<Vec3R>& points,
630 std::vector<float>& distances,
631 std::vector<unsigned char>& mask,
635 , mDistances(distances)
637 , mOverlapping(overlapping)
645 : mSphere(rhs.mSphere)
646 , mPoints(rhs.mPoints)
647 , mDistances(rhs.mDistances)
649 , mOverlapping(rhs.mOverlapping)
650 , mRadius(rhs.mRadius)
659 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
661 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
669 for (
size_t n = range.begin(); n != range.end(); ++n) {
670 if (mMask[n])
continue;
672 pos.
x() = float(mPoints[n].x()) - mSphere[0];
673 pos.y() = float(mPoints[n].y()) - mSphere[1];
674 pos.z() = float(mPoints[n].z()) - mSphere[2];
676 float dist = pos.length();
678 if (dist < mSphere[3]) {
684 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
687 if (mDistances[n] > mRadius) {
688 mRadius = mDistances[n];
701 template<
typename Gr
idT,
typename InterrupterT>
705 std::vector<openvdb::Vec4s>& spheres,
712 InterrupterT* interrupter)
715 spheres.reserve(maxSphereCount);
717 const bool addNBPoints = grid.activeVoxelCount() < 10000;
718 int instances =
std::max(instanceCount, maxSphereCount);
720 typedef typename GridT::TreeType TreeT;
721 typedef typename GridT::ValueType ValueT;
723 typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
724 typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
725 typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
732 typedef boost::mt11213b RandGen;
735 const TreeT& tree = grid.tree();
738 std::vector<Vec3R> instancePoints;
748 interiorMaskPtr->
tree().topologyUnion(tree);
751 if (interrupter && interrupter->wasInterrupted())
return;
755 instancePoints.reserve(instances);
759 ptnAcc,
Index64(addNBPoints ? (instances / 2) : instances), mtRand, interrupter);
761 scatter(*interiorMaskPtr);
764 if (interrupter && interrupter->wasInterrupted())
return;
766 std::vector<float> instanceRadius;
772 if (instancePoints.size() < size_t(instances)) {
773 const Int16TreeT& signTree = csp.
signTree();
774 typename Int16TreeT::LeafNodeType::ValueOnCIter it;
775 typename Int16TreeT::LeafCIter leafIt = signTree.cbeginLeaf();
777 for (; leafIt; ++leafIt) {
778 for (it = leafIt->cbeginValueOn(); it; ++it) {
779 const int flags = it.getValue();
780 if (!(0xE00 & flags) && (flags & 0x100)) {
781 instancePoints.push_back(transform.
indexToWorld(it.getCoord()));
784 if (instancePoints.size() == size_t(instances))
break;
786 if (instancePoints.size() == size_t(instances))
break;
791 if (interrupter && interrupter->wasInterrupted())
return;
793 if (!csp.
search(instancePoints, instanceRadius))
return;
795 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
796 float largestRadius = 0.0;
797 int largestRadiusIdx = 0;
799 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
800 if (instanceRadius[n] > largestRadius) {
801 largestRadius = instanceRadius[n];
802 largestRadiusIdx = int(n);
808 minRadius = float(minRadius * transform.
voxelSize()[0]);
809 maxRadius = float(maxRadius * transform.
voxelSize()[0]);
811 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
813 if (interrupter && interrupter->wasInterrupted())
return;
815 largestRadius =
std::min(maxRadius, largestRadius);
817 if (s != 0 && largestRadius < minRadius)
break;
819 sphere[0] = float(instancePoints[largestRadiusIdx].x());
820 sphere[1] = float(instancePoints[largestRadiusIdx].y());
821 sphere[2] = float(instancePoints[largestRadiusIdx].z());
822 sphere[3] = largestRadius;
824 spheres.push_back(sphere);
825 instanceMask[largestRadiusIdx] = 1;
828 sphere, instancePoints, instanceRadius, instanceMask, overlapping);
831 largestRadius = op.
radius();
832 largestRadiusIdx = op.
index();
839 template<
typename Gr
idT>
841 : mIsInitialized(false)
842 , mLeafBoundingSpheres(0)
843 , mNodeBoundingSpheres(0)
846 , mSurfacePointList()
854 template<
typename Gr
idT>
858 initialize<GridT, util::NullInterrupter>(grid, isovalue, NULL);
862 template<
typename Gr
idT>
863 template<
typename InterrupterT>
866 const GridT& grid,
float isovalue, InterrupterT* interrupter)
868 mIsInitialized =
false;
872 typedef typename GridT::ValueType ValueT;
880 LeafManagerT leafs(tree);
882 signDataOp(tree, leafs, ValueT(isovalue));
886 mSignTreePt = signDataOp.
signTree();
887 mIdxTreePt = signDataOp.
idxTree();
890 if (interrupter && interrupter->wasInterrupted())
return;
892 Int16LeafManagerT signLeafs(*mSignTreePt);
894 std::vector<size_t> regions(signLeafs.leafCount(), 0);
898 for (
size_t tmp = 0, n = 0, N = regions.size(); n < N; ++n) {
900 regions[n] = mPointListSize;
901 mPointListSize += tmp;
904 if (mPointListSize == 0)
return;
906 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
909 pointOp(signLeafs, tree, *mIdxTreePt, mSurfacePointList, regions, transform, isovalue);
913 mIdxTreePt->topologyUnion(*mSignTreePt);
916 if (interrupter && interrupter->wasInterrupted())
return;
919 CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
924 dim[0] = std::abs(dim[0]);
925 dim[1] = std::abs(dim[1]);
926 dim[2] = std::abs(dim[2]);
929 mMaxRadiusSqr *= 0.51f;
930 mMaxRadiusSqr *= mMaxRadiusSqr;
933 IntLeafManagerT idxLeafs(*mIdxTreePt);
936 typedef typename IntTreeT::RootNodeType IntRootNodeT;
937 typedef typename IntRootNodeT::NodeChainType IntNodeChainT;
938 BOOST_STATIC_ASSERT(boost::mpl::size<IntNodeChainT>::value > 1);
939 typedef typename boost::mpl::at<IntNodeChainT, boost::mpl::int_<1> >::type IntInternalNodeT;
942 typename IntTreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
943 nIt.setMinDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
944 nIt.setMaxDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
946 std::vector<const IntInternalNodeT*> internalNodes;
948 const IntInternalNodeT* node = NULL;
951 if (node) internalNodes.push_back(node);
954 std::vector<IndexRange>().swap(mLeafRanges);
955 mLeafRanges.resize(internalNodes.size());
957 std::vector<const IntLeafT*>().swap(mLeafNodes);
958 mLeafNodes.reserve(idxLeafs.leafCount());
960 typename IntInternalNodeT::ChildOnCIter leafIt;
962 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
964 mLeafRanges[n].first = mLeafNodes.size();
966 size_t leafCount = 0;
967 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
968 mLeafNodes.push_back(&(*leafIt));
972 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
974 mLeafRanges[n].second = mLeafNodes.size();
977 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
978 mLeafBoundingSpheres.resize(mLeafNodes.size());
981 mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
985 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
986 mNodeBoundingSpheres.resize(internalNodes.size());
988 internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
990 mIsInitialized =
true;
994 template<
typename Gr
idT>
997 std::vector<float>& distances,
bool transformPoints)
999 if (!mIsInitialized)
return false;
1002 distances.resize(points.size(), mMaxRadiusSqr);
1005 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1006 mMaxNodeLeafs, transformPoints);
1014 template<
typename Gr
idT>
1018 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
1022 template<
typename Gr
idT>
1025 std::vector<float>& distances)
1027 return search(points, distances,
true);
1035 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
uint64_t Index64
Definition: Types.h:58
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:721
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Definition: Exceptions.h:39
Vec4< float > Vec4s
Definition: Vec4.h:545
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:97
Vec3< double > Vec3d
Definition: Vec3.h:629
OPENVDB_IMPORT void initialize()
Global registration of basic types.
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:999
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:302
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Vec3< float > Vec3s
Definition: Vec3.h:628