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) 2012-2018 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] = std::sqrt(maxDist);
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 
355  double tmpDist = (pos - avg).length() + mLeafBoundingSpheres[i][3];
356  if (tmpDist > maxDist) maxDist = tmpDist;
357  }
358 
359  Vec4R& sphere = mNodeBoundingSpheres[n];
360 
361  sphere[0] = avg[0];
362  sphere[1] = avg[1];
363  sphere[2] = avg[2];
364  sphere[3] = maxDist;
365  }
366 }
367 
368 
369 ////////////////////////////////////////
370 
371 
372 template<typename Index32LeafT>
374 {
375 public:
376  using IndexRange = std::pair<size_t, size_t>;
377 
379  std::vector<Vec3R>& instancePoints,
380  std::vector<float>& instanceDistances,
381  const PointList& surfacePointList,
382  const std::vector<const Index32LeafT*>& leafNodes,
383  const std::vector<IndexRange>& leafRanges,
384  const std::vector<Vec4R>& leafBoundingSpheres,
385  const std::vector<Vec4R>& nodeBoundingSpheres,
386  size_t maxNodeLeafs,
387  bool transformPoints = false);
388 
389 
390  void run(bool threaded = true);
391 
392 
393  void operator()(const tbb::blocked_range<size_t>&) const;
394 
395 private:
396 
397  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
398  void evalNode(size_t pointIndex, size_t nodeIndex) const;
399 
400 
401  std::vector<Vec3R>& mInstancePoints;
402  std::vector<float>& mInstanceDistances;
403 
404  const PointList& mSurfacePointList;
405 
406  const std::vector<const Index32LeafT*>& mLeafNodes;
407  const std::vector<IndexRange>& mLeafRanges;
408  const std::vector<Vec4R>& mLeafBoundingSpheres;
409  const std::vector<Vec4R>& mNodeBoundingSpheres;
410 
411  std::vector<float> mLeafDistances, mNodeDistances;
412 
413  const bool mTransformPoints;
414  size_t mClosestPointIndex;
415 };// ClosestPointDist
416 
417 
418 template<typename Index32LeafT>
420  std::vector<Vec3R>& instancePoints,
421  std::vector<float>& instanceDistances,
422  const PointList& surfacePointList,
423  const std::vector<const Index32LeafT*>& leafNodes,
424  const std::vector<IndexRange>& leafRanges,
425  const std::vector<Vec4R>& leafBoundingSpheres,
426  const std::vector<Vec4R>& nodeBoundingSpheres,
427  size_t maxNodeLeafs,
428  bool transformPoints)
429  : mInstancePoints(instancePoints)
430  , mInstanceDistances(instanceDistances)
431  , mSurfacePointList(surfacePointList)
432  , mLeafNodes(leafNodes)
433  , mLeafRanges(leafRanges)
434  , mLeafBoundingSpheres(leafBoundingSpheres)
435  , mNodeBoundingSpheres(nodeBoundingSpheres)
436  , mLeafDistances(maxNodeLeafs, 0.0)
437  , mNodeDistances(leafRanges.size(), 0.0)
438  , mTransformPoints(transformPoints)
439  , mClosestPointIndex(0)
440 {
441 }
442 
443 
444 template<typename Index32LeafT>
445 void
447 {
448  if (threaded) {
449  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
450  } else {
451  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
452  }
453 }
454 
455 template<typename Index32LeafT>
456 void
457 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
458 {
459  typename Index32LeafT::ValueOnCIter iter;
460  const Vec3s center = mInstancePoints[index];
461  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
462 
463  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
464 
465  const Vec3s& point = mSurfacePointList[iter.getValue()];
466  float tmpDist = (point - center).lengthSqr();
467 
468  if (tmpDist < mInstanceDistances[index]) {
469  mInstanceDistances[index] = tmpDist;
470  closestPointIndex = iter.getValue();
471  }
472  }
473 }
474 
475 
476 template<typename Index32LeafT>
477 void
478 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
479 {
480  if (nodeIndex >= mLeafRanges.size()) return;
481 
482  const Vec3R& pos = mInstancePoints[pointIndex];
483  float minDist = mInstanceDistances[pointIndex];
484  size_t minDistIdx = 0;
485  Vec3R center;
486  bool updatedDist = false;
487 
488  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
489  i < mLeafRanges[nodeIndex].second; ++i, ++n)
490  {
491  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
492 
493  center[0] = mLeafBoundingSpheres[i][0];
494  center[1] = mLeafBoundingSpheres[i][1];
495  center[2] = mLeafBoundingSpheres[i][2];
496  const auto radius = mLeafBoundingSpheres[i][3];
497 
498  distToLeaf = float(std::max(0.0, (pos - center).length() - radius));
499 
500  if (distToLeaf < minDist) {
501  minDist = distToLeaf;
502  minDistIdx = i;
503  updatedDist = true;
504  }
505  }
506 
507  if (!updatedDist) return;
508 
509  evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
510 
511  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
512  i < mLeafRanges[nodeIndex].second; ++i, ++n)
513  {
514  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
515  evalLeaf(pointIndex, *mLeafNodes[i]);
516  }
517  }
518 }
519 
520 
521 template<typename Index32LeafT>
522 void
523 ClosestPointDist<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
524 {
525  Vec3R center;
526  for (size_t n = range.begin(); n != range.end(); ++n) {
527 
528  const Vec3R& pos = mInstancePoints[n];
529  float minDist = mInstanceDistances[n];
530  size_t minDistIdx = 0;
531 
532  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
533  float& distToNode = const_cast<float&>(mNodeDistances[i]);
534 
535  center[0] = mNodeBoundingSpheres[i][0];
536  center[1] = mNodeBoundingSpheres[i][1];
537  center[2] = mNodeBoundingSpheres[i][2];
538  const auto radius = mNodeBoundingSpheres[i][3];
539 
540  distToNode = float(std::max(0.0, (pos - center).length() - radius));
541 
542  if (distToNode < minDist) {
543  minDist = distToNode;
544  minDistIdx = i;
545  }
546  }
547 
548  evalNode(n, minDistIdx);
549 
550  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
551  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
552  evalNode(n, i);
553  }
554  }
555 
556  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
557 
558  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
559  }
560 }
561 
562 
564 {
565 public:
566  UpdatePoints(
567  const Vec4s& sphere,
568  const std::vector<Vec3R>& points,
569  std::vector<float>& distances,
570  std::vector<unsigned char>& mask,
571  bool overlapping);
572 
573  float radius() const { return mRadius; }
574  int index() const { return mIndex; }
575 
576  inline void run(bool threaded = true);
577 
578 
579  UpdatePoints(UpdatePoints&, tbb::split);
580  inline void operator()(const tbb::blocked_range<size_t>& range);
581  void join(const UpdatePoints& rhs)
582  {
583  if (rhs.mRadius > mRadius) {
584  mRadius = rhs.mRadius;
585  mIndex = rhs.mIndex;
586  }
587  }
588 
589 private:
590  const Vec4s& mSphere;
591  const std::vector<Vec3R>& mPoints;
592  std::vector<float>& mDistances;
593  std::vector<unsigned char>& mMask;
594  bool mOverlapping;
595  float mRadius;
596  int mIndex;
597 };
598 
599 inline
601  const Vec4s& sphere,
602  const std::vector<Vec3R>& points,
603  std::vector<float>& distances,
604  std::vector<unsigned char>& mask,
605  bool overlapping)
606  : mSphere(sphere)
607  , mPoints(points)
608  , mDistances(distances)
609  , mMask(mask)
610  , mOverlapping(overlapping)
611  , mRadius(0.0)
612  , mIndex(0)
613 {
614 }
615 
616 inline
618  : mSphere(rhs.mSphere)
619  , mPoints(rhs.mPoints)
620  , mDistances(rhs.mDistances)
621  , mMask(rhs.mMask)
622  , mOverlapping(rhs.mOverlapping)
623  , mRadius(rhs.mRadius)
624  , mIndex(rhs.mIndex)
625 {
626 }
627 
628 inline void
629 UpdatePoints::run(bool threaded)
630 {
631  if (threaded) {
632  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
633  } else {
634  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
635  }
636 }
637 
638 inline void
639 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
640 {
641  Vec3s pos;
642  for (size_t n = range.begin(); n != range.end(); ++n) {
643  if (mMask[n]) continue;
644 
645  pos.x() = float(mPoints[n].x()) - mSphere[0];
646  pos.y() = float(mPoints[n].y()) - mSphere[1];
647  pos.z() = float(mPoints[n].z()) - mSphere[2];
648 
649  float dist = pos.length();
650 
651  if (dist < mSphere[3]) {
652  mMask[n] = 1;
653  continue;
654  }
655 
656  if (!mOverlapping) {
657  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
658  }
659 
660  if (mDistances[n] > mRadius) {
661  mRadius = mDistances[n];
662  mIndex = int(n);
663  }
664  }
665 }
666 
667 
668 } // namespace v2s_internal
669 
670 
671 ////////////////////////////////////////
672 
673 
674 template<typename GridT, typename InterrupterT>
675 inline void
677  const GridT& grid,
678  std::vector<openvdb::Vec4s>& spheres,
679  int maxSphereCount,
680  bool overlapping,
681  float minRadius,
682  float maxRadius,
683  float isovalue,
684  int instanceCount,
685  InterrupterT* interrupter)
686 {
687  fillWithSpheres(grid, spheres, Vec2i(1, maxSphereCount), overlapping,
688  minRadius, maxRadius, isovalue, instanceCount, interrupter);
689 }
690 
691 
692 template<typename GridT, typename InterrupterT>
693 inline void
695  const GridT& grid,
696  std::vector<openvdb::Vec4s>& spheres,
697  const Vec2i& sphereCount,
698  bool overlapping,
699  float minRadius,
700  float maxRadius,
701  float isovalue,
702  int instanceCount,
703  InterrupterT* interrupter)
704 {
705  spheres.clear();
706 
707  if (grid.empty()) return;
708 
709  const int
710  minSphereCount = sphereCount[0],
711  maxSphereCount = sphereCount[1];
712  if ((minSphereCount > maxSphereCount) || (maxSphereCount < 1)) {
713  OPENVDB_LOG_WARN("fillWithSpheres: minimum sphere count ("
714  << minSphereCount << ") exceeds maximum count (" << maxSphereCount << ")");
715  return;
716  }
717  spheres.reserve(maxSphereCount);
718 
719  auto gridPtr = grid.copy(); // shallow copy
720 
721  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
722  // Clamp the isovalue to the level set's background value minus epsilon.
723  // (In a valid narrow-band level set, all voxels, including background voxels,
724  // have values less than or equal to the background value, so an isovalue
725  // greater than or equal to the background value would produce a mask with
726  // effectively infinite extent.)
727  isovalue = std::min(isovalue,
728  static_cast<float>(gridPtr->background() - math::Tolerance<float>::value()));
729  } else if (gridPtr->getGridClass() == GRID_FOG_VOLUME) {
730  // Clamp the isovalue of a fog volume between epsilon and one,
731  // again to avoid a mask with infinite extent. (Recall that
732  // fog volume voxel values vary from zero outside to one inside.)
733  isovalue = math::Clamp(isovalue, math::Tolerance<float>::value(), 1.f);
734  }
735 
736  // ClosestSurfacePoint is inaccurate for small grids.
737  // Resample the input grid if it is too small.
738  auto numVoxels = gridPtr->activeVoxelCount();
739  if (numVoxels < 10000) {
740  const auto scale = 1.0 / math::Cbrt(2.0 * 10000.0 / double(numVoxels));
741  auto scaledXform = gridPtr->transform().copy();
742  scaledXform->preScale(scale);
743 
744  auto newGridPtr = levelSetRebuild(*gridPtr, isovalue,
745  LEVEL_SET_HALF_WIDTH, LEVEL_SET_HALF_WIDTH, scaledXform.get(), interrupter);
746 
747  const auto newNumVoxels = newGridPtr->activeVoxelCount();
748  if (newNumVoxels > numVoxels) {
749  OPENVDB_LOG_DEBUG_RUNTIME("fillWithSpheres: resampled input grid from "
750  << numVoxels << " voxel" << (numVoxels == 1 ? "" : "s")
751  << " to " << newNumVoxels << " voxel" << (newNumVoxels == 1 ? "" : "s"));
752  gridPtr = newGridPtr;
753  numVoxels = newNumVoxels;
754  }
755  }
756 
757  const bool addNarrowBandPoints = (numVoxels < 10000);
758  int instances = std::max(instanceCount, maxSphereCount);
759 
760  using TreeT = typename GridT::TreeType;
761  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
762  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
763 
764  using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
765  0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>; // mt11213b
766  RandGen mtRand(/*seed=*/0);
767 
768  const TreeT& tree = gridPtr->tree();
769  math::Transform transform = gridPtr->transform();
770 
771  std::vector<Vec3R> instancePoints;
772  {
773  // Compute a mask of the voxels enclosed by the isosurface.
774  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
775  if (gridPtr->getGridClass() == GRID_LEVEL_SET) {
776  interiorMaskPtr = sdfInteriorMask(*gridPtr, isovalue);
777  } else {
778  // For non-level-set grids, the interior mask comprises the active voxels.
779  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
780  interiorMaskPtr->setTransform(transform.copy());
781  interiorMaskPtr->tree().topologyUnion(tree);
782  }
783 
784  if (interrupter && interrupter->wasInterrupted()) return;
785 
786  // If the interior mask is small and eroding it results in an empty grid,
787  // use the uneroded mask instead. (But if the minimum sphere count is zero,
788  // then eroding away the mask is acceptable.)
789  if (!addNarrowBandPoints || (minSphereCount <= 0)) {
790  erodeVoxels(interiorMaskPtr->tree(), 1);
791  } else {
792  auto& maskTree = interiorMaskPtr->tree();
793  auto copyOfTree = StaticPtrCast<BoolTreeT>(maskTree.copy());
794  erodeVoxels(maskTree, 1);
795  if (maskTree.empty()) { interiorMaskPtr->setTree(copyOfTree); }
796  }
797 
798  // Scatter candidate sphere centroids (instancePoints)
799  instancePoints.reserve(instances);
800  v2s_internal::PointAccessor ptnAcc(instancePoints);
801 
802  const auto scatterCount = Index64(addNarrowBandPoints ? (instances / 2) : instances);
803 
805  ptnAcc, scatterCount, mtRand, 1.0, interrupter);
806  scatter(*interiorMaskPtr);
807  }
808 
809  if (interrupter && interrupter->wasInterrupted()) return;
810 
811  auto csp = ClosestSurfacePoint<GridT>::create(*gridPtr, isovalue, interrupter);
812  if (!csp) return;
813 
814  // Add extra instance points in the interior narrow band.
815  if (instancePoints.size() < size_t(instances)) {
816  const Int16TreeT& signTree = csp->signTree();
817  for (auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
818  for (auto it = leafIt->cbeginValueOn(); it; ++it) {
819  const int flags = int(it.getValue());
820  if (!(volume_to_mesh_internal::EDGES & flags)
822  {
823  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
824  }
825  if (instancePoints.size() == size_t(instances)) break;
826  }
827  if (instancePoints.size() == size_t(instances)) break;
828  }
829  }
830 
831  if (interrupter && interrupter->wasInterrupted()) return;
832 
833  // Assign a radius to each candidate sphere. The radius is the world-space
834  // distance from the sphere's center to the closest surface point.
835  std::vector<float> instanceRadius;
836  if (!csp->search(instancePoints, instanceRadius)) return;
837 
838  float largestRadius = 0.0;
839  int largestRadiusIdx = 0;
840  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
841  if (instanceRadius[n] > largestRadius) {
842  largestRadius = instanceRadius[n];
843  largestRadiusIdx = int(n);
844  }
845  }
846 
847  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
848 
849  minRadius = float(minRadius * transform.voxelSize()[0]);
850  maxRadius = float(maxRadius * transform.voxelSize()[0]);
851 
852  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
853 
854  if (interrupter && interrupter->wasInterrupted()) return;
855 
856  largestRadius = std::min(maxRadius, largestRadius);
857 
858  if ((int(s) >= minSphereCount) && (largestRadius < minRadius)) break;
859 
860  const Vec4s sphere(
861  float(instancePoints[largestRadiusIdx].x()),
862  float(instancePoints[largestRadiusIdx].y()),
863  float(instancePoints[largestRadiusIdx].z()),
864  largestRadius);
865 
866  spheres.push_back(sphere);
867  instanceMask[largestRadiusIdx] = 1;
868 
870  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
871  op.run();
872 
873  largestRadius = op.radius();
874  largestRadiusIdx = op.index();
875  }
876 } // fillWithSpheres
877 
878 
879 ////////////////////////////////////////
880 
881 
882 template<typename GridT>
883 template<typename InterrupterT>
884 inline typename ClosestSurfacePoint<GridT>::Ptr
885 ClosestSurfacePoint<GridT>::create(const GridT& grid, float isovalue, InterrupterT* interrupter)
886 {
887  auto csp = Ptr{new ClosestSurfacePoint};
888  if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
889  return csp;
890 }
891 
892 
893 template<typename GridT>
894 template<typename InterrupterT>
895 inline bool
897  const GridT& grid, float isovalue, InterrupterT* interrupter)
898 {
899  using Index32LeafManagerT = tree::LeafManager<Index32TreeT>;
900  using ValueT = typename GridT::ValueType;
901 
902  const TreeT& tree = grid.tree();
903  const math::Transform& transform = grid.transform();
904 
905  { // Extract surface point cloud
906 
907  BoolTreeT mask(false);
909 
910  mSignTreePt.reset(new Int16TreeT(0));
911  mIdxTreePt.reset(new Index32TreeT(std::numeric_limits<Index32>::max()));
912 
913 
915  *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
916 
917  if (interrupter && interrupter->wasInterrupted()) return false;
918 
919  // count unique points
920 
921  using Int16LeafNodeType = typename Int16TreeT::LeafNodeType;
922  using Index32LeafNodeType = typename Index32TreeT::LeafNodeType;
923 
924  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
925  mSignTreePt->getNodes(signFlagsLeafNodes);
926 
927  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
928 
929  hboost::scoped_array<Index32> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
930 
931  tbb::parallel_for(auxiliaryLeafNodeRange,
933  (signFlagsLeafNodes, leafNodeOffsets));
934 
935  {
936  Index32 pointCount = 0;
937  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
938  const Index32 tmp = leafNodeOffsets[n];
939  leafNodeOffsets[n] = pointCount;
940  pointCount += tmp;
941  }
942 
943  mPointListSize = size_t(pointCount);
944  mSurfacePointList.reset(new Vec3s[mPointListSize]);
945  }
946 
947 
948  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
949  mIdxTreePt->getNodes(pointIndexLeafNodes);
950 
951  tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
952  mSurfacePointList.get(), tree, pointIndexLeafNodes,
953  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
954  }
955 
956  if (interrupter && interrupter->wasInterrupted()) return false;
957 
958  Index32LeafManagerT idxLeafs(*mIdxTreePt);
959 
960  using Index32RootNodeT = typename Index32TreeT::RootNodeType;
961  using Index32NodeChainT = typename Index32RootNodeT::NodeChainType;
963  "expected tree depth greater than one");
964  using Index32InternalNodeT =
965  typename hboost::mpl::at<Index32NodeChainT, hboost::mpl::int_<1> >::type;
966 
967  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
968  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
969  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
970 
971  std::vector<const Index32InternalNodeT*> internalNodes;
972 
973  const Index32InternalNodeT* node = nullptr;
974  for (; nIt; ++nIt) {
975  nIt.getNode(node);
976  if (node) internalNodes.push_back(node);
977  }
978 
979  std::vector<IndexRange>().swap(mLeafRanges);
980  mLeafRanges.resize(internalNodes.size());
981 
982  std::vector<const Index32LeafT*>().swap(mLeafNodes);
983  mLeafNodes.reserve(idxLeafs.leafCount());
984 
985  typename Index32InternalNodeT::ChildOnCIter leafIt;
986  mMaxNodeLeafs = 0;
987  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
988 
989  mLeafRanges[n].first = mLeafNodes.size();
990 
991  size_t leafCount = 0;
992  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
993  mLeafNodes.push_back(&(*leafIt));
994  ++leafCount;
995  }
996 
997  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
998 
999  mLeafRanges[n].second = mLeafNodes.size();
1000  }
1001 
1002  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
1003  mLeafBoundingSpheres.resize(mLeafNodes.size());
1004 
1005  v2s_internal::LeafOp<Index32LeafT> leafBS(
1006  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
1007  leafBS.run();
1008 
1009 
1010  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
1011  mNodeBoundingSpheres.resize(internalNodes.size());
1012 
1013  v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
1014  nodeBS.run();
1015  return true;
1016 } // ClosestSurfacePoint::initialize
1017 
1018 
1019 template<typename GridT>
1020 inline bool
1021 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
1022  std::vector<float>& distances, bool transformPoints)
1023 {
1024  distances.clear();
1025  distances.resize(points.size(), std::numeric_limits<float>::infinity());
1026 
1027  v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
1028  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1029  mMaxNodeLeafs, transformPoints);
1030 
1031  cpd.run();
1032 
1033  return true;
1034 }
1035 
1036 
1037 template<typename GridT>
1038 inline bool
1039 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
1040 {
1041  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
1042 }
1043 
1044 
1045 template<typename GridT>
1046 inline bool
1048  std::vector<float>& distances)
1049 {
1050  return search(points, distances, true);
1051 }
1052 
1053 } // namespace tools
1054 } // namespace OPENVDB_VERSION_NAME
1055 } // namespace openvdb
1056 
1057 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
1058 
1059 // Copyright (c) 2012-2018 DreamWorks Animation LLC
1060 // All rights reserved. This software is distributed under the
1061 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Overwrite each input point with its closest surface point.
GLint first
Definition: glcorearb.h:404
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
GLenum GLint * range
Definition: glcorearb.h:1924
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)
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:343
void operator()(const tbb::blocked_range< size_t > &) const
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:798
GLint GLuint mask
Definition: glcorearb.h:123
GLbitfield flags
Definition: glcorearb.h:1595
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
GLint y
Definition: glcorearb.h:102
png_uint_32 i
Definition: png.h:2877
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:723
GLsizeiptr size
Definition: glcorearb.h:663
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.
GA_API const UT_StringHolder scale
LeafOp(std::vector< Vec4R > &leafBoundingSpheres, const std::vector< const Index32LeafT * > &leafNodes, const math::Transform &transform, const PointList &surfacePointList)
GLdouble n
Definition: glcorearb.h:2007
GLfloat f
Definition: glcorearb.h:1925
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
exint length() const
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:49
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
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)
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
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
GA_API const UT_StringHolder transform
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
GLint GLsizei count
Definition: glcorearb.h:404
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:878
GLsizei const GLfloat * value
Definition: glcorearb.h:823
#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)
typedef int
Definition: png.h:1175
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 index
Definition: glcorearb.h:785
GLint GLenum GLint x
Definition: glcorearb.h:408
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...
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
typename TreeT::template ValueConverter< Int16 >::Type Int16TreeT
void swap(TfErrorTransport &l, TfErrorTransport &r)
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...
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1113
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
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)
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:794
void setTree(TreeBase::Ptr) override
Associate the given tree with this grid, in place of its existing tree.
Definition: Grid.h:1310