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