35 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 36 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED 38 #include <openvdb/tree/LeafManager.h> 44 #include <boost/scoped_array.hpp> 45 #include <tbb/blocked_range.h> 46 #include <tbb/parallel_for.h> 47 #include <tbb/parallel_reduce.h> 78 template<
typename Gr
idT,
typename InterrupterT = util::NullInterrupter>
82 std::vector<openvdb::Vec4s>& spheres,
84 bool overlapping =
false,
85 float minRadius = 1.0,
88 int instanceCount = 10000,
89 InterrupterT* interrupter =
nullptr);
98 template<
typename Gr
idT>
102 using Ptr = std::unique_ptr<ClosestSurfacePoint>;
103 using TreeT =
typename GridT::TreeType;
104 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
105 using Index32TreeT =
typename TreeT::template ValueConverter<Index32>::Type;
106 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
119 template<
typename InterrupterT = util::NullInterrupter>
120 static inline Ptr create(
const GridT& grid,
float isovalue = 0.0,
121 InterrupterT* interrupter =
nullptr);
126 inline bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
131 inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
139 using Index32LeafT =
typename Index32TreeT::LeafNodeType;
140 using IndexRange = std::pair<size_t, size_t>;
142 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
143 std::vector<IndexRange> mLeafRanges;
144 std::vector<const Index32LeafT*> mLeafNodes;
146 size_t mPointListSize = 0, mMaxNodeLeafs = 0;
147 typename Index32TreeT::Ptr mIdxTreePt;
148 typename Int16TreeT::Ptr mSignTreePt;
151 template<
typename InterrupterT = util::NullInterrupter>
152 inline bool initialize(
const GridT&,
float isovalue, InterrupterT*);
153 inline bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
162 namespace v2s_internal {
173 mPoints.push_back(pos);
176 std::vector<Vec3R>& mPoints;
180 template<
typename Index32LeafT>
185 LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
186 const std::vector<const Index32LeafT*>& leafNodes,
190 void run(
bool threaded =
true);
193 void operator()(
const tbb::blocked_range<size_t>&)
const;
196 std::vector<Vec4R>& mLeafBoundingSpheres;
197 const std::vector<const Index32LeafT*>& mLeafNodes;
202 template<
typename Index32LeafT>
204 std::vector<Vec4R>& leafBoundingSpheres,
205 const std::vector<const Index32LeafT*>& leafNodes,
208 : mLeafBoundingSpheres(leafBoundingSpheres)
209 , mLeafNodes(leafNodes)
210 , mTransform(transform)
211 , mSurfacePointList(surfacePointList)
215 template<
typename Index32LeafT>
220 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
222 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
226 template<
typename Index32LeafT>
230 typename Index32LeafT::ValueOnCIter iter;
233 for (
size_t n = range.begin(); n != range.end(); ++n) {
239 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
240 avg += mSurfacePointList[iter.getValue()];
243 if (count > 1) avg *= float(1.0 /
double(count));
246 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
247 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
248 if (tmpDist > maxDist) maxDist = tmpDist;
251 Vec4R& sphere = mLeafBoundingSpheres[n];
255 sphere[3] = std::sqrt(maxDist);
265 NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
266 const std::vector<IndexRange>& leafRanges,
267 const std::vector<Vec4R>& leafBoundingSpheres);
269 inline void run(
bool threaded =
true);
271 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
274 std::vector<Vec4R>& mNodeBoundingSpheres;
275 const std::vector<IndexRange>& mLeafRanges;
276 const std::vector<Vec4R>& mLeafBoundingSpheres;
281 const std::vector<IndexRange>& leafRanges,
282 const std::vector<Vec4R>& leafBoundingSpheres)
283 : mNodeBoundingSpheres(nodeBoundingSpheres)
284 , mLeafRanges(leafRanges)
285 , mLeafBoundingSpheres(leafBoundingSpheres)
293 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
295 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
304 for (
size_t n = range.begin(); n != range.end(); ++n) {
310 int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
312 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
313 avg[0] += mLeafBoundingSpheres[i][0];
314 avg[1] += mLeafBoundingSpheres[i][1];
315 avg[2] += mLeafBoundingSpheres[i][2];
318 if (count > 1) avg *= float(1.0 /
double(count));
321 double maxDist = 0.0;
323 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
324 pos[0] = mLeafBoundingSpheres[i][0];
325 pos[1] = mLeafBoundingSpheres[i][1];
326 pos[2] = mLeafBoundingSpheres[i][2];
328 double tmpDist = (pos - avg).length() + mLeafBoundingSpheres[i][3];
329 if (tmpDist > maxDist) maxDist = tmpDist;
332 Vec4R& sphere = mNodeBoundingSpheres[n];
345 template<
typename Index32LeafT>
352 std::vector<Vec3R>& instancePoints,
353 std::vector<float>& instanceDistances,
355 const std::vector<const Index32LeafT*>& leafNodes,
356 const std::vector<IndexRange>& leafRanges,
357 const std::vector<Vec4R>& leafBoundingSpheres,
358 const std::vector<Vec4R>& nodeBoundingSpheres,
360 bool transformPoints =
false);
363 void run(
bool threaded =
true);
366 void operator()(
const tbb::blocked_range<size_t>&)
const;
370 void evalLeaf(
size_t index,
const Index32LeafT& leaf)
const;
371 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
374 std::vector<Vec3R>& mInstancePoints;
375 std::vector<float>& mInstanceDistances;
379 const std::vector<const Index32LeafT*>& mLeafNodes;
380 const std::vector<IndexRange>& mLeafRanges;
381 const std::vector<Vec4R>& mLeafBoundingSpheres;
382 const std::vector<Vec4R>& mNodeBoundingSpheres;
384 std::vector<float> mLeafDistances, mNodeDistances;
386 const bool mTransformPoints;
387 size_t mClosestPointIndex;
391 template<
typename Index32LeafT>
393 std::vector<Vec3R>& instancePoints,
394 std::vector<float>& instanceDistances,
396 const std::vector<const Index32LeafT*>& leafNodes,
397 const std::vector<IndexRange>& leafRanges,
398 const std::vector<Vec4R>& leafBoundingSpheres,
399 const std::vector<Vec4R>& nodeBoundingSpheres,
401 bool transformPoints)
402 : mInstancePoints(instancePoints)
403 , mInstanceDistances(instanceDistances)
404 , mSurfacePointList(surfacePointList)
405 , mLeafNodes(leafNodes)
406 , mLeafRanges(leafRanges)
407 , mLeafBoundingSpheres(leafBoundingSpheres)
408 , mNodeBoundingSpheres(nodeBoundingSpheres)
409 , mLeafDistances(maxNodeLeafs, 0.0)
410 , mNodeDistances(leafRanges.size(), 0.0)
411 , mTransformPoints(transformPoints)
412 , mClosestPointIndex(0)
417 template<
typename Index32LeafT>
422 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
424 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
428 template<
typename Index32LeafT>
432 typename Index32LeafT::ValueOnCIter iter;
433 const Vec3s center = mInstancePoints[index];
434 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
436 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
438 const Vec3s& point = mSurfacePointList[iter.getValue()];
439 float tmpDist = (point - center).lengthSqr();
441 if (tmpDist < mInstanceDistances[index]) {
442 mInstanceDistances[index] = tmpDist;
443 closestPointIndex = iter.getValue();
449 template<
typename Index32LeafT>
453 if (nodeIndex >= mLeafRanges.size())
return;
455 const Vec3R& pos = mInstancePoints[pointIndex];
456 float minDist = mInstanceDistances[pointIndex];
457 size_t minDistIdx = 0;
459 bool updatedDist =
false;
461 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
462 i < mLeafRanges[nodeIndex].second; ++i, ++n)
464 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
466 center[0] = mLeafBoundingSpheres[i][0];
467 center[1] = mLeafBoundingSpheres[i][1];
468 center[2] = mLeafBoundingSpheres[i][2];
469 const auto radius = mLeafBoundingSpheres[i][3];
471 distToLeaf = float(
std::max(0.0, (pos - center).length() - radius));
473 if (distToLeaf < minDist) {
474 minDist = distToLeaf;
480 if (!updatedDist)
return;
482 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
484 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
485 i < mLeafRanges[nodeIndex].second; ++i, ++n)
487 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
488 evalLeaf(pointIndex, *mLeafNodes[i]);
494 template<
typename Index32LeafT>
499 for (
size_t n = range.begin(); n != range.end(); ++n) {
501 const Vec3R& pos = mInstancePoints[n];
502 float minDist = mInstanceDistances[n];
503 size_t minDistIdx = 0;
505 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
506 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
508 center[0] = mNodeBoundingSpheres[i][0];
509 center[1] = mNodeBoundingSpheres[i][1];
510 center[2] = mNodeBoundingSpheres[i][2];
511 const auto radius = mNodeBoundingSpheres[i][3];
513 distToNode = float(
std::max(0.0, (pos - center).length() - radius));
515 if (distToNode < minDist) {
516 minDist = distToNode;
521 evalNode(n, minDistIdx);
523 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
524 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
529 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
531 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
541 const std::vector<Vec3R>& points,
542 std::vector<float>& distances,
543 std::vector<unsigned char>& mask,
547 int index()
const {
return mIndex; }
549 inline void run(
bool threaded =
true);
553 inline void operator()(
const tbb::blocked_range<size_t>& range);
556 if (rhs.mRadius > mRadius) {
557 mRadius = rhs.mRadius;
563 const Vec4s& mSphere;
564 const std::vector<Vec3R>& mPoints;
565 std::vector<float>& mDistances;
566 std::vector<unsigned char>& mMask;
575 const std::vector<Vec3R>& points,
576 std::vector<float>& distances,
577 std::vector<unsigned char>& mask,
581 , mDistances(distances)
583 , mOverlapping(overlapping)
591 : mSphere(rhs.mSphere)
592 , mPoints(rhs.mPoints)
593 , mDistances(rhs.mDistances)
595 , mOverlapping(rhs.mOverlapping)
596 , mRadius(rhs.mRadius)
605 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
607 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
615 for (
size_t n = range.begin(); n != range.end(); ++n) {
616 if (mMask[n])
continue;
618 pos.x() = float(mPoints[n].x()) - mSphere[0];
619 pos.y() = float(mPoints[n].y()) - mSphere[1];
620 pos.z() = float(mPoints[n].z()) - mSphere[2];
622 float dist = pos.length();
624 if (dist < mSphere[3]) {
630 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
633 if (mDistances[n] > mRadius) {
634 mRadius = mDistances[n];
647 template<
typename Gr
idT,
typename InterrupterT>
651 std::vector<openvdb::Vec4s>& spheres,
658 InterrupterT* interrupter)
661 spheres.reserve(maxSphereCount);
663 const bool addNBPoints = grid.activeVoxelCount() < 10000;
664 int instances =
std::max(instanceCount, maxSphereCount);
666 using TreeT =
typename GridT::TreeType;
667 using ValueT =
typename GridT::ValueType;
669 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
670 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
672 using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
673 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>;
676 const TreeT& tree = grid.tree();
679 std::vector<Vec3R> instancePoints;
702 interiorMaskPtr->
tree().topologyUnion(tree);
705 if (interrupter && interrupter->wasInterrupted())
return;
710 instancePoints.reserve(instances);
714 ptnAcc,
Index64(addNBPoints ? (instances / 2) : instances), mtRand, 1.0, interrupter);
716 scatter(*interiorMaskPtr);
719 if (interrupter && interrupter->wasInterrupted())
return;
725 if (instancePoints.size() < size_t(instances)) {
726 const Int16TreeT& signTree = csp->signTree();
727 for (
auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
728 for (
auto it = leafIt->cbeginValueOn(); it; ++it) {
729 const int flags = int(it.getValue());
733 instancePoints.push_back(transform.
indexToWorld(it.getCoord()));
735 if (instancePoints.size() == size_t(instances))
break;
737 if (instancePoints.size() == size_t(instances))
break;
741 if (interrupter && interrupter->wasInterrupted())
return;
743 std::vector<float> instanceRadius;
744 if (!csp->search(instancePoints, instanceRadius))
return;
746 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
747 float largestRadius = 0.0;
748 int largestRadiusIdx = 0;
750 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
751 if (instanceRadius[n] > largestRadius) {
752 largestRadius = instanceRadius[n];
753 largestRadiusIdx = int(n);
759 minRadius = float(minRadius * transform.
voxelSize()[0]);
760 maxRadius = float(maxRadius * transform.
voxelSize()[0]);
762 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
764 if (interrupter && interrupter->wasInterrupted())
return;
766 largestRadius =
std::min(maxRadius, largestRadius);
768 if (s != 0 && largestRadius < minRadius)
break;
770 sphere[0] = float(instancePoints[largestRadiusIdx].x());
771 sphere[1] = float(instancePoints[largestRadiusIdx].y());
772 sphere[2] = float(instancePoints[largestRadiusIdx].z());
773 sphere[3] = largestRadius;
775 spheres.push_back(sphere);
776 instanceMask[largestRadiusIdx] = 1;
779 sphere, instancePoints, instanceRadius, instanceMask, overlapping);
782 largestRadius = op.
radius();
783 largestRadiusIdx = op.
index();
791 template<
typename Gr
idT>
792 template<
typename InterrupterT>
797 if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
802 template<
typename Gr
idT>
803 template<
typename InterrupterT>
806 const GridT& grid,
float isovalue, InterrupterT* interrupter)
809 using ValueT =
typename GridT::ValueType;
820 mIdxTreePt.reset(
new Index32TreeT(boost::integer_traits<Index32>::const_max));
824 *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
826 if (interrupter && interrupter->wasInterrupted())
return false;
830 using Int16LeafNodeType =
typename Int16TreeT::LeafNodeType;
831 using Index32LeafNodeType =
typename Index32TreeT::LeafNodeType;
833 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
834 mSignTreePt->getNodes(signFlagsLeafNodes);
836 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
838 boost::scoped_array<Index32> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
840 tbb::parallel_for(auxiliaryLeafNodeRange,
842 (signFlagsLeafNodes, leafNodeOffsets));
846 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
847 const Index32 tmp = leafNodeOffsets[n];
852 mPointListSize = size_t(pointCount);
853 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
857 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
858 mIdxTreePt->getNodes(pointIndexLeafNodes);
861 mSurfacePointList.get(), tree, pointIndexLeafNodes,
862 signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
865 if (interrupter && interrupter->wasInterrupted())
return false;
867 Index32LeafManagerT idxLeafs(*mIdxTreePt);
869 using Index32RootNodeT =
typename Index32TreeT::RootNodeType;
870 using Index32NodeChainT =
typename Index32RootNodeT::NodeChainType;
871 BOOST_STATIC_ASSERT(boost::mpl::size<Index32NodeChainT>::value > 1);
872 using Index32InternalNodeT =
873 typename boost::mpl::at<Index32NodeChainT, boost::mpl::int_<1> >::type;
875 typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
876 nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
877 nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
879 std::vector<const Index32InternalNodeT*> internalNodes;
881 const Index32InternalNodeT* node =
nullptr;
884 if (node) internalNodes.push_back(node);
887 std::vector<IndexRange>().swap(mLeafRanges);
888 mLeafRanges.resize(internalNodes.size());
890 std::vector<const Index32LeafT*>().swap(mLeafNodes);
891 mLeafNodes.reserve(idxLeafs.leafCount());
893 typename Index32InternalNodeT::ChildOnCIter leafIt;
895 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
897 mLeafRanges[n].first = mLeafNodes.size();
899 size_t leafCount = 0;
900 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
901 mLeafNodes.push_back(&(*leafIt));
905 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
907 mLeafRanges[n].second = mLeafNodes.size();
910 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
911 mLeafBoundingSpheres.resize(mLeafNodes.size());
914 mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
918 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
919 mNodeBoundingSpheres.resize(internalNodes.size());
927 template<
typename Gr
idT>
930 std::vector<float>& distances,
bool transformPoints)
933 distances.resize(points.size(), std::numeric_limits<float>::infinity());
936 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
937 mMaxNodeLeafs, transformPoints);
945 template<
typename Gr
idT>
949 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
953 template<
typename Gr
idT>
956 std::vector<float>& distances)
958 return search(points, distances,
true);
965 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
Index64 pointCount(const PointDataTreeT &tree, const bool inCoreOnly=false)
Total points in the PointDataTree.
Definition: PointCount.h:241
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1112
uint64_t Index64
Definition: Types.h:59
Vec3< double > Vec3d
Definition: Vec3.h:679
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:797
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:230
Tolerance for floating-point comparison.
Definition: Math.h:117
Extract polygonal surfaces from scalar volumes.
uint32_t Index32
Definition: Types.h:58
static T value()
Definition: Math.h:117
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
Vec4< float > Vec4s
Definition: Vec4.h:586
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:39
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:343
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
OPENVDB_API void initialize()
Global registration of basic types.
Vec3< float > Vec3s
Definition: Vec3.h:678
SharedPtr< Grid > Ptr
Definition: Grid.h:502