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