OpenVDB  3.1.0
MeshToVolume.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2015 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 //
42 
43 
44 #ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
45 #define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
46 
47 #include <openvdb/Types.h>
48 #include <openvdb/math/FiniteDifference.h> // for GudonovsNormSqrd
49 #include <openvdb/math/Proximity.h> // for closestPointOnTriangleToPoint()
51 #include <openvdb/util/Util.h>
52 
53 #include "ChangeBackground.h"
54 #include "Prune.h" // for pruneInactive and pruneLevelSet
55 #include "SignedFloodFill.h" // for signedFloodFillWithValues
56 
57 #include <tbb/blocked_range.h>
58 #include <tbb/enumerable_thread_specific.h>
59 #include <tbb/parallel_for.h>
60 #include <tbb/parallel_reduce.h>
61 #include <tbb/partitioner.h>
62 #include <tbb/task_group.h>
63 #include <tbb/task_scheduler_init.h>
64 #include <tbb/tick_count.h>
65 
66 #include <boost/integer_traits.hpp> // const_max
67 #include <boost/math/special_functions/fpclassify.hpp> // for isfinite()
68 #include <boost/scoped_array.hpp>
69 
70 #include <deque>
71 #include <limits>
72 #include <sstream>
73 
74 
75 namespace openvdb {
77 namespace OPENVDB_VERSION_NAME {
78 namespace tools {
79 
80 
82 
83 
86 
92 
96 
100 
104 };
105 
106 
137 template <typename GridType, typename MeshDataAdapter>
138 inline typename GridType::Ptr
140  const MeshDataAdapter& mesh,
141  const math::Transform& transform,
142  float exteriorBandWidth = 3.0f,
143  float interiorBandWidth = 3.0f,
144  int flags = 0,
145  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = NULL);
146 
147 
163 template <typename GridType, typename MeshDataAdapter, typename Interrupter>
164 inline typename GridType::Ptr
166  Interrupter& interrupter,
167  const MeshDataAdapter& mesh,
168  const math::Transform& transform,
169  float exteriorBandWidth = 3.0f,
170  float interiorBandWidth = 3.0f,
171  int flags = 0,
172  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = NULL);
173 
174 
176 
177 
188 template<typename PointType, typename PolygonType>
190 
191  QuadAndTriangleDataAdapter(const std::vector<PointType>& points,
192  const std::vector<PolygonType>& polygons)
193  : mPointArray(points.empty() ? NULL : &points[0])
194  , mPointArraySize(points.size())
195  , mPolygonArray(polygons.empty() ? NULL : &polygons[0])
196  , mPolygonArraySize(polygons.size())
197  {
198  }
199 
200  QuadAndTriangleDataAdapter(const PointType * pointArray, size_t pointArraySize,
201  const PolygonType* polygonArray, size_t polygonArraySize)
202  : mPointArray(pointArray)
203  , mPointArraySize(pointArraySize)
204  , mPolygonArray(polygonArray)
205  , mPolygonArraySize(polygonArraySize)
206  {
207  }
208 
209  size_t polygonCount() const { return mPolygonArraySize; }
210  size_t pointCount() const { return mPointArraySize; }
211 
213  size_t vertexCount(size_t n) const {
214  return (PolygonType::size == 3 || mPolygonArray[n][3] == util::INVALID_IDX) ? 3 : 4;
215  }
216 
219  void getIndexSpacePoint(size_t n, size_t v, Vec3d& pos) const {
220  const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
221  pos[0] = double(p[0]);
222  pos[1] = double(p[1]);
223  pos[2] = double(p[2]);
224  }
225 
226 private:
227  PointType const * const mPointArray;
228  size_t const mPointArraySize;
229  PolygonType const * const mPolygonArray;
230  size_t const mPolygonArraySize;
231 }; // struct QuadAndTriangleDataAdapter
232 
233 
235 
236 
237 // Wrapper functions for the mesh to volume converter
238 
239 
255 template<typename GridType>
256 inline typename GridType::Ptr
258  const openvdb::math::Transform& xform,
259  const std::vector<Vec3s>& points,
260  const std::vector<Vec3I>& triangles,
261  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
262 
263 
279 template<typename GridType>
280 inline typename GridType::Ptr
282  const openvdb::math::Transform& xform,
283  const std::vector<Vec3s>& points,
284  const std::vector<Vec4I>& quads,
285  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
286 
287 
304 template<typename GridType>
305 inline typename GridType::Ptr
307  const openvdb::math::Transform& xform,
308  const std::vector<Vec3s>& points,
309  const std::vector<Vec3I>& triangles,
310  const std::vector<Vec4I>& quads,
311  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
312 
313 
332 template<typename GridType>
333 inline typename GridType::Ptr
335  const openvdb::math::Transform& xform,
336  const std::vector<Vec3s>& points,
337  const std::vector<Vec3I>& triangles,
338  const std::vector<Vec4I>& quads,
339  float exBandWidth,
340  float inBandWidth);
341 
342 
357 template<typename GridType>
358 inline typename GridType::Ptr
360  const openvdb::math::Transform& xform,
361  const std::vector<Vec3s>& points,
362  const std::vector<Vec3I>& triangles,
363  const std::vector<Vec4I>& quads,
364  float bandWidth);
365 
366 
368 
369 
376 template<typename GridType, typename VecType>
377 inline typename GridType::Ptr
379  const openvdb::math::Transform& xform,
380  typename VecType::ValueType halfWidth = LEVEL_SET_HALF_WIDTH);
381 
382 
384 
385 
392 template <typename FloatTreeT>
393 inline void
394 traceExteriorBoundaries(FloatTreeT& tree);
395 
396 
398 
399 
402 {
403 public:
404 
406 
408  struct EdgeData {
409  EdgeData(float dist = 1.0)
410  : mXDist(dist), mYDist(dist), mZDist(dist)
411  , mXPrim(util::INVALID_IDX)
412  , mYPrim(util::INVALID_IDX)
413  , mZPrim(util::INVALID_IDX)
414  {
415  }
416 
418  bool operator< (const EdgeData&) const { return false; }
421  bool operator> (const EdgeData&) const { return false; }
422  template<class T> EdgeData operator+(const T&) const { return *this; }
423  template<class T> EdgeData operator-(const T&) const { return *this; }
424  EdgeData operator-() const { return *this; }
426 
427  bool operator==(const EdgeData& rhs) const
428  {
429  return mXPrim == rhs.mXPrim && mYPrim == rhs.mYPrim && mZPrim == rhs.mZPrim;
430  }
431 
432  float mXDist, mYDist, mZDist;
433  Index32 mXPrim, mYPrim, mZPrim;
434  };
435 
438 
439 
441 
442 
444 
445 
453  void convert(const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
454 
455 
458  void getEdgeData(Accessor& acc, const Coord& ijk,
459  std::vector<Vec3d>& points, std::vector<Index32>& primitives);
460 
463  Accessor getAccessor() { return Accessor(mTree); }
464 
465 private:
466  void operator=(const MeshToVoxelEdgeData&) {}
467  TreeType mTree;
468  class GenEdgeData;
469 };
470 
471 
474 
475 
476 // Internal utility objects and implementation details
477 
478 namespace mesh_to_volume_internal {
479 
480 template<typename PointType>
482 
483  TransformPoints(const PointType* pointsIn, PointType* pointsOut,
484  const math::Transform& xform)
485  : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
486  {
487  }
488 
489  void operator()(const tbb::blocked_range<size_t>& range) const {
490 
491  Vec3d pos;
492 
493  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
494 
495  const PointType& wsP = mPointsIn[n];
496  pos[0] = double(wsP[0]);
497  pos[1] = double(wsP[1]);
498  pos[2] = double(wsP[2]);
499 
500  pos = mXform->worldToIndex(pos);
501 
502  PointType& isP = mPointsOut[n];
503  isP[0] = typename PointType::value_type(pos[0]);
504  isP[1] = typename PointType::value_type(pos[1]);
505  isP[2] = typename PointType::value_type(pos[2]);
506  }
507  }
508 
509  PointType const * const mPointsIn;
510  PointType * const mPointsOut;
511  math::Transform const * const mXform;
512 }; // TransformPoints
513 
514 
515 template<typename ValueType>
516 struct Tolerance
517 {
518  static ValueType epsilon() { return ValueType(1e-7); }
519  static ValueType minNarrowBandWidth() { return ValueType(1.0 + 1e-6); }
520 };
521 
522 
524 
525 
526 template<typename TreeType>
528 {
529 public:
530 
531  typedef typename TreeType::template ValueConverter<Int32>::Type Int32TreeType;
532 
533  typedef typename TreeType::LeafNodeType LeafNodeType;
534  typedef typename Int32TreeType::LeafNodeType Int32LeafNodeType;
535 
536  CombineLeafNodes(TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
537  LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
538  : mDistTree(&lhsDistTree)
539  , mIdxTree(&lhsIdxTree)
540  , mRhsDistNodes(rhsDistNodes)
541  , mRhsIdxNodes(rhsIdxNodes)
542  {
543  }
544 
545  void operator()(const tbb::blocked_range<size_t>& range) const {
546 
547  tree::ValueAccessor<TreeType> distAcc(*mDistTree);
548  tree::ValueAccessor<Int32TreeType> idxAcc(*mIdxTree);
549 
550  typedef typename LeafNodeType::ValueType DistValueType;
551  typedef typename Int32LeafNodeType::ValueType IndexValueType;
552 
553  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
554 
555  const Coord& origin = mRhsDistNodes[n]->origin();
556 
557  LeafNodeType* lhsDistNode = distAcc.probeLeaf(origin);
558  Int32LeafNodeType* lhsIdxNode = idxAcc.probeLeaf(origin);
559 
560  DistValueType* lhsDistData = lhsDistNode->buffer().data();
561  IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
562 
563  const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
564  const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
565 
566 
567  for (Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
568 
569  if (rhsIdxData[offset] != Int32(util::INVALID_IDX)) {
570 
571  const DistValueType& lhsValue = lhsDistData[offset];
572  const DistValueType& rhsValue = rhsDistData[offset];
573 
574  if (rhsValue < lhsValue) {
575  lhsDistNode->setValueOn(offset, rhsValue);
576  lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
577  } else if (math::isExactlyEqual(rhsValue, lhsValue)) {
578  lhsIdxNode->setValueOn(offset,
579  std::min(lhsIdxData[offset], rhsIdxData[offset]));
580  }
581  }
582  }
583 
584  delete mRhsDistNodes[n];
585  delete mRhsIdxNodes[n];
586  }
587  }
588 
589 private:
590 
591  TreeType * const mDistTree;
592  Int32TreeType * const mIdxTree;
593 
594  LeafNodeType ** const mRhsDistNodes;
595  Int32LeafNodeType ** const mRhsIdxNodes;
596 }; // class CombineLeafNodes
597 
598 
600 
601 
602 template<typename TreeType>
604 {
605  typedef typename TreeType::LeafNodeType LeafNodeType;
606 
607  StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
608  : mNodes(nodes.empty() ? NULL : &nodes[0]), mCoordinates(coordinates)
609  {
610  }
611 
612  void operator()(const tbb::blocked_range<size_t>& range) const {
613  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
614  Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
615  mCoordinates[n] = origin;
616  origin[0] = static_cast<int>(n);
617  }
618  }
619 
620  LeafNodeType ** const mNodes;
621  Coord * const mCoordinates;
622 };
623 
624 
625 template<typename TreeType>
627 {
628  typedef typename TreeType::LeafNodeType LeafNodeType;
629 
630  RestoreOrigin(std::vector<LeafNodeType*>& nodes, const Coord* coordinates)
631  : mNodes(nodes.empty() ? NULL : &nodes[0]), mCoordinates(coordinates)
632  {
633  }
634 
635  void operator()(const tbb::blocked_range<size_t>& range) const {
636  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
637  Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
638  origin[0] = mCoordinates[n][0];
639  }
640  }
641 
642  LeafNodeType ** const mNodes;
643  Coord const * const mCoordinates;
644 };
645 
646 
647 template<typename TreeType>
649 {
650 public:
651  typedef typename TreeType::LeafNodeType LeafNodeType;
652 
653  ComputeNodeConnectivity(const TreeType& tree, const Coord* coordinates,
654  size_t* offsets, size_t numNodes, const CoordBBox& bbox)
655  : mTree(&tree)
656  , mCoordinates(coordinates)
657  , mOffsets(offsets)
658  , mNumNodes(numNodes)
659  , mBBox(bbox)
660  {
661  }
662 
663  void operator()(const tbb::blocked_range<size_t>& range) const {
664 
665  size_t* offsetsNextX = mOffsets;
666  size_t* offsetsPrevX = mOffsets + mNumNodes;
667  size_t* offsetsNextY = mOffsets + mNumNodes * 2;
668  size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
669  size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
670  size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
671 
673  Coord ijk;
674 
675  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
676  const Coord& origin = mCoordinates[n];
677  offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(LeafNodeType::DIM, 0, 0));
678  offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-LeafNodeType::DIM, 0, 0));
679  offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, LeafNodeType::DIM, 0));
680  offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -LeafNodeType::DIM, 0));
681  offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, LeafNodeType::DIM));
682  offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -LeafNodeType::DIM));
683  }
684  }
685 
686  size_t findNeighbourNode(tree::ValueAccessor<const TreeType>& acc, const Coord& start, const Coord& step) const {
687 
688  Coord ijk = start + step;
689  CoordBBox bbox(mBBox);
690 
691  while (bbox.isInside(ijk)) {
692  const LeafNodeType* node = acc.probeConstLeaf(ijk);
693  if (node) return static_cast<size_t>(node->origin()[0]);
694  ijk += step;
695  }
696 
697  return boost::integer_traits<size_t>::const_max;
698  }
699 
700 
701 private:
702  // Disallow assignment
704 
705  TreeType const * const mTree;
706  Coord const * const mCoordinates;
707  size_t * const mOffsets;
708 
709  const size_t mNumNodes;
710  const CoordBBox mBBox;
711 }; // class ComputeNodeConnectivity
712 
713 
714 template<typename TreeType>
716 
717  enum { INVALID_OFFSET = boost::integer_traits<size_t>::const_max };
718 
719  typedef typename TreeType::LeafNodeType LeafNodeType;
720 
722  : mLeafNodes()
723  , mOffsets(NULL)
724  {
725  mLeafNodes.reserve(tree.leafCount());
726  tree.getNodes(mLeafNodes);
727 
728  if (mLeafNodes.empty()) return;
729 
730  CoordBBox bbox;
731  tree.evalLeafBoundingBox(bbox);
732 
733  const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
734 
735  // stash the leafnode origin coordinate and temporarily store the
736  // linear offset in the origin.x variable.
737  boost::scoped_array<Coord> coordinates(new Coord[mLeafNodes.size()]);
738  tbb::parallel_for(range, StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
739 
740  // build the leafnode offset table
741  mOffsets.reset(new size_t[mLeafNodes.size() * 6]);
742 
743 
744  tbb::parallel_for(range,
745  ComputeNodeConnectivity<TreeType>(tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
746 
747  // restore the leafnode origin coordinate
748  tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
749  }
750 
751  size_t size() const { return mLeafNodes.size(); }
752 
753  std::vector<LeafNodeType*>& nodes() { return mLeafNodes; }
754  const std::vector<LeafNodeType*>& nodes() const { return mLeafNodes; }
755 
756 
757  const size_t* offsetsNextX() const { return mOffsets.get(); }
758  const size_t* offsetsPrevX() const { return mOffsets.get() + mLeafNodes.size(); }
759 
760  const size_t* offsetsNextY() const { return mOffsets.get() + mLeafNodes.size() * 2; }
761  const size_t* offsetsPrevY() const { return mOffsets.get() + mLeafNodes.size() * 3; }
762 
763  const size_t* offsetsNextZ() const { return mOffsets.get() + mLeafNodes.size() * 4; }
764  const size_t* offsetsPrevZ() const { return mOffsets.get() + mLeafNodes.size() * 5; }
765 
766 private:
767  std::vector<LeafNodeType*> mLeafNodes;
768  boost::scoped_array<size_t> mOffsets;
769 }; // struct LeafNodeConnectivityTable
770 
771 
772 template<typename TreeType>
774 {
775 public:
776 
777  enum Axis { X_AXIS = 0, Y_AXIS = 1, Z_AXIS = 2 };
778 
779  typedef typename TreeType::ValueType ValueType;
780  typedef typename TreeType::LeafNodeType LeafNodeType;
782 
783  SweepExteriorSign(Axis axis, const std::vector<size_t>& startNodeIndices, ConnectivityTable& connectivity)
784  : mStartNodeIndices(startNodeIndices.empty() ? NULL : &startNodeIndices[0])
785  , mConnectivity(&connectivity)
786  , mAxis(axis)
787  {
788  }
789 
790  void operator()(const tbb::blocked_range<size_t>& range) const {
791 
792  std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
793 
794  // Z Axis
795  size_t idxA = 0, idxB = 1;
796  Index step = 1;
797 
798  const size_t* nextOffsets = mConnectivity->offsetsNextZ();
799  const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
800 
801  if (mAxis == Y_AXIS) {
802 
803  idxA = 0;
804  idxB = 2;
805  step = LeafNodeType::DIM;
806 
807  nextOffsets = mConnectivity->offsetsNextY();
808  prevOffsets = mConnectivity->offsetsPrevY();
809 
810  } else if (mAxis == X_AXIS) {
811 
812  idxA = 1;
813  idxB = 2;
814  step = LeafNodeType::DIM * LeafNodeType::DIM;
815 
816  nextOffsets = mConnectivity->offsetsNextX();
817  prevOffsets = mConnectivity->offsetsPrevX();
818  }
819 
820  Coord ijk(0, 0, 0);
821 
822  int& a = ijk[idxA];
823  int& b = ijk[idxB];
824 
825  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
826 
827  size_t startOffset = mStartNodeIndices[n];
828  size_t lastOffset = startOffset;
829 
830  Index pos(0);
831 
832  for (a = 0; a < int(LeafNodeType::DIM); ++a) {
833  for (b = 0; b < int(LeafNodeType::DIM); ++b) {
834 
835  pos = LeafNodeType::coordToOffset(ijk);
836  size_t offset = startOffset;
837 
838  // sweep in +axis direction until a boundary voxel is hit.
839  while ( offset != ConnectivityTable::INVALID_OFFSET &&
840  traceVoxelLine(*nodes[offset], pos, step) ) {
841 
842  lastOffset = offset;
843  offset = nextOffsets[offset];
844  }
845 
846  // find last leafnode in +axis direction
847  offset = lastOffset;
848  while (offset != ConnectivityTable::INVALID_OFFSET) {
849  lastOffset = offset;
850  offset = nextOffsets[offset];
851  }
852 
853  // sweep in -axis direction until a boundary voxel is hit.
854  offset = lastOffset;
855  pos += step * (LeafNodeType::DIM - 1);
856  while ( offset != ConnectivityTable::INVALID_OFFSET &&
857  traceVoxelLine(*nodes[offset], pos, -step)) {
858  offset = prevOffsets[offset];
859  }
860  }
861  }
862  }
863  }
864 
865 
866  bool traceVoxelLine(LeafNodeType& node, Index pos, Index step) const {
867 
868  ValueType* data = node.buffer().data();
869 
870  bool isOutside = true;
871 
872  for (Index i = 0; i < LeafNodeType::DIM; ++i) {
873 
874  ValueType& dist = data[pos];
875 
876  if (dist < ValueType(0.0)) {
877  isOutside = true;
878  } else {
879  // Boundary voxel check. (Voxel that intersects the surface)
880  if (!(dist > ValueType(0.75))) isOutside = false;
881 
882  if (isOutside) dist = ValueType(-dist);
883  }
884 
885  pos += step;
886  }
887 
888  return isOutside;
889  }
890 
891 
892 private:
893  size_t const * const mStartNodeIndices;
894  ConnectivityTable * const mConnectivity;
895 
896  const Axis mAxis;
897 }; // class SweepExteriorSign
898 
899 
900 template<typename LeafNodeType>
901 inline void
902 seedFill(LeafNodeType& node)
903 {
904  typedef typename LeafNodeType::ValueType ValueType;
905  typedef std::deque<Index> Queue;
906 
907 
908  ValueType* data = node.buffer().data();
909 
910  // find seed points
911  Queue seedPoints;
912  for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
913  if (data[pos] < 0.0) seedPoints.push_back(pos);
914  }
915 
916  if (seedPoints.empty()) return;
917 
918  // clear sign information
919  for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
920  ValueType& dist = data[*it];
921  dist = -dist;
922  }
923 
924  // flood fill
925 
926  Coord ijk(0, 0, 0);
927  Index pos(0), nextPos(0);
928 
929  while (!seedPoints.empty()) {
930 
931  pos = seedPoints.back();
932  seedPoints.pop_back();
933 
934  ValueType& dist = data[pos];
935 
936  if (!(dist < ValueType(0.0))) {
937 
938  dist = -dist; // flip sign
939 
940  ijk = LeafNodeType::offsetToLocalCoord(pos);
941 
942  if (ijk[0] != 0) { // i - 1, j, k
943  nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
944  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
945  }
946 
947  if (ijk[0] != (LeafNodeType::DIM - 1)) { // i + 1, j, k
948  nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
949  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
950  }
951 
952  if (ijk[1] != 0) { // i, j - 1, k
953  nextPos = pos - LeafNodeType::DIM;
954  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
955  }
956 
957  if (ijk[1] != (LeafNodeType::DIM - 1)) { // i, j + 1, k
958  nextPos = pos + LeafNodeType::DIM;
959  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
960  }
961 
962  if (ijk[2] != 0) { // i, j, k - 1
963  nextPos = pos - 1;
964  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
965  }
966 
967  if (ijk[2] != (LeafNodeType::DIM - 1)) { // i, j, k + 1
968  nextPos = pos + 1;
969  if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
970  }
971  }
972  }
973 } // seedFill()
974 
975 
976 template<typename LeafNodeType>
977 inline bool
978 scanFill(LeafNodeType& node)
979 {
980  bool updatedNode = false;
981 
982  typedef typename LeafNodeType::ValueType ValueType;
983  ValueType* data = node.buffer().data();
984 
985  Coord ijk(0, 0, 0);
986 
987  bool updatedSign = true;
988  while (updatedSign) {
989 
990  updatedSign = false;
991 
992  for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
993 
994  ValueType& dist = data[pos];
995 
996  if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
997 
998  ijk = LeafNodeType::offsetToLocalCoord(pos);
999 
1000  // i, j, k - 1
1001  if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1002  updatedSign = true;
1003  dist = ValueType(-dist);
1004 
1005  // i, j, k + 1
1006  } else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1007  updatedSign = true;
1008  dist = ValueType(-dist);
1009 
1010  // i, j - 1, k
1011  } else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1012  updatedSign = true;
1013  dist = ValueType(-dist);
1014 
1015  // i, j + 1, k
1016  } else if (ijk[1] != (LeafNodeType::DIM - 1) && data[pos + LeafNodeType::DIM] < ValueType(0.0)) {
1017  updatedSign = true;
1018  dist = ValueType(-dist);
1019 
1020  // i - 1, j, k
1021  } else if (ijk[0] != 0 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0)) {
1022  updatedSign = true;
1023  dist = ValueType(-dist);
1024 
1025  // i + 1, j, k
1026  } else if (ijk[0] != (LeafNodeType::DIM - 1) && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0)) {
1027  updatedSign = true;
1028  dist = ValueType(-dist);
1029  }
1030  }
1031  } // end value loop
1032 
1033  updatedNode |= updatedSign;
1034  } // end update loop
1035 
1036  return updatedNode;
1037 } // scanFill()
1038 
1039 
1040 template<typename TreeType>
1042 {
1043 public:
1044  typedef typename TreeType::ValueType ValueType;
1045  typedef typename TreeType::LeafNodeType LeafNodeType;
1046 
1047  SeedFillExteriorSign(std::vector<LeafNodeType*>& nodes, bool* changedNodeMask)
1048  : mNodes(nodes.empty() ? NULL : &nodes[0])
1049  , mChangedNodeMask(changedNodeMask)
1050  {
1051  }
1052 
1053  void operator()(const tbb::blocked_range<size_t>& range) const {
1054  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1055  if (mChangedNodeMask[n]) {
1056  //seedFill(*mNodes[n]);
1057  mChangedNodeMask[n] = scanFill(*mNodes[n]);
1058  }
1059  }
1060  }
1061 
1062  LeafNodeType ** const mNodes;
1063  bool * const mChangedNodeMask;
1064 };
1065 
1066 
1067 template<typename ValueType>
1069 {
1070  FillArray(ValueType* array, const ValueType v) : mArray(array), mValue(v) { }
1071 
1072  void operator()(const tbb::blocked_range<size_t>& range) const {
1073  const ValueType v = mValue;
1074  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1075  mArray[n] = v;
1076  }
1077  }
1078 
1079  ValueType * const mArray;
1080  const ValueType mValue;
1081 };
1082 
1083 
1084 template<typename ValueType>
1085 inline void
1086 fillArray(ValueType* array, const ValueType val, const size_t length)
1087 {
1088  const size_t grainSize = length / tbb::task_scheduler_init::default_num_threads();
1089  const tbb::blocked_range<size_t> range(0, length, grainSize);
1090  tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1091 }
1092 
1093 
1094 template<typename TreeType>
1096 {
1097 public:
1098  typedef typename TreeType::ValueType ValueType;
1099  typedef typename TreeType::LeafNodeType LeafNodeType;
1100 
1101  SyncVoxelMask(std::vector<LeafNodeType*>& nodes, const bool* changedNodeMask, bool* changedVoxelMask)
1102  : mNodes(nodes.empty() ? NULL : &nodes[0])
1103  , mChangedNodeMask(changedNodeMask)
1104  , mChangedVoxelMask(changedVoxelMask)
1105  {
1106  }
1107 
1108  void operator()(const tbb::blocked_range<size_t>& range) const {
1109  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1110 
1111  if (mChangedNodeMask[n]) {
1112  bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1113 
1114  ValueType* data = mNodes[n]->buffer().data();
1115 
1116  for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1117  if (mask[pos]) {
1118  data[pos] = ValueType(-data[pos]);
1119  mask[pos] = false;
1120  }
1121  }
1122  }
1123  }
1124  }
1125 
1126  LeafNodeType ** const mNodes;
1127  bool const * const mChangedNodeMask;
1128  bool * const mChangedVoxelMask;
1129 };
1130 
1131 
1132 template<typename TreeType>
1134 {
1135 public:
1136  typedef typename TreeType::ValueType ValueType;
1137  typedef typename TreeType::LeafNodeType LeafNodeType;
1139 
1140  SeedPoints(ConnectivityTable& connectivity, bool* changedNodeMask, bool* nodeMask, bool* changedVoxelMask)
1141  : mConnectivity(&connectivity)
1142  , mChangedNodeMask(changedNodeMask)
1143  , mNodeMask(nodeMask)
1144  , mChangedVoxelMask(changedVoxelMask)
1145  {
1146  }
1147 
1148  void operator()(const tbb::blocked_range<size_t>& range) const {
1149 
1150  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1151 
1152  if (!mChangedNodeMask[n]) {
1153 
1154  bool changedValue = false;
1155 
1156  changedValue |= processZ(n, /*firstFace=*/true);
1157  changedValue |= processZ(n, /*firstFace=*/false);
1158 
1159  changedValue |= processY(n, /*firstFace=*/true);
1160  changedValue |= processY(n, /*firstFace=*/false);
1161 
1162  changedValue |= processX(n, /*firstFace=*/true);
1163  changedValue |= processX(n, /*firstFace=*/false);
1164 
1165  mNodeMask[n] = changedValue;
1166  }
1167  }
1168  }
1169 
1170 
1171  bool processZ(const size_t n, bool firstFace) const
1172  {
1173  const size_t offset = firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1174  if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1175 
1176  bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1177 
1178  const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1179  const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1180 
1181  const Index lastOffset = LeafNodeType::DIM - 1;
1182  const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1183 
1184  Index tmpPos(0), pos(0);
1185  bool changedValue = false;
1186 
1187  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
1188  tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1189  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
1190  pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1191 
1192  if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1193  if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1194  changedValue = true;
1195  mask[pos + lhsOffset] = true;
1196  }
1197  }
1198  }
1199  }
1200 
1201  return changedValue;
1202  }
1203 
1204  return false;
1205  }
1206 
1207  bool processY(const size_t n, bool firstFace) const
1208  {
1209  const size_t offset = firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1210  if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1211 
1212  bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1213 
1214  const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1215  const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1216 
1217  const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1218  const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1219 
1220  Index tmpPos(0), pos(0);
1221  bool changedValue = false;
1222 
1223  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
1224  tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1225  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
1226  pos = tmpPos + z;
1227 
1228  if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1229  if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1230  changedValue = true;
1231  mask[pos + lhsOffset] = true;
1232  }
1233  }
1234  }
1235  }
1236 
1237  return changedValue;
1238  }
1239 
1240  return false;
1241  }
1242 
1243  bool processX(const size_t n, bool firstFace) const
1244  {
1245  const size_t offset = firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1246  if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1247 
1248  bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1249 
1250  const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1251  const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1252 
1253  const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1254  const Index lhsOffset = firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1255 
1256  Index tmpPos(0), pos(0);
1257  bool changedValue = false;
1258 
1259  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
1260  tmpPos = y << LeafNodeType::LOG2DIM;
1261  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
1262  pos = tmpPos + z;
1263 
1264  if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1265  if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1266  changedValue = true;
1267  mask[pos + lhsOffset] = true;
1268  }
1269  }
1270  }
1271  }
1272 
1273  return changedValue;
1274  }
1275 
1276  return false;
1277  }
1278 
1279  ConnectivityTable * const mConnectivity;
1280  bool * const mChangedNodeMask;
1281  bool * const mNodeMask;
1282  bool * const mChangedVoxelMask;
1283 };
1284 
1285 
1287 
1288 template<typename TreeType, typename MeshDataAdapter>
1290 {
1291  typedef typename TreeType::ValueType ValueType;
1292  typedef typename TreeType::LeafNodeType LeafNodeType;
1293  typedef typename TreeType::template ValueConverter<Int32>::Type Int32TreeType;
1294  typedef typename Int32TreeType::LeafNodeType Int32LeafNodeType;
1295 
1296  typedef std::pair<boost::shared_array<Vec3d>, boost::shared_array<bool> > LocalData;
1297  typedef tbb::enumerable_thread_specific<LocalData> LocalDataTable;
1298 
1300  std::vector<LeafNodeType*>& distNodes,
1301  const TreeType& distTree,
1302  const Int32TreeType& indexTree,
1303  const MeshDataAdapter& mesh)
1304  : mDistNodes(distNodes.empty() ? NULL : &distNodes[0])
1305  , mDistTree(&distTree)
1306  , mIndexTree(&indexTree)
1307  , mMesh(&mesh)
1308  , mLocalDataTable(new LocalDataTable())
1309  {
1310  }
1311 
1312 
1313  void operator()(const tbb::blocked_range<size_t>& range) const {
1314 
1315  tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1316  tree::ValueAccessor<const Int32TreeType> idxAcc(*mIndexTree);
1317 
1318  ValueType nval;
1319  CoordBBox bbox;
1320  Index xPos(0), yPos(0);
1321  Coord ijk, nijk, nodeMin, nodeMax;
1322  Vec3d cp, xyz, nxyz, dir1, dir2;
1323 
1324  LocalData& localData = mLocalDataTable->local();
1325 
1326  boost::shared_array<Vec3d>& points = localData.first;
1327  if (!points) points.reset(new Vec3d[LeafNodeType::SIZE * 2]);
1328 
1329  boost::shared_array<bool>& mask = localData.second;
1330  if (!mask) mask.reset(new bool[LeafNodeType::SIZE]);
1331 
1332 
1333  typename LeafNodeType::ValueOnCIter it;
1334 
1335  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1336 
1337  LeafNodeType& node = *mDistNodes[n];
1338  ValueType* data = node.buffer().data();
1339 
1340  const Int32LeafNodeType* idxNode = idxAcc.probeConstLeaf(node.origin());
1341  const Int32* idxData = idxNode->buffer().data();
1342 
1343  nodeMin = node.origin();
1344  nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1345 
1346  // reset computed voxel mask.
1347  memset(mask.get(), 0, sizeof(bool) * LeafNodeType::SIZE);
1348 
1349  for (it = node.cbeginValueOn(); it; ++it) {
1350  Index pos = it.pos();
1351 
1352  ValueType& dist = data[pos];
1353  if (dist < 0.0 || dist > 0.75) continue;
1354 
1355  ijk = node.offsetToGlobalCoord(pos);
1356 
1357  xyz[0] = double(ijk[0]);
1358  xyz[1] = double(ijk[1]);
1359  xyz[2] = double(ijk[2]);
1360 
1361 
1362  bbox.min() = Coord::maxComponent(ijk.offsetBy(-1), nodeMin);
1363  bbox.max() = Coord::minComponent(ijk.offsetBy(1), nodeMax);
1364 
1365  bool flipSign = false;
1366 
1367  for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1368  xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1369  for (nijk[1] = bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1370  yPos = xPos + ((nijk[1] & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
1371  for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1372  pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1373 
1374  const Int32& polyIdx = idxData[pos];
1375 
1376  if (polyIdx == Int32(util::INVALID_IDX) || !(data[pos] < -0.75)) continue;
1377 
1378  const Index pointIndex = pos * 2;
1379 
1380  if (!mask[pos]) {
1381 
1382  mask[pos] = true;
1383 
1384  nxyz[0] = double(nijk[0]);
1385  nxyz[1] = double(nijk[1]);
1386  nxyz[2] = double(nijk[2]);
1387 
1388  Vec3d& point = points[pointIndex];
1389 
1390  point = closestPoint(nxyz, polyIdx);
1391 
1392  Vec3d& direction = points[pointIndex + 1];
1393  direction = nxyz - point;
1394  direction.normalize();
1395  }
1396 
1397  dir1 = xyz - points[pointIndex];
1398  dir1.normalize();
1399 
1400  if (points[pointIndex + 1].dot(dir1) > 0.0) {
1401  flipSign = true;
1402  break;
1403  }
1404  }
1405  }
1406  }
1407 
1408  if (flipSign) {
1409  dist = -dist;
1410  } else {
1411  for (Int32 m = 0; m < 26; ++m) {
1412  nijk = ijk + util::COORD_OFFSETS[m];
1413 
1414  if (!bbox.isInside(nijk) && distAcc.probeValue(nijk, nval) && nval < -0.75) {
1415  nxyz[0] = double(nijk[0]);
1416  nxyz[1] = double(nijk[1]);
1417  nxyz[2] = double(nijk[2]);
1418 
1419  cp = closestPoint(nxyz, idxAcc.getValue(nijk));
1420 
1421  dir1 = xyz - cp;
1422  dir1.normalize();
1423 
1424  dir2 = nxyz - cp;
1425  dir2.normalize();
1426 
1427  if (dir2.dot(dir1) > 0.0) {
1428  dist = -dist;
1429  break;
1430  }
1431  }
1432  }
1433  }
1434 
1435  } // active voxel loop
1436  } // leaf node loop
1437  }
1438 
1439 private:
1440 
1441  Vec3d closestPoint(const Vec3d& center, Int32 polyIdx) const
1442  {
1443  Vec3d a, b, c, cp, uvw;
1444 
1445  const size_t polygon = size_t(polyIdx);
1446  mMesh->getIndexSpacePoint(polygon, 0, a);
1447  mMesh->getIndexSpacePoint(polygon, 1, b);
1448  mMesh->getIndexSpacePoint(polygon, 2, c);
1449 
1450  cp = closestPointOnTriangleToPoint(a, c, b, center, uvw);
1451 
1452  if (4 == mMesh->vertexCount(polygon)) {
1453 
1454  mMesh->getIndexSpacePoint(polygon, 3, b);
1455 
1456  c = closestPointOnTriangleToPoint(a, b, c, center, uvw);
1457 
1458  if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1459  cp = c;
1460  }
1461  }
1462 
1463  return cp;
1464  }
1465 
1466 
1467  LeafNodeType ** const mDistNodes;
1468  TreeType const * const mDistTree;
1469  Int32TreeType const * const mIndexTree;
1470  MeshDataAdapter const * const mMesh;
1471 
1472  boost::shared_ptr<LocalDataTable> mLocalDataTable;
1473 }; // ComputeIntersectingVoxelSign
1474 
1475 
1477 
1478 
1479 template<typename LeafNodeType>
1480 inline void
1481 maskNodeInternalNeighbours(const Index pos, bool (&mask)[26])
1482 {
1483  typedef LeafNodeType NodeT;
1484 
1485  const Coord ijk = NodeT::offsetToLocalCoord(pos);
1486 
1487  // Face adjacent neighbours
1488  // i+1, j, k
1489  mask[0] = ijk[0] != (NodeT::DIM - 1);
1490  // i-1, j, k
1491  mask[1] = ijk[0] != 0;
1492  // i, j+1, k
1493  mask[2] = ijk[1] != (NodeT::DIM - 1);
1494  // i, j-1, k
1495  mask[3] = ijk[1] != 0;
1496  // i, j, k+1
1497  mask[4] = ijk[2] != (NodeT::DIM - 1);
1498  // i, j, k-1
1499  mask[5] = ijk[2] != 0;
1500 
1501  // Edge adjacent neighbour
1502  // i+1, j, k-1
1503  mask[6] = mask[0] && mask[5];
1504  // i-1, j, k-1
1505  mask[7] = mask[1] && mask[5];
1506  // i+1, j, k+1
1507  mask[8] = mask[0] && mask[4];
1508  // i-1, j, k+1
1509  mask[9] = mask[1] && mask[4];
1510  // i+1, j+1, k
1511  mask[10] = mask[0] && mask[2];
1512  // i-1, j+1, k
1513  mask[11] = mask[1] && mask[2];
1514  // i+1, j-1, k
1515  mask[12] = mask[0] && mask[3];
1516  // i-1, j-1, k
1517  mask[13] = mask[1] && mask[3];
1518  // i, j-1, k+1
1519  mask[14] = mask[3] && mask[4];
1520  // i, j-1, k-1
1521  mask[15] = mask[3] && mask[5];
1522  // i, j+1, k+1
1523  mask[16] = mask[2] && mask[4];
1524  // i, j+1, k-1
1525  mask[17] = mask[2] && mask[5];
1526 
1527  // Corner adjacent neighbours
1528  // i-1, j-1, k-1
1529  mask[18] = mask[1] && mask[3] && mask[5];
1530  // i-1, j-1, k+1
1531  mask[19] = mask[1] && mask[3] && mask[4];
1532  // i+1, j-1, k+1
1533  mask[20] = mask[0] && mask[3] && mask[4];
1534  // i+1, j-1, k-1
1535  mask[21] = mask[0] && mask[3] && mask[5];
1536  // i-1, j+1, k-1
1537  mask[22] = mask[1] && mask[2] && mask[5];
1538  // i-1, j+1, k+1
1539  mask[23] = mask[1] && mask[2] && mask[4];
1540  // i+1, j+1, k+1
1541  mask[24] = mask[0] && mask[2] && mask[4];
1542  // i+1, j+1, k-1
1543  mask[25] = mask[0] && mask[2] && mask[5];
1544 }
1545 
1546 
1547 template<typename Compare, typename LeafNodeType>
1548 inline bool
1549 checkNeighbours(const Index pos, const typename LeafNodeType::ValueType * data, bool (&mask)[26])
1550 {
1551  typedef LeafNodeType NodeT;
1552 
1553  // i, j, k - 1
1554  if (mask[5] && Compare::check(data[pos - 1])) return true;
1555  // i, j, k + 1
1556  if (mask[4] && Compare::check(data[pos + 1])) return true;
1557  // i, j - 1, k
1558  if (mask[3] && Compare::check(data[pos - NodeT::DIM])) return true;
1559  // i, j + 1, k
1560  if (mask[2] && Compare::check(data[pos + NodeT::DIM])) return true;
1561  // i - 1, j, k
1562  if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM])) return true;
1563  // i + 1, j, k
1564  if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM])) return true;
1565  // i+1, j, k-1
1566  if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM])) return true;
1567  // i-1, j, k-1
1568  if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1])) return true;
1569  // i+1, j, k+1
1570  if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1])) return true;
1571  // i-1, j, k+1
1572  if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1])) return true;
1573  // i+1, j+1, k
1574  if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM])) return true;
1575  // i-1, j+1, k
1576  if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM])) return true;
1577  // i+1, j-1, k
1578  if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM])) return true;
1579  // i-1, j-1, k
1580  if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM])) return true;
1581  // i, j-1, k+1
1582  if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1])) return true;
1583  // i, j-1, k-1
1584  if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1])) return true;
1585  // i, j+1, k+1
1586  if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1])) return true;
1587  // i, j+1, k-1
1588  if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1])) return true;
1589  // i-1, j-1, k-1
1590  if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1])) return true;
1591  // i-1, j-1, k+1
1592  if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1])) return true;
1593  // i+1, j-1, k+1
1594  if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1])) return true;
1595  // i+1, j-1, k-1
1596  if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1])) return true;
1597  // i-1, j+1, k-1
1598  if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1])) return true;
1599  // i-1, j+1, k+1
1600  if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1])) return true;
1601  // i+1, j+1, k+1
1602  if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1])) return true;
1603  // i+1, j+1, k-1
1604  if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1])) return true;
1605 
1606  return false;
1607 }
1608 
1609 
1610 template<typename Compare, typename AccessorType>
1611 inline bool
1612 checkNeighbours(const Coord& ijk, AccessorType& acc, bool (&mask)[26])
1613 {
1614  for (Int32 m = 0; m < 26; ++m) {
1615  if (!mask[m] && Compare::check(acc.getValue(ijk + util::COORD_OFFSETS[m]))) {
1616  return true;
1617  }
1618  }
1619 
1620  return false;
1621 }
1622 
1623 
1624 template<typename TreeType>
1626 {
1627  typedef typename TreeType::ValueType ValueType;
1628  typedef typename TreeType::LeafNodeType LeafNodeType;
1629 
1630  struct IsNegative { static bool check(const ValueType v) { return v < ValueType(0.0); } };
1631 
1632  ValidateIntersectingVoxels(TreeType& tree, std::vector<LeafNodeType*>& nodes)
1633  : mTree(&tree)
1634  , mNodes(nodes.empty() ? NULL : &nodes[0])
1635  {
1636  }
1637 
1638  void operator()(const tbb::blocked_range<size_t>& range) const
1639  {
1641  bool neighbourMask[26];
1642 
1643  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1644 
1645  LeafNodeType& node = *mNodes[n];
1646  ValueType* data = node.buffer().data();
1647 
1648  typename LeafNodeType::ValueOnCIter it;
1649  for (it = node.cbeginValueOn(); it; ++it) {
1650 
1651  const Index pos = it.pos();
1652 
1653  ValueType& dist = data[pos];
1654  if (dist < 0.0 || dist > 0.75) continue;
1655 
1656  // Mask node internal neighbours
1657  maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1658 
1659  const bool hasNegativeNeighbour =
1660  checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1661  checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1662 
1663  if (!hasNegativeNeighbour) {
1664  // push over boundary voxel distance
1665  dist = ValueType(0.75) + Tolerance<ValueType>::epsilon();
1666  }
1667  }
1668  }
1669  }
1670 
1671  TreeType * const mTree;
1672  LeafNodeType ** const mNodes;
1673 }; // ValidateIntersectingVoxels
1674 
1675 
1676 template<typename TreeType>
1678 {
1679  typedef typename TreeType::ValueType ValueType;
1680  typedef typename TreeType::LeafNodeType LeafNodeType;
1681  typedef typename TreeType::template ValueConverter<Int32>::Type Int32TreeType;
1682 
1683  struct Comp { static bool check(const ValueType v) { return !(v > ValueType(0.75)); } };
1684 
1685  RemoveSelfIntersectingSurface(std::vector<LeafNodeType*>& nodes,
1686  TreeType& distTree, Int32TreeType& indexTree)
1687  : mNodes(nodes.empty() ? NULL : &nodes[0])
1688  , mDistTree(&distTree)
1689  , mIndexTree(&indexTree)
1690  {
1691  }
1692 
1693  void operator()(const tbb::blocked_range<size_t>& range) const
1694  {
1695  tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1696  tree::ValueAccessor<Int32TreeType> idxAcc(*mIndexTree);
1697  bool neighbourMask[26];
1698 
1699  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1700 
1701  LeafNodeType& distNode = *mNodes[n];
1702  ValueType* data = distNode.buffer().data();
1703 
1704  typename Int32TreeType::LeafNodeType* idxNode =
1705  idxAcc.probeLeaf(distNode.origin());
1706 
1707  typename LeafNodeType::ValueOnCIter it;
1708  for (it = distNode.cbeginValueOn(); it; ++it) {
1709 
1710  const Index pos = it.pos();
1711 
1712  if (!(data[pos] > 0.75)) continue;
1713 
1714  // Mask node internal neighbours
1715  maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1716 
1717  const bool hasBoundaryNeighbour =
1718  checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1719  checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos), distAcc, neighbourMask);
1720 
1721  if (!hasBoundaryNeighbour) {
1722  distNode.setValueOff(pos);
1723  idxNode->setValueOff(pos);
1724  }
1725  }
1726  }
1727  }
1728 
1729  LeafNodeType * * const mNodes;
1730  TreeType * const mDistTree;
1731  Int32TreeType * const mIndexTree;
1732 }; // RemoveSelfIntersectingSurface
1733 
1734 
1736 
1737 
1738 template<typename NodeType>
1740 {
1741  ReleaseChildNodes(NodeType ** nodes) : mNodes(nodes) {}
1742 
1743  void operator()(const tbb::blocked_range<size_t>& range) const {
1744 
1745  typedef typename NodeType::NodeMaskType NodeMaskType;
1746 
1747  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1748  const_cast<NodeMaskType&>(mNodes[n]->getChildMask()).setOff();
1749  }
1750  }
1751 
1752  NodeType ** const mNodes;
1753 };
1754 
1755 
1756 template<typename TreeType>
1757 inline void
1758 releaseLeafNodes(TreeType& tree)
1759 {
1760  typedef typename TreeType::RootNodeType RootNodeType;
1761  typedef typename RootNodeType::NodeChainType NodeChainType;
1762  typedef typename boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
1763 
1764  std::vector<InternalNodeType*> nodes;
1765  tree.getNodes(nodes);
1766 
1767  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1768  ReleaseChildNodes<InternalNodeType>(nodes.empty() ? NULL : &nodes[0]));
1769 }
1770 
1771 
1772 template<typename TreeType>
1774 {
1775  typedef typename TreeType::LeafNodeType LeafNodeType;
1776 
1777  StealUniqueLeafNodes(TreeType& lhsTree, TreeType& rhsTree,
1778  std::vector<LeafNodeType*>& overlappingNodes)
1779  : mLhsTree(&lhsTree)
1780  , mRhsTree(&rhsTree)
1781  , mNodes(&overlappingNodes)
1782  {
1783  }
1784 
1785  void operator()() const {
1786 
1787  std::vector<LeafNodeType*> rhsLeafNodes;
1788 
1789  rhsLeafNodes.reserve(mRhsTree->leafCount());
1790  //mRhsTree->getNodes(rhsLeafNodes);
1791  //releaseLeafNodes(*mRhsTree);
1792  mRhsTree->stealNodes(rhsLeafNodes);
1793 
1794  tree::ValueAccessor<TreeType> acc(*mLhsTree);
1795 
1796  for (size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1797  if (!acc.probeLeaf(rhsLeafNodes[n]->origin())) {
1798  acc.addLeaf(rhsLeafNodes[n]);
1799  } else {
1800  mNodes->push_back(rhsLeafNodes[n]);
1801  }
1802  }
1803  }
1804 
1805 private:
1806  TreeType * const mLhsTree;
1807  TreeType * const mRhsTree;
1808  std::vector<LeafNodeType*> * const mNodes;
1809 };
1810 
1811 
1812 template<typename DistTreeType, typename IndexTreeType>
1813 inline void
1814 combineData(DistTreeType& lhsDist, IndexTreeType& lhsIdx,
1815  DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1816 {
1817  typedef typename DistTreeType::LeafNodeType DistLeafNodeType;
1818  typedef typename IndexTreeType::LeafNodeType IndexLeafNodeType;
1819 
1820  std::vector<DistLeafNodeType*> overlappingDistNodes;
1821  std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1822 
1823  // Steal unique leafnodes
1824  tbb::task_group tasks;
1825  tasks.run(StealUniqueLeafNodes<DistTreeType>(lhsDist, rhsDist, overlappingDistNodes));
1826  tasks.run(StealUniqueLeafNodes<IndexTreeType>(lhsIdx, rhsIdx, overlappingIdxNodes));
1827  tasks.wait();
1828 
1829  // Combine overlapping leaf nodes
1830  if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1831  tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1832  CombineLeafNodes<DistTreeType>(lhsDist, lhsIdx, &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1833  }
1834 }
1835 
1836 
1843 template<typename TreeType>
1845 
1846  typedef boost::scoped_ptr<VoxelizationData> Ptr;
1847  typedef typename TreeType::ValueType ValueType;
1848 
1849  typedef typename TreeType::template ValueConverter<Int32>::Type Int32TreeType;
1850  typedef typename TreeType::template ValueConverter<unsigned char>::Type UCharTreeType;
1851 
1855 
1856 
1858  : distTree(std::numeric_limits<ValueType>::max())
1859  , distAcc(distTree)
1860  , indexTree(Int32(util::INVALID_IDX))
1861  , indexAcc(indexTree)
1862  , primIdTree(MaxPrimId)
1863  , primIdAcc(primIdTree)
1864  , mPrimCount(0)
1865  {
1866  }
1867 
1868  TreeType distTree;
1869  FloatTreeAcc distAcc;
1870 
1871  Int32TreeType indexTree;
1872  Int32TreeAcc indexAcc;
1873 
1874  UCharTreeType primIdTree;
1875  UCharTreeAcc primIdAcc;
1876 
1877  unsigned char getNewPrimId() {
1878 
1879  if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1880  mPrimCount = 0;
1881  primIdTree.clear();
1882  }
1883 
1884  return mPrimCount++;
1885  }
1886 
1887 private:
1888 
1889  enum { MaxPrimId = 100 };
1890 
1891  unsigned char mPrimCount;
1892 };
1893 
1894 
1895 template<typename TreeType, typename MeshDataAdapter, typename Interrupter = util::NullInterrupter>
1897 {
1898 public:
1899 
1901  typedef tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr> DataTable;
1902 
1903  VoxelizePolygons(DataTable& dataTable,
1904  const MeshDataAdapter& mesh,
1905  Interrupter* interrupter = NULL)
1906  : mDataTable(&dataTable)
1907  , mMesh(&mesh)
1908  , mInterrupter(interrupter)
1909  {
1910  }
1911 
1912  void operator()(const tbb::blocked_range<size_t>& range) const {
1913 
1914  typename VoxelizationDataType::Ptr& dataPtr = mDataTable->local();
1915  if (!dataPtr) dataPtr.reset(new VoxelizationDataType());
1916 
1917  Triangle prim;
1918 
1919  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1920 
1921  if (this->wasInterrupted()) {
1922  tbb::task::self().cancel_group_execution();
1923  break;
1924  }
1925 
1926  const size_t numVerts = mMesh->vertexCount(n);
1927 
1928  // rasterize triangles and quads.
1929  if (numVerts == 3 || numVerts == 4) {
1930 
1931  prim.index = Int32(n);
1932 
1933  mMesh->getIndexSpacePoint(n, 0, prim.a);
1934  mMesh->getIndexSpacePoint(n, 1, prim.b);
1935  mMesh->getIndexSpacePoint(n, 2, prim.c);
1936 
1937  evalTriangle(prim, *dataPtr);
1938 
1939  if (numVerts == 4) {
1940  mMesh->getIndexSpacePoint(n, 3, prim.b);
1941  evalTriangle(prim, *dataPtr);
1942  }
1943  }
1944  }
1945  }
1946 
1947 private:
1948 
1949  bool wasInterrupted() const { return mInterrupter && mInterrupter->wasInterrupted(); }
1950 
1951  struct Triangle { Vec3d a, b, c; Int32 index; };
1952 
1953  struct SubTask
1954  {
1955  SubTask(const Triangle& prim, DataTable& dataTable, size_t polygonCount)
1956  : mLocalDataTable(&dataTable)
1957  , mPrim(prim)
1958  , mPolygonCount(polygonCount)
1959  {
1960  }
1961 
1962  void operator()() const
1963  {
1964  const size_t minNumTask = size_t(tbb::task_scheduler_init::default_num_threads() * 10);
1965 
1966  if (mPolygonCount > minNumTask) {
1967 
1968  typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
1969  if (!dataPtr) dataPtr.reset(new VoxelizationDataType());
1970 
1971  voxelizeTriangle(mPrim, *dataPtr);
1972 
1973  } else {
1974  spawnTasks(mPrim, *mLocalDataTable, mPolygonCount);
1975  }
1976  }
1977 
1978  DataTable * const mLocalDataTable;
1979  const Triangle mPrim;
1980  const size_t mPolygonCount;
1981  };
1982 
1983 
1984  void evalTriangle(const Triangle& prim, VoxelizationDataType& data) const
1985  {
1986  const size_t minNumTask = size_t(tbb::task_scheduler_init::default_num_threads() * 10);
1987 
1988  if (mMesh->polygonCount() > minNumTask) {
1989 
1990  voxelizeTriangle(prim, data);
1991 
1992  } else {
1993 
1994  spawnTasks(prim, *mDataTable, mMesh->polygonCount());
1995  }
1996  }
1997 
1998  static void spawnTasks(const Triangle& mainPrim, DataTable& dataTable, size_t primCount)
1999  {
2000  const size_t newPrimCount = primCount * 4;
2001 
2002  tbb::task_group tasks;
2003 
2004  const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2005  const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2006  const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2007 
2008  Triangle prim;
2009  prim.index = mainPrim.index;
2010 
2011  prim.a = mainPrim.a;
2012  prim.b = ab;
2013  prim.c = ac;
2014  tasks.run(SubTask(prim, dataTable, newPrimCount));
2015 
2016  prim.a = ab;
2017  prim.b = bc;
2018  prim.c = ac;
2019  tasks.run(SubTask(prim, dataTable, newPrimCount));
2020 
2021  prim.a = ab;
2022  prim.b = mainPrim.b;
2023  prim.c = bc;
2024  tasks.run(SubTask(prim, dataTable, newPrimCount));
2025 
2026  prim.a = ac;
2027  prim.b = bc;
2028  prim.c = mainPrim.c;
2029  tasks.run(SubTask(prim, dataTable, newPrimCount));
2030 
2031  tasks.wait();
2032  }
2033 
2034  static void voxelizeTriangle(const Triangle& prim, VoxelizationDataType& data)
2035  {
2036  std::deque<Coord> coordList;
2037  Coord ijk, nijk;
2038 
2039  ijk = Coord::floor(prim.a);
2040  coordList.push_back(ijk);
2041 
2042  computeDistance(ijk, prim, data);
2043 
2044  unsigned char primId = data.getNewPrimId();
2045  data.primIdAcc.setValueOnly(ijk, primId);
2046 
2047  while (!coordList.empty()) {
2048  ijk = coordList.back();
2049  coordList.pop_back();
2050 
2051  for (Int32 i = 0; i < 26; ++i) {
2052  nijk = ijk + util::COORD_OFFSETS[i];
2053  if (primId != data.primIdAcc.getValue(nijk)) {
2054  data.primIdAcc.setValueOnly(nijk, primId);
2055  if(computeDistance(nijk, prim, data)) coordList.push_back(nijk);
2056  }
2057  }
2058  }
2059  }
2060 
2061  static bool computeDistance(const Coord& ijk, const Triangle& prim, VoxelizationDataType& data)
2062  {
2063  Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2064 
2065  typedef typename TreeType::ValueType ValueType;
2066 
2067  const ValueType dist = ValueType((voxelCenter -
2068  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, voxelCenter, uvw)).lengthSqr());
2069 
2070  const ValueType oldDist = data.distAcc.getValue(ijk);
2071 
2072  if (dist < oldDist) {
2073  data.distAcc.setValue(ijk, dist);
2074  data.indexAcc.setValue(ijk, prim.index);
2075  } else if (math::isExactlyEqual(dist, oldDist)) {
2076  // makes reduction deterministic when different polygons
2077  // produce the same distance value.
2078  data.indexAcc.setValueOnly(ijk, std::min(prim.index, data.indexAcc.getValue(ijk)));
2079  }
2080 
2081  return !(dist > 0.75); // true if the primitive intersects the voxel.
2082  }
2083 
2084  DataTable * const mDataTable;
2085  MeshDataAdapter const * const mMesh;
2086  Interrupter * const mInterrupter;
2087 }; // VoxelizePolygons
2088 
2089 
2091 
2092 
2093 template<typename TreeType>
2095 {
2097  typedef typename TreeType::LeafNodeType LeafNodeType;
2098 
2099  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2100  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
2101 
2102  DiffLeafNodeMask(const TreeType& rhsTree,
2103  std::vector<BoolLeafNodeType*>& lhsNodes)
2104  : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? NULL : &lhsNodes[0])
2105  {
2106  }
2107 
2108  void operator()(const tbb::blocked_range<size_t>& range) const {
2109 
2110  tree::ValueAccessor<const TreeType> acc(*mRhsTree);
2111 
2112  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2113 
2114  BoolLeafNodeType* lhsNode = mLhsNodes[n];
2115  const LeafNodeType* rhsNode = acc.probeConstLeaf(lhsNode->origin());
2116 
2117  if (rhsNode) lhsNode->topologyDifference(*rhsNode, false);
2118  }
2119  }
2120 
2121 private:
2122  TreeType const * const mRhsTree;
2123  BoolLeafNodeType ** const mLhsNodes;
2124 };
2125 
2126 
2127 template<typename LeafNodeTypeA, typename LeafNodeTypeB>
2129 {
2130  UnionValueMasks(std::vector<LeafNodeTypeA*>& nodesA, std::vector<LeafNodeTypeB*>& nodesB)
2131  : mNodesA(nodesA.empty() ? NULL : &nodesA[0])
2132  , mNodesB(nodesB.empty() ? NULL : &nodesB[0])
2133  {
2134  }
2135 
2136  void operator()(const tbb::blocked_range<size_t>& range) const {
2137  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2138  mNodesA[n]->topologyUnion(*mNodesB[n]);
2139  }
2140  }
2141 
2142 private:
2143  LeafNodeTypeA ** const mNodesA;
2144  LeafNodeTypeB ** const mNodesB;
2145 };
2146 
2147 
2148 template<typename TreeType>
2150 {
2151  typedef typename TreeType::LeafNodeType LeafNodeType;
2152 
2153  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2154  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
2155 
2156  ConstructVoxelMask(BoolTreeType& maskTree, const TreeType& tree, std::vector<LeafNodeType*>& nodes)
2157  : mTree(&tree)
2158  , mNodes(nodes.empty() ? NULL : &nodes[0])
2159  , mLocalMaskTree(false)
2160  , mMaskTree(&maskTree)
2161  {
2162  }
2163 
2165  : mTree(rhs.mTree)
2166  , mNodes(rhs.mNodes)
2167  , mLocalMaskTree(false)
2168  , mMaskTree(&mLocalMaskTree)
2169  {
2170  }
2171 
2172  void operator()(const tbb::blocked_range<size_t>& range)
2173  {
2174  typedef typename LeafNodeType::ValueOnCIter Iterator;
2175 
2177  tree::ValueAccessor<BoolTreeType> maskAcc(*mMaskTree);
2178 
2179  Coord ijk, nijk, localCorod;
2180  Index pos, npos;
2181 
2182  for (size_t n = range.begin(); n != range.end(); ++n) {
2183 
2184  LeafNodeType& node = *mNodes[n];
2185 
2186  CoordBBox bbox = node.getNodeBoundingBox();
2187  bbox.expand(-1);
2188 
2189  BoolLeafNodeType& maskNode = *maskAcc.touchLeaf(node.origin());
2190 
2191  for (Iterator it = node.cbeginValueOn(); it; ++it) {
2192  ijk = it.getCoord();
2193  pos = it.pos();
2194 
2195  localCorod = LeafNodeType::offsetToLocalCoord(pos);
2196 
2197  if (localCorod[2] < int(LeafNodeType::DIM - 1)) {
2198  npos = pos + 1;
2199  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2200  } else {
2201  nijk = ijk.offsetBy(0, 0, 1);
2202  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2203  }
2204 
2205  if (localCorod[2] > 0) {
2206  npos = pos - 1;
2207  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2208  } else {
2209  nijk = ijk.offsetBy(0, 0, -1);
2210  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2211  }
2212 
2213  if (localCorod[1] < int(LeafNodeType::DIM - 1)) {
2214  npos = pos + LeafNodeType::DIM;
2215  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2216  } else {
2217  nijk = ijk.offsetBy(0, 1, 0);
2218  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2219  }
2220 
2221  if (localCorod[1] > 0) {
2222  npos = pos - LeafNodeType::DIM;
2223  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2224  } else {
2225  nijk = ijk.offsetBy(0, -1, 0);
2226  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2227  }
2228 
2229  if (localCorod[0] < int(LeafNodeType::DIM - 1)) {
2230  npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2231  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2232  } else {
2233  nijk = ijk.offsetBy(1, 0, 0);
2234  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2235  }
2236 
2237  if (localCorod[0] > 0) {
2238  npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2239  if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2240  } else {
2241  nijk = ijk.offsetBy(-1, 0, 0);
2242  if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2243  }
2244  }
2245  }
2246  }
2247 
2248  void join(ConstructVoxelMask& rhs) { mMaskTree->merge(*rhs.mMaskTree); }
2249 
2250 private:
2251  TreeType const * const mTree;
2252  LeafNodeType ** const mNodes;
2253 
2254  BoolTreeType mLocalMaskTree;
2255  BoolTreeType * const mMaskTree;
2256 };
2257 
2258 
2260 template<typename TreeType, typename MeshDataAdapter>
2262 {
2263  typedef typename TreeType::ValueType ValueType;
2264  typedef typename TreeType::LeafNodeType LeafNodeType;
2265  typedef typename TreeType::template ValueConverter<Int32>::Type Int32TreeType;
2266  typedef typename Int32TreeType::LeafNodeType Int32LeafNodeType;
2267  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2268  typedef typename BoolTreeType::LeafNodeType BoolLeafNodeType;
2269 
2271  std::vector<BoolLeafNodeType*>& maskNodes,
2272  BoolTreeType& maskTree,
2273  TreeType& distTree,
2274  Int32TreeType& indexTree,
2275  const MeshDataAdapter& mesh,
2276  ValueType exteriorBandWidth,
2277  ValueType interiorBandWidth,
2278  ValueType voxelSize)
2279  : mMaskNodes(maskNodes.empty() ? NULL : &maskNodes[0])
2280  , mMaskTree(&maskTree)
2281  , mDistTree(&distTree)
2282  , mIndexTree(&indexTree)
2283  , mMesh(&mesh)
2284  , mNewMaskTree(false)
2285  , mDistNodes()
2286  , mUpdatedDistNodes()
2287  , mIndexNodes()
2288  , mUpdatedIndexNodes()
2289  , mExteriorBandWidth(exteriorBandWidth)
2290  , mInteriorBandWidth(interiorBandWidth)
2291  , mVoxelSize(voxelSize)
2292  {
2293  }
2294 
2295  ExpandNarrowband(const ExpandNarrowband& rhs, tbb::split)
2296  : mMaskNodes(rhs.mMaskNodes)
2297  , mMaskTree(rhs.mMaskTree)
2298  , mDistTree(rhs.mDistTree)
2299  , mIndexTree(rhs.mIndexTree)
2300  , mMesh(rhs.mMesh)
2301  , mNewMaskTree(false)
2302  , mDistNodes()
2303  , mUpdatedDistNodes()
2304  , mIndexNodes()
2305  , mUpdatedIndexNodes()
2306  , mExteriorBandWidth(rhs.mExteriorBandWidth)
2307  , mInteriorBandWidth(rhs.mInteriorBandWidth)
2308  , mVoxelSize(rhs.mVoxelSize)
2309  {
2310  }
2311 
2312  void join(ExpandNarrowband& rhs) {
2313  mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2314  mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2315 
2316  mUpdatedDistNodes.insert(mUpdatedDistNodes.end(), rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2317  mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(), rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2318 
2319  mNewMaskTree.merge(rhs.mNewMaskTree);
2320  }
2321 
2322 
2323  void operator()(const tbb::blocked_range<size_t>& range) {
2324 
2325  tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
2326  tree::ValueAccessor<TreeType> distAcc(*mDistTree);
2327  tree::ValueAccessor<Int32TreeType> indexAcc(*mIndexTree);
2328 
2329  std::vector<Int32> primitives;
2330  primitives.reserve(26);
2331 
2332  LeafNodeType * newDistNodePt = NULL;
2333  Int32LeafNodeType * newIndexNodePt = NULL;
2334 
2335  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2336 
2337  BoolLeafNodeType& maskNode = *mMaskNodes[n];
2338  if (maskNode.isEmpty()) continue;
2339 
2340  Coord ijk = maskNode.origin();
2341 
2342  bool usingNewNodes = false;
2343 
2344  LeafNodeType * distNodePt = distAcc.probeLeaf(ijk);
2345  Int32LeafNodeType * indexNodePt = indexAcc.probeLeaf(ijk);
2346 
2347  assert(!distNodePt == !indexNodePt);
2348 
2349  if (!distNodePt && !indexNodePt) {
2350 
2351  const ValueType backgroundDist = distAcc.getValue(ijk);
2352 
2353  if (!newDistNodePt && !newIndexNodePt) {
2354  newDistNodePt = new LeafNodeType(ijk, backgroundDist);
2355  newIndexNodePt = new Int32LeafNodeType(ijk, indexAcc.getValue(ijk));
2356  } else {
2357 
2358  if ((backgroundDist < ValueType(0.0)) != (newDistNodePt->getValue(0) < ValueType(0.0))) {
2359  newDistNodePt->buffer().fill(backgroundDist);
2360  }
2361 
2362  newDistNodePt->setOrigin(ijk);
2363  newIndexNodePt->setOrigin(ijk);
2364  }
2365 
2366  distNodePt = newDistNodePt;
2367  indexNodePt = newIndexNodePt;
2368 
2369  usingNewNodes = true;
2370  }
2371 
2372  bool updatedValues = false;
2373 
2374  for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2375 
2376  ijk = it.getCoord();
2377  const Index pos = it.pos();
2378 
2379  Int32 closestPrimIdx = 0;
2380  const ValueType distance =
2381  computeDistance(ijk, distAcc, indexAcc, primitives, closestPrimIdx);
2382 
2383  const bool inside = distNodePt->getValue(pos) < ValueType(0.0);
2384 
2385  if (!inside && distance < mExteriorBandWidth) {
2386  distNodePt->setValueOnly(pos, distance);
2387  indexNodePt->setValueOn(pos, closestPrimIdx);
2388  } else if (inside && distance < mInteriorBandWidth) {
2389  distNodePt->setValueOnly(pos, -distance);
2390  indexNodePt->setValueOn(pos, closestPrimIdx);
2391  } else {
2392  continue;
2393  }
2394 
2395  for (Int32 i = 0; i < 6; ++i) {
2396  newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[i]);
2397  }
2398 
2399  updatedValues = true;
2400  }
2401 
2402 
2403  if (updatedValues && usingNewNodes) {
2404 
2405  distNodePt->topologyUnion(*indexNodePt);
2406 
2407  mDistNodes.push_back(distNodePt);
2408  mIndexNodes.push_back(indexNodePt);
2409 
2410  newDistNodePt = NULL;
2411  newIndexNodePt = NULL;
2412 
2413  } else if (updatedValues) {
2414 
2415  mUpdatedDistNodes.push_back(distNodePt);
2416  mUpdatedIndexNodes.push_back(indexNodePt);
2417  }
2418  }
2419 
2420 
2421  }
2422 
2424 
2425  BoolTreeType& newMaskTree() { return mNewMaskTree; }
2426 
2427  std::vector<LeafNodeType*>& newDistNodes() { return mDistNodes; }
2428  std::vector<LeafNodeType*>& updatedDistNodes() { return mUpdatedDistNodes; }
2429 
2430  std::vector<Int32LeafNodeType*>& newIndexNodes() { return mIndexNodes; }
2431  std::vector<Int32LeafNodeType*>& updatedIndexNodes() { return mUpdatedIndexNodes; }
2432 
2433 private:
2434 
2435  ValueType
2436  computeDistance(const Coord& ijk,
2438  std::vector<Int32>& primitives, Int32& closestPrimIdx) const
2439  {
2440  ValueType minDist = std::numeric_limits<ValueType>::max();
2441  primitives.clear();
2442 
2443  const Coord ijkMin = ijk.offsetBy(-1);
2444  const Coord ijkMax = ijk.offsetBy(1);
2445  const Coord nodeMin = ijkMin & ~(LeafNodeType::DIM - 1);
2446  const Coord nodeMax = ijkMax & ~(LeafNodeType::DIM - 1);
2447 
2448  CoordBBox bbox;
2449  Coord nijk;
2450 
2451  for (nijk[0] = nodeMin[0]; nijk[0] <= nodeMax[0]; nijk[0] += LeafNodeType::DIM) {
2452  for (nijk[1] = nodeMin[1]; nijk[1] <= nodeMax[1]; nijk[1] += LeafNodeType::DIM) {
2453  for (nijk[2] = nodeMin[2]; nijk[2] <= nodeMax[2]; nijk[2] += LeafNodeType::DIM) {
2454 
2455  if (LeafNodeType* distleaf = distAcc.probeLeaf(nijk)) {
2456 
2457  bbox.min() = Coord::maxComponent(ijkMin, nijk);
2458  bbox.max() = Coord::minComponent(ijkMax, nijk.offsetBy(LeafNodeType::DIM - 1));
2459 
2460  evalLeafNode(bbox, *distleaf, *idxAcc.probeLeaf(nijk), primitives, minDist);
2461  }
2462  }
2463  }
2464  }
2465 
2466  const ValueType tmpDist = evalPrimitives(ijk, primitives, closestPrimIdx);
2467  return tmpDist > minDist ? tmpDist : minDist + mVoxelSize;
2468  }
2469 
2470  void
2471  evalLeafNode(const CoordBBox& bbox, LeafNodeType& distLeaf,
2472  Int32LeafNodeType& idxLeaf, std::vector<Int32>& primitives, ValueType& minNeighbourDist) const
2473  {
2474  ValueType tmpDist;
2475  Index xPos(0), yPos(0), pos(0);
2476 
2477  for (int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2478  xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2479  for (int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2480  yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2481  for (int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2482  pos = yPos + (z & (LeafNodeType::DIM - 1u));
2483  if (distLeaf.probeValue(pos, tmpDist)) {
2484  primitives.push_back(idxLeaf.getValue(pos));
2485  minNeighbourDist = std::min(std::abs(tmpDist), minNeighbourDist);
2486  }
2487  }
2488  }
2489  }
2490  }
2491 
2492  ValueType
2493  evalPrimitives(const Coord& ijk, std::vector<Int32>& primitives, Int32& closestPrimIdx) const
2494  {
2495  std::sort(primitives.begin(), primitives.end());
2496 
2497  Int32 lastPrim = -1;
2498  Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2499  double primDist, tmpDist, dist = std::numeric_limits<double>::max();
2500  for (size_t n = 0, N = primitives.size(); n < N; ++n) {
2501 
2502  if (primitives[n] == lastPrim) continue;
2503 
2504  lastPrim = primitives[n];
2505 
2506  const size_t polygon = size_t(lastPrim);
2507 
2508  mMesh->getIndexSpacePoint(polygon, 0, a);
2509  mMesh->getIndexSpacePoint(polygon, 1, b);
2510  mMesh->getIndexSpacePoint(polygon, 2, c);
2511 
2512  primDist = (voxelCenter -
2513  closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
2514 
2515  // Split-up quad into a second triangle
2516  if (4 == mMesh->vertexCount(polygon)) {
2517 
2518  mMesh->getIndexSpacePoint(polygon, 3, b);
2519 
2520  tmpDist = (voxelCenter - closestPointOnTriangleToPoint(
2521  a, b, c, voxelCenter, uvw)).lengthSqr();
2522 
2523  if (tmpDist < primDist) primDist = tmpDist;
2524  }
2525 
2526  if (primDist < dist) {
2527  dist = primDist;
2528  closestPrimIdx = lastPrim;
2529  }
2530  }
2531 
2532  return ValueType(std::sqrt(dist)) * mVoxelSize;
2533  }
2534 
2535 
2537 
2538 
2539  BoolLeafNodeType ** const mMaskNodes;
2540  BoolTreeType * const mMaskTree;
2541  TreeType * const mDistTree;
2542  Int32TreeType * const mIndexTree;
2543 
2544  MeshDataAdapter const * const mMesh;
2545 
2546  BoolTreeType mNewMaskTree;
2547 
2548  std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2549  std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2550 
2551  const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2552 }; // ExpandNarrowband
2553 
2554 
2555 template<typename TreeType, typename Int32TreeType, typename BoolTreeType, typename MeshDataAdapter>
2556 inline void
2558  TreeType& distTree,
2559  Int32TreeType& indexTree,
2560  BoolTreeType& maskTree,
2561  std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2562  const MeshDataAdapter& mesh,
2563  typename TreeType::ValueType exteriorBandWidth,
2564  typename TreeType::ValueType interiorBandWidth,
2565  typename TreeType::ValueType voxelSize)
2566 {
2567  typedef typename TreeType::LeafNodeType LeafNodeType;
2568  typedef typename Int32TreeType::LeafNodeType Int32LeafNodeType;
2569 
2571  expandOp(maskNodes, maskTree, distTree, indexTree,
2572  mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2573 
2574  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2575 
2576  {
2577  tree::ValueAccessor<TreeType> acc(distTree);
2578  typedef typename std::vector<LeafNodeType*> LeafNodePtVec;
2579  LeafNodePtVec& nodes = expandOp.newDistNodes();
2580  for (typename LeafNodePtVec::iterator it = nodes.begin(), end = nodes.end(); it != end; ++it) {
2581  acc.addLeaf(*it);
2582  }
2583  }
2584 
2585  {
2586  tree::ValueAccessor<Int32TreeType> acc(indexTree);
2587  typedef typename std::vector<Int32LeafNodeType*> LeafNodePtVec;
2588  LeafNodePtVec& nodes = expandOp.newIndexNodes();
2589  for (typename LeafNodePtVec::iterator it = nodes.begin(), end = nodes.end(); it != end; ++it) {
2590  acc.addLeaf(*it);
2591  }
2592  }
2593 
2594  tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.updatedIndexNodes().size()),
2596 
2597  maskTree.clear();
2598  maskTree.merge(expandOp.newMaskTree());
2599 }
2600 
2601 
2603 
2604 
2605 // Transform values (sqrt, world space scaling and sign flip if sdf)
2606 template<typename TreeType>
2608 {
2609  typedef typename TreeType::LeafNodeType LeafNodeType;
2610  typedef typename TreeType::ValueType ValueType;
2611 
2612  TransformValues(std::vector<LeafNodeType*>& nodes,
2613  ValueType voxelSize, bool unsignedDist)
2614  : mNodes(&nodes[0])
2615  , mVoxelSize(voxelSize)
2616  , mUnsigned(unsignedDist)
2617  {
2618  }
2619 
2620  void operator()(const tbb::blocked_range<size_t>& range) const {
2621 
2622  typename LeafNodeType::ValueOnIter iter;
2623 
2624  const bool udf = mUnsigned;
2625  const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2626 
2627  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2628 
2629  for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2630  ValueType& val = const_cast<ValueType&>(iter.getValue());
2631  val = w[udf || (val < ValueType(0.0))] * std::sqrt(std::abs(val));
2632  }
2633  }
2634  }
2635 
2636 private:
2637  LeafNodeType * * const mNodes;
2638  const ValueType mVoxelSize;
2639  const bool mUnsigned;
2640 };
2641 
2642 
2643 // Inactivate values outside the (exBandWidth, inBandWidth) range.
2644 template<typename TreeType>
2646 {
2647  typedef typename TreeType::LeafNodeType LeafNodeType;
2648  typedef typename TreeType::ValueType ValueType;
2649 
2650  InactivateValues(std::vector<LeafNodeType*>& nodes,
2651  ValueType exBandWidth, ValueType inBandWidth)
2652  : mNodes(nodes.empty() ? NULL : &nodes[0])
2653  , mExBandWidth(exBandWidth)
2654  , mInBandWidth(inBandWidth)
2655  {
2656  }
2657 
2658  void operator()(const tbb::blocked_range<size_t>& range) const {
2659 
2660  typename LeafNodeType::ValueOnIter iter;
2661  const ValueType exVal = mExBandWidth;
2662  const ValueType inVal = -mInBandWidth;
2663 
2664  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2665 
2666  for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2667 
2668  ValueType& val = const_cast<ValueType&>(iter.getValue());
2669 
2670  const bool inside = val < ValueType(0.0);
2671 
2672  if (inside && !(val > inVal)) {
2673  val = inVal;
2674  iter.setValueOff();
2675  } else if (!inside && !(val < exVal)) {
2676  val = exVal;
2677  iter.setValueOff();
2678  }
2679  }
2680  }
2681  }
2682 
2683 private:
2684  LeafNodeType * * const mNodes;
2685  const ValueType mExBandWidth, mInBandWidth;
2686 };
2687 
2688 
2689 template<typename TreeType>
2691 {
2692  typedef typename TreeType::LeafNodeType LeafNodeType;
2693  typedef typename TreeType::ValueType ValueType;
2694 
2695  OffsetValues(std::vector<LeafNodeType*>& nodes, ValueType offset)
2696  : mNodes(nodes.empty() ? NULL : &nodes[0]), mOffset(offset)
2697  {
2698  }
2699 
2700  void operator()(const tbb::blocked_range<size_t>& range) const {
2701 
2702  const ValueType offset = mOffset;
2703 
2704  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2705 
2706  typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2707 
2708  for (; iter; ++iter) {
2709  ValueType& val = const_cast<ValueType&>(iter.getValue());
2710  val += offset;
2711  }
2712  }
2713  }
2714 
2715 private:
2716  LeafNodeType * * const mNodes;
2717  const ValueType mOffset;
2718 };
2719 
2720 
2721 template<typename TreeType>
2723 {
2724  typedef typename TreeType::LeafNodeType LeafNodeType;
2725  typedef typename TreeType::ValueType ValueType;
2726 
2727  Renormalize(const TreeType& tree, const std::vector<LeafNodeType*>& nodes, ValueType* buffer, ValueType voxelSize)
2728  : mTree(&tree)
2729  , mNodes(nodes.empty() ? NULL : &nodes[0])
2730  , mBuffer(buffer)
2731  , mVoxelSize(voxelSize)
2732  {
2733  }
2734 
2735  void operator()(const tbb::blocked_range<size_t>& range) const
2736  {
2737  typedef math::Vec3<ValueType> Vec3Type;
2738 
2740 
2741  Coord ijk;
2742  Vec3Type up, down;
2743 
2744  const ValueType dx = mVoxelSize, invDx = ValueType(1.0) / mVoxelSize;
2745 
2746  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2747 
2748  ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2749 
2750  typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2751  for (; iter; ++iter) {
2752 
2753  const ValueType phi0 = *iter;
2754 
2755  ijk = iter.getCoord();
2756 
2757  up[0] = acc.getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2758  up[1] = acc.getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2759  up[2] = acc.getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2760 
2761  down[0] = phi0 - acc.getValue(ijk.offsetBy(-1, 0, 0));
2762  down[1] = phi0 - acc.getValue(ijk.offsetBy(0, -1, 0));
2763  down[2] = phi0 - acc.getValue(ijk.offsetBy(0, 0, -1));
2764 
2765  const ValueType normSqGradPhi = math::GudonovsNormSqrd(phi0 > 0.0, down, up);
2766 
2767  const ValueType diff = math::Sqrt(normSqGradPhi) * invDx - ValueType(1.0);
2768  const ValueType S = phi0 / (math::Sqrt(math::Pow2(phi0) + normSqGradPhi));
2769 
2770  bufferData[iter.pos()] = phi0 - dx * S * diff;
2771  }
2772  }
2773  }
2774 
2775 private:
2776  TreeType const * const mTree;
2777  LeafNodeType const * const * const mNodes;
2778  ValueType * const mBuffer;
2779 
2780  const ValueType mVoxelSize;
2781 };
2782 
2783 
2784 template<typename TreeType>
2786 {
2787  typedef typename TreeType::LeafNodeType LeafNodeType;
2788  typedef typename TreeType::ValueType ValueType;
2789 
2790  MinCombine(std::vector<LeafNodeType*>& nodes, const ValueType* buffer)
2791  : mNodes(nodes.empty() ? NULL : &nodes[0]), mBuffer(buffer)
2792  {
2793  }
2794 
2795  void operator()(const tbb::blocked_range<size_t>& range) const {
2796 
2797  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2798 
2799  const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2800 
2801  typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2802 
2803  for (; iter; ++iter) {
2804  ValueType& val = const_cast<ValueType&>(iter.getValue());
2805  val = std::min(val, bufferData[iter.pos()]);
2806  }
2807  }
2808  }
2809 
2810 private:
2811  LeafNodeType * * const mNodes;
2812  ValueType const * const mBuffer;
2813 };
2814 
2815 
2816 } // mesh_to_volume_internal namespace
2817 
2818 
2820 
2821 // Utility method implementation
2822 
2823 
2824 template <typename FloatTreeT>
2825 inline void
2826 traceExteriorBoundaries(FloatTreeT& tree)
2827 {
2829 
2830  ConnectivityTable nodeConnectivity(tree);
2831 
2832  std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
2833 
2834  for (size_t n = 0; n < nodeConnectivity.size(); ++n) {
2835  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
2836  xStartNodes.push_back(n);
2837  }
2838 
2839  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
2840  yStartNodes.push_back(n);
2841  }
2842 
2843  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
2844  zStartNodes.push_back(n);
2845  }
2846  }
2847 
2849 
2850  tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
2851  SweepingOp(SweepingOp::Z_AXIS, zStartNodes, nodeConnectivity));
2852 
2853  tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
2854  SweepingOp(SweepingOp::Y_AXIS, yStartNodes, nodeConnectivity));
2855 
2856  tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
2857  SweepingOp(SweepingOp::X_AXIS, xStartNodes, nodeConnectivity));
2858 
2859  const size_t numLeafNodes = nodeConnectivity.size();
2860  const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
2861 
2862  boost::scoped_array<bool> changedNodeMaskA(new bool[numLeafNodes]);
2863  boost::scoped_array<bool> changedNodeMaskB(new bool[numLeafNodes]);
2864  boost::scoped_array<bool> changedVoxelMask(new bool[numVoxels]);
2865 
2866  memset(changedNodeMaskA.get(), 1, sizeof(bool) * numLeafNodes);
2867  mesh_to_volume_internal::fillArray(changedVoxelMask.get(), false, numVoxels);
2868 
2869  const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
2870 
2871  bool nodesUpdated = false;
2872  do {
2873  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
2874  nodeConnectivity.nodes(), changedNodeMaskA.get()));
2875 
2876  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedPoints<FloatTreeT>(nodeConnectivity,
2877  changedNodeMaskA.get(), changedNodeMaskB.get(), changedVoxelMask.get()));
2878 
2879  changedNodeMaskA.swap(changedNodeMaskB);
2880 
2881  nodesUpdated = false;
2882  for (size_t n = 0; n < numLeafNodes; ++n) {
2883  nodesUpdated |= changedNodeMaskA[n];
2884  if (nodesUpdated) break;
2885  }
2886 
2887  if (nodesUpdated) {
2888  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
2889  nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
2890  }
2891  } while (nodesUpdated);
2892 
2893 } // void traceExteriorBoundaries()
2894 
2895 
2897 
2898 
2899 template <typename GridType, typename MeshDataAdapter, typename Interrupter>
2900 inline typename GridType::Ptr
2902  Interrupter& interrupter,
2903  const MeshDataAdapter& mesh,
2904  const math::Transform& transform,
2905  float exteriorBandWidth,
2906  float interiorBandWidth,
2907  int flags,
2908  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid)
2909 {
2910  typedef typename GridType::Ptr GridTypePtr;
2911  typedef typename GridType::TreeType TreeType;
2912  typedef typename TreeType::LeafNodeType LeafNodeType;
2913  typedef typename GridType::ValueType ValueType;
2914 
2915  typedef typename GridType::template ValueConverter<Int32>::Type Int32GridType;
2916  typedef typename Int32GridType::TreeType Int32TreeType;
2917 
2918  typedef typename TreeType::template ValueConverter<bool>::Type BoolTreeType;
2919 
2921 
2922  // Setup
2923 
2924  GridTypePtr distGrid(new GridType(std::numeric_limits<ValueType>::max()));
2925  distGrid->setTransform(transform.copy());
2926 
2927  ValueType exteriorWidth = ValueType(exteriorBandWidth);
2928  ValueType interiorWidth = ValueType(interiorBandWidth);
2929 
2930  // inf interior width is all right, this value makes the converter fill
2931  // interior regions with distance values.
2932  if (!boost::math::isfinite(exteriorWidth) || boost::math::isnan(interiorWidth)) {
2933  std::stringstream msg;
2934  msg << "Illegal narrow band width: exterior = " << exteriorWidth
2935  << ", interior = " << interiorWidth;
2936  OPENVDB_LOG_DEBUG(msg.str());
2937  return distGrid;
2938  }
2939 
2940  const ValueType voxelSize = ValueType(transform.voxelSize()[0]);
2941 
2942  if (!boost::math::isfinite(voxelSize) || math::isZero(voxelSize)) {
2943  std::stringstream msg;
2944  msg << "Illegal transform, voxel size = " << voxelSize;
2945  OPENVDB_LOG_DEBUG(msg.str());
2946  return distGrid;
2947  }
2948 
2949  // convert narrow band width from voxel units to world space units.
2950  exteriorWidth *= voxelSize;
2951  // avoid the unit conversion if the interior band width is set to
2952  // inf or std::numeric_limits<float>::max()
2953  if (interiorWidth < std::numeric_limits<ValueType>::max()) {
2954  interiorWidth *= voxelSize;
2955  }
2956 
2957  const bool computeSignedDistanceField = (flags & UNSIGNED_DISTANCE_FIELD) == 0;
2958  const bool removeIntersectingVoxels = (flags & DISABLE_INTERSECTING_VOXEL_REMOVAL) == 0;
2959  const bool renormalizeValues = (flags & DISABLE_RENORMALIZATION) == 0;
2960  const bool trimNarrowBand = (flags & DISABLE_NARROW_BAND_TRIMMING) == 0;
2961 
2962  Int32GridType* indexGrid = NULL;
2963 
2964  typename Int32GridType::Ptr temporaryIndexGrid;
2965 
2966  if (polygonIndexGrid) {
2967  indexGrid = polygonIndexGrid;
2968  } else {
2969  temporaryIndexGrid.reset(new Int32GridType(Int32(util::INVALID_IDX)));
2970  indexGrid = temporaryIndexGrid.get();
2971  }
2972 
2973  indexGrid->newTree();
2974  indexGrid->setTransform(transform.copy());
2975 
2976  if (computeSignedDistanceField) {
2977  distGrid->setGridClass(GRID_LEVEL_SET);
2978  } else {
2979  distGrid->setGridClass(GRID_UNKNOWN);
2980  interiorWidth = ValueType(0.0);
2981  }
2982 
2983  TreeType& distTree = distGrid->tree();
2984  Int32TreeType& indexTree = indexGrid->tree();
2985 
2986 
2988 
2989  // Voxelize mesh
2990 
2991  {
2992  typedef mesh_to_volume_internal::VoxelizationData<TreeType> VoxelizationDataType;
2993  typedef tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr> DataTable;
2994 
2995  DataTable data;
2997 
2998  const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
2999 
3000  tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3001 
3002  for (typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3003  VoxelizationDataType& dataItem = **i;
3004  mesh_to_volume_internal::combineData(distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3005  }
3006  }
3007 
3008  // the progress estimates are based on the observed average time for a few different
3009  // test cases and is only intended to provide some rough progression feedback to the user.
3010  if (interrupter.wasInterrupted(30)) return distGrid;
3011 
3012 
3014 
3015  // Classify interior and exterior regions
3016 
3017  if (computeSignedDistanceField) {
3018 
3019  // determines the inside/outside state for the narrow band of voxels.
3020  traceExteriorBoundaries(distTree);
3021 
3022  std::vector<LeafNodeType*> nodes;
3023  nodes.reserve(distTree.leafCount());
3024  distTree.getNodes(nodes);
3025 
3026  const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3027 
3029 
3030  tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3031 
3032  if (interrupter.wasInterrupted(45)) return distGrid;
3033 
3034  // remove voxels created by self intersecting portions of the mesh
3035  if (removeIntersectingVoxels) {
3036 
3037  tbb::parallel_for(nodeRange,
3039 
3040  tbb::parallel_for(nodeRange,
3042 
3043  tools::pruneInactive(distTree, /*threading=*/true);
3044  tools::pruneInactive(indexTree, /*threading=*/true);
3045  }
3046  }
3047 
3048  if (interrupter.wasInterrupted(50)) return distGrid;
3049 
3050  if (distTree.activeVoxelCount() == 0) {
3051  distGrid.reset((new GridType(ValueType(0.0))));
3052  return distGrid;
3053  }
3054 
3055  // transform values (world space scaling etc.)
3056  {
3057  std::vector<LeafNodeType*> nodes;
3058  nodes.reserve(distTree.leafCount());
3059  distTree.getNodes(nodes);
3060 
3061  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3062  mesh_to_volume_internal::TransformValues<TreeType>(nodes, voxelSize, !computeSignedDistanceField));
3063  }
3064 
3065  // propagate sign information into tile regions
3066  if (computeSignedDistanceField) {
3067  distTree.root().setBackground(exteriorWidth, /*updateChildNodes=*/false);
3068  tools::signedFloodFillWithValues(distTree, exteriorWidth, -interiorWidth);
3069  } else {
3070  tools::changeBackground(distTree, exteriorWidth);
3071  }
3072 
3073  if (interrupter.wasInterrupted(54)) return distGrid;
3074 
3075 
3077 
3078  // Expand the narrow band region
3079 
3080  const ValueType minBandWidth = voxelSize * ValueType(2.0);
3081 
3082  if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3083 
3084  // create the initial voxel mask.
3085  BoolTreeType maskTree(false);
3086 
3087  {
3088  std::vector<LeafNodeType*> nodes;
3089  nodes.reserve(distTree.leafCount());
3090  distTree.getNodes(nodes);
3091 
3092  mesh_to_volume_internal::ConstructVoxelMask<TreeType> op(maskTree, distTree, nodes);
3093  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3094  }
3095 
3096  // progress estimation
3097  unsigned maxIterations = std::numeric_limits<unsigned>::max();
3098 
3099  float progress = 54.0f, step = 0.0f;
3100  double estimated =
3101  2.0 * std::ceil((std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3102 
3103  if (estimated < double(maxIterations)) {
3104  maxIterations = unsigned(estimated);
3105  step = 40.0f / float(maxIterations);
3106  }
3107 
3108  std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3109 
3110  // expand
3111  unsigned count = 0;
3112  while (true) {
3113 
3114  if (interrupter.wasInterrupted(int(progress))) return distGrid;
3115 
3116  const size_t maskNodeCount = maskTree.leafCount();
3117  if (maskNodeCount == 0) break;
3118 
3119  maskNodes.clear();
3120  maskNodes.reserve(maskNodeCount);
3121  maskTree.getNodes(maskNodes);
3122 
3123  const tbb::blocked_range<size_t> range(0, maskNodes.size());
3124 
3125  tbb::parallel_for(range,
3127 
3128  mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3129  mesh, exteriorWidth, interiorWidth, voxelSize);
3130 
3131  if ((++count) >= maxIterations) break;
3132  progress += step;
3133  }
3134  }
3135 
3136  if (interrupter.wasInterrupted(94)) return distGrid;
3137 
3138  if (!polygonIndexGrid) indexGrid->clear();
3139 
3140 
3142 
3143  // Renormalize distances to smooth out bumps caused by self intersecting
3144  // and overlapping portions of the mesh and renormalize the level set.
3145 
3146  if (computeSignedDistanceField && renormalizeValues) {
3147 
3148  std::vector<LeafNodeType*> nodes;
3149  nodes.reserve(distTree.leafCount());
3150  distTree.getNodes(nodes);
3151 
3152  boost::scoped_array<ValueType> buffer(new ValueType[LeafNodeType::SIZE * nodes.size()]);
3153 
3154  const ValueType offset = ValueType(0.8 * voxelSize);
3155 
3156  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3158 
3159  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3160  mesh_to_volume_internal::Renormalize<TreeType>(distTree, nodes, buffer.get(), voxelSize));
3161 
3162  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3163  mesh_to_volume_internal::MinCombine<TreeType>(nodes, buffer.get()));
3164 
3165  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3167  }
3168 
3169  if (interrupter.wasInterrupted(99)) return distGrid;
3170 
3171 
3173 
3174  // Remove active voxels that exceed the narrow band limits
3175 
3176  if (trimNarrowBand && std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3177 
3178  std::vector<LeafNodeType*> nodes;
3179  nodes.reserve(distTree.leafCount());
3180  distTree.getNodes(nodes);
3181 
3182  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3183  mesh_to_volume_internal::InactivateValues<TreeType>(nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3184 
3185  tools::pruneLevelSet(distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3186  }
3187 
3188  return distGrid;
3189 }
3190 
3191 
3192 template <typename GridType, typename MeshDataAdapter>
3193 inline typename GridType::Ptr
3195  const MeshDataAdapter& mesh,
3196  const math::Transform& transform,
3197  float exteriorBandWidth,
3198  float interiorBandWidth,
3199  int flags,
3200  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid)
3201 {
3202  util::NullInterrupter nullInterrupter;
3203  return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3204  exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3205 }
3206 
3207 
3209 
3210 
3212 template<typename GridType>
3213 inline typename boost::enable_if<boost::is_floating_point<typename GridType::ValueType>,
3214 typename GridType::Ptr>::type
3216  const openvdb::math::Transform& xform,
3217  const std::vector<Vec3s>& points,
3218  const std::vector<Vec3I>& triangles,
3219  const std::vector<Vec4I>& quads,
3220  float exBandWidth,
3221  float inBandWidth,
3222  bool unsignedDistanceField = false)
3223 {
3224  if (points.empty()) {
3225  return typename GridType::Ptr(new GridType(typename GridType::ValueType(exBandWidth)));
3226  }
3227 
3228  const size_t numPoints = points.size();
3229  boost::scoped_array<Vec3s> indexSpacePoints(new Vec3s[numPoints]);
3230 
3231  // transform points to local grid index space
3232  tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3234  &points[0], indexSpacePoints.get(), xform));
3235 
3236  const int conversionFlags = unsignedDistanceField ? UNSIGNED_DISTANCE_FIELD : 0;
3237 
3238  if (quads.empty()) {
3239 
3241  mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3242 
3243  return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3244 
3245  } else if (triangles.empty()) {
3246 
3248  mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3249 
3250  return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3251  }
3252 
3253  // pack primitives
3254 
3255  const size_t numPrimitives = triangles.size() + quads.size();
3256  boost::scoped_array<Vec4I> prims(new Vec4I[numPrimitives]);
3257 
3258  for (size_t n = 0, N = triangles.size(); n < N; ++n) {
3259  const Vec3I& triangle = triangles[n];
3260  Vec4I& prim = prims[n];
3261  prim[0] = triangle[0];
3262  prim[1] = triangle[1];
3263  prim[2] = triangle[2];
3264  prim[3] = util::INVALID_IDX;
3265  }
3266 
3267  const size_t offset = triangles.size();
3268  for (size_t n = 0, N = quads.size(); n < N; ++n) {
3269  prims[offset + n] = quads[n];
3270  }
3271 
3273  mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3274 
3275  return meshToVolume<GridType>(mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3276 }
3277 
3278 
3281 template<typename GridType>
3282 inline typename boost::disable_if<boost::is_floating_point<typename GridType::ValueType>,
3283 typename GridType::Ptr>::type
3285  const math::Transform& /*xform*/,
3286  const std::vector<Vec3s>& /*points*/,
3287  const std::vector<Vec3I>& /*triangles*/,
3288  const std::vector<Vec4I>& /*quads*/,
3289  float /*exBandWidth*/,
3290  float /*inBandWidth*/,
3291  bool /*unsignedDistanceField*/ = false)
3292 {
3294  "mesh to volume conversion is supported only for scalar floating-point grids");
3295 }
3296 
3297 
3299 
3300 
3301 template<typename GridType>
3302 inline typename GridType::Ptr
3304  const openvdb::math::Transform& xform,
3305  const std::vector<Vec3s>& points,
3306  const std::vector<Vec3I>& triangles,
3307  float halfWidth)
3308 {
3309  std::vector<Vec4I> quads(0);
3310  return doMeshConversion<GridType>(xform, points, triangles, quads,
3311  halfWidth, halfWidth);
3312 }
3313 
3314 
3315 template<typename GridType>
3316 inline typename GridType::Ptr
3318  const openvdb::math::Transform& xform,
3319  const std::vector<Vec3s>& points,
3320  const std::vector<Vec4I>& quads,
3321  float halfWidth)
3322 {
3323  std::vector<Vec3I> triangles(0);
3324  return doMeshConversion<GridType>(xform, points, triangles, quads,
3325  halfWidth, halfWidth);
3326 }
3327 
3328 
3329 template<typename GridType>
3330 inline typename GridType::Ptr
3332  const openvdb::math::Transform& xform,
3333  const std::vector<Vec3s>& points,
3334  const std::vector<Vec3I>& triangles,
3335  const std::vector<Vec4I>& quads,
3336  float halfWidth)
3337 {
3338  return doMeshConversion<GridType>(xform, points, triangles, quads,
3339  halfWidth, halfWidth);
3340 }
3341 
3342 
3343 template<typename GridType>
3344 inline typename GridType::Ptr
3346  const openvdb::math::Transform& xform,
3347  const std::vector<Vec3s>& points,
3348  const std::vector<Vec3I>& triangles,
3349  const std::vector<Vec4I>& quads,
3350  float exBandWidth,
3351  float inBandWidth)
3352 {
3353  return doMeshConversion<GridType>(xform, points, triangles,
3354  quads, exBandWidth, inBandWidth);
3355 }
3356 
3357 
3358 template<typename GridType>
3359 inline typename GridType::Ptr
3361  const openvdb::math::Transform& xform,
3362  const std::vector<Vec3s>& points,
3363  const std::vector<Vec3I>& triangles,
3364  const std::vector<Vec4I>& quads,
3365  float bandWidth)
3366 {
3367  return doMeshConversion<GridType>(xform, points, triangles, quads,
3368  bandWidth, bandWidth, true);
3369 }
3370 
3371 
3373 
3374 
3375 // Required by several of the tree nodes
3376 inline std::ostream&
3377 operator<<(std::ostream& ostr, const MeshToVoxelEdgeData::EdgeData& rhs)
3378 {
3379  ostr << "{[ " << rhs.mXPrim << ", " << rhs.mXDist << "]";
3380  ostr << " [ " << rhs.mYPrim << ", " << rhs.mYDist << "]";
3381  ostr << " [ " << rhs.mZPrim << ", " << rhs.mZDist << "]}";
3382  return ostr;
3383 }
3384 
3385 // Required by math::Abs
3386 inline MeshToVoxelEdgeData::EdgeData
3388 {
3389  return x;
3390 }
3391 
3392 
3394 
3395 
3397 {
3398 public:
3399 
3400  GenEdgeData(
3401  const std::vector<Vec3s>& pointList,
3402  const std::vector<Vec4I>& polygonList);
3403 
3404  void run(bool threaded = true);
3405 
3406  GenEdgeData(GenEdgeData& rhs, tbb::split);
3407  inline void operator() (const tbb::blocked_range<size_t> &range);
3408  inline void join(GenEdgeData& rhs);
3409 
3410  inline TreeType& tree() { return mTree; }
3411 
3412 private:
3413  void operator=(const GenEdgeData&) {}
3414 
3415  struct Primitive { Vec3d a, b, c, d; Int32 index; };
3416 
3417  template<bool IsQuad>
3418  inline void voxelize(const Primitive&);
3419 
3420  template<bool IsQuad>
3421  inline bool evalPrimitive(const Coord&, const Primitive&);
3422 
3423  inline bool rayTriangleIntersection( const Vec3d& origin, const Vec3d& dir,
3424  const Vec3d& a, const Vec3d& b, const Vec3d& c, double& t);
3425 
3426 
3427  TreeType mTree;
3428  Accessor mAccessor;
3429 
3430  const std::vector<Vec3s>& mPointList;
3431  const std::vector<Vec4I>& mPolygonList;
3432 
3433  // Used internally for acceleration
3434  typedef TreeType::ValueConverter<Int32>::Type IntTreeT;
3435  IntTreeT mLastPrimTree;
3436  tree::ValueAccessor<IntTreeT> mLastPrimAccessor;
3437 }; // class MeshToVoxelEdgeData::GenEdgeData
3438 
3439 
3440 inline
3441 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3442  const std::vector<Vec3s>& pointList,
3443  const std::vector<Vec4I>& polygonList)
3444  : mTree(EdgeData())
3445  , mAccessor(mTree)
3446  , mPointList(pointList)
3447  , mPolygonList(polygonList)
3448  , mLastPrimTree(Int32(util::INVALID_IDX))
3449  , mLastPrimAccessor(mLastPrimTree)
3450 {
3451 }
3452 
3453 
3454 inline
3456  : mTree(EdgeData())
3457  , mAccessor(mTree)
3458  , mPointList(rhs.mPointList)
3459  , mPolygonList(rhs.mPolygonList)
3460  , mLastPrimTree(Int32(util::INVALID_IDX))
3461  , mLastPrimAccessor(mLastPrimTree)
3462 {
3463 }
3464 
3465 
3466 inline void
3468 {
3469  if (threaded) {
3470  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
3471  } else {
3472  (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3473  }
3474 }
3475 
3476 
3477 inline void
3479 {
3480  typedef TreeType::RootNodeType RootNodeType;
3481  typedef RootNodeType::NodeChainType NodeChainType;
3482  BOOST_STATIC_ASSERT(boost::mpl::size<NodeChainType>::value > 1);
3483  typedef boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type InternalNodeType;
3484 
3485  Coord ijk;
3486  Index offset;
3487 
3488  rhs.mTree.clearAllAccessors();
3489 
3490  TreeType::LeafIter leafIt = rhs.mTree.beginLeaf();
3491  for ( ; leafIt; ++leafIt) {
3492  ijk = leafIt->origin();
3493 
3494  TreeType::LeafNodeType* lhsLeafPt = mTree.probeLeaf(ijk);
3495 
3496  if (!lhsLeafPt) {
3497 
3498  mAccessor.addLeaf(rhs.mAccessor.probeLeaf(ijk));
3499  InternalNodeType* node = rhs.mAccessor.getNode<InternalNodeType>();
3500  node->stealNode<TreeType::LeafNodeType>(ijk, EdgeData(), false);
3501  rhs.mAccessor.clear();
3502 
3503  } else {
3504 
3505  TreeType::LeafNodeType::ValueOnCIter it = leafIt->cbeginValueOn();
3506  for ( ; it; ++it) {
3507 
3508  offset = it.pos();
3509  const EdgeData& rhsValue = it.getValue();
3510 
3511  if (!lhsLeafPt->isValueOn(offset)) {
3512  lhsLeafPt->setValueOn(offset, rhsValue);
3513  } else {
3514 
3515  EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
3516 
3517  if (rhsValue.mXDist < lhsValue.mXDist) {
3518  lhsValue.mXDist = rhsValue.mXDist;
3519  lhsValue.mXPrim = rhsValue.mXPrim;
3520  }
3521 
3522  if (rhsValue.mYDist < lhsValue.mYDist) {
3523  lhsValue.mYDist = rhsValue.mYDist;
3524  lhsValue.mYPrim = rhsValue.mYPrim;
3525  }
3526 
3527  if (rhsValue.mZDist < lhsValue.mZDist) {
3528  lhsValue.mZDist = rhsValue.mZDist;
3529  lhsValue.mZPrim = rhsValue.mZPrim;
3530  }
3531 
3532  }
3533  } // end value iteration
3534  }
3535  } // end leaf iteration
3536 }
3537 
3538 
3539 inline void
3540 MeshToVoxelEdgeData::GenEdgeData::operator()(const tbb::blocked_range<size_t> &range)
3541 {
3542  Primitive prim;
3543 
3544  for (size_t n = range.begin(); n < range.end(); ++n) {
3545 
3546  const Vec4I& verts = mPolygonList[n];
3547 
3548  prim.index = Int32(n);
3549  prim.a = Vec3d(mPointList[verts[0]]);
3550  prim.b = Vec3d(mPointList[verts[1]]);
3551  prim.c = Vec3d(mPointList[verts[2]]);
3552 
3553  if (util::INVALID_IDX != verts[3]) {
3554  prim.d = Vec3d(mPointList[verts[3]]);
3555  voxelize<true>(prim);
3556  } else {
3557  voxelize<false>(prim);
3558  }
3559  }
3560 }
3561 
3562 
3563 template<bool IsQuad>
3564 inline void
3565 MeshToVoxelEdgeData::GenEdgeData::voxelize(const Primitive& prim)
3566 {
3567  std::deque<Coord> coordList;
3568  Coord ijk, nijk;
3569 
3570  ijk = Coord::floor(prim.a);
3571  coordList.push_back(ijk);
3572 
3573  evalPrimitive<IsQuad>(ijk, prim);
3574 
3575  while (!coordList.empty()) {
3576 
3577  ijk = coordList.back();
3578  coordList.pop_back();
3579 
3580  for (Int32 i = 0; i < 26; ++i) {
3581  nijk = ijk + util::COORD_OFFSETS[i];
3582 
3583  if (prim.index != mLastPrimAccessor.getValue(nijk)) {
3584  mLastPrimAccessor.setValue(nijk, prim.index);
3585  if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3586  }
3587  }
3588  }
3589 }
3590 
3591 
3592 template<bool IsQuad>
3593 inline bool
3594 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(const Coord& ijk, const Primitive& prim)
3595 {
3596  Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3597  bool intersecting = false;
3598  double t;
3599 
3600  EdgeData edgeData;
3601  mAccessor.probeValue(ijk, edgeData);
3602 
3603  // Evaluate first triangle
3604  double dist = (org -
3605  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, org, uvw)).lengthSqr();
3606 
3607  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3608  if (t < edgeData.mXDist) {
3609  edgeData.mXDist = float(t);
3610  edgeData.mXPrim = prim.index;
3611  intersecting = true;
3612  }
3613  }
3614 
3615  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3616  if (t < edgeData.mYDist) {
3617  edgeData.mYDist = float(t);
3618  edgeData.mYPrim = prim.index;
3619  intersecting = true;
3620  }
3621  }
3622 
3623  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3624  if (t < edgeData.mZDist) {
3625  edgeData.mZDist = float(t);
3626  edgeData.mZPrim = prim.index;
3627  intersecting = true;
3628  }
3629  }
3630 
3631  if (IsQuad) {
3632  // Split quad into a second triangle and calculate distance.
3633  double secondDist = (org -
3634  closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, org, uvw)).lengthSqr();
3635 
3636  if (secondDist < dist) dist = secondDist;
3637 
3638  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3639  if (t < edgeData.mXDist) {
3640  edgeData.mXDist = float(t);
3641  edgeData.mXPrim = prim.index;
3642  intersecting = true;
3643  }
3644  }
3645 
3646  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3647  if (t < edgeData.mYDist) {
3648  edgeData.mYDist = float(t);
3649  edgeData.mYPrim = prim.index;
3650  intersecting = true;
3651  }
3652  }
3653 
3654  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3655  if (t < edgeData.mZDist) {
3656  edgeData.mZDist = float(t);
3657  edgeData.mZPrim = prim.index;
3658  intersecting = true;
3659  }
3660  }
3661  }
3662 
3663  if (intersecting) mAccessor.setValue(ijk, edgeData);
3664 
3665  return (dist < 0.86602540378443861);
3666 }
3667 
3668 
3669 inline bool
3670 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3671  const Vec3d& origin, const Vec3d& dir,
3672  const Vec3d& a, const Vec3d& b, const Vec3d& c,
3673  double& t)
3674 {
3675  // Check if ray is parallel with triangle
3676 
3677  Vec3d e1 = b - a;
3678  Vec3d e2 = c - a;
3679  Vec3d s1 = dir.cross(e2);
3680 
3681  double divisor = s1.dot(e1);
3682  if (!(std::abs(divisor) > 0.0)) return false;
3683 
3684  // Compute barycentric coordinates
3685 
3686  double inv_divisor = 1.0 / divisor;
3687  Vec3d d = origin - a;
3688  double b1 = d.dot(s1) * inv_divisor;
3689 
3690  if (b1 < 0.0 || b1 > 1.0) return false;
3691 
3692  Vec3d s2 = d.cross(e1);
3693  double b2 = dir.dot(s2) * inv_divisor;
3694 
3695  if (b2 < 0.0 || (b1 + b2) > 1.0) return false;
3696 
3697  // Compute distance to intersection point
3698 
3699  t = e2.dot(s2) * inv_divisor;
3700  return (t < 0.0) ? false : true;
3701 }
3702 
3703 
3705 
3706 
3707 inline
3709  : mTree(EdgeData())
3710 {
3711 }
3712 
3713 
3714 inline void
3716  const std::vector<Vec3s>& pointList,
3717  const std::vector<Vec4I>& polygonList)
3718 {
3719  GenEdgeData converter(pointList, polygonList);
3720  converter.run();
3721 
3722  mTree.clear();
3723  mTree.merge(converter.tree());
3724 }
3725 
3726 
3727 inline void
3729  Accessor& acc,
3730  const Coord& ijk,
3731  std::vector<Vec3d>& points,
3732  std::vector<Index32>& primitives)
3733 {
3734  EdgeData data;
3735  Vec3d point;
3736 
3737  Coord coord = ijk;
3738 
3739  if (acc.probeValue(coord, data)) {
3740 
3741  if (data.mXPrim != util::INVALID_IDX) {
3742  point[0] = double(coord[0]) + data.mXDist;
3743  point[1] = double(coord[1]);
3744  point[2] = double(coord[2]);
3745 
3746  points.push_back(point);
3747  primitives.push_back(data.mXPrim);
3748  }
3749 
3750  if (data.mYPrim != util::INVALID_IDX) {
3751  point[0] = double(coord[0]);
3752  point[1] = double(coord[1]) + data.mYDist;
3753  point[2] = double(coord[2]);
3754 
3755  points.push_back(point);
3756  primitives.push_back(data.mYPrim);
3757  }
3758 
3759  if (data.mZPrim != util::INVALID_IDX) {
3760  point[0] = double(coord[0]);
3761  point[1] = double(coord[1]);
3762  point[2] = double(coord[2]) + data.mZDist;
3763 
3764  points.push_back(point);
3765  primitives.push_back(data.mZPrim);
3766  }
3767 
3768  }
3769 
3770  coord[0] += 1;
3771 
3772  if (acc.probeValue(coord, data)) {
3773 
3774  if (data.mYPrim != util::INVALID_IDX) {
3775  point[0] = double(coord[0]);
3776  point[1] = double(coord[1]) + data.mYDist;
3777  point[2] = double(coord[2]);
3778 
3779  points.push_back(point);
3780  primitives.push_back(data.mYPrim);
3781  }
3782 
3783  if (data.mZPrim != util::INVALID_IDX) {
3784  point[0] = double(coord[0]);
3785  point[1] = double(coord[1]);
3786  point[2] = double(coord[2]) + data.mZDist;
3787 
3788  points.push_back(point);
3789  primitives.push_back(data.mZPrim);
3790  }
3791  }
3792 
3793  coord[2] += 1;
3794 
3795  if (acc.probeValue(coord, data)) {
3796  if (data.mYPrim != util::INVALID_IDX) {
3797  point[0] = double(coord[0]);
3798  point[1] = double(coord[1]) + data.mYDist;
3799  point[2] = double(coord[2]);
3800 
3801  points.push_back(point);
3802  primitives.push_back(data.mYPrim);
3803  }
3804  }
3805 
3806  coord[0] -= 1;
3807 
3808  if (acc.probeValue(coord, data)) {
3809 
3810  if (data.mXPrim != util::INVALID_IDX) {
3811  point[0] = double(coord[0]) + data.mXDist;
3812  point[1] = double(coord[1]);
3813  point[2] = double(coord[2]);
3814 
3815  points.push_back(point);
3816  primitives.push_back(data.mXPrim);
3817  }
3818 
3819  if (data.mYPrim != util::INVALID_IDX) {
3820  point[0] = double(coord[0]);
3821  point[1] = double(coord[1]) + data.mYDist;
3822  point[2] = double(coord[2]);
3823 
3824  points.push_back(point);
3825  primitives.push_back(data.mYPrim);
3826  }
3827  }
3828 
3829 
3830  coord[1] += 1;
3831 
3832  if (acc.probeValue(coord, data)) {
3833 
3834  if (data.mXPrim != util::INVALID_IDX) {
3835  point[0] = double(coord[0]) + data.mXDist;
3836  point[1] = double(coord[1]);
3837  point[2] = double(coord[2]);
3838 
3839  points.push_back(point);
3840  primitives.push_back(data.mXPrim);
3841  }
3842  }
3843 
3844  coord[2] -= 1;
3845 
3846  if (acc.probeValue(coord, data)) {
3847 
3848  if (data.mXPrim != util::INVALID_IDX) {
3849  point[0] = double(coord[0]) + data.mXDist;
3850  point[1] = double(coord[1]);
3851  point[2] = double(coord[2]);
3852 
3853  points.push_back(point);
3854  primitives.push_back(data.mXPrim);
3855  }
3856 
3857  if (data.mZPrim != util::INVALID_IDX) {
3858  point[0] = double(coord[0]);
3859  point[1] = double(coord[1]);
3860  point[2] = double(coord[2]) + data.mZDist;
3861 
3862  points.push_back(point);
3863  primitives.push_back(data.mZPrim);
3864  }
3865  }
3866 
3867  coord[0] += 1;
3868 
3869  if (acc.probeValue(coord, data)) {
3870 
3871  if (data.mZPrim != util::INVALID_IDX) {
3872  point[0] = double(coord[0]);
3873  point[1] = double(coord[1]);
3874  point[2] = double(coord[2]) + data.mZDist;
3875 
3876  points.push_back(point);
3877  primitives.push_back(data.mZPrim);
3878  }
3879  }
3880 }
3881 
3882 
3883 template<typename GridType, typename VecType>
3884 inline typename GridType::Ptr
3886  const openvdb::math::Transform& xform,
3887  typename VecType::ValueType halfWidth)
3888 {
3889  const Vec3s pmin = Vec3s(xform.worldToIndex(bbox.min()));
3890  const Vec3s pmax = Vec3s(xform.worldToIndex(bbox.max()));
3891 
3892  Vec3s points[8];
3893  points[0] = Vec3s(pmin[0], pmin[1], pmin[2]);
3894  points[1] = Vec3s(pmin[0], pmin[1], pmax[2]);
3895  points[2] = Vec3s(pmax[0], pmin[1], pmax[2]);
3896  points[3] = Vec3s(pmax[0], pmin[1], pmin[2]);
3897  points[4] = Vec3s(pmin[0], pmax[1], pmin[2]);
3898  points[5] = Vec3s(pmin[0], pmax[1], pmax[2]);
3899  points[6] = Vec3s(pmax[0], pmax[1], pmax[2]);
3900  points[7] = Vec3s(pmax[0], pmax[1], pmin[2]);
3901 
3902  Vec4I faces[6];
3903  faces[0] = Vec4I(0, 1, 2, 3); // bottom
3904  faces[1] = Vec4I(7, 6, 5, 4); // top
3905  faces[2] = Vec4I(4, 5, 1, 0); // front
3906  faces[3] = Vec4I(6, 7, 3, 2); // back
3907  faces[4] = Vec4I(0, 3, 7, 4); // left
3908  faces[5] = Vec4I(1, 5, 6, 2); // right
3909 
3910  QuadAndTriangleDataAdapter<Vec3s, Vec4I> mesh(points, 8, faces, 6);
3911 
3912  return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
3913 }
3914 
3915 
3916 } // namespace tools
3917 } // namespace OPENVDB_VERSION_NAME
3918 } // namespace openvdb
3919 
3920 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
3921 
3922 // Copyright (c) 2012-2015 DreamWorks Animation LLC
3923 // All rights reserved. This software is distributed under the
3924 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void releaseLeafNodes(TreeType &tree)
Definition: MeshToVolume.h:1758
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:374
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2647
std::vector< Int32LeafNodeType * > & updatedIndexNodes()
Definition: MeshToVolume.h:2431
Index32 Index
Definition: Types.h:58
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2795
EdgeData(float dist=1.0)
Definition: MeshToVolume.h:409
GridType::Ptr meshToVolume(Interrupter &interrupter, const MeshDataAdapter &mesh, const math::Transform &transform, float exteriorBandWidth=3.0f, float interiorBandWidth=3.0f, int flags=0, typename GridType::template ValueConverter< Int32 >::Type *polygonIndexGrid=NULL)
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
Definition: MeshToVolume.h:2901
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1137
OPENVDB_API const Index32 INVALID_IDX
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:545
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:424
const size_t * offsetsPrevX() const
Definition: MeshToVolume.h:758
static ValueType epsilon()
Definition: MeshToVolume.h:518
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1627
UCharTreeType primIdTree
Definition: MeshToVolume.h:1874
tree::ValueAccessor< Int32TreeType > Int32TreeAcc
Definition: MeshToVolume.h:1853
float mYDist
Definition: MeshToVolume.h:432
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2700
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1679
const ValueType mValue
Definition: MeshToVolume.h:1080
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:349
GridType::Ptr createLevelSetBox(const math::BBox< VecType > &bbox, const openvdb::math::Transform &xform, typename VecType::ValueType halfWidth=LEVEL_SET_HALF_WIDTH)
Return a grid of type GridType containing a narrow-band level set representation of a box...
Definition: MeshToVolume.h:3885
void fillArray(ValueType *array, const ValueType val, const size_t length)
Definition: MeshToVolume.h:1086
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:658
DiffLeafNodeMask(const TreeType &rhsTree, std::vector< BoolLeafNodeType * > &lhsNodes)
Definition: MeshToVolume.h:2102
Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:534
const size_t * offsetsPrevY() const
Definition: MeshToVolume.h:761
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1126
tree::ValueAccessor< TreeType > Accessor
Definition: MeshToVolume.h:437
math::Transform const *const mXform
Definition: MeshToVolume.h:511
boost::disable_if< boost::is_floating_point< typename GridType::ValueType >, typename GridType::Ptr >::type doMeshConversion(const math::Transform &, const std::vector< Vec3s > &, const std::vector< Vec3I > &, const std::vector< Vec4I > &, float, float, bool=false)
Definition: MeshToVolume.h:3284
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:348
void join(ExpandNarrowband &rhs)
Definition: MeshToVolume.h:2312
Int32TreeType *const mIndexTree
Definition: MeshToVolume.h:1731
size_t polygonCount() const
Definition: MeshToVolume.h:209
StealUniqueLeafNodes(TreeType &lhsTree, TreeType &rhsTree, std::vector< LeafNodeType * > &overlappingNodes)
Definition: MeshToVolume.h:1777
bool *const mChangedNodeMask
Definition: MeshToVolume.h:1063
Index32 mXPrim
Definition: MeshToVolume.h:433
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:610
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << ...'.
Definition: logging.h:45
SeedPoints(ConnectivityTable &connectivity, bool *changedNodeMask, bool *nodeMask, bool *changedVoxelMask)
Definition: MeshToVolume.h:1140
size_t pointCount() const
Definition: MeshToVolume.h:210
FillArray(ValueType *array, const ValueType v)
Definition: MeshToVolume.h:1070
QuadAndTriangleDataAdapter(const std::vector< PointType > &points, const std::vector< PolygonType > &polygons)
Definition: MeshToVolume.h:191
LeafNodeType **const mNodes
Definition: MeshToVolume.h:642
VoxelizationData< TreeType > VoxelizationDataType
Definition: MeshToVolume.h:1900
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2263
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1143
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
Renormalize(const TreeType &tree, const std::vector< LeafNodeType * > &nodes, ValueType *buffer, ValueType voxelSize)
Definition: MeshToVolume.h:2727
Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:2266
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2787
bool *const mChangedVoxelMask
Definition: MeshToVolume.h:1128
ReleaseChildNodes(NodeType **nodes)
Definition: MeshToVolume.h:1741
Definition: Math.h:841
static ValueType minNarrowBandWidth()
Definition: MeshToVolume.h:519
void join(GenEdgeData &rhs)
Definition: MeshToVolume.h:3478
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2154
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:651
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1477
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
NodeType **const mNodes
Definition: MeshToVolume.h:1752
Int32TreeAcc indexAcc
Definition: MeshToVolume.h:1872
Index32 mYPrim
Definition: MeshToVolume.h:433
Efficient multi-threaded replacement of the background values in tree.
std::vector< Int32LeafNodeType * > & newIndexNodes()
Definition: MeshToVolume.h:2430
void expandNarrowband(TreeType &distTree, Int32TreeType &indexTree, BoolTreeType &maskTree, std::vector< typename BoolTreeType::LeafNodeType * > &maskNodes, const MeshDataAdapter &mesh, typename TreeType::ValueType exteriorBandWidth, typename TreeType::ValueType interiorBandWidth, typename TreeType::ValueType voxelSize)
Definition: MeshToVolume.h:2557
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:520
Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:1294
uint32_t Index32
Definition: Types.h:56
StashOriginAndStoreOffset(std::vector< LeafNodeType * > &nodes, Coord *coordinates)
Definition: MeshToVolume.h:607
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:3345
bool operator==(const EdgeData &rhs) const
Definition: MeshToVolume.h:427
tree::Tree4< EdgeData, 5, 4, 3 >::Type TreeType
Definition: MeshToVolume.h:436
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1108
ValueType *const mArray
Definition: MeshToVolume.h:1079
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:3728
Definition: Types.h:205
ConstructVoxelMask(BoolTreeType &maskTree, const TreeType &tree, std::vector< LeafNodeType * > &nodes)
Definition: MeshToVolume.h:2156
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1313
size_t findNeighbourNode(tree::ValueAccessor< const TreeType > &acc, const Coord &start, const Coord &step) const
Definition: MeshToVolume.h:686
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...
std::vector< LeafNodeType * > & nodes()
Definition: MeshToVolume.h:753
TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1681
const size_t * offsetsNextZ() const
Definition: MeshToVolume.h:763
void changeBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &background, bool threaded=true, size_t grainSize=32)
Replace the background value in all the nodes of a tree.
Definition: ChangeBackground.h:230
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:158
void operator()() const
Definition: MeshToVolume.h:1785
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1847
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:267
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:3715
Definition: Mat4.h:51
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2610
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:105
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:605
PointType const *const mPointsIn
Definition: MeshToVolume.h:509
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1148
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:170
Definition: Tree.h:203
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:292
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:232
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:780
Defined various multi-threaded utility functions for trees.
MeshToVolumeFlags
Mesh to volume conversion flags.
Definition: MeshToVolume.h:85
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1775
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1729
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:3331
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2172
const Vec3T & min() const
Return a const reference to the minimum point of the BBox.
Definition: BBox.h:81
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:1664
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:511
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:256
TreeType::template ValueConverter< unsigned char >::Type UCharTreeType
Definition: MeshToVolume.h:1850
MeshToVoxelEdgeData()
Definition: MeshToVolume.h:3708
void seedFill(LeafNodeType &node)
Definition: MeshToVolume.h:902
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:223
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1099
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2097
const Vec3T & max() const
Return a const reference to the maximum point of the BBox.
Definition: BBox.h:84
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:612
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2099
bool scanFill(LeafNodeType &node)
Definition: MeshToVolume.h:978
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:533
EdgeData operator-() const
Definition: MeshToVolume.h:424
Real GudonovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:353
TreeType & tree()
Definition: MeshToVolume.h:3410
#define OPENVDB_VERSION_NAME
Definition: version.h:43
bool processZ(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1171
Int32TreeType indexTree
Definition: MeshToVolume.h:1871
void run(bool threaded=true)
Definition: MeshToVolume.h:3467
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2658
size_t vertexCount(size_t n) const
Vertex count for polygon n.
Definition: MeshToVolume.h:213
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:263
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:120
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1774
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2620
void combineData(DistTreeType &lhsDist, IndexTreeType &lhsIdx, DistTreeType &rhsDist, IndexTreeType &rhsIdx)
Definition: MeshToVolume.h:1814
LeafNodeConnectivityTable(TreeType &tree)
Definition: MeshToVolume.h:721
MinCombine(std::vector< LeafNodeType * > &nodes, const ValueType *buffer)
Definition: MeshToVolume.h:2790
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1062
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2153
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:3540
Propagates the sign of distance values from the active voxels in the narrow band to the inactive valu...
static bool check(const ValueType v)
Definition: MeshToVolume.h:1683
OffsetValues(std::vector< LeafNodeType * > &nodes, ValueType offset)
Definition: MeshToVolume.h:2695
FloatTreeAcc distAcc
Definition: MeshToVolume.h:1869
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
ConstructVoxelMask(ConstructVoxelMask &rhs, tbb::split)
Definition: MeshToVolume.h:2164
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1693
EdgeData operator-(const T &) const
Definition: MeshToVolume.h:423
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:628
Definition: Mat.h:146
BoolTreeType & newMaskTree()
Definition: MeshToVolume.h:2425
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:429
float mXDist
Definition: MeshToVolume.h:432
void pruneLevelSet(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing nodes whose values are all inactive with inactive ...
Definition: Prune.h:402
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2264
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2268
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1628
Definition: Exceptions.h:39
ValidateIntersectingVoxels(TreeType &tree, std::vector< LeafNodeType * > &nodes)
Definition: MeshToVolume.h:1632
Ptr copy() const
Definition: Transform.h:77
SyncVoxelMask(std::vector< LeafNodeType * > &nodes, const bool *changedNodeMask, bool *changedVoxelMask)
Definition: MeshToVolume.h:1101
BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2100
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1045
tree::ValueAccessor< TreeType > AccessorType
Definition: MeshToVolume.h:2096
Coord const *const mCoordinates
Definition: MeshToVolume.h:643
tree::ValueAccessor< UCharTreeType > UCharTreeAcc
Definition: MeshToVolume.h:1854
TreeType::ValueType ValueType
Definition: MeshToVolume.h:779
Definition: Math.h:840
SweepExteriorSign(Axis axis, const std::vector< size_t > &startNodeIndices, ConnectivityTable &connectivity)
Definition: MeshToVolume.h:783
const std::vector< LeafNodeType * > & nodes() const
Definition: MeshToVolume.h:754
TBB body object to voxelize a mesh of triangles and/or quads into a collection of VDB grids...
Definition: MeshToVolume.h:1844
RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:211
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1098
const size_t * offsetsNextY() const
Definition: MeshToVolume.h:760
Vec3< double > Vec3d
Definition: Vec3.h:643
Definition: Types.h:206
Axis
Definition: Math.h:838
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2725
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2693
virtual void clear()
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:438
Type Pow2(Type x)
Return .
Definition: Math.h:514
ConnectivityTable *const mConnectivity
Definition: MeshToVolume.h:1279
void join(ConstructVoxelMask &rhs)
Definition: MeshToVolume.h:2248
LeafNodeConnectivityTable< TreeType > ConnectivityTable
Definition: MeshToVolume.h:1138
Axis-aligned bounding box.
Definition: BBox.h:47
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:393
SeedFillExteriorSign(std::vector< LeafNodeType * > &nodes, bool *changedNodeMask)
Definition: MeshToVolume.h:1047
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1680
bool traceVoxelLine(LeafNodeType &node, Index pos, Index step) const
Definition: MeshToVolume.h:866
TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2267
ExpandNarrowband(const ExpandNarrowband &rhs, tbb::split)
Definition: MeshToVolume.h:2295
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:790
bool checkNeighbours(const Coord &ijk, AccessorType &acc, bool(&mask)[26])
Definition: MeshToVolume.h:1612
bool processX(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1243
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition: MeshToVolume.h:2826
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1053
const size_t * offsetsPrevZ() const
Definition: MeshToVolume.h:764
TransformValues(std::vector< LeafNodeType * > &nodes, ValueType voxelSize, bool unsignedDist)
Definition: MeshToVolume.h:2612
tbb::enumerable_thread_specific< LocalData > LocalDataTable
Definition: MeshToVolume.h:1297
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:734
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1293
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:367
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:489
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2648
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1912
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1136
EdgeData operator+(const T &) const
Definition: MeshToVolume.h:422
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:109
_RootNodeType RootNodeType
Definition: Tree.h:209
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:719
Extracts and stores voxel edge intersection data from a mesh.
Definition: MeshToVolume.h:401
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1072
int32_t Int32
Definition: Types.h:60
ExpandNarrowband(std::vector< BoolLeafNodeType * > &maskNodes, BoolTreeType &maskTree, TreeType &distTree, Int32TreeType &indexTree, const MeshDataAdapter &mesh, ValueType exteriorBandWidth, ValueType interiorBandWidth, ValueType voxelSize)
Definition: MeshToVolume.h:2270
std::vector< LeafNodeType * > & updatedDistNodes()
Definition: MeshToVolume.h:2428
Index32 mZPrim
Definition: MeshToVolume.h:433
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3387
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2136
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:635
RestoreOrigin(std::vector< LeafNodeType * > &nodes, const Coord *coordinates)
Definition: MeshToVolume.h:630
TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:531
bool *const mNodeMask
Definition: MeshToVolume.h:1281
float mZDist
Definition: MeshToVolume.h:432
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:3360
math::Vec4< Index32 > Vec4I
Definition: Types.h:90
std::ostream & operator<<(std::ostream &ostr, const MeshToVoxelEdgeData::EdgeData &rhs)
Definition: MeshToVolume.h:3377
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2108
RemoveSelfIntersectingSurface(std::vector< LeafNodeType * > &nodes, TreeType &distTree, Int32TreeType &indexTree)
Definition: MeshToVolume.h:1685
bool processY(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1207
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2724
LeafNodeType **const mNodes
Definition: MeshToVolume.h:620
boost::scoped_ptr< VoxelizationData > Ptr
Definition: MeshToVolume.h:1846
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2609
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2735
bool const *const mChangedNodeMask
Definition: MeshToVolume.h:1127
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2323
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
void getIndexSpacePoint(size_t n, size_t v, Vec3d &pos) const
Returns position pos in local grid index space for polygon n and vertex v.
Definition: MeshToVolume.h:219
static bool check(const ValueType v)
Definition: MeshToVolume.h:1630
TreeType::ValueType ValueType
Definition: MeshToVolume.h:2788
const size_t * offsetsNextX() const
Definition: MeshToVolume.h:757
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1291
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:709
tree::ValueAccessor< TreeType > FloatTreeAcc
Definition: MeshToVolume.h:1852
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2692
PointType *const mPointsOut
Definition: MeshToVolume.h:510
bool *const mChangedNodeMask
Definition: MeshToVolume.h:1280
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:663
UnionValueMasks(std::vector< LeafNodeTypeA * > &nodesA, std::vector< LeafNodeTypeB * > &nodesB)
Definition: MeshToVolume.h:2130
TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1849
LeafNodeConnectivityTable< TreeType > ConnectivityTable
Definition: MeshToVolume.h:781
InactivateValues(std::vector< LeafNodeType * > &nodes, ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2650
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
Definition: Exceptions.h:87
Vec3< float > Vec3s
Definition: Vec3.h:642
unsigned char getNewPrimId()
Definition: MeshToVolume.h:1877
Definition: Math.h:839
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:266
ComputeIntersectingVoxelSign(std::vector< LeafNodeType * > &distNodes, const TreeType &distTree, const Int32TreeType &indexTree, const MeshDataAdapter &mesh)
Definition: MeshToVolume.h:1299
bool *const mChangedVoxelMask
Definition: MeshToVolume.h:1282
TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:2265
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1292
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1743
std::pair< boost::shared_array< Vec3d >, boost::shared_array< bool > > LocalData
Definition: MeshToVolume.h:1296
QuadAndTriangleDataAdapter(const PointType *pointArray, size_t pointArraySize, const PolygonType *polygonArray, size_t polygonArraySize)
Definition: MeshToVolume.h:200
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1672
VoxelizePolygons(DataTable &dataTable, const MeshDataAdapter &mesh, Interrupter *interrupter=NULL)
Definition: MeshToVolume.h:1903
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:212
GenEdgeData(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition: MeshToVolume.h:3441
TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2151
UCharTreeAcc primIdAcc
Definition: MeshToVolume.h:1875
tbb::enumerable_thread_specific< typename VoxelizationDataType::Ptr > DataTable
Definition: MeshToVolume.h:1901
TransformPoints(const PointType *pointsIn, PointType *pointsOut, const math::Transform &xform)
Definition: MeshToVolume.h:483
std::vector< LeafNodeType * > & newDistNodes()
Definition: MeshToVolume.h:2427
Contiguous quad and triangle data adapter class.
Definition: MeshToVolume.h:189
CombineLeafNodes(TreeType &lhsDistTree, Int32TreeType &lhsIdxTree, LeafNodeType **rhsDistNodes, Int32LeafNodeType **rhsIdxNodes)
Definition: MeshToVolume.h:536
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1228
void maskNodeInternalNeighbours(const Index pos, bool(&mask)[26])
Definition: MeshToVolume.h:1481
Internal edge data type.
Definition: MeshToVolume.h:408
Accessor getAccessor()
Definition: MeshToVolume.h:463
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1638
TreeType::ValueType ValueType
Definition: MeshToVolume.h:1044
ComputeNodeConnectivity(const TreeType &tree, const Coord *coordinates, size_t *offsets, size_t numNodes, const CoordBBox &bbox)
Definition: MeshToVolume.h:653