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