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