OpenVDB  3.0.0
VolumeToSpheres.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2014 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 tree::LeafManager<const TreeT> LeafManagerT;
728  typedef tree::LeafManager<IntTreeT> IntLeafManagerT;
729  typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
730 
731 
732  typedef boost::mt11213b RandGen;
733  RandGen mtRand(/*seed=*/0);
734 
735  const TreeT& tree = grid.tree();
736  const math::Transform& transform = grid.transform();
737 
738  std::vector<Vec3R> instancePoints;
739 
740  { // Scatter candidate sphere centroids (instancePoints)
741  typename Grid<BoolTreeT>::Ptr interiorMaskPtr;
742 
743  if (grid.getGridClass() == GRID_LEVEL_SET) {
744  interiorMaskPtr = sdfInteriorMask(grid, ValueT(isovalue));
745  } else {
746  interiorMaskPtr = typename Grid<BoolTreeT>::Ptr(Grid<BoolTreeT>::create(false));
747  interiorMaskPtr->setTransform(transform.copy());
748  interiorMaskPtr->tree().topologyUnion(tree);
749  }
750 
751  if (interrupter && interrupter->wasInterrupted()) return;
752 
753  erodeVoxels(interiorMaskPtr->tree(), 1);
754 
755  instancePoints.reserve(instances);
756  internal::PointAccessor ptnAcc(instancePoints);
757 
759  ptnAcc, Index64(addNBPoints ? (instances / 2) : instances), mtRand, interrupter);
760 
761  scatter(*interiorMaskPtr);
762  }
763 
764  if (interrupter && interrupter->wasInterrupted()) return;
765 
766  std::vector<float> instanceRadius;
767 
769  csp.initialize(grid, isovalue, interrupter);
770 
771  // add extra instance points in the interior narrow band.
772  if (instancePoints.size() < size_t(instances)) {
773  const Int16TreeT& signTree = csp.signTree();
774  typename Int16TreeT::LeafNodeType::ValueOnCIter it;
775  typename Int16TreeT::LeafCIter leafIt = signTree.cbeginLeaf();
776 
777  for (; leafIt; ++leafIt) {
778  for (it = leafIt->cbeginValueOn(); it; ++it) {
779  const int flags = it.getValue();
780  if (!(0xE00 & flags) && (flags & 0x100)) {
781  instancePoints.push_back(transform.indexToWorld(it.getCoord()));
782  }
783 
784  if (instancePoints.size() == size_t(instances)) break;
785  }
786  if (instancePoints.size() == size_t(instances)) break;
787  }
788  }
789 
790 
791  if (interrupter && interrupter->wasInterrupted()) return;
792 
793  if (!csp.search(instancePoints, instanceRadius)) return;
794 
795  std::vector<unsigned char> instanceMask(instancePoints.size(), 0);
796  float largestRadius = 0.0;
797  int largestRadiusIdx = 0;
798 
799  for (size_t n = 0, N = instancePoints.size(); n < N; ++n) {
800  if (instanceRadius[n] > largestRadius) {
801  largestRadius = instanceRadius[n];
802  largestRadiusIdx = int(n);
803  }
804  }
805 
806  Vec3s pos;
807  Vec4s sphere;
808  minRadius = float(minRadius * transform.voxelSize()[0]);
809  maxRadius = float(maxRadius * transform.voxelSize()[0]);
810 
811  for (size_t s = 0, S = std::min(size_t(maxSphereCount), instancePoints.size()); s < S; ++s) {
812 
813  if (interrupter && interrupter->wasInterrupted()) return;
814 
815  largestRadius = std::min(maxRadius, largestRadius);
816 
817  if (s != 0 && largestRadius < minRadius) break;
818 
819  sphere[0] = float(instancePoints[largestRadiusIdx].x());
820  sphere[1] = float(instancePoints[largestRadiusIdx].y());
821  sphere[2] = float(instancePoints[largestRadiusIdx].z());
822  sphere[3] = largestRadius;
823 
824  spheres.push_back(sphere);
825  instanceMask[largestRadiusIdx] = 1;
826 
828  sphere, instancePoints, instanceRadius, instanceMask, overlapping);
829  op.run();
830 
831  largestRadius = op.radius();
832  largestRadiusIdx = op.index();
833  }
834 }
835 
837 
838 
839 template<typename GridT>
841  : mIsInitialized(false)
842  , mLeafBoundingSpheres(0)
843  , mNodeBoundingSpheres(0)
844  , mLeafRanges(0)
845  , mLeafNodes(0)
846  , mSurfacePointList()
847  , mPointListSize(0)
848  , mMaxNodeLeafs(0)
849  , mMaxRadiusSqr(0.0)
850  , mIdxTreePt()
851 {
852 }
853 
854 template<typename GridT>
855 void
856 ClosestSurfacePoint<GridT>::initialize(const GridT& grid, float isovalue)
857 {
858  initialize<GridT, util::NullInterrupter>(grid, isovalue, NULL);
859 }
860 
861 
862 template<typename GridT>
863 template<typename InterrupterT>
864 void
866  const GridT& grid, float isovalue, InterrupterT* interrupter)
867 {
868  mIsInitialized = false;
869  typedef tree::LeafManager<const TreeT> LeafManagerT;
870  typedef tree::LeafManager<IntTreeT> IntLeafManagerT;
871  typedef tree::LeafManager<Int16TreeT> Int16LeafManagerT;
872  typedef typename GridT::ValueType ValueT;
873 
874  const TreeT& tree = grid.tree();
875  const math::Transform& transform = grid.transform();
876 
877  { // Extract surface point cloud
878 
879  {
880  LeafManagerT leafs(tree);
882  signDataOp(tree, leafs, ValueT(isovalue));
883 
884  signDataOp.run();
885 
886  mSignTreePt = signDataOp.signTree();
887  mIdxTreePt = signDataOp.idxTree();
888  }
889 
890  if (interrupter && interrupter->wasInterrupted()) return;
891 
892  Int16LeafManagerT signLeafs(*mSignTreePt);
893 
894  std::vector<size_t> regions(signLeafs.leafCount(), 0);
895  signLeafs.foreach(internal::CountPoints(regions));
896 
897  mPointListSize = 0;
898  for (size_t tmp = 0, n = 0, N = regions.size(); n < N; ++n) {
899  tmp = regions[n];
900  regions[n] = mPointListSize;
901  mPointListSize += tmp;
902  }
903 
904  if (mPointListSize == 0) return;
905 
906  mSurfacePointList.reset(new Vec3s[mPointListSize]);
907 
909  pointOp(signLeafs, tree, *mIdxTreePt, mSurfacePointList, regions, transform, isovalue);
910 
911  pointOp.run();
912 
913  mIdxTreePt->topologyUnion(*mSignTreePt);
914  }
915 
916  if (interrupter && interrupter->wasInterrupted()) return;
917 
918  // estimate max sphere radius (sqr dist)
919  CoordBBox bbox = grid.evalActiveVoxelBoundingBox();
920 
921  Vec3s dim = transform.indexToWorld(bbox.min()) -
922  transform.indexToWorld(bbox.max());
923 
924  dim[0] = std::abs(dim[0]);
925  dim[1] = std::abs(dim[1]);
926  dim[2] = std::abs(dim[2]);
927 
928  mMaxRadiusSqr = std::min(std::min(dim[0], dim[1]), dim[2]);
929  mMaxRadiusSqr *= 0.51f;
930  mMaxRadiusSqr *= mMaxRadiusSqr;
931 
932 
933  IntLeafManagerT idxLeafs(*mIdxTreePt);
934 
935 
936  typedef typename IntTreeT::RootNodeType IntRootNodeT;
937  typedef typename IntRootNodeT::NodeChainType IntNodeChainT;
938  BOOST_STATIC_ASSERT(boost::mpl::size<IntNodeChainT>::value > 1);
939  typedef typename boost::mpl::at<IntNodeChainT, boost::mpl::int_<1> >::type IntInternalNodeT;
940 
941 
942  typename IntTreeT::NodeCIter nIt = mIdxTreePt->cbeginNode();
943  nIt.setMinDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
944  nIt.setMaxDepth(IntTreeT::NodeCIter::LEAF_DEPTH - 1);
945 
946  std::vector<const IntInternalNodeT*> internalNodes;
947 
948  const IntInternalNodeT* node = NULL;
949  for (; nIt; ++nIt) {
950  nIt.getNode(node);
951  if (node) internalNodes.push_back(node);
952  }
953 
954  std::vector<IndexRange>().swap(mLeafRanges);
955  mLeafRanges.resize(internalNodes.size());
956 
957  std::vector<const IntLeafT*>().swap(mLeafNodes);
958  mLeafNodes.reserve(idxLeafs.leafCount());
959 
960  typename IntInternalNodeT::ChildOnCIter leafIt;
961  mMaxNodeLeafs = 0;
962  for (size_t n = 0, N = internalNodes.size(); n < N; ++n) {
963 
964  mLeafRanges[n].first = mLeafNodes.size();
965 
966  size_t leafCount = 0;
967  for (leafIt = internalNodes[n]->cbeginChildOn(); leafIt; ++leafIt) {
968  mLeafNodes.push_back(&(*leafIt));
969  ++leafCount;
970  }
971 
972  mMaxNodeLeafs = std::max(leafCount, mMaxNodeLeafs);
973 
974  mLeafRanges[n].second = mLeafNodes.size();
975  }
976 
977  std::vector<Vec4R>().swap(mLeafBoundingSpheres);
978  mLeafBoundingSpheres.resize(mLeafNodes.size());
979 
981  mLeafBoundingSpheres, mLeafNodes, transform, mSurfacePointList);
982  leafBS.run();
983 
984 
985  std::vector<Vec4R>().swap(mNodeBoundingSpheres);
986  mNodeBoundingSpheres.resize(internalNodes.size());
987 
988  internal::NodeBS nodeBS(mNodeBoundingSpheres, mLeafRanges, mLeafBoundingSpheres);
989  nodeBS.run();
990  mIsInitialized = true;
991 }
992 
993 
994 template<typename GridT>
995 bool
996 ClosestSurfacePoint<GridT>::search(std::vector<Vec3R>& points,
997  std::vector<float>& distances, bool transformPoints)
998 {
999  if (!mIsInitialized) return false;
1000 
1001  distances.clear();
1002  distances.resize(points.size(), mMaxRadiusSqr);
1003 
1004  internal::ClosestPointDist<IntLeafT> cpd(points, distances, mSurfacePointList,
1005  mLeafNodes, mLeafRanges, mLeafBoundingSpheres, mNodeBoundingSpheres,
1006  mMaxNodeLeafs, transformPoints);
1007 
1008  cpd.run();
1009 
1010  return true;
1011 }
1012 
1013 
1014 template<typename GridT>
1015 bool
1016 ClosestSurfacePoint<GridT>::search(const std::vector<Vec3R>& points, std::vector<float>& distances)
1017 {
1018  return search(const_cast<std::vector<Vec3R>& >(points), distances, false);
1019 }
1020 
1021 
1022 template<typename GridT>
1023 bool
1025  std::vector<float>& distances)
1026 {
1027  return search(points, distances, true);
1028 }
1029 
1030 
1031 } // namespace tools
1032 } // namespace OPENVDB_VERSION_NAME
1033 } // namespace openvdb
1034 
1035 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
1036 
1037 // Copyright (c) 2012-2014 DreamWorks Animation LLC
1038 // All rights reserved. This software is distributed under the
1039 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
Definition: VolumeToSpheres.h:214
Definition: VolumeToMesh.h:1645
void operator()(const tbb::blocked_range< size_t > &range)
Definition: VolumeToSpheres.h:666
void join(const UpdatePoints &rhs)
Definition: VolumeToSpheres.h:605
const Int16TreeT & signTree() const
Tree accessors.
Definition: VolumeToSpheres.h:183
void run(bool threaded=true)
Definition: VolumeToSpheres.h:344
NodeBS(std::vector< Vec4R > &nodeBoundingSpheres, const std::vector< IndexRange > &leafRanges, const std::vector< Vec4R > &leafBoundingSpheres)
Definition: VolumeToSpheres.h:334
ClosestSurfacePoint()
Definition: VolumeToSpheres.h:840
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
Definition: VolumeToSpheres.h:314
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:358
uint64_t Index64
Definition: Types.h:58
GridT::TreeType TreeT
Definition: VolumeToSpheres.h:132
Definition: VolumeToSpheres.h:587
const IntTreeT & indexTree() const
Tree accessors.
Definition: VolumeToSpheres.h:182
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
TreeT::template ValueConverter< int >::Type IntTreeT
Definition: VolumeToSpheres.h:133
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:404
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
boost::scoped_array< openvdb::Vec3s > PointList
Point and primitive list types.
Definition: VolumeToMesh.h:174
bool search(const std::vector< Vec3R > &points, std::vector< float > &distances)
Computes distance to closest surface.
Definition: VolumeToSpheres.h:1016
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:721
TreeT::template ValueConverter< Int16 >::Type Int16TreeT
Definition: VolumeToSpheres.h:134
IntTreeT::Ptr idxTree() const
Definition: VolumeToMesh.h:925
void run(bool threaded=true)
Definition: VolumeToSpheres.h:656
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:112
Int16TreeT::Ptr signTree() const
Definition: VolumeToMesh.h:924
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
Definition: Types.h:210
void run(bool threaded=true)
Definition: VolumeToSpheres.h:267
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:135
#define OPENVDB_VERSION_NAME
Definition: version.h:43
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:354
std::pair< size_t, size_t > IndexRange
Definition: VolumeToSpheres.h:317
void run(bool threaded=true)
Definition: VolumeToMesh.h:985
Ptr copy() const
Definition: Transform.h:77
Definition: Exceptions.h:39
void run(bool threaded=true)
Definition: VolumeToSpheres.h:474
Vec4< float > Vec4s
Definition: Vec4.h:545
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Definition: VolumeToSpheres.h:231
UpdatePoints(const Vec4s &sphere, const std::vector< Vec3R > &points, std::vector< float > &distances, std::vector< unsigned char > &mask, bool overlapping)
Definition: VolumeToSpheres.h:627
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:97
Vec3< double > Vec3d
Definition: Vec3.h:629
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:548
OPENVDB_IMPORT void initialize()
Global registration of basic types.
Counts the total number of points per leaf, accounts for cells with multiple points.
Definition: VolumeToMesh.h:1062
int index() const
Definition: VolumeToSpheres.h:598
void operator()(const tbb::blocked_range< size_t > &) const
Definition: VolumeToSpheres.h:278
Definition: VolumeToMesh.h:905
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
void run(bool threaded=true)
Definition: VolumeToMesh.h:1719
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
void setTransform(math::Transform::Ptr)
Associate the given transform with this grid, in place of its existing transform. ...
Definition: Grid.h:999
void add(const Vec3R &pos)
Definition: VolumeToSpheres.h:221
bool searchAndReplace(std::vector< Vec3R > &points, std::vector< float > &distances)
Performs closest point searches.
Definition: VolumeToSpheres.h:1024
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:302
PointAccessor(std::vector< Vec3R > &points)
Definition: VolumeToSpheres.h:216
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
Definition: Mat4.h:51
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Vec3< float > Vec3s
Definition: Vec3.h:628
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 initialize(const GridT &grid, float isovalue=0.0, InterrupterT *interrupter=NULL)
Extracts the surface points and constructs a spatial acceleration structure.
Definition: VolumeToSpheres.h:865
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
float radius() const
Definition: VolumeToSpheres.h:597
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:794
Accelerated closest surface point queries for narrow band level sets. Supports queries that originate...
Definition: VolumeToSpheres.h:129