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