OpenVDB  2.3.0
VolumeToSpheres.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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  void run(bool threaded = true);
324 
325 
326  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  Vec3s 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 = mLeafRanges[n].second - 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  float 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  float 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 IntLeafT>
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 IntLeafT*>& 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 IntLeafT& 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 IntLeafT*>& 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 IntLeafT>
449  std::vector<Vec3R>& instancePoints,
450  std::vector<float>& instanceDistances,
451  const PointList& surfacePointList,
452  const std::vector<const IntLeafT*>& 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 IntLeafT>
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 IntLeafT>
485 void
486 ClosestPointDist<IntLeafT>::evalLeaf(size_t index, const IntLeafT& leaf) const
487 {
488  typename IntLeafT::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 IntLeafT>
506 void
507 ClosestPointDist<IntLeafT>::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; 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 = (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; i < mLeafRanges[nodeIndex].second; ++i, ++n) {
537  if (mLeafDistances[n] < mInstanceDistances[pointIndex] && i != minDistIdx) {
538  evalLeaf(pointIndex, *mLeafNodes[i]);
539  }
540  }
541 }
542 
543 
544 template<typename IntLeafT>
545 void
546 ClosestPointDist<IntLeafT>::operator()(const tbb::blocked_range<size_t>& range) const
547 {
548  Vec3R center;
549  for (size_t n = range.begin(); n != range.end(); ++n) {
550 
551  const Vec3R& pos = mInstancePoints[n];
552  float minDist = mInstanceDistances[n];
553  size_t minDistIdx = 0;
554 
555  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
556  float& distToNode = const_cast<float&>(mNodeDistances[i]);
557 
558  center[0] = mNodeBoundingSpheres[i][0];
559  center[1] = mNodeBoundingSpheres[i][1];
560  center[2] = mNodeBoundingSpheres[i][2];
561 
562  distToNode = (pos - center).lengthSqr() - mNodeBoundingSpheres[i][3];
563 
564  if (distToNode < minDist) {
565  minDist = distToNode;
566  minDistIdx = i;
567  }
568  }
569 
570  evalNode(n, minDistIdx);
571 
572  for (size_t i = 0, I = mNodeDistances.size(); i < I; ++i) {
573  if (mNodeDistances[i] < mInstanceDistances[n] && i != minDistIdx) {
574  evalNode(n, i);
575  }
576  }
577 
578  mInstanceDistances[n] = std::sqrt(mInstanceDistances[n]);
579 
580  if (mTransformPoints) mInstancePoints[n] = mSurfacePointList[mClosestPointIndex];
581  }
582 }
583 
584 
586 {
587 public:
588  UpdatePoints(
589  const Vec4s& sphere,
590  const std::vector<Vec3R>& points,
591  std::vector<float>& distances,
592  std::vector<unsigned char>& mask,
593  bool overlapping);
594 
595  float radius() const { return mRadius; }
596  int index() const { return mIndex; };
597 
598  void run(bool threaded = true);
599 
600 
601  UpdatePoints(UpdatePoints&, tbb::split);
602  void operator()(const tbb::blocked_range<size_t>& range);
603  void join(const UpdatePoints& rhs)
604  {
605  if (rhs.mRadius > mRadius) {
606  mRadius = rhs.mRadius;
607  mIndex = rhs.mIndex;
608  }
609  }
610 
611 private:
612 
613  const Vec4s& mSphere;
614  const std::vector<Vec3R>& mPoints;
615 
616  std::vector<float>& mDistances;
617  std::vector<unsigned char>& mMask;
618 
619  bool mOverlapping;
620  float mRadius;
621  int mIndex;
622 };
623 
624 inline
626  const Vec4s& sphere,
627  const std::vector<Vec3R>& points,
628  std::vector<float>& distances,
629  std::vector<unsigned char>& mask,
630  bool overlapping)
631  : mSphere(sphere)
632  , mPoints(points)
633  , mDistances(distances)
634  , mMask(mask)
635  , mOverlapping(overlapping)
636  , mRadius(0.0)
637  , mIndex(0)
638 {
639 }
640 
641 inline
643  : mSphere(rhs.mSphere)
644  , mPoints(rhs.mPoints)
645  , mDistances(rhs.mDistances)
646  , mMask(rhs.mMask)
647  , mOverlapping(rhs.mOverlapping)
648  , mRadius(rhs.mRadius)
649  , mIndex(rhs.mIndex)
650 {
651 }
652 
653 inline void
654 UpdatePoints::run(bool threaded)
655 {
656  if (threaded) {
657  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPoints.size()), *this);
658  } else {
659  (*this)(tbb::blocked_range<size_t>(0, mPoints.size()));
660  }
661 }
662 
663 inline void
664 UpdatePoints::operator()(const tbb::blocked_range<size_t>& range)
665 {
666  Vec3s pos;
667  for (size_t n = range.begin(); n != range.end(); ++n) {
668  if (mMask[n]) continue;
669 
670  pos.x() = float(mPoints[n].x()) - mSphere[0];
671  pos.y() = float(mPoints[n].y()) - mSphere[1];
672  pos.z() = float(mPoints[n].z()) - mSphere[2];
673 
674  float dist = pos.length();
675 
676  if (dist < mSphere[3]) {
677  mMask[n] = 1;
678  continue;
679  }
680 
681  if (!mOverlapping) {
682  mDistances[n] = std::min(mDistances[n], (dist - mSphere[3]));
683  }
684 
685  if (mDistances[n] > mRadius) {
686  mRadius = mDistances[n];
687  mIndex = n;
688  }
689  }
690 }
691 
692 
693 } // namespace internal
694 
695 
697 
698 
699 template<typename GridT, typename InterrupterT>
700 inline void
702  const GridT& grid,
703  std::vector<openvdb::Vec4s>& spheres,
704  int maxSphereCount,
705  bool overlapping,
706  float minRadius,
707  float maxRadius,
708  float isovalue,
709  int instanceCount,
710  InterrupterT* interrupter)
711 {
712  spheres.clear();
713  spheres.reserve(maxSphereCount);
714 
715  const bool addNBPoints = grid.activeVoxelCount() < 10000;
716  int instances = std::max(instanceCount, maxSphereCount);
717 
718  typedef typename GridT::TreeType TreeT;
719  typedef typename GridT::ValueType ValueT;
720 
721  typedef typename TreeT::template ValueConverter<bool>::Type BoolTreeT;
722  typedef typename TreeT::template ValueConverter<int>::Type IntTreeT;
723  typedef typename TreeT::template ValueConverter<Int16>::Type Int16TreeT;
724 
725  typedef tree::LeafManager<const TreeT> LeafManagerT;
726  typedef tree::LeafManager<IntTreeT> IntLeafManagerT;
727  typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
728 
729 
730  typedef boost::mt11213b RandGen;
731  RandGen mtRand(/*seed=*/0);
732 
733  const TreeT& tree = grid.tree();
734  const math::Transform& transform = grid.transform();
735 
736  std::vector<Vec3R> instancePoints;
737 
738  { // Scatter candidate sphere centroids (instancePoints)
739  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
740 
741  if (grid.getGridClass() == GRID_LEVEL_SET) {
742  interiorMaskPtr = sdfInteriorMask(grid, ValueT(isovalue));
743  } else {
744  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
745  interiorMaskPtr->setTransform(transform.copy());
746  interiorMaskPtr->tree().topologyUnion(tree);
747  }
748 
749  if (interrupter && interrupter->wasInterrupted()) return;
750 
751  erodeVoxels(interiorMaskPtr->tree(), 1);
752 
753  instancePoints.reserve(instances);
754  internal::PointAccessor ptnAcc(instancePoints);
755 
757  scatter(ptnAcc, (addNBPoints ? (instances / 2) : instances), mtRand, interrupter);
758 
759  scatter(*interiorMaskPtr);
760  }
761 
762  if (interrupter && interrupter->wasInterrupted()) return;
763 
764  std::vector<float> instanceRadius;
765 
767  csp.initialize(grid, isovalue, interrupter);
768 
769  // add extra instance points in the interior narrow band.
770  if (instancePoints.size() < instances) {
771  const Int16TreeT& signTree = csp.signTree();
772  typename Int16TreeT::LeafNodeType::ValueOnCIter it;
773  typename Int16TreeT::LeafCIter leafIt = signTree.cbeginLeaf();
774 
775  for (; leafIt; ++leafIt) {
776  for (it = leafIt->cbeginValueOn(); it; ++it) {
777  const int flags = it.getValue();
778  if (!(0xE00 & flags) && (flags & 0x100)) {
779  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
780  }
781 
782  if (instancePoints.size() == instances) break;
783  }
784  if (instancePoints.size() == instances) break;
785  }
786  }
787 
788 
789  if (interrupter && interrupter->wasInterrupted()) return;
790 
791  if (!csp.search(instancePoints, instanceRadius)) return;
792 
793  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
794  float largestRadius = 0.0;
795  int largestRadiusIdx = 0;
796 
797  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
798  if (instanceRadius[n] > largestRadius) {
799  largestRadius = instanceRadius[n];
800  largestRadiusIdx = n;
801  }
802  }
803 
804  Vec3s pos;
805  Vec4s sphere;
806  minRadius *= transform.voxelSize()[0];
807  maxRadius *= transform.voxelSize()[0];
808 
809  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
810 
811  if (interrupter && interrupter->wasInterrupted()) return;
812 
813  largestRadius = std::min(maxRadius, largestRadius);
814 
815  if (s != 0 && largestRadius < minRadius) break;
816 
817  sphere[0] = float(instancePoints[largestRadiusIdx].x());
818  sphere[1] = float(instancePoints[largestRadiusIdx].y());
819  sphere[2] = float(instancePoints[largestRadiusIdx].z());
820  sphere[3] = largestRadius;
821 
822  spheres.push_back(sphere);
823  instanceMask[largestRadiusIdx] = 1;
824 
825  internal::UpdatePoints op(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<IntTreeT> IntLeafManagerT;
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  {
877  LeafManagerT leafs(tree);
879  signDataOp(tree, leafs, ValueT(isovalue));
880 
881  signDataOp.run();
882 
883  mSignTreePt = signDataOp.signTree();
884  mIdxTreePt = signDataOp.idxTree();
885  }
886 
887  if (interrupter && interrupter->wasInterrupted()) return;
888 
889  Int16LeafManagerT signLeafs(*mSignTreePt);
890 
891  std::vector<size_t> regions(signLeafs.leafCount(), 0);
892  signLeafs.foreach(internal::CountPoints(regions));
893 
894  mPointListSize = 0;
895  for (size_t tmp = 0, n = 0, N = regions.size(); n < N; ++n) {
896  tmp = regions[n];
897  regions[n] = mPointListSize;
898  mPointListSize += tmp;
899  }
900 
901  if (mPointListSize == 0) return;
902 
903  mSurfacePointList.reset(new Vec3s[mPointListSize]);
904 
906  pointOp(signLeafs, tree, *mIdxTreePt, mSurfacePointList, regions, transform, isovalue);
907 
908  pointOp.run();
909 
910  mIdxTreePt->topologyUnion(*mSignTreePt);
911  }
912 
913  if (interrupter && interrupter->wasInterrupted()) return;
914 
915  // estimate max sphere radius (sqr dist)
916  CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
917 
918  Vec3s dim = transform.indexToWorld(bbox.min()) -
919  transform.indexToWorld(bbox.max());
920 
921  dim[0] = std::abs(dim[0]);
922  dim[1] = std::abs(dim[1]);
923  dim[2] = std::abs(dim[2]);
924 
925  mMaxRadiusSqr = std::min(std::min(dim[0], dim[1]), dim[2]);
926  mMaxRadiusSqr *= 0.51;
927  mMaxRadiusSqr *= mMaxRadiusSqr;
928 
929 
930  IntLeafManagerT idxLeafs(*mIdxTreePt);
931 
932 
933  typedef typename IntTreeT::RootNodeType IntRootNodeT;
934  typedef typename IntRootNodeT::NodeChainType IntNodeChainT;
935  BOOST_STATIC_ASSERT(boost::mpl::size<IntNodeChainT>::value > 1);
936  typedef typename boost::mpl::at<IntNodeChainT, boost::mpl::int_<1> >::type IntInternalNodeT;
937 
938 
939  typename IntTreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
940  nIt.setMinDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
941  nIt.setMaxDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
942 
943  std::vector<const IntInternalNodeT*> internalNodes;
944 
945  const IntInternalNodeT* node = NULL;
946  for (; nIt; ++nIt) {
947  nIt.getNode(node);
948  if (node) internalNodes.push_back(node);
949  }
950 
951  std::vector<IndexRange>().swap(mLeafRanges);
952  mLeafRanges.resize(internalNodes.size());
953 
954  std::vector<const IntLeafT*>().swap(mLeafNodes);
955  mLeafNodes.reserve(idxLeafs.leafCount());
956 
957  typename IntInternalNodeT::ChildOnCIter leafIt;
958  mMaxNodeLeafs = 0;
959  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
960 
961  mLeafRanges[n].first = mLeafNodes.size();
962 
963  size_t leafCount = 0;
964  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
965  mLeafNodes.push_back(&(*leafIt));
966  ++leafCount;
967  }
968 
969  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
970 
971  mLeafRanges[n].second = mLeafNodes.size();
972  }
973 
974  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
975  mLeafBoundingSpheres.resize(mLeafNodes.size());
976 
977  internal::LeafBS<IntLeafT> leafBS(mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
978  leafBS.run();
979 
980 
981  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
982  mNodeBoundingSpheres.resize(internalNodes.size());
983 
984  internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
985  nodeBS.run();
986  mIsInitialized = true;
987 }
988 
989 
990 template<typename GridT>
991 bool
992 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
993  std::vector<float>& distances, bool transformPoints)
994 {
995  if (!mIsInitialized) return false;
996 
997  distances.clear();
998  distances.resize(points.size(), mMaxRadiusSqr);
999 
1000  internal::ClosestPointDist<IntLeafT> cpd(points, distances, mSurfacePointList,
1001  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1002  mMaxNodeLeafs, transformPoints);
1003 
1004  cpd.run();
1005 
1006  return true;
1007 }
1008 
1009 
1010 template<typename GridT>
1011 bool
1012 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
1013 {
1014  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
1015 }
1016 
1017 
1018 template<typename GridT>
1019 bool
1020 ClosestSurfacePoint<GridT>::searchAndReplace(std::vector<Vec3R>& points, 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-2013 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/ )
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:546
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:317
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToSpheres.h:133
Definition: VolumeToMesh.h:894
OPENVDB_IMPORT void initialize()
Global registration of basic types.
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:978
boost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:167
void add(const Vec3R &pos)
Definition: VolumeToSpheres.h:221
void join(const UpdatePoints &rhs)
Definition: VolumeToSpheres.h:603
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)
fillWithSpheres method variant that automatically infers the util::NullInterrupter.
Definition: VolumeToSpheres.h:107
Definition: VolumeToSpheres.h:314
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:448
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:278
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToSpheres.h:134
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:710
void run(bool threaded=true)
Definition: VolumeToMesh.h:974
IntTreeT::Ptr idxTree() const
Definition: VolumeToMesh.h:914
void run(bool threaded=true)
Definition: VolumeToSpheres.h:475
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Computes distance to closest surface.
Definition: VolumeToSpheres.h:1012
Counts the total number of points per leaf, accounts for cells with multiple points.
Definition: VolumeToMesh.h:1051
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:461
Definition: Types.h:170
Definition: VolumeToSpheres.h:214
Grid< typename GridType::TreeType::template ValueConverter< bool >::Type >::Ptr sdfInteriorMask(const GridType &grid, typename GridType::ValueType iso=lsutilGridZero< GridType >())
Threaded method to extract an interior region mask from a level set/SDF grid.
Definition: LevelSetUtil.h:357
const IntTreeT & indexTree() const
Tree accessors.
Definition: VolumeToSpheres.h:182
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
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:701
#define OPENVDB_VERSION_NAME
Definition: version.h:45
Int16TreeT::Ptr signTree() const
Definition: VolumeToMesh.h:913
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Performs closest point searches.
Definition: VolumeToSpheres.h:1020
NodeBS(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
Definition: VolumeToSpheres.h:335
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
float radius() const
Definition: VolumeToSpheres.h:595
ClosestSurfacePoint()
Definition: VolumeToSpheres.h:837
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
Definition: VolumeToSpheres.h:625
Accelerated closest surface point queries for narrow band level sets. Supports queries that originate...
Definition: VolumeToSpheres.h:129
Ptr copy() const
Definition: Transform.h:76
PointAccessor(std::vector< Vec3R > &points)
Definition: VolumeToSpheres.h:216
void operator()(const tbb::blocked_range< size_t > &range)
Definition: VolumeToSpheres.h:664
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:355
GridT::TreeType TreeT
Definition: VolumeToSpheres.h:132
const Int16TreeT & signTree() const
Tree accessors.
Definition: VolumeToSpheres.h:183
Vec3< float > Vec3s
Definition: Vec3.h:624
void run(bool threaded=true)
Definition: VolumeToSpheres.h:345
TreeType & tree()
Return the tree associated with this manager.
Definition: LeafManager.h:298
Vec4< float > Vec4s
Definition: Vec4.h:545
Definition: VolumeToSpheres.h:231
Definition: VolumeToMesh.h:1634
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:119
void run(bool threaded=true)
Definition: VolumeToMesh.h:1708
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:87
int index() const
Definition: VolumeToSpheres.h:596
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:405
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:94
void run(bool threaded=true)
Definition: VolumeToSpheres.h:267
void run(bool threaded=true)
Definition: VolumeToSpheres.h:654
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:134
OPENVDB_STATIC_SPECIALIZATION void erodeVoxels(TreeType &tree, int count=1)
Definition: Morphology.h:384
Definition: Mat4.h:51
Definition: VolumeToSpheres.h:585
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:65