OpenVDB  8.1.1
VolumeToSpheres.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
7 
8 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
10 
12 #include <openvdb/math/Math.h>
13 #include "Morphology.h" // for erodeActiveValues()
14 #include "PointScatter.h"
15 #include "LevelSetRebuild.h"
16 #include "LevelSetUtil.h"
17 #include "VolumeToMesh.h"
18 
19 #include <tbb/blocked_range.h>
20 #include <tbb/parallel_for.h>
21 #include <tbb/parallel_reduce.h>
22 
23 #include <algorithm> // for std::min(), std::max()
24 #include <cmath> // for std::sqrt()
25 #include <limits> // for std::numeric_limits
26 #include <memory>
27 #include <random>
28 #include <utility> // for std::pair
29 #include <vector>
30 
31 
32 namespace openvdb {
34 namespace OPENVDB_VERSION_NAME {
35 namespace tools {
36 
59 template<typename GridT, typename InterrupterT = util::NullInterrupter>
60 inline void
62  const GridT& grid,
63  std::vector<openvdb::Vec4s>& spheres,
64  const Vec2i& sphereCount = Vec2i(1, 50),
65  bool overlapping = false,
66  float minRadius = 1.0,
67  float maxRadius = std::numeric_limits<float>::max(),
68  float isovalue = 0.0,
69  int instanceCount = 10000,
70  InterrupterT* interrupter = nullptr);
71 
72 
74 
75 
79 template<typename GridT>
81 {
82 public:
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;
88 
100  template<typename InterrupterT = util::NullInterrupter>
101  static inline Ptr create(const GridT& grid, float isovalue = 0.0,
102  InterrupterT* interrupter = nullptr);
103 
107  inline bool search(const std::vector<Vec3R>& points, std::vector<float>& distances);
108 
112  inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
113 
115  const Index32TreeT& indexTree() const { return *mIdxTreePt; }
117  const Int16TreeT& signTree() const { return *mSignTreePt; }
118 
119 private:
120  using Index32LeafT = typename Index32TreeT::LeafNodeType;
121  using IndexRange = std::pair<size_t, size_t>;
122 
123  std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
124  std::vector<IndexRange> mLeafRanges;
125  std::vector<const Index32LeafT*> mLeafNodes;
126  PointList mSurfacePointList;
127  size_t mPointListSize = 0, mMaxNodeLeafs = 0;
128  typename Index32TreeT::Ptr mIdxTreePt;
129  typename Int16TreeT::Ptr mSignTreePt;
130 
131  ClosestSurfacePoint() = default;
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);
135 };
136 
137 
139 
140 
141 // Internal utility methods
142 
143 namespace v2s_internal {
144 
146 {
147  PointAccessor(std::vector<Vec3R>& points)
148  : mPoints(points)
149  {
150  }
151 
152  void add(const Vec3R &pos)
153  {
154  mPoints.push_back(pos);
155  }
156 private:
157  std::vector<Vec3R>& mPoints;
158 };
159 
160 
161 template<typename Index32LeafT>
162 class LeafOp
163 {
164 public:
165 
166  LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
167  const std::vector<const Index32LeafT*>& leafNodes,
168  const math::Transform& transform,
169  const PointList& surfacePointList);
170 
171  void run(bool threaded = true);
172 
173 
174  void operator()(const tbb::blocked_range<size_t>&) const;
175 
176 private:
177  std::vector<Vec4R>& mLeafBoundingSpheres;
178  const std::vector<const Index32LeafT*>& mLeafNodes;
179  const math::Transform& mTransform;
180  const PointList& mSurfacePointList;
181 };
182 
183 template<typename Index32LeafT>
185  std::vector<Vec4R>& leafBoundingSpheres,
186  const std::vector<const Index32LeafT*>& leafNodes,
187  const math::Transform& transform,
188  const PointList& surfacePointList)
189  : mLeafBoundingSpheres(leafBoundingSpheres)
190  , mLeafNodes(leafNodes)
191  , mTransform(transform)
192  , mSurfacePointList(surfacePointList)
193 {
194 }
195 
196 template<typename Index32LeafT>
197 void
199 {
200  if (threaded) {
201  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *this);
202  } else {
203  (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
204  }
205 }
206 
207 template<typename Index32LeafT>
208 void
209 LeafOp<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
210 {
211  typename Index32LeafT::ValueOnCIter iter;
212  Vec3s avg;
213 
214  for (size_t n = range.begin(); n != range.end(); ++n) {
215  avg[0] = 0.0;
216  avg[1] = 0.0;
217  avg[2] = 0.0;
218 
219  int count = 0;
220  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
221  avg += mSurfacePointList[iter.getValue()];
222  ++count;
223  }
224  if (count > 1) avg *= float(1.0 / double(count));
225 
226  float maxDist = 0.0;
227  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
228  float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
229  if (tmpDist > maxDist) maxDist = tmpDist;
230  }
231 
232  Vec4R& sphere = mLeafBoundingSpheres[n];
233  sphere[0] = avg[0];
234  sphere[1] = avg[1];
235  sphere[2] = avg[2];
236  sphere[3] = maxDist * 2.0; // padded radius
237  }
238 }
239 
240 
241 class NodeOp
242 {
243 public:
244  using IndexRange = std::pair<size_t, size_t>;
245 
246  NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
247  const std::vector<IndexRange>& leafRanges,
248  const std::vector<Vec4R>& leafBoundingSpheres);
249 
250  inline void run(bool threaded = true);
251 
252  inline void operator()(const tbb::blocked_range<size_t>&) const;
253 
254 private:
255  std::vector<Vec4R>& mNodeBoundingSpheres;
256  const std::vector<IndexRange>& mLeafRanges;
257  const std::vector<Vec4R>& mLeafBoundingSpheres;
258 };
259 
260 inline
261 NodeOp::NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
262  const std::vector<IndexRange>& leafRanges,
263  const std::vector<Vec4R>& leafBoundingSpheres)
264  : mNodeBoundingSpheres(nodeBoundingSpheres)
265  , mLeafRanges(leafRanges)
266  , mLeafBoundingSpheres(leafBoundingSpheres)
267 {
268 }
269 
270 inline void
271 NodeOp::run(bool threaded)
272 {
273  if (threaded) {
274  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *this);
275  } else {
276  (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
277  }
278 }
279 
280 inline void
281 NodeOp::operator()(const tbb::blocked_range<size_t>& range) const
282 {
283  Vec3d avg, pos;
284 
285  for (size_t n = range.begin(); n != range.end(); ++n) {
286 
287  avg[0] = 0.0;
288  avg[1] = 0.0;
289  avg[2] = 0.0;
290 
291  int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
292 
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];
297  }
298 
299  if (count > 1) avg *= float(1.0 / double(count));
300 
301 
302  double maxDist = 0.0;
303 
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];
309 
310  double tmpDist = (pos - avg).lengthSqr() + radiusSqr;
311  if (tmpDist > maxDist) maxDist = tmpDist;
312  }
313 
314  Vec4R& sphere = mNodeBoundingSpheres[n];
315 
316  sphere[0] = avg[0];
317  sphere[1] = avg[1];
318  sphere[2] = avg[2];
319  sphere[3] = maxDist * 2.0; // padded radius
320  }
321 }
322 
323 
325 
326 
327 template<typename Index32LeafT>
329 {
330 public:
331  using IndexRange = std::pair<size_t, size_t>;
332 
334  std::vector<Vec3R>& instancePoints,
335  std::vector<float>& instanceDistances,
336  const PointList& surfacePointList,
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,
341  size_t maxNodeLeafs,
342  bool transformPoints = false);
343 
344 
345  void run(bool threaded = true);
346 
347 
348  void operator()(const tbb::blocked_range<size_t>&) const;
349 
350 private:
351 
352  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
353  void evalNode(size_t pointIndex, size_t nodeIndex) const;
354 
355 
356  std::vector<Vec3R>& mInstancePoints;
357  std::vector<float>& mInstanceDistances;
358 
359  const PointList& mSurfacePointList;
360 
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;
365 
366  std::vector<float> mLeafDistances, mNodeDistances;
367 
368  const bool mTransformPoints;
369  size_t mClosestPointIndex;
370 };// ClosestPointDist
371 
372 
373 template<typename Index32LeafT>
375  std::vector<Vec3R>& instancePoints,
376  std::vector<float>& instanceDistances,
377  const PointList& surfacePointList,
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,
382  size_t maxNodeLeafs,
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)
395 {
396 }
397 
398 
399 template<typename Index32LeafT>
400 void
402 {
403  if (threaded) {
404  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
405  } else {
406  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
407  }
408 }
409 
410 template<typename Index32LeafT>
411 void
412 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
413 {
414  typename Index32LeafT::ValueOnCIter iter;
415  const Vec3s center = mInstancePoints[index];
416  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
417 
418  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
419 
420  const Vec3s& point = mSurfacePointList[iter.getValue()];
421  float tmpDist = (point - center).lengthSqr();
422 
423  if (tmpDist < mInstanceDistances[index]) {
424  mInstanceDistances[index] = tmpDist;
425  closestPointIndex = iter.getValue();
426  }
427  }
428 }
429 
430 
431 template<typename Index32LeafT>
432 void
433 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
434 {
435  if (nodeIndex >= mLeafRanges.size()) return;
436 
437  const Vec3R& pos = mInstancePoints[pointIndex];
438  float minDist = mInstanceDistances[pointIndex];
439  size_t minDistIdx = 0;
440  Vec3R center;
441  bool updatedDist = false;
442 
443  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
444  i < mLeafRanges[nodeIndex].second; ++i, ++n)
445  {
446  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
447 
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];
452 
453  distToLeaf = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
454 
455  if (distToLeaf < minDist) {
456  minDist = distToLeaf;
457  minDistIdx = i;
458  updatedDist = true;
459  }
460  }
461 
462  if (!updatedDist) return;
463 
464  evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
465 
466  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
467  i < mLeafRanges[nodeIndex].second; ++i, ++n)
468  {
469  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
470  evalLeaf(pointIndex, *mLeafNodes[i]);
471  }
472  }
473 }
474 
475 
476 template<typename Index32LeafT>
477 void
478 ClosestPointDist<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
479 {
480  Vec3R center;
481  for (size_t n = range.begin(); n != range.end(); ++n) {
482 
483  const Vec3R& pos = mInstancePoints[n];
484  float minDist = mInstanceDistances[n];
485  size_t minDistIdx = 0;
486 
487  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
488  float& distToNode = const_cast<float&>(mNodeDistances[i]);
489 
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];
494 
495  distToNode = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
496 
497  if (distToNode < minDist) {
498  minDist = distToNode;
499  minDistIdx = i;
500  }
501  }
502 
503  evalNode(n, minDistIdx);
504 
505  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
506  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
507  evalNode(n, i);
508  }
509  }
510 
511  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
512 
513  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
514  }
515 }
516 
517 
519 {
520 public:
521  UpdatePoints(
522  const Vec4s& sphere,
523  const std::vector<Vec3R>& points,
524  std::vector<float>& distances,
525  std::vector<unsigned char>& mask,
526  bool overlapping);
527 
528  float radius() const { return mRadius; }
529  int index() const { return mIndex; }
530 
531  inline void run(bool threaded = true);
532 
533 
534  UpdatePoints(UpdatePoints&, tbb::split);
535  inline void operator()(const tbb::blocked_range<size_t>& range);
536  void join(const UpdatePoints& rhs)
537  {
538  if (rhs.mRadius > mRadius) {
539  mRadius = rhs.mRadius;
540  mIndex = rhs.mIndex;
541  }
542  }
543 
544 private:
545  const Vec4s& mSphere;
546  const std::vector<Vec3R>& mPoints;
547  std::vector<float>& mDistances;
548  std::vector<unsigned char>& mMask;
549  bool mOverlapping;
550  float mRadius;
551  int mIndex;
552 };
553 
554 inline
556  const Vec4s& sphere,
557  const std::vector<Vec3R>& points,
558  std::vector<float>& distances,
559  std::vector<unsigned char>& mask,
560  bool overlapping)
561  : mSphere(sphere)
562  , mPoints(points)
563  , mDistances(distances)
564  , mMask(mask)
565  , mOverlapping(overlapping)
566  , mRadius(0.0)
567  , mIndex(0)
568 {
569 }
570 
571 inline
573  : mSphere(rhs.mSphere)
574  , mPoints(rhs.mPoints)
575  , mDistances(rhs.mDistances)
576  , mMask(rhs.mMask)
577  , mOverlapping(rhs.mOverlapping)
578  , mRadius(rhs.mRadius)
579  , mIndex(rhs.mIndex)
580 {
581 }
582 
583 inline void
584 UpdatePoints::run(bool threaded)
585 {
586  if (threaded) {
587  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
588  } else {
589  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
590  }
591 }
592 
593 inline void
594 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
595 {
596  Vec3s pos;
597  for (size_t n = range.begin(); n != range.end(); ++n) {
598  if (mMask[n]) continue;
599 
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];
603 
604  float dist = pos.length();
605 
606  if (dist < mSphere[3]) {
607  mMask[n] = 1;
608  continue;
609  }
610 
611  if (!mOverlapping) {
612  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
613  }
614 
615  if (mDistances[n] > mRadius) {
616  mRadius = mDistances[n];
617  mIndex = int(n);
618  }
619  }
620 }
621 
622 
623 } // namespace v2s_internal
624 
625 
627 
628 
629 template<typename GridT, typename InterrupterT>
630 inline void
632  const GridT& grid,
633  std::vector<openvdb::Vec4s>& spheres,
634  const Vec2i& sphereCount,
635  bool overlapping,
636  float minRadius,
637  float maxRadius,
638  float isovalue,
639  int instanceCount,
640  InterrupterT* interrupter)
641 {
642  spheres.clear();
643 
644  if (grid.empty()) return;
645 
646  const int
647  minSphereCount = sphereCount[0],
648  maxSphereCount = sphereCount[1];
649  if ((minSphereCount > maxSphereCount) || (maxSphereCount < 1)) {
650  OPENVDB_LOG_WARN("fillWithSpheres: minimum sphere count ("
651  << minSphereCount << ") exceeds maximum count (" << maxSphereCount << ")");
652  return;
653  }
654  spheres.reserve(maxSphereCount);
655 
656  auto gridPtr = grid.copy(); // shallow copy
657 
658  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
659  // Clamp the isovalue to the level set's background value minus epsilon.
660  // (In a valid narrow-band level set, all voxels, including background voxels,
661  // have values less than or equal to the background value, so an isovalue
662  // greater than or equal to the background value would produce a mask with
663  // effectively infinite extent.)
664  isovalue = std::min(isovalue,
665  static_cast<float>(gridPtr->background() - math::Tolerance<float>::value()));
666  } else if (gridPtr->getGridClass() == GRID_FOG_VOLUME) {
667  // Clamp the isovalue of a fog volume between epsilon and one,
668  // again to avoid a mask with infinite extent. (Recall that
669  // fog volume voxel values vary from zero outside to one inside.)
670  isovalue = math::Clamp(isovalue, math::Tolerance<float>::value(), 1.f);
671  }
672 
673  // ClosestSurfacePoint is inaccurate for small grids.
674  // Resample the input grid if it is too small.
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);
680 
681  auto newGridPtr = levelSetRebuild(*gridPtr, isovalue,
682  LEVEL_SET_HALF_WIDTH, LEVEL_SET_HALF_WIDTH, scaledXform.get(), interrupter);
683 
684  const auto newNumVoxels = newGridPtr->activeVoxelCount();
685  if (newNumVoxels > numVoxels) {
686  OPENVDB_LOG_DEBUG_RUNTIME("fillWithSpheres: resampled input grid from "
687  << numVoxels << " voxel" << (numVoxels == 1 ? "" : "s")
688  << " to " << newNumVoxels << " voxel" << (newNumVoxels == 1 ? "" : "s"));
689  gridPtr = newGridPtr;
690  numVoxels = newNumVoxels;
691  }
692  }
693 
694  const bool addNarrowBandPoints = (numVoxels < 10000);
695  int instances = std::max(instanceCount, maxSphereCount);
696 
697  using TreeT = typename GridT::TreeType;
698  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
699  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
700 
701  using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
702  0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>; // mt11213b
703  RandGen mtRand(/*seed=*/0);
704 
705  const TreeT& tree = gridPtr->tree();
706  math::Transform transform = gridPtr->transform();
707 
708  std::vector<Vec3R> instancePoints;
709  {
710  // Compute a mask of the voxels enclosed by the isosurface.
711  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
712  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
713  interiorMaskPtr = sdfInteriorMask(*gridPtr, isovalue);
714  } else {
715  // For non-level-set grids, the interior mask comprises the active voxels.
716  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
717  interiorMaskPtr->setTransform(transform.copy());
718  interiorMaskPtr->tree().topologyUnion(tree);
719  }
720 
721  if (interrupter && interrupter->wasInterrupted()) return;
722 
723  // If the interior mask is small and eroding it results in an empty grid,
724  // use the uneroded mask instead. (But if the minimum sphere count is zero,
725  // then eroding away the mask is acceptable.)
726  if (!addNarrowBandPoints || (minSphereCount <= 0)) {
727  tools::erodeActiveValues(interiorMaskPtr->tree(), /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
728  tools::pruneInactive(interiorMaskPtr->tree());
729  } else {
730  auto& maskTree = interiorMaskPtr->tree();
731  auto copyOfTree = StaticPtrCast<BoolTreeT>(maskTree.copy());
732  tools::erodeActiveValues(maskTree, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
733  tools::pruneInactive(maskTree);
734  if (maskTree.empty()) { interiorMaskPtr->setTree(copyOfTree); }
735  }
736 
737  // Scatter candidate sphere centroids (instancePoints)
738  instancePoints.reserve(instances);
739  v2s_internal::PointAccessor ptnAcc(instancePoints);
740 
741  const auto scatterCount = Index64(addNarrowBandPoints ? (instances / 2) : instances);
742 
744  ptnAcc, scatterCount, mtRand, 1.0, interrupter);
745  scatter(*interiorMaskPtr);
746  }
747 
748  if (interrupter && interrupter->wasInterrupted()) return;
749 
750  auto csp = ClosestSurfacePoint<GridT>::create(*gridPtr, isovalue, interrupter);
751  if (!csp) return;
752 
753  // Add extra instance points in the interior narrow band.
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());
759  if (!(volume_to_mesh_internal::EDGES & flags)
760  && (volume_to_mesh_internal::INSIDE & flags))
761  {
762  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
763  }
764  if (instancePoints.size() == size_t(instances)) break;
765  }
766  if (instancePoints.size() == size_t(instances)) break;
767  }
768  }
769 
770  if (interrupter && interrupter->wasInterrupted()) return;
771 
772  // Assign a radius to each candidate sphere. The radius is the world-space
773  // distance from the sphere's center to the closest surface point.
774  std::vector<float> instanceRadius;
775  if (!csp->search(instancePoints, instanceRadius)) return;
776 
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);
783  }
784  }
785 
786  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
787 
788  minRadius = float(minRadius * transform.voxelSize()[0]);
789  maxRadius = float(maxRadius * transform.voxelSize()[0]);
790 
791  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
792 
793  if (interrupter && interrupter->wasInterrupted()) return;
794 
795  largestRadius = std::min(maxRadius, largestRadius);
796 
797  if ((int(s) >= minSphereCount) && (largestRadius < minRadius)) break;
798 
799  const Vec4s sphere(
800  float(instancePoints[largestRadiusIdx].x()),
801  float(instancePoints[largestRadiusIdx].y()),
802  float(instancePoints[largestRadiusIdx].z()),
803  largestRadius);
804 
805  spheres.push_back(sphere);
806  instanceMask[largestRadiusIdx] = 1;
807 
809  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
810  op.run();
811 
812  largestRadius = op.radius();
813  largestRadiusIdx = op.index();
814  }
815 } // fillWithSpheres
816 
817 
819 
820 
821 template<typename GridT>
822 template<typename InterrupterT>
823 inline typename ClosestSurfacePoint<GridT>::Ptr
824 ClosestSurfacePoint<GridT>::create(const GridT& grid, float isovalue, InterrupterT* interrupter)
825 {
826  auto csp = Ptr{new ClosestSurfacePoint};
827  if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
828  return csp;
829 }
830 
831 
832 template<typename GridT>
833 template<typename InterrupterT>
834 inline bool
836  const GridT& grid, float isovalue, InterrupterT* interrupter)
837 {
838  using Index32LeafManagerT = tree::LeafManager<Index32TreeT>;
839  using ValueT = typename GridT::ValueType;
840 
841  const TreeT& tree = grid.tree();
842  const math::Transform& transform = grid.transform();
843 
844  { // Extract surface point cloud
845 
846  BoolTreeT mask(false);
848 
849  mSignTreePt.reset(new Int16TreeT(0));
850  mIdxTreePt.reset(new Index32TreeT(std::numeric_limits<Index32>::max()));
851 
852 
854  *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
855 
856  if (interrupter && interrupter->wasInterrupted()) return false;
857 
858  // count unique points
859 
860  using Int16LeafNodeType = typename Int16TreeT::LeafNodeType;
861  using Index32LeafNodeType = typename Index32TreeT::LeafNodeType;
862 
863  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
864  mSignTreePt->getNodes(signFlagsLeafNodes);
865 
866  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
867 
868  std::unique_ptr<Index32[]> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
869 
870  tbb::parallel_for(auxiliaryLeafNodeRange,
872  (signFlagsLeafNodes, leafNodeOffsets));
873 
874  {
875  Index32 pointCount = 0;
876  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
877  const Index32 tmp = leafNodeOffsets[n];
878  leafNodeOffsets[n] = pointCount;
879  pointCount += tmp;
880  }
881 
882  mPointListSize = size_t(pointCount);
883  mSurfacePointList.reset(new Vec3s[mPointListSize]);
884  }
885 
886 
887  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
888  mIdxTreePt->getNodes(pointIndexLeafNodes);
889 
890  tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
891  mSurfacePointList.get(), tree, pointIndexLeafNodes,
892  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
893  }
894 
895  if (interrupter && interrupter->wasInterrupted()) return false;
896 
897  Index32LeafManagerT idxLeafs(*mIdxTreePt);
898 
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>;
904 
905  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
906  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
907  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
908 
909  std::vector<const Index32InternalNodeT*> internalNodes;
910 
911  const Index32InternalNodeT* node = nullptr;
912  for (; nIt; ++nIt) {
913  nIt.getNode(node);
914  if (node) internalNodes.push_back(node);
915  }
916 
917  std::vector<IndexRange>().swap(mLeafRanges);
918  mLeafRanges.resize(internalNodes.size());
919 
920  std::vector<const Index32LeafT*>().swap(mLeafNodes);
921  mLeafNodes.reserve(idxLeafs.leafCount());
922 
923  typename Index32InternalNodeT::ChildOnCIter leafIt;
924  mMaxNodeLeafs = 0;
925  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
926 
927  mLeafRanges[n].first = mLeafNodes.size();
928 
929  size_t leafCount = 0;
930  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
931  mLeafNodes.push_back(&(*leafIt));
932  ++leafCount;
933  }
934 
935  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
936 
937  mLeafRanges[n].second = mLeafNodes.size();
938  }
939 
940  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
941  mLeafBoundingSpheres.resize(mLeafNodes.size());
942 
944  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
945  leafBS.run();
946 
947 
948  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
949  mNodeBoundingSpheres.resize(internalNodes.size());
950 
951  v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
952  nodeBS.run();
953  return true;
954 } // ClosestSurfacePoint::initialize
955 
956 
957 template<typename GridT>
958 inline bool
959 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
960  std::vector<float>& distances, bool transformPoints)
961 {
962  distances.clear();
963  distances.resize(points.size(), std::numeric_limits<float>::infinity());
964 
965  v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
966  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
967  mMaxNodeLeafs, transformPoints);
968 
969  cpd.run();
970 
971  return true;
972 }
973 
974 
975 template<typename GridT>
976 inline bool
977 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
978 {
979  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
980 }
981 
982 
983 template<typename GridT>
984 inline bool
986  std::vector<float>& distances)
987 {
988  return search(points, distances, true);
989 }
990 
991 } // namespace tools
992 } // namespace OPENVDB_VERSION_NAME
993 } // namespace openvdb
994 
995 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
typename TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToSpheres.h:85
void identifySurfaceIntersectingVoxels(typename InputTreeType::template ValueConverter< bool >::Type &intersectionTree, const InputTreeType &inputTree, typename InputTreeType::ValueType isovalue)
Definition: VolumeToMesh.h:3322
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: tools/PointScatter.h:87
Vec2< int32_t > Vec2i
Definition: Vec2.h:534
void join(const UpdatePoints &rhs)
Definition: VolumeToSpheres.h:536
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Overwrite each input point with its closest surface point.
Definition: VolumeToSpheres.h:985
ClosestPointDist(std::vector< Vec3R > &instancePoints, std::vector< float > &instanceDistances, const PointList &surfacePointList, const std::vector< const Index32LeafT * > &leafNodes, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres, const std::vector< Vec4R > &nodeBoundingSpheres, size_t maxNodeLeafs, bool transformPoints=false)
Definition: VolumeToSpheres.h:374
Definition: Morphology.h:56
uint32_t Index32
Definition: openvdb/Types.h:48
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2270
void operator()(const tbb::blocked_range< size_t > &range)
Definition: VolumeToSpheres.h:594
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
#define OPENVDB_LOG_DEBUG_RUNTIME(message)
Log a debugging message in both debug and optimized builds.
Definition: logging.h:267
static Ptr create(const GridT &grid, float isovalue=0.0, InterrupterT *interrupter=nullptr)
Extract surface points and construct a spatial acceleration structure.
Definition: VolumeToSpheres.h:824
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
We offer three different algorithms (each in its own class) for scattering of points in active voxels...
typename GridT::TreeType TreeT
Definition: VolumeToSpheres.h:84
OPENVDB_IMPORT void initialize()
Global registration of basic types.
void run(bool threaded=true)
Definition: VolumeToSpheres.h:198
Definition: VolumeToSpheres.h:145
Definition: Mat4.h:24
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:108
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:244
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:93
Tolerance for floating-point comparison.
Definition: Math.h:147
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Accelerated closest surface point queries for narrow band level sets.
Definition: VolumeToSpheres.h:80
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:260
void run(bool threaded=true)
Definition: VolumeToSpheres.h:584
Extract polygonal surfaces from scalar volumes.
void run(bool threaded=true)
Definition: VolumeToSpheres.h:401
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
Vec3< double > Vec3d
Definition: Vec3.h:668
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:354
void add(const Vec3R &pos)
Definition: VolumeToSpheres.h:152
float radius() const
Definition: VolumeToSpheres.h:528
float Cbrt(float x)
Return the cube root of a floating-point value.
Definition: Math.h:772
std::unique_ptr< openvdb::Vec3s[]> PointList
Point and primitive list types.
Definition: VolumeToMesh.h:150
GridType::Ptr levelSetRebuild(const GridType &grid, float isovalue=0, float halfWidth=float(LEVEL_SET_HALF_WIDTH), const math::Transform *xform=nullptr)
Return a new grid of type GridType that contains a narrow-band level set representation of an isosurf...
Definition: LevelSetRebuild.h:311
Implementation of morphological dilation and erosion.
typename TreeT::template ValueConverter< Index32 >::Type Index32TreeT
Definition: VolumeToSpheres.h:86
Definition: openvdb/Exceptions.h:13
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:637
Definition: Morphology.h:78
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Definition: PointCount.h:88
uint64_t Index64
Definition: openvdb/Types.h:49
const Index32TreeT & indexTree() const
Tree accessor.
Definition: VolumeToSpheres.h:115
PointAccessor(std::vector< Vec3R > &points)
Definition: VolumeToSpheres.h:147
std::unique_ptr< ClosestSurfacePoint > Ptr
Definition: VolumeToSpheres.h:83
NodeOp(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
Definition: VolumeToSpheres.h:261
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
Definition: VolumeToSpheres.h:555
typename TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToSpheres.h:87
Definition: openvdb/Types.h:334
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Compute the distance from each input point to its closest surface point.
Definition: VolumeToSpheres.h:977
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:478
Vec4< float > Vec4s
Definition: Vec4.h:565
Vec3< float > Vec3s
Definition: Vec3.h:667
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
static const Real LEVEL_SET_HALF_WIDTH
Definition: openvdb/Types.h:339
Definition: openvdb/Types.h:333
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1247
Definition: VolumeToSpheres.h:518
TreeType & tree()
Return a reference to this grid&#39;s tree, which might be shared with other grids.
Definition: Grid.h:917
#define OPENVDB_LOG_WARN(message)
Log a warning message of the form &#39;someVar << "some text" << ...&#39;.
Definition: logging.h:253
Ptr copy() const
Definition: Transform.h:50
Definition: VolumeToSpheres.h:162
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
int index() const
Definition: VolumeToSpheres.h:529
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:281
SharedPtr< Grid > Ptr
Definition: Grid.h:579
void erodeActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically erode all active values (i.e. both voxels and tiles) in a tree using one of three neare...
Definition: Morphology.h:1126
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:209
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:302
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
void setTree(TreeBase::Ptr) override
Associate the given tree with this grid, in place of its existing tree.
Definition: Grid.h:1473
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
void computeAuxiliaryData(typename InputTreeType::template ValueConverter< Int16 >::Type &signFlagsTree, typename InputTreeType::template ValueConverter< Index32 >::Type &pointIndexTree, const typename InputTreeType::template ValueConverter< bool >::Type &intersectionTree, const InputTreeType &inputTree, typename InputTreeType::ValueType isovalue)
Definition: VolumeToMesh.h:3864
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:331
static T value()
Definition: Math.h:147
void run(bool threaded=true)
Definition: VolumeToSpheres.h:271
void fillWithSpheres(const GridT &grid, std::vector< openvdb::Vec4s > &spheres, const Vec2i &sphereCount=Vec2i(1, 50), bool overlapping=false, float minRadius=1.0, float maxRadius=std::numeric_limits< float >::max(), float isovalue=0.0, int instanceCount=10000, InterrupterT *interrupter=nullptr)
Fill a closed level set or fog volume with adaptively-sized spheres.
Definition: VolumeToSpheres.h:631
const Int16TreeT & signTree() const
Tree accessor.
Definition: VolumeToSpheres.h:117
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
Definition: Transform.h:39
Definition: VolumeToSpheres.h:241