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