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