OpenVDB  7.0.0
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 <boost/mpl/at.hpp>
39 #include <boost/mpl/int.hpp>
40 #include <boost/mpl/size.hpp>
41 
42 #include <algorithm> // for std::sort()
43 #include <cmath> // for std::isfinite(), std::isnan()
44 #include <deque>
45 #include <limits>
46 #include <memory>
47 #include <sstream>
48 #include <type_traits>
49 #include <vector>
50 
51 namespace openvdb {
53 namespace OPENVDB_VERSION_NAME {
54 namespace tools {
55 
56 
58 
59 
62 
68 
72 
76 
80 };
81 
82 
113 template <typename GridType, typename MeshDataAdapter>
114 inline typename GridType::Ptr
116  const MeshDataAdapter& mesh,
117  const math::Transform& transform,
118  float exteriorBandWidth = 3.0f,
119  float interiorBandWidth = 3.0f,
120  int flags = 0,
121  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = nullptr);
122 
123 
139 template <typename GridType, typename MeshDataAdapter, typename Interrupter>
140 inline typename GridType::Ptr
142  Interrupter& interrupter,
143  const MeshDataAdapter& mesh,
144  const math::Transform& transform,
145  float exteriorBandWidth = 3.0f,
146  float interiorBandWidth = 3.0f,
147  int flags = 0,
148  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = nullptr);
149 
150 
152 
153 
164 template<typename PointType, typename PolygonType>
166 
167  QuadAndTriangleDataAdapter(const std::vector<PointType>& points,
168  const std::vector<PolygonType>& polygons)
169  : mPointArray(points.empty() ? nullptr : &points[0])
170  , mPointArraySize(points.size())
171  , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
172  , mPolygonArraySize(polygons.size())
173  {
174  }
175 
176  QuadAndTriangleDataAdapter(const PointType * pointArray, size_t pointArraySize,
177  const PolygonType* polygonArray, size_t polygonArraySize)
178  : mPointArray(pointArray)
179  , mPointArraySize(pointArraySize)
180  , mPolygonArray(polygonArray)
181  , mPolygonArraySize(polygonArraySize)
182  {
183  }
184 
185  size_t polygonCount() const { return mPolygonArraySize; }
186  size_t pointCount() const { return mPointArraySize; }
187 
189  size_t vertexCount(size_t n) const {
190  return (PolygonType::size == 3 || mPolygonArray[n][3] == util::INVALID_IDX) ? 3 : 4;
191  }
192 
195  void getIndexSpacePoint(size_t n, size_t v, Vec3d& pos) const {
196  const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
197  pos[0] = double(p[0]);
198  pos[1] = double(p[1]);
199  pos[2] = double(p[2]);
200  }
201 
202 private:
203  PointType const * const mPointArray;
204  size_t const mPointArraySize;
205  PolygonType const * const mPolygonArray;
206  size_t const mPolygonArraySize;
207 }; // struct QuadAndTriangleDataAdapter
208 
209 
211 
212 
213 // Convenience functions for the mesh to volume converter that wrap stl containers.
214 //
215 // Note the meshToVolume() method declared above is more flexible and better suited
216 // for arbitrary data structures.
217 
218 
234 template<typename GridType>
235 inline typename GridType::Ptr
237  const openvdb::math::Transform& xform,
238  const std::vector<Vec3s>& points,
239  const std::vector<Vec3I>& triangles,
240  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
241 
243 template<typename GridType, typename Interrupter>
244 inline typename GridType::Ptr
246  Interrupter& interrupter,
247  const openvdb::math::Transform& xform,
248  const std::vector<Vec3s>& points,
249  const std::vector<Vec3I>& triangles,
250  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
251 
252 
268 template<typename GridType>
269 inline typename GridType::Ptr
271  const openvdb::math::Transform& xform,
272  const std::vector<Vec3s>& points,
273  const std::vector<Vec4I>& quads,
274  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
275 
277 template<typename GridType, typename Interrupter>
278 inline typename GridType::Ptr
280  Interrupter& interrupter,
281  const openvdb::math::Transform& xform,
282  const std::vector<Vec3s>& points,
283  const std::vector<Vec4I>& quads,
284  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
285 
286 
303 template<typename GridType>
304 inline typename GridType::Ptr
306  const openvdb::math::Transform& xform,
307  const std::vector<Vec3s>& points,
308  const std::vector<Vec3I>& triangles,
309  const std::vector<Vec4I>& quads,
310  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
311 
313 template<typename GridType, typename Interrupter>
314 inline typename GridType::Ptr
316  Interrupter& interrupter,
317  const openvdb::math::Transform& xform,
318  const std::vector<Vec3s>& points,
319  const std::vector<Vec3I>& triangles,
320  const std::vector<Vec4I>& quads,
321  float halfWidth = float(LEVEL_SET_HALF_WIDTH));
322 
323 
342 template<typename GridType>
343 inline typename GridType::Ptr
345  const openvdb::math::Transform& xform,
346  const std::vector<Vec3s>& points,
347  const std::vector<Vec3I>& triangles,
348  const std::vector<Vec4I>& quads,
349  float exBandWidth,
350  float inBandWidth);
351 
353 template<typename GridType, typename Interrupter>
354 inline typename GridType::Ptr
356  Interrupter& interrupter,
357  const openvdb::math::Transform& xform,
358  const std::vector<Vec3s>& points,
359  const std::vector<Vec3I>& triangles,
360  const std::vector<Vec4I>& quads,
361  float exBandWidth,
362  float inBandWidth);
363 
364 
379 template<typename GridType>
380 inline typename GridType::Ptr
382  const openvdb::math::Transform& xform,
383  const std::vector<Vec3s>& points,
384  const std::vector<Vec3I>& triangles,
385  const std::vector<Vec4I>& quads,
386  float bandWidth);
387 
389 template<typename GridType, typename Interrupter>
390 inline typename GridType::Ptr
392  Interrupter& interrupter,
393  const openvdb::math::Transform& xform,
394  const std::vector<Vec3s>& points,
395  const std::vector<Vec3I>& triangles,
396  const std::vector<Vec4I>& quads,
397  float bandWidth);
398 
399 
401 
402 
409 template<typename GridType, typename VecType>
410 inline typename GridType::Ptr
412  const openvdb::math::Transform& xform,
413  typename VecType::ValueType halfWidth = LEVEL_SET_HALF_WIDTH);
414 
415 
417 
418 
425 template <typename FloatTreeT>
426 inline void
427 traceExteriorBoundaries(FloatTreeT& tree);
428 
429 
431 
432 
435 {
436 public:
437 
439 
441  struct EdgeData {
442  EdgeData(float dist = 1.0)
443  : mXDist(dist), mYDist(dist), mZDist(dist)
444  , mXPrim(util::INVALID_IDX)
445  , mYPrim(util::INVALID_IDX)
446  , mZPrim(util::INVALID_IDX)
447  {
448  }
449 
451  bool operator< (const EdgeData&) const { return false; }
454  bool operator> (const EdgeData&) const { return false; }
455  template<class T> EdgeData operator+(const T&) const { return *this; }
456  template<class T> EdgeData operator-(const T&) const { return *this; }
457  EdgeData operator-() const { return *this; }
459 
460  bool operator==(const EdgeData& rhs) const
461  {
462  return mXPrim == rhs.mXPrim && mYPrim == rhs.mYPrim && mZPrim == rhs.mZPrim;
463  }
464 
465  float mXDist, mYDist, mZDist;
466  Index32 mXPrim, mYPrim, mZPrim;
467  };
468 
471 
472 
474 
475 
477 
478 
486  void convert(const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
487 
488 
491  void getEdgeData(Accessor& acc, const Coord& ijk,
492  std::vector<Vec3d>& points, std::vector<Index32>& primitives);
493 
496  Accessor getAccessor() { return Accessor(mTree); }
497 
498 private:
499  void operator=(const MeshToVoxelEdgeData&) {}
500  TreeType mTree;
501  class GenEdgeData;
502 };
503 
504 
507 
508 
509 // Internal utility objects and implementation details
510 
511 namespace mesh_to_volume_internal {
512 
513 template<typename PointType>
515 
516  TransformPoints(const PointType* pointsIn, PointType* pointsOut,
517  const math::Transform& xform)
518  : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
519  {
520  }
521 
522  void operator()(const tbb::blocked_range<size_t>& range) const {
523 
524  Vec3d pos;
525 
526  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
527 
528  const PointType& wsP = mPointsIn[n];
529  pos[0] = double(wsP[0]);
530  pos[1] = double(wsP[1]);
531  pos[2] = double(wsP[2]);
532 
533  pos = mXform->worldToIndex(pos);
534 
535  PointType& isP = mPointsOut[n];
536  isP[0] = typename PointType::value_type(pos[0]);
537  isP[1] = typename PointType::value_type(pos[1]);
538  isP[2] = typename PointType::value_type(pos[2]);
539  }
540  }
541 
542  PointType const * const mPointsIn;
543  PointType * const mPointsOut;
544  math::Transform const * const mXform;
545 }; // TransformPoints
546 
547 
548 template<typename ValueType>
549 struct Tolerance
550 {
551  static ValueType epsilon() { return ValueType(1e-7); }
552  static ValueType minNarrowBandWidth() { return ValueType(1.0 + 1e-6); }
553 };
554 
555 
557 
558 
559 template<typename TreeType>
561 {
562 public:
563 
564  using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
565 
566  using LeafNodeType = typename TreeType::LeafNodeType;
567  using Int32LeafNodeType = typename Int32TreeType::LeafNodeType;
568 
569  CombineLeafNodes(TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
570  LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
571  : mDistTree(&lhsDistTree)
572  , mIdxTree(&lhsIdxTree)
573  , mRhsDistNodes(rhsDistNodes)
574  , mRhsIdxNodes(rhsIdxNodes)
575  {
576  }
577 
578  void operator()(const tbb::blocked_range<size_t>& range) const {
579 
580  tree::ValueAccessor<TreeType> distAcc(*mDistTree);
581  tree::ValueAccessor<Int32TreeType> idxAcc(*mIdxTree);
582 
583  using DistValueType = typename LeafNodeType::ValueType;
584  using IndexValueType = typename Int32LeafNodeType::ValueType;
585 
586  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
587 
588  const Coord& origin = mRhsDistNodes[n]->origin();
589 
590  LeafNodeType* lhsDistNode = distAcc.probeLeaf(origin);
591  Int32LeafNodeType* lhsIdxNode = idxAcc.probeLeaf(origin);
592 
593  DistValueType* lhsDistData = lhsDistNode->buffer().data();
594  IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
595 
596  const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
597  const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
598 
599 
600  for (Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
601 
602  if (rhsIdxData[offset] != Int32(util::INVALID_IDX)) {
603 
604  const DistValueType& lhsValue = lhsDistData[offset];
605  const DistValueType& rhsValue = rhsDistData[offset];
606 
607  if (rhsValue < lhsValue) {
608  lhsDistNode->setValueOn(offset, rhsValue);
609  lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
610  } else if (math::isExactlyEqual(rhsValue, lhsValue)) {
611  lhsIdxNode->setValueOn(offset,
612  std::min(lhsIdxData[offset], rhsIdxData[offset]));
613  }
614  }
615  }
616 
617  delete mRhsDistNodes[n];
618  delete mRhsIdxNodes[n];
619  }
620  }
621 
622 private:
623 
624  TreeType * const mDistTree;
625  Int32TreeType * const mIdxTree;
626 
627  LeafNodeType ** const mRhsDistNodes;
628  Int32LeafNodeType ** const mRhsIdxNodes;
629 }; // class CombineLeafNodes
630 
631 
633 
634 
635 template<typename TreeType>
637 {
638  using LeafNodeType = typename TreeType::LeafNodeType;
639 
640  StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
641  : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
642  {
643  }
644 
645  void operator()(const tbb::blocked_range<size_t>& range) const {
646  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
647  Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
648  mCoordinates[n] = origin;
649  origin[0] = static_cast<int>(n);
650  }
651  }
652 
654  Coord * const mCoordinates;
655 };
656 
657 
658 template<typename TreeType>
660 {
661  using LeafNodeType = typename TreeType::LeafNodeType;
662 
663  RestoreOrigin(std::vector<LeafNodeType*>& nodes, const Coord* coordinates)
664  : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
665  {
666  }
667 
668  void operator()(const tbb::blocked_range<size_t>& range) const {
669  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
670  Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
671  origin[0] = mCoordinates[n][0];
672  }
673  }
674 
676  Coord const * const mCoordinates;
677 };
678 
679 
680 template<typename TreeType>
682 {
683 public:
684  using LeafNodeType = typename TreeType::LeafNodeType;
685 
686  ComputeNodeConnectivity(const TreeType& tree, const Coord* coordinates,
687  size_t* offsets, size_t numNodes, const CoordBBox& bbox)
688  : mTree(&tree)
689  , mCoordinates(coordinates)
690  , mOffsets(offsets)
691  , mNumNodes(numNodes)
692  , mBBox(bbox)
693  {
694  }
695 
697 
698  // Disallow assignment
699  ComputeNodeConnectivity& operator=(const ComputeNodeConnectivity&) = delete;
700 
701  void operator()(const tbb::blocked_range<size_t>& range) const {
702 
703  size_t* offsetsNextX = mOffsets;
704  size_t* offsetsPrevX = mOffsets + mNumNodes;
705  size_t* offsetsNextY = mOffsets + mNumNodes * 2;
706  size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
707  size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
708  size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
709 
711  Coord ijk;
712 
713  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
714  const Coord& origin = mCoordinates[n];
715  offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(LeafNodeType::DIM, 0, 0));
716  offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-LeafNodeType::DIM, 0, 0));
717  offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, LeafNodeType::DIM, 0));
718  offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -LeafNodeType::DIM, 0));
719  offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, LeafNodeType::DIM));
720  offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -LeafNodeType::DIM));
721  }
722  }
723 
725  const Coord& start, const Coord& step) const
726  {
727  Coord ijk = start + step;
728  CoordBBox bbox(mBBox);
729 
730  while (bbox.isInside(ijk)) {
731  const LeafNodeType* node = acc.probeConstLeaf(ijk);
732  if (node) return static_cast<size_t>(node->origin()[0]);
733  ijk += step;
734  }
735 
737  }
738 
739 
740 private:
741  TreeType const * const mTree;
742  Coord const * const mCoordinates;
743  size_t * const mOffsets;
744 
745  const size_t mNumNodes;
746  const CoordBBox mBBox;
747 }; // class ComputeNodeConnectivity
748 
749 
750 template<typename TreeType>
752 {
753  enum { INVALID_OFFSET = std::numeric_limits<size_t>::max() };
754 
755  using LeafNodeType = typename TreeType::LeafNodeType;
756 
758  {
759  mLeafNodes.reserve(tree.leafCount());
760  tree.getNodes(mLeafNodes);
761 
762  if (mLeafNodes.empty()) return;
763 
764  CoordBBox bbox;
765  tree.evalLeafBoundingBox(bbox);
766 
767  const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
768 
769  // stash the leafnode origin coordinate and temporarily store the
770  // linear offset in the origin.x variable.
771  std::unique_ptr<Coord[]> coordinates{new Coord[mLeafNodes.size()]};
772  tbb::parallel_for(range,
773  StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
774 
775  // build the leafnode offset table
776  mOffsets.reset(new size_t[mLeafNodes.size() * 6]);
777 
778 
779  tbb::parallel_for(range, ComputeNodeConnectivity<TreeType>(
780  tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
781 
782  // restore the leafnode origin coordinate
783  tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
784  }
785 
786  size_t size() const { return mLeafNodes.size(); }
787 
788  std::vector<LeafNodeType*>& nodes() { return mLeafNodes; }
789  const std::vector<LeafNodeType*>& nodes() const { return mLeafNodes; }
790 
791 
792  const size_t* offsetsNextX() const { return mOffsets.get(); }
793  const size_t* offsetsPrevX() const { return mOffsets.get() + mLeafNodes.size(); }
794 
795  const size_t* offsetsNextY() const { return mOffsets.get() + mLeafNodes.size() * 2; }
796  const size_t* offsetsPrevY() const { return mOffsets.get() + mLeafNodes.size() * 3; }
797 
798  const size_t* offsetsNextZ() const { return mOffsets.get() + mLeafNodes.size() * 4; }
799  const size_t* offsetsPrevZ() const { return mOffsets.get() + mLeafNodes.size() * 5; }
800 
801 private:
802  std::vector<LeafNodeType*> mLeafNodes;
803  std::unique_ptr<size_t[]> mOffsets;
804 }; // struct LeafNodeConnectivityTable
805 
806 
807 template<typename TreeType>
809 {
810 public:
811 
812  enum Axis { X_AXIS = 0, Y_AXIS = 1, Z_AXIS = 2 };
813 
814  using ValueType = typename TreeType::ValueType;
815  using LeafNodeType = typename TreeType::LeafNodeType;
817 
818  SweepExteriorSign(Axis axis, const std::vector<size_t>& startNodeIndices,
819  ConnectivityTable& connectivity)
820  : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
821  , mConnectivity(&connectivity)
822  , mAxis(axis)
823  {
824  }
825 
826  void operator()(const tbb::blocked_range<size_t>& range) const {
827 
828  std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
829 
830  // Z Axis
831  size_t idxA = 0, idxB = 1;
832  Index step = 1;
833 
834  const size_t* nextOffsets = mConnectivity->offsetsNextZ();
835  const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
836 
837  if (mAxis == Y_AXIS) {
838 
839  idxA = 0;
840  idxB = 2;
841  step = LeafNodeType::DIM;
842 
843  nextOffsets = mConnectivity->offsetsNextY();
844  prevOffsets = mConnectivity->offsetsPrevY();
845 
846  } else if (mAxis == X_AXIS) {
847 
848  idxA = 1;
849  idxB = 2;
850  step = LeafNodeType::DIM * LeafNodeType::DIM;
851 
852  nextOffsets = mConnectivity->offsetsNextX();
853  prevOffsets = mConnectivity->offsetsPrevX();
854  }
855 
856  Coord ijk(0, 0, 0);
857 
858  int& a = ijk[idxA];
859  int& b = ijk[idxB];
860 
861  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
862 
863  size_t startOffset = mStartNodeIndices[n];
864  size_t lastOffset = startOffset;
865 
866  Index pos(0);
867 
868  for (a = 0; a < int(LeafNodeType::DIM); ++a) {
869  for (b = 0; b < int(LeafNodeType::DIM); ++b) {
870 
871  pos = LeafNodeType::coordToOffset(ijk);
872  size_t offset = startOffset;
873 
874  // sweep in +axis direction until a boundary voxel is hit.
875  while ( offset != ConnectivityTable::INVALID_OFFSET &&
876  traceVoxelLine(*nodes[offset], pos, step) ) {
877 
878  lastOffset = offset;
879  offset = nextOffsets[offset];
880  }
881 
882  // find last leafnode in +axis direction
883  offset = lastOffset;
884  while (offset != ConnectivityTable::INVALID_OFFSET) {
885  lastOffset = offset;
886  offset = nextOffsets[offset];
887  }
888 
889  // sweep in -axis direction until a boundary voxel is hit.
890  offset = lastOffset;
891  pos += step * (LeafNodeType::DIM - 1);
892  while ( offset != ConnectivityTable::INVALID_OFFSET &&
893  traceVoxelLine(*nodes[offset], pos, -step)) {
894  offset = prevOffsets[offset];
895  }
896  }
897  }
898  }
899  }
900 
901 
902  bool traceVoxelLine(LeafNodeType& node, Index pos, Index step) const {
903 
904  ValueType* data = node.buffer().data();
905 
906  bool isOutside = true;
907 
908  for (Index i = 0; i < LeafNodeType::DIM; ++i) {
909 
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, 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  mChangedNodeMask[n] = scanFill(*mNodes[n]);
1100  }
1101  }
1102  }
1103 
1105  bool * const mChangedNodeMask;
1106 };
1107 
1108 
1109 template<typename ValueType>
1111 {
1112  FillArray(ValueType* array, const ValueType v) : mArray(array), mValue(v) { }
1113 
1114  void operator()(const tbb::blocked_range<size_t>& range) const {
1115  const ValueType v = mValue;
1116  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1117  mArray[n] = v;
1118  }
1119  }
1120 
1121  ValueType * const mArray;
1122  const ValueType mValue;
1123 };
1124 
1125 
1126 template<typename ValueType>
1127 inline void
1128 fillArray(ValueType* array, const ValueType val, const size_t length)
1129 {
1130  const auto grainSize = std::max<size_t>(
1131  length / tbb::task_scheduler_init::default_num_threads(), 1024);
1132  const tbb::blocked_range<size_t> range(0, length, grainSize);
1133  tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1134 }
1135 
1136 
1137 template<typename TreeType>
1139 {
1140 public:
1141  using ValueType = typename TreeType::ValueType;
1142  using LeafNodeType = typename TreeType::LeafNodeType;
1143 
1144  SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
1145  const bool* changedNodeMask, bool* changedVoxelMask)
1146  : mNodes(nodes.empty() ? nullptr : &nodes[0])
1147  , mChangedNodeMask(changedNodeMask)
1148  , mChangedVoxelMask(changedVoxelMask)
1149  {
1150  }
1151 
1152  void operator()(const tbb::blocked_range<size_t>& range) const {
1153  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1154 
1155  if (mChangedNodeMask[n]) {
1156  bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1157 
1158  ValueType* data = mNodes[n]->buffer().data();
1159 
1160  for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1161  if (mask[pos]) {
1162  data[pos] = ValueType(-data[pos]);
1163  mask[pos] = false;
1164  }
1165  }
1166  }
1167  }
1168  }
1169 
1171  bool const * const mChangedNodeMask;
1172  bool * const mChangedVoxelMask;
1173 };
1174 
1175 
1176 template<typename TreeType>
1178 {
1179 public:
1180  using ValueType = typename TreeType::ValueType;
1181  using LeafNodeType = typename TreeType::LeafNodeType;
1183 
1185  bool* changedNodeMask, bool* nodeMask, bool* changedVoxelMask)
1186  : mConnectivity(&connectivity)
1187  , mChangedNodeMask(changedNodeMask)
1188  , mNodeMask(nodeMask)
1189  , mChangedVoxelMask(changedVoxelMask)
1190  {
1191  }
1192 
1193  void operator()(const tbb::blocked_range<size_t>& range) const {
1194 
1195  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1196 
1197  if (!mChangedNodeMask[n]) {
1198 
1199  bool changedValue = false;
1200 
1201  changedValue |= processZ(n, /*firstFace=*/true);
1202  changedValue |= processZ(n, /*firstFace=*/false);
1203 
1204  changedValue |= processY(n, /*firstFace=*/true);
1205  changedValue |= processY(n, /*firstFace=*/false);
1206 
1207  changedValue |= processX(n, /*firstFace=*/true);
1208  changedValue |= processX(n, /*firstFace=*/false);
1209 
1210  mNodeMask[n] = changedValue;
1211  }
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 boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
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  ConnectivityTable nodeConnectivity(tree);
3024 
3025  std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3026 
3027  for (size_t n = 0; n < nodeConnectivity.size(); ++n) {
3028  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3029  xStartNodes.push_back(n);
3030  }
3031 
3032  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3033  yStartNodes.push_back(n);
3034  }
3035 
3036  if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3037  zStartNodes.push_back(n);
3038  }
3039  }
3040 
3042 
3043  tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3044  SweepingOp(SweepingOp::Z_AXIS, zStartNodes, nodeConnectivity));
3045 
3046  tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3047  SweepingOp(SweepingOp::Y_AXIS, yStartNodes, nodeConnectivity));
3048 
3049  tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3050  SweepingOp(SweepingOp::X_AXIS, xStartNodes, nodeConnectivity));
3051 
3052  const size_t numLeafNodes = nodeConnectivity.size();
3053  const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3054 
3055  std::unique_ptr<bool[]> changedNodeMaskA{new bool[numLeafNodes]};
3056  std::unique_ptr<bool[]> changedNodeMaskB{new bool[numLeafNodes]};
3057  std::unique_ptr<bool[]> changedVoxelMask{new bool[numVoxels]};
3058 
3059  mesh_to_volume_internal::fillArray(changedNodeMaskA.get(), true, numLeafNodes);
3060  mesh_to_volume_internal::fillArray(changedNodeMaskB.get(), false, numLeafNodes);
3061  mesh_to_volume_internal::fillArray(changedVoxelMask.get(), false, numVoxels);
3062 
3063  const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3064 
3065  bool nodesUpdated = false;
3066  do {
3067  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
3068  nodeConnectivity.nodes(), changedNodeMaskA.get()));
3069 
3070  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedPoints<FloatTreeT>(
3071  nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3072  changedVoxelMask.get()));
3073 
3074  changedNodeMaskA.swap(changedNodeMaskB);
3075 
3076  nodesUpdated = false;
3077  for (size_t n = 0; n < numLeafNodes; ++n) {
3078  nodesUpdated |= changedNodeMaskA[n];
3079  if (nodesUpdated) break;
3080  }
3081 
3082  if (nodesUpdated) {
3083  tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
3084  nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3085  }
3086  } while (nodesUpdated);
3087 
3088 } // void traceExteriorBoundaries()
3089 
3090 
3092 
3093 
3094 template <typename GridType, typename MeshDataAdapter, typename Interrupter>
3095 inline typename GridType::Ptr
3097  Interrupter& interrupter,
3098  const MeshDataAdapter& mesh,
3099  const math::Transform& transform,
3100  float exteriorBandWidth,
3101  float interiorBandWidth,
3102  int flags,
3103  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid)
3104 {
3105  using GridTypePtr = typename GridType::Ptr;
3106  using TreeType = typename GridType::TreeType;
3107  using LeafNodeType = typename TreeType::LeafNodeType;
3108  using ValueType = typename GridType::ValueType;
3109 
3110  using Int32GridType = typename GridType::template ValueConverter<Int32>::Type;
3111  using Int32TreeType = typename Int32GridType::TreeType;
3112 
3113  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
3114 
3116 
3117  // Setup
3118 
3119  GridTypePtr distGrid(new GridType(std::numeric_limits<ValueType>::max()));
3120  distGrid->setTransform(transform.copy());
3121 
3122  ValueType exteriorWidth = ValueType(exteriorBandWidth);
3123  ValueType interiorWidth = ValueType(interiorBandWidth);
3124 
3125  // Note: inf interior width is all right, this value makes the converter fill
3126  // interior regions with distance values.
3127  if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3128  std::stringstream msg;
3129  msg << "Illegal narrow band width: exterior = " << exteriorWidth
3130  << ", interior = " << interiorWidth;
3131  OPENVDB_LOG_DEBUG(msg.str());
3132  return distGrid;
3133  }
3134 
3135  const ValueType voxelSize = ValueType(transform.voxelSize()[0]);
3136 
3137  if (!std::isfinite(voxelSize) || math::isZero(voxelSize)) {
3138  std::stringstream msg;
3139  msg << "Illegal transform, voxel size = " << voxelSize;
3140  OPENVDB_LOG_DEBUG(msg.str());
3141  return distGrid;
3142  }
3143 
3144  // Convert narrow band width from voxel units to world space units.
3145  exteriorWidth *= voxelSize;
3146  // Avoid the unit conversion if the interior band width is set to
3147  // inf or std::numeric_limits<float>::max().
3148  if (interiorWidth < std::numeric_limits<ValueType>::max()) {
3149  interiorWidth *= voxelSize;
3150  }
3151 
3152  const bool computeSignedDistanceField = (flags & UNSIGNED_DISTANCE_FIELD) == 0;
3153  const bool removeIntersectingVoxels = (flags & DISABLE_INTERSECTING_VOXEL_REMOVAL) == 0;
3154  const bool renormalizeValues = (flags & DISABLE_RENORMALIZATION) == 0;
3155  const bool trimNarrowBand = (flags & DISABLE_NARROW_BAND_TRIMMING) == 0;
3156 
3157  Int32GridType* indexGrid = nullptr;
3158 
3159  typename Int32GridType::Ptr temporaryIndexGrid;
3160 
3161  if (polygonIndexGrid) {
3162  indexGrid = polygonIndexGrid;
3163  } else {
3164  temporaryIndexGrid.reset(new Int32GridType(Int32(util::INVALID_IDX)));
3165  indexGrid = temporaryIndexGrid.get();
3166  }
3167 
3168  indexGrid->newTree();
3169  indexGrid->setTransform(transform.copy());
3170 
3171  if (computeSignedDistanceField) {
3172  distGrid->setGridClass(GRID_LEVEL_SET);
3173  } else {
3174  distGrid->setGridClass(GRID_UNKNOWN);
3175  interiorWidth = ValueType(0.0);
3176  }
3177 
3178  TreeType& distTree = distGrid->tree();
3179  Int32TreeType& indexTree = indexGrid->tree();
3180 
3181 
3183 
3184  // Voxelize mesh
3185 
3186  {
3187  using VoxelizationDataType = mesh_to_volume_internal::VoxelizationData<TreeType>;
3188  using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3189 
3190  DataTable data;
3191  using Voxelizer =
3193 
3194  const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3195 
3196  tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3197 
3198  for (typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3199  VoxelizationDataType& dataItem = **i;
3201  distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3202  }
3203  }
3204 
3205  // The progress estimates are based on the observed average time for a few different
3206  // test cases and is only intended to provide some rough progression feedback to the user.
3207  if (interrupter.wasInterrupted(30)) return distGrid;
3208 
3209 
3211 
3212  // Classify interior and exterior regions
3213 
3214  if (computeSignedDistanceField) {
3215 
3216  // Determines the inside/outside state for the narrow band of voxels.
3217  traceExteriorBoundaries(distTree);
3218 
3219  std::vector<LeafNodeType*> nodes;
3220  nodes.reserve(distTree.leafCount());
3221  distTree.getNodes(nodes);
3222 
3223  const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3224 
3225  using SignOp =
3227 
3228  tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3229 
3230  if (interrupter.wasInterrupted(45)) return distGrid;
3231 
3232  // Remove voxels created by self intersecting portions of the mesh.
3233  if (removeIntersectingVoxels) {
3234 
3235  tbb::parallel_for(nodeRange,
3237 
3238  tbb::parallel_for(nodeRange,
3240  nodes, distTree, indexTree));
3241 
3242  tools::pruneInactive(distTree, /*threading=*/true);
3243  tools::pruneInactive(indexTree, /*threading=*/true);
3244  }
3245  }
3246 
3247  if (interrupter.wasInterrupted(50)) return distGrid;
3248 
3249  if (distTree.activeVoxelCount() == 0) {
3250  distTree.clear();
3251  distTree.root().setBackground(exteriorWidth, /*updateChildNodes=*/false);
3252  return distGrid;
3253  }
3254 
3255  // Transform values (world space scaling etc.).
3256  {
3257  std::vector<LeafNodeType*> nodes;
3258  nodes.reserve(distTree.leafCount());
3259  distTree.getNodes(nodes);
3260 
3261  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3263  nodes, voxelSize, !computeSignedDistanceField));
3264  }
3265 
3266  // Propagate sign information into tile regions.
3267  if (computeSignedDistanceField) {
3268  distTree.root().setBackground(exteriorWidth, /*updateChildNodes=*/false);
3269  tools::signedFloodFillWithValues(distTree, exteriorWidth, -interiorWidth);
3270  } else {
3271  tools::changeBackground(distTree, exteriorWidth);
3272  }
3273 
3274  if (interrupter.wasInterrupted(54)) return distGrid;
3275 
3276 
3278 
3279  // Expand the narrow band region
3280 
3281  const ValueType minBandWidth = voxelSize * ValueType(2.0);
3282 
3283  if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3284 
3285  // Create the initial voxel mask.
3286  BoolTreeType maskTree(false);
3287 
3288  {
3289  std::vector<LeafNodeType*> nodes;
3290  nodes.reserve(distTree.leafCount());
3291  distTree.getNodes(nodes);
3292 
3293  mesh_to_volume_internal::ConstructVoxelMask<TreeType> op(maskTree, distTree, nodes);
3294  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3295  }
3296 
3297  // Progress estimation
3298  unsigned maxIterations = std::numeric_limits<unsigned>::max();
3299 
3300  float progress = 54.0f, step = 0.0f;
3301  double estimated =
3302  2.0 * std::ceil((std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3303 
3304  if (estimated < double(maxIterations)) {
3305  maxIterations = unsigned(estimated);
3306  step = 40.0f / float(maxIterations);
3307  }
3308 
3309  std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3310 
3311  unsigned count = 0;
3312  while (true) {
3313 
3314  if (interrupter.wasInterrupted(int(progress))) return distGrid;
3315 
3316  const size_t maskNodeCount = maskTree.leafCount();
3317  if (maskNodeCount == 0) break;
3318 
3319  maskNodes.clear();
3320  maskNodes.reserve(maskNodeCount);
3321  maskTree.getNodes(maskNodes);
3322 
3323  const tbb::blocked_range<size_t> range(0, maskNodes.size());
3324 
3325  tbb::parallel_for(range,
3327 
3328  mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3329  mesh, exteriorWidth, interiorWidth, voxelSize);
3330 
3331  if ((++count) >= maxIterations) break;
3332  progress += step;
3333  }
3334  }
3335 
3336  if (interrupter.wasInterrupted(94)) return distGrid;
3337 
3338  if (!polygonIndexGrid) indexGrid->clear();
3339 
3340 
3342 
3343  // Renormalize distances to smooth out bumps caused by self intersecting
3344  // and overlapping portions of the mesh and renormalize the level set.
3345 
3346  if (computeSignedDistanceField && renormalizeValues) {
3347 
3348  std::vector<LeafNodeType*> nodes;
3349  nodes.reserve(distTree.leafCount());
3350  distTree.getNodes(nodes);
3351 
3352  std::unique_ptr<ValueType[]> buffer{new ValueType[LeafNodeType::SIZE * nodes.size()]};
3353 
3354  const ValueType offset = ValueType(0.8 * voxelSize);
3355 
3356  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3358 
3359  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3361  distTree, nodes, buffer.get(), voxelSize));
3362 
3363  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3364  mesh_to_volume_internal::MinCombine<TreeType>(nodes, buffer.get()));
3365 
3366  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3369  }
3370 
3371  if (interrupter.wasInterrupted(99)) return distGrid;
3372 
3373 
3375 
3376  // Remove active voxels that exceed the narrow band limits
3377 
3378  if (trimNarrowBand && std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3379 
3380  std::vector<LeafNodeType*> nodes;
3381  nodes.reserve(distTree.leafCount());
3382  distTree.getNodes(nodes);
3383 
3384  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3386  nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3387 
3389  distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3390  }
3391 
3392  return distGrid;
3393 }
3394 
3395 
3396 template <typename GridType, typename MeshDataAdapter>
3397 inline typename GridType::Ptr
3399  const MeshDataAdapter& mesh,
3400  const math::Transform& transform,
3401  float exteriorBandWidth,
3402  float interiorBandWidth,
3403  int flags,
3404  typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid)
3405 {
3406  util::NullInterrupter nullInterrupter;
3407  return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3408  exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3409 }
3410 
3411 
3413 
3414 
3415 //{
3417 
3419 template<typename GridType, typename Interrupter>
3420 inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3421  typename GridType::Ptr>::type
3422 doMeshConversion(
3423  Interrupter& interrupter,
3424  const openvdb::math::Transform& xform,
3425  const std::vector<Vec3s>& points,
3426  const std::vector<Vec3I>& triangles,
3427  const std::vector<Vec4I>& quads,
3428  float exBandWidth,
3429  float inBandWidth,
3430  bool unsignedDistanceField = false)
3431 {
3432  if (points.empty()) {
3433  return typename GridType::Ptr(new GridType(typename GridType::ValueType(exBandWidth)));
3434  }
3435 
3436  const size_t numPoints = points.size();
3437  std::unique_ptr<Vec3s[]> indexSpacePoints{new Vec3s[numPoints]};
3438 
3439  // transform points to local grid index space
3440  tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3442  &points[0], indexSpacePoints.get(), xform));
3443 
3444  const int conversionFlags = unsignedDistanceField ? UNSIGNED_DISTANCE_FIELD : 0;
3445 
3446  if (quads.empty()) {
3447 
3449  mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3450 
3451  return meshToVolume<GridType>(
3452  interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3453 
3454  } else if (triangles.empty()) {
3455 
3457  mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3458 
3459  return meshToVolume<GridType>(
3460  interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3461  }
3462 
3463  // pack primitives
3464 
3465  const size_t numPrimitives = triangles.size() + quads.size();
3466  std::unique_ptr<Vec4I[]> prims{new Vec4I[numPrimitives]};
3467 
3468  for (size_t n = 0, N = triangles.size(); n < N; ++n) {
3469  const Vec3I& triangle = triangles[n];
3470  Vec4I& prim = prims[n];
3471  prim[0] = triangle[0];
3472  prim[1] = triangle[1];
3473  prim[2] = triangle[2];
3474  prim[3] = util::INVALID_IDX;
3475  }
3476 
3477  const size_t offset = triangles.size();
3478  for (size_t n = 0, N = quads.size(); n < N; ++n) {
3479  prims[offset + n] = quads[n];
3480  }
3481 
3483  mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3484 
3485  return meshToVolume<GridType>(interrupter, mesh, xform,
3486  exBandWidth, inBandWidth, conversionFlags);
3487 }
3488 
3489 
3492 template<typename GridType, typename Interrupter>
3493 inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3494  typename GridType::Ptr>::type
3495 doMeshConversion(
3496  Interrupter&,
3497  const math::Transform& /*xform*/,
3498  const std::vector<Vec3s>& /*points*/,
3499  const std::vector<Vec3I>& /*triangles*/,
3500  const std::vector<Vec4I>& /*quads*/,
3501  float /*exBandWidth*/,
3502  float /*inBandWidth*/,
3503  bool /*unsignedDistanceField*/ = false)
3504 {
3506  "mesh to volume conversion is supported only for scalar floating-point grids");
3507 }
3508 
3510 //}
3511 
3512 
3514 
3515 
3516 template<typename GridType>
3517 inline typename GridType::Ptr
3519  const openvdb::math::Transform& xform,
3520  const std::vector<Vec3s>& points,
3521  const std::vector<Vec3I>& triangles,
3522  float halfWidth)
3523 {
3524  util::NullInterrupter nullInterrupter;
3525  std::vector<Vec4I> quads(0);
3526  return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3527  halfWidth, halfWidth);
3528 }
3529 
3530 
3531 template<typename GridType, typename Interrupter>
3532 inline typename GridType::Ptr
3534  Interrupter& interrupter,
3535  const openvdb::math::Transform& xform,
3536  const std::vector<Vec3s>& points,
3537  const std::vector<Vec3I>& triangles,
3538  float halfWidth)
3539 {
3540  std::vector<Vec4I> quads(0);
3541  return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3542  halfWidth, halfWidth);
3543 }
3544 
3545 
3546 template<typename GridType>
3547 inline typename GridType::Ptr
3549  const openvdb::math::Transform& xform,
3550  const std::vector<Vec3s>& points,
3551  const std::vector<Vec4I>& quads,
3552  float halfWidth)
3553 {
3554  util::NullInterrupter nullInterrupter;
3555  std::vector<Vec3I> triangles(0);
3556  return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3557  halfWidth, halfWidth);
3558 }
3559 
3560 
3561 template<typename GridType, typename Interrupter>
3562 inline typename GridType::Ptr
3564  Interrupter& interrupter,
3565  const openvdb::math::Transform& xform,
3566  const std::vector<Vec3s>& points,
3567  const std::vector<Vec4I>& quads,
3568  float halfWidth)
3569 {
3570  std::vector<Vec3I> triangles(0);
3571  return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3572  halfWidth, halfWidth);
3573 }
3574 
3575 
3576 template<typename GridType>
3577 inline typename GridType::Ptr
3579  const openvdb::math::Transform& xform,
3580  const std::vector<Vec3s>& points,
3581  const std::vector<Vec3I>& triangles,
3582  const std::vector<Vec4I>& quads,
3583  float halfWidth)
3584 {
3585  util::NullInterrupter nullInterrupter;
3586  return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3587  halfWidth, halfWidth);
3588 }
3589 
3590 
3591 template<typename GridType, typename Interrupter>
3592 inline typename GridType::Ptr
3594  Interrupter& interrupter,
3595  const openvdb::math::Transform& xform,
3596  const std::vector<Vec3s>& points,
3597  const std::vector<Vec3I>& triangles,
3598  const std::vector<Vec4I>& quads,
3599  float halfWidth)
3600 {
3601  return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3602  halfWidth, halfWidth);
3603 }
3604 
3605 
3606 template<typename GridType>
3607 inline typename GridType::Ptr
3609  const openvdb::math::Transform& xform,
3610  const std::vector<Vec3s>& points,
3611  const std::vector<Vec3I>& triangles,
3612  const std::vector<Vec4I>& quads,
3613  float exBandWidth,
3614  float inBandWidth)
3615 {
3616  util::NullInterrupter nullInterrupter;
3617  return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles,
3618  quads, exBandWidth, inBandWidth);
3619 }
3620 
3621 
3622 template<typename GridType, typename Interrupter>
3623 inline typename GridType::Ptr
3625  Interrupter& interrupter,
3626  const openvdb::math::Transform& xform,
3627  const std::vector<Vec3s>& points,
3628  const std::vector<Vec3I>& triangles,
3629  const std::vector<Vec4I>& quads,
3630  float exBandWidth,
3631  float inBandWidth)
3632 {
3633  return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3634  quads, exBandWidth, inBandWidth);
3635 }
3636 
3637 
3638 template<typename GridType>
3639 inline typename GridType::Ptr
3641  const openvdb::math::Transform& xform,
3642  const std::vector<Vec3s>& points,
3643  const std::vector<Vec3I>& triangles,
3644  const std::vector<Vec4I>& quads,
3645  float bandWidth)
3646 {
3647  util::NullInterrupter nullInterrupter;
3648  return doMeshConversion<GridType>(nullInterrupter, xform, points, triangles, quads,
3649  bandWidth, bandWidth, true);
3650 }
3651 
3652 
3653 template<typename GridType, typename Interrupter>
3654 inline typename GridType::Ptr
3656  Interrupter& interrupter,
3657  const openvdb::math::Transform& xform,
3658  const std::vector<Vec3s>& points,
3659  const std::vector<Vec3I>& triangles,
3660  const std::vector<Vec4I>& quads,
3661  float bandWidth)
3662 {
3663  return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3664  bandWidth, bandWidth, true);
3665 }
3666 
3667 
3669 
3670 
3671 // Required by several of the tree nodes
3672 inline std::ostream&
3673 operator<<(std::ostream& ostr, const MeshToVoxelEdgeData::EdgeData& rhs)
3674 {
3675  ostr << "{[ " << rhs.mXPrim << ", " << rhs.mXDist << "]";
3676  ostr << " [ " << rhs.mYPrim << ", " << rhs.mYDist << "]";
3677  ostr << " [ " << rhs.mZPrim << ", " << rhs.mZDist << "]}";
3678  return ostr;
3679 }
3680 
3681 // Required by math::Abs
3684 {
3685  return x;
3686 }
3687 
3688 
3690 
3691 
3693 {
3694 public:
3695 
3696  GenEdgeData(
3697  const std::vector<Vec3s>& pointList,
3698  const std::vector<Vec4I>& polygonList);
3699 
3700  void run(bool threaded = true);
3701 
3702  GenEdgeData(GenEdgeData& rhs, tbb::split);
3703  inline void operator() (const tbb::blocked_range<size_t> &range);
3704  inline void join(GenEdgeData& rhs);
3705 
3706  inline TreeType& tree() { return mTree; }
3707 
3708 private:
3709  void operator=(const GenEdgeData&) {}
3710 
3711  struct Primitive { Vec3d a, b, c, d; Int32 index; };
3712 
3713  template<bool IsQuad>
3714  inline void voxelize(const Primitive&);
3715 
3716  template<bool IsQuad>
3717  inline bool evalPrimitive(const Coord&, const Primitive&);
3718 
3719  inline bool rayTriangleIntersection( const Vec3d& origin, const Vec3d& dir,
3720  const Vec3d& a, const Vec3d& b, const Vec3d& c, double& t);
3721 
3722 
3723  TreeType mTree;
3724  Accessor mAccessor;
3725 
3726  const std::vector<Vec3s>& mPointList;
3727  const std::vector<Vec4I>& mPolygonList;
3728 
3729  // Used internally for acceleration
3730  using IntTreeT = TreeType::ValueConverter<Int32>::Type;
3731  IntTreeT mLastPrimTree;
3732  tree::ValueAccessor<IntTreeT> mLastPrimAccessor;
3733 }; // class MeshToVoxelEdgeData::GenEdgeData
3734 
3735 
3736 inline
3737 MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3738  const std::vector<Vec3s>& pointList,
3739  const std::vector<Vec4I>& polygonList)
3740  : mTree(EdgeData())
3741  , mAccessor(mTree)
3742  , mPointList(pointList)
3743  , mPolygonList(polygonList)
3744  , mLastPrimTree(Int32(util::INVALID_IDX))
3745  , mLastPrimAccessor(mLastPrimTree)
3746 {
3747 }
3748 
3749 
3750 inline
3752  : mTree(EdgeData())
3753  , mAccessor(mTree)
3754  , mPointList(rhs.mPointList)
3755  , mPolygonList(rhs.mPolygonList)
3756  , mLastPrimTree(Int32(util::INVALID_IDX))
3757  , mLastPrimAccessor(mLastPrimTree)
3758 {
3759 }
3760 
3761 
3762 inline void
3764 {
3765  if (threaded) {
3766  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
3767  } else {
3768  (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
3769  }
3770 }
3771 
3772 
3773 inline void
3775 {
3776  using RootNodeType = TreeType::RootNodeType;
3777  using NodeChainType = RootNodeType::NodeChainType;
3778  static_assert(boost::mpl::size<NodeChainType>::value > 1, "expected tree height > 1");
3779  using InternalNodeType = boost::mpl::at<NodeChainType, boost::mpl::int_<1> >::type;
3780 
3781  Coord ijk;
3782  Index offset;
3783 
3784  rhs.mTree.clearAllAccessors();
3785 
3786  TreeType::LeafIter leafIt = rhs.mTree.beginLeaf();
3787  for ( ; leafIt; ++leafIt) {
3788  ijk = leafIt->origin();
3789 
3790  TreeType::LeafNodeType* lhsLeafPt = mTree.probeLeaf(ijk);
3791 
3792  if (!lhsLeafPt) {
3793 
3794  mAccessor.addLeaf(rhs.mAccessor.probeLeaf(ijk));
3795  InternalNodeType* node = rhs.mAccessor.getNode<InternalNodeType>();
3796  node->stealNode<TreeType::LeafNodeType>(ijk, EdgeData(), false);
3797  rhs.mAccessor.clear();
3798 
3799  } else {
3800 
3801  TreeType::LeafNodeType::ValueOnCIter it = leafIt->cbeginValueOn();
3802  for ( ; it; ++it) {
3803 
3804  offset = it.pos();
3805  const EdgeData& rhsValue = it.getValue();
3806 
3807  if (!lhsLeafPt->isValueOn(offset)) {
3808  lhsLeafPt->setValueOn(offset, rhsValue);
3809  } else {
3810 
3811  EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
3812 
3813  if (rhsValue.mXDist < lhsValue.mXDist) {
3814  lhsValue.mXDist = rhsValue.mXDist;
3815  lhsValue.mXPrim = rhsValue.mXPrim;
3816  }
3817 
3818  if (rhsValue.mYDist < lhsValue.mYDist) {
3819  lhsValue.mYDist = rhsValue.mYDist;
3820  lhsValue.mYPrim = rhsValue.mYPrim;
3821  }
3822 
3823  if (rhsValue.mZDist < lhsValue.mZDist) {
3824  lhsValue.mZDist = rhsValue.mZDist;
3825  lhsValue.mZPrim = rhsValue.mZPrim;
3826  }
3827 
3828  }
3829  } // end value iteration
3830  }
3831  } // end leaf iteration
3832 }
3833 
3834 
3835 inline void
3836 MeshToVoxelEdgeData::GenEdgeData::operator()(const tbb::blocked_range<size_t> &range)
3837 {
3838  Primitive prim;
3839 
3840  for (size_t n = range.begin(); n < range.end(); ++n) {
3841 
3842  const Vec4I& verts = mPolygonList[n];
3843 
3844  prim.index = Int32(n);
3845  prim.a = Vec3d(mPointList[verts[0]]);
3846  prim.b = Vec3d(mPointList[verts[1]]);
3847  prim.c = Vec3d(mPointList[verts[2]]);
3848 
3849  if (util::INVALID_IDX != verts[3]) {
3850  prim.d = Vec3d(mPointList[verts[3]]);
3851  voxelize<true>(prim);
3852  } else {
3853  voxelize<false>(prim);
3854  }
3855  }
3856 }
3857 
3858 
3859 template<bool IsQuad>
3860 inline void
3861 MeshToVoxelEdgeData::GenEdgeData::voxelize(const Primitive& prim)
3862 {
3863  std::deque<Coord> coordList;
3864  Coord ijk, nijk;
3865 
3866  ijk = Coord::floor(prim.a);
3867  coordList.push_back(ijk);
3868 
3869  evalPrimitive<IsQuad>(ijk, prim);
3870 
3871  while (!coordList.empty()) {
3872 
3873  ijk = coordList.back();
3874  coordList.pop_back();
3875 
3876  for (Int32 i = 0; i < 26; ++i) {
3877  nijk = ijk + util::COORD_OFFSETS[i];
3878 
3879  if (prim.index != mLastPrimAccessor.getValue(nijk)) {
3880  mLastPrimAccessor.setValue(nijk, prim.index);
3881  if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
3882  }
3883  }
3884  }
3885 }
3886 
3887 
3888 template<bool IsQuad>
3889 inline bool
3890 MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(const Coord& ijk, const Primitive& prim)
3891 {
3892  Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
3893  bool intersecting = false;
3894  double t;
3895 
3896  EdgeData edgeData;
3897  mAccessor.probeValue(ijk, edgeData);
3898 
3899  // Evaluate first triangle
3900  double dist = (org -
3901  closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, org, uvw)).lengthSqr();
3902 
3903  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
3904  if (t < edgeData.mXDist) {
3905  edgeData.mXDist = float(t);
3906  edgeData.mXPrim = prim.index;
3907  intersecting = true;
3908  }
3909  }
3910 
3911  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
3912  if (t < edgeData.mYDist) {
3913  edgeData.mYDist = float(t);
3914  edgeData.mYPrim = prim.index;
3915  intersecting = true;
3916  }
3917  }
3918 
3919  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
3920  if (t < edgeData.mZDist) {
3921  edgeData.mZDist = float(t);
3922  edgeData.mZPrim = prim.index;
3923  intersecting = true;
3924  }
3925  }
3926 
3927  if (IsQuad) {
3928  // Split quad into a second triangle and calculate distance.
3929  double secondDist = (org -
3930  closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, org, uvw)).lengthSqr();
3931 
3932  if (secondDist < dist) dist = secondDist;
3933 
3934  if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
3935  if (t < edgeData.mXDist) {
3936  edgeData.mXDist = float(t);
3937  edgeData.mXPrim = prim.index;
3938  intersecting = true;
3939  }
3940  }
3941 
3942  if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
3943  if (t < edgeData.mYDist) {
3944  edgeData.mYDist = float(t);
3945  edgeData.mYPrim = prim.index;
3946  intersecting = true;
3947  }
3948  }
3949 
3950  if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
3951  if (t < edgeData.mZDist) {
3952  edgeData.mZDist = float(t);
3953  edgeData.mZPrim = prim.index;
3954  intersecting = true;
3955  }
3956  }
3957  }
3958 
3959  if (intersecting) mAccessor.setValue(ijk, edgeData);
3960 
3961  return (dist < 0.86602540378443861);
3962 }
3963 
3964 
3965 inline bool
3966 MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
3967  const Vec3d& origin, const Vec3d& dir,
3968  const Vec3d& a, const Vec3d& b, const Vec3d& c,
3969  double& t)
3970 {
3971  // Check if ray is parallel with triangle
3972 
3973  Vec3d e1 = b - a;
3974  Vec3d e2 = c - a;
3975  Vec3d s1 = dir.cross(e2);
3976 
3977  double divisor = s1.dot(e1);
3978  if (!(std::abs(divisor) > 0.0)) return false;
3979 
3980  // Compute barycentric coordinates
3981 
3982  double inv_divisor = 1.0 / divisor;
3983  Vec3d d = origin - a;
3984  double b1 = d.dot(s1) * inv_divisor;
3985 
3986  if (b1 < 0.0 || b1 > 1.0) return false;
3987 
3988  Vec3d s2 = d.cross(e1);
3989  double b2 = dir.dot(s2) * inv_divisor;
3990 
3991  if (b2 < 0.0 || (b1 + b2) > 1.0) return false;
3992 
3993  // Compute distance to intersection point
3994 
3995  t = e2.dot(s2) * inv_divisor;
3996  return (t < 0.0) ? false : true;
3997 }
3998 
3999 
4001 
4002 
4003 inline
4005  : mTree(EdgeData())
4006 {
4007 }
4008 
4009 
4010 inline void
4012  const std::vector<Vec3s>& pointList,
4013  const std::vector<Vec4I>& polygonList)
4014 {
4015  GenEdgeData converter(pointList, polygonList);
4016  converter.run();
4017 
4018  mTree.clear();
4019  mTree.merge(converter.tree());
4020 }
4021 
4022 
4023 inline void
4025  Accessor& acc,
4026  const Coord& ijk,
4027  std::vector<Vec3d>& points,
4028  std::vector<Index32>& primitives)
4029 {
4030  EdgeData data;
4031  Vec3d point;
4032 
4033  Coord coord = ijk;
4034 
4035  if (acc.probeValue(coord, data)) {
4036 
4037  if (data.mXPrim != util::INVALID_IDX) {
4038  point[0] = double(coord[0]) + data.mXDist;
4039  point[1] = double(coord[1]);
4040  point[2] = double(coord[2]);
4041 
4042  points.push_back(point);
4043  primitives.push_back(data.mXPrim);
4044  }
4045 
4046  if (data.mYPrim != util::INVALID_IDX) {
4047  point[0] = double(coord[0]);
4048  point[1] = double(coord[1]) + data.mYDist;
4049  point[2] = double(coord[2]);
4050 
4051  points.push_back(point);
4052  primitives.push_back(data.mYPrim);
4053  }
4054 
4055  if (data.mZPrim != util::INVALID_IDX) {
4056  point[0] = double(coord[0]);
4057  point[1] = double(coord[1]);
4058  point[2] = double(coord[2]) + data.mZDist;
4059 
4060  points.push_back(point);
4061  primitives.push_back(data.mZPrim);
4062  }
4063 
4064  }
4065 
4066  coord[0] += 1;
4067 
4068  if (acc.probeValue(coord, data)) {
4069 
4070  if (data.mYPrim != util::INVALID_IDX) {
4071  point[0] = double(coord[0]);
4072  point[1] = double(coord[1]) + data.mYDist;
4073  point[2] = double(coord[2]);
4074 
4075  points.push_back(point);
4076  primitives.push_back(data.mYPrim);
4077  }
4078 
4079  if (data.mZPrim != util::INVALID_IDX) {
4080  point[0] = double(coord[0]);
4081  point[1] = double(coord[1]);
4082  point[2] = double(coord[2]) + data.mZDist;
4083 
4084  points.push_back(point);
4085  primitives.push_back(data.mZPrim);
4086  }
4087  }
4088 
4089  coord[2] += 1;
4090 
4091  if (acc.probeValue(coord, data)) {
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 
4102  coord[0] -= 1;
4103 
4104  if (acc.probeValue(coord, data)) {
4105 
4106  if (data.mXPrim != util::INVALID_IDX) {
4107  point[0] = double(coord[0]) + data.mXDist;
4108  point[1] = double(coord[1]);
4109  point[2] = double(coord[2]);
4110 
4111  points.push_back(point);
4112  primitives.push_back(data.mXPrim);
4113  }
4114 
4115  if (data.mYPrim != util::INVALID_IDX) {
4116  point[0] = double(coord[0]);
4117  point[1] = double(coord[1]) + data.mYDist;
4118  point[2] = double(coord[2]);
4119 
4120  points.push_back(point);
4121  primitives.push_back(data.mYPrim);
4122  }
4123  }
4124 
4125 
4126  coord[1] += 1;
4127 
4128  if (acc.probeValue(coord, data)) {
4129 
4130  if (data.mXPrim != util::INVALID_IDX) {
4131  point[0] = double(coord[0]) + data.mXDist;
4132  point[1] = double(coord[1]);
4133  point[2] = double(coord[2]);
4134 
4135  points.push_back(point);
4136  primitives.push_back(data.mXPrim);
4137  }
4138  }
4139 
4140  coord[2] -= 1;
4141 
4142  if (acc.probeValue(coord, data)) {
4143 
4144  if (data.mXPrim != util::INVALID_IDX) {
4145  point[0] = double(coord[0]) + data.mXDist;
4146  point[1] = double(coord[1]);
4147  point[2] = double(coord[2]);
4148 
4149  points.push_back(point);
4150  primitives.push_back(data.mXPrim);
4151  }
4152 
4153  if (data.mZPrim != util::INVALID_IDX) {
4154  point[0] = double(coord[0]);
4155  point[1] = double(coord[1]);
4156  point[2] = double(coord[2]) + data.mZDist;
4157 
4158  points.push_back(point);
4159  primitives.push_back(data.mZPrim);
4160  }
4161  }
4162 
4163  coord[0] += 1;
4164 
4165  if (acc.probeValue(coord, data)) {
4166 
4167  if (data.mZPrim != util::INVALID_IDX) {
4168  point[0] = double(coord[0]);
4169  point[1] = double(coord[1]);
4170  point[2] = double(coord[2]) + data.mZDist;
4171 
4172  points.push_back(point);
4173  primitives.push_back(data.mZPrim);
4174  }
4175  }
4176 }
4177 
4178 
4179 template<typename GridType, typename VecType>
4180 inline typename GridType::Ptr
4182  const openvdb::math::Transform& xform,
4183  typename VecType::ValueType halfWidth)
4184 {
4185  const Vec3s pmin = Vec3s(xform.worldToIndex(bbox.min()));
4186  const Vec3s pmax = Vec3s(xform.worldToIndex(bbox.max()));
4187 
4188  Vec3s points[8];
4189  points[0] = Vec3s(pmin[0], pmin[1], pmin[2]);
4190  points[1] = Vec3s(pmin[0], pmin[1], pmax[2]);
4191  points[2] = Vec3s(pmax[0], pmin[1], pmax[2]);
4192  points[3] = Vec3s(pmax[0], pmin[1], pmin[2]);
4193  points[4] = Vec3s(pmin[0], pmax[1], pmin[2]);
4194  points[5] = Vec3s(pmin[0], pmax[1], pmax[2]);
4195  points[6] = Vec3s(pmax[0], pmax[1], pmax[2]);
4196  points[7] = Vec3s(pmax[0], pmax[1], pmin[2]);
4197 
4198  Vec4I faces[6];
4199  faces[0] = Vec4I(0, 1, 2, 3); // bottom
4200  faces[1] = Vec4I(7, 6, 5, 4); // top
4201  faces[2] = Vec4I(4, 5, 1, 0); // front
4202  faces[3] = Vec4I(6, 7, 3, 2); // back
4203  faces[4] = Vec4I(0, 3, 7, 4); // left
4204  faces[5] = Vec4I(1, 5, 6, 2); // right
4205 
4206  QuadAndTriangleDataAdapter<Vec3s, Vec4I> mesh(points, 8, faces, 6);
4207 
4208  return meshToVolume<GridType>(mesh, xform, halfWidth, halfWidth);
4209 }
4210 
4211 
4212 } // namespace tools
4213 } // namespace OPENVDB_VERSION_NAME
4214 } // namespace openvdb
4215 
4216 #endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:755
bool *const mNodeMask
Definition: MeshToVolume.h:1332
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:390
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:237
SeedPoints(ConnectivityTable &connectivity, bool *changedNodeMask, bool *nodeMask, bool *changedVoxelMask)
Definition: MeshToVolume.h:1184
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:4181
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2928
float mZDist
Definition: MeshToVolume.h:465
std::vector< LeafNodeType * > & nodes()
Definition: MeshToVolume.h:788
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2261
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:281
typename tree::ValueAccessor< TreeType > AccessorType
Definition: MeshToVolume.h:2204
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1141
OPENVDB_API const Index32 INVALID_IDX
static ValueType epsilon()
Definition: MeshToVolume.h:551
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2892
static bool check(const ValueType v)
Definition: MeshToVolume.h:1737
UnionValueMasks(std::vector< LeafNodeTypeA * > &nodesA, std::vector< LeafNodeTypeB * > &nodesB)
Definition: MeshToVolume.h:2238
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:460
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:668
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
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition: Vec3.h:218
bool checkNeighbours(const Coord &ijk, AccessorType &acc, bool(&mask)[26])
Definition: MeshToVolume.h:1666
tbb::enumerable_thread_specific< typename VoxelizationDataType::Ptr > DataTable
Definition: MeshToVolume.h:1971
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
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition: Tree.h:1458
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1087
Base class for tree-traversal iterators over tile and voxel values.
Definition: TreeIterator.h:621
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1193
math::Transform const *const mXform
Definition: MeshToVolume.h:544
void expand(ElementType padding)
Pad this bounding box.
Definition: BBox.h:321
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1692
Index32 mXPrim
Definition: MeshToVolume.h:466
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2980
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1682
Vec3< double > Vec3d
Definition: Vec3.h:662
CombineLeafNodes(TreeType &lhsDistTree, Int32TreeType &lhsIdxTree, LeafNodeType **rhsDistNodes, Int32LeafNodeType **rhsIdxNodes)
Definition: MeshToVolume.h:569
Contiguous quad and triangle data adapter class.
Definition: MeshToVolume.h:165
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:195
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form &#39;someVar << "text" << ...&#39;.
Definition: logging.h:263
SweepExteriorSign(Axis axis, const std::vector< size_t > &startNodeIndices, ConnectivityTable &connectivity)
Definition: MeshToVolume.h:818
void join(GenEdgeData &rhs)
Definition: MeshToVolume.h:3774
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:266
std::vector< LeafNodeType * > & updatedDistNodes()
Definition: MeshToVolume.h:2583
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:512
EdgeData(float dist=1.0)
Definition: MeshToVolume.h:442
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1114
uint32_t Index32
Definition: Types.h:29
Index32 mYPrim
Definition: MeshToVolume.h:466
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
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:1152
std::unique_ptr< bool[]> MaskArray
Definition: MeshToVolume.h:1348
const size_t * offsetsPrevX() const
Definition: MeshToVolume.h:793
Int32TreeType indexTree
Definition: MeshToVolume.h:1925
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1783
std::shared_ptr< T > SharedPtr
Definition: Types.h:91
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3683
const size_t * offsetsPrevY() const
Definition: MeshToVolume.h:796
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:638
Efficient multi-threaded replacement of the background values in tree.
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2207
PointType *const mPointsOut
Definition: MeshToVolume.h:543
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:230
std::vector< Int32LeafNodeType * > & updatedIndexNodes()
Definition: MeshToVolume.h:2586
typename TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1344
Definition: Math.h:851
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:522
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2981
Definition: Coord.h:587
ReleaseChildNodes(NodeType **nodes)
Definition: MeshToVolume.h:1795
void fillArray(ValueType *array, const ValueType val, const size_t length)
Definition: MeshToVolume.h:1128
Axis-aligned bounding box.
Definition: BBox.h:23
ConstructVoxelMask(BoolTreeType &maskTree, const TreeType &tree, std::vector< LeafNodeType * > &nodes)
Definition: MeshToVolume.h:2264
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition: BBox.h:62
Type Pow2(Type x)
Return x2.
Definition: Math.h:495
TBB body object to voxelize a mesh of triangles and/or quads into a collection of VDB grids...
Definition: MeshToVolume.h:1898
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2840
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1681
bool *const mChangedVoxelMask
Definition: MeshToVolume.h:1172
PointType const *const mPointsIn
Definition: MeshToVolume.h:542
Ptr copy() const
Definition: Transform.h:50
tbb::enumerable_thread_specific< LocalData > LocalDataTable
Definition: MeshToVolume.h:1350
void operator()() const
Definition: MeshToVolume.h:1839
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2372
std::vector< LeafNodeType * > *const mNodes
Definition: MeshToVolume.h:2759
Int32TreeType *const mIndexTree
Definition: MeshToVolume.h:1785
InactivateValues(std::vector< LeafNodeType * > &nodes, ValueType exBandWidth, ValueType inBandWidth)
Definition: MeshToVolume.h:2842
typename TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:2375
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition: Tree.h:1155
TransformValues(std::vector< LeafNodeType * > &nodes, ValueType voxelSize, bool unsignedDist)
Definition: MeshToVolume.h:2804
OPENVDB_API const Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
LeafNodeType **const mNodes
Definition: MeshToVolume.h:653
TreeType & tree()
Definition: MeshToVolume.h:3706
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:3655
void join(ConstructVoxelMask &rhs)
Definition: MeshToVolume.h:2357
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1901
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1733
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:3096
LeafNodeType **const mNodes
Definition: MeshToVolume.h:675
Renormalize(const TreeType &tree, const std::vector< LeafNodeType * > &nodes, ValueType *buffer, ValueType voxelSize)
Definition: MeshToVolume.h:2919
Defined various multi-threaded utility functions for trees.
ComputeIntersectingVoxelSign(std::vector< LeafNodeType * > &distNodes, const TreeType &distTree, const Int32TreeType &indexTree, const MeshDataAdapter &mesh)
Definition: MeshToVolume.h:1352
bool operator<(const Fragment &rhs) const
Definition: MeshToVolume.h:2392
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:815
ExpandNarrowband(const ExpandNarrowband &rhs, tbb::split)
Definition: MeshToVolume.h:2422
std::unique_ptr< VoxelizationData > Ptr
Definition: MeshToVolume.h:1900
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:4011
typename Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:1345
void operator()() const
Definition: MeshToVolume.h:2750
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2916
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:661
Definition: Math.h:850
bool operator>(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:192
Internal edge data type.
Definition: MeshToVolume.h:441
std::vector< LeafNodeType * > & newDistNodes()
Definition: MeshToVolume.h:2582
void seedFill(LeafNodeType &node)
Definition: MeshToVolume.h:938
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:3624
Tree< typename RootNodeType::template ValueConverter< Int32 >::Type > Type
Definition: Tree.h:197
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2839
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1142
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:814
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1104
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1181
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition: TreeIterator.h:697
void clear() override
Remove all nodes from this cache, then reinsert the root node.
Definition: ValueAccessor.h:404
DiffLeafNodeMask(const TreeType &rhsTree, std::vector< BoolLeafNodeType * > &lhsNodes)
Definition: MeshToVolume.h:2210
std::vector< Int32LeafNodeType * > & newIndexNodes()
Definition: MeshToVolume.h:2585
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
float mXDist
Definition: MeshToVolume.h:465
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:189
static ValueType minNarrowBandWidth()
Definition: MeshToVolume.h:552
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2373
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1982
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Accessor getAccessor()
Definition: MeshToVolume.h:496
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:261
void clearAllAccessors()
Clear all registered accessors.
Definition: Tree.h:1517
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1342
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:3593
int32_t Int32
Definition: Types.h:33
typename TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1903
Fragment(Int32 idx_, Int32 x_, Int32 y_, Int32 z_, ValueType dist_)
Definition: MeshToVolume.h:2387
size_t pointCount() const
Definition: MeshToVolume.h:186
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2850
bool *const mChangedNodeMask
Definition: MeshToVolume.h:1331
_RootNodeType RootNodeType
Definition: Tree.h:182
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: MeshToVolume.h:2377
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:578
RestoreOrigin(std::vector< LeafNodeType * > &nodes, const Coord *coordinates)
Definition: MeshToVolume.h:663
MeshToVoxelEdgeData()
Definition: MeshToVolume.h:4004
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Extracts and stores voxel edge intersection data from a mesh.
Definition: MeshToVolume.h:434
MinCombine(std::vector< LeafNodeType * > &nodes, const ValueType *buffer)
Definition: MeshToVolume.h:2983
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2917
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1366
std::ostream & operator<<(std::ostream &ostr, const MeshToVoxelEdgeData::EdgeData &rhs)
Definition: MeshToVolume.h:3673
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2208
UCharTreeType primIdTree
Definition: MeshToVolume.h:1928
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1734
size_t vertexCount(size_t n) const
Vertex count for polygon n.
Definition: MeshToVolume.h:189
bool processY(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1254
StealUniqueLeafNodes(TreeType &lhsTree, TreeType &rhsTree, std::vector< LeafNodeType * > &overlappingNodes)
Definition: MeshToVolume.h:1831
Definition: Mat4.h:24
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1726
NodeType **const mNodes
Definition: MeshToVolume.h:1806
const ValueType mValue
Definition: MeshToVolume.h:1122
Definition: Exceptions.h:13
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition: BBox.h:64
bool *const mChangedVoxelMask
Definition: MeshToVolume.h:1333
bool processX(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1292
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
bool operator==(const EdgeData &rhs) const
Definition: MeshToVolume.h:460
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:395
void maskNodeInternalNeighbours(const Index pos, bool(&mask)[26])
Definition: MeshToVolume.h:1535
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2205
TreeType *const mTree
Definition: MeshToVolume.h:2758
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...
VoxelizePolygons(DataTable &dataTable, const MeshDataAdapter &mesh, Interrupter *interrupter=nullptr)
Definition: MeshToVolume.h:1973
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:645
size_t polygonCount() const
Definition: MeshToVolume.h:185
bool scanFill(LeafNodeType &node)
Definition: MeshToVolume.h:1014
bool traceVoxelLine(LeafNodeType &node, Index pos, Index step) const
Definition: MeshToVolume.h:902
void run(bool threaded=true)
Definition: MeshToVolume.h:3763
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
LeafNodeType **const mNodes
Definition: MeshToVolume.h:1170
AddNodes(TreeType &tree, std::vector< LeafNodeType * > &nodes)
Definition: MeshToVolume.h:2745
SyncVoxelMask(std::vector< LeafNodeType * > &nodes, const bool *changedNodeMask, bool *changedVoxelMask)
Definition: MeshToVolume.h:1144
size_t findNeighbourNode(tree::ValueAccessor< const TreeType > &acc, const Coord &start, const Coord &step) const
Definition: MeshToVolume.h:724
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:708
EdgeData operator-(const T &) const
Definition: MeshToVolume.h:456
Definition: Types.h:453
std::pair< PointArray, MaskArray > LocalData
Definition: MeshToVolume.h:1349
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1747
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels) ...
Definition: TreeIterator.h:1191
typename TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:1735
UCharTreeAcc primIdAcc
Definition: MeshToVolume.h:1929
ExpandNarrowband(std::vector< BoolLeafNodeType * > &maskNodes, BoolTreeType &maskTree, TreeType &distTree, Int32TreeType &indexTree, const MeshDataAdapter &mesh, ValueType exteriorBandWidth, ValueType interiorBandWidth, ValueType voxelSize)
Definition: MeshToVolume.h:2397
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2802
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:826
Definition: Transform.h:39
bool processZ(const size_t n, bool firstFace) const
Definition: MeshToVolume.h:1216
unsigned char getNewPrimId()
Definition: MeshToVolume.h:1931
GenEdgeData(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition: MeshToVolume.h:3737
ValueType *const mArray
Definition: MeshToVolume.h:1121
MeshToVolumeFlags
Mesh to volume conversion flags.
Definition: MeshToVolume.h:61
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2259
Index32 mZPrim
Definition: MeshToVolume.h:466
typename RootNodeType::LeafNodeType LeafNodeType
Definition: Tree.h:185
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1829
Int32TreeAcc indexAcc
Definition: MeshToVolume.h:1926
Index32 Index
Definition: Types.h:31
ComputeNodeConnectivity(const TreeType &tree, const Coord *coordinates, size_t *offsets, size_t numNodes, const CoordBBox &bbox)
Definition: MeshToVolume.h:686
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
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:359
math::Vec4< Index32 > Vec4I
Definition: Types.h:65
Vec3< float > Vec3s
Definition: Vec3.h:661
ValidateIntersectingVoxels(TreeType &tree, std::vector< LeafNodeType * > &nodes)
Definition: MeshToVolume.h:1686
Definition: Mat.h:170
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition: Tuple.h:180
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
QuadAndTriangleDataAdapter(const std::vector< PointType > &points, const std::vector< PolygonType > &polygons)
Definition: MeshToVolume.h:167
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2884
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2801
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2244
EdgeData operator+(const T &) const
Definition: MeshToVolume.h:455
StashOriginAndStoreOffset(std::vector< LeafNodeType * > &nodes, Coord *coordinates)
Definition: MeshToVolume.h:640
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
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1086
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:102
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:701
typename TreeType::template ValueConverter< Int32 >::Type Int32TreeType
Definition: MeshToVolume.h:564
float mYDist
Definition: MeshToVolume.h:465
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:4024
void join(ExpandNarrowband &rhs)
Definition: MeshToVolume.h:2439
RemoveSelfIntersectingSurface(std::vector< LeafNodeType * > &nodes, TreeType &distTree, Int32TreeType &indexTree)
Definition: MeshToVolume.h:1739
Definition: Tree.h:176
const std::vector< LeafNodeType * > & nodes() const
Definition: MeshToVolume.h:789
bool const *const mChangedNodeMask
Definition: MeshToVolume.h:1171
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1797
const size_t * offsetsNextX() const
Definition: MeshToVolume.h:792
Coord const *const mCoordinates
Definition: MeshToVolume.h:676
typename TreeType::template ValueConverter< unsigned char >::Type UCharTreeType
Definition: MeshToVolume.h:1904
void releaseLeafNodes(TreeType &tree)
Definition: MeshToVolume.h:1812
QuadAndTriangleDataAdapter(const PointType *pointArray, size_t pointArraySize, const PolygonType *polygonArray, size_t polygonArraySize)
Definition: MeshToVolume.h:176
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2281
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the voxel as well as its value.
Definition: ValueAccessor.h:240
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:388
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition: Transform.h:93
SeedFillExteriorSign(std::vector< LeafNodeType * > &nodes, bool *changedNodeMask)
Definition: MeshToVolume.h:1089
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:340
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2216
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:270
const size_t * offsetsNextY() const
Definition: MeshToVolume.h:795
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:2453
ConstructVoxelMask(ConstructVoxelMask &rhs, tbb::split)
Definition: MeshToVolume.h:2273
BoolTreeType & newMaskTree()
Definition: MeshToVolume.h:2580
EdgeData operator-() const
Definition: MeshToVolume.h:457
typename Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:567
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2262
LeafNodeConnectivityTable(TreeType &tree)
Definition: MeshToVolume.h:757
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:1180
static bool check(const ValueType v)
Definition: MeshToVolume.h:1684
FloatTreeAcc distAcc
Definition: MeshToVolume.h:1923
OffsetValues(std::vector< LeafNodeType * > &nodes, ValueType offset)
Definition: MeshToVolume.h:2887
const size_t * offsetsPrevZ() const
Definition: MeshToVolume.h:799
FillArray(ValueType *array, const ValueType v)
Definition: MeshToVolume.h:1112
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:1704
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2812
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition: Tree.h:1832
typename LeafNodeType::NodeMaskType NodeMaskType
Definition: MeshToVolume.h:2374
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
bool *const mChangedNodeMask
Definition: MeshToVolume.h:1105
Definition: Math.h:852
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:503
const size_t * offsetsNextZ() const
Definition: MeshToVolume.h:798
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:1095
typename BoolTreeType::LeafNodeType BoolLeafNodeType
Definition: MeshToVolume.h:2378
Definition: Exceptions.h:64
typename Int32TreeType::LeafNodeType Int32LeafNodeType
Definition: MeshToVolume.h:2376
Definition: Types.h:454
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:1343
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:684
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: MeshToVolume.h:2988
ValueConverter<T>::Type is the type of a tree having the same hierarchy as this tree but a different ...
Definition: Tree.h:196
ConnectivityTable *const mConnectivity
Definition: MeshToVolume.h:1330
NodeType * getNode()
Return the cached node of type NodeType. [Mainly for internal use].
Definition: ValueAccessor.h:315
Partitions points into BucketLog2Dim aligned buckets using a parallel radix-based sorting algorithm...
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:566
typename TreeType::ValueType ValueType
Definition: MeshToVolume.h:2885
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106
TransformPoints(const PointType *pointsIn, PointType *pointsOut, const math::Transform &xform)
Definition: MeshToVolume.h:516
typename TreeType::LeafNodeType LeafNodeType
Definition: MeshToVolume.h:2743
void operator()(const tbb::blocked_range< size_t > &range)
Definition: MeshToVolume.h:3836