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