OpenVDB  3.2.0
VolumeToSpheres.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2016 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 
31 #ifndef OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_SPHERES_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/tree/ValueAccessor.h>
35 #include <openvdb/tree/LeafManager.h>
36 #include <openvdb/tools/Morphology.h> // for erodeVoxels()
37 
38 #include <openvdb/tools/PointScatter.h>
39 #include <openvdb/tools/LevelSetUtil.h>
40 #include <openvdb/tools/VolumeToMesh.h>
41 
42 #include <boost/scoped_array.hpp>
43 #include <boost/scoped_ptr.hpp>
44 #include <tbb/blocked_range.h>
45 #include <tbb/parallel_for.h>
46 #include <tbb/parallel_reduce.h>
47 
48 #include <vector>
49 #include <limits> // std::numeric_limits
50 
52 
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
59 
89 template<typename GridT, typename InterrupterT>
90 inline void
92  const GridT& grid,
93  std::vector<openvdb::Vec4s>& spheres,
94  int maxSphereCount,
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 = NULL);
101 
102 
105 template<typename GridT>
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 {
117  fillWithSpheres<GridT, util::NullInterrupter>(grid, spheres,
118  maxSphereCount, overlapping, minRadius, maxRadius, isovalue, instanceCount);
119 }
120 
121 
123 
124 
128 template<typename GridT>
130 {
131 public:
132  typedef typename GridT::TreeType TreeT;
133  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
134  typedef typename TreeT::template ValueConverter<Index32>::Type Index32TreeT;
135  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
136 
137 
139 
140 
152  template<typename InterrupterT>
153  void initialize(const GridT& grid, float isovalue = 0.0, InterrupterT* interrupter = NULL);
154 
155 
158  void initialize(const GridT& grid, float isovalue = 0.0);
159 
160 
161 
168  bool search(const std::vector<Vec3R>& points, std::vector<float>& distances);
169 
170 
178  bool searchAndReplace(std::vector<Vec3R>& points, std::vector<float>& distances);
179 
180 
183  const Index32TreeT& indexTree() const { return *mIdxTreePt; }
184  const Int16TreeT& signTree() const { return *mSignTreePt; }
186 
187 private:
188  typedef typename Index32TreeT::LeafNodeType Index32LeafT;
189  typedef std::pair<size_t, size_t> IndexRange;
190 
191  bool mIsInitialized;
192  std::vector<Vec4R> mLeafBoundingSpheres, mNodeBoundingSpheres;
193  std::vector<IndexRange> mLeafRanges;
194  std::vector<const Index32LeafT*> mLeafNodes;
195  PointList mSurfacePointList;
196  size_t mPointListSize, mMaxNodeLeafs;
197  float mMaxRadiusSqr;
198  typename Index32TreeT::Ptr mIdxTreePt;
199  typename Int16TreeT::Ptr mSignTreePt;
200 
201  bool search(std::vector<Vec3R>&, std::vector<float>&, bool transformPoints);
202 };
203 
204 
206 
207 
208 
209 
210 // Internal utility methods
211 
212 
213 namespace internal {
214 
216 {
217  PointAccessor(std::vector<Vec3R>& points)
218  : mPoints(points)
219  {
220  }
221 
222  void add(const Vec3R &pos)
223  {
224  mPoints.push_back(pos);
225  }
226 private:
227  std::vector<Vec3R>& mPoints;
228 };
229 
230 
231 template<typename Index32LeafT>
232 class LeafBS
233 {
234 public:
235 
236  LeafBS(std::vector<Vec4R>& leafBoundingSpheres,
237  const std::vector<const Index32LeafT*>& leafNodes,
238  const math::Transform& transform,
239  const PointList& surfacePointList);
240 
241  void run(bool threaded = true);
242 
243 
244  void operator()(const tbb::blocked_range<size_t>&) const;
245 
246 private:
247  std::vector<Vec4R>& mLeafBoundingSpheres;
248  const std::vector<const Index32LeafT*>& mLeafNodes;
249  const math::Transform& mTransform;
250  const PointList& mSurfacePointList;
251 };
252 
253 template<typename Index32LeafT>
255  std::vector<Vec4R>& leafBoundingSpheres,
256  const std::vector<const Index32LeafT*>& leafNodes,
257  const math::Transform& transform,
258  const PointList& surfacePointList)
259  : mLeafBoundingSpheres(leafBoundingSpheres)
260  , mLeafNodes(leafNodes)
261  , mTransform(transform)
262  , mSurfacePointList(surfacePointList)
263 {
264 }
265 
266 template<typename Index32LeafT>
267 void
269 {
270  if (threaded) {
271  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafNodes.size()), *this);
272  } else {
273  (*this)(tbb::blocked_range<size_t>(0, mLeafNodes.size()));
274  }
275 }
276 
277 template<typename Index32LeafT>
278 void
279 LeafBS<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
280 {
281  typename Index32LeafT::ValueOnCIter iter;
282  Vec3s avg;
283 
284  for (size_t n = range.begin(); n != range.end(); ++n) {
285 
286  avg[0] = 0.0;
287  avg[1] = 0.0;
288  avg[2] = 0.0;
289 
290  int count = 0;
291  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
292  avg += mSurfacePointList[iter.getValue()];
293  ++count;
294  }
295 
296  if (count > 1) avg *= float(1.0 / double(count));
297 
298  float maxDist = 0.0;
299 
300  for (iter = mLeafNodes[n]->cbeginValueOn(); iter; ++iter) {
301  float tmpDist = (mSurfacePointList[iter.getValue()] - avg).lengthSqr();
302  if (tmpDist > maxDist) maxDist = tmpDist;
303  }
304 
305  Vec4R& sphere = mLeafBoundingSpheres[n];
306 
307  sphere[0] = avg[0];
308  sphere[1] = avg[1];
309  sphere[2] = avg[2];
310  sphere[3] = maxDist * 2.0; // padded radius
311  }
312 }
313 
314 
315 class NodeBS
316 {
317 public:
318  typedef std::pair<size_t, size_t> IndexRange;
319 
320  NodeBS(std::vector<Vec4R>& nodeBoundingSpheres,
321  const std::vector<IndexRange>& leafRanges,
322  const std::vector<Vec4R>& leafBoundingSpheres);
323 
324  inline void run(bool threaded = true);
325 
326  inline void operator()(const tbb::blocked_range<size_t>&) const;
327 
328 private:
329  std::vector<Vec4R>& mNodeBoundingSpheres;
330  const std::vector<IndexRange>& mLeafRanges;
331  const std::vector<Vec4R>& mLeafBoundingSpheres;
332 };
333 
334 inline
335 NodeBS::NodeBS(std::vector<Vec4R>& nodeBoundingSpheres,
336  const std::vector<IndexRange>& leafRanges,
337  const std::vector<Vec4R>& leafBoundingSpheres)
338  : mNodeBoundingSpheres(nodeBoundingSpheres)
339  , mLeafRanges(leafRanges)
340  , mLeafBoundingSpheres(leafBoundingSpheres)
341 {
342 }
343 
344 inline void
345 NodeBS::run(bool threaded)
346 {
347  if (threaded) {
348  tbb::parallel_for(tbb::blocked_range<size_t>(0, mLeafRanges.size()), *this);
349  } else {
350  (*this)(tbb::blocked_range<size_t>(0, mLeafRanges.size()));
351  }
352 }
353 
354 inline void
355 NodeBS::operator()(const tbb::blocked_range<size_t>& range) const
356 {
357  Vec3d avg, pos;
358 
359  for (size_t n = range.begin(); n != range.end(); ++n) {
360 
361  avg[0] = 0.0;
362  avg[1] = 0.0;
363  avg[2] = 0.0;
364 
365  int count = int(mLeafRanges[n].second) - int(mLeafRanges[n].first);
366 
367  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
368  avg[0] += mLeafBoundingSpheres[i][0];
369  avg[1] += mLeafBoundingSpheres[i][1];
370  avg[2] += mLeafBoundingSpheres[i][2];
371  }
372 
373  if (count > 1) avg *= float(1.0 / double(count));
374 
375 
376  double maxDist = 0.0;
377 
378  for (size_t i = mLeafRanges[n].first; i < mLeafRanges[n].second; ++i) {
379  pos[0] = mLeafBoundingSpheres[i][0];
380  pos[1] = mLeafBoundingSpheres[i][1];
381  pos[2] = mLeafBoundingSpheres[i][2];
382 
383  double tmpDist = (pos - avg).lengthSqr() + mLeafBoundingSpheres[i][3];
384  if (tmpDist > maxDist) maxDist = tmpDist;
385  }
386 
387  Vec4R& sphere = mNodeBoundingSpheres[n];
388 
389  sphere[0] = avg[0];
390  sphere[1] = avg[1];
391  sphere[2] = avg[2];
392  sphere[3] = maxDist * 2.0; // padded radius
393  }
394 }
395 
396 
397 
399 
400 
401 template<typename Index32LeafT>
403 {
404 public:
405  typedef std::pair<size_t, size_t> IndexRange;
406 
408  std::vector<Vec3R>& instancePoints,
409  std::vector<float>& instanceDistances,
410  const PointList& surfacePointList,
411  const std::vector<const Index32LeafT*>& leafNodes,
412  const std::vector<IndexRange>& leafRanges,
413  const std::vector<Vec4R>& leafBoundingSpheres,
414  const std::vector<Vec4R>& nodeBoundingSpheres,
415  size_t maxNodeLeafs,
416  bool transformPoints = false);
417 
418 
419  void run(bool threaded = true);
420 
421 
422  void operator()(const tbb::blocked_range<size_t>&) const;
423 
424 private:
425 
426  void evalLeaf(size_t index, const Index32LeafT& leaf) const;
427  void evalNode(size_t pointIndex, size_t nodeIndex) const;
428 
429 
430  std::vector<Vec3R>& mInstancePoints;
431  std::vector<float>& mInstanceDistances;
432 
433  const PointList& mSurfacePointList;
434 
435  const std::vector<const Index32LeafT*>& mLeafNodes;
436  const std::vector<IndexRange>& mLeafRanges;
437  const std::vector<Vec4R>& mLeafBoundingSpheres;
438  const std::vector<Vec4R>& mNodeBoundingSpheres;
439 
440  std::vector<float> mLeafDistances, mNodeDistances;
441 
442  const bool mTransformPoints;
443  size_t mClosestPointIndex;
444 };
445 
446 
447 template<typename Index32LeafT>
449  std::vector<Vec3R>& instancePoints,
450  std::vector<float>& instanceDistances,
451  const PointList& surfacePointList,
452  const std::vector<const Index32LeafT*>& leafNodes,
453  const std::vector<IndexRange>& leafRanges,
454  const std::vector<Vec4R>& leafBoundingSpheres,
455  const std::vector<Vec4R>& nodeBoundingSpheres,
456  size_t maxNodeLeafs,
457  bool transformPoints)
458  : mInstancePoints(instancePoints)
459  , mInstanceDistances(instanceDistances)
460  , mSurfacePointList(surfacePointList)
461  , mLeafNodes(leafNodes)
462  , mLeafRanges(leafRanges)
463  , mLeafBoundingSpheres(leafBoundingSpheres)
464  , mNodeBoundingSpheres(nodeBoundingSpheres)
465  , mLeafDistances(maxNodeLeafs, 0.0)
466  , mNodeDistances(leafRanges.size(), 0.0)
467  , mTransformPoints(transformPoints)
468  , mClosestPointIndex(0)
469 {
470 }
471 
472 
473 template<typename Index32LeafT>
474 void
476 {
477  if (threaded) {
478  tbb::parallel_for(tbb::blocked_range<size_t>(0, mInstancePoints.size()), *this);
479  } else {
480  (*this)(tbb::blocked_range<size_t>(0, mInstancePoints.size()));
481  }
482 }
483 
484 template<typename Index32LeafT>
485 void
486 ClosestPointDist<Index32LeafT>::evalLeaf(size_t index, const Index32LeafT& leaf) const
487 {
488  typename Index32LeafT::ValueOnCIter iter;
489  const Vec3s center = mInstancePoints[index];
490  size_t& closestPointIndex = const_cast<size_t&>(mClosestPointIndex);
491 
492  for (iter = leaf.cbeginValueOn(); iter; ++iter) {
493 
494  const Vec3s& point = mSurfacePointList[iter.getValue()];
495  float tmpDist = (point - center).lengthSqr();
496 
497  if (tmpDist < mInstanceDistances[index]) {
498  mInstanceDistances[index] = tmpDist;
499  closestPointIndex = iter.getValue();
500  }
501  }
502 }
503 
504 
505 template<typename Index32LeafT>
506 void
507 ClosestPointDist<Index32LeafT>::evalNode(size_t pointIndex, size_t nodeIndex) const
508 {
509  const Vec3R& pos = mInstancePoints[pointIndex];
510  float minDist = mInstanceDistances[pointIndex];
511  size_t minDistIdx = 0;
512  Vec3R center;
513  bool updatedDist = false;
514 
515  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
516  i < mLeafRanges[nodeIndex].second; ++i, ++n)
517  {
518  float& distToLeaf = const_cast<float&>(mLeafDistances[n]);
519 
520  center[0] = mLeafBoundingSpheres[i][0];
521  center[1] = mLeafBoundingSpheres[i][1];
522  center[2] = mLeafBoundingSpheres[i][2];
523 
524  distToLeaf = float((pos - center).lengthSqr() - mLeafBoundingSpheres[i][3]);
525 
526  if (distToLeaf < minDist) {
527  minDist = distToLeaf;
528  minDistIdx = i;
529  updatedDist = true;
530  }
531  }
532 
533  if (!updatedDist) return;
534 
535  evalLeaf(pointIndex, *mLeafNodes[minDistIdx]);
536 
537  for (size_t i = mLeafRanges[nodeIndex].first, n = 0;
538  i < mLeafRanges[nodeIndex].second; ++i, ++n)
539  {
540  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
541  evalLeaf(pointIndex, *mLeafNodes[i]);
542  }
543  }
544 }
545 
546 
547 template<typename Index32LeafT>
548 void
549 ClosestPointDist<Index32LeafT>::operator()(const tbb::blocked_range<size_t>& range) const
550 {
551  Vec3R center;
552  for (size_t n = range.begin(); n != range.end(); ++n) {
553 
554  const Vec3R& pos = mInstancePoints[n];
555  float minDist = mInstanceDistances[n];
556  size_t minDistIdx = 0;
557 
558  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
559  float& distToNode = const_cast<float&>(mNodeDistances[i]);
560 
561  center[0] = mNodeBoundingSpheres[i][0];
562  center[1] = mNodeBoundingSpheres[i][1];
563  center[2] = mNodeBoundingSpheres[i][2];
564 
565  distToNode = float((pos - center).lengthSqr() - mNodeBoundingSpheres[i][3]);
566 
567  if (distToNode < minDist) {
568  minDist = distToNode;
569  minDistIdx = i;
570  }
571  }
572 
573  evalNode(n, minDistIdx);
574 
575  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
576  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
577  evalNode(n, i);
578  }
579  }
580 
581  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
582 
583  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
584  }
585 }
586 
587 
589 {
590 public:
591  UpdatePoints(
592  const Vec4s& sphere,
593  const std::vector<Vec3R>& points,
594  std::vector<float>& distances,
595  std::vector<unsigned char>& mask,
596  bool overlapping);
597 
598  float radius() const { return mRadius; }
599  int index() const { return mIndex; }
600 
601  inline void run(bool threaded = true);
602 
603 
604  UpdatePoints(UpdatePoints&, tbb::split);
605  inline void operator()(const tbb::blocked_range<size_t>& range);
606  void join(const UpdatePoints& rhs)
607  {
608  if (rhs.mRadius > mRadius) {
609  mRadius = rhs.mRadius;
610  mIndex = rhs.mIndex;
611  }
612  }
613 
614 private:
615 
616  const Vec4s& mSphere;
617  const std::vector<Vec3R>& mPoints;
618 
619  std::vector<float>& mDistances;
620  std::vector<unsigned char>& mMask;
621 
622  bool mOverlapping;
623  float mRadius;
624  int mIndex;
625 };
626 
627 inline
629  const Vec4s& sphere,
630  const std::vector<Vec3R>& points,
631  std::vector<float>& distances,
632  std::vector<unsigned char>& mask,
633  bool overlapping)
634  : mSphere(sphere)
635  , mPoints(points)
636  , mDistances(distances)
637  , mMask(mask)
638  , mOverlapping(overlapping)
639  , mRadius(0.0)
640  , mIndex(0)
641 {
642 }
643 
644 inline
646  : mSphere(rhs.mSphere)
647  , mPoints(rhs.mPoints)
648  , mDistances(rhs.mDistances)
649  , mMask(rhs.mMask)
650  , mOverlapping(rhs.mOverlapping)
651  , mRadius(rhs.mRadius)
652  , mIndex(rhs.mIndex)
653 {
654 }
655 
656 inline void
657 UpdatePoints::run(bool threaded)
658 {
659  if (threaded) {
660  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
661  } else {
662  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
663  }
664 }
665 
666 inline void
667 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
668 {
669  Vec3s pos;
670  for (size_t n = range.begin(); n != range.end(); ++n) {
671  if (mMask[n]) continue;
672 
673  pos.x() = float(mPoints[n].x()) - mSphere[0];
674  pos.y() = float(mPoints[n].y()) - mSphere[1];
675  pos.z() = float(mPoints[n].z()) - mSphere[2];
676 
677  float dist = pos.length();
678 
679  if (dist < mSphere[3]) {
680  mMask[n] = 1;
681  continue;
682  }
683 
684  if (!mOverlapping) {
685  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
686  }
687 
688  if (mDistances[n] > mRadius) {
689  mRadius = mDistances[n];
690  mIndex = int(n);
691  }
692  }
693 }
694 
695 
696 } // namespace internal
697 
698 
700 
701 
702 template<typename GridT, typename InterrupterT>
703 inline void
705  const GridT& grid,
706  std::vector<openvdb::Vec4s>& spheres,
707  int maxSphereCount,
708  bool overlapping,
709  float minRadius,
710  float maxRadius,
711  float isovalue,
712  int instanceCount,
713  InterrupterT* interrupter)
714 {
715  spheres.clear();
716  spheres.reserve(maxSphereCount);
717 
718  const bool addNBPoints = grid.activeVoxelCount() < 10000;
719  int instances = std::max(instanceCount, maxSphereCount);
720 
721  typedef typename GridT::TreeType TreeT;
722  typedef typename GridT::ValueType ValueT;
723 
724  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
725  typedef typename TreeT::template ValueConverter<Index32>::Type Index32TreeT;
726  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
727 
728  typedef boost::mt11213b RandGen;
729  RandGen mtRand(/*seed=*/0);
730 
731  const TreeT& tree = grid.tree();
732  const math::Transform& transform = grid.transform();
733 
734  std::vector<Vec3R> instancePoints;
735 
736  { // Scatter candidate sphere centroids (instancePoints)
737  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
738 
739  if (grid.getGridClass() == GRID_LEVEL_SET) {
740  interiorMaskPtr = sdfInteriorMask(grid, ValueT(isovalue));
741  } else {
742  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
743  interiorMaskPtr->setTransform(transform.copy());
744  interiorMaskPtr->tree().topologyUnion(tree);
745  }
746 
747  if (interrupter && interrupter->wasInterrupted()) return;
748 
749  erodeVoxels(interiorMaskPtr->tree(), 1);
750 
751  instancePoints.reserve(instances);
752  internal::PointAccessor ptnAcc(instancePoints);
753 
755  ptnAcc, Index64(addNBPoints ? (instances / 2) : instances), mtRand, interrupter);
756 
757  scatter(*interiorMaskPtr);
758  }
759 
760  if (interrupter && interrupter->wasInterrupted()) return;
761 
762  std::vector<float> instanceRadius;
763 
765  csp.initialize(grid, isovalue, interrupter);
766 
767  // add extra instance points in the interior narrow band.
768  if (instancePoints.size() < size_t(instances)) {
769  const Int16TreeT& signTree = csp.signTree();
770  typename Int16TreeT::LeafNodeType::ValueOnCIter it;
771  typename Int16TreeT::LeafCIter leafIt = signTree.cbeginLeaf();
772 
773  for (; leafIt; ++leafIt) {
774  for (it = leafIt->cbeginValueOn(); it; ++it) {
775 
776  const int flags = int(it.getValue());
778  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
779  }
780 
781  if (instancePoints.size() == size_t(instances)) break;
782  }
783  if (instancePoints.size() == size_t(instances)) break;
784  }
785  }
786 
787 
788  if (interrupter && interrupter->wasInterrupted()) return;
789 
790  if (!csp.search(instancePoints, instanceRadius)) return;
791 
792  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
793  float largestRadius = 0.0;
794  int largestRadiusIdx = 0;
795 
796  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
797  if (instanceRadius[n] > largestRadius) {
798  largestRadius = instanceRadius[n];
799  largestRadiusIdx = int(n);
800  }
801  }
802 
803  Vec3s pos;
804  Vec4s sphere;
805  minRadius = float(minRadius * transform.voxelSize()[0]);
806  maxRadius = float(maxRadius * transform.voxelSize()[0]);
807 
808  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
809 
810  if (interrupter && interrupter->wasInterrupted()) return;
811 
812  largestRadius = std::min(maxRadius, largestRadius);
813 
814  if (s != 0 && largestRadius < minRadius) break;
815 
816  sphere[0] = float(instancePoints[largestRadiusIdx].x());
817  sphere[1] = float(instancePoints[largestRadiusIdx].y());
818  sphere[2] = float(instancePoints[largestRadiusIdx].z());
819  sphere[3] = largestRadius;
820 
821  spheres.push_back(sphere);
822  instanceMask[largestRadiusIdx] = 1;
823 
825  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
826  op.run();
827 
828  largestRadius = op.radius();
829  largestRadiusIdx = op.index();
830  }
831 }
832 
834 
835 
836 template<typename GridT>
838  : mIsInitialized(false)
839  , mLeafBoundingSpheres(0)
840  , mNodeBoundingSpheres(0)
841  , mLeafRanges(0)
842  , mLeafNodes(0)
843  , mSurfacePointList()
844  , mPointListSize(0)
845  , mMaxNodeLeafs(0)
846  , mMaxRadiusSqr(0.0)
847  , mIdxTreePt()
848 {
849 }
850 
851 template<typename GridT>
852 void
853 ClosestSurfacePoint<GridT>::initialize(const GridT& grid, float isovalue)
854 {
855  initialize<GridT, util::NullInterrupter>(grid, isovalue, NULL);
856 }
857 
858 
859 template<typename GridT>
860 template<typename InterrupterT>
861 void
863  const GridT& grid, float isovalue, InterrupterT* interrupter)
864 {
865  mIsInitialized = false;
866  typedef tree::LeafManager<const TreeT> LeafManagerT;
867  typedef tree::LeafManager<Index32TreeT> Index32LeafManagerT;
868  typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
869  typedef typename GridT::ValueType ValueT;
870 
871  const TreeT& tree = grid.tree();
872  const math::Transform& transform = grid.transform();
873 
874  { // Extract surface point cloud
875 
876  mSignTreePt.reset(new Int16TreeT(0));
877  mIdxTreePt.reset(new Index32TreeT(boost::integer_traits<Index32>::const_max));
878 
879  BoolTreeT mask(false);
881  volume_to_mesh_internal::computeAuxiliaryData(*mSignTreePt, *mIdxTreePt, mask, tree, ValueT(isovalue));
882 
883  if (interrupter && interrupter->wasInterrupted()) return;
884 
885  // count unique points
886 
887  typedef typename Int16TreeT::LeafNodeType Int16LeafNodeType;
888  typedef typename Index32TreeT::LeafNodeType Index32LeafNodeType;
889 
890  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
891  mSignTreePt->getNodes(signFlagsLeafNodes);
892 
893  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
894 
895  boost::scoped_array<Index32> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
896 
897  tbb::parallel_for(auxiliaryLeafNodeRange,
899  (signFlagsLeafNodes, leafNodeOffsets));
900 
901  {
902  Index32 pointCount = 0;
903  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
904  const Index32 tmp = leafNodeOffsets[n];
905  leafNodeOffsets[n] = pointCount;
906  pointCount += tmp;
907  }
908 
909  mPointListSize = size_t(pointCount);
910  mSurfacePointList.reset(new Vec3s[mPointListSize]);
911  }
912 
913 
914  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
915  mIdxTreePt->getNodes(pointIndexLeafNodes);
916 
917  tbb::parallel_for(auxiliaryLeafNodeRange,
918  volume_to_mesh_internal::ComputePoints<TreeT>(mSurfacePointList.get(), tree, pointIndexLeafNodes,
919  signFlagsLeafNodes, leafNodeOffsets, transform, ValueT(isovalue)));
920  }
921 
922  if (interrupter && interrupter->wasInterrupted()) return;
923 
924  // estimate max sphere radius (sqr dist)
925  CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
926 
927  Vec3s dim = transform.indexToWorld(bbox.min()) -
928  transform.indexToWorld(bbox.max());
929 
930  dim[0] = std::abs(dim[0]);
931  dim[1] = std::abs(dim[1]);
932  dim[2] = std::abs(dim[2]);
933 
934  mMaxRadiusSqr = std::min(std::min(dim[0], dim[1]), dim[2]);
935  mMaxRadiusSqr *= 0.51f;
936  mMaxRadiusSqr *= mMaxRadiusSqr;
937 
938 
939  Index32LeafManagerT idxLeafs(*mIdxTreePt);
940 
941  typedef typename Index32TreeT::RootNodeType Index32RootNodeT;
942  typedef typename Index32RootNodeT::NodeChainType Index32NodeChainT;
943  BOOST_STATIC_ASSERT(boost::mpl::size<Index32NodeChainT>::value > 1);
944  typedef typename boost::mpl::at<Index32NodeChainT, boost::mpl::int_<1> >::type Index32InternalNodeT;
945 
946  typename Index32TreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
947  nIt.setMinDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
948  nIt.setMaxDepth(Index32TreeT::NodeCIter::LEAF_DEPTH - 1);
949 
950  std::vector<const Index32InternalNodeT*> internalNodes;
951 
952  const Index32InternalNodeT* node = NULL;
953  for (; nIt; ++nIt) {
954  nIt.getNode(node);
955  if (node) internalNodes.push_back(node);
956  }
957 
958  std::vector<IndexRange>().swap(mLeafRanges);
959  mLeafRanges.resize(internalNodes.size());
960 
961  std::vector<const Index32LeafT*>().swap(mLeafNodes);
962  mLeafNodes.reserve(idxLeafs.leafCount());
963 
964  typename Index32InternalNodeT::ChildOnCIter leafIt;
965  mMaxNodeLeafs = 0;
966  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
967 
968  mLeafRanges[n].first = mLeafNodes.size();
969 
970  size_t leafCount = 0;
971  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
972  mLeafNodes.push_back(&(*leafIt));
973  ++leafCount;
974  }
975 
976  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
977 
978  mLeafRanges[n].second = mLeafNodes.size();
979  }
980 
981  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
982  mLeafBoundingSpheres.resize(mLeafNodes.size());
983 
985  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
986  leafBS.run();
987 
988 
989  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
990  mNodeBoundingSpheres.resize(internalNodes.size());
991 
992  internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
993  nodeBS.run();
994  mIsInitialized = true;
995 }
996 
997 
998 template<typename GridT>
999 bool
1000 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
1001  std::vector<float>& distances, bool transformPoints)
1002 {
1003  if (!mIsInitialized) return false;
1004 
1005  distances.clear();
1006  distances.resize(points.size(), mMaxRadiusSqr);
1007 
1008  internal::ClosestPointDist<Index32LeafT> cpd(points, distances, mSurfacePointList,
1009  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1010  mMaxNodeLeafs, transformPoints);
1011 
1012  cpd.run();
1013 
1014  return true;
1015 }
1016 
1017 
1018 template<typename GridT>
1019 bool
1020 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
1021 {
1022  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
1023 }
1024 
1025 
1026 template<typename GridT>
1027 bool
1029  std::vector<float>& distances)
1030 {
1031  return search(points, distances, true);
1032 }
1033 
1034 
1035 } // namespace tools
1036 } // namespace OPENVDB_VERSION_NAME
1037 } // namespace openvdb
1038 
1039 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
1040 
1041 // Copyright (c) 2012-2016 DreamWorks Animation LLC
1042 // All rights reserved. This software is distributed under the
1043 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
ClosestSurfacePoint()
Definition: VolumeToSpheres.h:837
Definition: VolumeToSpheres.h:588
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)
fillWithSpheres method variant that automatically infers the util::NullInterrupter.
Definition: VolumeToSpheres.h:107
Vec3< float > Vec3s
Definition: Vec3.h:650
void operator()(const tbb::blocked_range< size_t > &range)
Definition: VolumeToSpheres.h:667
void join(const UpdatePoints &rhs)
Definition: VolumeToSpheres.h:606
void run(bool threaded=true)
Definition: VolumeToSpheres.h:345
NodeBS(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
Definition: VolumeToSpheres.h:335
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:872
void run(bool threaded=true)
Definition: VolumeToSpheres.h:475
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:549
const Int16TreeT & signTree() const
Tree accessors.
Definition: VolumeToSpheres.h:184
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:115
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Computes distance to closest surface.
Definition: VolumeToSpheres.h:1020
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:112
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
Definition: VolumeToSpheres.h:315
GridT::TreeType TreeT
Definition: VolumeToSpheres.h:132
void run(bool threaded=true)
Definition: VolumeToSpheres.h:657
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
Ptr copy() const
Definition: Transform.h:77
TreeType & tree()
Return a reference to this grid&#39;s tree, which might be shared with other grids.
Definition: Grid.h:761
TreeT::template ValueConverter< Index32 >::Type Index32TreeT
Definition: VolumeToSpheres.h:134
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:338
Definition: VolumeToSpheres.h:232
void run(bool threaded=true)
Definition: VolumeToSpheres.h:268
uint64_t Index64
Definition: Types.h:57
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3< double > Vec3d
Definition: Vec3.h:651
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToSpheres.h:135
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:3883
boost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:182
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:318
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:355
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
const Index32TreeT & indexTree() const
Tree accessors.
Definition: VolumeToSpheres.h:183
Vec4< float > Vec4s
Definition: Vec4.h:583
Definition: Exceptions.h:39
uint32_t Index32
Definition: Types.h:56
OPENVDB_IMPORT void initialize()
Global registration of basic types.
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
Definition: VolumeToSpheres.h:628
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:2244
int index() const
Definition: VolumeToSpheres.h:599
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:405
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:448
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:279
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:135
void initialize(const GridT &grid, float isovalue=0.0, InterrupterT *interrupter=NULL)
Extracts the surface points and constructs a spatial acceleration structure.
Definition: VolumeToSpheres.h:862
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=NULL)
Threaded method to fill a closed level set or fog volume with adaptively sized spheres.
Definition: VolumeToSpheres.h:704
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Performs closest point searches.
Definition: VolumeToSpheres.h:1028
Definition: Types.h:213
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
void identifySurfaceIntersectingVoxels(typename InputTreeType::template ValueConverter< bool >::Type &intersectionTree, const InputTreeType &inputTree, typename InputTreeType::ValueType isovalue)
Definition: VolumeToMesh.h:3340
void add(const Vec3R &pos)
Definition: VolumeToSpheres.h:222
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:1039
Accelerated closest surface point queries for narrow band level sets. Supports queries that originate...
Definition: VolumeToSpheres.h:129
PointAccessor(std::vector< Vec3R > &points)
Definition: VolumeToSpheres.h:217
float radius() const
Definition: VolumeToSpheres.h:598
Definition: VolumeToSpheres.h:215
TreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: VolumeToSpheres.h:133
Definition: Mat4.h:51
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:128