OpenVDB  5.0.0
VolumeToSpheres.h
Go to the documentation of this file.
1 //
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 //
30 
34 
35 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
37 
38 #include <openvdb/tree/LeafManager.h>
39 #include "Morphology.h" // for erodeVoxels()
40 #include "PointScatter.h"
41 #include "LevelSetUtil.h"
42 #include "VolumeToMesh.h"
43 
44 #include <boost/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 <memory>
51 #include <vector>
52 
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
78 template<typename GridT, typename InterrupterT = util::NullInterrupter>
79 inline void
81  const GridT& grid,
82  std::vector<openvdb::Vec4s>& spheres,
83  int maxSphereCount,
84  bool overlapping = false,
85  float minRadius = 1.0,
86  float maxRadius = std::numeric_limits<float>::max(),
87  float isovalue = 0.0,
88  int instanceCount = 10000,
89  InterrupterT* interrupter = nullptr);
90 
91 
93 
94 
98 template<typename GridT>
100 {
101 public:
102  using Ptr = std::unique_ptr<ClosestSurfacePoint>;
103  using TreeT = typename GridT::TreeType;
104  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
105  using Index32TreeT = typename TreeT::template ValueConverter<Index32>::Type;
106  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
107 
119  template<typename InterrupterT = util::NullInterrupter>
120  static inline Ptr create(const GridT& grid, float isovalue = 0.0,
121  InterrupterT* interrupter = nullptr);
122 
126  inline bool search(const std::vector<Vec3R>& points, std::vector<float>& distances);
127 
131  inline bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
132 
134  const Index32TreeT& indexTree() const { return *mIdxTreePt; }
136  const Int16TreeT& signTree() const { return *mSignTreePt; }
137 
138 private:
139  using Index32LeafT = typename Index32TreeT::LeafNodeType;
140  using IndexRange = std::pair<size_t, size_t>;
141 
142  std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
143  std::vector<IndexRange> mLeafRanges;
144  std::vector<const Index32LeafT*> mLeafNodes;
145  PointList mSurfacePointList;
146  size_t mPointListSize = 0, mMaxNodeLeafs = 0;
147  typename Index32TreeT::Ptr mIdxTreePt;
148  typename Int16TreeT::Ptr mSignTreePt;
149 
150  ClosestSurfacePoint() = default;
151  template<typename InterrupterT = util::NullInterrupter>
152  inline bool initialize(const GridT&, float isovalue, InterrupterT*);
153  inline bool search(std::vector<Vec3R>&, std::vector<float>&, bool transformPoints);
154 };
155 
156 
158 
159 
160 // Internal utility methods
161 
162 namespace v2s_internal {
163 
165 {
166  PointAccessor(std::vector<Vec3R>& points)
167  : mPoints(points)
168  {
169  }
170 
171  void add(const Vec3R &pos)
172  {
173  mPoints.push_back(pos);
174  }
175 private:
176  std::vector<Vec3R>& mPoints;
177 };
178 
179 
180 template<typename Index32LeafT>
181 class LeafOp
182 {
183 public:
184 
185  LeafOp(std::vector<Vec4R>& leafBoundingSpheres,
186  const std::vector<const Index32LeafT*>& leafNodes,
187  const math::Transform& transform,
188  const PointList& surfacePointList);
189 
190  void run(bool threaded = true);
191 
192 
193  void operator()(const tbb::blocked_range<size_t>&) const;
194 
195 private:
196  std::vector<Vec4R>& mLeafBoundingSpheres;
197  const std::vector<const Index32LeafT*>& mLeafNodes;
198  const math::Transform& mTransform;
199  const PointList& mSurfacePointList;
200 };
201 
202 template<typename Index32LeafT>
204  std::vector<Vec4R>& leafBoundingSpheres,
205  const std::vector<const Index32LeafT*>& leafNodes,
206  const math::Transform& transform,
207  const PointList& surfacePointList)
208  : mLeafBoundingSpheres(leafBoundingSpheres)
209  , mLeafNodes(leafNodes)
210  , mTransform(transform)
211  , mSurfacePointList(surfacePointList)
212 {
213 }
214 
215 template<typename Index32LeafT>
216 void
218 {
219  if (threaded) {
220  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *this);
221  } else {
222  (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
223  }
224 }
225 
226 template<typename Index32LeafT>
227 void
228 LeafOp<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
229 {
230  typename Index32LeafT::ValueOnCIter iter;
231  Vec3s avg;
232 
233  for (size_t n = range.begin(); n != range.end(); ++n) {
234  avg[0] = 0.0;
235  avg[1] = 0.0;
236  avg[2] = 0.0;
237 
238  int count = 0;
239  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
240  avg += mSurfacePointList[iter.getValue()];
241  ++count;
242  }
243  if (count > 1) avg *= float(1.0 / double(count));
244 
245  float maxDist = 0.0;
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  sphere[0] = avg[0];
253  sphere[1] = avg[1];
254  sphere[2] = avg[2];
255  sphere[3] = std::sqrt(maxDist);
256  }
257 }
258 
259 
260 class NodeOp
261 {
262 public:
263  using IndexRange = std::pair<size_t, size_t>;
264 
265  NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
266  const std::vector<IndexRange>& leafRanges,
267  const std::vector<Vec4R>& leafBoundingSpheres);
268 
269  inline void run(bool threaded = true);
270 
271  inline void operator()(const tbb::blocked_range<size_t>&) const;
272 
273 private:
274  std::vector<Vec4R>& mNodeBoundingSpheres;
275  const std::vector<IndexRange>& mLeafRanges;
276  const std::vector<Vec4R>& mLeafBoundingSpheres;
277 };
278 
279 inline
280 NodeOp::NodeOp(std::vector<Vec4R>& nodeBoundingSpheres,
281  const std::vector<IndexRange>& leafRanges,
282  const std::vector<Vec4R>& leafBoundingSpheres)
283  : mNodeBoundingSpheres(nodeBoundingSpheres)
284  , mLeafRanges(leafRanges)
285  , mLeafBoundingSpheres(leafBoundingSpheres)
286 {
287 }
288 
289 inline void
290 NodeOp::run(bool threaded)
291 {
292  if (threaded) {
293  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *this);
294  } else {
295  (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
296  }
297 }
298 
299 inline void
300 NodeOp::operator()(const tbb::blocked_range<size_t>& range) const
301 {
302  Vec3d avg, pos;
303 
304  for (size_t n = range.begin(); n != range.end(); ++n) {
305 
306  avg[0] = 0.0;
307  avg[1] = 0.0;
308  avg[2] = 0.0;
309 
310  int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
311 
312  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
313  avg[0] += mLeafBoundingSpheres[i][0];
314  avg[1] += mLeafBoundingSpheres[i][1];
315  avg[2] += mLeafBoundingSpheres[i][2];
316  }
317 
318  if (count > 1) avg *= float(1.0 / double(count));
319 
320 
321  double maxDist = 0.0;
322 
323  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
324  pos[0] = mLeafBoundingSpheres[i][0];
325  pos[1] = mLeafBoundingSpheres[i][1];
326  pos[2] = mLeafBoundingSpheres[i][2];
327 
328  double tmpDist = (pos - avg).length() + mLeafBoundingSpheres[i][3];
329  if (tmpDist > maxDist) maxDist = tmpDist;
330  }
331 
332  Vec4R& sphere = mNodeBoundingSpheres[n];
333 
334  sphere[0] = avg[0];
335  sphere[1] = avg[1];
336  sphere[2] = avg[2];
337  sphere[3] = maxDist;
338  }
339 }
340 
341 
343 
344 
345 template<typename Index32LeafT>
347 {
348 public:
349  using IndexRange = std::pair<size_t, size_t>;
350 
352  std::vector<Vec3R>& instancePoints,
353  std::vector<float>& instanceDistances,
354  const PointList& surfacePointList,
355  const std::vector<const Index32LeafT*>& leafNodes,
356  const std::vector<IndexRange>& leafRanges,
357  const std::vector<Vec4R>& leafBoundingSpheres,
358  const std::vector<Vec4R>& nodeBoundingSpheres,
359  size_t maxNodeLeafs,
360  bool transformPoints = false);
361 
362 
363  void run(bool threaded = true);
364 
365 
366  void operator()(const tbb::blocked_range<size_t>&) const;
367 
368 private:
369 
370  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
371  void evalNode(size_t pointIndex, size_t nodeIndex) const;
372 
373 
374  std::vector<Vec3R>& mInstancePoints;
375  std::vector<float>& mInstanceDistances;
376 
377  const PointList& mSurfacePointList;
378 
379  const std::vector<const Index32LeafT*>& mLeafNodes;
380  const std::vector<IndexRange>& mLeafRanges;
381  const std::vector<Vec4R>& mLeafBoundingSpheres;
382  const std::vector<Vec4R>& mNodeBoundingSpheres;
383 
384  std::vector<float> mLeafDistances, mNodeDistances;
385 
386  const bool mTransformPoints;
387  size_t mClosestPointIndex;
388 };// ClosestPointDist
389 
390 
391 template<typename Index32LeafT>
393  std::vector<Vec3R>& instancePoints,
394  std::vector<float>& instanceDistances,
395  const PointList& surfacePointList,
396  const std::vector<const Index32LeafT*>& leafNodes,
397  const std::vector<IndexRange>& leafRanges,
398  const std::vector<Vec4R>& leafBoundingSpheres,
399  const std::vector<Vec4R>& nodeBoundingSpheres,
400  size_t maxNodeLeafs,
401  bool transformPoints)
402  : mInstancePoints(instancePoints)
403  , mInstanceDistances(instanceDistances)
404  , mSurfacePointList(surfacePointList)
405  , mLeafNodes(leafNodes)
406  , mLeafRanges(leafRanges)
407  , mLeafBoundingSpheres(leafBoundingSpheres)
408  , mNodeBoundingSpheres(nodeBoundingSpheres)
409  , mLeafDistances(maxNodeLeafs, 0.0)
410  , mNodeDistances(leafRanges.size(), 0.0)
411  , mTransformPoints(transformPoints)
412  , mClosestPointIndex(0)
413 {
414 }
415 
416 
417 template<typename Index32LeafT>
418 void
420 {
421  if (threaded) {
422  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
423  } else {
424  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
425  }
426 }
427 
428 template<typename Index32LeafT>
429 void
430 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
431 {
432  typename Index32LeafT::ValueOnCIter iter;
433  const Vec3s center = mInstancePoints[index];
434  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
435 
436  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
437 
438  const Vec3s& point = mSurfacePointList[iter.getValue()];
439  float tmpDist = (point - center).lengthSqr();
440 
441  if (tmpDist < mInstanceDistances[index]) {
442  mInstanceDistances[index] = tmpDist;
443  closestPointIndex = iter.getValue();
444  }
445  }
446 }
447 
448 
449 template<typename Index32LeafT>
450 void
451 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
452 {
453  if (nodeIndex >= mLeafRanges.size()) return;
454 
455  const Vec3R& pos = mInstancePoints[pointIndex];
456  float minDist = mInstanceDistances[pointIndex];
457  size_t minDistIdx = 0;
458  Vec3R center;
459  bool updatedDist = false;
460 
461  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
462  i < mLeafRanges[nodeIndex].second; ++i, ++n)
463  {
464  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
465 
466  center[0] = mLeafBoundingSpheres[i][0];
467  center[1] = mLeafBoundingSpheres[i][1];
468  center[2] = mLeafBoundingSpheres[i][2];
469  const auto radius = mLeafBoundingSpheres[i][3];
470 
471  distToLeaf = float(std::max(0.0, (pos - center).length() - radius));
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  const auto radius = mNodeBoundingSpheres[i][3];
512 
513  distToNode = float(std::max(0.0, (pos - center).length() - radius));
514 
515  if (distToNode < minDist) {
516  minDist = distToNode;
517  minDistIdx = i;
518  }
519  }
520 
521  evalNode(n, minDistIdx);
522 
523  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
524  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
525  evalNode(n, i);
526  }
527  }
528 
529  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
530 
531  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
532  }
533 }
534 
535 
537 {
538 public:
539  UpdatePoints(
540  const Vec4s& sphere,
541  const std::vector<Vec3R>& points,
542  std::vector<float>& distances,
543  std::vector<unsigned char>& mask,
544  bool overlapping);
545 
546  float radius() const { return mRadius; }
547  int index() const { return mIndex; }
548 
549  inline void run(bool threaded = true);
550 
551 
552  UpdatePoints(UpdatePoints&, tbb::split);
553  inline void operator()(const tbb::blocked_range<size_t>& range);
554  void join(const UpdatePoints& rhs)
555  {
556  if (rhs.mRadius > mRadius) {
557  mRadius = rhs.mRadius;
558  mIndex = rhs.mIndex;
559  }
560  }
561 
562 private:
563  const Vec4s& mSphere;
564  const std::vector<Vec3R>& mPoints;
565  std::vector<float>& mDistances;
566  std::vector<unsigned char>& mMask;
567  bool mOverlapping;
568  float mRadius;
569  int mIndex;
570 };
571 
572 inline
574  const Vec4s& sphere,
575  const std::vector<Vec3R>& points,
576  std::vector<float>& distances,
577  std::vector<unsigned char>& mask,
578  bool overlapping)
579  : mSphere(sphere)
580  , mPoints(points)
581  , mDistances(distances)
582  , mMask(mask)
583  , mOverlapping(overlapping)
584  , mRadius(0.0)
585  , mIndex(0)
586 {
587 }
588 
589 inline
591  : mSphere(rhs.mSphere)
592  , mPoints(rhs.mPoints)
593  , mDistances(rhs.mDistances)
594  , mMask(rhs.mMask)
595  , mOverlapping(rhs.mOverlapping)
596  , mRadius(rhs.mRadius)
597  , mIndex(rhs.mIndex)
598 {
599 }
600 
601 inline void
602 UpdatePoints::run(bool threaded)
603 {
604  if (threaded) {
605  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
606  } else {
607  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
608  }
609 }
610 
611 inline void
612 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
613 {
614  Vec3s pos;
615  for (size_t n = range.begin(); n != range.end(); ++n) {
616  if (mMask[n]) continue;
617 
618  pos.x() = float(mPoints[n].x()) - mSphere[0];
619  pos.y() = float(mPoints[n].y()) - mSphere[1];
620  pos.z() = float(mPoints[n].z()) - mSphere[2];
621 
622  float dist = pos.length();
623 
624  if (dist < mSphere[3]) {
625  mMask[n] = 1;
626  continue;
627  }
628 
629  if (!mOverlapping) {
630  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
631  }
632 
633  if (mDistances[n] > mRadius) {
634  mRadius = mDistances[n];
635  mIndex = int(n);
636  }
637  }
638 }
639 
640 
641 } // namespace v2s_internal
642 
643 
645 
646 
647 template<typename GridT, typename InterrupterT>
648 inline void
650  const GridT& grid,
651  std::vector<openvdb::Vec4s>& spheres,
652  int maxSphereCount,
653  bool overlapping,
654  float minRadius,
655  float maxRadius,
656  float isovalue,
657  int instanceCount,
658  InterrupterT* interrupter)
659 {
660  spheres.clear();
661  spheres.reserve(maxSphereCount);
662 
663  const bool addNBPoints = grid.activeVoxelCount() < 10000;
664  int instances = std::max(instanceCount, maxSphereCount);
665 
666  using TreeT = typename GridT::TreeType;
667  using ValueT = typename GridT::ValueType;
668 
669  using BoolTreeT = typename TreeT::template ValueConverter<bool>::Type;
670  using Int16TreeT = typename TreeT::template ValueConverter<Int16>::Type;
671 
672  using RandGen = std::mersenne_twister_engine<uint32_t, 32, 351, 175, 19,
673  0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253>; // mt11213b
674  RandGen mtRand(/*seed=*/0);
675 
676  const TreeT& tree = grid.tree();
677  const math::Transform& transform = grid.transform();
678 
679  std::vector<Vec3R> instancePoints;
680  {
681  // Compute a mask of the voxels enclosed by the isosurface.
682  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
683  if (grid.getGridClass() == GRID_LEVEL_SET) {
684  // Clamp the isovalue to the level set's background value minus epsilon.
685  // (In a valid narrow-band level set, all voxels, including background voxels,
686  // have values less than or equal to the background value, so an isovalue
687  // greater than or equal to the background value would produce a mask with
688  // effectively infinite extent.)
689  isovalue = std::min(isovalue,
690  static_cast<float>(tree.background() - math::Tolerance<ValueT>::value()));
691  interiorMaskPtr = sdfInteriorMask(grid, ValueT(isovalue));
692  } else {
693  if (grid.getGridClass() == GRID_FOG_VOLUME) {
694  // Clamp the isovalue of a fog volume between epsilon and one,
695  // again to avoid a mask with infinite extent. (Recall that
696  // fog volume voxel values vary from zero outside to one inside.)
697  isovalue = math::Clamp(isovalue, math::Tolerance<float>::value(), 1.f);
698  }
699  // For non-level-set grids, the interior mask comprises the active voxels.
700  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
701  interiorMaskPtr->setTransform(transform.copy());
702  interiorMaskPtr->tree().topologyUnion(tree);
703  }
704 
705  if (interrupter && interrupter->wasInterrupted()) return;
706 
707  erodeVoxels(interiorMaskPtr->tree(), 1);
708 
709  // Scatter candidate sphere centroids (instancePoints)
710  instancePoints.reserve(instances);
711  v2s_internal::PointAccessor ptnAcc(instancePoints);
712 
714  ptnAcc, Index64(addNBPoints ? (instances / 2) : instances), mtRand, 1.0, interrupter);
715 
716  scatter(*interiorMaskPtr);
717  }
718 
719  if (interrupter && interrupter->wasInterrupted()) return;
720 
721  auto csp = ClosestSurfacePoint<GridT>::create(grid, isovalue, interrupter);
722  if (!csp) return;
723 
724  // Add extra instance points in the interior narrow band.
725  if (instancePoints.size() < size_t(instances)) {
726  const Int16TreeT& signTree = csp->signTree();
727  for (auto leafIt = signTree.cbeginLeaf(); leafIt; ++leafIt) {
728  for (auto it = leafIt->cbeginValueOn(); it; ++it) {
729  const int flags = int(it.getValue());
730  if (!(volume_to_mesh_internal::EDGES & flags)
731  && (volume_to_mesh_internal::INSIDE & flags))
732  {
733  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
734  }
735  if (instancePoints.size() == size_t(instances)) break;
736  }
737  if (instancePoints.size() == size_t(instances)) break;
738  }
739  }
740 
741  if (interrupter && interrupter->wasInterrupted()) return;
742 
743  std::vector<float> instanceRadius;
744  if (!csp->search(instancePoints, instanceRadius)) return;
745 
746  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
747  float largestRadius = 0.0;
748  int largestRadiusIdx = 0;
749 
750  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
751  if (instanceRadius[n] > largestRadius) {
752  largestRadius = instanceRadius[n];
753  largestRadiusIdx = int(n);
754  }
755  }
756 
757  Vec3s pos;
758  Vec4s sphere;
759  minRadius = float(minRadius * transform.voxelSize()[0]);
760  maxRadius = float(maxRadius * transform.voxelSize()[0]);
761 
762  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
763 
764  if (interrupter && interrupter->wasInterrupted()) return;
765 
766  largestRadius = std::min(maxRadius, largestRadius);
767 
768  if (s != 0 && largestRadius < minRadius) break;
769 
770  sphere[0] = float(instancePoints[largestRadiusIdx].x());
771  sphere[1] = float(instancePoints[largestRadiusIdx].y());
772  sphere[2] = float(instancePoints[largestRadiusIdx].z());
773  sphere[3] = largestRadius;
774 
775  spheres.push_back(sphere);
776  instanceMask[largestRadiusIdx] = 1;
777 
779  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
780  op.run();
781 
782  largestRadius = op.radius();
783  largestRadiusIdx = op.index();
784  }
785 } // fillWithSpheres
786 
787 
789 
790 
791 template<typename GridT>
792 template<typename InterrupterT>
793 inline typename ClosestSurfacePoint<GridT>::Ptr
794 ClosestSurfacePoint<GridT>::create(const GridT& grid, float isovalue, InterrupterT* interrupter)
795 {
796  auto csp = Ptr{new ClosestSurfacePoint};
797  if (!csp->initialize(grid, isovalue, interrupter)) csp.reset();
798  return csp;
799 }
800 
801 
802 template<typename GridT>
803 template<typename InterrupterT>
804 inline bool
806  const GridT& grid, float isovalue, InterrupterT* interrupter)
807 {
808  using Index32LeafManagerT = tree::LeafManager<Index32TreeT>;
809  using ValueT = typename GridT::ValueType;
810 
811  const TreeT& tree = grid.tree();
812  const math::Transform& transform = grid.transform();
813 
814  { // Extract surface point cloud
815 
816  BoolTreeT mask(false);
818 
819  mSignTreePt.reset(new Int16TreeT(0));
820  mIdxTreePt.reset(new Index32TreeT(boost::integer_traits<Index32>::const_max));
821 
822 
824  *mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
825 
826  if (interrupter && interrupter->wasInterrupted()) return false;
827 
828  // count unique points
829 
830  using Int16LeafNodeType = typename Int16TreeT::LeafNodeType;
831  using Index32LeafNodeType = typename Index32TreeT::LeafNodeType;
832 
833  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
834  mSignTreePt->getNodes(signFlagsLeafNodes);
835 
836  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
837 
838  boost::scoped_array<Index32> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
839 
840  tbb::parallel_for(auxiliaryLeafNodeRange,
842  (signFlagsLeafNodes, leafNodeOffsets));
843 
844  {
845  Index32 pointCount = 0;
846  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
847  const Index32 tmp = leafNodeOffsets[n];
848  leafNodeOffsets[n] = pointCount;
849  pointCount += tmp;
850  }
851 
852  mPointListSize = size_t(pointCount);
853  mSurfacePointList.reset(new Vec3s[mPointListSize]);
854  }
855 
856 
857  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
858  mIdxTreePt->getNodes(pointIndexLeafNodes);
859 
860  tbb::parallel_for(auxiliaryLeafNodeRange, volume_to_mesh_internal::ComputePoints<TreeT>(
861  mSurfacePointList.get(), tree, pointIndexLeafNodes,
862  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
863  }
864 
865  if (interrupter && interrupter->wasInterrupted()) return false;
866 
867  Index32LeafManagerT idxLeafs(*mIdxTreePt);
868 
869  using Index32RootNodeT = typename Index32TreeT::RootNodeType;
870  using Index32NodeChainT = typename Index32RootNodeT::NodeChainType;
871  BOOST_STATIC_ASSERT(boost::mpl::size<Index32NodeChainT>::value > 1);
872  using Index32InternalNodeT =
873  typename boost::mpl::at<Index32NodeChainT, boost::mpl::int_<1> >::type;
874 
875  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
876  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
877  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
878 
879  std::vector<const Index32InternalNodeT*> internalNodes;
880 
881  const Index32InternalNodeT* node = nullptr;
882  for (; nIt; ++nIt) {
883  nIt.getNode(node);
884  if (node) internalNodes.push_back(node);
885  }
886 
887  std::vector<IndexRange>().swap(mLeafRanges);
888  mLeafRanges.resize(internalNodes.size());
889 
890  std::vector<const Index32LeafT*>().swap(mLeafNodes);
891  mLeafNodes.reserve(idxLeafs.leafCount());
892 
893  typename Index32InternalNodeT::ChildOnCIter leafIt;
894  mMaxNodeLeafs = 0;
895  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
896 
897  mLeafRanges[n].first = mLeafNodes.size();
898 
899  size_t leafCount = 0;
900  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
901  mLeafNodes.push_back(&(*leafIt));
902  ++leafCount;
903  }
904 
905  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
906 
907  mLeafRanges[n].second = mLeafNodes.size();
908  }
909 
910  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
911  mLeafBoundingSpheres.resize(mLeafNodes.size());
912 
914  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
915  leafBS.run();
916 
917 
918  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
919  mNodeBoundingSpheres.resize(internalNodes.size());
920 
921  v2s_internal::NodeOp nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
922  nodeBS.run();
923  return true;
924 } // ClosestSurfacePoint::initialize
925 
926 
927 template<typename GridT>
928 inline bool
929 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
930  std::vector<float>& distances, bool transformPoints)
931 {
932  distances.clear();
933  distances.resize(points.size(), std::numeric_limits<float>::infinity());
934 
935  v2s_internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
936  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
937  mMaxNodeLeafs, transformPoints);
938 
939  cpd.run();
940 
941  return true;
942 }
943 
944 
945 template<typename GridT>
946 inline bool
947 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
948 {
949  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
950 }
951 
952 
953 template<typename GridT>
954 inline bool
956  std::vector<float>& distances)
957 {
958  return search(points, distances, true);
959 }
960 
961 } // namespace tools
962 } // namespace OPENVDB_VERSION_NAME
963 } // namespace openvdb
964 
965 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
966 
967 // Copyright (c) 2012-2017 DreamWorks Animation LLC
968 // All rights reserved. This software is distributed under the
969 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
Definition: VolumeToSpheres.h:164
Index64 pointCount(const PointDataTreeT &tree, const bool inCoreOnly=false)
Total points in the PointDataTree.
Definition: PointCount.h:241
static Ptr create(const GridT &grid, float isovalue=0.0, InterrupterT *interrupter=nullptr)
Extract surface points and construct a spatial acceleration structure.
Definition: VolumeToSpheres.h:794
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1112
Definition: Types.h:269
void operator()(const tbb::blocked_range< size_t > &range)
Definition: VolumeToSpheres.h:612
void run(bool threaded=true)
Definition: VolumeToSpheres.h:602
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)
Definition: VolumeToSpheres.h:392
const Int16TreeT & signTree() const
Tree accessor.
Definition: VolumeToSpheres.h:136
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
uint64_t Index64
Definition: Types.h:59
void identifySurfaceIntersectingVoxels(typename InputTreeType::template ValueConverter< bool >::Type &intersectionTree, const InputTreeType &inputTree, typename InputTreeType::ValueType isovalue)
Definition: VolumeToMesh.h:3352
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Compute the distance from each input point to its closest surface point.
Definition: VolumeToSpheres.h:947
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
boost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:180
PointAccessor(std::vector< Vec3R > &points)
Definition: VolumeToSpheres.h:166
We offer three different algorithms (each in its own class) for scattering of points in active voxels...
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:349
int index() const
Definition: VolumeToSpheres.h:547
Vec3< double > Vec3d
Definition: Vec3.h:679
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
void run(bool threaded=true)
Definition: VolumeToSpheres.h:217
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:228
TreeType & tree()
Return a reference to this grid&#39;s tree, which might be shared with other grids.
Definition: Grid.h:797
void add(const Vec3R &pos)
Definition: VolumeToSpheres.h:171
void run(bool threaded=true)
Definition: VolumeToSpheres.h:290
Definition: Transform.h:66
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:230
Tolerance for floating-point comparison.
Definition: Math.h:117
Extract polygonal surfaces from scalar volumes.
uint32_t Index32
Definition: Types.h:58
static T value()
Definition: Math.h:117
Definition: Types.h:270
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Overwrite each input point with its closest surface point.
Definition: VolumeToSpheres.h:955
NodeOp(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
Definition: VolumeToSpheres.h:280
Vec4< float > Vec4s
Definition: Vec4.h:586
Implementation of morphological dilation and erosion.
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
Definition: VolumeToSpheres.h:573
Definition: Exceptions.h:39
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:135
Definition: VolumeToSpheres.h:536
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:496
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: tools/PointScatter.h:114
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
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:877
Definition: Mat4.h:51
void join(const UpdatePoints &rhs)
Definition: VolumeToSpheres.h:554
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2291
float radius() const
Definition: VolumeToSpheres.h:546
const Index32TreeT & indexTree() const
Tree accessor.
Definition: VolumeToSpheres.h:134
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.
Definition: VolumeToSpheres.h:649
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:263
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:343
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
typename GridT::TreeType TreeT
Definition: VolumeToSpheres.h:103
typename TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToSpheres.h:106
typename TreeT::template ValueConverter< Index32 >::Type Index32TreeT
Definition: VolumeToSpheres.h:105
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:55
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
Ptr copy() const
Definition: Transform.h:77
std::unique_ptr< ClosestSurfacePoint > Ptr
Definition: VolumeToSpheres.h:102
OPENVDB_API void initialize()
Global registration of basic types.
Vec3< float > Vec3s
Definition: Vec3.h:678
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)
Definition: VolumeToMesh.h:3894
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:300
void run(bool threaded=true)
Definition: VolumeToSpheres.h:419
Accelerated closest surface point queries for narrow band level sets.
Definition: VolumeToSpheres.h:99
SharedPtr< Grid > Ptr
Definition: Grid.h:502
Definition: VolumeToSpheres.h:181
typename TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToSpheres.h:104
Definition: VolumeToSpheres.h:260