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>
147 template<
typename InterrupterT>
148 void initialize(
const GridT& grid,
float isovalue = 0.0, InterrupterT* interrupter = NULL);
153 void initialize(
const GridT& grid,
float isovalue = 0.0);
163 bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
173 bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
176 typedef typename GridT::TreeType TreeT;
177 typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
178 typedef typename IntTreeT::LeafNodeType IntLeafT;
179 typedef std::pair<size_t, size_t> IndexRange;
182 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
183 std::vector<IndexRange> mLeafRanges;
184 std::vector<const IntLeafT*> mLeafNodes;
186 size_t mPointListSize, mMaxNodeLeafs;
188 typename IntTreeT::Ptr mIdxTreePt;
190 bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
213 mPoints.push_back(pos);
216 std::vector<Vec3R>& mPoints;
220 template<
typename IntLeafT>
225 LeafBS(std::vector<Vec4R>& leafBoundingSpheres,
226 const std::vector<const IntLeafT*>& leafNodes,
230 void run(
bool threaded =
true);
233 void operator()(
const tbb::blocked_range<size_t>&)
const;
236 std::vector<Vec4R>& mLeafBoundingSpheres;
237 const std::vector<const IntLeafT*>& mLeafNodes;
242 template<
typename IntLeafT>
244 std::vector<Vec4R>& leafBoundingSpheres,
245 const std::vector<const IntLeafT*>& leafNodes,
248 : mLeafBoundingSpheres(leafBoundingSpheres)
249 , mLeafNodes(leafNodes)
250 , mTransform(transform)
251 , mSurfacePointList(surfacePointList)
255 template<
typename IntLeafT>
260 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
262 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
266 template<
typename IntLeafT>
270 typename IntLeafT::ValueOnCIter iter;
273 for (
size_t n = range.begin(); n != range.end(); ++n) {
280 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
281 avg += mSurfacePointList[iter.getValue()];
285 if (count > 1) avg *= float(1.0 /
double(count));
287 float maxDist = mTransform.voxelSize()[0];
290 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
291 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
292 if (tmpDist > maxDist) maxDist = tmpDist;
295 Vec4R& sphere = mLeafBoundingSpheres[n];
310 NodeBS(std::vector<Vec4R>& nodeBoundingSpheres,
311 const std::vector<IndexRange>& leafRanges,
312 const std::vector<Vec4R>& leafBoundingSpheres);
314 void run(
bool threaded =
true);
317 void operator()(
const tbb::blocked_range<size_t>&)
const;
320 std::vector<Vec4R>& mNodeBoundingSpheres;
321 const std::vector<IndexRange>& mLeafRanges;
322 const std::vector<Vec4R>& mLeafBoundingSpheres;
327 const std::vector<IndexRange>& leafRanges,
328 const std::vector<Vec4R>& leafBoundingSpheres)
329 : mNodeBoundingSpheres(nodeBoundingSpheres)
330 , mLeafRanges(leafRanges)
331 , mLeafBoundingSpheres(leafBoundingSpheres)
339 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
341 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
350 for (
size_t n = range.begin(); n != range.end(); ++n) {
356 int count = mLeafRanges[n].second - mLeafRanges[n].first;
358 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
359 avg[0] += mLeafBoundingSpheres[i][0];
360 avg[1] += mLeafBoundingSpheres[i][1];
361 avg[2] += mLeafBoundingSpheres[i][2];
364 if (count > 1) avg *= float(1.0 /
double(count));
369 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
370 pos[0] = mLeafBoundingSpheres[i][0];
371 pos[1] = mLeafBoundingSpheres[i][1];
372 pos[2] = mLeafBoundingSpheres[i][2];
374 float tmpDist = (pos - avg).lengthSqr() + mLeafBoundingSpheres[i][3];
375 if (tmpDist > maxDist) maxDist = tmpDist;
378 Vec4R& sphere = mNodeBoundingSpheres[n];
392 template<
typename IntLeafT>
399 std::vector<Vec3R>& instancePoints,
400 std::vector<float>& instanceDistances,
402 const std::vector<const IntLeafT*>& leafNodes,
403 const std::vector<IndexRange>& leafRanges,
404 const std::vector<Vec4R>& leafBoundingSpheres,
405 const std::vector<Vec4R>& nodeBoundingSpheres,
407 bool transformPoints =
false);
410 void run(
bool threaded =
true);
413 void operator()(
const tbb::blocked_range<size_t>&)
const;
417 void evalLeaf(
size_t index,
const IntLeafT& leaf)
const;
418 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
421 std::vector<Vec3R>& mInstancePoints;
422 std::vector<float>& mInstanceDistances;
426 const std::vector<const IntLeafT*>& mLeafNodes;
427 const std::vector<IndexRange>& mLeafRanges;
428 const std::vector<Vec4R>& mLeafBoundingSpheres;
429 const std::vector<Vec4R>& mNodeBoundingSpheres;
431 std::vector<float> mLeafDistances, mNodeDistances;
433 const bool mTransformPoints;
434 size_t mClosestPointIndex;
438 template<
typename IntLeafT>
440 std::vector<Vec3R>& instancePoints,
441 std::vector<float>& instanceDistances,
443 const std::vector<const IntLeafT*>& leafNodes,
444 const std::vector<IndexRange>& leafRanges,
445 const std::vector<Vec4R>& leafBoundingSpheres,
446 const std::vector<Vec4R>& nodeBoundingSpheres,
448 bool transformPoints)
449 : mInstancePoints(instancePoints)
450 , mInstanceDistances(instanceDistances)
451 , mSurfacePointList(surfacePointList)
452 , mLeafNodes(leafNodes)
453 , mLeafRanges(leafRanges)
454 , mLeafBoundingSpheres(leafBoundingSpheres)
455 , mNodeBoundingSpheres(nodeBoundingSpheres)
456 , mLeafDistances(maxNodeLeafs, 0.0)
457 , mNodeDistances(leafRanges.size(), 0.0)
458 , mTransformPoints(transformPoints)
459 , mClosestPointIndex(0)
464 template<
typename IntLeafT>
469 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
471 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
475 template<
typename IntLeafT>
479 typename IntLeafT::ValueOnCIter iter;
480 const Vec3s center = mInstancePoints[index];
481 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
483 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
485 const Vec3s& point = mSurfacePointList[iter.getValue()];
486 float tmpDist = (point - center).lengthSqr();
488 if (tmpDist < mInstanceDistances[index]) {
489 mInstanceDistances[index] = tmpDist;
490 closestPointIndex = iter.getValue();
496 template<
typename IntLeafT>
498 ClosestPointDist<IntLeafT>::evalNode(
size_t pointIndex,
size_t nodeIndex)
const
500 const Vec3R& pos = mInstancePoints[pointIndex];
501 float minDist = mInstanceDistances[pointIndex];
502 size_t minDistIdx = 0;
504 bool updatedDist =
false;
507 for (
size_t i = 0, I = mLeafDistances.size(); i < I; ++i) {
508 float& distToLeaf =
const_cast<float&
>(mLeafDistances[i]);
512 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0; i < mLeafRanges[nodeIndex].second; ++i, ++n) {
514 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
516 center[0] = mLeafBoundingSpheres[i][0];
517 center[1] = mLeafBoundingSpheres[i][1];
518 center[2] = mLeafBoundingSpheres[i][2];
520 distToLeaf = (pos - center).lengthSqr() - mLeafBoundingSpheres[i][3];
522 if (distToLeaf < minDist) {
523 minDist = distToLeaf;
529 if (!updatedDist)
return;
531 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
533 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0; i < mLeafRanges[nodeIndex].second; ++i, ++n) {
534 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
535 evalLeaf(pointIndex, *mLeafNodes[i]);
541 template<
typename IntLeafT>
546 for (
size_t n = range.begin(); n != range.end(); ++n) {
548 const Vec3R& pos = mInstancePoints[n];
549 float minDist = mInstanceDistances[n];
550 size_t minDistIdx = 0;
553 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
554 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
556 center[0] = mNodeBoundingSpheres[i][0];
557 center[1] = mNodeBoundingSpheres[i][1];
558 center[2] = mNodeBoundingSpheres[i][2];
560 distToNode = (pos - center).lengthSqr() - mNodeBoundingSpheres[i][3];
562 if (distToNode < minDist) {
563 minDist = distToNode;
568 evalNode(n, minDistIdx);
570 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
571 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
576 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
578 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
588 const std::vector<Vec3R>& points,
589 std::vector<float>& distances,
590 std::vector<unsigned char>& mask,
594 int index()
const {
return mIndex; };
596 void run(
bool threaded =
true);
600 void operator()(
const tbb::blocked_range<size_t>& range);
603 if (rhs.mRadius > mRadius) {
604 mRadius = rhs.mRadius;
611 const Vec4s& mSphere;
612 const std::vector<Vec3R>& mPoints;
614 std::vector<float>& mDistances;
615 std::vector<unsigned char>& mMask;
625 const std::vector<Vec3R>& points,
626 std::vector<float>& distances,
627 std::vector<unsigned char>& mask,
631 , mDistances(distances)
633 , mOverlapping(overlapping)
641 : mSphere(rhs.mSphere)
642 , mPoints(rhs.mPoints)
643 , mDistances(rhs.mDistances)
645 , mOverlapping(rhs.mOverlapping)
646 , mRadius(rhs.mRadius)
655 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
657 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
665 for (
size_t n = range.begin(); n != range.end(); ++n) {
666 if (mMask[n])
continue;
668 pos.
x() = float(mPoints[n].x()) - mSphere[0];
669 pos.y() = float(mPoints[n].y()) - mSphere[1];
670 pos.z() = float(mPoints[n].z()) - mSphere[2];
672 float dist = pos.length();
674 if (dist < mSphere[3]) {
680 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
683 if (mDistances[n] > mRadius) {
684 mRadius = mDistances[n];
697 template<
typename Gr
idT,
typename InterrupterT>
701 std::vector<openvdb::Vec4s>& spheres,
708 InterrupterT* interrupter)
711 spheres.reserve(maxSphereCount);
713 int instances =
std::max(instanceCount, maxSphereCount);
715 typedef typename GridT::TreeType TreeT;
716 typedef typename GridT::ValueType ValueT;
718 typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
719 typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
720 typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
727 typedef boost::mt11213b RandGen;
730 const TreeT& tree = grid.tree();
733 std::vector<Vec3R> instancePoints;
743 interiorMaskPtr->
tree().topologyUnion(tree);
746 if (interrupter && interrupter->wasInterrupted())
return;
750 instancePoints.reserve(instances);
754 scatter(ptnAcc, instances, mtRand, interrupter);
756 scatter(*interiorMaskPtr);
759 if (interrupter && interrupter->wasInterrupted())
return;
761 std::vector<float> instanceRadius;
766 if (interrupter && interrupter->wasInterrupted())
return;
768 if (!csp.
search(instancePoints, instanceRadius))
return;
770 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
771 float largestRadius = 0.0;
772 int largestRadiusIdx = 0;
774 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
775 if (instanceRadius[n] > largestRadius) {
776 largestRadius = instanceRadius[n];
777 largestRadiusIdx = n;
786 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
788 if (interrupter && interrupter->wasInterrupted())
return;
790 largestRadius =
std::min(maxRadius, largestRadius);
792 if (s != 0 && largestRadius < minRadius)
break;
794 sphere[0] = float(instancePoints[largestRadiusIdx].x());
795 sphere[1] = float(instancePoints[largestRadiusIdx].y());
796 sphere[2] = float(instancePoints[largestRadiusIdx].z());
797 sphere[3] = largestRadius;
799 spheres.push_back(sphere);
800 instanceMask[largestRadiusIdx] = 1;
805 largestRadius = op.
radius();
806 largestRadiusIdx = op.
index();
813 template<
typename Gr
idT>
815 : mIsInitialized(false)
816 , mLeafBoundingSpheres(0)
817 , mNodeBoundingSpheres(0)
820 , mSurfacePointList()
828 template<
typename Gr
idT>
832 initialize<GridT, util::NullInterrupter>(grid, isovalue, NULL);
836 template<
typename Gr
idT>
837 template<
typename InterrupterT>
840 const GridT& grid,
float isovalue, InterrupterT* interrupter)
842 mIsInitialized =
false;
843 typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
847 typedef typename GridT::ValueType ValueT;
849 const TreeT& tree = grid.
tree();
854 typename Int16TreeT::Ptr signTreePt;
857 LeafManagerT leafs(tree);
859 signDataOp(tree, leafs, ValueT(isovalue));
864 mIdxTreePt = signDataOp.
idxTree();
867 if (interrupter && interrupter->wasInterrupted())
return;
869 Int16LeafManagerT signLeafs(*signTreePt);
871 std::vector<size_t> regions(signLeafs.leafCount(), 0);
875 for (
size_t tmp = 0, n = 0, N = regions.size(); n < N; ++n) {
877 regions[n] = mPointListSize;
878 mPointListSize += tmp;
881 if (mPointListSize == 0)
return;
883 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
886 pointOp(signLeafs, tree, *mIdxTreePt, mSurfacePointList, regions, transform, isovalue);
890 mIdxTreePt->topologyUnion(*signTreePt);
893 if (interrupter && interrupter->wasInterrupted())
return;
896 CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
901 dim[0] = std::abs(dim[0]);
902 dim[1] = std::abs(dim[1]);
903 dim[2] = std::abs(dim[2]);
906 mMaxRadiusSqr *= 0.51;
907 mMaxRadiusSqr *= mMaxRadiusSqr;
910 IntLeafManagerT idxLeafs(*mIdxTreePt);
913 typedef typename IntTreeT::RootNodeType IntRootNodeT;
914 typedef typename IntRootNodeT::NodeChainType IntNodeChainT;
915 BOOST_STATIC_ASSERT(boost::mpl::size<IntNodeChainT>::value > 1);
916 typedef typename boost::mpl::at<IntNodeChainT, boost::mpl::int_<1> >::type IntInternalNodeT;
920 typename IntTreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
921 nIt.setMinDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
922 nIt.setMaxDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
924 std::vector<const IntInternalNodeT*> internalNodes;
926 const IntInternalNodeT* node = NULL;
929 if (node) internalNodes.push_back(node);
932 std::vector<IndexRange>().swap(mLeafRanges);
933 mLeafRanges.resize(internalNodes.size());
935 std::vector<const IntLeafT*>().swap(mLeafNodes);
936 mLeafNodes.reserve(idxLeafs.leafCount());
938 typename IntInternalNodeT::ChildOnCIter leafIt;
940 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
942 mLeafRanges[n].first = mLeafNodes.size();
944 size_t leafCount = 0;
945 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
946 mLeafNodes.push_back(&(*leafIt));
950 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
952 mLeafRanges[n].second = mLeafNodes.size();
955 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
956 mLeafBoundingSpheres.resize(mLeafNodes.size());
962 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
963 mNodeBoundingSpheres.resize(internalNodes.size());
965 internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
967 mIsInitialized =
true;
971 template<
typename Gr
idT>
974 std::vector<float>& distances,
bool transformPoints)
976 if (!mIsInitialized)
return false;
979 distances.resize(points.size(), mMaxRadiusSqr);
982 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
983 mMaxNodeLeafs, transformPoints);
991 template<
typename Gr
idT>
995 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
999 template<
typename Gr
idT>
1003 return search(points, distances,
true);
1011 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:696
TreeType & tree()
Return the tree associated with this manager.
Definition: LeafManager.h:298
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
Vec3< float > Vec3s
Definition: Vec3.h:604
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:461
#define OPENVDB_VERSION_NAME
Definition: version.h:45
const MaskGridType * mMask
Definition: GridOperators.h:386
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
OPENVDB_IMPORT void initialize()
Global registration of basic types.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:94
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:964
Vec4< float > Vec4s
Definition: Vec4.h:523