HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VolumeToSpheres.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
29 ///////////////////////////////////////////////////////////////////////////
30 
31 /// @file tools/VolumeToSpheres.h
32 ///
33 /// @brief Fill a closed level set or fog volume with adaptively-sized spheres.
34 
35 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
37 
39 #include <openvdb/math/Math.h>
40 #include "Morphology.h" // for erodeVoxels()
41 #include "PointScatter.h"
42 #include "LevelSetRebuild.h"
43 #include "LevelSetUtil.h"
44 #include "VolumeToMesh.h"
45 
46 #include <hboost/mpl/at.hpp>
47 #include <hboost/mpl/int.hpp>
48 #include <hboost/scoped_array.hpp>
49 #include <tbb/blocked_range.h>
50 #include <tbb/parallel_for.h>
51 #include <tbb/parallel_reduce.h>
52 
53 #include <algorithm> // for std::min(), std::max()
54 #include <cmath> // for std::sqrt()
55 #include <limits> // for std::numeric_limits
56 #include <memory>
57 #include <random>
58 #include <utility> // for std::pair
59 #include <vector>
60 
61 
62 namespace openvdb {
64 namespace OPENVDB_VERSION_NAME {
65 namespace tools {
66 
67 /// @brief Fill a closed level set or fog volume with adaptively-sized spheres.
68 ///
69 /// @param grid a scalar grid that defines the surface to be filled with spheres
70 /// @param spheres an output array of 4-tuples representing the fitted spheres<BR>
71 /// The first three components of each tuple specify the sphere center,
72 /// and the fourth specifies the radius.
73 /// The spheres are ordered by radius, from largest to smallest.
74 /// @param sphereCount lower and upper bounds on the number of spheres to be generated<BR>
75 /// The actual number will be somewhere within the bounds.
76 /// @param overlapping toggle to allow spheres to overlap/intersect
77 /// @param minRadius the smallest allowable sphere size, in voxel units<BR>
78 /// @param maxRadius the largest allowable sphere size, in voxel units
79 /// @param isovalue the voxel value that determines the surface of the volume<BR>
80 /// The default value of zero works for signed distance fields,
81 /// while fog volumes require a larger positive value
82 /// (0.5 is a good initial guess).
83 /// @param instanceCount the number of interior points to consider for the sphere placement<BR>
84 /// Increasing this count increases the chances of finding optimal
85 /// sphere sizes.
86 /// @param interrupter pointer to an object adhering to the util::NullInterrupter interface
87 ///
88 /// @note The minimum sphere count takes precedence over the minimum radius.
89 template<typename GridT, typename InterrupterT = util::NullInterrupter>
90 inline void
92  const GridT& grid,
93  std::vector<openvdb::Vec4s>& spheres,
94  const Vec2i& sphereCount = Vec2i(1, 50),
95  bool overlapping = false,
96  float minRadius = 1.0,
97  float maxRadius = std::numeric_limits<float>::max(),
98  float isovalue = 0.0,
99  int instanceCount = 10000,
100  InterrupterT* interrupter = nullptr);
101 
102 
103 /// @deprecated Use the @a sphereCount overload instead.
104 template<typename GridT, typename InterrupterT = util::NullInterrupter>
106 inline void
108  const GridT& grid,
109  std::vector<openvdb::Vec4s>& spheres,
110  int maxSphereCount,
111  bool overlapping = false,
112  float minRadius = 1.0,
113  float maxRadius = std::numeric_limits<float>::max(),
114  float isovalue = 0.0,
115  int instanceCount = 10000,
116  InterrupterT* interrupter = nullptr);
117 
118 
119 ////////////////////////////////////////
120 
121 
122 /// @brief Accelerated closest surface point queries for narrow band level sets
123 /// @details Supports queries that originate at arbitrary world-space locations,
124 /// is not confined to the narrow band region of the input volume geometry.
125 template<typename GridT>
127 {
128 public:
129  using Ptr = std::unique_ptr<ClosestSurfacePoint>;
130  using TreeT = typename GridT::TreeType;
131  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
132  using Index32TreeT = typename TreeT::template ValueConverter<Index32>::Type;
133  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
134 
135  /// @brief Extract surface points and construct a spatial acceleration structure.
136  ///
137  /// @return a null pointer if the initialization fails for any reason,
138  /// otherwise a unique pointer to a newly-allocated ClosestSurfacePoint object.
139  ///
140  /// @param grid a scalar level set or fog volume
141  /// @param isovalue the voxel value that determines the surface of the volume
142  /// The default value of zero works for signed distance fields,
143  /// while fog volumes require a larger positive value
144  /// (0.5 is a good initial guess).
145  /// @param interrupter pointer to an object adhering to the util::NullInterrupter interface.
146  template<typename InterrupterT = util::NullInterrupter>
147  static inline Ptr create(const GridT& grid, float isovalue = 0.0,
148  InterrupterT* interrupter = nullptr);
149 
150  /// @brief Compute the distance from each input point to its closest surface point.
151  /// @param points input list of points in world space
152  /// @param distances output list of closest surface point distances
153  inline bool search(const std::vector<Vec3R>& points, std::vector<float>& distances);
154 
155  /// @brief Overwrite each input point with its closest surface point.
156  /// @param points input/output list of points in world space
157  /// @param distances output list of closest surface point distances
158  inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
159 
160  /// @brief Tree accessor
161  const Index32TreeT& indexTree() const { return *mIdxTreePt; }
162  /// @brief Tree accessor
163  const Int16TreeT& signTree() const { return *mSignTreePt; }
164 
165 private:
166  using Index32LeafT = typename Index32TreeT::LeafNodeType;
167  using IndexRange = std::pair<size_t, size_t>;
168 
169  std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
170  std::vector<IndexRange> mLeafRanges;
171  std::vector<const Index32LeafT*> mLeafNodes;
172  PointList mSurfacePointList;
173  size_t mPointListSize = 0, mMaxNodeLeafs = 0;
174  typename Index32TreeT::Ptr mIdxTreePt;
175  typename Int16TreeT::Ptr mSignTreePt;
176 
177  ClosestSurfacePoint() = default;
178  template<typename InterrupterT = util::NullInterrupter>
179  inline bool initialize(const GridT&, float isovalue, InterrupterT*);
180  inline bool search(std::vector<Vec3R>&, std::vector<float>&, bool transformPoints);
181 };
182 
183 
184 ////////////////////////////////////////
185 
186 
187 // Internal utility methods
188 
189 namespace v2s_internal {
190 
192 {
193  PointAccessor(std::vector<Vec3R>& points)
194  : mPoints(points)
195  {
196  }
197 
198  void add(const Vec3R &pos)
199  {
200  mPoints.push_back(pos);
201  }
202 private:
203  std::vector<Vec3R>& mPoints;
204 };
205 
206 
207 template<typename Index32LeafT>
208 class LeafOp
209 {
210 public:
211 
212  LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
213  const std::vector<const Index32LeafT*>& leafNodes,
214  const math::Transform& transform,
215  const PointList& surfacePointList);
216 
217  void run(bool threaded = true);
218 
219 
220  void operator()(const tbb::blocked_range<size_t>&) const;
221 
222 private:
223  std::vector<Vec4R>& mLeafBoundingSpheres;
224  const std::vector<const Index32LeafT*>& mLeafNodes;
225  const math::Transform& mTransform;
226  const PointList& mSurfacePointList;
227 };
228 
229 template<typename Index32LeafT>
231  std::vector<Vec4R>& leafBoundingSpheres,
232  const std::vector<const Index32LeafT*>& leafNodes,
233  const math::Transform& transform,
234  const PointList& surfacePointList)
235  : mLeafBoundingSpheres(leafBoundingSpheres)
236  , mLeafNodes(leafNodes)
237  , mTransform(transform)
238  , mSurfacePointList(surfacePointList)
239 {
240 }
241 
242 template<typename Index32LeafT>
243 void
245 {
246  if (threaded) {
247  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *this);
248  } else {
249  (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
250  }
251 }
252 
253 template<typename Index32LeafT>
254 void
255 LeafOp<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
256 {
257  typename Index32LeafT::ValueOnCIter iter;
258  Vec3s avg;
259 
260  for (size_t n = range.begin(); n != range.end(); ++n) {
261  avg[0] = 0.0;
262  avg[1] = 0.0;
263  avg[2] = 0.0;
264 
265  int count = 0;
266  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
267  avg += mSurfacePointList[iter.getValue()];
268  ++count;
269  }
270  if (count > 1) avg *= float(1.0 / double(count));
271 
272  float maxDist = 0.0;
273  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
274  float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
275  if (tmpDist > maxDist) maxDist = tmpDist;
276  }
277 
278  Vec4R& sphere = mLeafBoundingSpheres[n];
279  sphere[0] = avg[0];
280  sphere[1] = avg[1];
281  sphere[2] = avg[2];
282  sphere[3] = maxDist * 2.0; // padded radius
283  }
284 }
285 
286 
287 class NodeOp
288 {
289 public:
290  using IndexRange = std::pair<size_t, size_t>;
291 
292  NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
293  const std::vector<IndexRange>& leafRanges,
294  const std::vector<Vec4R>& leafBoundingSpheres);
295 
296  inline void run(bool threaded = true);
297 
298  inline void operator()(const tbb::blocked_range<size_t>&) const;
299 
300 private:
301  std::vector<Vec4R>& mNodeBoundingSpheres;
302  const std::vector<IndexRange>& mLeafRanges;
303  const std::vector<Vec4R>& mLeafBoundingSpheres;
304 };
305 
306 inline
307 NodeOp::NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
308  const std::vector<IndexRange>& leafRanges,
309  const std::vector<Vec4R>& leafBoundingSpheres)
310  : mNodeBoundingSpheres(nodeBoundingSpheres)
311  , mLeafRanges(leafRanges)
312  , mLeafBoundingSpheres(leafBoundingSpheres)
313 {
314 }
315 
316 inline void
317 NodeOp::run(bool threaded)
318 {
319  if (threaded) {
320  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *this);
321  } else {
322  (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
323  }
324 }
325 
326 inline void
327 NodeOp::operator()(const tbb::blocked_range<size_t>& range) const
328 {
329  Vec3d avg, pos;
330 
331  for (size_t n = range.begin(); n != range.end(); ++n) {
332 
333  avg[0] = 0.0;
334  avg[1] = 0.0;
335  avg[2] = 0.0;
336 
337  int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
338 
339  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
340  avg[0] += mLeafBoundingSpheres[i][0];
341  avg[1] += mLeafBoundingSpheres[i][1];
342  avg[2] += mLeafBoundingSpheres[i][2];
343  }
344 
345  if (count > 1) avg *= float(1.0 / double(count));
346 
347 
348  double maxDist = 0.0;
349 
350  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
351  pos[0] = mLeafBoundingSpheres[i][0];
352  pos[1] = mLeafBoundingSpheres[i][1];
353  pos[2] = mLeafBoundingSpheres[i][2];
354  const auto radiusSqr = mLeafBoundingSpheres[i][3];
355 
356  double tmpDist = (pos - avg).lengthSqr() + radiusSqr;
357  if (tmpDist > maxDist) maxDist = tmpDist;
358  }
359 
360  Vec4R& sphere = mNodeBoundingSpheres[n];
361 
362  sphere[0] = avg[0];
363  sphere[1] = avg[1];
364  sphere[2] = avg[2];
365  sphere[3] = maxDist * 2.0; // padded radius
366  }
367 }
368 
369 
370 ////////////////////////////////////////
371 
372 
373 template<typename Index32LeafT>
375 {
376 public:
377  using IndexRange = std::pair<size_t, size_t>;
378 
380  std::vector<Vec3R>& instancePoints,
381  std::vector<float>& instanceDistances,
382  const PointList& surfacePointList,
383  const std::vector<const Index32LeafT*>& leafNodes,
384  const std::vector<IndexRange>& leafRanges,
385  const std::vector<Vec4R>& leafBoundingSpheres,
386  const std::vector<Vec4R>& nodeBoundingSpheres,
387  size_t maxNodeLeafs,
388  bool transformPoints = false);
389 
390 
391  void run(bool threaded = true);
392 
393 
394  void operator()(const tbb::blocked_range<size_t>&) const;
395 
396 private:
397 
398  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
399  void evalNode(size_t pointIndex, size_t nodeIndex) const;
400 
401 
402  std::vector<Vec3R>& mInstancePoints;
403  std::vector<float>& mInstanceDistances;
404 
405  const PointList& mSurfacePointList;
406 
407  const std::vector<const Index32LeafT*>& mLeafNodes;
408  const std::vector<IndexRange>& mLeafRanges;
409  const std::vector<Vec4R>& mLeafBoundingSpheres;
410  const std::vector<Vec4R>& mNodeBoundingSpheres;
411 
412  std::vector<float> mLeafDistances, mNodeDistances;
413 
414  const bool mTransformPoints;
415  size_t mClosestPointIndex;
416 };// ClosestPointDist
417 
418 
419 template<typename Index32LeafT>
421  std::vector<Vec3R>& instancePoints,
422  std::vector<float>& instanceDistances,
423  const PointList& surfacePointList,
424  const std::vector<const Index32LeafT*>& leafNodes,
425  const std::vector<IndexRange>& leafRanges,
426  const std::vector<Vec4R>& leafBoundingSpheres,
427  const std::vector<Vec4R>& nodeBoundingSpheres,
428  size_t maxNodeLeafs,
429  bool transformPoints)
430  : mInstancePoints(instancePoints)
431  , mInstanceDistances(instanceDistances)
432  , mSurfacePointList(surfacePointList)
433  , mLeafNodes(leafNodes)
434  , mLeafRanges(leafRanges)
435  , mLeafBoundingSpheres(leafBoundingSpheres)
436  , mNodeBoundingSpheres(nodeBoundingSpheres)
437  , mLeafDistances(maxNodeLeafs, 0.0)
438  , mNodeDistances(leafRanges.size(), 0.0)
439  , mTransformPoints(transformPoints)
440  , mClosestPointIndex(0)
441 {
442 }
443 
444 
445 template<typename Index32LeafT>
446 void
448 {
449  if (threaded) {
450  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
451  } else {
452  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
453  }
454 }
455 
456 template<typename Index32LeafT>
457 void
458 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
459 {
460  typename Index32LeafT::ValueOnCIter iter;
461  const Vec3s center = mInstancePoints[index];
462  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
463 
464  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
465 
466  const Vec3s& point = mSurfacePointList[iter.getValue()];
467  float tmpDist = (point - center).lengthSqr();
468 
469  if (tmpDist < mInstanceDistances[index]) {
470  mInstanceDistances[index] = tmpDist;
471  closestPointIndex = iter.getValue();
472  }
473  }
474 }
475 
476 
477 template<typename Index32LeafT>
478 void
479 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
480 {
481  if (nodeIndex >= mLeafRanges.size()) return;
482 
483  const Vec3R& pos = mInstancePoints[pointIndex];
484  float minDist = mInstanceDistances[pointIndex];
485  size_t minDistIdx = 0;
486  Vec3R center;
487  bool updatedDist = false;
488 
489  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
490  i < mLeafRanges[nodeIndex].second; ++i, ++n)
491  {
492  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
493 
494  center[0] = mLeafBoundingSpheres[i][0];
495  center[1] = mLeafBoundingSpheres[i][1];
496  center[2] = mLeafBoundingSpheres[i][2];
497  const auto radiusSqr = mLeafBoundingSpheres[i][3];
498 
499  distToLeaf = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
500 
501  if (distToLeaf < minDist) {
502  minDist = distToLeaf;
503  minDistIdx = i;
504  updatedDist = true;
505  }
506  }
507 
508  if (!updatedDist) return;
509 
510  evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
511 
512  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
513  i < mLeafRanges[nodeIndex].second; ++i, ++n)
514  {
515  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
516  evalLeaf(pointIndex, *mLeafNodes[i]);
517  }
518  }
519 }
520 
521 
522 template<typename Index32LeafT>
523 void
524 ClosestPointDist<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
525 {
526  Vec3R center;
527  for (size_t n = range.begin(); n != range.end(); ++n) {
528 
529  const Vec3R& pos = mInstancePoints[n];
530  float minDist = mInstanceDistances[n];
531  size_t minDistIdx = 0;
532 
533  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
534  float& distToNode = const_cast<float&>(mNodeDistances[i]);
535 
536  center[0] = mNodeBoundingSpheres[i][0];
537  center[1] = mNodeBoundingSpheres[i][1];
538  center[2] = mNodeBoundingSpheres[i][2];
539  const auto radiusSqr = mNodeBoundingSpheres[i][3];
540 
541  distToNode = float(std::max(0.0, (pos - center).lengthSqr() - radiusSqr));
542 
543  if (distToNode < minDist) {
544  minDist = distToNode;
545  minDistIdx = i;
546  }
547  }
548 
549  evalNode(n, minDistIdx);
550 
551  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
552  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
553  evalNode(n, i);
554  }
555  }
556 
557  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
558 
559  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
560  }
561 }
562 
563 
565 {
566 public:
567  UpdatePoints(
568  const Vec4s& sphere,
569  const std::vector<Vec3R>& points,
570  std::vector<float>& distances,
571  std::vector<unsigned char>& mask,
572  bool overlapping);
573 
574  float radius() const { return mRadius; }
575  int index() const { return mIndex; }
576 
577  inline void run(bool threaded = true);
578 
579 
581  inline void operator()(const tbb::blocked_range<size_t>& range);
582  void join(const UpdatePoints& rhs)
583  {
584  if (rhs.mRadius > mRadius) {
585  mRadius = rhs.mRadius;
586  mIndex = rhs.mIndex;
587  }
588  }
589 
590 private:
591  const Vec4s& mSphere;
592  const std::vector<Vec3R>& mPoints;
593  std::vector<float>& mDistances;
594  std::vector<unsigned char>& mMask;
595  bool mOverlapping;
596  float mRadius;
597  int mIndex;
598 };
599 
600 inline
602  const Vec4s& sphere,
603  const std::vector<Vec3R>& points,
604  std::vector<float>& distances,
605  std::vector<unsigned char>& mask,
606  bool overlapping)
607  : mSphere(sphere)
608  , mPoints(points)
609  , mDistances(distances)
610  , mMask(mask)
611  , mOverlapping(overlapping)
612  , mRadius(0.0)
613  , mIndex(0)
614 {
615 }
616 
617 inline
619  : mSphere(rhs.mSphere)
620  , mPoints(rhs.mPoints)
621  , mDistances(rhs.mDistances)
622  , mMask(rhs.mMask)
623  , mOverlapping(rhs.mOverlapping)
624  , mRadius(rhs.mRadius)
625  , mIndex(rhs.mIndex)
626 {
627 }
628 
629 inline void
630 UpdatePoints::run(bool threaded)
631 {
632  if (threaded) {
633  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
634  } else {
635  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
636  }
637 }
638 
639 inline void
640 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
641 {
642  Vec3s pos;
643  for (size_t n = range.begin(); n != range.end(); ++n) {
644  if (mMask[n]) continue;
645 
646  pos.x() = float(mPoints[n].x()) - mSphere[0];
647  pos.y() = float(mPoints[n].y()) - mSphere[1];
648  pos.z() = float(mPoints[n].z()) - mSphere[2];
649 
650  float dist = pos.length();
651 
652  if (dist < mSphere[3]) {
653  mMask[n] = 1;
654  continue;
655  }
656 
657  if (!mOverlapping) {
658  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
659  }
660 
661  if (mDistances[n] > mRadius) {
662  mRadius = mDistances[n];
663  mIndex = int(n);
664  }
665  }
666 }
667 
668 
669 } // namespace v2s_internal
670 
671 
672 ////////////////////////////////////////
673 
674 
675 template<typename GridT, typename InterrupterT>
676 inline void
678  const GridT& grid,
679  std::vector<openvdb::Vec4s>& spheres,
680  int maxSphereCount,
681  bool overlapping,
682  float minRadius,
683  float maxRadius,
684  float isovalue,
685  int instanceCount,
686  InterrupterT* interrupter)
687 {
688  fillWithSpheres(grid, spheres, Vec2i(1, maxSphereCount), overlapping,
689  minRadius, maxRadius, isovalue, instanceCount, interrupter);
690 }
691 
692 
693 template<typename GridT, typename InterrupterT>
694 inline void
696  const GridT& grid,
697  std::vector<openvdb::Vec4s>& spheres,
698  const Vec2i& sphereCount,
699  bool overlapping,
700  float minRadius,
701  float maxRadius,
702  float isovalue,
703  int instanceCount,
704  InterrupterT* interrupter)
705 {
706  spheres.clear();
707 
708  if (grid.empty()) return;
709 
710  const int
711  minSphereCount = sphereCount[0],
712  maxSphereCount = sphereCount[1];
713  if ((minSphereCount > maxSphereCount) || (maxSphereCount < 1)) {
714  OPENVDB_LOG_WARN("fillWithSpheres: minimum sphere count ("
715  << minSphereCount << ") exceeds maximum count (" << maxSphereCount << ")");
716  return;
717  }
718  spheres.reserve(maxSphereCount);
719 
720  auto gridPtr = grid.copy(); // shallow copy
721 
722  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
723  // Clamp the isovalue to the level set's background value minus epsilon.
724  // (In a valid narrow-band level set, all voxels, including background voxels,
725  // have values less than or equal to the background value, so an isovalue
726  // greater than or equal to the background value would produce a mask with
727  // effectively infinite extent.)
728  isovalue = std::min(isovalue,
729  static_cast<float>(gridPtr->background() - math::Tolerance<float>::value()));
730  } else if (gridPtr->getGridClass() == GRID_FOG_VOLUME) {
731  // Clamp the isovalue of a fog volume between epsilon and one,
732  // again to avoid a mask with infinite extent. (Recall that
733  // fog volume voxel values vary from zero outside to one inside.)
734  isovalue = math::Clamp(isovalue, math::Tolerance<float>::value(), 1.f);
735  }
736 
737  // ClosestSurfacePoint is inaccurate for small grids.
738  // Resample the input grid if it is too small.
739  auto numVoxels = gridPtr->activeVoxelCount();
740  if (numVoxels < 10000) {
741  const auto scale = 1.0 / math::Cbrt(2.0 * 10000.0 / double(numVoxels));
742  auto scaledXform = gridPtr->transform().copy();
743  scaledXform->preScale(scale);
744 
745  auto newGridPtr = levelSetRebuild(*gridPtr, isovalue,
746  LEVEL_SET_HALF_WIDTH, LEVEL_SET_HALF_WIDTH, scaledXform.get(), interrupter);
747 
748  const auto newNumVoxels = newGridPtr->activeVoxelCount();
749  if (newNumVoxels > numVoxels) {
750  OPENVDB_LOG_DEBUG_RUNTIME("fillWithSpheres: resampled input grid from "
751  << numVoxels << " voxel" << (numVoxels == 1 ? "" : "s")
752  << " to " << newNumVoxels << " voxel" << (newNumVoxels == 1 ? "" : "s"));
753  gridPtr = newGridPtr;
754  numVoxels = newNumVoxels;
755  }
756  }
757 
758  const bool addNarrowBandPoints = (numVoxels < 10000);
759  int instances = std::max(instanceCount, maxSphereCount);
760 
761  using TreeT = typename GridT::TreeType;
762  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
763  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
764 
765  using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
766  0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>; // mt11213b
767  RandGen mtRand(/*seed=*/0);
768 
769  const TreeT& tree = gridPtr->tree();
770  math::Transform transform = gridPtr->transform();
771 
772  std::vector<Vec3R> instancePoints;
773  {
774  // Compute a mask of the voxels enclosed by the isosurface.
775  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
776  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
777  interiorMaskPtr = sdfInteriorMask(*gridPtr, isovalue);
778  } else {
779  // For non-level-set grids, the interior mask comprises the active voxels.
780  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
781  interiorMaskPtr->setTransform(transform.copy());
782  interiorMaskPtr->tree().topologyUnion(tree);
783  }
784 
785  if (interrupter && interrupter->wasInterrupted()) return;
786 
787  // If the interior mask is small and eroding it results in an empty grid,
788  // use the uneroded mask instead. (But if the minimum sphere count is zero,
789  // then eroding away the mask is acceptable.)
790  if (!addNarrowBandPoints || (minSphereCount <= 0)) {
791  erodeVoxels(interiorMaskPtr->tree(), 1);
792  } else {
793  auto& maskTree = interiorMaskPtr->tree();
794  auto copyOfTree = StaticPtrCast<BoolTreeT>(maskTree.copy());
795  erodeVoxels(maskTree, 1);
796  if (maskTree.empty()) { interiorMaskPtr->setTree(copyOfTree); }
797  }
798 
799  // Scatter candidate sphere centroids (instancePoints)
800  instancePoints.reserve(instances);
801  v2s_internal::PointAccessor ptnAcc(instancePoints);
802 
803  const auto scatterCount = Index64(addNarrowBandPoints ? (instances / 2) : instances);
804 
806  ptnAcc, scatterCount, mtRand, 1.0, interrupter);
807  scatter(*interiorMaskPtr);
808  }
809 
810  if (interrupter && interrupter->wasInterrupted()) return;
811 
812  auto csp = ClosestSurfacePoint<GridT>::create(*gridPtr, isovalue, interrupter);
813  if (!csp) return;
814 
815  // Add extra instance points in the interior narrow band.
816  if (instancePoints.size() < size_t(instances)) {
817  const Int16TreeT& signTree = csp->signTree();
818  for (auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
819  for (auto it = leafIt->cbeginValueOn(); it; ++it) {
820  const int flags = int(it.getValue());
821  if (!(volume_to_mesh_internal::EDGES & flags)
823  {
824  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
825  }
826  if (instancePoints.size() == size_t(instances)) break;
827  }
828  if (instancePoints.size() == size_t(instances)) break;
829  }
830  }
831 
832  if (interrupter && interrupter->wasInterrupted()) return;
833 
834  // Assign a radius to each candidate sphere. The radius is the world-space
835  // distance from the sphere's center to the closest surface point.
836  std::vector<float> instanceRadius;
837  if (!csp->search(instancePoints, instanceRadius)) return;
838 
839  float largestRadius = 0.0;
840  int largestRadiusIdx = 0;
841  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
842  if (instanceRadius[n] > largestRadius) {
843  largestRadius = instanceRadius[n];
844  largestRadiusIdx = int(n);
845  }
846  }
847 
848  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
849 
850  minRadius = float(minRadius * transform.voxelSize()[0]);
851  maxRadius = float(maxRadius * transform.voxelSize()[0]);
852 
853  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
854 
855  if (interrupter && interrupter->wasInterrupted()) return;
856 
857  largestRadius = std::min(maxRadius, largestRadius);
858 
859  if ((int(s) >= minSphereCount) && (largestRadius < minRadius)) break;
860 
861  const Vec4s sphere(
862  float(instancePoints[largestRadiusIdx].x()),
863  float(instancePoints[largestRadiusIdx].y()),
864  float(instancePoints[largestRadiusIdx].z()),
865  largestRadius);
866 
867  spheres.push_back(sphere);
868  instanceMask[largestRadiusIdx] = 1;
869 
871  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
872  op.run();
873 
874  largestRadius = op.radius();
875  largestRadiusIdx = op.index();
876  }
877 } // fillWithSpheres
878 
879 
880 ////////////////////////////////////////
881 
882 
883 template<typename GridT>
884 template<typename InterrupterT>
885 inline typename ClosestSurfacePoint<GridT>::Ptr
886 ClosestSurfacePoint<GridT>::create(const GridT& grid, float isovalue, InterrupterT* interrupter)
887 {
888  auto csp = Ptr{new ClosestSurfacePoint};
889  if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
890  return csp;
891 }
892 
893 
894 template<typename GridT>
895 template<typename InterrupterT>
896 inline bool
898  const GridT& grid, float isovalue, InterrupterT* interrupter)
899 {
900  using Index32LeafManagerT = tree::LeafManager<Index32TreeT>;
901  using ValueT = typename GridT::ValueType;
902 
903  const TreeT& tree = grid.tree();
904  const math::Transform& transform = grid.transform();
905 
906  { // Extract surface point cloud
907 
908  BoolTreeT mask(false);
910 
911  mSignTreePt.reset(new Int16TreeT(0));
912  mIdxTreePt.reset(new Index32TreeT(std::numeric_limits<Index32>::max()));
913 
914 
916  *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
917 
918  if (interrupter && interrupter->wasInterrupted()) return false;
919 
920  // count unique points
921 
922  using Int16LeafNodeType = typename Int16TreeT::LeafNodeType;
923  using Index32LeafNodeType = typename Index32TreeT::LeafNodeType;
924 
925  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
926  mSignTreePt->getNodes(signFlagsLeafNodes);
927 
928  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
929 
930  hboost::scoped_array<Index32> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
931 
932  tbb::parallel_for(auxiliaryLeafNodeRange,
934  (signFlagsLeafNodes, leafNodeOffsets));
935 
936  {
937  Index32 pointCount = 0;
938  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
939  const Index32 tmp = leafNodeOffsets[n];
940  leafNodeOffsets[n] = pointCount;
941  pointCount += tmp;
942  }
943 
944  mPointListSize = size_t(pointCount);
945  mSurfacePointList.reset(new Vec3s[mPointListSize]);
946  }
947 
948 
949  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
950  mIdxTreePt->getNodes(pointIndexLeafNodes);
951 
952  tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
953  mSurfacePointList.get(), tree, pointIndexLeafNodes,
954  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
955  }
956 
957  if (interrupter && interrupter->wasInterrupted()) return false;
958 
959  Index32LeafManagerT idxLeafs(*mIdxTreePt);
960 
961  using Index32RootNodeT = typename Index32TreeT::RootNodeType;
962  using Index32NodeChainT = typename Index32RootNodeT::NodeChainType;
964  "expected tree depth greater than one");
965  using Index32InternalNodeT =
966  typename hboost::mpl::at<Index32NodeChainT, hboost::mpl::int_<1> >::type;
967 
968  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
969  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
970  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
971 
972  std::vector<const Index32InternalNodeT*> internalNodes;
973 
974  const Index32InternalNodeT* node = nullptr;
975  for (; nIt; ++nIt) {
976  nIt.getNode(node);
977  if (node) internalNodes.push_back(node);
978  }
979 
980  std::vector<IndexRange>().swap(mLeafRanges);
981  mLeafRanges.resize(internalNodes.size());
982 
983  std::vector<const Index32LeafT*>().swap(mLeafNodes);
984  mLeafNodes.reserve(idxLeafs.leafCount());
985 
986  typename Index32InternalNodeT::ChildOnCIter leafIt;
987  mMaxNodeLeafs = 0;
988  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
989 
990  mLeafRanges[n].first = mLeafNodes.size();
991 
992  size_t leafCount = 0;
993  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
994  mLeafNodes.push_back(&(*leafIt));
995  ++leafCount;
996  }
997 
998  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
999 
1000  mLeafRanges[n].second = mLeafNodes.size();
1001  }
1002 
1003  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
1004  mLeafBoundingSpheres.resize(mLeafNodes.size());
1005 
1006  v2s_internal::LeafOp<Index32LeafT> leafBS(
1007  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
1008  leafBS.run();
1009 
1010 
1011  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
1012  mNodeBoundingSpheres.resize(internalNodes.size());
1013 
1014  v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
1015  nodeBS.run();
1016  return true;
1017 } // ClosestSurfacePoint::initialize
1018 
1019 
1020 template<typename GridT>
1021 inline bool
1022 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
1023  std::vector<float>& distances, bool transformPoints)
1024 {
1025  distances.clear();
1026  distances.resize(points.size(), std::numeric_limits<float>::infinity());
1027 
1028  v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
1029  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1030  mMaxNodeLeafs, transformPoints);
1031 
1032  cpd.run();
1033 
1034  return true;
1035 }
1036 
1037 
1038 template<typename GridT>
1039 inline bool
1040 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
1041 {
1042  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
1043 }
1044 
1045 
1046 template<typename GridT>
1047 inline bool
1049  std::vector<float>& distances)
1050 {
1051  return search(points, distances, true);
1052 }
1053 
1054 } // namespace tools
1055 } // namespace OPENVDB_VERSION_NAME
1056 } // namespace openvdb
1057 
1058 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
1059 
1060 // Copyright (c) DreamWorks Animation LLC
1061 // All rights reserved. This software is distributed under the
1062 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
GLdouble s
Definition: glew.h:1390
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Overwrite each input point with its closest surface point.
vint4 max(const vint4 &a, const vint4 &b)
Definition: simd.h:4703
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:153
GLenum GLint * range
Definition: glew.h:3500
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.
void operator()(const tbb::blocked_range< size_t > &) const
GLsizeiptr size
Definition: glew.h:1681
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:135
#define OPENVDB_LOG_WARN(mesg)
Definition: logging.h:301
math::Vec3< Real > Vec3R
Definition: Types.h:79
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)
GLuint GLenum GLenum transform
Definition: glew.h:14742
GLuint index
Definition: glew.h:1814
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:342
void operator()(const tbb::blocked_range< size_t > &) const
GLenum GLenum GLenum GLenum GLenum scale
Definition: glew.h:13880
GLenum GLsizei GLsizei GLsizei GLsizei GLbitfield flags
Definition: glew.h:2864
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7231
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:919
void swap(T &lhs, T &rhs)
Definition: pugixml.cpp:7172
const GLint * first
Definition: glew.h:1528
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:200
GLenum GLint GLuint mask
Definition: glew.h:1845
Tolerance for floating-point comparison.
Definition: Math.h:117
float Cbrt(float x)
Return the cube root of a floating-point value.
Definition: Math.h:743
GLdouble GLdouble z
Definition: glew.h:1559
typename TreeT::template ValueConverter< bool >::Type BoolTreeT
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.
LeafOp(std::vector< Vec4R > &leafBoundingSpheres, const std::vector< const Index32LeafT * > &leafNodes, const math::Transform &transform, const PointList &surfacePointList)
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
void OIIO_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
GLclampf f
Definition: glew.h:3499
GLint GLint GLint GLint GLint x
Definition: glew.h:1252
exint length() const
GLint GLint GLint GLint GLint GLint y
Definition: glew.h:1252
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:114
Extract polygonal surfaces from scalar volumes.
#define OPENVDB_DEPRECATED
Definition: Platform.h:69
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...
void operator()(const tbb::blocked_range< size_t > &range)
GLsizei n
Definition: glew.h:4040
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
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:115
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:55
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:230
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...
OPENVDB_STATIC_SPECIALIZATION 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:880
#define OPENVDB_LOG_DEBUG_RUNTIME(mesg)
Definition: logging.h:305
const Index32TreeT & indexTree() const
Tree accessor.
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)
const Int16TreeT & signTree() const
Tree accessor.
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
GLuint GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glew.h:3446
GLuint GLuint GLsizei count
Definition: glew.h:1253
GA_API const UT_StringHolder N
Implementation of morphological dilation and erosion.
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...
typename TreeT::template ValueConverter< Int16 >::Type Int16TreeT
hboost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:179
NodeOp(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
vint4 min(const vint4 &a, const vint4 &b)
Definition: simd.h:4694
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1245
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:146
GLsizei const GLfloat * value
Definition: glew.h:1849
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.
void operator()(const tbb::blocked_range< size_t > &) const
void identifySurfaceIntersectingVoxels(typename InputTreeType::template ValueConverter< bool >::Type &intersectionTree, const InputTreeType &inputTree, typename InputTreeType::ValueType isovalue)
type
Definition: core.h:528
void setTree(TreeBase::Ptr) override
Associate the given tree with this grid, in place of its existing tree.
Definition: Grid.h:1442