OpenVDB  3.0.0
MeshToVolume.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_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/Types.h>
35 #include <openvdb/math/FiniteDifference.h>
36 #include <openvdb/math/Operators.h> // for ISGradientNormSqrd
37 #include <openvdb/math/Proximity.h> // for closestPointOnTriangleToPoint()
38 #include <openvdb/tools/Morphology.h> // for dilateVoxels()
39 #include <openvdb/util/NullInterrupter.h>
40 #include <openvdb/util/Util.h> // for nearestCoord()
41 #include "ChangeBackground.h"
42 #include "Prune.h"// for pruneInactive and pruneLevelSet
43 #include "SignedFloodFill.h" // for signedFloodFillWithValues
44 
45 #include <tbb/blocked_range.h>
46 #include <tbb/parallel_for.h>
47 #include <tbb/parallel_reduce.h>
48 
49 #include <boost/math/special_functions/fpclassify.hpp> // for isfinite()
50 
51 #include <deque>
52 #include <limits>
53 #include <sstream>
54 
55 namespace openvdb {
57 namespace OPENVDB_VERSION_NAME {
58 namespace tools {
59 
60 
62 
63 
64 // Wrapper functions for the MeshToVolume converter
65 
66 
82 template<typename GridType>
83 inline typename GridType::Ptr
85  const openvdb::math::Transform& xform,
86  const std::vector<Vec3s>& points,
87  const std::vector<Vec3I>& triangles,
88  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
89 
90 
106 template<typename GridType>
107 inline typename GridType::Ptr
109  const openvdb::math::Transform& xform,
110  const std::vector<Vec3s>& points,
111  const std::vector<Vec4I>& quads,
112  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
113 
114 
131 template<typename GridType>
132 inline typename GridType::Ptr
134  const openvdb::math::Transform& xform,
135  const std::vector<Vec3s>& points,
136  const std::vector<Vec3I>& triangles,
137  const std::vector<Vec4I>& quads,
138  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
139 
140 
159 template<typename GridType>
160 inline typename GridType::Ptr
162  const openvdb::math::Transform& xform,
163  const std::vector<Vec3s>& points,
164  const std::vector<Vec3I>& triangles,
165  const std::vector<Vec4I>& quads,
166  float exBandWidth,
167  float inBandWidth);
168 
169 
184 template<typename GridType>
185 inline typename GridType::Ptr
187  const openvdb::math::Transform& xform,
188  const std::vector<Vec3s>& points,
189  const std::vector<Vec3I>& triangles,
190  const std::vector<Vec4I>& quads,
191  float bandWidth);
192 
193 
195 
196 
199 
200 
201 // MeshToVolume
202 template<typename FloatGridT, typename InterruptT = util::NullInterrupter>
204 {
205 public:
206  typedef typename FloatGridT::TreeType FloatTreeT;
207  typedef typename FloatTreeT::ValueType FloatValueT;
208  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
210  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
212 
213  MeshToVolume(openvdb::math::Transform::Ptr&, int conversionFlags = 0,
214  InterruptT *interrupter = NULL, int signSweeps = 1);
215 
227  void convertToLevelSet(
228  const std::vector<Vec3s>& pointList,
229  const std::vector<Vec4I>& polygonList,
230  FloatValueT exBandWidth = FloatValueT(LEVEL_SET_HALF_WIDTH),
231  FloatValueT inBandWidth = FloatValueT(LEVEL_SET_HALF_WIDTH));
232 
241  void convertToUnsignedDistanceField(const std::vector<Vec3s>& pointList,
242  const std::vector<Vec4I>& polygonList, FloatValueT exBandWidth);
243 
244  void clear();
245 
247  typename FloatGridT::Ptr distGridPtr() const { return mDistGrid; }
248 
251  typename IntGridT::Ptr indexGridPtr() const { return mIndexGrid; }
252 
253 private:
254  // disallow copy by assignment
255  void operator=(const MeshToVolume<FloatGridT, InterruptT>&) {}
256 
257  void doConvert(const std::vector<Vec3s>&, const std::vector<Vec4I>&,
258  FloatValueT exBandWidth, FloatValueT inBandWidth, bool unsignedDistField = false);
259 
260  bool wasInterrupted(int percent = -1) const {
261  return mInterrupter && mInterrupter->wasInterrupted(percent);
262  }
263 
264  openvdb::math::Transform::Ptr mTransform;
265  int mConversionFlags, mSignSweeps;
266 
267  typename FloatGridT::Ptr mDistGrid;
268  typename IntGridT::Ptr mIndexGrid;
269  typename BoolGridT::Ptr mIntersectingVoxelsGrid;
270 
271  InterruptT *mInterrupter;
272 };
273 
274 
276 
277 
280 {
281 public:
282 
284 
286  struct EdgeData {
287  EdgeData(float dist = 1.0)
288  : mXDist(dist), mYDist(dist), mZDist(dist)
289  , mXPrim(util::INVALID_IDX)
290  , mYPrim(util::INVALID_IDX)
291  , mZPrim(util::INVALID_IDX)
292  {
293  }
294 
296  bool operator< (const EdgeData&) const { return false; }
299  bool operator> (const EdgeData&) const { return false; }
300  template<class T> EdgeData operator+(const T&) const { return *this; }
301  template<class T> EdgeData operator-(const T&) const { return *this; }
302  EdgeData operator-() const { return *this; }
304 
305  bool operator==(const EdgeData& rhs) const
306  {
307  return mXPrim == rhs.mXPrim && mYPrim == rhs.mYPrim && mZPrim == rhs.mZPrim;
308  }
309 
310  float mXDist, mYDist, mZDist;
311  Index32 mXPrim, mYPrim, mZPrim;
312  };
313 
316 
317 
319 
320 
322 
323 
331  void convert(const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
332 
333 
336  void getEdgeData(Accessor& acc, const Coord& ijk,
337  std::vector<Vec3d>& points, std::vector<Index32>& primitives);
338 
341  Accessor getAccessor() { return Accessor(mTree); }
342 
343 private:
344  void operator=(const MeshToVoxelEdgeData&) {}
345  TreeType mTree;
346  class GenEdgeData;
347 };
348 
349 
352 
353 
354 // Internal utility objects and implementation details
355 
356 namespace internal {
357 
358 
360 {
361 public:
362  PointTransform(const std::vector<Vec3s>& pointsIn, std::vector<Vec3s>& pointsOut,
363  const math::Transform& xform)
364  : mPointsIn(pointsIn)
365  , mPointsOut(&pointsOut)
366  , mXform(xform)
367  {
368  }
369 
370  void run(bool threaded = true)
371  {
372  if (threaded) {
373  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPointsOut->size()), *this);
374  } else {
375  (*this)(tbb::blocked_range<size_t>(0, mPointsOut->size()));
376  }
377  }
378 
379  inline void operator()(const tbb::blocked_range<size_t>& range) const
380  {
381  for (size_t n = range.begin(); n < range.end(); ++n) {
382  (*mPointsOut)[n] = Vec3s(mXform.worldToIndex(mPointsIn[n]));
383  }
384  }
385 
386 private:
387  const std::vector<Vec3s>& mPointsIn;
388  std::vector<Vec3s> * const mPointsOut;
389  const math::Transform& mXform;
390 };
391 
392 
393 template<typename ValueType>
394 struct Tolerance
395 {
396  static ValueType epsilon() { return ValueType(1e-7); }
397  static ValueType minNarrowBandWidth() { return ValueType(1.0 + 1e-6); }
398 };
399 
400 
401 template<typename FloatTreeT, typename IntTreeT>
402 inline void
403 combine(FloatTreeT& lhsDist, IntTreeT& lhsIndex, FloatTreeT& rhsDist, IntTreeT& rhsIndex)
404 {
405  typedef typename FloatTreeT::ValueType FloatValueT;
406  typename tree::ValueAccessor<FloatTreeT> lhsDistAccessor(lhsDist);
407  typename tree::ValueAccessor<IntTreeT> lhsIndexAccessor(lhsIndex);
408  typename tree::ValueAccessor<IntTreeT> rhsIndexAccessor(rhsIndex);
409  typename FloatTreeT::LeafCIter iter = rhsDist.cbeginLeaf();
410 
411  FloatValueT rhsValue;
412  Coord ijk;
413 
414  for ( ; iter; ++iter) {
415  typename FloatTreeT::LeafNodeType::ValueOnCIter it = iter->cbeginValueOn();
416 
417  for ( ; it; ++it) {
418 
419  ijk = it.getCoord();
420  rhsValue = it.getValue();
421  FloatValueT& lhsValue = const_cast<FloatValueT&>(lhsDistAccessor.getValue(ijk));
422 
423  if (-rhsValue < std::abs(lhsValue)) {
424  lhsValue = rhsValue;
425  lhsIndexAccessor.setValue(ijk, rhsIndexAccessor.getValue(ijk));
426  }
427  }
428  }
429 }
430 
431 
433 
434 
442 template<typename FloatTreeT, typename InterruptT = util::NullInterrupter>
444 {
445 public:
446  typedef typename FloatTreeT::ValueType FloatValueT;
447  typedef typename FloatTreeT::LeafNodeType FloatLeafT;
449  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
450  typedef typename IntTreeT::LeafNodeType IntLeafT;
452  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
453  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
455 
456 
457  MeshVoxelizer(const std::vector<Vec3s>& pointList,
458  const std::vector<Vec4I>& polygonList, InterruptT *interrupter = NULL);
459 
461 
462  void run(bool threaded = true);
463 
465  void operator() (const tbb::blocked_range<size_t> &range);
467 
468  FloatTreeT& sqrDistTree() { return mSqrDistTree; }
469  IntTreeT& primIndexTree() { return mPrimIndexTree; }
470  BoolTreeT& intersectionTree() { return mIntersectionTree; }
471 
472 private:
473  // disallow copy by assignment
474  void operator=(const MeshVoxelizer<FloatTreeT, InterruptT>&) {}
475  bool wasInterrupted() const { return mInterrupter && mInterrupter->wasInterrupted(); }
476 
477  bool evalVoxel(const Coord& ijk, const Int32 polyIdx);
478 
479  const std::vector<Vec3s>& mPointList;
480  const std::vector<Vec4I>& mPolygonList;
481 
482  FloatTreeT mSqrDistTree;
483  FloatAccessorT mSqrDistAccessor;
484 
485  IntTreeT mPrimIndexTree;
486  IntAccessorT mPrimIndexAccessor;
487 
488  BoolTreeT mIntersectionTree;
489  BoolAccessorT mIntersectionAccessor;
490 
491  // Used internally for acceleration
492  IntTreeT mLastPrimTree;
493  IntAccessorT mLastPrimAccessor;
494 
495  InterruptT *mInterrupter;
496 
497 
498  struct Primitive { Vec3d a, b, c, d; Int32 index; };
499 
500  template<bool IsQuad>
501  bool evalPrimitive(const Coord&, const Primitive&);
502 
503  template<bool IsQuad>
504  void voxelize(const Primitive&);
505 };
506 
507 
508 template<typename FloatTreeT, typename InterruptT>
509 void
511 {
512  if (threaded) {
513  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
514  } else {
515  (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
516  }
517 }
518 
519 template<typename FloatTreeT, typename InterruptT>
521  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
522  InterruptT *interrupter)
523  : mPointList(pointList)
524  , mPolygonList(polygonList)
525  , mSqrDistTree(std::numeric_limits<FloatValueT>::max())
526  , mSqrDistAccessor(mSqrDistTree)
527  , mPrimIndexTree(Int32(util::INVALID_IDX))
528  , mPrimIndexAccessor(mPrimIndexTree)
529  , mIntersectionTree(false)
530  , mIntersectionAccessor(mIntersectionTree)
531  , mLastPrimTree(Int32(util::INVALID_IDX))
532  , mLastPrimAccessor(mLastPrimTree)
533  , mInterrupter(interrupter)
534 {
535 }
536 
537 template<typename FloatTreeT, typename InterruptT>
539  MeshVoxelizer<FloatTreeT, InterruptT>& rhs, tbb::split)
540  : mPointList(rhs.mPointList)
541  , mPolygonList(rhs.mPolygonList)
542  , mSqrDistTree(std::numeric_limits<FloatValueT>::max())
543  , mSqrDistAccessor(mSqrDistTree)
544  , mPrimIndexTree(Int32(util::INVALID_IDX))
545  , mPrimIndexAccessor(mPrimIndexTree)
546  , mIntersectionTree(false)
547  , mIntersectionAccessor(mIntersectionTree)
548  , mLastPrimTree(Int32(util::INVALID_IDX))
549  , mLastPrimAccessor(mLastPrimTree)
550  , mInterrupter(rhs.mInterrupter)
551 {
552 }
553 
554 
555 template<typename FloatTreeT, typename InterruptT>
556 void
557 MeshVoxelizer<FloatTreeT, InterruptT>::operator()(const tbb::blocked_range<size_t> &range)
558 {
559  Primitive prim;
560 
561  for (size_t n = range.begin(); n < range.end(); ++n) {
562 
563  if (mInterrupter && mInterrupter->wasInterrupted()) {
564  tbb::task::self().cancel_group_execution();
565  break;
566  }
567 
568  const Vec4I& verts = mPolygonList[n];
569 
570  prim.index = Int32(n);
571  prim.a = Vec3d(mPointList[verts[0]]);
572  prim.b = Vec3d(mPointList[verts[1]]);
573  prim.c = Vec3d(mPointList[verts[2]]);
574 
575  if (util::INVALID_IDX != verts[3]) {
576  prim.d = Vec3d(mPointList[verts[3]]);
577  voxelize<true>(prim);
578  } else {
579  voxelize<false>(prim);
580  }
581  }
582 }
583 
584 
585 template<typename FloatTreeT, typename InterruptT>
586 template<bool IsQuad>
587 void
589 {
590  std::deque<Coord> coordList;
591  Coord ijk, nijk;
592 
593  ijk = util::nearestCoord(prim.a);
594  coordList.push_back(ijk);
595 
596  evalPrimitive<IsQuad>(ijk, prim);
597 
598  while (!coordList.empty()) {
599  if(wasInterrupted()) break;
600 
601  ijk = coordList.back();
602  coordList.pop_back();
603 
604  mIntersectionAccessor.setActiveState(ijk, true);
605 
606  for (Int32 i = 0; i < 26; ++i) {
607  nijk = ijk + util::COORD_OFFSETS[i];
608  if (prim.index != mLastPrimAccessor.getValue(nijk)) {
609  mLastPrimAccessor.setValue(nijk, prim.index);
610  if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
611  }
612  }
613  }
614 }
615 
616 
617 template<typename FloatTreeT, typename InterruptT>
618 template<bool IsQuad>
619 bool
620 MeshVoxelizer<FloatTreeT, InterruptT>::evalPrimitive(const Coord& ijk, const Primitive& prim)
621 {
622  Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
623 
624  // Evaluate first triangle
625  FloatValueT dist = FloatValueT((voxelCenter -
626  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, voxelCenter, uvw)).lengthSqr());
627 
628  if (IsQuad) {
629  // Split quad into a second triangle and calculate distance.
630  FloatValueT secondDist = FloatValueT((voxelCenter -
631  closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, voxelCenter, uvw)).lengthSqr());
632 
633  if (secondDist < dist) dist = secondDist;
634  }
635 
636  FloatValueT oldDist = std::abs(mSqrDistAccessor.getValue(ijk));
637 
638  if (dist < oldDist) {
639  mSqrDistAccessor.setValue(ijk, -dist);
640  mPrimIndexAccessor.setValue(ijk, prim.index);
641  } else if (math::isExactlyEqual(dist, oldDist)) {
642  // makes reduction deterministic when different polygons
643  // produce the same distance value.
644  mPrimIndexAccessor.setValue(ijk, std::min(prim.index, mPrimIndexAccessor.getValue(ijk)));
645  }
646 
647  return (dist < 0.86602540378443861);
648 }
649 
650 
651 template<typename FloatTreeT, typename InterruptT>
652 void
654 {
655  typedef typename FloatTreeT::RootNodeType FloatRootNodeT;
656  typedef typename FloatRootNodeT::NodeChainType FloatNodeChainT;
657  BOOST_STATIC_ASSERT(boost::mpl::size<FloatNodeChainT>::value > 1);
658  typedef typename boost::mpl::at<FloatNodeChainT, boost::mpl::int_<1> >::type FloatInternalNodeT;
659 
660  typedef typename IntTreeT::RootNodeType IntRootNodeT;
661  typedef typename IntRootNodeT::NodeChainType IntNodeChainT;
662  BOOST_STATIC_ASSERT(boost::mpl::size<IntNodeChainT>::value > 1);
663  typedef typename boost::mpl::at<IntNodeChainT, boost::mpl::int_<1> >::type IntInternalNodeT;
664 
665  const FloatValueT background = std::numeric_limits<FloatValueT>::max();
666 
667  Coord ijk;
668  Index offset;
669 
670  rhs.mSqrDistTree.clearAllAccessors();
671  rhs.mPrimIndexTree.clearAllAccessors();
672 
673  typename FloatTreeT::LeafIter leafIt = rhs.mSqrDistTree.beginLeaf();
674  for ( ; leafIt; ++leafIt) {
675 
676  ijk = leafIt->origin();
677  FloatLeafT* lhsDistLeafPt = mSqrDistAccessor.probeLeaf(ijk);
678 
679  if (!lhsDistLeafPt) {
680 
681  // Steals leaf nodes through their parent, always the last internal-node
682  // stored in the ValueAccessor's node chain, avoiding the overhead of
683  // the root node. This is significantly faster than going through the
684  // tree or root node.
685  mSqrDistAccessor.addLeaf(rhs.mSqrDistAccessor.probeLeaf(ijk));
686  FloatInternalNodeT* floatNode =
687  rhs.mSqrDistAccessor.template getNode<FloatInternalNodeT>();
688  floatNode->template stealNode<FloatLeafT>(ijk, background, false);
689 
690  mPrimIndexAccessor.addLeaf(rhs.mPrimIndexAccessor.probeLeaf(ijk));
691  IntInternalNodeT* intNode =
692  rhs.mPrimIndexAccessor.template getNode<IntInternalNodeT>();
693  intNode->template stealNode<IntLeafT>(ijk, util::INVALID_IDX, false);
694 
695  } else {
696 
697  IntLeafT* lhsIdxLeafPt = mPrimIndexAccessor.probeLeaf(ijk);
698  IntLeafT* rhsIdxLeafPt = rhs.mPrimIndexAccessor.probeLeaf(ijk);
699  FloatValueT lhsValue, rhsValue;
700 
701  typename FloatLeafT::ValueOnCIter it = leafIt->cbeginValueOn();
702  for ( ; it; ++it) {
703 
704  offset = it.pos();
705 
706  lhsValue = std::abs(lhsDistLeafPt->getValue(offset));
707  rhsValue = std::abs(it.getValue());
708 
709  if (rhsValue < lhsValue) {
710  lhsDistLeafPt->setValueOn(offset, it.getValue());
711  lhsIdxLeafPt->setValueOn(offset, rhsIdxLeafPt->getValue(offset));
712  } else if (math::isExactlyEqual(rhsValue, lhsValue)) {
713  lhsIdxLeafPt->setValueOn(offset,
714  std::min(lhsIdxLeafPt->getValue(offset), rhsIdxLeafPt->getValue(offset)));
715  }
716  }
717  }
718  }
719 
720  mIntersectionTree.merge(rhs.mIntersectionTree);
721 
722  rhs.mSqrDistTree.clear();
723  rhs.mPrimIndexTree.clear();
724  rhs.mIntersectionTree.clear();
725 }
726 
727 
729 
730 
731 // ContourTracer
734 template<typename FloatTreeT, typename InterruptT = util::NullInterrupter>
736 {
737 public:
738  typedef typename FloatTreeT::ValueType FloatValueT;
740  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
742 
743  ContourTracer(FloatTreeT&, const BoolTreeT&, InterruptT *interrupter = NULL);
745 
746  void run(bool threaded = true);
747 
749  void operator()(const tbb::blocked_range<int> &range) const;
750 
751 private:
752  void operator=(const ContourTracer<FloatTreeT, InterruptT>&) {}
753 
754  int sparseScan(int slice) const;
755 
756  FloatTreeT& mDistTree;
757  DistAccessorT mDistAccessor;
758 
759  const BoolTreeT& mIntersectionTree;
760  BoolAccessorT mIntersectionAccessor;
761 
762  CoordBBox mBBox;
763 
765  std::vector<Index> mStepSize;
766 
767  InterruptT *mInterrupter;
768 };
769 
770 
771 template<typename FloatTreeT, typename InterruptT>
772 void
774 {
775  if (threaded) {
776  tbb::parallel_for(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1), *this);
777  } else {
778  (*this)(tbb::blocked_range<int>(mBBox.min()[0], mBBox.max()[0]+1));
779  }
780 }
781 
782 
783 template<typename FloatTreeT, typename InterruptT>
785  FloatTreeT& distTree, const BoolTreeT& intersectionTree, InterruptT *interrupter)
786  : mDistTree(distTree)
787  , mDistAccessor(mDistTree)
788  , mIntersectionTree(intersectionTree)
789  , mIntersectionAccessor(mIntersectionTree)
790  , mBBox(CoordBBox())
791  , mStepSize(0)
792  , mInterrupter(interrupter)
793 {
794  // Build the step size table for different tree value depths.
795  std::vector<Index> dims;
796  mDistTree.getNodeLog2Dims(dims);
797 
798  mStepSize.resize(dims.size()+1, 1);
799  Index exponent = 0;
800  for (int idx = static_cast<int>(dims.size()) - 1; idx > -1; --idx) {
801  exponent += dims[idx];
802  mStepSize[idx] = 1 << exponent;
803  }
804 
805  mDistTree.evalLeafBoundingBox(mBBox);
806 
807  // Make sure that mBBox coincides with the min and max corners of the internal nodes.
808  const int tileDim = mStepSize[0];
809 
810  for (size_t i = 0; i < 3; ++i) {
811 
812  int n;
813  double diff = std::abs(double(mBBox.min()[i])) / double(tileDim);
814 
815  if (mBBox.min()[i] <= tileDim) {
816  n = int(std::ceil(diff));
817  mBBox.min()[i] = - n * tileDim;
818  } else {
819  n = int(std::floor(diff));
820  mBBox.min()[i] = n * tileDim;
821  }
822 
823  n = int(std::ceil(std::abs(double(mBBox.max()[i] - mBBox.min()[i])) / double(tileDim)));
824  mBBox.max()[i] = mBBox.min()[i] + n * tileDim;
825  }
826 }
827 
828 
829 template<typename FloatTreeT, typename InterruptT>
832  : mDistTree(rhs.mDistTree)
833  , mDistAccessor(mDistTree)
834  , mIntersectionTree(rhs.mIntersectionTree)
835  , mIntersectionAccessor(mIntersectionTree)
836  , mBBox(rhs.mBBox)
837  , mStepSize(rhs.mStepSize)
838  , mInterrupter(rhs.mInterrupter)
839 {
840 }
841 
842 
843 template<typename FloatTreeT, typename InterruptT>
844 void
845 ContourTracer<FloatTreeT, InterruptT>::operator()(const tbb::blocked_range<int> &range) const
846 {
847  // Slice up the volume and trace contours.
848  int iStep = 1;
849  for (int n = range.begin(); n < range.end(); n += iStep) {
850 
851  if (mInterrupter && mInterrupter->wasInterrupted()) {
852  tbb::task::self().cancel_group_execution();
853  break;
854  }
855 
856  iStep = sparseScan(n);
857  }
858 }
859 
860 
861 template<typename FloatTreeT, typename InterruptT>
862 int
864 {
865  bool lastVoxelWasOut = true;
866  int last_k = mBBox.min()[2];
867 
868  Coord ijk(slice, mBBox.min()[1], mBBox.min()[2]);
869  Coord step(mStepSize[mDistAccessor.getValueDepth(ijk) + 1]);
870  Coord n_ijk;
871 
872  for (ijk[1] = mBBox.min()[1]; ijk[1] <= mBBox.max()[1]; ijk[1] += step[1]) { // j
873 
874  if (mInterrupter && mInterrupter->wasInterrupted()) {
875  break;
876  }
877 
878  step[1] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
879  step[0] = std::min(step[0], step[1]);
880 
881  for (ijk[2] = mBBox.min()[2]; ijk[2] <= mBBox.max()[2]; ijk[2] += step[2]) { // k
882 
883  step[2] = mStepSize[mDistAccessor.getValueDepth(ijk) + 1];
884  step[1] = std::min(step[1], step[2]);
885  step[0] = std::min(step[0], step[2]);
886 
887  // If the current voxel is set?
888  if (mDistAccessor.isValueOn(ijk)) {
889 
890  // Is this a boundary voxel?
891  if (mIntersectionAccessor.isValueOn(ijk)) {
892 
893  lastVoxelWasOut = false;
894  last_k = ijk[2];
895 
896  } else if (lastVoxelWasOut) {
897 
898  FloatValueT& val = const_cast<FloatValueT&>(mDistAccessor.getValue(ijk));
899  val = -val; // flip sign
900 
901  } else {
902 
903  FloatValueT val;
904  for (Int32 n = 3; n < 6; n += 2) {
905  n_ijk = ijk + util::COORD_OFFSETS[n];
906 
907  if (mDistAccessor.probeValue(n_ijk, val) && val > 0) {
908  lastVoxelWasOut = true;
909  break;
910  }
911  }
912 
913  if (lastVoxelWasOut) {
914 
915  FloatValueT& v = const_cast<FloatValueT&>(mDistAccessor.getValue(ijk));
916  v = -v; // flip sign
917 
918  const int tmp_k = ijk[2];
919 
920  // backtrace
921  for (--ijk[2]; ijk[2] >= last_k; --ijk[2]) {
922  if (mIntersectionAccessor.isValueOn(ijk)) break;
923  FloatValueT& vb =
924  const_cast<FloatValueT&>(mDistAccessor.getValue(ijk));
925  if (vb < FloatValueT(0.0)) vb = -vb; // flip sign
926  }
927 
928  last_k = tmp_k;
929  ijk[2] = tmp_k;
930 
931  } else {
932  last_k = std::min(ijk[2], last_k);
933  }
934 
935  }
936 
937  } // end isValueOn check
938  } // end k
939  } // end j
940  return step[0];
941 }
942 
943 
945 
946 
948 template<typename FloatTreeT, typename InterruptT = util::NullInterrupter>
949 class SignMask
950 {
951 public:
952  typedef typename FloatTreeT::ValueType FloatValueT;
953  typedef typename FloatTreeT::LeafNodeType FloatLeafT;
956  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
957  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
960 
961 
962  SignMask(const FloatLeafManager&, const FloatTreeT&, const BoolTreeT&,
963  InterruptT *interrupter = NULL);
964 
966 
967  void run(bool threaded = true);
968 
969  SignMask(SignMask<FloatTreeT, InterruptT>& rhs, tbb::split);
970  void operator() (const tbb::blocked_range<size_t> &range);
972 
973  BoolTreeT& signMaskTree() { return mSignMaskTree; }
974 
975 private:
976  // disallow copy by assignment
977  void operator=(const SignMask<FloatTreeT, InterruptT>&) {}
978  bool wasInterrupted() const { return mInterrupter && mInterrupter->wasInterrupted(); }
979 
980  const FloatLeafManager& mDistLeafs;
981  const FloatTreeT& mDistTree;
982  const BoolTreeT& mIntersectionTree;
983 
984  BoolTreeT mSignMaskTree;
985 
986  InterruptT *mInterrupter;
987 }; // class SignMask
988 
989 
990 template<typename FloatTreeT, typename InterruptT>
992  const FloatLeafManager& distLeafs, const FloatTreeT& distTree,
993  const BoolTreeT& intersectionTree, InterruptT *interrupter)
994  : mDistLeafs(distLeafs)
995  , mDistTree(distTree)
996  , mIntersectionTree(intersectionTree)
997  , mSignMaskTree(false)
998  , mInterrupter(interrupter)
999 {
1000 }
1001 
1002 
1003 template<typename FloatTreeT, typename InterruptT>
1005  SignMask<FloatTreeT, InterruptT>& rhs, tbb::split)
1006  : mDistLeafs(rhs.mDistLeafs)
1007  , mDistTree(rhs.mDistTree)
1008  , mIntersectionTree(rhs.mIntersectionTree)
1009  , mSignMaskTree(false)
1010  , mInterrupter(rhs.mInterrupter)
1011 {
1012 }
1013 
1014 
1015 template<typename FloatTreeT, typename InterruptT>
1016 void
1018 {
1019  if (threaded) tbb::parallel_reduce(mDistLeafs.getRange(), *this);
1020  else (*this)(mDistLeafs.getRange());
1021 }
1022 
1023 
1024 template<typename FloatTreeT, typename InterruptT>
1025 void
1026 SignMask<FloatTreeT, InterruptT>::operator()(const tbb::blocked_range<size_t> &range)
1027 {
1028  FloatAccessorT distAcc(mDistTree);
1029  BoolConstAccessorT intersectionAcc(mIntersectionTree);
1030  BoolAccessorT maskAcc(mSignMaskTree);
1031 
1032  FloatValueT value;
1033  CoordBBox bbox;
1034  Coord& maxCoord = bbox.max();
1035  Coord& minCoord = bbox.min();
1036  Coord ijk;
1037  const int extent = BoolLeafT::DIM - 1;
1038 
1039  for (size_t n = range.begin(); n < range.end(); ++n) {
1040 
1041  const FloatLeafT& distLeaf = mDistLeafs.leaf(n);
1042 
1043  minCoord = distLeaf.origin();
1044  maxCoord[0] = minCoord[0] + extent;
1045  maxCoord[1] = minCoord[1] + extent;
1046  maxCoord[2] = minCoord[2] + extent;
1047 
1048  const BoolLeafT* intersectionLeaf = intersectionAcc.probeConstLeaf(minCoord);
1049 
1050  BoolLeafT* maskLeafPt = new BoolLeafT(minCoord, false);
1051  BoolLeafT& maskLeaf = *maskLeafPt;
1052  bool addLeaf = false;
1053 
1054  bbox.expand(-1);
1055 
1056  typename FloatLeafT::ValueOnCIter it = distLeaf.cbeginValueOn();
1057  for (; it; ++it) {
1058  if (intersectionLeaf && intersectionLeaf->isValueOn(it.pos())) continue;
1059  if (it.getValue() < FloatValueT(0.0)) {
1060  ijk = it.getCoord();
1061  if (bbox.isInside(ijk)) {
1062  for (size_t i = 0; i < 6; ++i) {
1063  if (distLeaf.probeValue(ijk+util::COORD_OFFSETS[i], value) && value>0.0) {
1064  maskLeaf.setValueOn(ijk);
1065  addLeaf = true;
1066  break;
1067  }
1068  }
1069  } else {
1070  for (size_t i = 0; i < 6; ++i) {
1071  if (distAcc.probeValue(ijk+util::COORD_OFFSETS[i], value) && value>0.0) {
1072  maskLeaf.setValueOn(ijk);
1073  addLeaf = true;
1074  break;
1075  }
1076  }
1077  }
1078  }
1079  }
1080 
1081  if (addLeaf) maskAcc.addLeaf(maskLeafPt);
1082  else delete maskLeafPt;
1083  }
1084 }
1085 
1086 
1087 template<typename FloatTreeT, typename InterruptT>
1088 void
1090 {
1091  mSignMaskTree.merge(rhs.mSignMaskTree);
1092 }
1093 
1094 
1096 
1097 
1099 template<typename FloatTreeT, typename InterruptT = util::NullInterrupter>
1101 {
1102 public:
1103  typedef typename FloatTreeT::ValueType FloatValueT;
1104  typedef typename FloatTreeT::LeafNodeType FloatLeafT;
1106  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1107  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1111 
1112  PropagateSign(BoolLeafManager&, FloatTreeT&, const BoolTreeT&, InterruptT *interrupter = NULL);
1113 
1115 
1116  void run(bool threaded = true);
1117 
1119  void operator() (const tbb::blocked_range<size_t> &range);
1121 
1122  BoolTreeT& signMaskTree() { return mSignMaskTree; }
1123 
1124 private:
1125  // disallow copy by assignment
1126  void operator=(const PropagateSign<FloatTreeT, InterruptT>&);
1127  bool wasInterrupted() const { return mInterrupter && mInterrupter->wasInterrupted(); }
1128 
1129  BoolLeafManager& mOldSignMaskLeafs;
1130  FloatTreeT& mDistTree;
1131  const BoolTreeT& mIntersectionTree;
1132 
1133  BoolTreeT mSignMaskTree;
1134  InterruptT *mInterrupter;
1135 };
1136 
1137 
1138 template<typename FloatTreeT, typename InterruptT>
1140  FloatTreeT& distTree, const BoolTreeT& intersectionTree, InterruptT *interrupter)
1141  : mOldSignMaskLeafs(signMaskLeafs)
1142  , mDistTree(distTree)
1143  , mIntersectionTree(intersectionTree)
1144  , mSignMaskTree(false)
1145  , mInterrupter(interrupter)
1146 {
1147 }
1148 
1149 
1150 template<typename FloatTreeT, typename InterruptT>
1152  PropagateSign<FloatTreeT, InterruptT>& rhs, tbb::split)
1153  : mOldSignMaskLeafs(rhs.mOldSignMaskLeafs)
1154  , mDistTree(rhs.mDistTree)
1155  , mIntersectionTree(rhs.mIntersectionTree)
1156  , mSignMaskTree(false)
1157  , mInterrupter(rhs.mInterrupter)
1158 {
1159 }
1160 
1161 
1162 template<typename FloatTreeT, typename InterruptT>
1163 void
1165 {
1166  if (threaded) tbb::parallel_reduce(mOldSignMaskLeafs.getRange(), *this);
1167  else (*this)(mOldSignMaskLeafs.getRange());
1168 }
1169 
1170 
1171 template<typename FloatTreeT, typename InterruptT>
1172 void
1173 PropagateSign<FloatTreeT, InterruptT>::operator()(const tbb::blocked_range<size_t> &range)
1174 {
1175  FloatAccessorT distAcc(mDistTree);
1176  BoolConstAccessorT intersectionAcc(mIntersectionTree);
1177  BoolAccessorT maskAcc(mSignMaskTree);
1178 
1179  std::deque<Coord> coordList;
1180 
1181  FloatValueT value;
1182  CoordBBox bbox;
1183  Coord& maxCoord = bbox.max();
1184  Coord& minCoord = bbox.min();
1185  Coord ijk, nijk;
1186  const int extent = BoolLeafT::DIM - 1;
1187 
1188  for (size_t n = range.begin(); n < range.end(); ++n) {
1189  BoolLeafT& oldMaskLeaf = mOldSignMaskLeafs.leaf(n);
1190 
1191  minCoord = oldMaskLeaf.origin();
1192  maxCoord[0] = minCoord[0] + extent;
1193  maxCoord[1] = minCoord[1] + extent;
1194  maxCoord[2] = minCoord[2] + extent;
1195 
1196  FloatLeafT& distLeaf = *distAcc.probeLeaf(minCoord);
1197  const BoolLeafT* intersectionLeaf = intersectionAcc.probeConstLeaf(minCoord);
1198 
1199  typename BoolLeafT::ValueOnCIter it = oldMaskLeaf.cbeginValueOn();
1200  for (; it; ++it) {
1201  coordList.push_back(it.getCoord());
1202 
1203  while (!coordList.empty()) {
1204 
1205  ijk = coordList.back();
1206  coordList.pop_back();
1207 
1208  FloatValueT& dist = const_cast<FloatValueT&>(distLeaf.getValue(ijk));
1209  if (dist < FloatValueT(0.0)) {
1210  dist = -dist; // flip sign
1211 
1212  for (size_t i = 0; i < 6; ++i) {
1213  nijk = ijk + util::COORD_OFFSETS[i];
1214  if (bbox.isInside(nijk)) {
1215  if (intersectionLeaf && intersectionLeaf->isValueOn(nijk)) continue;
1216 
1217  if (distLeaf.probeValue(nijk, value) && value < 0.0) {
1218  coordList.push_back(nijk);
1219  }
1220 
1221  } else {
1222  if(!intersectionAcc.isValueOn(nijk) &&
1223  distAcc.probeValue(nijk, value) && value < 0.0) {
1224  maskAcc.setValueOn(nijk);
1225  }
1226  }
1227  }
1228  }
1229  }
1230  }
1231  }
1232 }
1233 
1234 
1235 template<typename FloatTreeT, typename InterruptT>
1236 void
1238 {
1239  mSignMaskTree.merge(rhs.mSignMaskTree);
1240 }
1241 
1242 
1244 
1245 
1246 // IntersectingVoxelSign
1250 template<typename FloatTreeT>
1252 {
1253 public:
1254  typedef typename FloatTreeT::ValueType FloatValueT;
1256  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
1258  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1261 
1263  const std::vector<Vec3s>& pointList,
1264  const std::vector<Vec4I>& polygonList,
1265  FloatTreeT& distTree,
1266  IntTreeT& indexTree,
1267  BoolTreeT& intersectionTree,
1268  BoolLeafManager& leafs);
1269 
1271 
1272  void run(bool threaded = true);
1273 
1275  void operator()(const tbb::blocked_range<size_t>&) const;
1276 
1277 private:
1278  void operator=(const IntersectingVoxelSign<FloatTreeT>&) {}
1279 
1280  Vec3d getClosestPoint(const Coord& ijk, const Vec4I& prim) const;
1281 
1282  std::vector<Vec3s> const * const mPointList;
1283  std::vector<Vec4I> const * const mPolygonList;
1284 
1285  FloatTreeT& mDistTree;
1286  IntTreeT& mIndexTree;
1287  BoolTreeT& mIntersectionTree;
1288 
1289  BoolLeafManager& mLeafs;
1290 };
1291 
1292 
1293 template<typename FloatTreeT>
1294 void
1296 {
1297  if (threaded) tbb::parallel_for(mLeafs.getRange(), *this);
1298  else (*this)(mLeafs.getRange());
1299 }
1300 
1301 
1302 template<typename FloatTreeT>
1304  const std::vector<Vec3s>& pointList,
1305  const std::vector<Vec4I>& polygonList,
1306  FloatTreeT& distTree,
1307  IntTreeT& indexTree,
1308  BoolTreeT& intersectionTree,
1309  BoolLeafManager& leafs)
1310  : mPointList(&pointList)
1311  , mPolygonList(&polygonList)
1312  , mDistTree(distTree)
1313  , mIndexTree(indexTree)
1314  , mIntersectionTree(intersectionTree)
1315  , mLeafs(leafs)
1316 {
1317 }
1318 
1319 
1320 template<typename FloatTreeT>
1323  : mPointList(rhs.mPointList)
1324  , mPolygonList(rhs.mPolygonList)
1325  , mDistTree(rhs.mDistTree)
1326  , mIndexTree(rhs.mIndexTree)
1327  , mIntersectionTree(rhs.mIntersectionTree)
1328  , mLeafs(rhs.mLeafs)
1329 {
1330 }
1331 
1332 
1333 template<typename FloatTreeT>
1334 void
1336  const tbb::blocked_range<size_t>& range) const
1337 {
1338  Coord ijk, nijk;
1339 
1340  FloatAccessorT distAcc(mDistTree);
1341  BoolAccessorT maskAcc(mIntersectionTree);
1342  IntAccessorT idxAcc(mIndexTree);
1343 
1344  FloatValueT tmpValue;
1345  Vec3d cpt, center, dir1, dir2;
1346 
1347  typename BoolTreeT::LeafNodeType::ValueOnCIter iter;
1348  for (size_t n = range.begin(); n < range.end(); ++n) {
1349  iter = mLeafs.leaf(n).cbeginValueOn();
1350  for (; iter; ++iter) {
1351 
1352  ijk = iter.getCoord();
1353 
1354  FloatValueT value = distAcc.getValue(ijk);
1355 
1356  if (!(value < FloatValueT(0.0))) continue;
1357 
1358  center = Vec3d(ijk[0], ijk[1], ijk[2]);
1359 
1360  for (Int32 i = 0; i < 26; ++i) {
1361  nijk = ijk + util::COORD_OFFSETS[i];
1362 
1363  if (!maskAcc.isValueOn(nijk) && distAcc.probeValue(nijk, tmpValue)) {
1364  if (tmpValue < FloatValueT(0.0)) continue;
1365 
1366  const Vec4I& prim = (*mPolygonList)[idxAcc.getValue(nijk)];
1367 
1368  cpt = getClosestPoint(nijk, prim);
1369 
1370  dir1 = center - cpt;
1371  dir1.normalize();
1372 
1373  dir2 = Vec3d(nijk[0], nijk[1], nijk[2]) - cpt;
1374  dir2.normalize();
1375 
1376  if (dir2.dot(dir1) > 0.0) {
1377  distAcc.setValue(ijk, -value);
1378  break;
1379  }
1380  }
1381  }
1382  }
1383  }
1384 }
1385 
1386 
1387 template<typename FloatTreeT>
1388 Vec3d
1389 IntersectingVoxelSign<FloatTreeT>::getClosestPoint(const Coord& ijk, const Vec4I& prim) const
1390 {
1391  Vec3d voxelCenter(ijk[0], ijk[1], ijk[2]);
1392 
1393  // Evaluate first triangle
1394  const Vec3d a((*mPointList)[prim[0]]);
1395  const Vec3d b((*mPointList)[prim[1]]);
1396  const Vec3d c((*mPointList)[prim[2]]);
1397 
1398  Vec3d uvw;
1399  Vec3d cpt1 = closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw);
1400 
1401  // Evaluate second triangle if quad.
1402  if (prim[3] != util::INVALID_IDX) {
1403 
1404  Vec3d diff1 = voxelCenter - cpt1;
1405 
1406  const Vec3d d((*mPointList)[prim[3]]);
1407 
1408  Vec3d cpt2 = closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw);
1409  Vec3d diff2 = voxelCenter - cpt2;
1410 
1411  if (diff2.lengthSqr() < diff1.lengthSqr()) {
1412  return cpt2;
1413  }
1414  }
1415 
1416  return cpt1;
1417 }
1418 
1419 
1421 
1422 
1423 // IntersectingVoxelCleaner
1426 template<typename FloatTreeT>
1428 {
1429 public:
1430  typedef typename FloatTreeT::ValueType FloatValueT;
1432  typedef typename FloatTreeT::LeafNodeType DistLeafT;
1433  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
1435  typedef typename IntTreeT::LeafNodeType IntLeafT;
1436  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1438  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1440 
1441  IntersectingVoxelCleaner(FloatTreeT& distTree, IntTreeT& indexTree,
1442  BoolTreeT& intersectionTree, BoolLeafManager& leafs);
1443 
1445 
1446  void run(bool threaded = true);
1447 
1449  void operator()(const tbb::blocked_range<size_t>&) const;
1450 
1451 private:
1452  void operator=(const IntersectingVoxelCleaner<FloatTreeT>&) {}
1453 
1454  FloatTreeT& mDistTree;
1455  IntTreeT& mIndexTree;
1456  BoolTreeT& mIntersectionTree;
1457  BoolLeafManager& mLeafs;
1458 };
1459 
1460 
1461 template<typename FloatTreeT>
1462 void
1464 {
1465  if (threaded) tbb::parallel_for(mLeafs.getRange(), *this);
1466  else (*this)(mLeafs.getRange());
1467 
1468  tools::pruneInactive(mIntersectionTree, threaded);
1469 }
1470 
1471 
1472 template<typename FloatTreeT>
1474  FloatTreeT& distTree,
1475  IntTreeT& indexTree,
1476  BoolTreeT& intersectionTree,
1477  BoolLeafManager& leafs)
1478  : mDistTree(distTree)
1479  , mIndexTree(indexTree)
1480  , mIntersectionTree(intersectionTree)
1481  , mLeafs(leafs)
1482 {
1483 }
1484 
1485 
1486 template<typename FloatTreeT>
1489  : mDistTree(rhs.mDistTree)
1490  , mIndexTree(rhs.mIndexTree)
1491  , mIntersectionTree(rhs.mIntersectionTree)
1492  , mLeafs(rhs.mLeafs)
1493 {
1494 }
1495 
1496 
1497 template<typename FloatTreeT>
1498 void
1500  const tbb::blocked_range<size_t>& range) const
1501 {
1502  Coord ijk, m_ijk;
1503  bool turnOff;
1504  FloatValueT value;
1505  Index offset;
1506 
1507  typename BoolLeafT::ValueOnCIter iter;
1508 
1509  IntAccessorT indexAcc(mIndexTree);
1510  DistAccessorT distAcc(mDistTree);
1511  BoolAccessorT maskAcc(mIntersectionTree);
1512 
1513  for (size_t n = range.begin(); n < range.end(); ++n) {
1514 
1515  BoolLeafT& maskLeaf = mLeafs.leaf(n);
1516 
1517  ijk = maskLeaf.origin();
1518 
1519  DistLeafT * distLeaf = distAcc.probeLeaf(ijk);
1520  if (distLeaf) {
1521  iter = maskLeaf.cbeginValueOn();
1522  for (; iter; ++iter) {
1523 
1524  offset = iter.pos();
1525 
1526  if(distLeaf->getValue(offset) > 0.0) continue;
1527 
1528  ijk = iter.getCoord();
1529  turnOff = true;
1530  for (Int32 m = 0; m < 26; ++m) {
1531  m_ijk = ijk + util::COORD_OFFSETS[m];
1532  if (distAcc.probeValue(m_ijk, value)) {
1533  if (value > 0.0) {
1534  turnOff = false;
1535  break;
1536  }
1537  }
1538  }
1539 
1540  if (turnOff) {
1541  maskLeaf.setValueOff(offset);
1542  distLeaf->setValueOn(offset, FloatValueT(-0.86602540378443861));
1543  }
1544  }
1545  }
1546  }
1547 }
1548 
1549 
1551 
1552 
1553 // ShellVoxelCleaner
1556 template<typename FloatTreeT>
1558 {
1559 public:
1560  typedef typename FloatTreeT::ValueType FloatValueT;
1562  typedef typename FloatTreeT::LeafNodeType DistLeafT;
1564  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
1566  typedef typename IntTreeT::LeafNodeType IntLeafT;
1567  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1569  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1570 
1571  ShellVoxelCleaner(FloatTreeT& distTree, DistArrayT& leafs, IntTreeT& indexTree,
1572  BoolTreeT& intersectionTree);
1573 
1575 
1576  void run(bool threaded = true);
1577 
1579  void operator()(const tbb::blocked_range<size_t>&) const;
1580 
1581 private:
1582  void operator=(const ShellVoxelCleaner<FloatTreeT>&) {}
1583 
1584  FloatTreeT& mDistTree;
1585  DistArrayT& mLeafs;
1586  IntTreeT& mIndexTree;
1587  BoolTreeT& mIntersectionTree;
1588 };
1589 
1590 
1591 template<typename FloatTreeT>
1592 void
1594 {
1595  if (threaded) tbb::parallel_for(mLeafs.getRange(), *this);
1596  else (*this)(mLeafs.getRange());
1597 
1598  tools::pruneInactive(mDistTree, threaded);
1599  tools::pruneInactive(mIndexTree, threaded);
1600 }
1601 
1602 
1603 template<typename FloatTreeT>
1605  FloatTreeT& distTree,
1606  DistArrayT& leafs,
1607  IntTreeT& indexTree,
1608  BoolTreeT& intersectionTree)
1609  : mDistTree(distTree)
1610  , mLeafs(leafs)
1611  , mIndexTree(indexTree)
1612  , mIntersectionTree(intersectionTree)
1613 {
1614 }
1615 
1616 
1617 template<typename FloatTreeT>
1619  const ShellVoxelCleaner<FloatTreeT> &rhs)
1620  : mDistTree(rhs.mDistTree)
1621  , mLeafs(rhs.mLeafs)
1622  , mIndexTree(rhs.mIndexTree)
1623  , mIntersectionTree(rhs.mIntersectionTree)
1624 {
1625 }
1626 
1627 
1628 template<typename FloatTreeT>
1629 void
1631  const tbb::blocked_range<size_t>& range) const
1632 {
1633  Coord ijk, m_ijk;
1634  bool turnOff;
1635  FloatValueT value;
1636  Index offset;
1637 
1638  typename DistLeafT::ValueOnCIter iter;
1639  const FloatValueT distBG = mDistTree.background();
1640  const Int32 indexBG = mIntersectionTree.background();
1641 
1642  IntAccessorT indexAcc(mIndexTree);
1643  DistAccessorT distAcc(mDistTree);
1644  BoolAccessorT maskAcc(mIntersectionTree);
1645 
1646 
1647  for (size_t n = range.begin(); n < range.end(); ++n) {
1648 
1649  DistLeafT& distLeaf = mLeafs.leaf(n);
1650 
1651  ijk = distLeaf.origin();
1652 
1653  const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
1654  IntLeafT& indexLeaf = *indexAcc.probeLeaf(ijk);
1655 
1656  iter = distLeaf.cbeginValueOn();
1657  for (; iter; ++iter) {
1658 
1659  value = iter.getValue();
1660  if(value > 0.0) continue;
1661 
1662  offset = iter.pos();
1663  if (maskLeaf && maskLeaf->isValueOn(offset)) continue;
1664 
1665  ijk = iter.getCoord();
1666  turnOff = true;
1667  for (Int32 m = 0; m < 26; ++m) {
1668  m_ijk = ijk + util::COORD_OFFSETS[m];
1669  if (maskAcc.isValueOn(m_ijk)) {
1670  turnOff = false;
1671  break;
1672  }
1673  }
1674 
1675  if (turnOff) {
1676  distLeaf.setValueOff(offset, distBG);
1677  indexLeaf.setValueOff(offset, indexBG);
1678  }
1679  }
1680  }
1681 }
1682 
1683 
1685 
1686 
1687 template<typename TreeType>
1689 {
1691 
1692  CopyActiveVoxelsOp(TreeType& tree) : mAcc(tree) { }
1693 
1694  template <typename LeafNodeType>
1695  void operator()(LeafNodeType &leaf, size_t) const
1696  {
1697  LeafNodeType* rhsLeaf = const_cast<LeafNodeType*>(mAcc.probeLeaf(leaf.origin()));
1698  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1699  for (; iter; ++iter) {
1700  rhsLeaf->setValueOnly(iter.pos(), iter.getValue());
1701  }
1702  }
1703 
1704 private:
1705  AccessorT mAcc;
1706 };
1707 
1708 
1709 // ExpandNB
1712 template<typename FloatTreeT>
1714 {
1715 public:
1716  typedef typename FloatTreeT::ValueType FloatValueT;
1717  typedef typename FloatTreeT::LeafNodeType FloatLeafT;
1719  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
1720  typedef typename IntTreeT::LeafNodeType IntLeafT;
1722  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1723  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1726 
1727  ExpandNB(BoolLeafManager& leafs,
1728  FloatTreeT& distTree, IntTreeT& indexTree, BoolTreeT& maskTree,
1729  FloatValueT exteriorBandWidth, FloatValueT interiorBandWidth, FloatValueT voxelSize,
1730  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
1731 
1732  void run(bool threaded = true);
1733 
1734  void operator()(const tbb::blocked_range<size_t>&);
1735  void join(ExpandNB<FloatTreeT>&);
1736  ExpandNB(const ExpandNB<FloatTreeT>&, tbb::split);
1738 
1739 private:
1740  void operator=(const ExpandNB<FloatTreeT>&) {}
1741 
1742  double evalVoxelDist(const Coord&, FloatAccessorT&, IntAccessorT&,
1743  BoolAccessorT&, std::vector<Int32>&, Int32&) const;
1744 
1745  double evalVoxelDist(const Coord&, FloatLeafT&, IntLeafT&,
1746  BoolLeafT&, std::vector<Int32>&, Int32&) const;
1747 
1748  double closestPrimDist(const Coord&, std::vector<Int32>&, Int32&) const;
1749 
1750  BoolLeafManager& mMaskLeafs;
1751 
1752  FloatTreeT& mDistTree;
1753  IntTreeT& mIndexTree;
1754  BoolTreeT& mMaskTree;
1755 
1756  const FloatValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1757  const std::vector<Vec3s>& mPointList;
1758  const std::vector<Vec4I>& mPolygonList;
1759 
1760  FloatTreeT mNewDistTree;
1761  IntTreeT mNewIndexTree;
1762  BoolTreeT mNewMaskTree;
1763 };
1764 
1765 
1766 template<typename FloatTreeT>
1768  BoolLeafManager& leafs,
1769  FloatTreeT& distTree,
1770  IntTreeT& indexTree,
1771  BoolTreeT& maskTree,
1772  FloatValueT exteriorBandWidth,
1773  FloatValueT interiorBandWidth,
1774  FloatValueT voxelSize,
1775  const std::vector<Vec3s>& pointList,
1776  const std::vector<Vec4I>& polygonList)
1777  : mMaskLeafs(leafs)
1778  , mDistTree(distTree)
1779  , mIndexTree(indexTree)
1780  , mMaskTree(maskTree)
1781  , mExteriorBandWidth(exteriorBandWidth)
1782  , mInteriorBandWidth(interiorBandWidth)
1783  , mVoxelSize(voxelSize)
1784  , mPointList(pointList)
1785  , mPolygonList(polygonList)
1786  , mNewDistTree(std::numeric_limits<FloatValueT>::max())
1787  , mNewIndexTree(Int32(util::INVALID_IDX))
1788  , mNewMaskTree(false)
1789 {
1790 }
1791 
1792 
1793 template<typename FloatTreeT>
1795  : mMaskLeafs(rhs.mMaskLeafs)
1796  , mDistTree(rhs.mDistTree)
1797  , mIndexTree(rhs.mIndexTree)
1798  , mMaskTree(rhs.mMaskTree)
1799  , mExteriorBandWidth(rhs.mExteriorBandWidth)
1800  , mInteriorBandWidth(rhs.mInteriorBandWidth)
1801  , mVoxelSize(rhs.mVoxelSize)
1802  , mPointList(rhs.mPointList)
1803  , mPolygonList(rhs.mPolygonList)
1804  , mNewDistTree(std::numeric_limits<FloatValueT>::max())
1805  , mNewIndexTree(Int32(util::INVALID_IDX))
1806  , mNewMaskTree(false)
1807 {
1808 }
1809 
1810 
1811 template<typename FloatTreeT>
1812 void
1814 {
1815  if (threaded) tbb::parallel_reduce(mMaskLeafs.getRange(), *this);
1816  else (*this)(mMaskLeafs.getRange());
1817 
1818  // Copy only the active voxels (tree::merge does branch stealing
1819  // which also moves indicative values).
1820  mDistTree.topologyUnion(mNewDistTree);
1821  tree::LeafManager<FloatTreeT> leafs(mNewDistTree);
1822  leafs.foreach(CopyActiveVoxelsOp<FloatTreeT>(mDistTree));
1823 
1824  mIndexTree.merge(mNewIndexTree);
1825 
1826  mMaskTree.clear();
1827  mMaskTree.merge(mNewMaskTree);
1828 
1829 }
1830 
1831 
1832 template<typename FloatTreeT>
1833 void
1834 ExpandNB<FloatTreeT>::operator()(const tbb::blocked_range<size_t>& range)
1835 {
1836  Coord ijk;
1837  Int32 closestPrim = 0;
1838  Index pos = 0;
1839  FloatValueT distance;
1840  bool inside;
1841 
1842  FloatAccessorT newDistAcc(mNewDistTree);
1843  IntAccessorT newIndexAcc(mNewIndexTree);
1844  BoolAccessorT newMaskAcc(mNewMaskTree);
1845 
1846  FloatAccessorT distAcc(mDistTree);
1847  IntAccessorT indexAcc(mIndexTree);
1848  BoolAccessorT maskAcc(mMaskTree);
1849 
1850  CoordBBox bbox;
1851  std::vector<Int32> primitives(18);
1852 
1853  for (size_t n = range.begin(); n < range.end(); ++n) {
1854 
1855  BoolLeafT& maskLeaf = mMaskLeafs.leaf(n);
1856 
1857  if (maskLeaf.isEmpty()) continue;
1858 
1859  ijk = maskLeaf.origin();
1860 
1861  FloatLeafT* distLeafPt = distAcc.probeLeaf(ijk);
1862 
1863  if (!distLeafPt) {
1864  distLeafPt = new FloatLeafT(ijk, distAcc.getValue(ijk));
1865  newDistAcc.addLeaf(distLeafPt);
1866  }
1867 
1868  IntLeafT* indexLeafPt = indexAcc.probeLeaf(ijk);
1869  if (!indexLeafPt) indexLeafPt = newIndexAcc.touchLeaf(ijk);
1870 
1871  bbox = maskLeaf.getNodeBoundingBox();
1872  bbox.expand(-1);
1873 
1874  typename BoolLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1875  for (; iter; ++iter) {
1876 
1877  ijk = iter.getCoord();
1878 
1879  if (bbox.isInside(ijk)) {
1880  distance = FloatValueT(evalVoxelDist(ijk, *distLeafPt, *indexLeafPt, maskLeaf,
1881  primitives, closestPrim));
1882  } else {
1883  distance = FloatValueT(evalVoxelDist(ijk, distAcc, indexAcc, maskAcc,
1884  primitives, closestPrim));
1885  }
1886 
1887  pos = iter.pos();
1888 
1889  inside = distLeafPt->getValue(pos) < FloatValueT(0.0);
1890 
1891  if (!inside && distance < mExteriorBandWidth) {
1892  distLeafPt->setValueOn(pos, distance);
1893  indexLeafPt->setValueOn(pos, closestPrim);
1894  } else if (inside && distance < mInteriorBandWidth) {
1895  distLeafPt->setValueOn(pos, -distance);
1896  indexLeafPt->setValueOn(pos, closestPrim);
1897  } else {
1898  continue;
1899  }
1900 
1901  for (Int32 i = 0; i < 6; ++i) {
1902  newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[i]);
1903  }
1904  }
1905  }
1906 }
1907 
1908 
1909 template<typename FloatTreeT>
1910 double
1912  const Coord& ijk,
1913  FloatAccessorT& distAcc,
1914  IntAccessorT& indexAcc,
1915  BoolAccessorT& maskAcc,
1916  std::vector<Int32>& prims,
1917  Int32& closestPrim) const
1918 {
1919  FloatValueT tmpDist, minDist = std::numeric_limits<FloatValueT>::max();
1920  prims.clear();
1921 
1922  // Collect primitive indices from active neighbors and min distance.
1923  Coord n_ijk;
1924  for (Int32 n = 0; n < 18; ++n) {
1925  n_ijk = ijk + util::COORD_OFFSETS[n];
1926  if (!maskAcc.isValueOn(n_ijk) && distAcc.probeValue(n_ijk, tmpDist)) {
1927  prims.push_back(indexAcc.getValue(n_ijk));
1928  tmpDist = std::abs(tmpDist);
1929  if (tmpDist < minDist) minDist = tmpDist;
1930  }
1931  }
1932 
1933  // Calc. this voxels distance to the closest primitive.
1934  tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1935 
1936  // Forces the gradient to be monotonic for non-manifold
1937  // polygonal models with self-intersections.
1938  return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1939 }
1940 
1941 
1942 // Leaf specialized version.
1943 template<typename FloatTreeT>
1944 double
1945 ExpandNB<FloatTreeT>::evalVoxelDist(
1946  const Coord& ijk,
1947  FloatLeafT& distLeaf,
1948  IntLeafT& indexLeaf,
1949  BoolLeafT& maskLeaf,
1950  std::vector<Int32>& prims,
1951  Int32& closestPrim) const
1952 {
1953  FloatValueT tmpDist, minDist = std::numeric_limits<FloatValueT>::max();
1954  prims.clear();
1955 
1956  Index pos;
1957  for (Int32 n = 0; n < 18; ++n) {
1958  pos = FloatLeafT::coordToOffset(ijk + util::COORD_OFFSETS[n]);
1959  if (!maskLeaf.isValueOn(pos) && distLeaf.probeValue(pos, tmpDist)) {
1960  prims.push_back(indexLeaf.getValue(pos));
1961  tmpDist = std::abs(tmpDist);
1962  if (tmpDist < minDist) minDist = tmpDist;
1963  }
1964  }
1965 
1966  tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1967  return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1968 }
1969 
1970 
1971 template<typename FloatTreeT>
1972 double
1973 ExpandNB<FloatTreeT>::closestPrimDist(const Coord& ijk,
1974  std::vector<Int32>& prims, Int32& closestPrim) const
1975 {
1976  std::sort(prims.begin(), prims.end());
1977 
1978  Int32 lastPrim = -1;
1979  Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
1980  double primDist, tmpDist, dist = std::numeric_limits<double>::max();
1981 
1982  for (size_t n = 0, N = prims.size(); n < N; ++n) {
1983  if (prims[n] == lastPrim) continue;
1984 
1985  lastPrim = prims[n];
1986 
1987  const Vec4I& verts = mPolygonList[lastPrim];
1988 
1989  // Evaluate first triangle
1990  const Vec3d a(mPointList[verts[0]]);
1991  const Vec3d b(mPointList[verts[1]]);
1992  const Vec3d c(mPointList[verts[2]]);
1993 
1994  primDist = (voxelCenter -
1995  closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
1996 
1997  // Split-up quad into a second triangle and calac distance.
1998  if (util::INVALID_IDX != verts[3]) {
1999  const Vec3d d(mPointList[verts[3]]);
2000 
2001  tmpDist = (voxelCenter -
2002  closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw)).lengthSqr();
2003 
2004  if (tmpDist < primDist) primDist = tmpDist;
2005  }
2006 
2007  if (primDist < dist) {
2008  dist = primDist;
2009  closestPrim = lastPrim;
2010  }
2011  }
2012 
2013  return std::sqrt(dist) * double(mVoxelSize);
2014 }
2015 
2016 
2017 template<typename FloatTreeT>
2018 void
2020 {
2021  mNewDistTree.merge(rhs.mNewDistTree);
2022  mNewIndexTree.merge(rhs.mNewIndexTree);
2023  mNewMaskTree.merge(rhs.mNewMaskTree);
2024 }
2025 
2026 
2028 
2029 
2030 template<typename ValueType>
2032 {
2033  SqrtAndScaleOp(ValueType voxelSize, bool unsignedDist = false)
2034  : mVoxelSize(voxelSize)
2035  , mUnsigned(unsignedDist)
2036  {
2037  }
2038 
2039  template <typename LeafNodeType>
2040  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2041  {
2042  ValueType w[2];
2043  w[0] = mVoxelSize;
2044  w[1] = -mVoxelSize;
2045 
2046  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2047  for (; iter; ++iter) {
2048  ValueType& val = const_cast<ValueType&>(iter.getValue());
2049  val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
2050  }
2051  }
2052 
2053 private:
2054  ValueType mVoxelSize;
2055  const bool mUnsigned;
2056 };
2057 
2058 
2059 template<typename ValueType>
2061 {
2062  VoxelSignOp(ValueType exBandWidth, ValueType inBandWidth)
2063  : mExBandWidth(exBandWidth)
2064  , mInBandWidth(inBandWidth)
2065  {
2066  }
2067 
2068  template <typename LeafNodeType>
2069  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2070  {
2071  ValueType bgValues[2];
2072  bgValues[0] = mExBandWidth;
2073  bgValues[1] = -mInBandWidth;
2074 
2075  typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
2076 
2077  for (; iter; ++iter) {
2078  ValueType& val = const_cast<ValueType&>(iter.getValue());
2079  val = bgValues[int(val < ValueType(0.0))];
2080  }
2081  }
2082 
2083 private:
2084  ValueType mExBandWidth, mInBandWidth;
2085 };
2086 
2087 
2088 template<typename ValueType>
2089 struct TrimOp
2090 {
2091  TrimOp(ValueType exBandWidth, ValueType inBandWidth)
2092  : mExBandWidth(exBandWidth)
2093  , mInBandWidth(inBandWidth)
2094  {
2095  }
2096 
2097  template <typename LeafNodeType>
2098  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2099  {
2100  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2101 
2102  for (; iter; ++iter) {
2103  ValueType& val = const_cast<ValueType&>(iter.getValue());
2104  const bool inside = val < ValueType(0.0);
2105 
2106  if (inside && !(val > -mInBandWidth)) {
2107  val = -mInBandWidth;
2108  iter.setValueOff();
2109  } else if (!inside && !(val < mExBandWidth)) {
2110  val = mExBandWidth;
2111  iter.setValueOff();
2112  }
2113  }
2114  }
2115 
2116 private:
2117  ValueType mExBandWidth, mInBandWidth;
2118 };
2119 
2120 
2121 template<typename ValueType>
2122 struct OffsetOp
2123 {
2124  OffsetOp(ValueType offset): mOffset(offset) {}
2125 
2126  void resetOffset(ValueType offset) { mOffset = offset; }
2127 
2128  template <typename LeafNodeType>
2129  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2130  {
2131  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2132  for (; iter; ++iter) {
2133  ValueType& val = const_cast<ValueType&>(iter.getValue());
2134  val += mOffset;
2135  }
2136  }
2137 
2138 private:
2139  ValueType mOffset;
2140 };
2141 
2142 
2143 template<typename GridType, typename ValueType>
2144 struct RenormOp
2145 {
2147  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
2150 
2151  RenormOp(GridType& grid, LeafManagerType& leafs, ValueType voxelSize, ValueType cfl = 1.0)
2152  : mGrid(grid)
2153  , mLeafs(leafs)
2154  , mVoxelSize(voxelSize)
2155  , mCFL(cfl)
2156  {
2157  }
2158 
2159  void resetCFL(ValueType cfl) { mCFL = cfl; }
2160 
2161  template <typename LeafNodeType>
2162  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2163  {
2164  const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
2165  Stencil stencil(mGrid);
2166  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
2167 
2168  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2169  for (; iter; ++iter) {
2170  stencil.moveTo(iter);
2171 
2172  const ValueType normSqGradPhi =
2174 
2175  const ValueType phi0 = iter.getValue();
2176  const ValueType diff = math::Sqrt(normSqGradPhi) * invDx - one;
2177  const ValueType S = phi0 / (math::Sqrt(math::Pow2(phi0) + normSqGradPhi));
2178 
2179  buffer.setValue(iter.pos(), phi0 - dt * S * diff);
2180  }
2181  }
2182 
2183 private:
2184  GridType& mGrid;
2185  LeafManagerType& mLeafs;
2186  ValueType mVoxelSize, mCFL;
2187 };
2188 
2189 
2190 template<typename TreeType, typename ValueType>
2191 struct MinOp
2192 {
2195 
2196  MinOp(LeafManagerType& leafs): mLeafs(leafs) {}
2197 
2198  template <typename LeafNodeType>
2199  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2200  {
2201  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
2202  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2203 
2204  for (; iter; ++iter) {
2205  ValueType& val = const_cast<ValueType&>(iter.getValue());
2206  val = std::min(val, buffer.getValue(iter.pos()));
2207  }
2208  }
2209 
2210 private:
2211  LeafManagerType& mLeafs;
2212 };
2213 
2214 
2215 template<typename TreeType, typename ValueType>
2217 {
2220 
2221  MergeBufferOp(LeafManagerType& leafs, size_t bufferIndex = 1)
2222  : mLeafs(leafs)
2223  , mBufferIndex(bufferIndex)
2224  {
2225  }
2226 
2227  template <typename LeafNodeType>
2228  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2229  {
2230  BufferType& buffer = mLeafs.getBuffer(leafIndex, mBufferIndex);
2231  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2232  Index offset;
2233 
2234  for (; iter; ++iter) {
2235  offset = iter.pos();
2236  leaf.setValueOnly(offset, buffer.getValue(offset));
2237  }
2238  }
2239 
2240 private:
2241  LeafManagerType& mLeafs;
2242  const size_t mBufferIndex;
2243 };
2244 
2245 
2246 template<typename TreeType>
2248 {
2250  typedef typename TreeType::LeafNodeType LeafNodeT;
2251 
2252  LeafTopologyDiffOp(TreeType& tree) : mAcc(tree) { }
2253 
2254  template <typename LeafNodeType>
2255  void operator()(LeafNodeType &leaf, size_t) const
2256  {
2257  const LeafNodeT* rhsLeaf = mAcc.probeConstLeaf(leaf.origin());
2258  if (rhsLeaf) leaf.topologyDifference(*rhsLeaf, false);
2259  }
2260 
2261 private:
2262  AccessorT mAcc;
2263 };
2264 
2265 } // internal namespace
2266 
2267 
2269 
2270 
2271 // MeshToVolume
2272 
2273 template<typename FloatGridT, typename InterruptT>
2275  openvdb::math::Transform::Ptr& transform, int conversionFlags,
2276  InterruptT *interrupter, int signSweeps)
2277  : mTransform(transform)
2278  , mConversionFlags(conversionFlags)
2279  , mSignSweeps(signSweeps)
2280  , mInterrupter(interrupter)
2281 {
2282  clear();
2283  mSignSweeps = std::min(mSignSweeps, 1);
2284 }
2285 
2286 
2287 template<typename FloatGridT, typename InterruptT>
2288 void
2290 {
2291  mDistGrid = FloatGridT::create(std::numeric_limits<FloatValueT>::max());
2292  mIndexGrid = IntGridT::create(Int32(util::INVALID_IDX));
2293  mIntersectingVoxelsGrid = BoolGridT::create(false);
2294 }
2295 
2296 
2297 template<typename FloatGridT, typename InterruptT>
2298 inline void
2300  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2301  FloatValueT exBandWidth, FloatValueT inBandWidth)
2302 {
2303  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
2304  exBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), exBandWidth);
2305  inBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), inBandWidth);
2306  const FloatValueT vs = FloatValueT(mTransform->voxelSize()[0]);
2307 
2308  // Convert from index space units to world-space units. To fill the
2309  // interior, inBandWidth is passed FLOAT_MAX. Don't multiply with vs if so.
2310  exBandWidth *= vs;
2311  if (inBandWidth < std::numeric_limits<FloatValueT>::max()) {
2312  inBandWidth *= vs;
2313  }
2314 
2315  doConvert(pointList, polygonList, exBandWidth, inBandWidth);
2316  mDistGrid->setGridClass(GRID_LEVEL_SET);
2317 }
2318 
2319 
2320 template<typename FloatGridT, typename InterruptT>
2321 inline void
2323  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2324  FloatValueT exBandWidth)
2325 {
2326  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
2327  exBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), exBandWidth);
2328  const FloatValueT vs = FloatValueT(mTransform->voxelSize()[0]);
2329  doConvert(pointList, polygonList, vs * exBandWidth, 0.0, true);
2330  mDistGrid->setGridClass(GRID_UNKNOWN);
2331 }
2332 
2333 
2334 template<typename FloatGridT, typename InterruptT>
2335 void
2337  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2338  FloatValueT exBandWidth, FloatValueT inBandWidth, bool unsignedDistField)
2339 {
2340  mDistGrid->setTransform(mTransform);
2341  mIndexGrid->setTransform(mTransform);
2342  const bool rawData = OUTPUT_RAW_DATA & mConversionFlags;
2343 
2344  // Note that inBandWidth is allowed to be infinite when filling the interior.
2345  if (!boost::math::isfinite(exBandWidth) || boost::math::isnan(inBandWidth)) {
2346  std::stringstream msg;
2347  msg << "Illegal narrow band width: exterior = " << exBandWidth
2348  << ", interior = " << inBandWidth;
2349  OPENVDB_THROW(ValueError, msg.str());
2350  }
2351 
2352 
2353  // The progress estimates given to the interrupter are based on the
2354  // observed average time for each stage and therefore not alway
2355  // accurate. The goal is to give some progression feedback to the user.
2356 
2357  if (wasInterrupted(1)) return;
2358 
2359  // Voxelize mesh
2360  {
2361  internal::MeshVoxelizer<FloatTreeT, InterruptT>
2362  voxelizer(pointList, polygonList, mInterrupter);
2363 
2364  voxelizer.run();
2365 
2366  if (wasInterrupted(18)) return;
2367 
2368  mDistGrid->tree().merge(voxelizer.sqrDistTree());
2369  mIndexGrid->tree().merge(voxelizer.primIndexTree());
2370  mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
2371  }
2372 
2373  if (!unsignedDistField) {
2374  // Determine the inside/outside state for the narrow band of voxels.
2375  {
2376  // Slices up the volume and label the exterior contour of each slice in parallel.
2377  internal::ContourTracer<FloatTreeT, InterruptT> trace(
2378  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2379  for (int i = 0; i < mSignSweeps; ++i) {
2380 
2381  if (wasInterrupted(19)) return;
2382 
2383  trace.run();
2384 
2385  if (wasInterrupted(24)) return;
2386 
2387  // Propagate sign information between the slices.
2388  BoolTreeT signMaskTree(false);
2389  {
2390  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2391  internal::SignMask<FloatTreeT, InterruptT> signMaskOp(leafs,
2392  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2393  signMaskOp.run();
2394  signMaskTree.merge(signMaskOp.signMaskTree());
2395  }
2396 
2397  if (wasInterrupted(25)) return;
2398 
2399  while (true) {
2400  tree::LeafManager<BoolTreeT> leafs(signMaskTree);
2401  if(leafs.leafCount() == 0) break;
2402 
2403  internal::PropagateSign<FloatTreeT, InterruptT> sign(leafs,
2404  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2405 
2406  sign.run();
2407 
2408  signMaskTree.clear();
2409  signMaskTree.merge(sign.signMaskTree());
2410  }
2411  }
2412  }
2413 
2414 
2415  if (wasInterrupted(28)) return;
2416  {
2417  tree::LeafManager<BoolTreeT> leafs(mIntersectingVoxelsGrid->tree());
2418 
2419  // Determine the sign of the mesh intersecting voxels.
2420  internal::IntersectingVoxelSign<FloatTreeT> sign(pointList, polygonList,
2421  mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2422 
2423  sign.run();
2424 
2425  if (wasInterrupted(34)) return;
2426 
2427  // Remove mesh intersecting voxels that where set by rasterizing
2428  // self-intersecting portions of the mesh.
2429  internal::IntersectingVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2430  mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2431  cleaner.run();
2432  }
2433 
2434  // Remove shell voxels that where set by rasterizing
2435  // self-intersecting portions of the mesh.
2436  {
2437  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2438 
2439  internal::ShellVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2440  leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
2441 
2442  cleaner.run();
2443  }
2444 
2445  if (wasInterrupted(38)) return;
2446 
2447  } else { // if unsigned dist. field
2448  inBandWidth = FloatValueT(0.0);
2449  }
2450 
2451  if (mDistGrid->activeVoxelCount() == 0) {
2452  tools::changeBackground(mDistGrid->tree(), exBandWidth);
2453  return;
2454  }
2455 
2456  mIntersectingVoxelsGrid->clear();
2457  const FloatValueT voxelSize = FloatValueT(mTransform->voxelSize()[0]);
2458 
2459  { // Transform values (world space scaling etc.)
2460  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2461  leafs.foreach(internal::SqrtAndScaleOp<FloatValueT>(voxelSize, unsignedDistField));
2462  }
2463 
2464  if (wasInterrupted(40)) return;
2465 
2466  if (!unsignedDistField) { // Propagate sign information to inactive values.
2467  mDistGrid->tree().root().setBackground(exBandWidth, /*updateChildNodes=*/false);
2468  tools::signedFloodFillWithValues(mDistGrid->tree(), exBandWidth, -inBandWidth);
2469  }
2470 
2471  if (wasInterrupted(46)) return;
2472 
2473  // Narrow-band dilation
2474  const FloatValueT minWidth = FloatValueT(voxelSize * 2.0);
2475  if (inBandWidth > minWidth || exBandWidth > minWidth) {
2476 
2477  // Create the initial voxel mask.
2478  BoolTreeT maskTree(false);
2479  maskTree.topologyUnion(mDistGrid->tree());
2480 
2481  if (wasInterrupted(48)) return;
2482 
2483  internal::LeafTopologyDiffOp<FloatTreeT> diffOp(mDistGrid->tree());
2484  openvdb::tools::dilateVoxels(maskTree);
2485 
2486  unsigned maxIterations = std::numeric_limits<unsigned>::max();
2487  float progress = 48, step = 0.0;
2488  // progress estimation..
2489  double estimated =
2490  2.0 * std::ceil((std::max(inBandWidth, exBandWidth) - minWidth) / voxelSize);
2491  if (estimated < double(maxIterations)) {
2492  maxIterations = unsigned(estimated);
2493  step = 42.f / float(maxIterations);
2494  }
2495 
2496  unsigned count = 0;
2497  while (true) {
2498 
2499  if (wasInterrupted(int(progress))) return;
2500 
2501  tree::LeafManager<BoolTreeT> leafs(maskTree);
2502 
2503  if (leafs.leafCount() == 0) break;
2504 
2505  leafs.foreach(diffOp);
2506 
2507  internal::ExpandNB<FloatTreeT> expand(
2508  leafs, mDistGrid->tree(), mIndexGrid->tree(), maskTree,
2509  exBandWidth, inBandWidth, voxelSize, pointList, polygonList);
2510 
2511  expand.run();
2512 
2513  if ((++count) >= maxIterations) break;
2514  progress += step;
2515  }
2516  }
2517 
2518  if (!bool(GENERATE_PRIM_INDEX_GRID & mConversionFlags)) mIndexGrid->clear();
2519 
2520  if (wasInterrupted(80)) return;
2521 
2522  // Renormalize distances to smooth out bumps caused by self-intersecting
2523  // and overlapping portions of the mesh and renormalize the level set.
2524  if (!unsignedDistField && !rawData) {
2525 
2526  tools::pruneLevelSet(mDistGrid->tree(), exBandWidth, -inBandWidth);
2527 
2528  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree(), 1);
2529 
2530  const FloatValueT offset = FloatValueT(0.8 * voxelSize);
2531  if (wasInterrupted(82)) return;
2532 
2533  internal::OffsetOp<FloatValueT> offsetOp(-offset);
2534 
2535  leafs.foreach(offsetOp);
2536 
2537  if (wasInterrupted(84)) return;
2538 
2539  leafs.foreach(internal::RenormOp<FloatGridT, FloatValueT>(*mDistGrid, leafs, voxelSize));
2540 
2541  leafs.foreach(internal::MinOp<FloatTreeT, FloatValueT>(leafs));
2542 
2543  if (wasInterrupted(95)) return;
2544 
2545  offsetOp.resetOffset(offset - internal::Tolerance<FloatValueT>::epsilon());
2546  leafs.foreach(offsetOp);
2547  }
2548 
2549  if (wasInterrupted(98)) return;
2550 
2551  const FloatValueT minTrimWidth = FloatValueT(voxelSize * 4.0);
2552  if (inBandWidth < minTrimWidth || exBandWidth < minTrimWidth) {
2553 
2554  // If the narrow band was not expanded, we might need to trim off
2555  // some of the active voxels in order to respect the narrow band limits.
2556  // (The mesh voxelization step generates some extra 'shell' voxels)
2557 
2558  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2559  leafs.foreach(internal::TrimOp<FloatValueT>(
2560  exBandWidth, unsignedDistField ? exBandWidth : inBandWidth));
2561 
2562  tools::pruneLevelSet(mDistGrid->tree(), exBandWidth, unsignedDistField ? -exBandWidth : -inBandWidth);
2563  }
2564 }
2565 
2566 
2568 
2569 
2571 template<typename GridType>
2572 inline typename boost::enable_if<boost::is_floating_point<typename GridType::ValueType>,
2573 typename GridType::Ptr>::type
2575  const openvdb::math::Transform& xform,
2576  const std::vector<Vec3s>& points,
2577  const std::vector<Vec3I>& triangles,
2578  const std::vector<Vec4I>& quads,
2579  float exBandWidth,
2580  float inBandWidth,
2581  bool unsignedDistanceField = false)
2582 {
2583  std::vector<Vec3s> indexSpacePoints(points.size());
2584 
2585  { // Copy and transform (required for MeshToVolume) points to grid space.
2586  internal::PointTransform ptnXForm(points, indexSpacePoints, xform);
2587  ptnXForm.run();
2588  }
2589 
2590  // Copy primitives
2591  std::vector<Vec4I> primitives(triangles.size() + quads.size());
2592 
2593  for (size_t n = 0, N = triangles.size(); n < N; ++n) {
2594  Vec4I& prim = primitives[n];
2595  const Vec3I& triangle = triangles[n];
2596  prim[0] = triangle[0];
2597  prim[1] = triangle[1];
2598  prim[2] = triangle[2];
2599  prim[3] = util::INVALID_IDX;
2600  }
2601 
2602  for (size_t n = 0, N = quads.size(); n < N; ++n) {
2603  primitives[n + triangles.size()] = quads[n];
2604  }
2605 
2606  typename GridType::ValueType exWidth(exBandWidth);
2607  typename GridType::ValueType inWidth(inBandWidth);
2608 
2609 
2610  math::Transform::Ptr transform = xform.copy();
2611  MeshToVolume<GridType> vol(transform);
2612 
2613  if (!unsignedDistanceField) {
2614  vol.convertToLevelSet(indexSpacePoints, primitives, exWidth, inWidth);
2615  } else {
2616  vol.convertToUnsignedDistanceField(indexSpacePoints, primitives, exWidth);
2617  }
2618 
2619  return vol.distGridPtr();
2620 }
2621 
2622 
2625 template<typename GridType>
2626 inline typename boost::disable_if<boost::is_floating_point<typename GridType::ValueType>,
2627 typename GridType::Ptr>::type
2629  const math::Transform& /*xform*/,
2630  const std::vector<Vec3s>& /*points*/,
2631  const std::vector<Vec3I>& /*triangles*/,
2632  const std::vector<Vec4I>& /*quads*/,
2633  float /*exBandWidth*/,
2634  float /*inBandWidth*/,
2635  bool /*unsignedDistanceField*/ = false)
2636 {
2638  "mesh to volume conversion is supported only for scalar, floating-point grids");
2639 }
2640 
2641 
2643 
2644 
2645 template<typename GridType>
2646 inline typename GridType::Ptr
2648  const openvdb::math::Transform& xform,
2649  const std::vector<Vec3s>& points,
2650  const std::vector<Vec3I>& triangles,
2651  float halfWidth)
2652 {
2653  std::vector<Vec4I> quads(0);
2654  return doMeshConversion<GridType>(xform, points, triangles, quads,
2655  halfWidth, halfWidth);
2656 }
2657 
2658 
2659 template<typename GridType>
2660 inline typename GridType::Ptr
2662  const openvdb::math::Transform& xform,
2663  const std::vector<Vec3s>& points,
2664  const std::vector<Vec4I>& quads,
2665  float halfWidth)
2666 {
2667  std::vector<Vec3I> triangles(0);
2668  return doMeshConversion<GridType>(xform, points, triangles, quads,
2669  halfWidth, halfWidth);
2670 }
2671 
2672 
2673 template<typename GridType>
2674 inline typename GridType::Ptr
2676  const openvdb::math::Transform& xform,
2677  const std::vector<Vec3s>& points,
2678  const std::vector<Vec3I>& triangles,
2679  const std::vector<Vec4I>& quads,
2680  float halfWidth)
2681 {
2682  return doMeshConversion<GridType>(xform, points, triangles, quads,
2683  halfWidth, halfWidth);
2684 }
2685 
2686 
2687 template<typename GridType>
2688 inline typename GridType::Ptr
2690  const openvdb::math::Transform& xform,
2691  const std::vector<Vec3s>& points,
2692  const std::vector<Vec3I>& triangles,
2693  const std::vector<Vec4I>& quads,
2694  float exBandWidth,
2695  float inBandWidth)
2696 {
2697  return doMeshConversion<GridType>(xform, points, triangles,
2698  quads, exBandWidth, inBandWidth);
2699 }
2700 
2701 
2702 template<typename GridType>
2703 inline typename GridType::Ptr
2705  const openvdb::math::Transform& xform,
2706  const std::vector<Vec3s>& points,
2707  const std::vector<Vec3I>& triangles,
2708  const std::vector<Vec4I>& quads,
2709  float bandWidth)
2710 {
2711  return doMeshConversion<GridType>(xform, points, triangles, quads,
2712  bandWidth, bandWidth, true);
2713 }
2714 
2715 
2717 
2718 
2719 // Required by several of the tree nodes
2720 inline std::ostream&
2721 operator<<(std::ostream& ostr, const MeshToVoxelEdgeData::EdgeData& rhs)
2722 {
2723  ostr << "{[ " << rhs.mXPrim << ", " << rhs.mXDist << "]";
2724  ostr << " [ " << rhs.mYPrim << ", " << rhs.mYDist << "]";
2725  ostr << " [ " << rhs.mZPrim << ", " << rhs.mZDist << "]}";
2726  return ostr;
2727 }
2728 
2729 // Required by math::Abs
2730 inline MeshToVoxelEdgeData::EdgeData
2732 {
2733  return x;
2734 }
2735 
2736 
2738 
2739 
2741 {
2742 public:
2743 
2744  GenEdgeData(
2745  const std::vector<Vec3s>& pointList,
2746  const std::vector<Vec4I>& polygonList);
2747 
2748  void run(bool threaded = true);
2749 
2750  GenEdgeData(GenEdgeData& rhs, tbb::split);
2751  inline void operator() (const tbb::blocked_range<size_t> &range);
2752  inline void join(GenEdgeData& rhs);
2753 
2754  inline TreeType& tree() { return mTree; }
2755 
2756 private:
2757  void operator=(const GenEdgeData&) {}
2758 
2759  struct Primitive { Vec3d a, b, c, d; Int32 index; };
2760 
2761  template<bool IsQuad>
2762  inline void voxelize(const Primitive&);
2763 
2764  template<bool IsQuad>
2765  inline bool evalPrimitive(const Coord&, const Primitive&);
2766 
2767  inline bool rayTriangleIntersection( const Vec3d& origin, const Vec3d& dir,
2768  const Vec3d& a, const Vec3d& b, const Vec3d& c, double& t);
2769 
2770 
2771  TreeType mTree;
2772  Accessor mAccessor;
2773 
2774  const std::vector<Vec3s>& mPointList;
2775  const std::vector<Vec4I>& mPolygonList;
2776 
2777  // Used internally for acceleration
2778  typedef TreeType::ValueConverter<Int32>::Type IntTreeT;
2779  IntTreeT mLastPrimTree;
2780  tree::ValueAccessor<IntTreeT> mLastPrimAccessor;
2781 }; // class MeshToVoxelEdgeData::GenEdgeData
2782 
2783 
2784 inline
2786  const std::vector<Vec3s>& pointList,
2787  const std::vector<Vec4I>& polygonList)
2788  : mTree(EdgeData())
2789  , mAccessor(mTree)
2790  , mPointList(pointList)
2791  , mPolygonList(polygonList)
2792  , mLastPrimTree(Int32(util::INVALID_IDX))
2793  , mLastPrimAccessor(mLastPrimTree)
2794 {
2795 }
2796 
2797 
2798 inline
2800  : mTree(EdgeData())
2801  , mAccessor(mTree)
2802  , mPointList(rhs.mPointList)
2803  , mPolygonList(rhs.mPolygonList)
2804  , mLastPrimTree(Int32(util::INVALID_IDX))
2805  , mLastPrimAccessor(mLastPrimTree)
2806 {
2807 }
2808 
2809 
2810 inline void
2812 {
2813  if (threaded) {
2814  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
2815  } else {
2816  (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
2817  }
2818 }
2819 
2820 
2821 inline void
2823 {
2824  typedef TreeType::RootNodeType RootNodeType;
2825  typedef RootNodeType::NodeChainType NodeChainType;
2826  BOOST_STATIC_ASSERT(boost::mpl::size<NodeChainType>::value > 1);
2827  typedef boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
2828 
2829  Coord ijk;
2830  Index offset;
2831 
2832  rhs.mTree.clearAllAccessors();
2833 
2834  TreeType::LeafIter leafIt = rhs.mTree.beginLeaf();
2835  for ( ; leafIt; ++leafIt) {
2836  ijk = leafIt->origin();
2837 
2838  TreeType::LeafNodeType* lhsLeafPt = mTree.probeLeaf(ijk);
2839 
2840  if (!lhsLeafPt) {
2841 
2842  mAccessor.addLeaf(rhs.mAccessor.probeLeaf(ijk));
2843  InternalNodeType* node = rhs.mAccessor.getNode<InternalNodeType>();
2844  node->stealNode<TreeType::LeafNodeType>(ijk, EdgeData(), false);
2845  rhs.mAccessor.clear();
2846 
2847  } else {
2848 
2849  TreeType::LeafNodeType::ValueOnCIter it = leafIt->cbeginValueOn();
2850  for ( ; it; ++it) {
2851 
2852  offset = it.pos();
2853  const EdgeData& rhsValue = it.getValue();
2854 
2855  if (!lhsLeafPt->isValueOn(offset)) {
2856  lhsLeafPt->setValueOn(offset, rhsValue);
2857  } else {
2858 
2859  EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
2860 
2861  if (rhsValue.mXDist < lhsValue.mXDist) {
2862  lhsValue.mXDist = rhsValue.mXDist;
2863  lhsValue.mXPrim = rhsValue.mXPrim;
2864  }
2865 
2866  if (rhsValue.mYDist < lhsValue.mYDist) {
2867  lhsValue.mYDist = rhsValue.mYDist;
2868  lhsValue.mYPrim = rhsValue.mYPrim;
2869  }
2870 
2871  if (rhsValue.mZDist < lhsValue.mZDist) {
2872  lhsValue.mZDist = rhsValue.mZDist;
2873  lhsValue.mZPrim = rhsValue.mZPrim;
2874  }
2875 
2876  }
2877  } // end value iteration
2878  }
2879  } // end leaf iteration
2880 }
2881 
2882 
2883 inline void
2884 MeshToVoxelEdgeData::GenEdgeData::operator()(const tbb::blocked_range<size_t> &range)
2885 {
2886  Primitive prim;
2887 
2888  for (size_t n = range.begin(); n < range.end(); ++n) {
2889 
2890  const Vec4I& verts = mPolygonList[n];
2891 
2892  prim.index = Int32(n);
2893  prim.a = Vec3d(mPointList[verts[0]]);
2894  prim.b = Vec3d(mPointList[verts[1]]);
2895  prim.c = Vec3d(mPointList[verts[2]]);
2896 
2897  if (util::INVALID_IDX != verts[3]) {
2898  prim.d = Vec3d(mPointList[verts[3]]);
2899  voxelize<true>(prim);
2900  } else {
2901  voxelize<false>(prim);
2902  }
2903  }
2904 }
2905 
2906 
2907 template<bool IsQuad>
2908 inline void
2909 MeshToVoxelEdgeData::GenEdgeData::voxelize(const Primitive& prim)
2910 {
2911  std::deque<Coord> coordList;
2912  Coord ijk, nijk;
2913 
2914  ijk = util::nearestCoord(prim.a);
2915  coordList.push_back(ijk);
2916 
2917  evalPrimitive<IsQuad>(ijk, prim);
2918 
2919  while (!coordList.empty()) {
2920 
2921  ijk = coordList.back();
2922  coordList.pop_back();
2923 
2924  for (Int32 i = 0; i < 26; ++i) {
2925  nijk = ijk + util::COORD_OFFSETS[i];
2926 
2927  if (prim.index != mLastPrimAccessor.getValue(nijk)) {
2928  mLastPrimAccessor.setValue(nijk, prim.index);
2929  if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
2930  }
2931  }
2932  }
2933 }
2934 
2935 
2936 template<bool IsQuad>
2937 inline bool
2938 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(const Coord& ijk, const Primitive& prim)
2939 {
2940  Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
2941  bool intersecting = false;
2942  double t;
2943 
2944  EdgeData edgeData;
2945  mAccessor.probeValue(ijk, edgeData);
2946 
2947  // Evaluate first triangle
2948  double dist = (org -
2949  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, org, uvw)).lengthSqr();
2950 
2951  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
2952  if (t < edgeData.mXDist) {
2953  edgeData.mXDist = float(t);
2954  edgeData.mXPrim = prim.index;
2955  intersecting = true;
2956  }
2957  }
2958 
2959  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
2960  if (t < edgeData.mYDist) {
2961  edgeData.mYDist = float(t);
2962  edgeData.mYPrim = prim.index;
2963  intersecting = true;
2964  }
2965  }
2966 
2967  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
2968  if (t < edgeData.mZDist) {
2969  edgeData.mZDist = float(t);
2970  edgeData.mZPrim = prim.index;
2971  intersecting = true;
2972  }
2973  }
2974 
2975  if (IsQuad) {
2976  // Split quad into a second triangle and calculate distance.
2977  double secondDist = (org -
2978  closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, org, uvw)).lengthSqr();
2979 
2980  if (secondDist < dist) dist = secondDist;
2981 
2982  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
2983  if (t < edgeData.mXDist) {
2984  edgeData.mXDist = float(t);
2985  edgeData.mXPrim = prim.index;
2986  intersecting = true;
2987  }
2988  }
2989 
2990  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
2991  if (t < edgeData.mYDist) {
2992  edgeData.mYDist = float(t);
2993  edgeData.mYPrim = prim.index;
2994  intersecting = true;
2995  }
2996  }
2997 
2998  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
2999  if (t < edgeData.mZDist) {
3000  edgeData.mZDist = float(t);
3001  edgeData.mZPrim = prim.index;
3002  intersecting = true;
3003  }
3004  }
3005  }
3006 
3007  if (intersecting) mAccessor.setValue(ijk, edgeData);
3008 
3009  return (dist < 0.86602540378443861);
3010 }
3011 
3012 
3013 inline bool
3014 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3015  const Vec3d& origin, const Vec3d& dir,
3016  const Vec3d& a, const Vec3d& b, const Vec3d& c,
3017  double& t)
3018 {
3019  // Check if ray is parallel with triangle
3020 
3021  Vec3d e1 = b - a;
3022  Vec3d e2 = c - a;
3023  Vec3d s1 = dir.cross(e2);
3024 
3025  double divisor = s1.dot(e1);
3026  if (!(std::abs(divisor) > 0.0)) return false;
3027 
3028  // Compute barycentric coordinates
3029 
3030  double inv_divisor = 1.0 / divisor;
3031  Vec3d d = origin - a;
3032  double b1 = d.dot(s1) * inv_divisor;
3033 
3034  if (b1 < 0.0 || b1 > 1.0) return false;
3035 
3036  Vec3d s2 = d.cross(e1);
3037  double b2 = dir.dot(s2) * inv_divisor;
3038 
3039  if (b2 < 0.0 || (b1 + b2) > 1.0) return false;
3040 
3041  // Compute distance to intersection point
3042 
3043  t = e2.dot(s2) * inv_divisor;
3044  return (t < 0.0) ? false : true;
3045 }
3046 
3047 
3049 
3050 
3051 inline
3053  : mTree(EdgeData())
3054 {
3055 }
3056 
3057 
3058 inline void
3060  const std::vector<Vec3s>& pointList,
3061  const std::vector<Vec4I>& polygonList)
3062 {
3063  GenEdgeData converter(pointList, polygonList);
3064  converter.run();
3065 
3066  mTree.clear();
3067  mTree.merge(converter.tree());
3068 }
3069 
3070 
3071 inline void
3073  Accessor& acc,
3074  const Coord& ijk,
3075  std::vector<Vec3d>& points,
3076  std::vector<Index32>& primitives)
3077 {
3078  EdgeData data;
3079  Vec3d point;
3080 
3081  Coord coord = ijk;
3082 
3083  if (acc.probeValue(coord, data)) {
3084 
3085  if (data.mXPrim != util::INVALID_IDX) {
3086  point[0] = double(coord[0]) + data.mXDist;
3087  point[1] = double(coord[1]);
3088  point[2] = double(coord[2]);
3089 
3090  points.push_back(point);
3091  primitives.push_back(data.mXPrim);
3092  }
3093 
3094  if (data.mYPrim != util::INVALID_IDX) {
3095  point[0] = double(coord[0]);
3096  point[1] = double(coord[1]) + data.mYDist;
3097  point[2] = double(coord[2]);
3098 
3099  points.push_back(point);
3100  primitives.push_back(data.mYPrim);
3101  }
3102 
3103  if (data.mZPrim != util::INVALID_IDX) {
3104  point[0] = double(coord[0]);
3105  point[1] = double(coord[1]);
3106  point[2] = double(coord[2]) + data.mZDist;
3107 
3108  points.push_back(point);
3109  primitives.push_back(data.mZPrim);
3110  }
3111 
3112  }
3113 
3114  coord[0] += 1;
3115 
3116  if (acc.probeValue(coord, data)) {
3117 
3118  if (data.mYPrim != util::INVALID_IDX) {
3119  point[0] = double(coord[0]);
3120  point[1] = double(coord[1]) + data.mYDist;
3121  point[2] = double(coord[2]);
3122 
3123  points.push_back(point);
3124  primitives.push_back(data.mYPrim);
3125  }
3126 
3127  if (data.mZPrim != util::INVALID_IDX) {
3128  point[0] = double(coord[0]);
3129  point[1] = double(coord[1]);
3130  point[2] = double(coord[2]) + data.mZDist;
3131 
3132  points.push_back(point);
3133  primitives.push_back(data.mZPrim);
3134  }
3135  }
3136 
3137  coord[2] += 1;
3138 
3139  if (acc.probeValue(coord, data)) {
3140  if (data.mYPrim != util::INVALID_IDX) {
3141  point[0] = double(coord[0]);
3142  point[1] = double(coord[1]) + data.mYDist;
3143  point[2] = double(coord[2]);
3144 
3145  points.push_back(point);
3146  primitives.push_back(data.mYPrim);
3147  }
3148  }
3149 
3150  coord[0] -= 1;
3151 
3152  if (acc.probeValue(coord, data)) {
3153 
3154  if (data.mXPrim != util::INVALID_IDX) {
3155  point[0] = double(coord[0]) + data.mXDist;
3156  point[1] = double(coord[1]);
3157  point[2] = double(coord[2]);
3158 
3159  points.push_back(point);
3160  primitives.push_back(data.mXPrim);
3161  }
3162 
3163  if (data.mYPrim != util::INVALID_IDX) {
3164  point[0] = double(coord[0]);
3165  point[1] = double(coord[1]) + data.mYDist;
3166  point[2] = double(coord[2]);
3167 
3168  points.push_back(point);
3169  primitives.push_back(data.mYPrim);
3170  }
3171  }
3172 
3173 
3174  coord[1] += 1;
3175 
3176  if (acc.probeValue(coord, data)) {
3177 
3178  if (data.mXPrim != util::INVALID_IDX) {
3179  point[0] = double(coord[0]) + data.mXDist;
3180  point[1] = double(coord[1]);
3181  point[2] = double(coord[2]);
3182 
3183  points.push_back(point);
3184  primitives.push_back(data.mXPrim);
3185  }
3186  }
3187 
3188  coord[2] -= 1;
3189 
3190  if (acc.probeValue(coord, data)) {
3191 
3192  if (data.mXPrim != util::INVALID_IDX) {
3193  point[0] = double(coord[0]) + data.mXDist;
3194  point[1] = double(coord[1]);
3195  point[2] = double(coord[2]);
3196 
3197  points.push_back(point);
3198  primitives.push_back(data.mXPrim);
3199  }
3200 
3201  if (data.mZPrim != util::INVALID_IDX) {
3202  point[0] = double(coord[0]);
3203  point[1] = double(coord[1]);
3204  point[2] = double(coord[2]) + data.mZDist;
3205 
3206  points.push_back(point);
3207  primitives.push_back(data.mZPrim);
3208  }
3209  }
3210 
3211  coord[0] += 1;
3212 
3213  if (acc.probeValue(coord, data)) {
3214 
3215  if (data.mZPrim != util::INVALID_IDX) {
3216  point[0] = double(coord[0]);
3217  point[1] = double(coord[1]);
3218  point[2] = double(coord[2]) + data.mZDist;
3219 
3220  points.push_back(point);
3221  primitives.push_back(data.mZPrim);
3222  }
3223  }
3224 }
3225 
3226 
3227 } // namespace tools
3228 } // namespace OPENVDB_VERSION_NAME
3229 } // namespace openvdb
3230 
3231 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
3232 
3233 // Copyright (c) 2012-2014 DreamWorks Animation LLC
3234 // All rights reserved. This software is distributed under the
3235 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1630
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:1695
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1107
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
Definition: ValueAccessor.h:328
tree::ValueAccessor< TreeType > Accessor
Definition: MeshToVolume.h:315
~IntersectingVoxelCleaner()
Definition: MeshToVolume.h:1444
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:956
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1566
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1430
MeshToVolume(openvdb::math::Transform::Ptr &, int conversionFlags=0, InterruptT *interrupter=NULL, int signSweeps=1)
Definition: MeshToVolume.h:2274
void clear()
Definition: MeshToVolume.h:2289
GridType::Ptr meshToUnsignedDistanceField(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float bandWidth)
Convert a triangle and quad mesh to an unsigned distance field.
Definition: MeshToVolume.h:2704
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1568
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2194
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1724
IntersectingVoxelSign(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, FloatTreeT &distTree, IntTreeT &indexTree, BoolTreeT &intersectionTree, BoolLeafManager &leafs)
Definition: MeshToVolume.h:1303
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:378
TBB body object that removes intersecting voxels that were set via voxelization of self-intersecting ...
Definition: MeshToVolume.h:1427
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:260
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:2731
OPENVDB_STATIC_SPECIALIZATION void dilateVoxels(TreeType &tree, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically dilate all leaf-level active voxels in a tree using one of three nearest neighbor conne...
Definition: Morphology.h:774
EdgeData(float dist=1.0)
Definition: MeshToVolume.h:287
FloatTreeT::LeafNodeType DistLeafT
Definition: MeshToVolume.h:1432
OPENVDB_API const Index32 INVALID_IDX
Accessor getAccessor()
Definition: MeshToVolume.h:341
Definition: Mat.h:146
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:269
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:303
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:170
boost::enable_if< boost::is_floating_point< typename GridType::ValueType >, typename GridType::Ptr >::type doMeshConversion(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float exBandWidth, float inBandWidth, bool unsignedDistanceField=false)
Definition: MeshToVolume.h:2574
Index32 mXPrim
Definition: MeshToVolume.h:311
_RootNodeType RootNodeType
Definition: Tree.h:209
TBB body object that traversers all intersecting voxels (defined by the intersectingVoxelsGrid) and p...
Definition: MeshToVolume.h:1251
Extracts and stores voxel edge intersection data from a mesh.
Definition: MeshToVolume.h:279
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:208
boost::shared_ptr< Transform > Ptr
Definition: Transform.h:69
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:557
tree::LeafManager< typename GridType::TreeType > LeafManagerType
Definition: MeshToVolume.h:2148
void operator()(const tbb::blocked_range< int > &range) const
Definition: MeshToVolume.h:845
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:1431
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1434
float mYDist
Definition: MeshToVolume.h:310
void operator()(const tbb::blocked_range< size_t > &)
Definition: MeshToVolume.h:1834
tree::LeafManager< FloatTreeT > FloatLeafManager
Definition: MeshToVolume.h:954
void run(bool threaded=true)
Definition: MeshToVolume.h:1813
tree::Tree4< EdgeData, 5, 4, 3 >::Type TreeType
Definition: MeshToVolume.h:314
tree::ValueAccessor< TreeType > AccessorT
Definition: MeshToVolume.h:2249
GenEdgeData(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition: MeshToVolume.h:2785
math::Vec4< Index32 > Vec4I
Definition: Types.h:91
tree::ValueAccessor< const BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:741
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:447
BoolTreeT & signMaskTree()
Definition: MeshToVolume.h:1122
BoolTreeT & intersectionTree()
Definition: MeshToVolume.h:470
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
static ValueType epsilon()
Definition: MeshToVolume.h:396
Scheme::template ISStencil< GridType >::StencilType Stencil
Definition: MeshToVolume.h:2147
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:453
Efficient multi-threaded replacement of the background values in tree.
math::BIAS_SCHEME< math::FIRST_BIAS > Scheme
Definition: MeshToVolume.h:2146
Type Pow2(Type x)
Return .
Definition: Math.h:498
IntGridT::Ptr indexGridPtr() const
Definition: MeshToVolume.h:251
TBB body object that performs a parallel flood fill.
Definition: MeshToVolume.h:1100
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1439
Definition: MeshToVolume.h:198
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:379
tree::ValueAccessor< TreeType > AccessorT
Definition: MeshToVolume.h:1690
Definition: MeshToVolume.h:2191
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:210
IntersectingVoxelCleaner(FloatTreeT &distTree, IntTreeT &indexTree, BoolTreeT &intersectionTree, BoolLeafManager &leafs)
Definition: MeshToVolume.h:1473
Definition: Operators.h:152
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:1104
TreeType::LeafNodeType LeafNodeT
Definition: MeshToVolume.h:2250
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1567
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2228
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:739
Grid< BoolTreeT > BoolGridT
Definition: MeshToVolume.h:211
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2129
~ShellVoxelCleaner()
Definition: MeshToVolume.h:1574
ExpandNB(BoolLeafManager &leafs, FloatTreeT &distTree, IntTreeT &indexTree, BoolTreeT &maskTree, FloatValueT exteriorBandWidth, FloatValueT interiorBandWidth, FloatValueT voxelSize, const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition: MeshToVolume.h:1767
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1564
Index32 mYPrim
Definition: MeshToVolume.h:311
FloatTreeT & sqrDistTree()
Definition: MeshToVolume.h:468
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2884
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:738
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1254
Grid< IntTreeT > IntGridT
Definition: MeshToVolume.h:209
bool operator==(const EdgeData &rhs) const
Definition: MeshToVolume.h:305
int32_t Int32
Definition: Types.h:61
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:734
void getEdgeData(Accessor &acc, const Coord &ijk, std::vector< Vec3d > &points, std::vector< Index32 > &primitives)
Returns intersection points with corresponding primitive indices for the given ijk voxel...
Definition: MeshToVolume.h:3072
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:334
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:1718
Internal edge data type.
Definition: MeshToVolume.h:286
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:449
void convert(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Threaded method to extract voxel edge data, the closest intersection point and corresponding primitiv...
Definition: MeshToVolume.h:3059
void convertToLevelSet(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, FloatValueT exBandWidth=FloatValueT(LEVEL_SET_HALF_WIDTH), FloatValueT inBandWidth=FloatValueT(LEVEL_SET_HALF_WIDTH))
Mesh to Level Set / Signed Distance Field conversion.
Definition: MeshToVolume.h:2299
void run(bool threaded=true)
Definition: MeshToVolume.h:1463
TBB body object to voxelize a mesh of triangles and/or quads into a collection of VDB grids...
Definition: MeshToVolume.h:443
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:246
~MeshVoxelizer()
Definition: MeshToVolume.h:460
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2162
PropagateSign(BoolLeafManager &, FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:1139
Defined various multi-threaded utility functions for trees.
EdgeData operator-() const
Definition: MeshToVolume.h:302
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1723
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
Definition: LeafManager.h:476
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:210
void join(GenEdgeData &rhs)
Definition: MeshToVolume.h:2822
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:1717
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:485
void combine(FloatTreeT &lhsDist, IntTreeT &lhsIndex, FloatTreeT &rhsDist, IntTreeT &rhsIndex)
Definition: MeshToVolume.h:403
Definition: MeshToVolume.h:2122
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:1561
Definition: Types.h:210
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1565
TBB body object that that finds seed points for the parallel flood fill.
Definition: MeshToVolume.h:949
MeshToVoxelEdgeData()
Definition: MeshToVolume.h:3052
FloatTreeT::LeafNodeType DistLeafT
Definition: MeshToVolume.h:1562
MinOp(LeafManagerType &leafs)
Definition: MeshToVolume.h:2196
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:693
void run(bool threaded=true)
Definition: MeshToVolume.h:1295
ContourTracer(FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:784
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:446
#define OPENVDB_VERSION_NAME
Definition: version.h:43
SqrtAndScaleOp(ValueType voxelSize, bool unsignedDist=false)
Definition: MeshToVolume.h:2033
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2219
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1499
void resetCFL(ValueType cfl)
Definition: MeshToVolume.h:2159
TBB body object that partitions a volume into 2D slices that can be processed in parallel and marks t...
Definition: MeshToVolume.h:735
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:217
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1258
Index32 Index
Definition: Types.h:59
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1721
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:216
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:953
Definition: Exceptions.h:87
EdgeData operator-(const T &) const
Definition: MeshToVolume.h:301
~SignMask()
Definition: MeshToVolume.h:965
~ExpandNB()
Definition: MeshToVolume.h:1737
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:1173
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:740
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1435
Definition: MeshToVolume.h:2060
void changeBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &background, bool threaded=true, size_t grainSize=32)
Replace the background value in all the nodes of a tree.
Definition: ChangeBackground.h:230
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1228
SignMask(const FloatLeafManager &, const FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:991
ScalarToVectorConverter< GridType >::Type::Ptr cpt(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Closest-Point Transform (CPT) from a distance field.
Definition: GridOperators.h:950
TBB body object that removes non-intersecting voxels that where set by rasterizing self-intersecting ...
Definition: MeshToVolume.h:1557
void convertToUnsignedDistanceField(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, FloatValueT exBandWidth)
Mesh to Unsigned Distance Field conversion.
Definition: MeshToVolume.h:2322
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or NULL if no such node exists...
Definition: ValueAccessor.h:383
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
void join(PropagateSign< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:1237
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1720
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:223
~IntersectingVoxelSign()
Definition: MeshToVolume.h:1270
Definition: Exceptions.h:39
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2069
BufferType & getBuffer(size_t leafIdx, size_t bufferIdx) const
Return the leaf or auxiliary buffer for the leaf node at index leafIdx. If bufferIdx is zero...
Definition: LeafManager.h:330
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:580
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1722
tree::LeafManager< FloatTreeT > DistArrayT
Definition: MeshToVolume.h:1563
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1436
tree::LeafManager< TreeType > LeafManagerType
Definition: MeshToVolume.h:2193
Definition: MeshToVolume.h:203
LeafTopologyDiffOp(TreeType &tree)
Definition: MeshToVolume.h:2252
Vec3< double > Vec3d
Definition: Vec3.h:629
float mXDist
Definition: MeshToVolume.h:310
CopyActiveVoxelsOp(TreeType &tree)
Definition: MeshToVolume.h:1692
OffsetOp(ValueType offset)
Definition: MeshToVolume.h:2124
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1256
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1103
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1435
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2255
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:454
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:1105
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2149
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:952
virtual void clear()
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:392
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:450
Definition: Types.h:209
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1437
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:158
LeafNodeT * touchLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, create one, but preserve the values and active states of all voxels.
Definition: ValueAccessor.h:347
void run(bool threaded=true)
Definition: MeshToVolume.h:773
TreeType & tree()
Definition: MeshToVolume.h:2754
void run(bool threaded=true)
Definition: MeshToVolume.h:2811
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2098
Definition: MeshToVolume.h:2089
Coord nearestCoord(const Vec3d &voxelCoord)
Return voxelCoord rounded to the closest integer coordinates.
Definition: Util.h:56
Definition: Exceptions.h:88
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:957
tree::ValueAccessor< const BoolTreeT > BoolConstAccessorT
Definition: MeshToVolume.h:1110
void run(bool threaded=true)
Definition: MeshToVolume.h:1593
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1438
~ContourTracer()
Definition: MeshToVolume.h:744
EdgeData operator+(const T &) const
Definition: MeshToVolume.h:300
IntTreeT & primIndexTree()
Definition: MeshToVolume.h:469
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:658
ShellVoxelCleaner(FloatTreeT &distTree, DistArrayT &leafs, IntTreeT &indexTree, BoolTreeT &intersectionTree)
Definition: MeshToVolume.h:1604
void run(bool threaded=true)
Definition: MeshToVolume.h:370
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1569
MergeBufferOp(LeafManagerType &leafs, size_t bufferIndex=1)
Definition: MeshToVolume.h:2221
Index32 mZPrim
Definition: MeshToVolume.h:311
~PropagateSign()
Definition: MeshToVolume.h:1114
PointTransform(const std::vector< Vec3s > &pointsIn, std::vector< Vec3s > &pointsOut, const math::Transform &xform)
Definition: MeshToVolume.h:362
void pruneLevelSet(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing nodes whose values are all inactive with inactive ...
Definition: Prune.h:369
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1335
Definition: MeshToVolume.h:394
tree::ValueAccessor< const BoolTreeT > BoolConstAccessorT
Definition: MeshToVolume.h:959
void join(ExpandNB< FloatTreeT > &)
Definition: MeshToVolume.h:2019
void resetOffset(ValueType offset)
Definition: MeshToVolume.h:2126
tree::ValueAccessor< const FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:955
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1719
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1106
FloatGridT::TreeType FloatTreeT
Definition: MeshToVolume.h:206
Definition: Tree.h:203
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1100
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:391
std::ostream & operator<<(std::ostream &ostr, const MeshToVoxelEdgeData::EdgeData &rhs)
Definition: MeshToVolume.h:2721
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:452
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1109
VoxelSignOp(ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2062
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1731
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, return NULL.
Definition: Tree.h:1622
GridType::Ptr meshToLevelSet(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, float halfWidth=float(LEVEL_SET_HALF_WIDTH))
Convert a triangle mesh to a level set volume.
Definition: MeshToVolume.h:2647
TBB body object to expand the level set narrow band.
Definition: MeshToVolume.h:1713
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2199
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:451
Definition: Mat4.h:51
uint32_t Index32
Definition: Types.h:57
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1259
MeshVoxelizer(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:520
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
RenormOp(GridType &grid, LeafManagerType &leafs, ValueType voxelSize, ValueType cfl=1.0)
Definition: MeshToVolume.h:2151
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
BoolTreeT & signMaskTree()
Definition: MeshToVolume.h:973
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:1026
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1560
Vec3< float > Vec3s
Definition: Vec3.h:628
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1716
float mZDist
Definition: MeshToVolume.h:310
void join(MeshVoxelizer< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:653
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:220
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2040
Definition: MeshToVolume.h:2144
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:1255
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:448
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:207
void join(SignMask< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:1089
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1108
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1257
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:120
OPENVDB_API Vec3d closestPointOnTriangleToPoint(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p, Vec3d &uvw)
Closest Point on Triangle to Point. Given a triangle abc and a point p, return the point on abc close...
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1260
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1433
void run(bool threaded=true)
Definition: MeshToVolume.h:1017
void run(bool threaded=true)
Definition: MeshToVolume.h:1164
GridType::Ptr meshToLevelSet(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float halfWidth=float(LEVEL_SET_HALF_WIDTH))
Convert a triangle and quad mesh to a level set volume.
Definition: MeshToVolume.h:2675
void setValue(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
Definition: ValueAccessor.h:241
tree::LeafManager< TreeType > LeafManagerType
Definition: MeshToVolume.h:2218
RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:211
FloatGridT::Ptr distGridPtr() const
Returns a narrow-band (signed) distance field / level set grid.
Definition: MeshToVolume.h:247
static ValueType minNarrowBandWidth()
Definition: MeshToVolume.h:397
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1725
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:958
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:232
TrimOp(ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2091
GridType::Ptr meshToSignedDistanceField(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float exBandWidth, float inBandWidth)
Convert a triangle and quad mesh to a signed distance field with an asymmetrical narrow band...
Definition: MeshToVolume.h:2689