35 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 36 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 44 #include <boost/mpl/at.hpp> 45 #include <boost/mpl/int.hpp> 46 #include <boost/scoped_array.hpp> 47 #include <tbb/blocked_range.h> 48 #include <tbb/parallel_for.h> 49 #include <tbb/parallel_reduce.h> 84 template<
typename Gr
idT,
typename InterrupterT = util::NullInterrupter>
88 std::vector<openvdb::Vec4s>& spheres,
90 bool overlapping =
false,
91 float minRadius = 1.0,
94 int instanceCount = 10000,
95 InterrupterT* interrupter =
nullptr);
104 template<
typename Gr
idT>
108 using Ptr = std::unique_ptr<ClosestSurfacePoint>;
109 using TreeT =
typename GridT::TreeType;
110 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
111 using Index32TreeT =
typename TreeT::template ValueConverter<Index32>::Type;
112 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
125 template<
typename InterrupterT = util::NullInterrupter>
126 static inline Ptr create(
const GridT& grid,
float isovalue = 0.0,
127 InterrupterT* interrupter =
nullptr);
132 inline bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
137 inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
145 using Index32LeafT =
typename Index32TreeT::LeafNodeType;
146 using IndexRange = std::pair<size_t, size_t>;
148 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
149 std::vector<IndexRange> mLeafRanges;
150 std::vector<const Index32LeafT*> mLeafNodes;
152 size_t mPointListSize = 0, mMaxNodeLeafs = 0;
153 typename Index32TreeT::Ptr mIdxTreePt;
154 typename Int16TreeT::Ptr mSignTreePt;
157 template<
typename InterrupterT = util::NullInterrupter>
158 inline bool initialize(
const GridT&,
float isovalue, InterrupterT*);
159 inline bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
168 namespace v2s_internal {
179 mPoints.push_back(pos);
182 std::vector<Vec3R>& mPoints;
186 template<
typename Index32LeafT>
191 LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
192 const std::vector<const Index32LeafT*>& leafNodes,
196 void run(
bool threaded =
true);
199 void operator()(
const tbb::blocked_range<size_t>&)
const;
202 std::vector<Vec4R>& mLeafBoundingSpheres;
203 const std::vector<const Index32LeafT*>& mLeafNodes;
208 template<
typename Index32LeafT>
210 std::vector<Vec4R>& leafBoundingSpheres,
211 const std::vector<const Index32LeafT*>& leafNodes,
214 : mLeafBoundingSpheres(leafBoundingSpheres)
215 , mLeafNodes(leafNodes)
216 , mTransform(transform)
217 , mSurfacePointList(surfacePointList)
221 template<
typename Index32LeafT>
226 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
228 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
232 template<
typename Index32LeafT>
236 typename Index32LeafT::ValueOnCIter iter;
239 for (
size_t n = range.begin(); n != range.end(); ++n) {
245 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
246 avg += mSurfacePointList[iter.getValue()];
249 if (count > 1) avg *= float(1.0 /
double(count));
252 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
253 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
254 if (tmpDist > maxDist) maxDist = tmpDist;
257 Vec4R& sphere = mLeafBoundingSpheres[n];
261 sphere[3] = std::sqrt(maxDist);
271 NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
272 const std::vector<IndexRange>& leafRanges,
273 const std::vector<Vec4R>& leafBoundingSpheres);
275 inline void run(
bool threaded =
true);
277 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
280 std::vector<Vec4R>& mNodeBoundingSpheres;
281 const std::vector<IndexRange>& mLeafRanges;
282 const std::vector<Vec4R>& mLeafBoundingSpheres;
287 const std::vector<IndexRange>& leafRanges,
288 const std::vector<Vec4R>& leafBoundingSpheres)
289 : mNodeBoundingSpheres(nodeBoundingSpheres)
290 , mLeafRanges(leafRanges)
291 , mLeafBoundingSpheres(leafBoundingSpheres)
299 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
301 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
310 for (
size_t n = range.begin(); n != range.end(); ++n) {
316 int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
318 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
319 avg[0] += mLeafBoundingSpheres[i][0];
320 avg[1] += mLeafBoundingSpheres[i][1];
321 avg[2] += mLeafBoundingSpheres[i][2];
324 if (count > 1) avg *= float(1.0 /
double(count));
327 double maxDist = 0.0;
329 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
330 pos[0] = mLeafBoundingSpheres[i][0];
331 pos[1] = mLeafBoundingSpheres[i][1];
332 pos[2] = mLeafBoundingSpheres[i][2];
334 double tmpDist = (pos - avg).length() + mLeafBoundingSpheres[i][3];
335 if (tmpDist > maxDist) maxDist = tmpDist;
338 Vec4R& sphere = mNodeBoundingSpheres[n];
351 template<
typename Index32LeafT>
358 std::vector<Vec3R>& instancePoints,
359 std::vector<float>& instanceDistances,
361 const std::vector<const Index32LeafT*>& leafNodes,
362 const std::vector<IndexRange>& leafRanges,
363 const std::vector<Vec4R>& leafBoundingSpheres,
364 const std::vector<Vec4R>& nodeBoundingSpheres,
366 bool transformPoints =
false);
369 void run(
bool threaded =
true);
372 void operator()(
const tbb::blocked_range<size_t>&)
const;
376 void evalLeaf(
size_t index,
const Index32LeafT& leaf)
const;
377 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
380 std::vector<Vec3R>& mInstancePoints;
381 std::vector<float>& mInstanceDistances;
385 const std::vector<const Index32LeafT*>& mLeafNodes;
386 const std::vector<IndexRange>& mLeafRanges;
387 const std::vector<Vec4R>& mLeafBoundingSpheres;
388 const std::vector<Vec4R>& mNodeBoundingSpheres;
390 std::vector<float> mLeafDistances, mNodeDistances;
392 const bool mTransformPoints;
393 size_t mClosestPointIndex;
397 template<
typename Index32LeafT>
399 std::vector<Vec3R>& instancePoints,
400 std::vector<float>& instanceDistances,
402 const std::vector<const Index32LeafT*>& leafNodes,
403 const std::vector<IndexRange>& leafRanges,
404 const std::vector<Vec4R>& leafBoundingSpheres,
405 const std::vector<Vec4R>& nodeBoundingSpheres,
407 bool transformPoints)
408 : mInstancePoints(instancePoints)
409 , mInstanceDistances(instanceDistances)
410 , mSurfacePointList(surfacePointList)
411 , mLeafNodes(leafNodes)
412 , mLeafRanges(leafRanges)
413 , mLeafBoundingSpheres(leafBoundingSpheres)
414 , mNodeBoundingSpheres(nodeBoundingSpheres)
415 , mLeafDistances(maxNodeLeafs, 0.0)
416 , mNodeDistances(leafRanges.size(), 0.0)
417 , mTransformPoints(transformPoints)
418 , mClosestPointIndex(0)
423 template<
typename Index32LeafT>
428 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
430 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
434 template<
typename Index32LeafT>
438 typename Index32LeafT::ValueOnCIter iter;
439 const Vec3s center = mInstancePoints[index];
440 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
442 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
444 const Vec3s& point = mSurfacePointList[iter.getValue()];
445 float tmpDist = (point - center).lengthSqr();
447 if (tmpDist < mInstanceDistances[index]) {
448 mInstanceDistances[index] = tmpDist;
449 closestPointIndex = iter.getValue();
455 template<
typename Index32LeafT>
459 if (nodeIndex >= mLeafRanges.size())
return;
461 const Vec3R& pos = mInstancePoints[pointIndex];
462 float minDist = mInstanceDistances[pointIndex];
463 size_t minDistIdx = 0;
465 bool updatedDist =
false;
467 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
468 i < mLeafRanges[nodeIndex].second; ++i, ++n)
470 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
472 center[0] = mLeafBoundingSpheres[i][0];
473 center[1] = mLeafBoundingSpheres[i][1];
474 center[2] = mLeafBoundingSpheres[i][2];
475 const auto radius = mLeafBoundingSpheres[i][3];
477 distToLeaf = float(
std::max(0.0, (pos - center).length() - radius));
479 if (distToLeaf < minDist) {
480 minDist = distToLeaf;
486 if (!updatedDist)
return;
488 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
490 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
491 i < mLeafRanges[nodeIndex].second; ++i, ++n)
493 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
494 evalLeaf(pointIndex, *mLeafNodes[i]);
500 template<
typename Index32LeafT>
505 for (
size_t n = range.begin(); n != range.end(); ++n) {
507 const Vec3R& pos = mInstancePoints[n];
508 float minDist = mInstanceDistances[n];
509 size_t minDistIdx = 0;
511 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
512 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
514 center[0] = mNodeBoundingSpheres[i][0];
515 center[1] = mNodeBoundingSpheres[i][1];
516 center[2] = mNodeBoundingSpheres[i][2];
517 const auto radius = mNodeBoundingSpheres[i][3];
519 distToNode = float(
std::max(0.0, (pos - center).length() - radius));
521 if (distToNode < minDist) {
522 minDist = distToNode;
527 evalNode(n, minDistIdx);
529 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
530 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
535 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
537 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
547 const std::vector<Vec3R>& points,
548 std::vector<float>& distances,
549 std::vector<unsigned char>& mask,
553 int index()
const {
return mIndex; }
555 inline void run(
bool threaded =
true);
559 inline void operator()(
const tbb::blocked_range<size_t>& range);
562 if (rhs.mRadius > mRadius) {
563 mRadius = rhs.mRadius;
569 const Vec4s& mSphere;
570 const std::vector<Vec3R>& mPoints;
571 std::vector<float>& mDistances;
572 std::vector<unsigned char>& mMask;
581 const std::vector<Vec3R>& points,
582 std::vector<float>& distances,
583 std::vector<unsigned char>& mask,
587 , mDistances(distances)
589 , mOverlapping(overlapping)
597 : mSphere(rhs.mSphere)
598 , mPoints(rhs.mPoints)
599 , mDistances(rhs.mDistances)
601 , mOverlapping(rhs.mOverlapping)
602 , mRadius(rhs.mRadius)
611 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
613 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
621 for (
size_t n = range.begin(); n != range.end(); ++n) {
622 if (mMask[n])
continue;
624 pos.x() = float(mPoints[n].x()) - mSphere[0];
625 pos.y() = float(mPoints[n].y()) - mSphere[1];
626 pos.z() = float(mPoints[n].z()) - mSphere[2];
628 float dist = pos.length();
630 if (dist < mSphere[3]) {
636 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
639 if (mDistances[n] > mRadius) {
640 mRadius = mDistances[n];
653 template<
typename Gr
idT,
typename InterrupterT>
657 std::vector<openvdb::Vec4s>& spheres,
664 InterrupterT* interrupter)
667 spheres.reserve(maxSphereCount);
669 const bool addNBPoints = grid.activeVoxelCount() < 10000;
670 int instances =
std::max(instanceCount, maxSphereCount);
672 using TreeT =
typename GridT::TreeType;
673 using ValueT =
typename GridT::ValueType;
675 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
676 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
678 using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
679 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>;
682 const TreeT& tree = grid.tree();
685 std::vector<Vec3R> instancePoints;
708 interiorMaskPtr->
tree().topologyUnion(tree);
711 if (interrupter && interrupter->wasInterrupted())
return;
716 instancePoints.reserve(instances);
720 ptnAcc,
Index64(addNBPoints ? (instances / 2) : instances), mtRand, 1.0, interrupter);
722 scatter(*interiorMaskPtr);
725 if (interrupter && interrupter->wasInterrupted())
return;
731 if (instancePoints.size() < size_t(instances)) {
732 const Int16TreeT& signTree = csp->signTree();
733 for (
auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
734 for (
auto it = leafIt->cbeginValueOn(); it; ++it) {
735 const int flags = int(it.getValue());
739 instancePoints.push_back(transform.
indexToWorld(it.getCoord()));
741 if (instancePoints.size() == size_t(instances))
break;
743 if (instancePoints.size() == size_t(instances))
break;
747 if (interrupter && interrupter->wasInterrupted())
return;
749 std::vector<float> instanceRadius;
750 if (!csp->search(instancePoints, instanceRadius))
return;
752 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
753 float largestRadius = 0.0;
754 int largestRadiusIdx = 0;
756 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
757 if (instanceRadius[n] > largestRadius) {
758 largestRadius = instanceRadius[n];
759 largestRadiusIdx = int(n);
765 minRadius = float(minRadius * transform.
voxelSize()[0]);
766 maxRadius = float(maxRadius * transform.
voxelSize()[0]);
768 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
770 if (interrupter && interrupter->wasInterrupted())
return;
772 largestRadius =
std::min(maxRadius, largestRadius);
774 if (s != 0 && largestRadius < minRadius)
break;
776 sphere[0] = float(instancePoints[largestRadiusIdx].x());
777 sphere[1] = float(instancePoints[largestRadiusIdx].y());
778 sphere[2] = float(instancePoints[largestRadiusIdx].z());
779 sphere[3] = largestRadius;
781 spheres.push_back(sphere);
782 instanceMask[largestRadiusIdx] = 1;
785 sphere, instancePoints, instanceRadius, instanceMask, overlapping);
788 largestRadius = op.
radius();
789 largestRadiusIdx = op.
index();
797 template<
typename Gr
idT>
798 template<
typename InterrupterT>
803 if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
808 template<
typename Gr
idT>
809 template<
typename InterrupterT>
812 const GridT& grid,
float isovalue, InterrupterT* interrupter)
815 using ValueT =
typename GridT::ValueType;
830 *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
832 if (interrupter && interrupter->wasInterrupted())
return false;
836 using Int16LeafNodeType =
typename Int16TreeT::LeafNodeType;
837 using Index32LeafNodeType =
typename Index32TreeT::LeafNodeType;
839 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
840 mSignTreePt->getNodes(signFlagsLeafNodes);
842 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
844 boost::scoped_array<Index32> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
846 tbb::parallel_for(auxiliaryLeafNodeRange,
848 (signFlagsLeafNodes, leafNodeOffsets));
852 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
853 const Index32 tmp = leafNodeOffsets[n];
858 mPointListSize = size_t(pointCount);
859 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
863 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
864 mIdxTreePt->getNodes(pointIndexLeafNodes);
867 mSurfacePointList.get(), tree, pointIndexLeafNodes,
868 signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
871 if (interrupter && interrupter->wasInterrupted())
return false;
873 Index32LeafManagerT idxLeafs(*mIdxTreePt);
875 using Index32RootNodeT =
typename Index32TreeT::RootNodeType;
876 using Index32NodeChainT =
typename Index32RootNodeT::NodeChainType;
877 static_assert(boost::mpl::size<Index32NodeChainT>::value > 1,
878 "expected tree depth greater than one");
879 using Index32InternalNodeT =
880 typename boost::mpl::at<Index32NodeChainT, boost::mpl::int_<1> >::type;
882 typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
883 nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
884 nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
886 std::vector<const Index32InternalNodeT*> internalNodes;
888 const Index32InternalNodeT* node =
nullptr;
891 if (node) internalNodes.push_back(node);
894 std::vector<IndexRange>().swap(mLeafRanges);
895 mLeafRanges.resize(internalNodes.size());
897 std::vector<const Index32LeafT*>().swap(mLeafNodes);
898 mLeafNodes.reserve(idxLeafs.leafCount());
900 typename Index32InternalNodeT::ChildOnCIter leafIt;
902 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
904 mLeafRanges[n].first = mLeafNodes.size();
906 size_t leafCount = 0;
907 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
908 mLeafNodes.push_back(&(*leafIt));
912 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
914 mLeafRanges[n].second = mLeafNodes.size();
917 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
918 mLeafBoundingSpheres.resize(mLeafNodes.size());
921 mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
925 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
926 mNodeBoundingSpheres.resize(internalNodes.size());
934 template<
typename Gr
idT>
937 std::vector<float>& distances,
bool transformPoints)
940 distances.resize(points.size(), std::numeric_limits<float>::infinity());
943 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
944 mMaxNodeLeafs, transformPoints);
952 template<
typename Gr
idT>
956 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
960 template<
typename Gr
idT>
963 std::vector<float>& distances)
965 return search(points, distances,
true);
972 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
uint32_t Index32
Definition: Types.h:59
Index64 pointCount(const PointDataTreeT &tree, const bool inCoreOnly=false)
Total points in the PointDataTree.
Definition: PointCount.h:241
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:343
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:230
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1112
Tolerance for floating-point comparison.
Definition: Math.h:117
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
OPENVDB_IMPORT void initialize()
Global registration of basic types.
Extract polygonal surfaces from scalar volumes.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
Vec3< double > Vec3d
Definition: Vec3.h:679
uint64_t Index64
Definition: Types.h:60
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:40
Vec3< float > Vec3s
Definition: Vec3.h:678
static T value()
Definition: Math.h:117
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:797
Vec4< float > Vec4s
Definition: Vec4.h:586
SharedPtr< Grid > Ptr
Definition: Grid.h:502
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:55
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...