8 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
19 #include <tbb/blocked_range.h>
20 #include <tbb/parallel_for.h>
21 #include <tbb/parallel_reduce.h>
59 template<
typename Gr
idT,
typename InterrupterT = util::NullInterrupter>
63 std::vector<openvdb::Vec4s>& spheres,
65 bool overlapping =
false,
66 float minRadius = 1.0,
69 int instanceCount = 10000,
70 InterrupterT* interrupter =
nullptr);
79 template<
typename Gr
idT>
83 using Ptr = std::unique_ptr<ClosestSurfacePoint>;
84 using TreeT =
typename GridT::TreeType;
85 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
86 using Index32TreeT =
typename TreeT::template ValueConverter<Index32>::Type;
87 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
100 template<
typename InterrupterT = util::NullInterrupter>
101 static inline Ptr create(
const GridT& grid,
float isovalue = 0.0,
102 InterrupterT* interrupter =
nullptr);
107 inline bool search(
const std::vector<Vec3R>& points, std::vector<float>& distances);
112 inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
120 using Index32LeafT =
typename Index32TreeT::LeafNodeType;
121 using IndexRange = std::pair<size_t, size_t>;
123 std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
124 std::vector<IndexRange> mLeafRanges;
125 std::vector<const Index32LeafT*> mLeafNodes;
127 size_t mPointListSize = 0, mMaxNodeLeafs = 0;
128 typename Index32TreeT::Ptr mIdxTreePt;
129 typename Int16TreeT::Ptr mSignTreePt;
132 template<
typename InterrupterT = util::NullInterrupter>
133 inline bool initialize(
const GridT&,
float isovalue, InterrupterT*);
134 inline bool search(std::vector<Vec3R>&, std::vector<float>&,
bool transformPoints);
143 namespace v2s_internal {
154 mPoints.push_back(pos);
157 std::vector<Vec3R>& mPoints;
161 template<
typename Index32LeafT>
166 LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
167 const std::vector<const Index32LeafT*>& leafNodes,
171 void run(
bool threaded =
true);
174 void operator()(
const tbb::blocked_range<size_t>&)
const;
177 std::vector<Vec4R>& mLeafBoundingSpheres;
178 const std::vector<const Index32LeafT*>& mLeafNodes;
183 template<
typename Index32LeafT>
185 std::vector<Vec4R>& leafBoundingSpheres,
186 const std::vector<const Index32LeafT*>& leafNodes,
189 : mLeafBoundingSpheres(leafBoundingSpheres)
190 , mLeafNodes(leafNodes)
191 , mTransform(transform)
192 , mSurfacePointList(surfacePointList)
196 template<
typename Index32LeafT>
201 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *
this);
203 (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
207 template<
typename Index32LeafT>
211 typename Index32LeafT::ValueOnCIter iter;
214 for (
size_t n = range.begin(); n != range.end(); ++n) {
220 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
221 avg += mSurfacePointList[iter.getValue()];
224 if (count > 1) avg *= float(1.0 /
double(count));
227 for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
228 float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
229 if (tmpDist > maxDist) maxDist = tmpDist;
232 Vec4R& sphere = mLeafBoundingSpheres[n];
236 sphere[3] = maxDist * 2.0;
246 NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
247 const std::vector<IndexRange>& leafRanges,
248 const std::vector<Vec4R>& leafBoundingSpheres);
250 inline void run(
bool threaded =
true);
252 inline void operator()(
const tbb::blocked_range<size_t>&)
const;
255 std::vector<Vec4R>& mNodeBoundingSpheres;
256 const std::vector<IndexRange>& mLeafRanges;
257 const std::vector<Vec4R>& mLeafBoundingSpheres;
262 const std::vector<IndexRange>& leafRanges,
263 const std::vector<Vec4R>& leafBoundingSpheres)
264 : mNodeBoundingSpheres(nodeBoundingSpheres)
265 , mLeafRanges(leafRanges)
266 , mLeafBoundingSpheres(leafBoundingSpheres)
274 tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *
this);
276 (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
285 for (
size_t n = range.begin(); n != range.end(); ++n) {
291 int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
293 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
294 avg[0] += mLeafBoundingSpheres[i][0];
295 avg[1] += mLeafBoundingSpheres[i][1];
296 avg[2] += mLeafBoundingSpheres[i][2];
299 if (count > 1) avg *= float(1.0 /
double(count));
302 double maxDist = 0.0;
304 for (
size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
305 pos[0] = mLeafBoundingSpheres[i][0];
306 pos[1] = mLeafBoundingSpheres[i][1];
307 pos[2] = mLeafBoundingSpheres[i][2];
308 const auto radiusSqr = mLeafBoundingSpheres[i][3];
310 double tmpDist = (pos - avg).lengthSqr() + radiusSqr;
311 if (tmpDist > maxDist) maxDist = tmpDist;
314 Vec4R& sphere = mNodeBoundingSpheres[n];
319 sphere[3] = maxDist * 2.0;
327 template<
typename Index32LeafT>
334 std::vector<Vec3R>& instancePoints,
335 std::vector<float>& instanceDistances,
337 const std::vector<const Index32LeafT*>& leafNodes,
338 const std::vector<IndexRange>& leafRanges,
339 const std::vector<Vec4R>& leafBoundingSpheres,
340 const std::vector<Vec4R>& nodeBoundingSpheres,
342 bool transformPoints =
false);
345 void run(
bool threaded =
true);
348 void operator()(
const tbb::blocked_range<size_t>&)
const;
352 void evalLeaf(
size_t index,
const Index32LeafT& leaf)
const;
353 void evalNode(
size_t pointIndex,
size_t nodeIndex)
const;
356 std::vector<Vec3R>& mInstancePoints;
357 std::vector<float>& mInstanceDistances;
361 const std::vector<const Index32LeafT*>& mLeafNodes;
362 const std::vector<IndexRange>& mLeafRanges;
363 const std::vector<Vec4R>& mLeafBoundingSpheres;
364 const std::vector<Vec4R>& mNodeBoundingSpheres;
366 std::vector<float> mLeafDistances, mNodeDistances;
368 const bool mTransformPoints;
369 size_t mClosestPointIndex;
373 template<
typename Index32LeafT>
375 std::vector<Vec3R>& instancePoints,
376 std::vector<float>& instanceDistances,
378 const std::vector<const Index32LeafT*>& leafNodes,
379 const std::vector<IndexRange>& leafRanges,
380 const std::vector<Vec4R>& leafBoundingSpheres,
381 const std::vector<Vec4R>& nodeBoundingSpheres,
383 bool transformPoints)
384 : mInstancePoints(instancePoints)
385 , mInstanceDistances(instanceDistances)
386 , mSurfacePointList(surfacePointList)
387 , mLeafNodes(leafNodes)
388 , mLeafRanges(leafRanges)
389 , mLeafBoundingSpheres(leafBoundingSpheres)
390 , mNodeBoundingSpheres(nodeBoundingSpheres)
391 , mLeafDistances(maxNodeLeafs, 0.0)
392 , mNodeDistances(leafRanges.size(), 0.0)
393 , mTransformPoints(transformPoints)
394 , mClosestPointIndex(0)
399 template<
typename Index32LeafT>
404 tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *
this);
406 (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
410 template<
typename Index32LeafT>
414 typename Index32LeafT::ValueOnCIter iter;
415 const Vec3s center = mInstancePoints[index];
416 size_t& closestPointIndex =
const_cast<size_t&
>(mClosestPointIndex);
418 for (iter = leaf.cbeginValueOn(); iter; ++iter) {
420 const Vec3s& point = mSurfacePointList[iter.getValue()];
421 float tmpDist = (point - center).lengthSqr();
423 if (tmpDist < mInstanceDistances[index]) {
424 mInstanceDistances[index] = tmpDist;
425 closestPointIndex = iter.getValue();
431 template<
typename Index32LeafT>
433 ClosestPointDist<Index32LeafT>::evalNode(
size_t pointIndex,
size_t nodeIndex)
const
435 if (nodeIndex >= mLeafRanges.size())
return;
437 const Vec3R& pos = mInstancePoints[pointIndex];
438 float minDist = mInstanceDistances[pointIndex];
439 size_t minDistIdx = 0;
441 bool updatedDist =
false;
443 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
444 i < mLeafRanges[nodeIndex].second; ++i, ++n)
446 float& distToLeaf =
const_cast<float&
>(mLeafDistances[n]);
448 center[0] = mLeafBoundingSpheres[i][0];
449 center[1] = mLeafBoundingSpheres[i][1];
450 center[2] = mLeafBoundingSpheres[i][2];
451 const auto radiusSqr = mLeafBoundingSpheres[i][3];
453 distToLeaf = float(
std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
455 if (distToLeaf < minDist) {
456 minDist = distToLeaf;
462 if (!updatedDist)
return;
464 evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
466 for (
size_t i = mLeafRanges[nodeIndex].first, n = 0;
467 i < mLeafRanges[nodeIndex].second; ++i, ++n)
469 if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
470 evalLeaf(pointIndex, *mLeafNodes[i]);
476 template<
typename Index32LeafT>
481 for (
size_t n = range.begin(); n != range.end(); ++n) {
483 const Vec3R& pos = mInstancePoints[n];
484 float minDist = mInstanceDistances[n];
485 size_t minDistIdx = 0;
487 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
488 float& distToNode =
const_cast<float&
>(mNodeDistances[i]);
490 center[0] = mNodeBoundingSpheres[i][0];
491 center[1] = mNodeBoundingSpheres[i][1];
492 center[2] = mNodeBoundingSpheres[i][2];
493 const auto radiusSqr = mNodeBoundingSpheres[i][3];
495 distToNode = float(
std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
497 if (distToNode < minDist) {
498 minDist = distToNode;
503 evalNode(n, minDistIdx);
505 for (
size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
506 if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
511 mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
513 if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
523 const std::vector<Vec3R>& points,
524 std::vector<float>& distances,
525 std::vector<unsigned char>& mask,
529 int index()
const {
return mIndex; }
531 inline void run(
bool threaded =
true);
535 inline void operator()(
const tbb::blocked_range<size_t>& range);
538 if (rhs.mRadius > mRadius) {
539 mRadius = rhs.mRadius;
545 const Vec4s& mSphere;
546 const std::vector<Vec3R>& mPoints;
547 std::vector<float>& mDistances;
548 std::vector<unsigned char>& mMask;
557 const std::vector<Vec3R>& points,
558 std::vector<float>& distances,
559 std::vector<unsigned char>& mask,
563 , mDistances(distances)
565 , mOverlapping(overlapping)
573 : mSphere(rhs.mSphere)
574 , mPoints(rhs.mPoints)
575 , mDistances(rhs.mDistances)
577 , mOverlapping(rhs.mOverlapping)
578 , mRadius(rhs.mRadius)
587 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *
this);
589 (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
597 for (
size_t n = range.begin(); n != range.end(); ++n) {
598 if (mMask[n])
continue;
600 pos.x() = float(mPoints[n].x()) - mSphere[0];
601 pos.y() = float(mPoints[n].y()) - mSphere[1];
602 pos.z() = float(mPoints[n].z()) - mSphere[2];
604 float dist = pos.length();
606 if (dist < mSphere[3]) {
612 mDistances[n] =
std::min(mDistances[n], (dist - mSphere[3]));
615 if (mDistances[n] > mRadius) {
616 mRadius = mDistances[n];
629 template<
typename Gr
idT,
typename InterrupterT>
633 std::vector<openvdb::Vec4s>& spheres,
634 const Vec2i& sphereCount,
640 InterrupterT* interrupter)
644 if (grid.empty())
return;
647 minSphereCount = sphereCount[0],
648 maxSphereCount = sphereCount[1];
649 if ((minSphereCount > maxSphereCount) || (maxSphereCount < 1)) {
651 << minSphereCount <<
") exceeds maximum count (" << maxSphereCount <<
")");
654 spheres.reserve(maxSphereCount);
656 auto gridPtr = grid.copy();
675 auto numVoxels = gridPtr->activeVoxelCount();
676 if (numVoxels < 10000) {
677 const auto scale = 1.0 /
math::Cbrt(2.0 * 10000.0 /
double(numVoxels));
678 auto scaledXform = gridPtr->transform().copy();
679 scaledXform->preScale(
scale);
684 const auto newNumVoxels = newGridPtr->activeVoxelCount();
685 if (newNumVoxels > numVoxels) {
687 << numVoxels <<
" voxel" << (numVoxels == 1 ?
"" :
"s")
688 <<
" to " << newNumVoxels <<
" voxel" << (newNumVoxels == 1 ?
"" :
"s"));
689 gridPtr = newGridPtr;
690 numVoxels = newNumVoxels;
694 const bool addNarrowBandPoints = (numVoxels < 10000);
695 int instances =
std::max(instanceCount, maxSphereCount);
697 using TreeT =
typename GridT::TreeType;
698 using BoolTreeT =
typename TreeT::template ValueConverter<bool>::Type;
699 using Int16TreeT =
typename TreeT::template ValueConverter<Int16>::Type;
701 using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
702 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>;
705 const TreeT& tree = gridPtr->tree();
708 std::vector<Vec3R> instancePoints;
717 interiorMaskPtr->setTransform(transform.
copy());
718 interiorMaskPtr->
tree().topologyUnion(tree);
721 if (interrupter && interrupter->wasInterrupted())
return;
726 if (!addNarrowBandPoints || (minSphereCount <= 0)) {
730 auto& maskTree = interiorMaskPtr->
tree();
731 auto copyOfTree = StaticPtrCast<BoolTreeT>(maskTree.copy());
734 if (maskTree.empty()) { interiorMaskPtr->
setTree(copyOfTree); }
738 instancePoints.reserve(instances);
741 const auto scatterCount =
Index64(addNarrowBandPoints ? (instances / 2) : instances);
744 ptnAcc, scatterCount, mtRand, 1.0, interrupter);
745 scatter(*interiorMaskPtr);
748 if (interrupter && interrupter->wasInterrupted())
return;
754 if (instancePoints.size() <
size_t(instances)) {
755 const Int16TreeT& signTree = csp->signTree();
756 for (
auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
757 for (
auto it = leafIt->cbeginValueOn(); it; ++it) {
758 const int flags = int(it.getValue());
762 instancePoints.push_back(transform.
indexToWorld(it.getCoord()));
764 if (instancePoints.size() ==
size_t(instances))
break;
766 if (instancePoints.size() ==
size_t(instances))
break;
770 if (interrupter && interrupter->wasInterrupted())
return;
774 std::vector<float> instanceRadius;
775 if (!csp->search(instancePoints, instanceRadius))
return;
777 float largestRadius = 0.0;
778 int largestRadiusIdx = 0;
779 for (
size_t n = 0, N = instancePoints.size(); n < N; ++n) {
780 if (instanceRadius[n] > largestRadius) {
781 largestRadius = instanceRadius[n];
782 largestRadiusIdx = int(n);
786 std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
788 minRadius = float(minRadius * transform.
voxelSize()[0]);
789 maxRadius = float(maxRadius * transform.
voxelSize()[0]);
791 for (
size_t s = 0, S =
std::min(
size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
793 if (interrupter && interrupter->wasInterrupted())
return;
795 largestRadius =
std::min(maxRadius, largestRadius);
797 if ((
int(s) >= minSphereCount) && (largestRadius < minRadius))
break;
800 float(instancePoints[largestRadiusIdx].x()),
801 float(instancePoints[largestRadiusIdx].y()),
802 float(instancePoints[largestRadiusIdx].z()),
805 spheres.push_back(sphere);
806 instanceMask[largestRadiusIdx] = 1;
809 sphere, instancePoints, instanceRadius, instanceMask, overlapping);
812 largestRadius = op.
radius();
813 largestRadiusIdx = op.
index();
821 template<
typename Gr
idT>
822 template<
typename InterrupterT>
827 if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
832 template<
typename Gr
idT>
833 template<
typename InterrupterT>
836 const GridT& grid,
float isovalue, InterrupterT* interrupter)
839 using ValueT =
typename GridT::ValueType;
841 const TreeT& tree = grid.
tree();
846 BoolTreeT mask(
false);
849 mSignTreePt.reset(
new Int16TreeT(0));
854 *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
856 if (interrupter && interrupter->wasInterrupted())
return false;
860 using Int16LeafNodeType =
typename Int16TreeT::LeafNodeType;
861 using Index32LeafNodeType =
typename Index32TreeT::LeafNodeType;
863 std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
864 mSignTreePt->getNodes(signFlagsLeafNodes);
866 const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
868 std::unique_ptr<Index32[]> leafNodeOffsets(
new Index32[signFlagsLeafNodes.size()]);
870 tbb::parallel_for(auxiliaryLeafNodeRange,
872 (signFlagsLeafNodes, leafNodeOffsets));
876 for (
size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
877 const Index32 tmp = leafNodeOffsets[n];
883 mSurfacePointList.reset(
new Vec3s[mPointListSize]);
887 std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
888 mIdxTreePt->getNodes(pointIndexLeafNodes);
890 tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
891 mSurfacePointList.get(), tree, pointIndexLeafNodes,
892 signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
895 if (interrupter && interrupter->wasInterrupted())
return false;
897 Index32LeafManagerT idxLeafs(*mIdxTreePt);
899 using Index32RootNodeT =
typename Index32TreeT::RootNodeType;
900 using Index32NodeChainT =
typename Index32RootNodeT::NodeChainType;
901 static_assert(Index32NodeChainT::Size > 1,
902 "expected tree depth greater than one");
903 using Index32InternalNodeT =
typename Index32NodeChainT::template Get<1>;
905 typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
906 nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
907 nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
909 std::vector<const Index32InternalNodeT*> internalNodes;
911 const Index32InternalNodeT* node =
nullptr;
914 if (node) internalNodes.push_back(node);
917 std::vector<IndexRange>().swap(mLeafRanges);
918 mLeafRanges.resize(internalNodes.size());
920 std::vector<const Index32LeafT*>().swap(mLeafNodes);
921 mLeafNodes.reserve(idxLeafs.leafCount());
923 typename Index32InternalNodeT::ChildOnCIter leafIt;
925 for (
size_t n = 0, N = internalNodes.size(); n < N; ++n) {
927 mLeafRanges[n].first = mLeafNodes.size();
929 size_t leafCount = 0;
930 for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
931 mLeafNodes.push_back(&(*leafIt));
935 mMaxNodeLeafs =
std::max(leafCount, mMaxNodeLeafs);
937 mLeafRanges[n].second = mLeafNodes.size();
940 std::vector<Vec4R>().swap(mLeafBoundingSpheres);
941 mLeafBoundingSpheres.resize(mLeafNodes.size());
943 v2s_internal::LeafOp<Index32LeafT> leafBS(
944 mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
948 std::vector<Vec4R>().swap(mNodeBoundingSpheres);
949 mNodeBoundingSpheres.resize(internalNodes.size());
951 v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
957 template<
typename Gr
idT>
960 std::vector<float>& distances,
bool transformPoints)
963 distances.resize(points.size(), std::numeric_limits<float>::infinity());
965 v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
966 mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
967 mMaxNodeLeafs, transformPoints);
975 template<
typename Gr
idT>
979 return search(
const_cast<std::vector<Vec3R>&
>(points), distances,
false);
983 template<
typename Gr
idT>
986 std::vector<float>& distances)
988 return search(points, distances,
true);
995 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED