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