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