HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VolumeToSpheres.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 /// @file tools/VolumeToSpheres.h
5 ///
6 /// @brief Fill a closed level set or fog volume with adaptively-sized spheres.
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 
37 /// @brief Fill a closed level set or fog volume with adaptively-sized spheres.
38 ///
39 /// @param grid a scalar grid that defines the surface to be filled with spheres
40 /// @param spheres an output array of 4-tuples representing the fitted spheres<BR>
41 /// The first three components of each tuple specify the sphere center,
42 /// and the fourth specifies the radius.
43 /// The spheres are ordered by radius, from largest to smallest.
44 /// @param sphereCount lower and upper bounds on the number of spheres to be generated<BR>
45 /// The actual number will be somewhere within the bounds.
46 /// @param overlapping toggle to allow spheres to overlap/intersect
47 /// @param minRadius the smallest allowable sphere size, in voxel units<BR>
48 /// @param maxRadius the largest allowable sphere size, in voxel units
49 /// @param isovalue the voxel value that determines the surface of the volume<BR>
50 /// The default value of zero works for signed distance fields,
51 /// while fog volumes require a larger positive value
52 /// (0.5 is a good initial guess).
53 /// @param instanceCount the number of interior points to consider for the sphere placement<BR>
54 /// Increasing this count increases the chances of finding optimal
55 /// sphere sizes.
56 /// @param interrupter pointer to an object adhering to the util::NullInterrupter interface
57 ///
58 /// @note The minimum sphere count takes precedence over the minimum radius.
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 
73 ////////////////////////////////////////
74 
75 
76 /// @brief Accelerated closest surface point queries for narrow band level sets
77 /// @details Supports queries that originate at arbitrary world-space locations,
78 /// is not confined to the narrow band region of the input volume geometry.
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 
89  /// @brief Extract surface points and construct a spatial acceleration structure.
90  ///
91  /// @return a null pointer if the initialization fails for any reason,
92  /// otherwise a unique pointer to a newly-allocated ClosestSurfacePoint object.
93  ///
94  /// @param grid a scalar level set or fog volume
95  /// @param isovalue the voxel value that determines the surface of the volume
96  /// The default value of zero works for signed distance fields,
97  /// while fog volumes require a larger positive value
98  /// (0.5 is a good initial guess).
99  /// @param interrupter pointer to an object adhering to the util::NullInterrupter interface.
100  template<typename InterrupterT = util::NullInterrupter>
101  static inline Ptr create(const GridT& grid, float isovalue = 0.0,
102  InterrupterT* interrupter = nullptr);
103 
104  /// @brief Compute the distance from each input point to its closest surface point.
105  /// @param points input list of points in world space
106  /// @param distances output list of closest surface point distances
107  inline bool search(const std::vector<Vec3R>& points, std::vector<float>& distances);
108 
109  /// @brief Overwrite each input point with its closest surface point.
110  /// @param points input/output list of points in world space
111  /// @param distances output list of closest surface point distances
112  inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
113 
114  /// @brief Tree accessor
115  const Index32TreeT& indexTree() const { return *mIdxTreePt; }
116  /// @brief Tree accessor
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 
138 ////////////////////////////////////////
139 
140 
141 // Internal utility methods
142 
143 /// @cond OPENVDB_DOCS_INTERNAL
144 
145 namespace v2s_internal {
146 
147 struct PointAccessor
148 {
149  PointAccessor(std::vector<Vec3R>& points)
150  : mPoints(points)
151  {
152  }
153 
154  void add(const Vec3R &pos)
155  {
156  mPoints.push_back(pos);
157  }
158 private:
159  std::vector<Vec3R>& mPoints;
160 };
161 
162 
163 template<typename Index32LeafT>
164 class LeafOp
165 {
166 public:
167 
168  LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
169  const std::vector<const Index32LeafT*>& leafNodes,
170  const math::Transform& transform,
171  const PointList& surfacePointList);
172 
173  void run(bool threaded = true);
174 
175 
176  void operator()(const tbb::blocked_range<size_t>&) const;
177 
178 private:
179  std::vector<Vec4R>& mLeafBoundingSpheres;
180  const std::vector<const Index32LeafT*>& mLeafNodes;
181  const math::Transform& mTransform;
182  const PointList& mSurfacePointList;
183 };
184 
185 template<typename Index32LeafT>
186 LeafOp<Index32LeafT>::LeafOp(
187  std::vector<Vec4R>& leafBoundingSpheres,
188  const std::vector<const Index32LeafT*>& leafNodes,
189  const math::Transform& transform,
190  const PointList& surfacePointList)
191  : mLeafBoundingSpheres(leafBoundingSpheres)
192  , mLeafNodes(leafNodes)
193  , mTransform(transform)
194  , mSurfacePointList(surfacePointList)
195 {
196 }
197 
198 template<typename Index32LeafT>
199 void
200 LeafOp<Index32LeafT>::run(bool threaded)
201 {
202  if (threaded) {
203  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *this);
204  } else {
205  (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
206  }
207 }
208 
209 template<typename Index32LeafT>
210 void
211 LeafOp<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
212 {
213  typename Index32LeafT::ValueOnCIter iter;
214  Vec3s avg;
215 
216  for (size_t n = range.begin(); n != range.end(); ++n) {
217  avg[0] = 0.0;
218  avg[1] = 0.0;
219  avg[2] = 0.0;
220 
221  int count = 0;
222  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
223  avg += mSurfacePointList[iter.getValue()];
224  ++count;
225  }
226  if (count > 1) avg *= float(1.0 / double(count));
227 
228  float maxDist = 0.0;
229  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
230  float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
231  if (tmpDist > maxDist) maxDist = tmpDist;
232  }
233 
234  Vec4R& sphere = mLeafBoundingSpheres[n];
235  sphere[0] = avg[0];
236  sphere[1] = avg[1];
237  sphere[2] = avg[2];
238  sphere[3] = maxDist * 2.0; // padded radius
239  }
240 }
241 
242 
243 class NodeOp
244 {
245 public:
246  using IndexRange = std::pair<size_t, size_t>;
247 
248  NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
249  const std::vector<IndexRange>& leafRanges,
250  const std::vector<Vec4R>& leafBoundingSpheres);
251 
252  inline void run(bool threaded = true);
253 
254  inline void operator()(const tbb::blocked_range<size_t>&) const;
255 
256 private:
257  std::vector<Vec4R>& mNodeBoundingSpheres;
258  const std::vector<IndexRange>& mLeafRanges;
259  const std::vector<Vec4R>& mLeafBoundingSpheres;
260 };
261 
262 inline
263 NodeOp::NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
264  const std::vector<IndexRange>& leafRanges,
265  const std::vector<Vec4R>& leafBoundingSpheres)
266  : mNodeBoundingSpheres(nodeBoundingSpheres)
267  , mLeafRanges(leafRanges)
268  , mLeafBoundingSpheres(leafBoundingSpheres)
269 {
270 }
271 
272 inline void
273 NodeOp::run(bool threaded)
274 {
275  if (threaded) {
276  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *this);
277  } else {
278  (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
279  }
280 }
281 
282 inline void
283 NodeOp::operator()(const tbb::blocked_range<size_t>& range) const
284 {
285  Vec3d avg, pos;
286 
287  for (size_t n = range.begin(); n != range.end(); ++n) {
288 
289  avg[0] = 0.0;
290  avg[1] = 0.0;
291  avg[2] = 0.0;
292 
293  int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
294 
295  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
296  avg[0] += mLeafBoundingSpheres[i][0];
297  avg[1] += mLeafBoundingSpheres[i][1];
298  avg[2] += mLeafBoundingSpheres[i][2];
299  }
300 
301  if (count > 1) avg *= float(1.0 / double(count));
302 
303 
304  double maxDist = 0.0;
305 
306  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
307  pos[0] = mLeafBoundingSpheres[i][0];
308  pos[1] = mLeafBoundingSpheres[i][1];
309  pos[2] = mLeafBoundingSpheres[i][2];
310  const auto radiusSqr = mLeafBoundingSpheres[i][3];
311 
312  double tmpDist = (pos - avg).lengthSqr() + radiusSqr;
313  if (tmpDist > maxDist) maxDist = tmpDist;
314  }
315 
316  Vec4R& sphere = mNodeBoundingSpheres[n];
317 
318  sphere[0] = avg[0];
319  sphere[1] = avg[1];
320  sphere[2] = avg[2];
321  sphere[3] = maxDist * 2.0; // padded radius
322  }
323 }
324 
325 
326 ////////////////////////////////////////
327 
328 
329 template<typename Index32LeafT>
330 class ClosestPointDist
331 {
332 public:
333  using IndexRange = std::pair<size_t, size_t>;
334 
335  ClosestPointDist(
336  std::vector<Vec3R>& instancePoints,
337  std::vector<float>& instanceDistances,
338  const PointList& surfacePointList,
339  const std::vector<const Index32LeafT*>& leafNodes,
340  const std::vector<IndexRange>& leafRanges,
341  const std::vector<Vec4R>& leafBoundingSpheres,
342  const std::vector<Vec4R>& nodeBoundingSpheres,
343  size_t maxNodeLeafs,
344  bool transformPoints = false);
345 
346 
347  void run(bool threaded = true);
348 
349 
350  void operator()(const tbb::blocked_range<size_t>&) const;
351 
352 private:
353 
354  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
355  void evalNode(size_t pointIndex, size_t nodeIndex) const;
356 
357 
358  std::vector<Vec3R>& mInstancePoints;
359  std::vector<float>& mInstanceDistances;
360 
361  const PointList& mSurfacePointList;
362 
363  const std::vector<const Index32LeafT*>& mLeafNodes;
364  const std::vector<IndexRange>& mLeafRanges;
365  const std::vector<Vec4R>& mLeafBoundingSpheres;
366  const std::vector<Vec4R>& mNodeBoundingSpheres;
367 
368  std::vector<float> mLeafDistances, mNodeDistances;
369 
370  const bool mTransformPoints;
371  size_t mClosestPointIndex;
372 };// ClosestPointDist
373 
374 
375 template<typename Index32LeafT>
376 ClosestPointDist<Index32LeafT>::ClosestPointDist(
377  std::vector<Vec3R>& instancePoints,
378  std::vector<float>& instanceDistances,
379  const PointList& surfacePointList,
380  const std::vector<const Index32LeafT*>& leafNodes,
381  const std::vector<IndexRange>& leafRanges,
382  const std::vector<Vec4R>& leafBoundingSpheres,
383  const std::vector<Vec4R>& nodeBoundingSpheres,
384  size_t maxNodeLeafs,
385  bool transformPoints)
386  : mInstancePoints(instancePoints)
387  , mInstanceDistances(instanceDistances)
388  , mSurfacePointList(surfacePointList)
389  , mLeafNodes(leafNodes)
390  , mLeafRanges(leafRanges)
391  , mLeafBoundingSpheres(leafBoundingSpheres)
392  , mNodeBoundingSpheres(nodeBoundingSpheres)
393  , mLeafDistances(maxNodeLeafs, 0.0)
394  , mNodeDistances(leafRanges.size(), 0.0)
395  , mTransformPoints(transformPoints)
396  , mClosestPointIndex(0)
397 {
398 }
399 
400 
401 template<typename Index32LeafT>
402 void
403 ClosestPointDist<Index32LeafT>::run(bool threaded)
404 {
405  if (threaded) {
406  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
407  } else {
408  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
409  }
410 }
411 
412 template<typename Index32LeafT>
413 void
414 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
415 {
416  typename Index32LeafT::ValueOnCIter iter;
417  const Vec3s center = mInstancePoints[index];
418  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
419 
420  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
421 
422  const Vec3s& point = mSurfacePointList[iter.getValue()];
423  float tmpDist = (point - center).lengthSqr();
424 
425  if (tmpDist < mInstanceDistances[index]) {
426  mInstanceDistances[index] = tmpDist;
427  closestPointIndex = iter.getValue();
428  }
429  }
430 }
431 
432 
433 template<typename Index32LeafT>
434 void
435 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
436 {
437  if (nodeIndex >= mLeafRanges.size()) return;
438 
439  const Vec3R& pos = mInstancePoints[pointIndex];
440  float minDist = mInstanceDistances[pointIndex];
441  size_t minDistIdx = 0;
442  Vec3R center;
443  bool updatedDist = false;
444 
445  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
446  i < mLeafRanges[nodeIndex].second; ++i, ++n)
447  {
448  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
449 
450  center[0] = mLeafBoundingSpheres[i][0];
451  center[1] = mLeafBoundingSpheres[i][1];
452  center[2] = mLeafBoundingSpheres[i][2];
453  const auto radiusSqr = mLeafBoundingSpheres[i][3];
454 
455  distToLeaf = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
456 
457  if (distToLeaf < minDist) {
458  minDist = distToLeaf;
459  minDistIdx = i;
460  updatedDist = true;
461  }
462  }
463 
464  if (!updatedDist) return;
465 
466  evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
467 
468  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
469  i < mLeafRanges[nodeIndex].second; ++i, ++n)
470  {
471  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
472  evalLeaf(pointIndex, *mLeafNodes[i]);
473  }
474  }
475 }
476 
477 
478 template<typename Index32LeafT>
479 void
480 ClosestPointDist<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
481 {
482  Vec3R center;
483  for (size_t n = range.begin(); n != range.end(); ++n) {
484 
485  const Vec3R& pos = mInstancePoints[n];
486  float minDist = mInstanceDistances[n];
487  size_t minDistIdx = 0;
488 
489  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
490  float& distToNode = const_cast<float&>(mNodeDistances[i]);
491 
492  center[0] = mNodeBoundingSpheres[i][0];
493  center[1] = mNodeBoundingSpheres[i][1];
494  center[2] = mNodeBoundingSpheres[i][2];
495  const auto radiusSqr = mNodeBoundingSpheres[i][3];
496 
497  distToNode = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
498 
499  if (distToNode < minDist) {
500  minDist = distToNode;
501  minDistIdx = i;
502  }
503  }
504 
505  evalNode(n, minDistIdx);
506 
507  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
508  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
509  evalNode(n, i);
510  }
511  }
512 
513  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
514 
515  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
516  }
517 }
518 
519 
520 class UpdatePoints
521 {
522 public:
523  UpdatePoints(
524  const Vec4s& sphere,
525  const std::vector<Vec3R>& points,
526  std::vector<float>& distances,
527  std::vector<unsigned char>& mask,
528  bool overlapping);
529 
530  float radius() const { return mRadius; }
531  int index() const { return mIndex; }
532 
533  inline void run(bool threaded = true);
534 
535 
536  UpdatePoints(UpdatePoints&, tbb::split);
537  inline void operator()(const tbb::blocked_range<size_t>& range);
538  void join(const UpdatePoints& rhs)
539  {
540  if (rhs.mRadius > mRadius) {
541  mRadius = rhs.mRadius;
542  mIndex = rhs.mIndex;
543  }
544  }
545 
546 private:
547  const Vec4s& mSphere;
548  const std::vector<Vec3R>& mPoints;
549  std::vector<float>& mDistances;
550  std::vector<unsigned char>& mMask;
551  bool mOverlapping;
552  float mRadius;
553  int mIndex;
554 };
555 
556 inline
557 UpdatePoints::UpdatePoints(
558  const Vec4s& sphere,
559  const std::vector<Vec3R>& points,
560  std::vector<float>& distances,
561  std::vector<unsigned char>& mask,
562  bool overlapping)
563  : mSphere(sphere)
564  , mPoints(points)
565  , mDistances(distances)
566  , mMask(mask)
567  , mOverlapping(overlapping)
568  , mRadius(0.0)
569  , mIndex(0)
570 {
571 }
572 
573 inline
574 UpdatePoints::UpdatePoints(UpdatePoints& rhs, tbb::split)
575  : mSphere(rhs.mSphere)
576  , mPoints(rhs.mPoints)
577  , mDistances(rhs.mDistances)
578  , mMask(rhs.mMask)
579  , mOverlapping(rhs.mOverlapping)
580  , mRadius(rhs.mRadius)
581  , mIndex(rhs.mIndex)
582 {
583 }
584 
585 inline void
586 UpdatePoints::run(bool threaded)
587 {
588  if (threaded) {
589  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
590  } else {
591  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
592  }
593 }
594 
595 inline void
596 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
597 {
598  Vec3s pos;
599  for (size_t n = range.begin(); n != range.end(); ++n) {
600  if (mMask[n]) continue;
601 
602  pos.x() = float(mPoints[n].x()) - mSphere[0];
603  pos.y() = float(mPoints[n].y()) - mSphere[1];
604  pos.z() = float(mPoints[n].z()) - mSphere[2];
605 
606  float dist = pos.length();
607 
608  if (dist < mSphere[3]) {
609  mMask[n] = 1;
610  continue;
611  }
612 
613  if (!mOverlapping) {
614  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
615  }
616 
617  if (mDistances[n] > mRadius) {
618  mRadius = mDistances[n];
619  mIndex = int(n);
620  }
621  }
622 }
623 
624 
625 } // namespace v2s_internal
626 
627 /// @endcond
628 
629 ////////////////////////////////////////
630 
631 
632 template<typename GridT, typename InterrupterT>
633 inline void
635  const GridT& grid,
636  std::vector<openvdb::Vec4s>& spheres,
637  const Vec2i& sphereCount,
638  bool overlapping,
639  float minRadius,
640  float maxRadius,
641  float isovalue,
642  int instanceCount,
643  InterrupterT* interrupter)
644 {
645  spheres.clear();
646 
647  if (grid.empty()) return;
648 
649  const int
650  minSphereCount = sphereCount[0],
651  maxSphereCount = sphereCount[1];
652  if ((minSphereCount > maxSphereCount) || (maxSphereCount < 1)) {
653  OPENVDB_LOG_WARN("fillWithSpheres: minimum sphere count ("
654  << minSphereCount << ") exceeds maximum count (" << maxSphereCount << ")");
655  return;
656  }
657  spheres.reserve(maxSphereCount);
658 
659  auto gridPtr = grid.copy(); // shallow copy
660 
661  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
662  // Clamp the isovalue to the level set's background value minus epsilon.
663  // (In a valid narrow-band level set, all voxels, including background voxels,
664  // have values less than or equal to the background value, so an isovalue
665  // greater than or equal to the background value would produce a mask with
666  // effectively infinite extent.)
667  isovalue = std::min(isovalue,
668  static_cast<float>(gridPtr->background() - math::Tolerance<float>::value()));
669  } else if (gridPtr->getGridClass() == GRID_FOG_VOLUME) {
670  // Clamp the isovalue of a fog volume between epsilon and one,
671  // again to avoid a mask with infinite extent. (Recall that
672  // fog volume voxel values vary from zero outside to one inside.)
673  isovalue = math::Clamp(isovalue, math::Tolerance<float>::value(), 1.f);
674  }
675 
676  // ClosestSurfacePoint is inaccurate for small grids.
677  // Resample the input grid if it is too small.
678  auto numVoxels = gridPtr->activeVoxelCount();
679  if (numVoxels < 10000) {
680  const auto scale = 1.0 / math::Cbrt(2.0 * 10000.0 / double(numVoxels));
681  auto scaledXform = gridPtr->transform().copy();
682  scaledXform->preScale(scale);
683 
684  auto newGridPtr = levelSetRebuild(*gridPtr, isovalue,
685  LEVEL_SET_HALF_WIDTH, LEVEL_SET_HALF_WIDTH, scaledXform.get(), interrupter);
686 
687  const auto newNumVoxels = newGridPtr->activeVoxelCount();
688  if (newNumVoxels > numVoxels) {
689  OPENVDB_LOG_DEBUG_RUNTIME("fillWithSpheres: resampled input grid from "
690  << numVoxels << " voxel" << (numVoxels == 1 ? "" : "s")
691  << " to " << newNumVoxels << " voxel" << (newNumVoxels == 1 ? "" : "s"));
692  gridPtr = newGridPtr;
693  numVoxels = newNumVoxels;
694  }
695  }
696 
697  const bool addNarrowBandPoints = (numVoxels < 10000);
698  int instances = std::max(instanceCount, maxSphereCount);
699 
700  using TreeT = typename GridT::TreeType;
701  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
702  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
703 
704  using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
705  0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>; // mt11213b
706  RandGen mtRand(/*seed=*/0);
707 
708  const TreeT& tree = gridPtr->tree();
709  math::Transform transform = gridPtr->transform();
710 
711  std::vector<Vec3R> instancePoints;
712  {
713  // Compute a mask of the voxels enclosed by the isosurface.
714  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
715  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
716  interiorMaskPtr = sdfInteriorMask(*gridPtr, isovalue);
717  } else {
718  // For non-level-set grids, the interior mask comprises the active voxels.
719  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
720  interiorMaskPtr->setTransform(transform.copy());
721  interiorMaskPtr->tree().topologyUnion(tree);
722  }
723 
724  if (interrupter && interrupter->wasInterrupted()) return;
725 
726  // If the interior mask is small and eroding it results in an empty grid,
727  // use the uneroded mask instead. (But if the minimum sphere count is zero,
728  // then eroding away the mask is acceptable.)
729  if (!addNarrowBandPoints || (minSphereCount <= 0)) {
730  tools::erodeActiveValues(interiorMaskPtr->tree(), /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
731  tools::pruneInactive(interiorMaskPtr->tree());
732  } else {
733  auto& maskTree = interiorMaskPtr->tree();
734  auto copyOfTree = StaticPtrCast<BoolTreeT>(maskTree.copy());
735  tools::erodeActiveValues(maskTree, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
736  tools::pruneInactive(maskTree);
737  if (maskTree.empty()) { interiorMaskPtr->setTree(copyOfTree); }
738  }
739 
740  // Scatter candidate sphere centroids (instancePoints)
741  instancePoints.reserve(instances);
742  v2s_internal::PointAccessor ptnAcc(instancePoints);
743 
744  const auto scatterCount = Index64(addNarrowBandPoints ? (instances / 2) : instances);
745 
747  ptnAcc, scatterCount, mtRand, 1.0, interrupter);
748  scatter(*interiorMaskPtr);
749  }
750 
751  if (interrupter && interrupter->wasInterrupted()) return;
752 
753  auto csp = ClosestSurfacePoint<GridT>::create(*gridPtr, isovalue, interrupter);
754  if (!csp) return;
755 
756  // Add extra instance points in the interior narrow band.
757  if (instancePoints.size() < size_t(instances)) {
758  const Int16TreeT& signTree = csp->signTree();
759  for (auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
760  for (auto it = leafIt->cbeginValueOn(); it; ++it) {
761  const int flags = int(it.getValue());
762  if (!(volume_to_mesh_internal::EDGES & flags)
763  && (volume_to_mesh_internal::INSIDE & flags))
764  {
765  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
766  }
767  if (instancePoints.size() == size_t(instances)) break;
768  }
769  if (instancePoints.size() == size_t(instances)) break;
770  }
771  }
772 
773  if (interrupter && interrupter->wasInterrupted()) return;
774 
775  // Assign a radius to each candidate sphere. The radius is the world-space
776  // distance from the sphere's center to the closest surface point.
777  std::vector<float> instanceRadius;
778  if (!csp->search(instancePoints, instanceRadius)) return;
779 
780  float largestRadius = 0.0;
781  int largestRadiusIdx = 0;
782  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
783  if (instanceRadius[n] > largestRadius) {
784  largestRadius = instanceRadius[n];
785  largestRadiusIdx = int(n);
786  }
787  }
788 
789  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
790 
791  minRadius = float(minRadius * transform.voxelSize()[0]);
792  maxRadius = float(maxRadius * transform.voxelSize()[0]);
793 
794  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
795 
796  if (interrupter && interrupter->wasInterrupted()) return;
797 
798  largestRadius = std::min(maxRadius, largestRadius);
799 
800  if ((int(s) >= minSphereCount) && (largestRadius < minRadius)) break;
801 
802  const Vec4s sphere(
803  float(instancePoints[largestRadiusIdx].x()),
804  float(instancePoints[largestRadiusIdx].y()),
805  float(instancePoints[largestRadiusIdx].z()),
806  largestRadius);
807 
808  spheres.push_back(sphere);
809  instanceMask[largestRadiusIdx] = 1;
810 
811  v2s_internal::UpdatePoints op(
812  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
813  op.run();
814 
815  largestRadius = op.radius();
816  largestRadiusIdx = op.index();
817  }
818 } // fillWithSpheres
819 
820 
821 ////////////////////////////////////////
822 
823 
824 template<typename GridT>
825 template<typename InterrupterT>
826 inline typename ClosestSurfacePoint<GridT>::Ptr
827 ClosestSurfacePoint<GridT>::create(const GridT& grid, float isovalue, InterrupterT* interrupter)
828 {
829  auto csp = Ptr{new ClosestSurfacePoint};
830  if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
831  return csp;
832 }
833 
834 
835 template<typename GridT>
836 template<typename InterrupterT>
837 inline bool
839  const GridT& grid, float isovalue, InterrupterT* interrupter)
840 {
841  using Index32LeafManagerT = tree::LeafManager<Index32TreeT>;
842  using ValueT = typename GridT::ValueType;
843 
844  const TreeT& tree = grid.tree();
845  const math::Transform& transform = grid.transform();
846 
847  { // Extract surface point cloud
848 
849  BoolTreeT mask(false);
850  volume_to_mesh_internal::identifySurfaceIntersectingVoxels(mask, tree, ValueT(isovalue));
851 
852  mSignTreePt.reset(new Int16TreeT(0));
853  mIdxTreePt.reset(new Index32TreeT(std::numeric_limits<Index32>::max()));
854 
855 
856  volume_to_mesh_internal::computeAuxiliaryData(
857  *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
858 
859  if (interrupter && interrupter->wasInterrupted()) return false;
860 
861  // count unique points
862 
863  using Int16LeafNodeType = typename Int16TreeT::LeafNodeType;
864  using Index32LeafNodeType = typename Index32TreeT::LeafNodeType;
865 
866  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
867  mSignTreePt->getNodes(signFlagsLeafNodes);
868 
869  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
870 
871  std::unique_ptr<Index32[]> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
872 
873  tbb::parallel_for(auxiliaryLeafNodeRange,
874  volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
875  (signFlagsLeafNodes, leafNodeOffsets));
876 
877  {
878  Index32 pointCount = 0;
879  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
880  const Index32 tmp = leafNodeOffsets[n];
881  leafNodeOffsets[n] = pointCount;
882  pointCount += tmp;
883  }
884 
885  mPointListSize = size_t(pointCount);
886  mSurfacePointList.reset(new Vec3s[mPointListSize]);
887  }
888 
889 
890  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
891  mIdxTreePt->getNodes(pointIndexLeafNodes);
892 
893  tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
894  mSurfacePointList.get(), tree, pointIndexLeafNodes,
895  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
896  }
897 
898  if (interrupter && interrupter->wasInterrupted()) return false;
899 
900  Index32LeafManagerT idxLeafs(*mIdxTreePt);
901 
902  using Index32RootNodeT = typename Index32TreeT::RootNodeType;
903  using Index32NodeChainT = typename Index32RootNodeT::NodeChainType;
904  static_assert(Index32NodeChainT::Size > 1,
905  "expected tree depth greater than one");
906  using Index32InternalNodeT = typename Index32NodeChainT::template Get<1>;
907 
908  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
909  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
910  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
911 
912  std::vector<const Index32InternalNodeT*> internalNodes;
913 
914  const Index32InternalNodeT* node = nullptr;
915  for (; nIt; ++nIt) {
916  nIt.getNode(node);
917  if (node) internalNodes.push_back(node);
918  }
919 
920  std::vector<IndexRange>().swap(mLeafRanges);
921  mLeafRanges.resize(internalNodes.size());
922 
923  std::vector<const Index32LeafT*>().swap(mLeafNodes);
924  mLeafNodes.reserve(idxLeafs.leafCount());
925 
926  typename Index32InternalNodeT::ChildOnCIter leafIt;
927  mMaxNodeLeafs = 0;
928  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
929 
930  mLeafRanges[n].first = mLeafNodes.size();
931 
932  size_t leafCount = 0;
933  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
934  mLeafNodes.push_back(&(*leafIt));
935  ++leafCount;
936  }
937 
938  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
939 
940  mLeafRanges[n].second = mLeafNodes.size();
941  }
942 
943  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
944  mLeafBoundingSpheres.resize(mLeafNodes.size());
945 
946  v2s_internal::LeafOp<Index32LeafT> leafBS(
947  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
948  leafBS.run();
949 
950 
951  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
952  mNodeBoundingSpheres.resize(internalNodes.size());
953 
954  v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
955  nodeBS.run();
956  return true;
957 } // ClosestSurfacePoint::initialize
958 
959 
960 template<typename GridT>
961 inline bool
962 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
963  std::vector<float>& distances, bool transformPoints)
964 {
965  distances.clear();
966  distances.resize(points.size(), std::numeric_limits<float>::infinity());
967 
968  v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
969  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
970  mMaxNodeLeafs, transformPoints);
971 
972  cpd.run();
973 
974  return true;
975 }
976 
977 
978 template<typename GridT>
979 inline bool
980 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
981 {
982  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
983 }
984 
985 
986 template<typename GridT>
987 inline bool
989  std::vector<float>& distances)
990 {
991  return search(points, distances, true);
992 }
993 
994 } // namespace tools
995 } // namespace OPENVDB_VERSION_NAME
996 } // namespace openvdb
997 
998 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
std::unique_ptr< openvdb::Vec3s[]> PointList
Point and primitive list types.
Definition: VolumeToMesh.h:150
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Overwrite each input point with its closest surface point.
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
GA_API const UT_StringHolder dist
static Ptr create(const GridT &grid, float isovalue=0.0, InterrupterT *interrupter=nullptr)
Extract surface points and construct a spatial acceleration structure.
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:108
#define OPENVDB_LOG_WARN(mesg)
Definition: logging.h:277
GLenum GLenum GLenum GLenum GLenum scale
Definition: glew.h:14163
math::Vec3< Real > Vec3R
Definition: Types.h:68
GLint first
Definition: glcorearb.h:404
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:302
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7458
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:917
void swap(T &lhs, T &rhs)
Definition: pugixml.cpp:7172
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:178
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:1130
Tolerance for floating-point comparison.
Definition: Math.h:147
float Cbrt(float x)
Return the cube root of a floating-point value.
Definition: Math.h:772
GLint GLenum GLint x
Definition: glcorearb.h:408
typename TreeT::template ValueConverter< bool >::Type BoolTreeT
GLsizeiptr size
Definition: glcorearb.h:663
GLuint GLenum GLenum transform
Definition: glew.h:15055
Accelerated closest surface point queries for narrow band level sets.
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Compute the distance from each input point to its closest surface point.
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:93
void OIIO_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
exint length() const
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:87
Extract polygonal surfaces from scalar volumes.
GLint GLuint mask
Definition: glcorearb.h:123
typename TreeT::template ValueConverter< Index32 >::Type Index32TreeT
We offer three different algorithms (each in its own class) for scattering of points in active voxels...
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
GLint GLsizei count
Definition: glcorearb.h:404
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
OPENVDB_API void initialize()
Global registration of basic types.
Definition: logging.h:294
GLbitfield flags
Definition: glcorearb.h:1595
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
math::Vec4< Real > Vec4R
Definition: Types.h:83
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:260
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glew.h:3460
GLdouble n
Definition: glcorearb.h:2007
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...
#define OPENVDB_LOG_DEBUG_RUNTIME(mesg)
Definition: logging.h:281
const Index32TreeT & indexTree() const
Tree accessor.
const Int16TreeT & signTree() const
Tree accessor.
GLuint index
Definition: glcorearb.h:785
GA_API const UT_StringHolder N
Implementation of morphological dilation and erosion.
GLfloat f
Definition: glcorearb.h:1925
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...
GLenum GLint * range
Definition: glcorearb.h:1924
typename TreeT::template ValueConverter< Int16 >::Type Int16TreeT
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
ImageBuf OIIO_API add(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1247
GLdouble s
Definition: glew.h:1395
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:114
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.
GLint y
Definition: glcorearb.h:102
arg_join< It, Sentinel, char > join(It begin, Sentinel end, string_view sep)
Definition: format.h:3681
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 setTree(TreeBase::Ptr) override
Associate the given tree with this grid, in place of its existing tree.
Definition: Grid.h:1473