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