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