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