OpenVDB  2.3.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 = mBBox.min()[2];
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.0) 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.0) {
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 < 26; ++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 template<typename TreeType>
1682 {
1684 
1685  CopyActiveVoxelsOp(TreeType& tree) : mAcc(tree) { }
1686 
1687  template <typename LeafNodeType>
1688  void operator()(LeafNodeType &leaf, size_t) const
1689  {
1690  LeafNodeType* rhsLeaf = const_cast<LeafNodeType*>(mAcc.probeLeaf(leaf.origin()));
1691  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
1692  for (; iter; ++iter) {
1693  rhsLeaf->setValueOnly(iter.pos(), iter.getValue());
1694  }
1695  }
1696 
1697 private:
1698  AccessorT mAcc;
1699 };
1700 
1701 
1702 // ExpandNB
1705 template<typename FloatTreeT>
1707 {
1708 public:
1709  typedef typename FloatTreeT::ValueType FloatValueT;
1710  typedef typename FloatTreeT::LeafNodeType FloatLeafT;
1712  typedef typename FloatTreeT::template ValueConverter<Int32>::Type IntTreeT;
1713  typedef typename IntTreeT::LeafNodeType IntLeafT;
1715  typedef typename FloatTreeT::template ValueConverter<bool>::Type BoolTreeT;
1716  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1719 
1720  ExpandNB(BoolLeafManager& leafs,
1721  FloatTreeT& distTree, IntTreeT& indexTree, BoolTreeT& maskTree,
1722  FloatValueT exteriorBandWidth, FloatValueT interiorBandWidth, FloatValueT voxelSize,
1723  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
1724 
1725  void run(bool threaded = true);
1726 
1727  void operator()(const tbb::blocked_range<size_t>&);
1728  void join(ExpandNB<FloatTreeT>&);
1729  ExpandNB(const ExpandNB<FloatTreeT>&, tbb::split);
1731 
1732 private:
1733  void operator=(const ExpandNB<FloatTreeT>&) {}
1734 
1735  double evalVoxelDist(const Coord&, FloatAccessorT&, IntAccessorT&,
1736  BoolAccessorT&, std::vector<Int32>&, Int32&) const;
1737 
1738  double evalVoxelDist(const Coord&, FloatLeafT&, IntLeafT&,
1739  BoolLeafT&, std::vector<Int32>&, Int32&) const;
1740 
1741  double closestPrimDist(const Coord&, std::vector<Int32>&, Int32&) const;
1742 
1743  BoolLeafManager& mMaskLeafs;
1744 
1745  FloatTreeT& mDistTree;
1746  IntTreeT& mIndexTree;
1747  BoolTreeT& mMaskTree;
1748 
1749  const FloatValueT mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
1750  const std::vector<Vec3s>& mPointList;
1751  const std::vector<Vec4I>& mPolygonList;
1752 
1753  FloatTreeT mNewDistTree;
1754  IntTreeT mNewIndexTree;
1755  BoolTreeT mNewMaskTree;
1756 };
1757 
1758 
1759 template<typename FloatTreeT>
1761  BoolLeafManager& leafs,
1762  FloatTreeT& distTree,
1763  IntTreeT& indexTree,
1764  BoolTreeT& maskTree,
1765  FloatValueT exteriorBandWidth,
1766  FloatValueT interiorBandWidth,
1767  FloatValueT voxelSize,
1768  const std::vector<Vec3s>& pointList,
1769  const std::vector<Vec4I>& polygonList)
1770  : mMaskLeafs(leafs)
1771  , mDistTree(distTree)
1772  , mIndexTree(indexTree)
1773  , mMaskTree(maskTree)
1774  , mExteriorBandWidth(exteriorBandWidth)
1775  , mInteriorBandWidth(interiorBandWidth)
1776  , mVoxelSize(voxelSize)
1777  , mPointList(pointList)
1778  , mPolygonList(polygonList)
1779  , mNewDistTree(std::numeric_limits<FloatValueT>::max())
1780  , mNewIndexTree(Int32(util::INVALID_IDX))
1781  , mNewMaskTree(false)
1782 {
1783 }
1784 
1785 
1786 template<typename FloatTreeT>
1788  : mMaskLeafs(rhs.mMaskLeafs)
1789  , mDistTree(rhs.mDistTree)
1790  , mIndexTree(rhs.mIndexTree)
1791  , mMaskTree(rhs.mMaskTree)
1792  , mExteriorBandWidth(rhs.mExteriorBandWidth)
1793  , mInteriorBandWidth(rhs.mInteriorBandWidth)
1794  , mVoxelSize(rhs.mVoxelSize)
1795  , mPointList(rhs.mPointList)
1796  , mPolygonList(rhs.mPolygonList)
1797  , mNewDistTree(std::numeric_limits<FloatValueT>::max())
1798  , mNewIndexTree(Int32(util::INVALID_IDX))
1799  , mNewMaskTree(false)
1800 {
1801 }
1802 
1803 
1804 template<typename FloatTreeT>
1805 void
1807 {
1808  if (threaded) tbb::parallel_reduce(mMaskLeafs.getRange(), *this);
1809  else (*this)(mMaskLeafs.getRange());
1810 
1811  // Copy only the active voxels (tree::merge does branch stealing
1812  // which also moves indicative values).
1813  mDistTree.topologyUnion(mNewDistTree);
1814  tree::LeafManager<FloatTreeT> leafs(mNewDistTree);
1815  leafs.foreach(CopyActiveVoxelsOp<FloatTreeT>(mDistTree));
1816 
1817  mIndexTree.merge(mNewIndexTree);
1818 
1819  mMaskTree.clear();
1820  mMaskTree.merge(mNewMaskTree);
1821 
1822 }
1823 
1824 
1825 template<typename FloatTreeT>
1826 void
1827 ExpandNB<FloatTreeT>::operator()(const tbb::blocked_range<size_t>& range)
1828 {
1829  Coord ijk;
1830  Int32 closestPrim = 0;
1831  Index pos = 0;
1832  FloatValueT distance;
1833  bool inside;
1834 
1835  FloatAccessorT newDistAcc(mNewDistTree);
1836  IntAccessorT newIndexAcc(mNewIndexTree);
1837  BoolAccessorT newMaskAcc(mNewMaskTree);
1838 
1839  FloatAccessorT distAcc(mDistTree);
1840  IntAccessorT indexAcc(mIndexTree);
1841  BoolAccessorT maskAcc(mMaskTree);
1842 
1843  CoordBBox bbox;
1844  std::vector<Int32> primitives(18);
1845 
1846  for (size_t n = range.begin(); n < range.end(); ++n) {
1847 
1848  BoolLeafT& maskLeaf = mMaskLeafs.leaf(n);
1849 
1850  if (maskLeaf.isEmpty()) continue;
1851 
1852  ijk = maskLeaf.origin();
1853 
1854  FloatLeafT* distLeafPt = distAcc.probeLeaf(ijk);
1855 
1856  if (!distLeafPt) {
1857  distLeafPt = new FloatLeafT(ijk, distAcc.getValue(ijk));
1858  newDistAcc.addLeaf(distLeafPt);
1859  }
1860 
1861  IntLeafT* indexLeafPt = indexAcc.probeLeaf(ijk);
1862  if (!indexLeafPt) indexLeafPt = newIndexAcc.touchLeaf(ijk);
1863 
1864  bbox = maskLeaf.getNodeBoundingBox();
1865  bbox.expand(-1);
1866 
1867  typename BoolLeafT::ValueOnIter iter = maskLeaf.beginValueOn();
1868  for (; iter; ++iter) {
1869 
1870  ijk = iter.getCoord();
1871 
1872  if (bbox.isInside(ijk)) {
1873  distance = evalVoxelDist(ijk, *distLeafPt, *indexLeafPt, maskLeaf,
1874  primitives, closestPrim);
1875  } else {
1876  distance = evalVoxelDist(ijk, distAcc, indexAcc, maskAcc,
1877  primitives, closestPrim);
1878  }
1879 
1880  pos = iter.pos();
1881 
1882  inside = distLeafPt->getValue(pos) < FloatValueT(0.0);
1883 
1884  if (!inside && distance < mExteriorBandWidth) {
1885  distLeafPt->setValueOn(pos, distance);
1886  indexLeafPt->setValueOn(pos, closestPrim);
1887  } else if (inside && distance < mInteriorBandWidth) {
1888  distLeafPt->setValueOn(pos, -distance);
1889  indexLeafPt->setValueOn(pos, closestPrim);
1890  } else {
1891  continue;
1892  }
1893 
1894  for (Int32 i = 0; i < 6; ++i) {
1895  newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[i]);
1896  }
1897  }
1898  }
1899 }
1900 
1901 
1902 template<typename FloatTreeT>
1903 double
1905  const Coord& ijk,
1906  FloatAccessorT& distAcc,
1907  IntAccessorT& indexAcc,
1908  BoolAccessorT& maskAcc,
1909  std::vector<Int32>& prims,
1910  Int32& closestPrim) const
1911 {
1912  FloatValueT tmpDist, minDist = std::numeric_limits<FloatValueT>::max();
1913  prims.clear();
1914 
1915  // Collect primitive indices from active neighbors and min distance.
1916  Coord n_ijk;
1917  for (Int32 n = 0; n < 18; ++n) {
1918  n_ijk = ijk + util::COORD_OFFSETS[n];
1919  if (!maskAcc.isValueOn(n_ijk) && distAcc.probeValue(n_ijk, tmpDist)) {
1920  prims.push_back(indexAcc.getValue(n_ijk));
1921  tmpDist = std::abs(tmpDist);
1922  if (tmpDist < minDist) minDist = tmpDist;
1923  }
1924  }
1925 
1926  // Calc. this voxels distance to the closest primitive.
1927  tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1928 
1929  // Forces the gradient to be monotonic for non-manifold
1930  // polygonal models with self-intersections.
1931  return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1932 }
1933 
1934 
1935 // Leaf specialized version.
1936 template<typename FloatTreeT>
1937 double
1938 ExpandNB<FloatTreeT>::evalVoxelDist(
1939  const Coord& ijk,
1940  FloatLeafT& distLeaf,
1941  IntLeafT& indexLeaf,
1942  BoolLeafT& maskLeaf,
1943  std::vector<Int32>& prims,
1944  Int32& closestPrim) const
1945 {
1946  FloatValueT tmpDist, minDist = std::numeric_limits<FloatValueT>::max();
1947  prims.clear();
1948 
1949  Index pos;
1950  for (Int32 n = 0; n < 18; ++n) {
1951  pos = FloatLeafT::coordToOffset(ijk + util::COORD_OFFSETS[n]);
1952  if (!maskLeaf.isValueOn(pos) && distLeaf.probeValue(pos, tmpDist)) {
1953  prims.push_back(indexLeaf.getValue(pos));
1954  tmpDist = std::abs(tmpDist);
1955  if (tmpDist < minDist) minDist = tmpDist;
1956  }
1957  }
1958 
1959  tmpDist = FloatValueT(closestPrimDist(ijk, prims, closestPrim));
1960  return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
1961 }
1962 
1963 
1964 template<typename FloatTreeT>
1965 double
1966 ExpandNB<FloatTreeT>::closestPrimDist(const Coord& ijk,
1967  std::vector<Int32>& prims, Int32& closestPrim) const
1968 {
1969  std::sort(prims.begin(), prims.end());
1970 
1971  Int32 lastPrim = -1;
1972  Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
1973  double primDist, tmpDist, dist = std::numeric_limits<double>::max();
1974 
1975  for (size_t n = 0, N = prims.size(); n < N; ++n) {
1976  if (prims[n] == lastPrim) continue;
1977 
1978  lastPrim = prims[n];
1979 
1980  const Vec4I& verts = mPolygonList[lastPrim];
1981 
1982  // Evaluate first triangle
1983  const Vec3d a(mPointList[verts[0]]);
1984  const Vec3d b(mPointList[verts[1]]);
1985  const Vec3d c(mPointList[verts[2]]);
1986 
1987  primDist = (voxelCenter -
1988  closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
1989 
1990  // Split-up quad into a second triangle and calac distance.
1991  if (util::INVALID_IDX != verts[3]) {
1992  const Vec3d d(mPointList[verts[3]]);
1993 
1994  tmpDist = (voxelCenter -
1995  closestPointOnTriangleToPoint(a, d, c, voxelCenter, uvw)).lengthSqr();
1996 
1997  if (tmpDist < primDist) primDist = tmpDist;
1998  }
1999 
2000  if (primDist < dist) {
2001  dist = primDist;
2002  closestPrim = lastPrim;
2003  }
2004  }
2005 
2006  return std::sqrt(dist) * double(mVoxelSize);
2007 }
2008 
2009 
2010 template<typename FloatTreeT>
2011 void
2013 {
2014  mNewDistTree.merge(rhs.mNewDistTree);
2015  mNewIndexTree.merge(rhs.mNewIndexTree);
2016  mNewMaskTree.merge(rhs.mNewMaskTree);
2017 }
2018 
2019 
2021 
2022 
2023 template<typename ValueType>
2025 {
2026  SqrtAndScaleOp(ValueType voxelSize, bool unsignedDist = false)
2027  : mVoxelSize(voxelSize)
2028  , mUnsigned(unsignedDist)
2029  {
2030  }
2031 
2032  template <typename LeafNodeType>
2033  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2034  {
2035  ValueType w[2];
2036  w[0] = mVoxelSize;
2037  w[1] = -mVoxelSize;
2038 
2039  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2040  for (; iter; ++iter) {
2041  ValueType& val = const_cast<ValueType&>(iter.getValue());
2042  val = w[!mUnsigned && int(val < ValueType(0.0))] * std::sqrt(std::abs(val));
2043  }
2044  }
2045 
2046 private:
2047  ValueType mVoxelSize;
2048  const bool mUnsigned;
2049 };
2050 
2051 
2052 template<typename ValueType>
2054 {
2055  VoxelSignOp(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  ValueType bgValues[2];
2065  bgValues[0] = mExBandWidth;
2066  bgValues[1] = -mInBandWidth;
2067 
2068  typename LeafNodeType::ValueOffIter iter = leaf.beginValueOff();
2069 
2070  for (; iter; ++iter) {
2071  ValueType& val = const_cast<ValueType&>(iter.getValue());
2072  val = bgValues[int(val < ValueType(0.0))];
2073  }
2074  }
2075 
2076 private:
2077  ValueType mExBandWidth, mInBandWidth;
2078 };
2079 
2080 
2081 template<typename ValueType>
2082 struct TrimOp
2083 {
2084  TrimOp(ValueType exBandWidth, ValueType inBandWidth)
2085  : mExBandWidth(exBandWidth)
2086  , mInBandWidth(inBandWidth)
2087  {
2088  }
2089 
2090  template <typename LeafNodeType>
2091  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2092  {
2093  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2094 
2095  for (; iter; ++iter) {
2096  ValueType& val = const_cast<ValueType&>(iter.getValue());
2097  const bool inside = val < ValueType(0.0);
2098 
2099  if (inside && !(val > -mInBandWidth)) {
2100  val = -mInBandWidth;
2101  iter.setValueOff();
2102  } else if (!inside && !(val < mExBandWidth)) {
2103  val = mExBandWidth;
2104  iter.setValueOff();
2105  }
2106  }
2107  }
2108 
2109 private:
2110  ValueType mExBandWidth, mInBandWidth;
2111 };
2112 
2113 
2114 template<typename ValueType>
2115 struct OffsetOp
2116 {
2117  OffsetOp(ValueType offset): mOffset(offset) {}
2118 
2119  void resetOffset(ValueType offset) { mOffset = offset; }
2120 
2121  template <typename LeafNodeType>
2122  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2123  {
2124  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2125  for (; iter; ++iter) {
2126  ValueType& val = const_cast<ValueType&>(iter.getValue());
2127  val += mOffset;
2128  }
2129  }
2130 
2131 private:
2132  ValueType mOffset;
2133 };
2134 
2135 
2136 template<typename GridType, typename ValueType>
2137 struct RenormOp
2138 {
2140  typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
2143 
2144  RenormOp(GridType& grid, LeafManagerType& leafs, ValueType voxelSize, ValueType cfl = 1.0)
2145  : mGrid(grid)
2146  , mLeafs(leafs)
2147  , mVoxelSize(voxelSize)
2148  , mCFL(cfl)
2149  {
2150  }
2151 
2152  void resetCFL(ValueType cfl) { mCFL = cfl; }
2153 
2154  template <typename LeafNodeType>
2155  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2156  {
2157  const ValueType dt = mCFL * mVoxelSize, one(1.0), invDx = one / mVoxelSize;
2158  Stencil stencil(mGrid);
2159  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
2160 
2161  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2162  for (; iter; ++iter) {
2163  stencil.moveTo(iter);
2164 
2165  const ValueType normSqGradPhi =
2167 
2168  const ValueType phi0 = iter.getValue();
2169  const ValueType diff = math::Sqrt(normSqGradPhi) * invDx - one;
2170  const ValueType S = phi0 / (math::Sqrt(math::Pow2(phi0) + normSqGradPhi));
2171 
2172  buffer.setValue(iter.pos(), phi0 - dt * S * diff);
2173  }
2174  }
2175 
2176 private:
2177  GridType& mGrid;
2178  LeafManagerType& mLeafs;
2179  ValueType mVoxelSize, mCFL;
2180 };
2181 
2182 
2183 template<typename TreeType, typename ValueType>
2184 struct MinOp
2185 {
2188 
2189  MinOp(LeafManagerType& leafs): mLeafs(leafs) {}
2190 
2191  template <typename LeafNodeType>
2192  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2193  {
2194  BufferType& buffer = mLeafs.getBuffer(leafIndex, 1);
2195  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2196 
2197  for (; iter; ++iter) {
2198  ValueType& val = const_cast<ValueType&>(iter.getValue());
2199  val = std::min(val, buffer.getValue(iter.pos()));
2200  }
2201  }
2202 
2203 private:
2204  LeafManagerType& mLeafs;
2205 };
2206 
2207 
2208 template<typename TreeType, typename ValueType>
2210 {
2213 
2214  MergeBufferOp(LeafManagerType& leafs, size_t bufferIndex = 1)
2215  : mLeafs(leafs)
2216  , mBufferIndex(bufferIndex)
2217  {
2218  }
2219 
2220  template <typename LeafNodeType>
2221  void operator()(LeafNodeType &leaf, size_t leafIndex) const
2222  {
2223  BufferType& buffer = mLeafs.getBuffer(leafIndex, mBufferIndex);
2224  typename LeafNodeType::ValueOnIter iter = leaf.beginValueOn();
2225  Index offset;
2226 
2227  for (; iter; ++iter) {
2228  offset = iter.pos();
2229  leaf.setValueOnly(offset, buffer.getValue(offset));
2230  }
2231  }
2232 
2233 private:
2234  LeafManagerType& mLeafs;
2235  const size_t mBufferIndex;
2236 };
2237 
2238 
2239 template<typename TreeType>
2241 {
2243  typedef typename TreeType::LeafNodeType LeafNodeT;
2244 
2245  LeafTopologyDiffOp(TreeType& tree) : mAcc(tree) { }
2246 
2247  template <typename LeafNodeType>
2248  void operator()(LeafNodeType &leaf, size_t) const
2249  {
2250  const LeafNodeT* rhsLeaf = mAcc.probeConstLeaf(leaf.origin());
2251  if (rhsLeaf) leaf.topologyDifference(*rhsLeaf, false);
2252  }
2253 
2254 private:
2255  AccessorT mAcc;
2256 };
2257 
2258 
2259 template<typename ValueType>
2260 struct SDFPrune
2261 {
2262  SDFPrune(const ValueType& out, const ValueType& in)
2263  : outside(out)
2264  , inside(in)
2265  {
2266  }
2267 
2268  template <typename ChildType>
2269  bool operator()(ChildType& child)
2270  {
2271  child.pruneOp(*this);
2272  if (!child.isInactive()) return false;
2273  value = math::isNegative(child.getFirstValue()) ? inside : outside;
2274  return true;
2275  }
2276 
2277  static const bool state = false;
2278  const ValueType outside, inside;
2279  ValueType value;
2280 };
2281 
2282 
2283 } // internal namespace
2284 
2285 
2287 
2288 
2289 // MeshToVolume
2290 
2291 template<typename FloatGridT, typename InterruptT>
2293  openvdb::math::Transform::Ptr& transform, int conversionFlags,
2294  InterruptT *interrupter, int signSweeps)
2295  : mTransform(transform)
2296  , mConversionFlags(conversionFlags)
2297  , mSignSweeps(signSweeps)
2298  , mInterrupter(interrupter)
2299 {
2300  clear();
2301  mSignSweeps = std::min(mSignSweeps, 1);
2302 }
2303 
2304 
2305 template<typename FloatGridT, typename InterruptT>
2306 void
2308 {
2309  mDistGrid = FloatGridT::create(std::numeric_limits<FloatValueT>::max());
2310  mIndexGrid = IntGridT::create(Int32(util::INVALID_IDX));
2311  mIntersectingVoxelsGrid = BoolGridT::create(false);
2312 }
2313 
2314 
2315 template<typename FloatGridT, typename InterruptT>
2316 inline void
2318  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2319  FloatValueT exBandWidth, FloatValueT inBandWidth)
2320 {
2321  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
2322  exBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), exBandWidth);
2323  inBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), inBandWidth);
2324  const FloatValueT vs = mTransform->voxelSize()[0];
2325  doConvert(pointList, polygonList, vs * exBandWidth, vs * inBandWidth);
2326  mDistGrid->setGridClass(GRID_LEVEL_SET);
2327 }
2328 
2329 
2330 template<typename FloatGridT, typename InterruptT>
2331 inline void
2333  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2334  FloatValueT exBandWidth)
2335 {
2336  // The narrow band width is exclusive, the shortest valid distance has to be > 1 voxel
2337  exBandWidth = std::max(internal::Tolerance<FloatValueT>::minNarrowBandWidth(), exBandWidth);
2338  const FloatValueT vs = mTransform->voxelSize()[0];
2339  doConvert(pointList, polygonList, vs * exBandWidth, 0.0, true);
2340  mDistGrid->setGridClass(GRID_UNKNOWN);
2341 }
2342 
2343 
2344 template<typename FloatGridT, typename InterruptT>
2345 void
2347  const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList,
2348  FloatValueT exBandWidth, FloatValueT inBandWidth, bool unsignedDistField)
2349 {
2350  mDistGrid->setTransform(mTransform);
2351  mIndexGrid->setTransform(mTransform);
2352  const bool rawData = OUTPUT_RAW_DATA & mConversionFlags;
2353 
2354  // The progress estimates given to the interrupter are based on the
2355  // observed average time for each stage and therefore not alway
2356  // accurate. The goal is to give some progression feedback to the user.
2357 
2358  if (wasInterrupted(1)) return;
2359 
2360  // Voxelize mesh
2361  {
2363  voxelizer(pointList, polygonList, mInterrupter);
2364 
2365  voxelizer.run();
2366 
2367  if (wasInterrupted(18)) return;
2368 
2369  mDistGrid->tree().merge(voxelizer.sqrDistTree());
2370  mIndexGrid->tree().merge(voxelizer.primIndexTree());
2371  mIntersectingVoxelsGrid->tree().merge(voxelizer.intersectionTree());
2372  }
2373 
2374  if (!unsignedDistField) {
2375  // Determine the inside/outside state for the narrow band of voxels.
2376  {
2377  // Slices up the volume and label the exterior contour of each slice in parallel.
2378  internal::ContourTracer<FloatTreeT, InterruptT> trace(
2379  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2380  for (int i = 0; i < mSignSweeps; ++i) {
2381 
2382  if (wasInterrupted(19)) return;
2383 
2384  trace.run();
2385 
2386  if (wasInterrupted(24)) return;
2387 
2388  // Propagate sign information between the slices.
2389  BoolTreeT signMaskTree(false);
2390  {
2391  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2392  internal::SignMask<FloatTreeT, InterruptT> signMaskOp(leafs,
2393  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2394  signMaskOp.run();
2395  signMaskTree.merge(signMaskOp.signMaskTree());
2396  }
2397 
2398  if (wasInterrupted(25)) return;
2399 
2400  while (true) {
2401  tree::LeafManager<BoolTreeT> leafs(signMaskTree);
2402  if(leafs.leafCount() == 0) break;
2403 
2404  internal::PropagateSign<FloatTreeT, InterruptT> sign(leafs,
2405  mDistGrid->tree(), mIntersectingVoxelsGrid->tree(), mInterrupter);
2406 
2407  sign.run();
2408 
2409  signMaskTree.clear();
2410  signMaskTree.merge(sign.signMaskTree());
2411  }
2412  }
2413  }
2414 
2415 
2416  if (wasInterrupted(28)) return;
2417  {
2418  tree::LeafManager<BoolTreeT> leafs(mIntersectingVoxelsGrid->tree());
2419 
2420  // Determine the sign of the mesh intersecting voxels.
2421  internal::IntersectingVoxelSign<FloatTreeT> sign(pointList, polygonList,
2422  mDistGrid->tree(), mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2423 
2424  sign.run();
2425 
2426  if (wasInterrupted(34)) return;
2427 
2428  // Remove mesh intersecting voxels that where set by rasterizing
2429  // self-intersecting portions of the mesh.
2430  internal::IntersectingVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2431  mIndexGrid->tree(), mIntersectingVoxelsGrid->tree(), leafs);
2432  cleaner.run();
2433  }
2434 
2435  // Remove shell voxels that where set by rasterizing
2436  // self-intersecting portions of the mesh.
2437  {
2438  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2439 
2440  internal::ShellVoxelCleaner<FloatTreeT> cleaner(mDistGrid->tree(),
2441  leafs, mIndexGrid->tree(), mIntersectingVoxelsGrid->tree());
2442 
2443  cleaner.run();
2444  }
2445 
2446  if (wasInterrupted(38)) return;
2447 
2448  } else { // if unsigned dist. field
2449  inBandWidth = FloatValueT(0.0);
2450  }
2451 
2452  if (mDistGrid->activeVoxelCount() == 0) return;
2453 
2454  mIntersectingVoxelsGrid->clear();
2455  const FloatValueT voxelSize(mTransform->voxelSize()[0]);
2456 
2457  { // Transform values (world space scaling etc.)
2458  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2459  leafs.foreach(internal::SqrtAndScaleOp<FloatValueT>(voxelSize, unsignedDistField));
2460  }
2461 
2462  if (wasInterrupted(40)) return;
2463 
2464  if (!unsignedDistField) { // Propagate sign information to inactive values.
2465  mDistGrid->tree().root().setBackground(exBandWidth, /*updateChildNodes=*/false);
2466  mDistGrid->tree().signedFloodFill(exBandWidth, -inBandWidth);
2467  }
2468 
2469  if (wasInterrupted(46)) return;
2470 
2471  // Narrow-band dilation
2472  const FloatValueT minWidth = voxelSize * 2.0;
2473  if (inBandWidth > minWidth || exBandWidth > minWidth) {
2474 
2475  // Create the initial voxel mask.
2476  BoolTreeT maskTree(false);
2477  maskTree.topologyUnion(mDistGrid->tree());
2478 
2479  if (wasInterrupted(48)) return;
2480 
2481  internal::LeafTopologyDiffOp<FloatTreeT> diffOp(mDistGrid->tree());
2482  openvdb::tools::dilateVoxels(maskTree);
2483 
2484  unsigned maxIterations = std::numeric_limits<unsigned>::max();
2485  float progress = 48, step = 0.0;
2486  // progress estimation..
2487  double estimated =
2488  2.0 * std::ceil((std::max(inBandWidth, exBandWidth) - minWidth) / voxelSize);
2489  if (estimated < double(maxIterations)) {
2490  maxIterations = unsigned(estimated);
2491  step = 42.0 / float(maxIterations);
2492  }
2493 
2494  unsigned count = 0;
2495  while (true) {
2496 
2497  if (wasInterrupted(int(progress))) return;
2498 
2499  tree::LeafManager<BoolTreeT> leafs(maskTree);
2500 
2501  if (leafs.leafCount() == 0) break;
2502 
2503  leafs.foreach(diffOp);
2504 
2505  internal::ExpandNB<FloatTreeT> expand(
2506  leafs, mDistGrid->tree(), mIndexGrid->tree(), maskTree,
2507  exBandWidth, inBandWidth, voxelSize, pointList, polygonList);
2508 
2509  expand.run();
2510 
2511  if ((++count) >= maxIterations) break;
2512  progress += step;
2513  }
2514  }
2515 
2516  if (!bool(GENERATE_PRIM_INDEX_GRID & mConversionFlags)) mIndexGrid->clear();
2517 
2518  if (wasInterrupted(80)) return;
2519 
2520  // Renormalize distances to smooth out bumps caused by self-intersecting
2521  // and overlapping portions of the mesh and renormalize the level set.
2522  if (!unsignedDistField && !rawData) {
2523 
2524  mDistGrid->tree().pruneLevelSet();
2525  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree(), 1);
2526 
2527  const FloatValueT offset = 0.8 * voxelSize;
2528  if (wasInterrupted(82)) return;
2529 
2530  internal::OffsetOp<FloatValueT> offsetOp(-offset);
2531 
2532  leafs.foreach(offsetOp);
2533 
2534  if (wasInterrupted(84)) return;
2535 
2536  leafs.foreach(internal::RenormOp<FloatGridT, FloatValueT>(*mDistGrid, leafs, voxelSize));
2537 
2538  leafs.foreach(internal::MinOp<FloatTreeT, FloatValueT>(leafs));
2539 
2540  if (wasInterrupted(95)) return;
2541 
2542  offsetOp.resetOffset(offset - internal::Tolerance<FloatValueT>::epsilon());
2543  leafs.foreach(offsetOp);
2544  }
2545 
2546  if (wasInterrupted(98)) return;
2547 
2548  const FloatValueT minTrimWidth = voxelSize * 4.0;
2549  if (inBandWidth < minTrimWidth || exBandWidth < minTrimWidth) {
2550 
2551  // If the narrow band was not expanded, we might need to trim off
2552  // some of the active voxels in order to respect the narrow band limits.
2553  // (The mesh voxelization step generates some extra 'shell' voxels)
2554 
2555  tree::LeafManager<FloatTreeT> leafs(mDistGrid->tree());
2556  leafs.foreach(internal::TrimOp<FloatValueT>(
2557  exBandWidth, unsignedDistField ? exBandWidth : inBandWidth));
2558 
2559  internal::SDFPrune<FloatValueT> sdfPrune(exBandWidth, -inBandWidth);
2560  mDistGrid->tree().pruneOp(sdfPrune);
2561  }
2562 }
2563 
2564 
2566 
2567 
2569 template<typename GridType>
2570 inline typename boost::enable_if<boost::is_floating_point<typename GridType::ValueType>,
2571 typename GridType::Ptr>::type
2573  const openvdb::math::Transform& xform,
2574  const std::vector<Vec3s>& points,
2575  const std::vector<Vec3I>& triangles,
2576  const std::vector<Vec4I>& quads,
2577  float exBandWidth,
2578  float inBandWidth,
2579  bool unsignedDistanceField = false)
2580 {
2581  std::vector<Vec3s> indexSpacePoints(points.size());
2582 
2583  { // Copy and transform (required for MeshToVolume) points to grid space.
2584  internal::PointTransform ptnXForm(points, indexSpacePoints, xform);
2585  ptnXForm.run();
2586  }
2587 
2588  // Copy primitives
2589  std::vector<Vec4I> primitives(triangles.size() + quads.size());
2590 
2591  for (size_t n = 0, N = triangles.size(); n < N; ++n) {
2592  Vec4I& prim = primitives[n];
2593  const Vec3I& triangle = triangles[n];
2594  prim[0] = triangle[0];
2595  prim[1] = triangle[1];
2596  prim[2] = triangle[2];
2597  prim[3] = util::INVALID_IDX;
2598  }
2599 
2600  for (size_t n = 0, N = quads.size(); n < N; ++n) {
2601  primitives[n + triangles.size()] = quads[n];
2602  }
2603 
2604  typename GridType::ValueType exWidth(exBandWidth);
2605  typename GridType::ValueType inWidth(inBandWidth);
2606 
2607 
2608  math::Transform::Ptr transform = xform.copy();
2609  MeshToVolume<GridType> vol(transform);
2610 
2611  if (!unsignedDistanceField) {
2612  vol.convertToLevelSet(indexSpacePoints, primitives, exWidth, inWidth);
2613  } else {
2614  vol.convertToUnsignedDistanceField(indexSpacePoints, primitives, exWidth);
2615  }
2616 
2617  return vol.distGridPtr();
2618 }
2619 
2620 
2623 template<typename GridType>
2624 inline typename boost::disable_if<boost::is_floating_point<typename GridType::ValueType>,
2625 typename GridType::Ptr>::type
2627  const math::Transform& /*xform*/,
2628  const std::vector<Vec3s>& /*points*/,
2629  const std::vector<Vec3I>& /*triangles*/,
2630  const std::vector<Vec4I>& /*quads*/,
2631  float /*exBandWidth*/,
2632  float /*inBandWidth*/,
2633  bool unsignedDistanceField = false)
2634 {
2636  "mesh to volume conversion is supported only for scalar, floating-point grids");
2637 }
2638 
2639 
2641 
2642 
2643 template<typename GridType>
2644 inline typename GridType::Ptr
2646  const openvdb::math::Transform& xform,
2647  const std::vector<Vec3s>& points,
2648  const std::vector<Vec3I>& triangles,
2649  float halfWidth)
2650 {
2651  std::vector<Vec4I> quads(0);
2652  return doMeshConversion<GridType>(xform, points, triangles, quads,
2653  halfWidth, halfWidth);
2654 }
2655 
2656 
2657 template<typename GridType>
2658 inline typename GridType::Ptr
2660  const openvdb::math::Transform& xform,
2661  const std::vector<Vec3s>& points,
2662  const std::vector<Vec4I>& quads,
2663  float halfWidth)
2664 {
2665  std::vector<Vec3I> triangles(0);
2666  return doMeshConversion<GridType>(xform, points, triangles, quads,
2667  halfWidth, halfWidth);
2668 }
2669 
2670 
2671 template<typename GridType>
2672 inline typename GridType::Ptr
2674  const openvdb::math::Transform& xform,
2675  const std::vector<Vec3s>& points,
2676  const std::vector<Vec3I>& triangles,
2677  const std::vector<Vec4I>& quads,
2678  float halfWidth)
2679 {
2680  return doMeshConversion<GridType>(xform, points, triangles, quads,
2681  halfWidth, halfWidth);
2682 }
2683 
2684 
2685 template<typename GridType>
2686 inline typename GridType::Ptr
2688  const openvdb::math::Transform& xform,
2689  const std::vector<Vec3s>& points,
2690  const std::vector<Vec3I>& triangles,
2691  const std::vector<Vec4I>& quads,
2692  float exBandWidth,
2693  float inBandWidth)
2694 {
2695  return doMeshConversion<GridType>(xform, points, triangles,
2696  quads, exBandWidth, inBandWidth);
2697 }
2698 
2699 
2700 template<typename GridType>
2701 inline typename GridType::Ptr
2703  const openvdb::math::Transform& xform,
2704  const std::vector<Vec3s>& points,
2705  const std::vector<Vec3I>& triangles,
2706  const std::vector<Vec4I>& quads,
2707  float bandWidth)
2708 {
2709  return doMeshConversion<GridType>(xform, points, triangles, quads,
2710  bandWidth, bandWidth, true);
2711 }
2712 
2713 
2715 
2716 
2717 // Required by several of the tree nodes
2718 inline std::ostream&
2719 operator<<(std::ostream& ostr, const MeshToVoxelEdgeData::EdgeData& rhs)
2720 {
2721  ostr << "{[ " << rhs.mXPrim << ", " << rhs.mXDist << "]";
2722  ostr << " [ " << rhs.mYPrim << ", " << rhs.mYDist << "]";
2723  ostr << " [ " << rhs.mZPrim << ", " << rhs.mZDist << "]}";
2724  return ostr;
2725 }
2726 
2727 // Required by math::Abs
2728 inline MeshToVoxelEdgeData::EdgeData
2730 {
2731  return x;
2732 }
2733 
2734 
2736 
2737 
2739 {
2740 public:
2741 
2742  GenEdgeData(
2743  const std::vector<Vec3s>& pointList,
2744  const std::vector<Vec4I>& polygonList);
2745 
2746  void run(bool threaded = true);
2747 
2748  GenEdgeData(GenEdgeData& rhs, tbb::split);
2749  inline void operator() (const tbb::blocked_range<size_t> &range);
2750  inline void join(GenEdgeData& rhs);
2751 
2752  inline TreeType& tree() { return mTree; }
2753 
2754 private:
2755  void operator=(const GenEdgeData&) {}
2756 
2757  struct Primitive { Vec3d a, b, c, d; Int32 index; };
2758 
2759  template<bool IsQuad>
2760  inline void voxelize(const Primitive&);
2761 
2762  template<bool IsQuad>
2763  inline bool evalPrimitive(const Coord&, const Primitive&);
2764 
2765  inline bool rayTriangleIntersection( const Vec3d& origin, const Vec3d& dir,
2766  const Vec3d& a, const Vec3d& b, const Vec3d& c, double& t);
2767 
2768 
2769  TreeType mTree;
2770  Accessor mAccessor;
2771 
2772  const std::vector<Vec3s>& mPointList;
2773  const std::vector<Vec4I>& mPolygonList;
2774 
2775  // Used internally for acceleration
2776  typedef TreeType::ValueConverter<Int32>::Type IntTreeT;
2777  IntTreeT mLastPrimTree;
2778  tree::ValueAccessor<IntTreeT> mLastPrimAccessor;
2779 }; // class MeshToVoxelEdgeData::GenEdgeData
2780 
2781 
2782 inline
2784  const std::vector<Vec3s>& pointList,
2785  const std::vector<Vec4I>& polygonList)
2786  : mTree(EdgeData())
2787  , mAccessor(mTree)
2788  , mPointList(pointList)
2789  , mPolygonList(polygonList)
2790  , mLastPrimTree(Int32(util::INVALID_IDX))
2791  , mLastPrimAccessor(mLastPrimTree)
2792 {
2793 }
2794 
2795 
2796 inline
2798  : mTree(EdgeData())
2799  , mAccessor(mTree)
2800  , mPointList(rhs.mPointList)
2801  , mPolygonList(rhs.mPolygonList)
2802  , mLastPrimTree(Int32(util::INVALID_IDX))
2803  , mLastPrimAccessor(mLastPrimTree)
2804 {
2805 }
2806 
2807 
2808 inline void
2810 {
2811  if (threaded) {
2812  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
2813  } else {
2814  (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
2815  }
2816 }
2817 
2818 
2819 inline void
2821 {
2822  typedef TreeType::RootNodeType RootNodeType;
2823  typedef RootNodeType::NodeChainType NodeChainType;
2824  BOOST_STATIC_ASSERT(boost::mpl::size<NodeChainType>::value > 1);
2825  typedef boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
2826 
2827  Coord ijk;
2828  Index offset;
2829 
2830  rhs.mTree.clearAllAccessors();
2831 
2832  TreeType::LeafIter leafIt = rhs.mTree.beginLeaf();
2833  for ( ; leafIt; ++leafIt) {
2834  ijk = leafIt->origin();
2835 
2836  TreeType::LeafNodeType* lhsLeafPt = mTree.probeLeaf(ijk);
2837 
2838  if (!lhsLeafPt) {
2839 
2840  mAccessor.addLeaf(rhs.mAccessor.probeLeaf(ijk));
2841  InternalNodeType* node = rhs.mAccessor.getNode<InternalNodeType>();
2842  node->stealNode<TreeType::LeafNodeType>(ijk, EdgeData(), false);
2843  rhs.mAccessor.clear();
2844 
2845  } else {
2846 
2847  TreeType::LeafNodeType::ValueOnCIter it = leafIt->cbeginValueOn();
2848  for ( ; it; ++it) {
2849 
2850  offset = it.pos();
2851  const EdgeData& rhsValue = it.getValue();
2852 
2853  if (!lhsLeafPt->isValueOn(offset)) {
2854  lhsLeafPt->setValueOn(offset, rhsValue);
2855  } else {
2856 
2857  EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
2858 
2859  if (rhsValue.mXDist < lhsValue.mXDist) {
2860  lhsValue.mXDist = rhsValue.mXDist;
2861  lhsValue.mXPrim = rhsValue.mXPrim;
2862  }
2863 
2864  if (rhsValue.mYDist < lhsValue.mYDist) {
2865  lhsValue.mYDist = rhsValue.mYDist;
2866  lhsValue.mYPrim = rhsValue.mYPrim;
2867  }
2868 
2869  if (rhsValue.mZDist < lhsValue.mZDist) {
2870  lhsValue.mZDist = rhsValue.mZDist;
2871  lhsValue.mZPrim = rhsValue.mZPrim;
2872  }
2873 
2874  }
2875  } // end value iteration
2876  }
2877  } // end leaf iteration
2878 }
2879 
2880 
2881 inline void
2882 MeshToVoxelEdgeData::GenEdgeData::operator()(const tbb::blocked_range<size_t> &range)
2883 {
2884  Primitive prim;
2885 
2886  for (size_t n = range.begin(); n < range.end(); ++n) {
2887 
2888  const Vec4I& verts = mPolygonList[n];
2889 
2890  prim.index = Int32(n);
2891  prim.a = Vec3d(mPointList[verts[0]]);
2892  prim.b = Vec3d(mPointList[verts[1]]);
2893  prim.c = Vec3d(mPointList[verts[2]]);
2894 
2895  if (util::INVALID_IDX != verts[3]) {
2896  prim.d = Vec3d(mPointList[verts[3]]);
2897  voxelize<true>(prim);
2898  } else {
2899  voxelize<false>(prim);
2900  }
2901  }
2902 }
2903 
2904 
2905 template<bool IsQuad>
2906 inline void
2907 MeshToVoxelEdgeData::GenEdgeData::voxelize(const Primitive& prim)
2908 {
2909  std::deque<Coord> coordList;
2910  Coord ijk, nijk;
2911 
2912  ijk = util::nearestCoord(prim.a);
2913  coordList.push_back(ijk);
2914 
2915  evalPrimitive<IsQuad>(ijk, prim);
2916 
2917  while (!coordList.empty()) {
2918 
2919  ijk = coordList.back();
2920  coordList.pop_back();
2921 
2922  for (Int32 i = 0; i < 26; ++i) {
2923  nijk = ijk + util::COORD_OFFSETS[i];
2924 
2925  if (prim.index != mLastPrimAccessor.getValue(nijk)) {
2926  mLastPrimAccessor.setValue(nijk, prim.index);
2927  if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
2928  }
2929  }
2930  }
2931 }
2932 
2933 
2934 template<bool IsQuad>
2935 inline bool
2936 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(const Coord& ijk, const Primitive& prim)
2937 {
2938  Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
2939  bool intersecting = false;
2940  double t;
2941 
2942  EdgeData edgeData;
2943  mAccessor.probeValue(ijk, edgeData);
2944 
2945  // Evaluate first triangle
2946  double dist = (org -
2947  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, org, uvw)).lengthSqr();
2948 
2949  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
2950  if (t < edgeData.mXDist) {
2951  edgeData.mXDist = t;
2952  edgeData.mXPrim = prim.index;
2953  intersecting = true;
2954  }
2955  }
2956 
2957  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
2958  if (t < edgeData.mYDist) {
2959  edgeData.mYDist = t;
2960  edgeData.mYPrim = prim.index;
2961  intersecting = true;
2962  }
2963  }
2964 
2965  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
2966  if (t < edgeData.mZDist) {
2967  edgeData.mZDist = t;
2968  edgeData.mZPrim = prim.index;
2969  intersecting = true;
2970  }
2971  }
2972 
2973  if (IsQuad) {
2974  // Split quad into a second triangle and calculate distance.
2975  double secondDist = (org -
2976  closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, org, uvw)).lengthSqr();
2977 
2978  if (secondDist < dist) dist = secondDist;
2979 
2980  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
2981  if (t < edgeData.mXDist) {
2982  edgeData.mXDist = t;
2983  edgeData.mXPrim = prim.index;
2984  intersecting = true;
2985  }
2986  }
2987 
2988  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
2989  if (t < edgeData.mYDist) {
2990  edgeData.mYDist = t;
2991  edgeData.mYPrim = prim.index;
2992  intersecting = true;
2993  }
2994  }
2995 
2996  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
2997  if (t < edgeData.mZDist) {
2998  edgeData.mZDist = t;
2999  edgeData.mZPrim = prim.index;
3000  intersecting = true;
3001  }
3002  }
3003  }
3004 
3005  if (intersecting) mAccessor.setValue(ijk, edgeData);
3006 
3007  return (dist < 0.86602540378443861);
3008 }
3009 
3010 
3011 inline bool
3012 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3013  const Vec3d& origin, const Vec3d& dir,
3014  const Vec3d& a, const Vec3d& b, const Vec3d& c,
3015  double& t)
3016 {
3017  // Check if ray is parallel with triangle
3018 
3019  Vec3d e1 = b - a;
3020  Vec3d e2 = c - a;
3021  Vec3d s1 = dir.cross(e2);
3022 
3023  double divisor = s1.dot(e1);
3024  if (!(std::abs(divisor) > 0.0)) return false;
3025 
3026  // Compute barycentric coordinates
3027 
3028  double inv_divisor = 1.0 / divisor;
3029  Vec3d d = origin - a;
3030  double b1 = d.dot(s1) * inv_divisor;
3031 
3032  if (b1 < 0.0 || b1 > 1.0) return false;
3033 
3034  Vec3d s2 = d.cross(e1);
3035  double b2 = dir.dot(s2) * inv_divisor;
3036 
3037  if (b2 < 0.0 || (b1 + b2) > 1.0) return false;
3038 
3039  // Compute distance to intersection point
3040 
3041  t = e2.dot(s2) * inv_divisor;
3042  return (t < 0.0) ? false : true;
3043 }
3044 
3045 
3047 
3048 
3049 inline
3051  : mTree(EdgeData())
3052 {
3053 }
3054 
3055 
3056 inline void
3058  const std::vector<Vec3s>& pointList,
3059  const std::vector<Vec4I>& polygonList)
3060 {
3061  GenEdgeData converter(pointList, polygonList);
3062  converter.run();
3063 
3064  mTree.clear();
3065  mTree.merge(converter.tree());
3066 }
3067 
3068 
3069 inline void
3071  Accessor& acc,
3072  const Coord& ijk,
3073  std::vector<Vec3d>& points,
3074  std::vector<Index32>& primitives)
3075 {
3076  EdgeData data;
3077  Vec3d point;
3078 
3079  Coord coord = ijk;
3080 
3081  if (acc.probeValue(coord, data)) {
3082 
3083  if (data.mXPrim != util::INVALID_IDX) {
3084  point[0] = double(coord[0]) + data.mXDist;
3085  point[1] = double(coord[1]);
3086  point[2] = double(coord[2]);
3087 
3088  points.push_back(point);
3089  primitives.push_back(data.mXPrim);
3090  }
3091 
3092  if (data.mYPrim != util::INVALID_IDX) {
3093  point[0] = double(coord[0]);
3094  point[1] = double(coord[1]) + data.mYDist;
3095  point[2] = double(coord[2]);
3096 
3097  points.push_back(point);
3098  primitives.push_back(data.mYPrim);
3099  }
3100 
3101  if (data.mZPrim != util::INVALID_IDX) {
3102  point[0] = double(coord[0]);
3103  point[1] = double(coord[1]);
3104  point[2] = double(coord[2]) + data.mZDist;
3105 
3106  points.push_back(point);
3107  primitives.push_back(data.mZPrim);
3108  }
3109 
3110  }
3111 
3112  coord[0] += 1;
3113 
3114  if (acc.probeValue(coord, data)) {
3115 
3116  if (data.mYPrim != util::INVALID_IDX) {
3117  point[0] = double(coord[0]);
3118  point[1] = double(coord[1]) + data.mYDist;
3119  point[2] = double(coord[2]);
3120 
3121  points.push_back(point);
3122  primitives.push_back(data.mYPrim);
3123  }
3124 
3125  if (data.mZPrim != util::INVALID_IDX) {
3126  point[0] = double(coord[0]);
3127  point[1] = double(coord[1]);
3128  point[2] = double(coord[2]) + data.mZDist;
3129 
3130  points.push_back(point);
3131  primitives.push_back(data.mZPrim);
3132  }
3133  }
3134 
3135  coord[2] += 1;
3136 
3137  if (acc.probeValue(coord, data)) {
3138  if (data.mYPrim != util::INVALID_IDX) {
3139  point[0] = double(coord[0]);
3140  point[1] = double(coord[1]) + data.mYDist;
3141  point[2] = double(coord[2]);
3142 
3143  points.push_back(point);
3144  primitives.push_back(data.mYPrim);
3145  }
3146  }
3147 
3148  coord[0] -= 1;
3149 
3150  if (acc.probeValue(coord, data)) {
3151 
3152  if (data.mXPrim != util::INVALID_IDX) {
3153  point[0] = double(coord[0]) + data.mXDist;
3154  point[1] = double(coord[1]);
3155  point[2] = double(coord[2]);
3156 
3157  points.push_back(point);
3158  primitives.push_back(data.mXPrim);
3159  }
3160 
3161  if (data.mYPrim != util::INVALID_IDX) {
3162  point[0] = double(coord[0]);
3163  point[1] = double(coord[1]) + data.mYDist;
3164  point[2] = double(coord[2]);
3165 
3166  points.push_back(point);
3167  primitives.push_back(data.mYPrim);
3168  }
3169  }
3170 
3171 
3172  coord[1] += 1;
3173 
3174  if (acc.probeValue(coord, data)) {
3175 
3176  if (data.mXPrim != util::INVALID_IDX) {
3177  point[0] = double(coord[0]) + data.mXDist;
3178  point[1] = double(coord[1]);
3179  point[2] = double(coord[2]);
3180 
3181  points.push_back(point);
3182  primitives.push_back(data.mXPrim);
3183  }
3184  }
3185 
3186  coord[2] -= 1;
3187 
3188  if (acc.probeValue(coord, data)) {
3189 
3190  if (data.mXPrim != util::INVALID_IDX) {
3191  point[0] = double(coord[0]) + data.mXDist;
3192  point[1] = double(coord[1]);
3193  point[2] = double(coord[2]);
3194 
3195  points.push_back(point);
3196  primitives.push_back(data.mXPrim);
3197  }
3198 
3199  if (data.mZPrim != util::INVALID_IDX) {
3200  point[0] = double(coord[0]);
3201  point[1] = double(coord[1]);
3202  point[2] = double(coord[2]) + data.mZDist;
3203 
3204  points.push_back(point);
3205  primitives.push_back(data.mZPrim);
3206  }
3207  }
3208 
3209  coord[0] += 1;
3210 
3211  if (acc.probeValue(coord, data)) {
3212 
3213  if (data.mZPrim != util::INVALID_IDX) {
3214  point[0] = double(coord[0]);
3215  point[1] = double(coord[1]);
3216  point[2] = double(coord[2]) + data.mZDist;
3217 
3218  points.push_back(point);
3219  primitives.push_back(data.mZPrim);
3220  }
3221  }
3222 }
3223 
3224 
3225 } // namespace tools
3226 } // namespace OPENVDB_VERSION_NAME
3227 } // namespace openvdb
3228 
3229 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
3230 
3231 // Copyright (c) 2012-2013 DreamWorks Animation LLC
3232 // All rights reserved. This software is distributed under the
3233 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1102
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
void operator()(const tbb::blocked_range< size_t > &)
Definition: MeshToVolume.h:1827
BoolTreeT & signMaskTree()
Definition: MeshToVolume.h:1115
Definition: Tree.h:178
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1709
BoolTreeT & signMaskTree()
Definition: MeshToVolume.h:966
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1560
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
IntersectingVoxelSign(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, FloatTreeT &distTree, IntTreeT &indexTree, BoolTreeT &intersectionTree, BoolLeafManager &leafs)
Definition: MeshToVolume.h:1296
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:533
~IntersectingVoxelCleaner()
Definition: MeshToVolume.h:1437
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:1760
void combine(FloatTreeT &lhsDist, IntTreeT &lhsIndex, FloatTreeT &rhsDist, IntTreeT &rhsIndex)
Definition: MeshToVolume.h:396
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:1248
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:2702
~SignMask()
Definition: MeshToVolume.h:958
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:303
Definition: MeshToVolume.h:2053
void run(bool threaded=true)
Definition: MeshToVolume.h:1586
tree::LeafManager< FloatTreeT > DistArrayT
Definition: MeshToVolume.h:1556
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1250
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1247
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2062
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1715
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:951
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
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:442
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:372
SqrtAndScaleOp(ValueType voxelSize, bool unsignedDist=false)
Definition: MeshToVolume.h:2026
Scheme::template ISStencil< GridType >::StencilType Stencil
Definition: MeshToVolume.h:2140
tree::ValueAccessor< TreeType > Accessor
Definition: MeshToVolume.h:308
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:653
PointTransform(const std::vector< Vec3s > &pointsIn, std::vector< Vec3s > &pointsOut, const math::Transform &xform)
Definition: MeshToVolume.h:355
~PropagateSign()
Definition: MeshToVolume.h:1107
tree::ValueAccessor< const BoolTreeT > BoolConstAccessorT
Definition: MeshToVolume.h:1103
TBB body object that traversers all intersecting voxels (defined by the intersectingVoxelsGrid) and p...
Definition: MeshToVolume.h:1244
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1426
bool operator==(const EdgeData &rhs) const
Definition: MeshToVolume.h:298
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:2729
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1430
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1370
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1713
Internal edge data type.
Definition: MeshToVolume.h:279
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:440
TBB body object that removes intersecting voxels that were set via voxelization of self-intersecting ...
Definition: MeshToVolume.h:1420
void run(bool threaded=true)
Definition: MeshToVolume.h:1010
Index32 Index
Definition: Types.h:57
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1712
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
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
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1096
void convertToUnsignedDistanceField(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, FloatValueT exBandWidth)
Mesh to Unsigned Distance Field conversion.
Definition: MeshToVolume.h:2332
std::ostream & operator<<(std::ostream &ostr, const MeshToVoxelEdgeData::EdgeData &rhs)
Definition: MeshToVolume.h:2719
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1100
OPENVDB_API const Index32 INVALID_IDX
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1716
Definition: Mat.h:146
GridType::Ptr meshToLevelSet(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float halfWidth=float(LEVEL_SET_HALF_WIDTH))
Convert a triangle and quad mesh to a level set volume.
Definition: MeshToVolume.h:2673
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1251
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:445
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2192
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2033
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2091
FloatTreeT & sqrDistTree()
Definition: MeshToVolume.h:461
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:946
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1064
Definition: MeshToVolume.h:2184
MeshToVolume(openvdb::math::Transform::Ptr &, int conversionFlags=0, InterruptT *interrupter=NULL, int signSweeps=1)
Definition: MeshToVolume.h:2292
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:1554
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
tree::LeafManager< TreeType > LeafManagerType
Definition: MeshToVolume.h:2211
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1557
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1427
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1253
FloatGridT::TreeType FloatTreeT
Definition: MeshToVolume.h:199
TBB body object to expand the level set narrow band.
Definition: MeshToVolume.h:1706
IntersectingVoxelCleaner(FloatTreeT &distTree, IntTreeT &indexTree, BoolTreeT &intersectionTree, BoolLeafManager &leafs)
Definition: MeshToVolume.h:1466
void run(bool threaded=true)
Definition: MeshToVolume.h:1806
Type Pow2(Type x)
Return .
Definition: Math.h:458
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:176
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1428
tree::LeafManager< FloatTreeT > FloatLeafManager
Definition: MeshToVolume.h:947
~ContourTracer()
Definition: MeshToVolume.h:737
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:733
Definition: MeshToVolume.h:2260
TBB body object that that finds seed points for the parallel flood fill.
Definition: MeshToVolume.h:942
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2221
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1432
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:1098
SDFPrune(const ValueType &out, const ValueType &in)
Definition: MeshToVolume.h:2262
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1423
TrimOp(ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2084
EdgeData(float dist=1.0)
Definition: MeshToVolume.h:280
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:439
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:1019
void run(bool threaded=true)
Definition: MeshToVolume.h:1456
OPENVDB_API Vec3d closestPointOnTriangleToPoint(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p, Vec3d &uvw)
Closest Point on Triangle to Point. Given a triangle abc and a point p, return the point on abc close...
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:443
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:1602
boost::shared_ptr< Grid > Ptr
Definition: Grid.h:461
Definition: Types.h:170
static ValueType epsilon()
Definition: MeshToVolume.h:389
void join(MeshVoxelizer< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:646
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:200
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:2572
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:949
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:210
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:733
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:447
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:1249
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:54
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:228
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:444
Definition: Exceptions.h:87
const ValueType outside
Definition: MeshToVolume.h:2278
#define OPENVDB_VERSION_NAME
Definition: version.h:45
EdgeData operator-() const
Definition: MeshToVolume.h:295
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1693
Definition: MeshToVolume.h:2082
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:199
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 > FloatAccessorT
Definition: MeshToVolume.h:1711
static Accessor::ValueType result(const Accessor &grid, const Coord &ijk)
Definition: Operators.h:260
math::BIAS_SCHEME< math::FIRST_BIAS > Scheme
Definition: MeshToVolume.h:2139
void run(bool threaded=true)
Definition: MeshToVolume.h:1288
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
tree::ValueAccessor< TreeType > AccessorT
Definition: MeshToVolume.h:2242
TBB body object that removes non-intersecting voxels that where set by rasterizing self-intersecting ...
Definition: MeshToVolume.h:1550
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:351
~IntersectingVoxelSign()
Definition: MeshToVolume.h:1263
EdgeData operator-(const T &) const
Definition: MeshToVolume.h:294
Grid< IntTreeT > IntGridT
Definition: MeshToVolume.h:202
Coord nearestCoord(const Vec3d &voxelCoord)
Return voxelCoord rounded to the closest integer coordinates.
Definition: util/Util.h:55
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2187
OffsetOp(ValueType offset)
Definition: MeshToVolume.h:2117
~ExpandNB()
Definition: MeshToVolume.h:1730
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1101
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
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1714
TBB body object that partitions a volume into 2D slices that can be processed in parallel and marks t...
Definition: MeshToVolume.h:728
Index32 mYPrim
Definition: MeshToVolume.h:304
MeshToVoxelEdgeData()
Definition: MeshToVolume.h:3050
void join(SignMask< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:1082
~ShellVoxelCleaner()
Definition: MeshToVolume.h:1567
void run(bool threaded=true)
Definition: MeshToVolume.h:1157
RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:186
ScalarToVectorConverter< GridType >::Type::Ptr cpt(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Closest-Point Transform (CPT) from a distance field.
Definition: GridOperators.h:950
void join(GenEdgeData &rhs)
Definition: MeshToVolume.h:2820
GenEdgeData(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition: MeshToVolume.h:2783
static ValueType minNarrowBandWidth()
Definition: MeshToVolume.h:390
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1429
tree::LeafManager< typename GridType::TreeType > LeafManagerType
Definition: MeshToVolume.h:2141
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:158
TreeType & tree()
Definition: MeshToVolume.h:2752
tree::ValueAccessor< const BoolTreeT > BoolConstAccessorT
Definition: MeshToVolume.h:952
Definition: Operators.h:152
Extracts and stores voxel edge intersection data from a mesh.
Definition: MeshToVolume.h:272
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:732
BoolTreeT & intersectionTree()
Definition: MeshToVolume.h:463
uint32_t Index32
Definition: Types.h:55
void operator()(LeafNodeType &leaf, size_t leafIndex) const
Definition: MeshToVolume.h:2155
int32_t Int32
Definition: Types.h:59
MeshVoxelizer(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:513
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2882
math::Vec4< Index32 > Vec4I
Definition: Types.h:89
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:658
IntTreeT::LeafNodeType IntLeafT
Definition: MeshToVolume.h:1559
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1211
float mZDist
Definition: MeshToVolume.h:303
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
Definition: LeafManager.h:463
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:1097
tree::ValueAccessor< IntTreeT > IntAccessorT
Definition: MeshToVolume.h:1558
Vec3< double > Vec3d
Definition: Vec3.h:625
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:3070
_RootNodeType RootNodeType
Definition: Tree.h:184
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:2645
MergeBufferOp(LeafManagerType &leafs, size_t bufferIndex=1)
Definition: MeshToVolume.h:2214
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:203
Definition: Types.h:169
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:1688
void run(bool threaded=true)
Definition: MeshToVolume.h:766
Definition: MeshToVolume.h:191
void run(bool threaded=true)
Definition: MeshToVolume.h:2809
static const bool state
Definition: MeshToVolume.h:2277
SignMask(const FloatLeafManager &, const FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:984
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1623
VoxelSignOp(ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2055
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2212
Definition: MeshToVolume.h:196
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1562
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1718
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:3057
~MeshVoxelizer()
Definition: MeshToVolume.h:453
Definition: MeshToVolume.h:2137
float mXDist
Definition: MeshToVolume.h:303
Index32 mXPrim
Definition: MeshToVolume.h:304
tree::ValueAccessor< const FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:948
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:731
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:2317
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1328
const ValueType inside
Definition: MeshToVolume.h:2278
void operator()(const tbb::blocked_range< size_t > &) const
Definition: MeshToVolume.h:1492
IntTreeT & primIndexTree()
Definition: MeshToVolume.h:462
Definition: MeshToVolume.h:2115
void clear()
Definition: MeshToVolume.h:2307
CopyActiveVoxelsOp(TreeType &tree)
Definition: MeshToVolume.h:1685
tree::Tree4< EdgeData, 5, 4, 3 >::Type TreeType
Definition: MeshToVolume.h:307
RenormOp(GridType &grid, LeafManagerType &leafs, ValueType voxelSize, ValueType cfl=1.0)
Definition: MeshToVolume.h:2144
float mYDist
Definition: MeshToVolume.h:303
void resetOffset(ValueType offset)
Definition: MeshToVolume.h:2119
tree::ValueAccessor< FloatTreeT > FloatAccessorT
Definition: MeshToVolume.h:441
FloatTreeT::template ValueConverter< bool >::Type BoolTreeT
Definition: MeshToVolume.h:1099
EdgeData operator+(const T &) const
Definition: MeshToVolume.h:293
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
ShellVoxelCleaner(FloatTreeT &distTree, DistArrayT &leafs, IntTreeT &indexTree, BoolTreeT &intersectionTree)
Definition: MeshToVolume.h:1597
Index32 mZPrim
Definition: MeshToVolume.h:304
LeafManagerType::BufferType BufferType
Definition: MeshToVolume.h:2142
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
Definition: MeshToVolume.h:387
MinOp(LeafManagerType &leafs)
Definition: MeshToVolume.h:2189
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:170
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:217
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:2687
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
FloatTreeT::template ValueConverter< Int32 >::Type IntTreeT
Definition: MeshToVolume.h:201
bool isNegative(const Type &x)
Return true if x is less than zero.
Definition: Math.h:304
TBB body object to voxelize a mesh of triangles and/or quads into a collection of VDB grids...
Definition: MeshToVolume.h:436
TreeType::LeafNodeType LeafNodeT
Definition: MeshToVolume.h:2243
LeafTopologyDiffOp(TreeType &tree)
Definition: MeshToVolume.h:2245
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:945
tree::ValueAccessor< const BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:734
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1252
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:446
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:198
tree::LeafManager< BoolTreeT > BoolLeafManager
Definition: MeshToVolume.h:1717
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:550
FloatTreeT::LeafNodeType DistLeafT
Definition: MeshToVolume.h:1555
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
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:1431
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
PropagateSign(BoolLeafManager &, FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:1132
TBB body object that performs a parallel flood fill.
Definition: MeshToVolume.h:1093
void run(bool threaded=true)
Definition: MeshToVolume.h:363
boost::shared_ptr< Transform > Ptr
Definition: Transform.h:68
OPENVDB_STATIC_SPECIALIZATION void dilateVoxels(TreeType &tree, int count=1)
Definition: Morphology.h:364
void join(PropagateSign< FloatTreeT, InterruptT > &rhs)
Definition: MeshToVolume.h:1230
Accessor getAccessor()
Definition: MeshToVolume.h:334
tree::ValueAccessor< FloatTreeT > DistAccessorT
Definition: MeshToVolume.h:1424
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2248
ValueType value
Definition: MeshToVolume.h:2279
tree::LeafManager< TreeType > LeafManagerType
Definition: MeshToVolume.h:2186
IntGridT::Ptr indexGridPtr() const
Definition: MeshToVolume.h:244
bool operator()(ChildType &child)
Definition: MeshToVolume.h:2269
FloatTreeT::ValueType FloatValueT
Definition: MeshToVolume.h:1553
void resetCFL(ValueType cfl)
Definition: MeshToVolume.h:2152
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
ContourTracer(FloatTreeT &, const BoolTreeT &, InterruptT *interrupter=NULL)
Definition: MeshToVolume.h:777
tree::ValueAccessor< TreeType > AccessorT
Definition: MeshToVolume.h:1683
FloatTreeT::LeafNodeType FloatLeafT
Definition: MeshToVolume.h:1710
void operator()(const tbb::blocked_range< int > &range) const
Definition: MeshToVolume.h:838
void join(ExpandNB< FloatTreeT > &)
Definition: MeshToVolume.h:2012
BoolTreeT::LeafNodeType BoolLeafT
Definition: MeshToVolume.h:950
Definition: Mat4.h:51
tree::ValueAccessor< BoolTreeT > BoolAccessorT
Definition: MeshToVolume.h:1561
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:65
FloatGridT::Ptr distGridPtr() const
Returns a narrow-band (signed) distance field / level set grid.
Definition: MeshToVolume.h:240
Grid< BoolTreeT > BoolGridT
Definition: MeshToVolume.h:204
void operator()(LeafNodeType &leaf, size_t) const
Definition: MeshToVolume.h:2122
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:1166