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