OpenVDB  12.0.0
VolumeToMesh.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 
4 /// @file VolumeToMesh.h
5 ///
6 /// @brief Extract polygonal surfaces from scalar volumes.
7 ///
8 /// @author Mihai Alden
9 
10 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
12 
13 #include <openvdb/Platform.h>
14 #include <openvdb/math/Operators.h> // for ISGradient
16 #include <openvdb/util/Util.h> // for INVALID_IDX
17 #include <openvdb/util/Assert.h>
18 #include <openvdb/openvdb.h>
19 
20 #include <tbb/blocked_range.h>
21 #include <tbb/parallel_for.h>
22 #include <tbb/parallel_reduce.h>
23 #include <tbb/task_arena.h>
24 
25 #include <cmath> // for std::isfinite()
26 #include <cstring> // for std::memset
27 #include <map>
28 #include <memory>
29 #include <set>
30 #include <type_traits>
31 #include <vector>
32 
33 
34 namespace openvdb {
36 namespace OPENVDB_VERSION_NAME {
37 namespace tools {
38 
39 
40 ////////////////////////////////////////
41 
42 
43 // Wrapper functions for the VolumeToMesh converter
44 
45 
46 /// @brief Uniformly mesh any scalar grid that has a continuous isosurface.
47 ///
48 /// @param grid a scalar grid to mesh
49 /// @param points output list of world space points
50 /// @param quads output quad index list
51 /// @param isovalue determines which isosurface to mesh
52 ///
53 /// @throw TypeError if @a grid does not have a scalar value type
54 template<typename GridType>
55 void
57  const GridType& grid,
58  std::vector<Vec3s>& points,
59  std::vector<Vec4I>& quads,
60  double isovalue = 0.0);
61 
62 
63 /// @brief Adaptively mesh any scalar grid that has a continuous isosurface.
64 /// @details When converting to polygons, the adaptivity threshold determines
65 /// how closely the isosurface is matched by the resulting mesh. Higher
66 /// thresholds will allow more variation in polygon size, using fewer
67 /// polygons to express the surface. Triangles will only be created for
68 /// areas of the mesh which hit the adaptivity threshold and can't be
69 /// represented as quads.
70 /// @note Do not use this method just to get a triangle mesh - use the above
71 /// method and post process the quad index list.
72 ///
73 /// @param grid a scalar grid to mesh
74 /// @param points output list of world space points
75 /// @param triangles output triangle index list
76 /// @param quads output quad index list
77 /// @param isovalue determines which isosurface to mesh
78 /// @param adaptivity surface adaptivity threshold [0 to 1]
79 /// @param relaxDisorientedTriangles toggle relaxing disoriented triangles during
80 /// adaptive meshing.
81 ///
82 /// @throw TypeError if @a grid does not have a scalar value type
83 template<typename GridType>
84 void
86  const GridType& grid,
87  std::vector<Vec3s>& points,
88  std::vector<Vec3I>& triangles,
89  std::vector<Vec4I>& quads,
90  double isovalue = 0.0,
91  double adaptivity = 0.0,
92  bool relaxDisorientedTriangles = true);
93 
94 
95 ////////////////////////////////////////
96 
97 
98 /// @brief Polygon flags, used for reference based meshing.
100 
101 
102 /// @brief Collection of quads and triangles
104 {
105 public:
106 
107  inline PolygonPool();
108  inline PolygonPool(const size_t numQuads, const size_t numTriangles);
109 
110  inline void copy(const PolygonPool& rhs);
111 
112  inline void resetQuads(size_t size);
113  inline void clearQuads();
114 
115  inline void resetTriangles(size_t size);
116  inline void clearTriangles();
117 
118 
119  // polygon accessor methods
120 
121  const size_t& numQuads() const { return mNumQuads; }
122 
123  openvdb::Vec4I& quad(size_t n) { return mQuads[n]; }
124  const openvdb::Vec4I& quad(size_t n) const { return mQuads[n]; }
125 
126 
127  const size_t& numTriangles() const { return mNumTriangles; }
128 
129  openvdb::Vec3I& triangle(size_t n) { return mTriangles[n]; }
130  const openvdb::Vec3I& triangle(size_t n) const { return mTriangles[n]; }
131 
132 
133  // polygon flags accessor methods
134 
135  char& quadFlags(size_t n) { return mQuadFlags[n]; }
136  const char& quadFlags(size_t n) const { return mQuadFlags[n]; }
137 
138  char& triangleFlags(size_t n) { return mTriangleFlags[n]; }
139  const char& triangleFlags(size_t n) const { return mTriangleFlags[n]; }
140 
141 
142  // reduce the polygon containers, n has to
143  // be smaller than the current container size.
144 
145  inline bool trimQuads(const size_t n, bool reallocate = false);
146  inline bool trimTrinagles(const size_t n, bool reallocate = false);
147 
148 private:
149  // disallow copy by assignment
150  void operator=(const PolygonPool&) {}
151 
152  size_t mNumQuads, mNumTriangles;
153  std::unique_ptr<openvdb::Vec4I[]> mQuads;
154  std::unique_ptr<openvdb::Vec3I[]> mTriangles;
155  std::unique_ptr<char[]> mQuadFlags, mTriangleFlags;
156 };
157 
158 
159 /// @{
160 /// @brief Point and primitive list types.
161 using PointList = std::unique_ptr<openvdb::Vec3s[]>;
162 using PolygonPoolList = std::unique_ptr<PolygonPool[]>;
163 /// @}
164 
165 
166 ////////////////////////////////////////
167 
168 
169 /// @brief Mesh any scalar grid that has a continuous isosurface.
171 {
172 
173  /// @param isovalue Determines which isosurface to mesh.
174  /// @param adaptivity Adaptivity threshold [0 to 1]
175  /// @param relaxDisorientedTriangles Toggle relaxing disoriented triangles during
176  /// adaptive meshing.
177  VolumeToMesh(double isovalue = 0, double adaptivity = 0, bool relaxDisorientedTriangles = true);
178 
179  //////////
180 
181  /// @{
182  // Mesh data accessors
183 
184  size_t pointListSize() const { return mPointListSize; }
185  PointList& pointList() { return mPoints; }
186  const PointList& pointList() const { return mPoints; }
187 
188  size_t polygonPoolListSize() const { return mPolygonPoolListSize; }
189  PolygonPoolList& polygonPoolList() { return mPolygons; }
190  const PolygonPoolList& polygonPoolList() const { return mPolygons; }
191 
192  std::vector<uint8_t>& pointFlags() { return mPointFlags; }
193  const std::vector<uint8_t>& pointFlags() const { return mPointFlags; }
194  /// @}
195 
196 
197  //////////
198 
199 
200  /// @brief Main call
201  /// @note Call with scalar typed grid.
202  template<typename InputGridType>
203  void operator()(const InputGridType&);
204 
205 
206  //////////
207 
208 
209  /// @brief When surfacing fractured SDF fragments, the original unfractured
210  /// SDF grid can be used to eliminate seam lines and tag polygons that are
211  /// coincident with the reference surface with the @c POLYFLAG_EXTERIOR
212  /// flag and polygons that are in proximity to the seam lines with the
213  /// @c POLYFLAG_FRACTURE_SEAM flag. (The performance cost for using this
214  /// reference based scheme compared to the regular meshing scheme is
215  /// approximately 15% for the first fragment and neglect-able for
216  /// subsequent fragments.)
217  ///
218  /// @note Attributes from the original asset such as uv coordinates, normals etc.
219  /// are typically transferred to polygons that are marked with the
220  /// @c POLYFLAG_EXTERIOR flag. Polygons that are not marked with this flag
221  /// are interior to reference surface and might need projected UV coordinates
222  /// or a different material. Polygons marked as @c POLYFLAG_FRACTURE_SEAM can
223  /// be used to drive secondary elements such as debris and dust in a FX pipeline.
224  ///
225  /// @param grid reference surface grid of @c GridT type.
226  /// @param secAdaptivity Secondary adaptivity threshold [0 to 1]. Used in regions
227  /// that do not exist in the reference grid. (Parts of the
228  /// fragment surface that are not coincident with the
229  /// reference surface.)
230  void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0);
231 
232 
233  /// @param mask A boolean grid whose active topology defines the region to mesh.
234  /// @param invertMask Toggle to mesh the complement of the mask.
235  /// @note The mask's tree configuration has to match @c GridT's tree configuration.
236  void setSurfaceMask(const GridBase::ConstPtr& mask, bool invertMask = false);
237 
238  /// @param grid A scalar grid used as a spatial multiplier for the adaptivity threshold.
239  /// @note The grid's tree configuration has to match @c GridT's tree configuration.
240  void setSpatialAdaptivity(const GridBase::ConstPtr& grid);
241 
242 
243  /// @param tree A boolean tree whose active topology defines the adaptivity mask.
244  /// @note The tree configuration has to match @c GridT's tree configuration.
245  void setAdaptivityMask(const TreeBase::ConstPtr& tree);
246 
247 private:
248  // Disallow copying
249  VolumeToMesh(const VolumeToMesh&);
250  VolumeToMesh& operator=(const VolumeToMesh&);
251 
252 
253  PointList mPoints;
254  PolygonPoolList mPolygons;
255 
256  size_t mPointListSize, mSeamPointListSize, mPolygonPoolListSize;
257  double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
258 
259  GridBase::ConstPtr mRefGrid, mSurfaceMaskGrid, mAdaptivityGrid;
260  TreeBase::ConstPtr mAdaptivityMaskTree;
261 
262  TreeBase::Ptr mRefSignTree, mRefIdxTree;
263 
264  bool mInvertSurfaceMask, mRelaxDisorientedTriangles;
265 
266  std::unique_ptr<uint32_t[]> mQuantizedSeamPoints;
267  std::vector<uint8_t> mPointFlags;
268 }; // struct VolumeToMesh
269 
270 
271 ////////////////////////////////////////
272 
273 
274 /// @brief Given a set of tangent elements, @c points with corresponding @c normals,
275 /// this method returns the intersection point of all tangent elements.
276 ///
277 /// @note Used to extract surfaces with sharp edges and corners from volume data,
278 /// see the following paper for details: "Feature Sensitive Surface
279 /// Extraction from Volume Data, Kobbelt et al. 2001".
280 inline Vec3d findFeaturePoint(
281  const std::vector<Vec3d>& points,
282  const std::vector<Vec3d>& normals)
283 {
284  using Mat3d = math::Mat3d;
285 
286  Vec3d avgPos(0.0);
287 
288  if (points.empty()) return avgPos;
289 
290  for (size_t n = 0, N = points.size(); n < N; ++n) {
291  avgPos += points[n];
292  }
293 
294  avgPos /= double(points.size());
295 
296  // Unique components of the 3x3 A^TA matrix, where A is
297  // the matrix of normals.
298  double m00=0,m01=0,m02=0,
299  m11=0,m12=0,
300  m22=0;
301 
302  // The rhs vector, A^Tb, where b = n dot p
303  Vec3d rhs(0.0);
304 
305  for (size_t n = 0, N = points.size(); n < N; ++n) {
306 
307  const Vec3d& n_ref = normals[n];
308 
309  // A^TA
310  m00 += n_ref[0] * n_ref[0]; // diagonal
311  m11 += n_ref[1] * n_ref[1];
312  m22 += n_ref[2] * n_ref[2];
313 
314  m01 += n_ref[0] * n_ref[1]; // Upper-tri
315  m02 += n_ref[0] * n_ref[2];
316  m12 += n_ref[1] * n_ref[2];
317 
318  // A^Tb (centered around the origin)
319  rhs += n_ref * n_ref.dot(points[n] - avgPos);
320  }
321 
322  Mat3d A(m00,m01,m02,
323  m01,m11,m12,
324  m02,m12,m22);
325 
326  /*
327  // Inverse
328  const double det = A.det();
329  if (det > 0.01) {
330  Mat3d A_inv = A.adjoint();
331  A_inv *= (1.0 / det);
332 
333  return avgPos + A_inv * rhs;
334  }
335  */
336 
337  // Compute the pseudo inverse
338 
339  math::Mat3d eigenVectors;
340  Vec3d eigenValues;
341 
342  diagonalizeSymmetricMatrix(A, eigenVectors, eigenValues, 300);
343 
344  Mat3d D = Mat3d::identity();
345 
346 
347  double tolerance = std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
348  tolerance = std::max(tolerance, std::abs(eigenValues[2]));
349  tolerance *= 0.01;
350 
351  int clamped = 0;
352  for (int i = 0; i < 3; ++i ) {
353  if (std::abs(eigenValues[i]) < tolerance) {
354  D[i][i] = 0.0;
355  ++clamped;
356  } else {
357  D[i][i] = 1.0 / eigenValues[i];
358  }
359  }
360 
361  // Assemble the pseudo inverse and calc. the intersection point
362  if (clamped < 3) {
363  Mat3d pseudoInv = eigenVectors * D * eigenVectors.transpose();
364  return avgPos + pseudoInv * rhs;
365  }
366 
367  return avgPos;
368 }
369 
370 
371 ////////////////////////////////////////////////////////////////////////////////
372 ////////////////////////////////////////////////////////////////////////////////
373 
374 
375 // Internal utility objects and implementation details
376 
377 /// @cond OPENVDB_DOCS_INTERNAL
378 
379 namespace volume_to_mesh_internal {
380 
381 template<typename ValueType>
382 struct FillArray
383 {
384  FillArray(ValueType* array, const ValueType& v) : mArray(array), mValue(v) { }
385 
386  void operator()(const tbb::blocked_range<size_t>& range) const {
387  const ValueType v = mValue;
388  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
389  mArray[n] = v;
390  }
391  }
392 
393  ValueType * const mArray;
394  const ValueType mValue;
395 };
396 
397 
398 template<typename ValueType>
399 inline void
400 fillArray(ValueType* array, const ValueType& val, const size_t length)
401 {
402  const auto grainSize = std::max<size_t>(
403  length / tbb::this_task_arena::max_concurrency(), 1024);
404  const tbb::blocked_range<size_t> range(0, length, grainSize);
405  tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
406 }
407 
408 
409 /// @brief Bit-flags used to classify cells.
410 enum { SIGNS = 0xFF, EDGES = 0xE00, INSIDE = 0x100,
411  XEDGE = 0x200, YEDGE = 0x400, ZEDGE = 0x800, SEAM = 0x1000};
412 
413 
414 /// @brief Used to quickly determine if a given cell is adaptable.
415 const bool sAdaptable[256] = {
416  1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,
417  1,0,1,1,0,0,1,1,0,0,0,1,0,0,1,1,1,1,1,1,0,0,1,1,0,1,0,1,0,0,0,1,
418  1,0,0,0,1,0,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
419  1,0,1,1,1,0,1,1,0,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,0,0,0,0,0,0,0,1,
420  1,0,0,0,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,1,1,0,1,1,1,0,1,
421  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,0,1,0,0,0,1,
422  1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,1,1,1,0,0,1,0,0,0,1,1,0,0,1,1,0,1,
423  1,0,1,0,1,0,1,0,1,0,0,0,1,0,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1};
424 
425 
426 /// @brief Contains the ambiguous face index for certain cell configuration.
427 const unsigned char sAmbiguousFace[256] = {
428  0,0,0,0,0,5,0,0,0,0,5,0,0,0,0,0,0,0,1,0,0,5,1,0,4,0,0,0,4,0,0,0,
429  0,1,0,0,2,0,0,0,0,1,5,0,2,0,0,0,0,0,0,0,2,0,0,0,4,0,0,0,0,0,0,0,
430  0,0,2,2,0,5,0,0,3,3,0,0,0,0,0,0,6,6,0,0,6,0,0,0,0,0,0,0,0,0,0,0,
431  0,1,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
432  0,4,0,4,3,0,3,0,0,0,5,0,0,0,0,0,0,0,1,0,3,0,0,0,0,0,0,0,0,0,0,0,
433  6,0,6,0,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
434  0,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
435  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
436 
437 
438 /// @brief Lookup table for different cell sign configurations. The first entry specifies
439 /// the total number of points that need to be generated inside a cell and the
440 /// remaining 12 entries indicate different edge groups.
441 const unsigned char sEdgeGroupTable[256][13] = {
442  {0,0,0,0,0,0,0,0,0,0,0,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},
443  {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},{1,1,1,1,1,0,0,0,0,1,0,1,0},
444  {1,1,0,1,0,0,0,0,0,0,1,1,0},{1,0,0,1,1,0,0,0,0,1,1,1,0},{1,0,0,1,1,0,0,0,0,0,0,0,1},
445  {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,1,1,1,1,0,0,0,0,0,1,0,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},
446  {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},
447  {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,0,0,0,0,1,0,0,1,1,0,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},
448  {1,1,1,0,0,1,0,0,1,1,1,0,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},
449  {1,1,1,1,1,1,0,0,1,0,0,1,0},{1,1,0,1,0,1,0,0,1,1,1,1,0},{1,0,0,1,1,1,0,0,1,0,1,1,0},
450  {1,0,0,1,1,1,0,0,1,1,0,0,1},{1,1,0,1,0,1,0,0,1,0,0,0,1},{2,2,1,1,2,1,0,0,1,2,1,0,1},
451  {1,0,1,1,0,1,0,0,1,0,1,0,1},{1,0,1,0,1,1,0,0,1,1,0,1,1},{1,1,1,0,0,1,0,0,1,0,0,1,1},
452  {2,1,0,0,1,2,0,0,2,1,2,2,2},{1,0,0,0,0,1,0,0,1,0,1,1,1},{1,0,0,0,0,1,1,0,0,0,1,0,0},
453  {1,1,0,0,1,1,1,0,0,1,1,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},
454  {1,0,1,1,0,1,1,0,0,0,1,1,0},{2,2,2,1,1,1,1,0,0,1,2,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},
455  {1,0,0,1,1,1,1,0,0,1,0,1,0},{2,0,0,2,2,1,1,0,0,0,1,0,2},{1,1,0,1,0,1,1,0,0,1,1,0,1},
456  {1,1,1,1,1,1,1,0,0,0,0,0,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},{1,0,1,0,1,1,1,0,0,0,1,1,1},
457  {2,1,1,0,0,2,2,0,0,2,1,2,2},{1,1,0,0,1,1,1,0,0,0,0,1,1},{1,0,0,0,0,1,1,0,0,1,0,1,1},
458  {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},
459  {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,0,1,1,0,0,1,0,1,1,1,1,0},{2,1,1,2,2,0,2,0,2,0,1,2,0},
460  {1,1,0,1,0,0,1,0,1,1,0,1,0},{1,0,0,1,1,0,1,0,1,0,0,1,0},{1,0,0,1,1,0,1,0,1,1,1,0,1},
461  {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,1,2,2,1,0,2,0,2,1,0,0,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},
462  {2,0,2,0,2,0,1,0,1,2,2,1,1},{2,2,2,0,0,0,1,0,1,0,2,1,1},{2,2,0,0,2,0,1,0,1,2,0,1,1},
463  {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,0,0,0,0,0,1,1,0,0,0,1,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},
464  {1,1,1,0,0,0,1,1,0,0,1,1,0},{1,0,1,0,1,0,1,1,0,1,1,1,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},
465  {1,1,1,1,1,0,1,1,0,1,0,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},{1,0,0,1,1,0,1,1,0,1,1,0,0},
466  {1,0,0,1,1,0,1,1,0,0,0,1,1},{1,1,0,1,0,0,1,1,0,1,0,1,1},{2,1,2,2,1,0,1,1,0,0,1,2,1},
467  {2,0,1,1,0,0,2,2,0,2,2,1,2},{1,0,1,0,1,0,1,1,0,0,0,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},
468  {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,0,0,0,0,0,1,1,0,1,1,0,1},{1,0,0,0,0,1,1,1,1,1,0,1,0},
469  {1,1,0,0,1,1,1,1,1,0,0,1,0},{2,1,1,0,0,2,2,1,1,1,2,1,0},{2,0,2,0,2,1,1,2,2,0,1,2,0},
470  {1,0,1,1,0,1,1,1,1,1,0,0,0},{2,2,2,1,1,2,2,1,1,0,0,0,0},{2,2,0,2,0,1,1,2,2,2,1,0,0},
471  {2,0,0,1,1,2,2,1,1,0,2,0,0},{2,0,0,1,1,1,1,2,2,1,0,1,2},{2,2,0,2,0,2,2,1,1,0,0,2,1},
472  {4,3,2,2,3,4,4,1,1,3,4,2,1},{3,0,2,2,0,1,1,3,3,0,1,2,3},{2,0,2,0,2,2,2,1,1,2,0,0,1},
473  {2,1,1,0,0,1,1,2,2,0,0,0,2},{3,1,0,0,1,2,2,3,3,1,2,0,3},{2,0,0,0,0,1,1,2,2,0,1,0,2},
474  {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,1,0,0,1,1,0,1,0,1,1,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},
475  {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},{2,1,1,2,2,2,0,2,0,2,1,0,0},
476  {1,1,0,1,0,1,0,1,0,0,0,0,0},{1,0,0,1,1,1,0,1,0,1,0,0,0},{1,0,0,1,1,1,0,1,0,0,1,1,1},
477  {2,2,0,2,0,1,0,1,0,1,2,2,1},{2,2,1,1,2,2,0,2,0,0,0,1,2},{2,0,2,2,0,1,0,1,0,1,0,2,1},
478  {1,0,1,0,1,1,0,1,0,0,1,0,1},{2,2,2,0,0,1,0,1,0,1,2,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},
479  {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,0,0,0,0,0,0,1,1,1,1,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},
480  {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},
481  {2,2,2,1,1,0,0,1,1,0,2,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},{1,0,0,1,1,0,0,1,1,0,0,0,0},
482  {2,0,0,2,2,0,0,1,1,2,2,2,1},{2,1,0,1,0,0,0,2,2,0,1,1,2},{3,2,1,1,2,0,0,3,3,2,0,1,3},
483  {2,0,1,1,0,0,0,2,2,0,0,1,2},{2,0,1,0,1,0,0,2,2,1,1,0,2},{2,1,1,0,0,0,0,2,2,0,1,0,2},
484  {2,1,0,0,1,0,0,2,2,1,0,0,2},{1,0,0,0,0,0,0,1,1,0,0,0,1},{1,0,0,0,0,0,0,1,1,0,0,0,1},
485  {1,1,0,0,1,0,0,1,1,1,0,0,1},{2,1,1,0,0,0,0,2,2,0,1,0,2},{1,0,1,0,1,0,0,1,1,1,1,0,1},
486  {1,0,1,1,0,0,0,1,1,0,0,1,1},{2,1,1,2,2,0,0,1,1,1,0,1,2},{1,1,0,1,0,0,0,1,1,0,1,1,1},
487  {2,0,0,1,1,0,0,2,2,2,2,2,1},{1,0,0,1,1,0,0,1,1,0,0,0,0},{1,1,0,1,0,0,0,1,1,1,0,0,0},
488  {1,1,1,1,1,0,0,1,1,0,1,0,0},{1,0,1,1,0,0,0,1,1,1,1,0,0},{1,0,1,0,1,0,0,1,1,0,0,1,0},
489  {1,1,1,0,0,0,0,1,1,1,0,1,0},{1,1,0,0,1,0,0,1,1,0,1,1,0},{1,0,0,0,0,0,0,1,1,1,1,1,0},
490  {1,0,0,0,0,1,0,1,0,1,0,0,1},{1,1,0,0,1,1,0,1,0,0,0,0,1},{1,1,1,0,0,1,0,1,0,1,1,0,1},
491  {1,0,1,0,1,1,0,1,0,0,1,0,1},{1,0,1,1,0,1,0,1,0,1,0,1,1},{2,2,2,1,1,2,0,2,0,0,0,2,1},
492  {2,1,0,1,0,2,0,2,0,1,2,2,1},{2,0,0,2,2,1,0,1,0,0,1,1,2},{1,0,0,1,1,1,0,1,0,1,0,0,0},
493  {1,1,0,1,0,1,0,1,0,0,0,0,0},{2,1,2,2,1,2,0,2,0,1,2,0,0},{1,0,1,1,0,1,0,1,0,0,1,0,0},
494  {1,0,1,0,1,1,0,1,0,1,0,1,0},{1,1,1,0,0,1,0,1,0,0,0,1,0},{2,2,0,0,2,1,0,1,0,2,1,1,0},
495  {1,0,0,0,0,1,0,1,0,0,1,1,0},{1,0,0,0,0,1,1,1,1,0,1,0,1},{2,1,0,0,1,2,1,1,2,2,1,0,1},
496  {1,1,1,0,0,1,1,1,1,0,0,0,1},{2,0,2,0,2,1,2,2,1,1,0,0,2},{2,0,1,1,0,1,2,2,1,0,1,2,1},
497  {4,1,1,3,3,2,4,4,2,2,1,4,3},{2,2,0,2,0,2,1,1,2,0,0,1,2},{3,0,0,1,1,2,3,3,2,2,0,3,1},
498  {1,0,0,1,1,1,1,1,1,0,1,0,0},{2,2,0,2,0,1,2,2,1,1,2,0,0},{2,2,1,1,2,2,1,1,2,0,0,0,0},
499  {2,0,1,1,0,2,1,1,2,2,0,0,0},{2,0,2,0,2,2,1,1,2,0,2,1,0},{3,1,1,0,0,3,2,2,3,3,1,2,0},
500  {2,1,0,0,1,1,2,2,1,0,0,2,0},{2,0,0,0,0,2,1,1,2,2,0,1,0},{1,0,0,0,0,0,1,1,0,1,1,0,1},
501  {1,1,0,0,1,0,1,1,0,0,1,0,1},{1,1,1,0,0,0,1,1,0,1,0,0,1},{1,0,1,0,1,0,1,1,0,0,0,0,1},
502  {2,0,2,2,0,0,1,1,0,2,2,1,2},{3,1,1,2,2,0,3,3,0,0,1,3,2},{2,1,0,1,0,0,2,2,0,1,0,2,1},
503  {2,0,0,1,1,0,2,2,0,0,0,2,1},{1,0,0,1,1,0,1,1,0,1,1,0,0},{1,1,0,1,0,0,1,1,0,0,1,0,0},
504  {2,2,1,1,2,0,1,1,0,2,0,0,0},{1,0,1,1,0,0,1,1,0,0,0,0,0},{2,0,1,0,1,0,2,2,0,1,1,2,0},
505  {2,1,1,0,0,0,2,2,0,0,1,2,0},{2,1,0,0,1,0,2,2,0,1,0,2,0},{1,0,0,0,0,0,1,1,0,0,0,1,0},
506  {1,0,0,0,0,0,1,0,1,0,0,1,1},{1,1,0,0,1,0,1,0,1,1,0,1,1},{1,1,1,0,0,0,1,0,1,0,1,1,1},
507  {2,0,2,0,2,0,1,0,1,1,1,2,2},{1,0,1,1,0,0,1,0,1,0,0,0,1},{2,2,2,1,1,0,2,0,2,2,0,0,1},
508  {1,1,0,1,0,0,1,0,1,0,1,0,1},{2,0,0,2,2,0,1,0,1,1,1,0,2},{1,0,0,1,1,0,1,0,1,0,0,1,0},
509  {1,1,0,1,0,0,1,0,1,1,0,1,0},{2,2,1,1,2,0,2,0,2,0,2,1,0},{2,0,2,2,0,0,1,0,1,1,1,2,0},
510  {1,0,1,0,1,0,1,0,1,0,0,0,0},{1,1,1,0,0,0,1,0,1,1,0,0,0},{1,1,0,0,1,0,1,0,1,0,1,0,0},
511  {1,0,0,0,0,0,1,0,1,1,1,0,0},{1,0,0,0,0,1,1,0,0,1,0,1,1},{1,1,0,0,1,1,1,0,0,0,0,1,1},
512  {2,2,2,0,0,1,1,0,0,2,1,2,2},{2,0,1,0,1,2,2,0,0,0,2,1,1},{1,0,1,1,0,1,1,0,0,1,0,0,1},
513  {2,1,1,2,2,1,1,0,0,0,0,0,2},{2,1,0,1,0,2,2,0,0,1,2,0,1},{2,0,0,2,2,1,1,0,0,0,1,0,2},
514  {1,0,0,1,1,1,1,0,0,1,0,1,0},{1,1,0,1,0,1,1,0,0,0,0,1,0},{3,1,2,2,1,3,3,0,0,1,3,2,0},
515  {2,0,1,1,0,2,2,0,0,0,2,1,0},{1,0,1,0,1,1,1,0,0,1,0,0,0},{1,1,1,0,0,1,1,0,0,0,0,0,0},
516  {2,2,0,0,2,1,1,0,0,2,1,0,0},{1,0,0,0,0,1,1,0,0,0,1,0,0},{1,0,0,0,0,1,0,0,1,0,1,1,1},
517  {2,2,0,0,2,1,0,0,1,1,2,2,2},{1,1,1,0,0,1,0,0,1,0,0,1,1},{2,0,1,0,1,2,0,0,2,2,0,1,1},
518  {1,0,1,1,0,1,0,0,1,0,1,0,1},{3,1,1,3,3,2,0,0,2,2,1,0,3},{1,1,0,1,0,1,0,0,1,0,0,0,1},
519  {2,0,0,2,2,1,0,0,1,1,0,0,2},{1,0,0,1,1,1,0,0,1,0,1,1,0},{2,1,0,1,0,2,0,0,2,2,1,1,0},
520  {2,1,2,2,1,1,0,0,1,0,0,2,0},{2,0,1,1,0,2,0,0,2,2,0,1,0},{1,0,1,0,1,1,0,0,1,0,1,0,0},
521  {2,1,1,0,0,2,0,0,2,2,1,0,0},{1,1,0,0,1,1,0,0,1,0,0,0,0},{1,0,0,0,0,1,0,0,1,1,0,0,0},
522  {1,0,0,0,0,0,0,0,0,1,1,1,1},{1,1,0,0,1,0,0,0,0,0,1,1,1},{1,1,1,0,0,0,0,0,0,1,0,1,1},
523  {1,0,1,0,1,0,0,0,0,0,0,1,1},{1,0,1,1,0,0,0,0,0,1,1,0,1},{2,1,1,2,2,0,0,0,0,0,1,0,2},
524  {1,1,0,1,0,0,0,0,0,1,0,0,1},{1,0,0,1,1,0,0,0,0,0,0,0,1},{1,0,0,1,1,0,0,0,0,1,1,1,0},
525  {1,1,0,1,0,0,0,0,0,0,1,1,0},{2,1,2,2,1,0,0,0,0,1,0,2,0},{1,0,1,1,0,0,0,0,0,0,0,1,0},
526  {1,0,1,0,1,0,0,0,0,1,1,0,0},{1,1,1,0,0,0,0,0,0,0,1,0,0},{1,1,0,0,1,0,0,0,0,1,0,0,0},
527  {0,0,0,0,0,0,0,0,0,0,0,0,0}};
528 
529 ////////////////////////////////////////
530 
531 inline bool
532 isPlanarQuad(
533  const Vec3d& p0, const Vec3d& p1,
534  const Vec3d& p2, const Vec3d& p3,
535  const double epsilon = 0.001)
536 {
537  // compute representative plane
538  Vec3d normal = (p2-p0).cross(p1-p3);
539  normal.normalize();
540  const Vec3d centroid = (p0 + p1 + p2 + p3);
541  const double d = centroid.dot(normal) * 0.25;
542 
543 
544  // test vertice distance to plane
545  double absDist = std::abs(p0.dot(normal) - d);
546  if (absDist > epsilon) return false;
547 
548  absDist = std::abs(p1.dot(normal) - d);
549  if (absDist > epsilon) return false;
550 
551  absDist = std::abs(p2.dot(normal) - d);
552  if (absDist > epsilon) return false;
553 
554  absDist = std::abs(p3.dot(normal) - d);
555  if (absDist > epsilon) return false;
556 
557  return true;
558 }
559 
560 
561 ////////////////////////////////////////
562 
563 
564 /// @{
565 /// @brief Utility methods for point quantization.
566 
567 enum {
568  MASK_FIRST_10_BITS = 0x000003FF,
569  MASK_DIRTY_BIT = 0x80000000,
570  MASK_INVALID_BIT = 0x40000000
571 };
572 
573 inline uint32_t
574 packPoint(const Vec3d& v)
575 {
576  uint32_t data = 0;
577 
578  // values are expected to be in the [0.0 to 1.0] range.
579  OPENVDB_ASSERT(!(v.x() > 1.0) && !(v.y() > 1.0) && !(v.z() > 1.0));
580  OPENVDB_ASSERT(!(v.x() < 0.0) && !(v.y() < 0.0) && !(v.z() < 0.0));
581 
582  data |= (uint32_t(v.x() * 1023.0) & MASK_FIRST_10_BITS) << 20;
583  data |= (uint32_t(v.y() * 1023.0) & MASK_FIRST_10_BITS) << 10;
584  data |= (uint32_t(v.z() * 1023.0) & MASK_FIRST_10_BITS);
585 
586  return data;
587 }
588 
589 inline Vec3d
590 unpackPoint(uint32_t data)
591 {
592  Vec3d v;
593  v.z() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
594  data = data >> 10;
595  v.y() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
596  data = data >> 10;
597  v.x() = double(data & MASK_FIRST_10_BITS) * 0.0009775171;
598 
599  return v;
600 }
601 
602 /// @}
603 
604 ////////////////////////////////////////
605 
606 template<typename T>
607 inline bool isBoolValue() { return false; }
608 
609 template<>
610 inline bool isBoolValue<bool>() { return true; }
611 
612 template<typename T>
613 inline bool isInsideValue(T value, T isovalue) { return value < isovalue; }
614 
615 template<>
616 inline bool isInsideValue<bool>(bool value, bool /*isovalue*/) { return value; }
617 
618 
619 /// @brief Minor wrapper around the Leaf API to avoid atomic access with
620 /// delayed loading.
621 template <typename LeafT,
622  bool IsBool = std::is_same<typename LeafT::ValueType, bool>::value>
623 struct LeafBufferAccessor
624 {
625  using T = typename LeafT::ValueType;
626  LeafBufferAccessor(const LeafT& leaf) : mData(leaf.buffer().data()) {}
627  inline T get(const Index idx) const { return mData[idx]; }
628  const T* const mData;
629 };
630 
631 template <typename LeafT>
632 struct LeafBufferAccessor<LeafT, true>
633 {
634  using T = bool;
635  LeafBufferAccessor(const LeafT& leaf) : mLeaf(leaf) {}
636  inline T get(const Index idx) const { return mLeaf.getValue(idx); }
637  const LeafT& mLeaf;
638 };
639 
640 
641 /// @brief Whether a coordinate does not lie at the positive edge of a leaf node.
642 template <typename LeafT>
643 bool isInternalLeafCoord(const Coord& ijk)
644 {
645  return
646  ijk[0] < int(LeafT::DIM - 1) &&
647  ijk[1] < int(LeafT::DIM - 1) &&
648  ijk[2] < int(LeafT::DIM - 1);
649 }
650 
651 
652 /// @brief Extracts the eight corner values for a cell starting at the given @ijk coordinate.
653 template<typename AccessorT, typename ValueT>
654 inline void
655 getCellVertexValues(const AccessorT& accessor,
656  Coord ijk,
657  std::array<ValueT, 8>& values)
658 {
659  values[0] = ValueT(accessor.getValue(ijk)); // i, j, k
660  ++ijk[0];
661  values[1] = ValueT(accessor.getValue(ijk)); // i+1, j, k
662  ++ijk[2];
663  values[2] = ValueT(accessor.getValue(ijk)); // i+1, j, k+1
664  --ijk[0];
665  values[3] = ValueT(accessor.getValue(ijk)); // i, j, k+1
666  --ijk[2]; ++ijk[1];
667  values[4] = ValueT(accessor.getValue(ijk)); // i, j+1, k
668  ++ijk[0];
669  values[5] = ValueT(accessor.getValue(ijk)); // i+1, j+1, k
670  ++ijk[2];
671  values[6] = ValueT(accessor.getValue(ijk)); // i+1, j+1, k+1
672  --ijk[0];
673  values[7] = ValueT(accessor.getValue(ijk)); // i, j+1, k+1
674 }
675 
676 
677 /// @brief Extracts the eight corner values for a cell starting at the given @ijk coordinate.
678 template<typename LeafT, typename ValueT>
679 inline void
680 getCellVertexValues(const LeafT& leaf,
681  const Index offset,
682  std::array<ValueT, 8>& values)
683 {
684  const LeafBufferAccessor<LeafT> acc(leaf);
685 
686  values[0] = ValueT(acc.get(offset)); // i, j, k
687  values[3] = ValueT(acc.get(offset + 1)); // i, j, k+1
688  values[4] = ValueT(acc.get(offset + LeafT::DIM)); // i, j+1, k
689  values[7] = ValueT(acc.get(offset + LeafT::DIM + 1)); // i, j+1, k+1
690  values[1] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM))); // i+1, j, k
691  values[2] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1)); // i+1, j, k+1
692  values[5] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM)); // i+1, j+1, k
693  values[6] = ValueT(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1)); // i+1, j+1, k+1
694 }
695 
696 
697 template<typename ValueType>
698 inline uint8_t
699 computeSignFlags(const std::array<ValueType, 8>& values, const ValueType iso)
700 {
701  unsigned signs = 0;
702  signs |= isInsideValue(values[0], iso) ? 1u : 0u;
703  signs |= isInsideValue(values[1], iso) ? 2u : 0u;
704  signs |= isInsideValue(values[2], iso) ? 4u : 0u;
705  signs |= isInsideValue(values[3], iso) ? 8u : 0u;
706  signs |= isInsideValue(values[4], iso) ? 16u : 0u;
707  signs |= isInsideValue(values[5], iso) ? 32u : 0u;
708  signs |= isInsideValue(values[6], iso) ? 64u : 0u;
709  signs |= isInsideValue(values[7], iso) ? 128u : 0u;
710  return uint8_t(signs);
711 }
712 
713 
714 /// @brief General method that computes the cell-sign configuration at the given
715 /// @c ijk coordinate.
716 template<typename AccessorT>
717 inline uint8_t
718 evalCellSigns(const AccessorT& accessor, const Coord& ijk, typename AccessorT::ValueType iso)
719 {
720  unsigned signs = 0;
721  Coord coord = ijk; // i, j, k
722  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 1u;
723  coord[0] += 1; // i+1, j, k
724  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 2u;
725  coord[2] += 1; // i+1, j, k+1
726  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 4u;
727  coord[0] = ijk[0]; // i, j, k+1
728  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 8u;
729  coord[1] += 1; coord[2] = ijk[2]; // i, j+1, k
730  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 16u;
731  coord[0] += 1; // i+1, j+1, k
732  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 32u;
733  coord[2] += 1; // i+1, j+1, k+1
734  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 64u;
735  coord[0] = ijk[0]; // i, j+1, k+1
736  if (isInsideValue(accessor.getValue(coord), iso)) signs |= 128u;
737  return uint8_t(signs);
738 }
739 
740 
741 /// @brief Leaf node optimized method that computes the cell-sign configuration
742 /// at the given local @c offset
743 template<typename LeafT>
744 inline uint8_t
745 evalCellSigns(const LeafT& leaf, const Index offset, typename LeafT::ValueType iso)
746 {
747  const LeafBufferAccessor<LeafT> acc(leaf);
748 
749  unsigned signs = 0;
750  if (isInsideValue(acc.get(offset), iso)) signs |= 1u; // i, j, k
751  if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM)), iso)) signs |= 2u; // i+1, j, k
752  if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + 1), iso)) signs |= 4u; // i+1, j, k+1
753  if (isInsideValue(acc.get(offset + 1), iso)) signs |= 8u; // i, j, k+1
754  if (isInsideValue(acc.get(offset + LeafT::DIM), iso)) signs |= 16u; // i, j+1, k
755  if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM), iso)) signs |= 32u; // i+1, j+1, k
756  if (isInsideValue(acc.get(offset + (LeafT::DIM * LeafT::DIM) + LeafT::DIM + 1), iso)) signs |= 64u; // i+1, j+1, k+1
757  if (isInsideValue(acc.get(offset + LeafT::DIM + 1), iso)) signs |= 128u; // i, j+1, k+1
758  return uint8_t(signs);
759 }
760 
761 
762 /// @brief Used to correct topological ambiguities related to two adjacent cells
763 /// that share an ambiguous face.
764 template<class AccessorT>
765 inline void
766 correctCellSigns(uint8_t& signs,
767  const uint8_t face,
768  const AccessorT& acc,
769  Coord ijk,
770  const typename AccessorT::ValueType iso)
771 {
772  switch (int(face)) {
773  case 1:
774  ijk[2] -= 1;
775  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 3) signs = uint8_t(~signs);
776  break;
777  case 2:
778  ijk[0] += 1;
779  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 4) signs = uint8_t(~signs);
780  break;
781  case 3:
782  ijk[2] += 1;
783  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 1) signs = uint8_t(~signs);
784  break;
785  case 4:
786  ijk[0] -= 1;
787  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 2) signs = uint8_t(~signs);
788  break;
789  case 5:
790  ijk[1] -= 1;
791  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 6) signs = uint8_t(~signs);
792  break;
793  case 6:
794  ijk[1] += 1;
795  if (sAmbiguousFace[evalCellSigns(acc, ijk, iso)] == 5) signs = uint8_t(~signs);
796  break;
797  default:
798  break;
799  }
800 }
801 
802 
803 template<class AccessorT>
804 inline bool
805 isNonManifold(const AccessorT& accessor, const Coord& ijk,
806  typename AccessorT::ValueType isovalue, const int dim)
807 {
808  const int hDim = dim >> 1;
809  bool m, p[8]; // Corner signs
810 
811  Coord coord = ijk; // i, j, k
812  p[0] = isInsideValue(accessor.getValue(coord), isovalue);
813  coord[0] += dim; // i+dim, j, k
814  p[1] = isInsideValue(accessor.getValue(coord), isovalue);
815  coord[2] += dim; // i+dim, j, k+dim
816  p[2] = isInsideValue(accessor.getValue(coord), isovalue);
817  coord[0] = ijk[0]; // i, j, k+dim
818  p[3] = isInsideValue(accessor.getValue(coord), isovalue);
819  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
820  p[4] = isInsideValue(accessor.getValue(coord), isovalue);
821  coord[0] += dim; // i+dim, j+dim, k
822  p[5] = isInsideValue(accessor.getValue(coord), isovalue);
823  coord[2] += dim; // i+dim, j+dim, k+dim
824  p[6] = isInsideValue(accessor.getValue(coord), isovalue);
825  coord[0] = ijk[0]; // i, j+dim, k+dim
826  p[7] = isInsideValue(accessor.getValue(coord), isovalue);
827 
828  // Check if the corner sign configuration is ambiguous
829  unsigned signs = 0;
830  if (p[0]) signs |= 1u;
831  if (p[1]) signs |= 2u;
832  if (p[2]) signs |= 4u;
833  if (p[3]) signs |= 8u;
834  if (p[4]) signs |= 16u;
835  if (p[5]) signs |= 32u;
836  if (p[6]) signs |= 64u;
837  if (p[7]) signs |= 128u;
838  if (!sAdaptable[signs]) return true;
839 
840  // Manifold check
841 
842  // Evaluate edges
843  const int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
844  const int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
845  const int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
846 
847  // edge 1
848  coord.reset(ip, j, k);
849  m = isInsideValue(accessor.getValue(coord), isovalue);
850  if (p[0] != m && p[1] != m) return true;
851 
852  // edge 2
853  coord.reset(ipp, j, kp);
854  m = isInsideValue(accessor.getValue(coord), isovalue);
855  if (p[1] != m && p[2] != m) return true;
856 
857  // edge 3
858  coord.reset(ip, j, kpp);
859  m = isInsideValue(accessor.getValue(coord), isovalue);
860  if (p[2] != m && p[3] != m) return true;
861 
862  // edge 4
863  coord.reset(i, j, kp);
864  m = isInsideValue(accessor.getValue(coord), isovalue);
865  if (p[0] != m && p[3] != m) return true;
866 
867  // edge 5
868  coord.reset(ip, jpp, k);
869  m = isInsideValue(accessor.getValue(coord), isovalue);
870  if (p[4] != m && p[5] != m) return true;
871 
872  // edge 6
873  coord.reset(ipp, jpp, kp);
874  m = isInsideValue(accessor.getValue(coord), isovalue);
875  if (p[5] != m && p[6] != m) return true;
876 
877  // edge 7
878  coord.reset(ip, jpp, kpp);
879  m = isInsideValue(accessor.getValue(coord), isovalue);
880  if (p[6] != m && p[7] != m) return true;
881 
882  // edge 8
883  coord.reset(i, jpp, kp);
884  m = isInsideValue(accessor.getValue(coord), isovalue);
885  if (p[7] != m && p[4] != m) return true;
886 
887  // edge 9
888  coord.reset(i, jp, k);
889  m = isInsideValue(accessor.getValue(coord), isovalue);
890  if (p[0] != m && p[4] != m) return true;
891 
892  // edge 10
893  coord.reset(ipp, jp, k);
894  m = isInsideValue(accessor.getValue(coord), isovalue);
895  if (p[1] != m && p[5] != m) return true;
896 
897  // edge 11
898  coord.reset(ipp, jp, kpp);
899  m = isInsideValue(accessor.getValue(coord), isovalue);
900  if (p[2] != m && p[6] != m) return true;
901 
902 
903  // edge 12
904  coord.reset(i, jp, kpp);
905  m = isInsideValue(accessor.getValue(coord), isovalue);
906  if (p[3] != m && p[7] != m) return true;
907 
908 
909  // Evaluate faces
910 
911  // face 1
912  coord.reset(ip, jp, k);
913  m = isInsideValue(accessor.getValue(coord), isovalue);
914  if (p[0] != m && p[1] != m && p[4] != m && p[5] != m) return true;
915 
916  // face 2
917  coord.reset(ipp, jp, kp);
918  m = isInsideValue(accessor.getValue(coord), isovalue);
919  if (p[1] != m && p[2] != m && p[5] != m && p[6] != m) return true;
920 
921  // face 3
922  coord.reset(ip, jp, kpp);
923  m = isInsideValue(accessor.getValue(coord), isovalue);
924  if (p[2] != m && p[3] != m && p[6] != m && p[7] != m) return true;
925 
926  // face 4
927  coord.reset(i, jp, kp);
928  m = isInsideValue(accessor.getValue(coord), isovalue);
929  if (p[0] != m && p[3] != m && p[4] != m && p[7] != m) return true;
930 
931  // face 5
932  coord.reset(ip, j, kp);
933  m = isInsideValue(accessor.getValue(coord), isovalue);
934  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m) return true;
935 
936  // face 6
937  coord.reset(ip, jpp, kp);
938  m = isInsideValue(accessor.getValue(coord), isovalue);
939  if (p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
940 
941  // test cube center
942  coord.reset(ip, jp, kp);
943  m = isInsideValue(accessor.getValue(coord), isovalue);
944  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
945  p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
946 
947  return false;
948 }
949 
950 
951 ////////////////////////////////////////
952 
953 
954 template <class LeafType>
955 inline void
956 mergeVoxels(LeafType& leaf, const Coord& start, const int dim, const int regionId)
957 {
958  Coord ijk;
959  const Coord end = start.offsetBy(dim);
960 
961  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
962  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
963  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
964  leaf.setValueOnly(ijk, regionId);
965  }
966  }
967  }
968 }
969 
970 
971 // Note that we must use ValueType::value_type or else Visual C++ gets confused
972 // thinking that it is a constructor.
973 template <class LeafType>
974 inline bool
975 isMergeable(const LeafType& leaf,
976  const Coord& start,
977  const int dim,
978  typename LeafType::ValueType::value_type adaptivity)
979 {
980  if (adaptivity < 1e-6) return false;
981 
982  using VecT = typename LeafType::ValueType;
983  Coord ijk;
984  const Coord end = start.offsetBy(dim);
985 
986  std::vector<VecT> norms;
987  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
988  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
989  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
990  if (!leaf.isValueOn(ijk)) continue;
991  norms.push_back(leaf.getValue(ijk));
992  }
993  }
994  }
995 
996  const size_t N = norms.size();
997  for (size_t ni = 0; ni < N; ++ni) {
998  VecT n_i = norms[ni];
999  for (size_t nj = 0; nj < N; ++nj) {
1000  VecT n_j = norms[nj];
1001  if ((1.0 - n_i.dot(n_j)) > adaptivity) return false;
1002  }
1003  }
1004  return true;
1005 }
1006 
1007 
1008 ////////////////////////////////////////
1009 
1010 
1011 /// linear interpolation.
1012 inline double evalZeroCrossing(double v0, double v1, double iso) { return (iso - v0) / (v1 - v0); }
1013 
1014 
1015 /// @brief Computes the average cell point for a given edge group.
1016 inline Vec3d
1017 computePoint(const std::array<double, 8>& values,
1018  const unsigned char signs,
1019  const unsigned char edgeGroup,
1020  const double iso)
1021 {
1022  Vec3d avg(0.0, 0.0, 0.0);
1023  int samples = 0;
1024 
1025  if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
1026  avg[0] += evalZeroCrossing(values[0], values[1], iso);
1027  ++samples;
1028  }
1029 
1030  if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
1031  avg[0] += 1.0;
1032  avg[2] += evalZeroCrossing(values[1], values[2], iso);
1033  ++samples;
1034  }
1035 
1036  if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
1037  avg[0] += evalZeroCrossing(values[3], values[2], iso);
1038  avg[2] += 1.0;
1039  ++samples;
1040  }
1041 
1042  if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
1043  avg[2] += evalZeroCrossing(values[0], values[3], iso);
1044  ++samples;
1045  }
1046 
1047  if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
1048  avg[0] += evalZeroCrossing(values[4], values[5], iso);
1049  avg[1] += 1.0;
1050  ++samples;
1051  }
1052 
1053  if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
1054  avg[0] += 1.0;
1055  avg[1] += 1.0;
1056  avg[2] += evalZeroCrossing(values[5], values[6], iso);
1057  ++samples;
1058  }
1059 
1060  if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
1061  avg[0] += evalZeroCrossing(values[7], values[6], iso);
1062  avg[1] += 1.0;
1063  avg[2] += 1.0;
1064  ++samples;
1065  }
1066 
1067  if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
1068  avg[1] += 1.0;
1069  avg[2] += evalZeroCrossing(values[4], values[7], iso);
1070  ++samples;
1071  }
1072 
1073  if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
1074  avg[1] += evalZeroCrossing(values[0], values[4], iso);
1075  ++samples;
1076  }
1077 
1078  if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
1079  avg[0] += 1.0;
1080  avg[1] += evalZeroCrossing(values[1], values[5], iso);
1081  ++samples;
1082  }
1083 
1084  if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
1085  avg[0] += 1.0;
1086  avg[1] += evalZeroCrossing(values[2], values[6], iso);
1087  avg[2] += 1.0;
1088  ++samples;
1089  }
1090 
1091  if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
1092  avg[1] += evalZeroCrossing(values[3], values[7], iso);
1093  avg[2] += 1.0;
1094  ++samples;
1095  }
1096 
1097  if (samples > 1) {
1098  const double w = 1.0 / double(samples);
1099  avg *= w;
1100  }
1101 
1102  return avg;
1103 }
1104 
1105 
1106 /// @brief Computes the average cell point for a given edge group, ignoring edge
1107 /// samples present in the @c signsMask configuration.
1108 inline int
1109 computeMaskedPoint(Vec3d& avg,
1110  const std::array<double, 8>& values,
1111  const unsigned char signs,
1112  const unsigned char signsMask,
1113  const unsigned char edgeGroup,
1114  const double iso)
1115 {
1116  avg = Vec3d(0.0, 0.0, 0.0);
1117  int samples = 0;
1118 
1119  if (sEdgeGroupTable[signs][1] == edgeGroup &&
1120  sEdgeGroupTable[signsMask][1] == 0) { // Edged: 0 - 1
1121  avg[0] += evalZeroCrossing(values[0], values[1], iso);
1122  ++samples;
1123  }
1124 
1125  if (sEdgeGroupTable[signs][2] == edgeGroup &&
1126  sEdgeGroupTable[signsMask][2] == 0) { // Edged: 1 - 2
1127  avg[0] += 1.0;
1128  avg[2] += evalZeroCrossing(values[1], values[2], iso);
1129  ++samples;
1130  }
1131 
1132  if (sEdgeGroupTable[signs][3] == edgeGroup &&
1133  sEdgeGroupTable[signsMask][3] == 0) { // Edged: 3 - 2
1134  avg[0] += evalZeroCrossing(values[3], values[2], iso);
1135  avg[2] += 1.0;
1136  ++samples;
1137  }
1138 
1139  if (sEdgeGroupTable[signs][4] == edgeGroup &&
1140  sEdgeGroupTable[signsMask][4] == 0) { // Edged: 0 - 3
1141  avg[2] += evalZeroCrossing(values[0], values[3], iso);
1142  ++samples;
1143  }
1144 
1145  if (sEdgeGroupTable[signs][5] == edgeGroup &&
1146  sEdgeGroupTable[signsMask][5] == 0) { // Edged: 4 - 5
1147  avg[0] += evalZeroCrossing(values[4], values[5], iso);
1148  avg[1] += 1.0;
1149  ++samples;
1150  }
1151 
1152  if (sEdgeGroupTable[signs][6] == edgeGroup &&
1153  sEdgeGroupTable[signsMask][6] == 0) { // Edged: 5 - 6
1154  avg[0] += 1.0;
1155  avg[1] += 1.0;
1156  avg[2] += evalZeroCrossing(values[5], values[6], iso);
1157  ++samples;
1158  }
1159 
1160  if (sEdgeGroupTable[signs][7] == edgeGroup &&
1161  sEdgeGroupTable[signsMask][7] == 0) { // Edged: 7 - 6
1162  avg[0] += evalZeroCrossing(values[7], values[6], iso);
1163  avg[1] += 1.0;
1164  avg[2] += 1.0;
1165  ++samples;
1166  }
1167 
1168  if (sEdgeGroupTable[signs][8] == edgeGroup &&
1169  sEdgeGroupTable[signsMask][8] == 0) { // Edged: 4 - 7
1170  avg[1] += 1.0;
1171  avg[2] += evalZeroCrossing(values[4], values[7], iso);
1172  ++samples;
1173  }
1174 
1175  if (sEdgeGroupTable[signs][9] == edgeGroup &&
1176  sEdgeGroupTable[signsMask][9] == 0) { // Edged: 0 - 4
1177  avg[1] += evalZeroCrossing(values[0], values[4], iso);
1178  ++samples;
1179  }
1180 
1181  if (sEdgeGroupTable[signs][10] == edgeGroup &&
1182  sEdgeGroupTable[signsMask][10] == 0) { // Edged: 1 - 5
1183  avg[0] += 1.0;
1184  avg[1] += evalZeroCrossing(values[1], values[5], iso);
1185  ++samples;
1186  }
1187 
1188  if (sEdgeGroupTable[signs][11] == edgeGroup &&
1189  sEdgeGroupTable[signsMask][11] == 0) { // Edged: 2 - 6
1190  avg[0] += 1.0;
1191  avg[1] += evalZeroCrossing(values[2], values[6], iso);
1192  avg[2] += 1.0;
1193  ++samples;
1194  }
1195 
1196  if (sEdgeGroupTable[signs][12] == edgeGroup &&
1197  sEdgeGroupTable[signsMask][12] == 0) { // Edged: 3 - 7
1198  avg[1] += evalZeroCrossing(values[3], values[7], iso);
1199  avg[2] += 1.0;
1200  ++samples;
1201  }
1202 
1203  if (samples > 1) {
1204  const double w = 1.0 / double(samples);
1205  avg *= w;
1206  }
1207 
1208  return samples;
1209 }
1210 
1211 
1212 /// @brief Computes the average cell point for a given edge group, by computing
1213 /// convex weights based on the distance from the sample point @c p.
1214 inline Vec3d
1215 computeWeightedPoint(const Vec3d& p,
1216  const std::array<double, 8>& values,
1217  const unsigned char signs,
1218  const unsigned char edgeGroup,
1219  const double iso)
1220 {
1221  std::vector<Vec3d> samples;
1222  samples.reserve(8);
1223 
1224  Vec3d avg(0.0, 0.0, 0.0);
1225 
1226  if (sEdgeGroupTable[signs][1] == edgeGroup) { // Edged: 0 - 1
1227  avg[0] = evalZeroCrossing(values[0], values[1], iso);
1228  avg[1] = 0.0;
1229  avg[2] = 0.0;
1230 
1231  samples.push_back(avg);
1232  }
1233 
1234  if (sEdgeGroupTable[signs][2] == edgeGroup) { // Edged: 1 - 2
1235  avg[0] = 1.0;
1236  avg[1] = 0.0;
1237  avg[2] = evalZeroCrossing(values[1], values[2], iso);
1238 
1239  samples.push_back(avg);
1240  }
1241 
1242  if (sEdgeGroupTable[signs][3] == edgeGroup) { // Edged: 3 - 2
1243  avg[0] = evalZeroCrossing(values[3], values[2], iso);
1244  avg[1] = 0.0;
1245  avg[2] = 1.0;
1246 
1247  samples.push_back(avg);
1248  }
1249 
1250  if (sEdgeGroupTable[signs][4] == edgeGroup) { // Edged: 0 - 3
1251  avg[0] = 0.0;
1252  avg[1] = 0.0;
1253  avg[2] = evalZeroCrossing(values[0], values[3], iso);
1254 
1255  samples.push_back(avg);
1256  }
1257 
1258  if (sEdgeGroupTable[signs][5] == edgeGroup) { // Edged: 4 - 5
1259  avg[0] = evalZeroCrossing(values[4], values[5], iso);
1260  avg[1] = 1.0;
1261  avg[2] = 0.0;
1262 
1263  samples.push_back(avg);
1264  }
1265 
1266  if (sEdgeGroupTable[signs][6] == edgeGroup) { // Edged: 5 - 6
1267  avg[0] = 1.0;
1268  avg[1] = 1.0;
1269  avg[2] = evalZeroCrossing(values[5], values[6], iso);
1270 
1271  samples.push_back(avg);
1272  }
1273 
1274  if (sEdgeGroupTable[signs][7] == edgeGroup) { // Edged: 7 - 6
1275  avg[0] = evalZeroCrossing(values[7], values[6], iso);
1276  avg[1] = 1.0;
1277  avg[2] = 1.0;
1278 
1279  samples.push_back(avg);
1280  }
1281 
1282  if (sEdgeGroupTable[signs][8] == edgeGroup) { // Edged: 4 - 7
1283  avg[0] = 0.0;
1284  avg[1] = 1.0;
1285  avg[2] = evalZeroCrossing(values[4], values[7], iso);
1286 
1287  samples.push_back(avg);
1288  }
1289 
1290  if (sEdgeGroupTable[signs][9] == edgeGroup) { // Edged: 0 - 4
1291  avg[0] = 0.0;
1292  avg[1] = evalZeroCrossing(values[0], values[4], iso);
1293  avg[2] = 0.0;
1294 
1295  samples.push_back(avg);
1296  }
1297 
1298  if (sEdgeGroupTable[signs][10] == edgeGroup) { // Edged: 1 - 5
1299  avg[0] = 1.0;
1300  avg[1] = evalZeroCrossing(values[1], values[5], iso);
1301  avg[2] = 0.0;
1302 
1303  samples.push_back(avg);
1304  }
1305 
1306  if (sEdgeGroupTable[signs][11] == edgeGroup) { // Edged: 2 - 6
1307  avg[0] = 1.0;
1308  avg[1] = evalZeroCrossing(values[2], values[6], iso);
1309  avg[2] = 1.0;
1310 
1311  samples.push_back(avg);
1312  }
1313 
1314  if (sEdgeGroupTable[signs][12] == edgeGroup) { // Edged: 3 - 7
1315  avg[0] = 0.0;
1316  avg[1] = evalZeroCrossing(values[3], values[7], iso);
1317  avg[2] = 1.0;
1318 
1319  samples.push_back(avg);
1320  }
1321 
1322  OPENVDB_ASSERT(!samples.empty());
1323  if (samples.size() == 1) {
1324  return samples.front();
1325  }
1326 
1327  std::vector<double> weights;
1328  weights.reserve(samples.size());
1329 
1330  for (const Vec3d& s : samples) {
1331  weights.emplace_back((s-p).lengthSqr());
1332  }
1333 
1334  double minWeight = weights.front();
1335  double maxWeight = weights.front();
1336 
1337  for (size_t i = 1, I = weights.size(); i < I; ++i) {
1338  minWeight = std::min(minWeight, weights[i]);
1339  maxWeight = std::max(maxWeight, weights[i]);
1340  }
1341 
1342  const double offset = maxWeight + minWeight * 0.1;
1343  for (size_t i = 0, I = weights.size(); i < I; ++i) {
1344  weights[i] = offset - weights[i];
1345  }
1346 
1347  double weightSum = 0.0;
1348  for (size_t i = 0, I = weights.size(); i < I; ++i) {
1349  weightSum += weights[i];
1350  }
1351 
1352  avg.setZero();
1353  for (size_t i = 0, I = samples.size(); i < I; ++i) {
1354  avg += samples[i] * (weights[i] / weightSum);
1355  }
1356 
1357  return avg;
1358 }
1359 
1360 
1361 /// @brief Computes the average cell points defined by the sign configuration
1362 /// @c signs and the given corner values @c values.
1363 inline size_t
1364 computeCellPoints(std::array<Vec3d, 4>& points,
1365  const std::array<double, 8>& values,
1366  const unsigned char signs,
1367  const double iso)
1368 {
1369  size_t offset = 0;
1370  for (size_t n = 1, N = sEdgeGroupTable[signs][0] + 1; n < N; ++n, ++offset) {
1371  OPENVDB_ASSERT(offset < 4);
1372  points[offset] = computePoint(values, signs, uint8_t(n), iso);
1373  }
1374  return offset;
1375 }
1376 
1377 
1378 /// @brief Given a sign configuration @c lhsSigns and an edge group @c groupId,
1379 /// finds the corresponding edge group in a different sign configuration
1380 /// @c rhsSigns. Returns -1 if no match is found.
1381 inline int
1382 matchEdgeGroup(unsigned char groupId, unsigned char lhsSigns, unsigned char rhsSigns)
1383 {
1384  int id = -1;
1385  for (size_t i = 1; i <= 12; ++i) {
1386  if (sEdgeGroupTable[lhsSigns][i] == groupId && sEdgeGroupTable[rhsSigns][i] != 0) {
1387  id = sEdgeGroupTable[rhsSigns][i];
1388  break;
1389  }
1390  }
1391  return id;
1392 }
1393 
1394 
1395 /// @brief Computes the average cell points defined by the sign configuration
1396 /// @c signs and the given corner values @c values. Combines data from
1397 /// two different level sets to eliminate seam lines when meshing
1398 /// fractured segments.
1399 inline size_t
1400 computeCellPoints(std::array<Vec3d, 4>& points,
1401  std::array<bool, 4>& weightedPointMask,
1402  const std::array<double, 8>& lhsValues,
1403  const std::array<double, 8>& rhsValues,
1404  const unsigned char lhsSigns,
1405  const unsigned char rhsSigns,
1406  const double iso,
1407  const size_t pointIdx,
1408  const uint32_t * seamPointArray)
1409 {
1410  size_t offset = 0;
1411  for (size_t n = 1, N = sEdgeGroupTable[lhsSigns][0] + 1; n < N; ++n, ++offset)
1412  {
1413  OPENVDB_ASSERT(offset < 4);
1414  const int id = matchEdgeGroup(uint8_t(n), lhsSigns, rhsSigns);
1415 
1416  if (id != -1) {
1417 
1418  const unsigned char e = uint8_t(id);
1419  const uint32_t quantizedPoint = seamPointArray[pointIdx + (id - 1)];
1420 
1421  if ((quantizedPoint & MASK_DIRTY_BIT) && !(quantizedPoint & MASK_INVALID_BIT)) {
1422  const Vec3d p = unpackPoint(quantizedPoint);
1423  points[offset] = computeWeightedPoint(p, rhsValues, rhsSigns, e, iso);
1424  weightedPointMask[offset] = true;
1425  } else {
1426  points[offset] = computePoint(rhsValues, rhsSigns, e, iso);
1427  weightedPointMask[offset] = false;
1428  }
1429 
1430  } else {
1431  points[offset] = computePoint(lhsValues, lhsSigns, uint8_t(n), iso);
1432  weightedPointMask[offset] = false;
1433  }
1434  }
1435  return offset;
1436 }
1437 
1438 
1439 template <typename InputTreeType>
1440 struct ComputePoints
1441 {
1442  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
1443  using InputValueType = typename InputLeafNodeType::ValueType;
1444 
1445  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
1446  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
1447 
1448  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
1449  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
1450 
1451  ComputePoints(Vec3s * pointArray,
1452  const InputTreeType& inputTree,
1453  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1454  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1455  const std::unique_ptr<Index32[]>& leafNodeOffsets,
1456  const math::Transform& xform,
1457  const double iso);
1458 
1459  void setRefData(const InputTreeType& refInputTree,
1460  const Index32TreeType& refPointIndexTree,
1461  const Int16TreeType& refSignFlagsTree,
1462  const uint32_t * quantizedSeamLinePoints,
1463  uint8_t * seamLinePointsFlags);
1464 
1465  void operator()(const tbb::blocked_range<size_t>&) const;
1466 
1467 private:
1468  Vec3s * const mPoints;
1469  InputTreeType const * const mInputTree;
1470  Index32LeafNodeType * const * const mPointIndexNodes;
1471  Int16LeafNodeType const * const * const mSignFlagsNodes;
1472  Index32 const * const mNodeOffsets;
1473  math::Transform const mTransform;
1474  double const mIsovalue;
1475  // reference meshing data
1476  InputTreeType const * mRefInputTree;
1477  Index32TreeType const * mRefPointIndexTree;
1478  Int16TreeType const * mRefSignFlagsTree;
1479  uint32_t const * mQuantizedSeamLinePoints;
1480  uint8_t * mSeamLinePointsFlags;
1481 }; // struct ComputePoints
1482 
1483 
1484 template <typename InputTreeType>
1485 ComputePoints<InputTreeType>::ComputePoints(
1486  Vec3s * pointArray,
1487  const InputTreeType& inputTree,
1488  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
1489  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1490  const std::unique_ptr<Index32[]>& leafNodeOffsets,
1491  const math::Transform& xform,
1492  const double iso)
1493  : mPoints(pointArray)
1494  , mInputTree(&inputTree)
1495  , mPointIndexNodes(pointIndexLeafNodes.data())
1496  , mSignFlagsNodes(signFlagsLeafNodes.data())
1497  , mNodeOffsets(leafNodeOffsets.get())
1498  , mTransform(xform)
1499  , mIsovalue(iso)
1500  , mRefInputTree(nullptr)
1501  , mRefPointIndexTree(nullptr)
1502  , mRefSignFlagsTree(nullptr)
1503  , mQuantizedSeamLinePoints(nullptr)
1504  , mSeamLinePointsFlags(nullptr)
1505 {
1506 }
1507 
1508 template <typename InputTreeType>
1509 void
1510 ComputePoints<InputTreeType>::setRefData(
1511  const InputTreeType& refInputTree,
1512  const Index32TreeType& refPointIndexTree,
1513  const Int16TreeType& refSignFlagsTree,
1514  const uint32_t * quantizedSeamLinePoints,
1515  uint8_t * seamLinePointsFlags)
1516 {
1517  mRefInputTree = &refInputTree;
1518  mRefPointIndexTree = &refPointIndexTree;
1519  mRefSignFlagsTree = &refSignFlagsTree;
1520  mQuantizedSeamLinePoints = quantizedSeamLinePoints;
1521  mSeamLinePointsFlags = seamLinePointsFlags;
1522 }
1523 
1524 template <typename InputTreeType>
1525 void
1526 ComputePoints<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range) const
1527 {
1528  using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
1529  using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
1530  using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
1531 
1532  using IndexType = typename Index32TreeType::ValueType;
1533 
1534  using IndexArray = std::vector<Index>;
1535  using IndexArrayMap = std::map<IndexType, IndexArray>;
1536 
1537  InputTreeAccessor inputAcc(*mInputTree);
1538 
1539  Vec3d xyz;
1540  Coord ijk;
1541  std::array<Vec3d, 4> points;
1542  std::array<bool, 4> weightedPointMask;
1543  std::array<double, 8> values, refValues;
1544  const double iso = mIsovalue;
1545 
1546  // reference data accessors
1547 
1548  std::unique_ptr<InputTreeAccessor> refInputAcc;
1549  std::unique_ptr<Index32TreeAccessor> refPointIndexAcc;
1550  std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
1551 
1552  const bool hasReferenceData = mRefInputTree && mRefPointIndexTree && mRefSignFlagsTree;
1553 
1554  if (hasReferenceData) {
1555  refInputAcc.reset(new InputTreeAccessor(*mRefInputTree));
1556  refPointIndexAcc.reset(new Index32TreeAccessor(*mRefPointIndexTree));
1557  refSignFlagsAcc.reset(new Int16TreeAccessor(*mRefSignFlagsTree));
1558  }
1559 
1560  for (size_t n = range.begin(), N = range.end(); n != N; ++n)
1561  {
1562  Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[n];
1563  const Coord& origin = pointIndexNode.origin();
1564 
1565  const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1566  const InputLeafNodeType * inputNode = inputAcc.probeConstLeaf(origin);
1567 
1568  // get reference data
1569  const InputLeafNodeType * refInputNode = nullptr;
1570  const Index32LeafNodeType * refPointIndexNode = nullptr;
1571  const Int16LeafNodeType * refSignFlagsNode = nullptr;
1572 
1573  if (hasReferenceData) {
1574  refInputNode = refInputAcc->probeConstLeaf(origin);
1575  refPointIndexNode = refPointIndexAcc->probeConstLeaf(origin);
1576  refSignFlagsNode = refSignFlagsAcc->probeConstLeaf(origin);
1577  }
1578 
1579  IndexType pointOffset = IndexType(mNodeOffsets[n]);
1580  IndexArrayMap regions;
1581 
1582  auto* const pidxData = pointIndexNode.buffer().data();
1583  const auto* const sfData = signFlagsNode.buffer().data();
1584 
1585  for (auto it = pointIndexNode.beginValueOn(); it; ++it)
1586  {
1587  const Index offset = it.pos();
1588  IndexType& id = pidxData[offset];
1589 
1590  if (id != 0) {
1591  if (id != IndexType(util::INVALID_IDX)) {
1592  regions[id].push_back(offset);
1593  }
1594  continue;
1595  }
1596 
1597  id = pointOffset;
1598 
1599  const Int16 flags = sfData[offset];
1600  const uint8_t signs = uint8_t(SIGNS & flags);
1601  uint8_t refSigns = 0;
1602 
1603  if ((flags & SEAM) && refPointIndexNode && refSignFlagsNode) {
1604  if (refSignFlagsNode->isValueOn(offset)) {
1605  refSigns = uint8_t(SIGNS & refSignFlagsNode->getValue(offset));
1606  }
1607  }
1608 
1609  ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1610 
1611  const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1612 
1613  ijk += origin;
1614 
1615  if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1616  else getCellVertexValues(inputAcc, ijk, values);
1617 
1618  size_t count, weightcount;
1619 
1620  if (refSigns == 0) {
1621  count = computeCellPoints(points, values, signs, iso);
1622  weightcount = 0;
1623  } else {
1624  if (inclusiveCell && refInputNode) {
1625  getCellVertexValues(*refInputNode, offset, refValues);
1626  } else {
1627  getCellVertexValues(*refInputAcc, ijk, refValues);
1628  }
1629  count = computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1630  iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1631  weightcount = count;
1632  }
1633 
1634  xyz = ijk.asVec3d();
1635 
1636  for (size_t i = 0; i < count; ++i) {
1637 
1638  Vec3d& point = points[i];
1639 
1640  // Checks for both NaN and inf vertex positions, i.e. any value that is not finite.
1641  if (!std::isfinite(point[0]) ||
1642  !std::isfinite(point[1]) ||
1643  !std::isfinite(point[2]))
1644  {
1646  "VolumeToMesh encountered NaNs or infs in the input VDB!"
1647  " Hint: Check the input and consider using the \"Diagnostics\" tool "
1648  "to detect and resolve the NaNs.");
1649  }
1650 
1651  point += xyz;
1652  point = mTransform.indexToWorld(point);
1653 
1654  Vec3s& pos = mPoints[pointOffset];
1655  pos[0] = float(point[0]);
1656  pos[1] = float(point[1]);
1657  pos[2] = float(point[2]);
1658 
1659  if (mSeamLinePointsFlags && weightcount && weightedPointMask[i]) {
1660  mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1661  }
1662 
1663  ++pointOffset;
1664  }
1665  }
1666 
1667  // generate collapsed region points
1668  for (auto it = regions.begin(); it != regions.end(); ++it)
1669  {
1670  Vec3d avg(0.0);
1671 
1672  const IndexArray& voxels = it->second;
1673  for (size_t i = 0, I = voxels.size(); i < I; ++i) {
1674 
1675  const Index offset = voxels[i];
1676  ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1677 
1678  const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1679 
1680  ijk += origin;
1681 
1682  pidxData[offset] = pointOffset;
1683 
1684  const uint8_t signs = uint8_t(SIGNS & sfData[offset]);
1685 
1686  if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1687  else getCellVertexValues(inputAcc, ijk, values);
1688 
1689  computeCellPoints(points, values, signs, iso);
1690 
1691  avg[0] += double(ijk[0]) + points[0][0];
1692  avg[1] += double(ijk[1]) + points[0][1];
1693  avg[2] += double(ijk[2]) + points[0][2];
1694  }
1695 
1696  if (voxels.size() > 1) {
1697  const double w = 1.0 / double(voxels.size());
1698  avg *= w;
1699  }
1700 
1701  avg = mTransform.indexToWorld(avg);
1702 
1703  Vec3s& pos = mPoints[pointOffset];
1704  pos[0] = float(avg[0]);
1705  pos[1] = float(avg[1]);
1706  pos[2] = float(avg[2]);
1707 
1708  ++pointOffset;
1709  }
1710  }
1711 } // ComputePoints::operator()
1712 
1713 
1714 ////////////////////////////////////////
1715 
1716 
1717 template <typename InputTreeType>
1718 struct SeamLineWeights
1719 {
1720  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
1721  using InputValueType = typename InputLeafNodeType::ValueType;
1722 
1723  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
1724  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
1725 
1726  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
1727  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
1728 
1729  SeamLineWeights(const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1730  const InputTreeType& inputTree,
1731  const Index32TreeType& refPointIndexTree,
1732  const Int16TreeType& refSignFlagsTree,
1733  uint32_t * quantizedPoints,
1734  InputValueType iso)
1735  : mSignFlagsNodes(signFlagsLeafNodes.data())
1736  , mInputTree(&inputTree)
1737  , mRefPointIndexTree(&refPointIndexTree)
1738  , mRefSignFlagsTree(&refSignFlagsTree)
1739  , mQuantizedPoints(quantizedPoints)
1740  , mIsovalue(iso)
1741  {
1742  }
1743 
1744  void operator()(const tbb::blocked_range<size_t>& range) const
1745  {
1746  tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
1747  tree::ValueAccessor<const Index32TreeType> pointIndexTreeAcc(*mRefPointIndexTree);
1748  tree::ValueAccessor<const Int16TreeType> signFlagsTreeAcc(*mRefSignFlagsTree);
1749 
1750  std::array<double, 8> values;
1751  const double iso = double(mIsovalue);
1752  Coord ijk;
1753  Vec3d pos;
1754 
1755  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1756 
1757  const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1758  const Coord& origin = signFlagsNode.origin();
1759 
1760  const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.probeConstLeaf(origin);
1761  if (!refSignNode) continue;
1762 
1763  const Index32LeafNodeType* refPointIndexNode =
1764  pointIndexTreeAcc.probeConstLeaf(origin);
1765  if (!refPointIndexNode) continue;
1766 
1767  const InputLeafNodeType* inputNode = inputTreeAcc.probeConstLeaf(origin);
1768 
1769  const auto* const sfData = signFlagsNode.buffer().data();
1770  const auto* const rfIdxData = refPointIndexNode->buffer().data();
1771  const auto* const rsData = refSignNode->buffer().data();
1772 
1773  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it)
1774  {
1775  const Index offset = it.pos();
1776  const Int16 flags = sfData[offset];
1777 
1778  ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1779 
1780  const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1781 
1782  ijk += origin;
1783 
1784  if ((flags & SEAM) && refSignNode->isValueOn(offset)) {
1785 
1786  const uint8_t lhsSigns = uint8_t(SIGNS & flags);
1787  const uint8_t rhsSigns = uint8_t(SIGNS & rsData[offset]);
1788 
1789  if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1790  else getCellVertexValues(inputTreeAcc, ijk, values);
1791 
1792  for (unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1793 
1794  const int id = matchEdgeGroup(uint8_t(i), lhsSigns, rhsSigns);
1795 
1796  if (id != -1) {
1797 
1798  uint32_t& data = mQuantizedPoints[rfIdxData[offset] + (id - 1)];
1799 
1800  if (!(data & MASK_DIRTY_BIT)) {
1801 
1802  const int samples = computeMaskedPoint(
1803  pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1804 
1805  if (samples > 0) data = packPoint(pos);
1806  else data = MASK_INVALID_BIT;
1807 
1808  data |= MASK_DIRTY_BIT;
1809  }
1810  }
1811  } // end point group loop
1812  }
1813  } // end value on loop
1814  } // end leaf node loop
1815  }
1816 
1817 private:
1818  Int16LeafNodeType const * const * const mSignFlagsNodes;
1819  InputTreeType const * const mInputTree;
1820  Index32TreeType const * const mRefPointIndexTree;
1821  Int16TreeType const * const mRefSignFlagsTree;
1822  uint32_t * const mQuantizedPoints;
1823  InputValueType const mIsovalue;
1824 }; // struct SeamLineWeights
1825 
1826 
1827 template <typename TreeType>
1828 struct SetSeamLineFlags
1829 {
1830  using LeafNodeType = typename TreeType::LeafNodeType;
1831 
1832  SetSeamLineFlags(const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1833  const TreeType& refSignFlagsTree)
1834  : mSignFlagsNodes(signFlagsLeafNodes.data())
1835  , mRefSignFlagsTree(&refSignFlagsTree)
1836  {
1837  }
1838 
1839  void operator()(const tbb::blocked_range<size_t>& range) const
1840  {
1841  tree::ValueAccessor<const TreeType> refSignFlagsTreeAcc(*mRefSignFlagsTree);
1842 
1843  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1844 
1845  LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1846  const Coord& origin = signFlagsNode.origin();
1847 
1848  const LeafNodeType* refSignNode = refSignFlagsTreeAcc.probeConstLeaf(origin);
1849  if (!refSignNode) continue;
1850 
1851  const auto* const rsData = refSignNode->buffer().data();
1852  auto* const sfData = signFlagsNode.buffer().data();
1853 
1854  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1855  const Index offset = it.pos();
1856 
1857  const uint8_t rhsSigns = uint8_t(rsData[offset] & SIGNS);
1858 
1859  if (sEdgeGroupTable[rhsSigns][0] > 0) {
1860 
1861  typename LeafNodeType::ValueType& value = sfData[offset];
1862  const uint8_t lhsSigns = uint8_t(value & SIGNS);
1863 
1864  if (rhsSigns != lhsSigns) {
1865  value |= SEAM;
1866  }
1867  }
1868 
1869  } // end value on loop
1870 
1871  } // end leaf node loop
1872  }
1873 
1874 private:
1875  LeafNodeType * const * const mSignFlagsNodes;
1876  TreeType const * const mRefSignFlagsTree;
1877 }; // struct SetSeamLineFlags
1878 
1879 
1880 template <typename BoolTreeType, typename SignDataType>
1881 struct TransferSeamLineFlags
1882 {
1883  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1884 
1885  using SignDataTreeType = typename BoolTreeType::template ValueConverter<SignDataType>::Type;
1886  using SignDataLeafNodeType = typename SignDataTreeType::LeafNodeType;
1887 
1888  TransferSeamLineFlags(const std::vector<SignDataLeafNodeType*>& signFlagsLeafNodes,
1889  const BoolTreeType& maskTree)
1890  : mSignFlagsNodes(signFlagsLeafNodes.data())
1891  , mMaskTree(&maskTree)
1892  {
1893  }
1894 
1895  void operator()(const tbb::blocked_range<size_t>& range) const
1896  {
1897  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1898 
1899  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1900 
1901  SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1902  const Coord& origin = signFlagsNode.origin();
1903 
1904  const BoolLeafNodeType * maskNode = maskAcc.probeConstLeaf(origin);
1905  if (!maskNode) continue;
1906 
1907  auto* const sfData = signFlagsNode.buffer().data();
1908 
1909  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1910  const Index offset = it.pos();
1911 
1912  if (maskNode->isValueOn(offset)) {
1913  sfData[offset] |= SEAM;
1914  }
1915  } // end value on loop
1916  } // end leaf node loop
1917  }
1918 
1919 private:
1920  SignDataLeafNodeType * const * const mSignFlagsNodes;
1921  BoolTreeType const * const mMaskTree;
1922 }; // struct TransferSeamLineFlags
1923 
1924 
1925 template <typename TreeType>
1926 struct MaskSeamLineVoxels
1927 {
1928  using LeafNodeType = typename TreeType::LeafNodeType;
1929  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1930 
1931  MaskSeamLineVoxels(const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1932  const TreeType& signFlagsTree,
1933  BoolTreeType& mask)
1934  : mSignFlagsNodes(signFlagsLeafNodes.data())
1935  , mSignFlagsTree(&signFlagsTree)
1936  , mTempMask(false)
1937  , mMask(&mask)
1938  {
1939  }
1940 
1941  MaskSeamLineVoxels(MaskSeamLineVoxels& rhs, tbb::split)
1942  : mSignFlagsNodes(rhs.mSignFlagsNodes)
1943  , mSignFlagsTree(rhs.mSignFlagsTree)
1944  , mTempMask(false)
1945  , mMask(&mTempMask)
1946  {
1947  }
1948 
1949  void join(MaskSeamLineVoxels& rhs) { mMask->merge(*rhs.mMask); }
1950 
1951  void operator()(const tbb::blocked_range<size_t>& range)
1952  {
1953  using ValueOnCIter = typename LeafNodeType::ValueOnCIter;
1954  using ValueType = typename LeafNodeType::ValueType;
1955 
1956  tree::ValueAccessor<const TreeType> signFlagsAcc(*mSignFlagsTree);
1957  tree::ValueAccessor<BoolTreeType> maskAcc(*mMask);
1958  Coord ijk;
1959 
1960  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1961 
1962  LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1963  auto* const sfData = signFlagsNode.buffer().data();
1964 
1965  for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1966 
1967  const ValueType flags = sfData[it.pos()];
1968 
1969  if (!(flags & SEAM) && (flags & EDGES)) {
1970 
1971  ijk = it.getCoord();
1972 
1973  bool isSeamLineVoxel = false;
1974 
1975  if (flags & XEDGE) {
1976  ijk[1] -= 1;
1977  isSeamLineVoxel = (signFlagsAcc.getValue(ijk) & SEAM);
1978  ijk[2] -= 1;
1979  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1980  ijk[1] += 1;
1981  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1982  ijk[2] += 1;
1983  }
1984 
1985  if (!isSeamLineVoxel && flags & YEDGE) {
1986  ijk[2] -= 1;
1987  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1988  ijk[0] -= 1;
1989  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1990  ijk[2] += 1;
1991  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1992  ijk[0] += 1;
1993  }
1994 
1995  if (!isSeamLineVoxel && flags & ZEDGE) {
1996  ijk[1] -= 1;
1997  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1998  ijk[0] -= 1;
1999  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2000  ijk[1] += 1;
2001  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2002  ijk[0] += 1;
2003  }
2004 
2005  if (isSeamLineVoxel) {
2006  maskAcc.setValue(it.getCoord(), true);
2007  }
2008  }
2009  } // end value on loop
2010 
2011  } // end leaf node loop
2012  }
2013 
2014 private:
2015  LeafNodeType * const * const mSignFlagsNodes;
2016  TreeType const * const mSignFlagsTree;
2017  BoolTreeType mTempMask;
2018  BoolTreeType * const mMask;
2019 }; // struct MaskSeamLineVoxels
2020 
2021 
2022 template<typename SignDataTreeType>
2023 inline void
2024 markSeamLineData(SignDataTreeType& signFlagsTree, const SignDataTreeType& refSignFlagsTree)
2025 {
2026  using SignDataType = typename SignDataTreeType::ValueType;
2027  using SignDataLeafNodeType = typename SignDataTreeType::LeafNodeType;
2028  using BoolTreeType = typename SignDataTreeType::template ValueConverter<bool>::Type;
2029 
2030  std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2031  signFlagsTree.getNodes(signFlagsLeafNodes);
2032 
2033  const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2034 
2035  tbb::parallel_for(nodeRange,
2036  SetSeamLineFlags<SignDataTreeType>(signFlagsLeafNodes, refSignFlagsTree));
2037 
2038  BoolTreeType seamLineMaskTree(false);
2039 
2040  MaskSeamLineVoxels<SignDataTreeType>
2041  maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2042 
2043  tbb::parallel_reduce(nodeRange, maskSeamLine);
2044 
2045  tbb::parallel_for(nodeRange,
2046  TransferSeamLineFlags<BoolTreeType, SignDataType>(signFlagsLeafNodes, seamLineMaskTree));
2047 }
2048 
2049 
2050 ////////////////////////////////////////
2051 
2052 
2053 template <typename InputGridType>
2054 struct MergeVoxelRegions
2055 {
2056  using InputTreeType = typename InputGridType::TreeType;
2057  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
2058  using InputValueType = typename InputLeafNodeType::ValueType;
2059 
2060  using FloatTreeType = typename InputTreeType::template ValueConverter<float>::Type;
2061  using FloatLeafNodeType = typename FloatTreeType::LeafNodeType;
2062  using FloatGridType = Grid<FloatTreeType>;
2063 
2064  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
2065  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
2066 
2067  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
2068  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
2069 
2070  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
2071  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2072 
2073  using MaskTreeType = typename InputTreeType::template ValueConverter<ValueMask>::Type;
2074  using MaskLeafNodeType = typename MaskTreeType::LeafNodeType;
2075 
2076  MergeVoxelRegions(const InputGridType& inputGrid,
2077  const Index32TreeType& pointIndexTree,
2078  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2079  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2080  InputValueType iso,
2081  float adaptivity,
2082  bool invertSurfaceOrientation);
2083 
2084  void setSpatialAdaptivity(const FloatGridType& grid)
2085  {
2086  mSpatialAdaptivityTree = &grid.tree();
2087  mSpatialAdaptivityTransform = &grid.transform();
2088  }
2089 
2090  void setAdaptivityMask(const BoolTreeType& mask)
2091  {
2092  mMaskTree = &mask;
2093  }
2094 
2095  void setRefSignFlagsData(const Int16TreeType& signFlagsData, float internalAdaptivity)
2096  {
2097  mRefSignFlagsTree = &signFlagsData;
2098  mInternalAdaptivity = internalAdaptivity;
2099  }
2100 
2101  void operator()(const tbb::blocked_range<size_t>&) const;
2102 
2103 private:
2104  InputTreeType const * const mInputTree;
2105  math::Transform const * const mInputTransform;
2106 
2107  Index32TreeType const * const mPointIndexTree;
2108  Index32LeafNodeType * const * const mPointIndexNodes;
2109  Int16LeafNodeType const * const * const mSignFlagsNodes;
2110 
2111  InputValueType mIsovalue;
2112  float mSurfaceAdaptivity, mInternalAdaptivity;
2113  bool mInvertSurfaceOrientation;
2114 
2115  FloatTreeType const * mSpatialAdaptivityTree;
2116  BoolTreeType const * mMaskTree;
2117  Int16TreeType const * mRefSignFlagsTree;
2118  math::Transform const * mSpatialAdaptivityTransform;
2119 }; // struct MergeVoxelRegions
2120 
2121 
2122 template <typename InputGridType>
2123 MergeVoxelRegions<InputGridType>::MergeVoxelRegions(
2124  const InputGridType& inputGrid,
2125  const Index32TreeType& pointIndexTree,
2126  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2127  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2128  InputValueType iso,
2129  float adaptivity,
2130  bool invertSurfaceOrientation)
2131  : mInputTree(&inputGrid.tree())
2132  , mInputTransform(&inputGrid.transform())
2133  , mPointIndexTree(&pointIndexTree)
2134  , mPointIndexNodes(pointIndexLeafNodes.data())
2135  , mSignFlagsNodes(signFlagsLeafNodes.data())
2136  , mIsovalue(iso)
2137  , mSurfaceAdaptivity(adaptivity)
2138  , mInternalAdaptivity(adaptivity)
2139  , mInvertSurfaceOrientation(invertSurfaceOrientation)
2140  , mSpatialAdaptivityTree(nullptr)
2141  , mMaskTree(nullptr)
2142  , mRefSignFlagsTree(nullptr)
2143  , mSpatialAdaptivityTransform(nullptr)
2144 {
2145 }
2146 
2147 
2148 template <typename InputGridType>
2149 void
2150 MergeVoxelRegions<InputGridType>::operator()(const tbb::blocked_range<size_t>& range) const
2151 {
2152  using Vec3sType = math::Vec3<float>;
2153  using Vec3sLeafNodeType = typename InputLeafNodeType::template ValueConverter<Vec3sType>::Type;
2154 
2155  using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
2156  using FloatTreeAccessor = tree::ValueAccessor<const FloatTreeType>;
2157  using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
2158  using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
2159  using BoolTreeAccessor = tree::ValueAccessor<const BoolTreeType>;
2160 
2161  std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2162  if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2163  spatialAdaptivityAcc.reset(new FloatTreeAccessor(*mSpatialAdaptivityTree));
2164  }
2165 
2166  std::unique_ptr<BoolTreeAccessor> maskAcc;
2167  if (mMaskTree) {
2168  maskAcc.reset(new BoolTreeAccessor(*mMaskTree));
2169  }
2170 
2171  std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2172  if (mRefSignFlagsTree) {
2173  refSignFlagsAcc.reset(new Int16TreeAccessor(*mRefSignFlagsTree));
2174  }
2175 
2176  InputTreeAccessor inputAcc(*mInputTree);
2177  Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2178 
2179  MaskLeafNodeType mask;
2180 
2181  const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2182  std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2183 
2184  Coord ijk, end;
2185  const int LeafDim = InputLeafNodeType::DIM;
2186 
2187  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
2188 
2189  mask.setValuesOff();
2190 
2191  const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
2192  Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[n];
2193 
2194  const Coord& origin = pointIndexNode.origin();
2195  end = origin.offsetBy(LeafDim);
2196 
2197  // Mask off seam line adjacent voxels
2198  if (maskAcc) {
2199  const BoolLeafNodeType* maskLeaf = maskAcc->probeConstLeaf(origin);
2200  if (maskLeaf != nullptr) {
2201  for (auto it = maskLeaf->cbeginValueOn(); it; ++it)
2202  {
2203  mask.setActiveState(it.getCoord() & ~1u, true);
2204  }
2205  }
2206  }
2207 
2208  float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2209  mInternalAdaptivity : mSurfaceAdaptivity;
2210 
2211  bool useGradients = adaptivity < 1.0f;
2212 
2213  // Set region adaptivity
2214  FloatLeafNodeType adaptivityLeaf(origin, adaptivity);
2215 
2216  if (spatialAdaptivityAcc) {
2217  useGradients = false;
2218  for (Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++offset) {
2219  ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2220  ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2221  mInputTransform->indexToWorld(ijk));
2222  float weight = spatialAdaptivityAcc->getValue(ijk);
2223  float adaptivityValue = weight * adaptivity;
2224  if (adaptivityValue < 1.0f) useGradients = true;
2225  adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2226  }
2227  }
2228 
2229  // Mask off ambiguous voxels
2230  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2231  const Int16 flags = it.getValue();
2232  const unsigned char signs = static_cast<unsigned char>(SIGNS & int(flags));
2233 
2234  if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2235 
2236  mask.setActiveState(it.getCoord() & ~1u, true);
2237 
2238  } else if (flags & EDGES) {
2239 
2240  bool maskRegion = false;
2241 
2242  ijk = it.getCoord();
2243  if (!pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2244 
2245  if (!maskRegion && flags & XEDGE) {
2246  ijk[1] -= 1;
2247  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2248  ijk[2] -= 1;
2249  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2250  ijk[1] += 1;
2251  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2252  ijk[2] += 1;
2253  }
2254 
2255  if (!maskRegion && flags & YEDGE) {
2256  ijk[2] -= 1;
2257  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2258  ijk[0] -= 1;
2259  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2260  ijk[2] += 1;
2261  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2262  ijk[0] += 1;
2263  }
2264 
2265  if (!maskRegion && flags & ZEDGE) {
2266  ijk[1] -= 1;
2267  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2268  ijk[0] -= 1;
2269  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2270  ijk[1] += 1;
2271  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2272  ijk[0] += 1;
2273  }
2274 
2275  if (maskRegion) {
2276  mask.setActiveState(it.getCoord() & ~1u, true);
2277  }
2278  }
2279  }
2280 
2281  // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
2282  int dim = 2;
2283  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2284  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2285  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2286  if (!mask.isValueOn(ijk) && isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2287  mask.setActiveState(ijk, true);
2288  }
2289  }
2290  }
2291  }
2292 
2293  // Compute the gradient for the remaining voxels
2294 
2295  if (useGradients) {
2296 
2297  if (gradientNode) {
2298  gradientNode->setValuesOff();
2299  } else {
2300  gradientNode.reset(new Vec3sLeafNodeType());
2301  }
2302 
2303  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it)
2304  {
2305  ijk = it.getCoord();
2306  if (!mask.isValueOn(ijk & ~1u))
2307  {
2308  Vec3sType dir(math::ISGradient<math::CD_2ND>::result(inputAcc, ijk));
2309  dir.normalize();
2310 
2311  if (invertGradientDir) {
2312  dir = -dir;
2313  }
2314 
2315  gradientNode->setValueOn(it.pos(), dir);
2316  }
2317  }
2318  }
2319 
2320  // Merge regions
2321  int regionId = 1;
2322  for ( ; dim <= LeafDim; dim = dim << 1) {
2323  const unsigned coordMask = ~((dim << 1) - 1);
2324  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2325  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2326  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2327 
2328  adaptivity = adaptivityLeaf.getValue(ijk);
2329 
2330  if (mask.isValueOn(ijk)
2331  || isNonManifold(inputAcc, ijk, mIsovalue, dim)
2332  || (useGradients && !isMergeable(*gradientNode, ijk, dim, adaptivity)))
2333  {
2334  mask.setActiveState(ijk & coordMask, true);
2335  } else {
2336  mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2337  }
2338  }
2339  }
2340  }
2341  }
2342  }
2343 } // MergeVoxelRegions::operator()
2344 
2345 
2346 ////////////////////////////////////////
2347 
2348 
2349 // Constructs qudas
2350 struct UniformPrimBuilder
2351 {
2352  UniformPrimBuilder(): mIdx(0), mPolygonPool(nullptr) {}
2353 
2354  void init(const size_t upperBound, PolygonPool& quadPool)
2355  {
2356  mPolygonPool = &quadPool;
2357  mPolygonPool->resetQuads(upperBound);
2358  mIdx = 0;
2359  }
2360 
2361  template<typename IndexType>
2362  void addPrim(const math::Vec4<IndexType>& verts, bool reverse, char flags = 0)
2363  {
2364  if (!reverse) {
2365  mPolygonPool->quad(mIdx) = verts;
2366  } else {
2367  Vec4I& quad = mPolygonPool->quad(mIdx);
2368  quad[0] = verts[3];
2369  quad[1] = verts[2];
2370  quad[2] = verts[1];
2371  quad[3] = verts[0];
2372  }
2373  mPolygonPool->quadFlags(mIdx) = flags;
2374  ++mIdx;
2375  }
2376 
2377  void done()
2378  {
2379  mPolygonPool->trimQuads(mIdx);
2380  }
2381 
2382 private:
2383  size_t mIdx;
2384  PolygonPool* mPolygonPool;
2385 };
2386 
2387 
2388 // Constructs qudas and triangles
2389 struct AdaptivePrimBuilder
2390 {
2391  AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(nullptr) {}
2392 
2393  void init(const size_t upperBound, PolygonPool& polygonPool)
2394  {
2395  mPolygonPool = &polygonPool;
2396  mPolygonPool->resetQuads(upperBound);
2397  mPolygonPool->resetTriangles(upperBound);
2398 
2399  mQuadIdx = 0;
2400  mTriangleIdx = 0;
2401  }
2402 
2403  template<typename IndexType>
2404  void addPrim(const math::Vec4<IndexType>& verts, bool reverse, char flags = 0)
2405  {
2406  if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2407  && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2408  mPolygonPool->quadFlags(mQuadIdx) = flags;
2409  addQuad(verts, reverse);
2410  } else if (
2411  verts[0] == verts[3] &&
2412  verts[1] != verts[2] &&
2413  verts[1] != verts[0] &&
2414  verts[2] != verts[0]) {
2415  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2416  addTriangle(verts[0], verts[1], verts[2], reverse);
2417  } else if (
2418  verts[1] == verts[2] &&
2419  verts[0] != verts[3] &&
2420  verts[0] != verts[1] &&
2421  verts[3] != verts[1]) {
2422  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2423  addTriangle(verts[0], verts[1], verts[3], reverse);
2424  } else if (
2425  verts[0] == verts[1] &&
2426  verts[2] != verts[3] &&
2427  verts[2] != verts[0] &&
2428  verts[3] != verts[0]) {
2429  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2430  addTriangle(verts[0], verts[2], verts[3], reverse);
2431  } else if (
2432  verts[2] == verts[3] &&
2433  verts[0] != verts[1] &&
2434  verts[0] != verts[2] &&
2435  verts[1] != verts[2]) {
2436  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2437  addTriangle(verts[0], verts[1], verts[2], reverse);
2438  }
2439  }
2440 
2441 
2442  void done()
2443  {
2444  mPolygonPool->trimQuads(mQuadIdx, /*reallocate=*/true);
2445  mPolygonPool->trimTrinagles(mTriangleIdx, /*reallocate=*/true);
2446  }
2447 
2448 private:
2449 
2450  template<typename IndexType>
2451  void addQuad(const math::Vec4<IndexType>& verts, bool reverse)
2452  {
2453  if (!reverse) {
2454  mPolygonPool->quad(mQuadIdx) = verts;
2455  } else {
2456  Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2457  quad[0] = verts[3];
2458  quad[1] = verts[2];
2459  quad[2] = verts[1];
2460  quad[3] = verts[0];
2461  }
2462  ++mQuadIdx;
2463  }
2464 
2465  void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
2466  {
2467  Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2468 
2469  prim[1] = v1;
2470 
2471  if (!reverse) {
2472  prim[0] = v0;
2473  prim[2] = v2;
2474  } else {
2475  prim[0] = v2;
2476  prim[2] = v0;
2477  }
2478  ++mTriangleIdx;
2479  }
2480 
2481  size_t mQuadIdx, mTriangleIdx;
2482  PolygonPool *mPolygonPool;
2483 };
2484 
2485 
2486 template<typename SignAccT, typename IdxAccT, typename PrimBuilder>
2487 inline void
2488 constructPolygons(
2489  bool invertSurfaceOrientation,
2490  Int16 flags,
2491  Int16 refFlags,
2492  const Vec3i& offsets,
2493  const Coord& ijk,
2494  const SignAccT& signAcc,
2495  const IdxAccT& idxAcc,
2496  PrimBuilder& mesher)
2497 {
2498  using IndexType = typename IdxAccT::ValueType;
2499 
2500  IndexType v0 = IndexType(util::INVALID_IDX);
2501  const bool isActive = idxAcc.probeValue(ijk, v0);
2502  if (isActive == false || v0 == IndexType(util::INVALID_IDX)) return;
2503 
2504  char tag[2];
2505  tag[0] = (flags & SEAM) ? POLYFLAG_FRACTURE_SEAM : 0;
2506  tag[1] = tag[0] | char(POLYFLAG_EXTERIOR);
2507 
2508  bool isInside = flags & INSIDE;
2509 
2510  isInside = invertSurfaceOrientation ? !isInside : isInside;
2511 
2512  Coord coord = ijk;
2513  math::Vec4<IndexType> quad(0,0,0,0);
2514 
2515  if (flags & XEDGE) {
2516 
2517  quad[0] = v0 + offsets[0];
2518 
2519  // i, j-1, k
2520  coord[1]--;
2521  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2522  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2523  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2524 
2525  // i, j-1, k-1
2526  coord[2]--;
2527  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2528  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2529  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2530 
2531  // i, j, k-1
2532  coord[1]++;
2533  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2534  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2535  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2536 
2537  if (activeValues) {
2538  mesher.addPrim(quad, isInside, tag[bool(refFlags & XEDGE)]);
2539  }
2540 
2541  coord[2]++; // i, j, k
2542  }
2543 
2544 
2545  if (flags & YEDGE) {
2546 
2547  quad[0] = v0 + offsets[1];
2548 
2549  // i, j, k-1
2550  coord[2]--;
2551  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2552  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2553  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2554 
2555  // i-1, j, k-1
2556  coord[0]--;
2557  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2558  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2559  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2560 
2561  // i-1, j, k
2562  coord[2]++;
2563  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2564  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2565  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2566 
2567  if (activeValues) {
2568  mesher.addPrim(quad, isInside, tag[bool(refFlags & YEDGE)]);
2569  }
2570 
2571  coord[0]++; // i, j, k
2572  }
2573 
2574 
2575  if (flags & ZEDGE) {
2576 
2577  quad[0] = v0 + offsets[2];
2578 
2579  // i, j-1, k
2580  coord[1]--;
2581  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2582  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2583  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2584 
2585  // i-1, j-1, k
2586  coord[0]--;
2587  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2588  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2589  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2590 
2591  // i-1, j, k
2592  coord[1]++;
2593  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2594  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2595  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2596 
2597  if (activeValues) {
2598  mesher.addPrim(quad, !isInside, tag[bool(refFlags & ZEDGE)]);
2599  }
2600  }
2601 }
2602 
2603 
2604 ////////////////////////////////////////
2605 
2606 
2607 template<typename InputTreeType>
2608 struct MaskTileBorders
2609 {
2610  using InputValueType = typename InputTreeType::ValueType;
2611  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
2612 
2613 
2614  MaskTileBorders(const InputTreeType& inputTree, InputValueType iso,
2615  BoolTreeType& mask, const Vec4i* tileArray)
2616  : mInputTree(&inputTree)
2617  , mIsovalue(iso)
2618  , mTempMask(false)
2619  , mMask(&mask)
2620  , mTileArray(tileArray)
2621  {
2622  }
2623 
2624  MaskTileBorders(MaskTileBorders& rhs, tbb::split)
2625  : mInputTree(rhs.mInputTree)
2626  , mIsovalue(rhs.mIsovalue)
2627  , mTempMask(false)
2628  , mMask(&mTempMask)
2629  , mTileArray(rhs.mTileArray)
2630  {
2631  }
2632 
2633  void join(MaskTileBorders& rhs) { mMask->merge(*rhs.mMask); }
2634 
2635  void operator()(const tbb::blocked_range<size_t>&);
2636 
2637 private:
2638  InputTreeType const * const mInputTree;
2639  InputValueType const mIsovalue;
2640  BoolTreeType mTempMask;
2641  BoolTreeType * const mMask;
2642  Vec4i const * const mTileArray;
2643 }; // MaskTileBorders
2644 
2645 
2646 template<typename InputTreeType>
2647 void
2648 MaskTileBorders<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
2649 {
2650  tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
2651 
2652  CoordBBox region, bbox;
2653  Coord ijk, nijk;
2654 
2655  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
2656 
2657  const Vec4i& tile = mTileArray[n];
2658 
2659  bbox.min()[0] = tile[0];
2660  bbox.min()[1] = tile[1];
2661  bbox.min()[2] = tile[2];
2662  bbox.max() = bbox.min();
2663  bbox.max().offset(tile[3]);
2664 
2665  InputValueType value = mInputTree->background();
2666 
2667  const bool isInside = isInsideValue(inputTreeAcc.getValue(bbox.min()), mIsovalue);
2668  const int valueDepth = inputTreeAcc.getValueDepth(bbox.min());
2669 
2670  // eval x-edges
2671 
2672  ijk = bbox.max();
2673  nijk = ijk;
2674  ++nijk[0];
2675 
2676  bool processRegion = true;
2677  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2678  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2679  }
2680 
2681  if (processRegion) {
2682  region = bbox;
2683  region.expand(1);
2684  region.min()[0] = region.max()[0] = ijk[0];
2685  mMask->fill(region, false);
2686  }
2687 
2688 
2689  ijk = bbox.min();
2690  --ijk[0];
2691 
2692  processRegion = true;
2693  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2694  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2695  && isInside != isInsideValue(value, mIsovalue));
2696  }
2697 
2698  if (processRegion) {
2699  region = bbox;
2700  region.expand(1);
2701  region.min()[0] = region.max()[0] = ijk[0];
2702  mMask->fill(region, false);
2703  }
2704 
2705 
2706  // eval y-edges
2707 
2708  ijk = bbox.max();
2709  nijk = ijk;
2710  ++nijk[1];
2711 
2712  processRegion = true;
2713  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2714  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2715  }
2716 
2717  if (processRegion) {
2718  region = bbox;
2719  region.expand(1);
2720  region.min()[1] = region.max()[1] = ijk[1];
2721  mMask->fill(region, false);
2722  }
2723 
2724 
2725  ijk = bbox.min();
2726  --ijk[1];
2727 
2728  processRegion = true;
2729  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2730  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2731  && isInside != isInsideValue(value, mIsovalue));
2732  }
2733 
2734  if (processRegion) {
2735  region = bbox;
2736  region.expand(1);
2737  region.min()[1] = region.max()[1] = ijk[1];
2738  mMask->fill(region, false);
2739  }
2740 
2741 
2742  // eval z-edges
2743 
2744  ijk = bbox.max();
2745  nijk = ijk;
2746  ++nijk[2];
2747 
2748  processRegion = true;
2749  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2750  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2751  }
2752 
2753  if (processRegion) {
2754  region = bbox;
2755  region.expand(1);
2756  region.min()[2] = region.max()[2] = ijk[2];
2757  mMask->fill(region, false);
2758  }
2759 
2760  ijk = bbox.min();
2761  --ijk[2];
2762 
2763  processRegion = true;
2764  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2765  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2766  && isInside != isInsideValue(value, mIsovalue));
2767  }
2768 
2769  if (processRegion) {
2770  region = bbox;
2771  region.expand(1);
2772  region.min()[2] = region.max()[2] = ijk[2];
2773  mMask->fill(region, false);
2774  }
2775  }
2776 } // MaskTileBorders::operator()
2777 
2778 
2779 template<typename InputTreeType>
2780 inline void
2781 maskActiveTileBorders(const InputTreeType& inputTree,
2782  const typename InputTreeType::ValueType iso,
2783  typename InputTreeType::template ValueConverter<bool>::Type& mask)
2784 {
2785  typename InputTreeType::ValueOnCIter tileIter(inputTree);
2786  tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2787 
2788  size_t tileCount = 0;
2789  for ( ; tileIter; ++tileIter) {
2790  ++tileCount;
2791  }
2792 
2793  if (tileCount > 0) {
2794  std::unique_ptr<Vec4i[]> tiles(new Vec4i[tileCount]);
2795 
2796  CoordBBox bbox;
2797  size_t index = 0;
2798 
2799  tileIter = inputTree.cbeginValueOn();
2800  tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2801 
2802  for (; tileIter; ++tileIter) {
2803  Vec4i& tile = tiles[index++];
2804  tileIter.getBoundingBox(bbox);
2805  tile[0] = bbox.min()[0];
2806  tile[1] = bbox.min()[1];
2807  tile[2] = bbox.min()[2];
2808  tile[3] = bbox.max()[0] - bbox.min()[0];
2809  }
2810 
2811  MaskTileBorders<InputTreeType> op(inputTree, iso, mask, tiles.get());
2812  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2813  }
2814 }
2815 
2816 
2817 ////////////////////////////////////////
2818 
2819 
2820 // Utility class for the volumeToMesh wrapper
2821 class PointListCopy
2822 {
2823 public:
2824  PointListCopy(const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2825  : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2826  {
2827  }
2828 
2829  void operator()(const tbb::blocked_range<size_t>& range) const
2830  {
2831  for (size_t n = range.begin(); n < range.end(); ++n) {
2832  mPointsOut[n] = mPointsIn[n];
2833  }
2834  }
2835 
2836 private:
2837  const PointList& mPointsIn;
2838  std::vector<Vec3s>& mPointsOut;
2839 };
2840 
2841 
2842 ////////////////////////////////////////////////////////////////////////////////
2843 ////////////////////////////////////////////////////////////////////////////////
2844 
2845 
2846 struct LeafNodeVoxelOffsets
2847 {
2848  using IndexVector = std::vector<Index>;
2849 
2850  template<typename LeafNodeType>
2851  void constructOffsetList();
2852 
2853  /// Return internal core voxel offsets.
2854  const IndexVector& core() const { return mCore; }
2855 
2856 
2857  /// Return front face voxel offsets.
2858  const IndexVector& minX() const { return mMinX; }
2859 
2860  /// Return back face voxel offsets.
2861  const IndexVector& maxX() const { return mMaxX; }
2862 
2863 
2864  /// Return bottom face voxel offsets.
2865  const IndexVector& minY() const { return mMinY; }
2866 
2867  /// Return top face voxel offsets.
2868  const IndexVector& maxY() const { return mMaxY; }
2869 
2870 
2871  /// Return left face voxel offsets.
2872  const IndexVector& minZ() const { return mMinZ; }
2873 
2874  /// Return right face voxel offsets.
2875  const IndexVector& maxZ() const { return mMaxZ; }
2876 
2877 
2878  /// Return voxel offsets with internal neighbours in x + 1.
2879  const IndexVector& internalNeighborsX() const { return mInternalNeighborsX; }
2880 
2881  /// Return voxel offsets with internal neighbours in y + 1.
2882  const IndexVector& internalNeighborsY() const { return mInternalNeighborsY; }
2883 
2884  /// Return voxel offsets with internal neighbours in z + 1.
2885  const IndexVector& internalNeighborsZ() const { return mInternalNeighborsZ; }
2886 
2887 
2888 private:
2889  IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2890  mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2891 }; // struct LeafNodeOffsets
2892 
2893 
2894 template<typename LeafNodeType>
2895 inline void
2896 LeafNodeVoxelOffsets::constructOffsetList()
2897 {
2898  // internal core voxels
2899  mCore.clear();
2900  mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2901 
2902  for (Index x = 1; x < (LeafNodeType::DIM - 1); ++x) {
2903  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2904  for (Index y = 1; y < (LeafNodeType::DIM - 1); ++y) {
2905  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2906  for (Index z = 1; z < (LeafNodeType::DIM - 1); ++z) {
2907  mCore.push_back(offsetXY + z);
2908  }
2909  }
2910  }
2911 
2912  // internal neighbors in x + 1
2913  mInternalNeighborsX.clear();
2914  mInternalNeighborsX.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2915 
2916  for (Index x = 0; x < (LeafNodeType::DIM - 1); ++x) {
2917  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2918  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2919  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2920  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2921  mInternalNeighborsX.push_back(offsetXY + z);
2922  }
2923  }
2924  }
2925 
2926  // internal neighbors in y + 1
2927  mInternalNeighborsY.clear();
2928  mInternalNeighborsY.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2929 
2930  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2931  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2932  for (Index y = 0; y < (LeafNodeType::DIM - 1); ++y) {
2933  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2934  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2935  mInternalNeighborsY.push_back(offsetXY + z);
2936  }
2937  }
2938  }
2939 
2940  // internal neighbors in z + 1
2941  mInternalNeighborsZ.clear();
2942  mInternalNeighborsZ.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2943 
2944  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2945  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2946  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2947  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2948  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2949  mInternalNeighborsZ.push_back(offsetXY + z);
2950  }
2951  }
2952  }
2953 
2954  // min x
2955  mMinX.clear();
2956  mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2957  {
2958  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2959  const Index offsetXY = (y << LeafNodeType::LOG2DIM);
2960  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2961  mMinX.push_back(offsetXY + z);
2962  }
2963  }
2964  }
2965 
2966  // max x
2967  mMaxX.clear();
2968  mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2969  {
2970  const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
2971  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2972  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2973  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2974  mMaxX.push_back(offsetXY + z);
2975  }
2976  }
2977  }
2978 
2979  // min y
2980  mMinY.clear();
2981  mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2982  {
2983  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2984  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2985  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2986  mMinY.push_back(offsetX + z);
2987  }
2988  }
2989  }
2990 
2991  // max y
2992  mMaxY.clear();
2993  mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2994  {
2995  const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
2996  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2997  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2998  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2999  mMaxY.push_back(offsetX + offsetY + z);
3000  }
3001  }
3002  }
3003 
3004  // min z
3005  mMinZ.clear();
3006  mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3007  {
3008  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
3009  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3010  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
3011  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3012  mMinZ.push_back(offsetXY);
3013  }
3014  }
3015  }
3016 
3017  // max z
3018  mMaxZ.clear();
3019  mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3020  {
3021  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
3022  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3023  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
3024  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3025  mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3026  }
3027  }
3028  }
3029 }
3030 
3031 
3032 ////////////////////////////////////////
3033 
3034 
3035 /// Utility method to marks all voxels that share an edge.
3036 template<typename AccessorT, int _AXIS>
3037 struct VoxelEdgeAccessor
3038 {
3039  enum { AXIS = _AXIS };
3040  AccessorT& acc;
3041 
3042  VoxelEdgeAccessor(AccessorT& _acc) : acc(_acc) {}
3043 
3044  void set(Coord ijk) {
3045  if (_AXIS == 0) { // x + 1 edge
3046  acc.setActiveState(ijk);
3047  --ijk[1]; // set i, j-1, k
3048  acc.setActiveState(ijk);
3049  --ijk[2]; // set i, j-1, k-1
3050  acc.setActiveState(ijk);
3051  ++ijk[1]; // set i, j, k-1
3052  acc.setActiveState(ijk);
3053  } else if (_AXIS == 1) { // y + 1 edge
3054  acc.setActiveState(ijk);
3055  --ijk[2]; // set i, j, k-1
3056  acc.setActiveState(ijk);
3057  --ijk[0]; // set i-1, j, k-1
3058  acc.setActiveState(ijk);
3059  ++ijk[2]; // set i-1, j, k
3060  acc.setActiveState(ijk);
3061  } else { // z + 1 edge
3062  acc.setActiveState(ijk);
3063  --ijk[1]; // set i, j-1, k
3064  acc.setActiveState(ijk);
3065  --ijk[0]; // set i-1, j-1, k
3066  acc.setActiveState(ijk);
3067  ++ijk[1]; // set i-1, j, k
3068  acc.setActiveState(ijk);
3069  }
3070  }
3071 };
3072 
3073 
3074 /// Utility method to check for sign changes along the x + 1, y + 1 or z + 1 directions.
3075 /// The direction is determined by the @a edgeAcc parameter. Only voxels that have internal
3076 /// neighbours are evaluated.
3077 template<typename VoxelEdgeAcc, typename LeafNodeT>
3078 void
3079 evalInternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3080  const LeafNodeT& leafnode,
3081  const LeafNodeVoxelOffsets& voxels,
3082  const typename LeafNodeT::ValueType iso)
3083 {
3084  Index nvo = 1; // neighbour voxel offset, z + 1 direction assumed initially.
3085  const std::vector<Index>* offsets = &voxels.internalNeighborsZ();
3086 
3087  if (VoxelEdgeAcc::AXIS == 0) { // x + 1 direction
3088  nvo = LeafNodeT::DIM * LeafNodeT::DIM;
3089  offsets = &voxels.internalNeighborsX();
3090  } else if (VoxelEdgeAcc::AXIS == 1) { // y + 1 direction
3091  nvo = LeafNodeT::DIM;
3092  offsets = &voxels.internalNeighborsY();
3093  }
3094 
3095  const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3096 
3097  for (size_t n = 0, N = offsets->size(); n < N; ++n) {
3098  const Index& pos = (*offsets)[n];
3099  const bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3100  if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3101  isInsideValue(lhsAcc.get((pos + nvo)), iso))) {
3102  edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3103  }
3104  }
3105 }
3106 
3107 
3108 /// Utility method to check for sign changes along the x + 1, y + 1 or z + 1 directions.
3109 /// The direction is determined by the @a edgeAcc parameter. All voxels that reside in the
3110 /// specified leafnode face: back, top or right are evaluated.
3111 template<typename LeafNodeT, typename TreeAcc, typename VoxelEdgeAcc>
3112 void
3113 evalExternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3114  const TreeAcc& acc,
3115  const LeafNodeT& lhsNode,
3116  const LeafNodeVoxelOffsets& voxels,
3117  const typename LeafNodeT::ValueType iso)
3118 {
3119  const std::vector<Index>* lhsOffsets = &voxels.maxX();
3120  const std::vector<Index>* rhsOffsets = &voxels.minX();
3121  Coord ijk = lhsNode.origin();
3122 
3123  if (VoxelEdgeAcc::AXIS == 0) { // back leafnode face
3124  ijk[0] += LeafNodeT::DIM;
3125  } else if (VoxelEdgeAcc::AXIS == 1) { // top leafnode face
3126  ijk[1] += LeafNodeT::DIM;
3127  lhsOffsets = &voxels.maxY();
3128  rhsOffsets = &voxels.minY();
3129  } else if (VoxelEdgeAcc::AXIS == 2) { // right leafnode face
3130  ijk[2] += LeafNodeT::DIM;
3131  lhsOffsets = &voxels.maxZ();
3132  rhsOffsets = &voxels.minZ();
3133  }
3134 
3135  typename LeafNodeT::ValueType value;
3136  const LeafNodeT* rhsNodePt = acc.probeConstLeaf(ijk);
3137 
3138  const LeafBufferAccessor<LeafNodeT> lhsAcc(lhsNode);
3139 
3140  if (rhsNodePt) {
3141  const LeafBufferAccessor<LeafNodeT> rhsAcc(*rhsNodePt);
3142 
3143  for (size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3144  const Index& pos = (*lhsOffsets)[n];
3145  bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[n]);
3146  if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3147  isInsideValue(rhsAcc.get((*rhsOffsets)[n]), iso))) {
3148  edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3149  }
3150  }
3151  } else if (!acc.probeValue(ijk, value)) {
3152  const bool inside = isInsideValue(value, iso);
3153  for (size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3154  const Index& pos = (*lhsOffsets)[n];
3155  if (lhsNode.isValueOn(pos) && (inside != isInsideValue(lhsAcc.get(pos), iso))) {
3156  edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3157  }
3158  }
3159  }
3160 }
3161 
3162 
3163 /// Utility method to check for sign changes along the x - 1, y - 1 or z - 1 directions.
3164 /// The direction is determined by the @a edgeAcc parameter. All voxels that reside in the
3165 /// specified leafnode face: front, bottom or left are evaluated.
3166 template<typename LeafNodeT, typename TreeAcc, typename VoxelEdgeAcc>
3167 void
3168 evalExternalVoxelEdgesInv(VoxelEdgeAcc& edgeAcc,
3169  const TreeAcc& acc,
3170  const LeafNodeT& leafnode,
3171  const LeafNodeVoxelOffsets& voxels,
3172  const typename LeafNodeT::ValueType iso)
3173 {
3174  Coord ijk = leafnode.origin();
3175  if (VoxelEdgeAcc::AXIS == 0) --ijk[0]; // front leafnode face
3176  else if (VoxelEdgeAcc::AXIS == 1) --ijk[1]; // bottom leafnode face
3177  else if (VoxelEdgeAcc::AXIS == 2) --ijk[2]; // left leafnode face
3178 
3179  typename LeafNodeT::ValueType value;
3180  if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3181 
3182  const std::vector<Index>* offsets = &voxels.internalNeighborsX();
3183  if (VoxelEdgeAcc::AXIS == 1) offsets = &voxels.internalNeighborsY();
3184  else if (VoxelEdgeAcc::AXIS == 2) offsets = &voxels.internalNeighborsZ();
3185 
3186  const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3187 
3188  const bool inside = isInsideValue(value, iso);
3189  for (size_t n = 0, N = offsets->size(); n < N; ++n) {
3190 
3191  const Index& pos = (*offsets)[n];
3192  if (leafnode.isValueOn(pos)
3193  && (inside != isInsideValue(lhsAcc.get(pos), iso)))
3194  {
3195  ijk = leafnode.offsetToGlobalCoord(pos);
3196  if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3197  else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3198  else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3199 
3200  edgeAcc.set(ijk);
3201  }
3202  }
3203  }
3204 }
3205 
3206 
3207 template<typename InputTreeType>
3208 struct IdentifyIntersectingVoxels
3209 {
3210  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3211  using InputValueType = typename InputLeafNodeType::ValueType;
3212 
3213  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3214 
3215  IdentifyIntersectingVoxels(
3216  const InputTreeType& inputTree,
3217  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3218  const LeafNodeVoxelOffsets& offsets,
3219  BoolTreeType& intersectionTree,
3220  InputValueType iso);
3221 
3222  IdentifyIntersectingVoxels(IdentifyIntersectingVoxels&, tbb::split);
3223  void operator()(const tbb::blocked_range<size_t>&);
3224  void join(const IdentifyIntersectingVoxels& rhs) {
3225  mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3226  }
3227 
3228 private:
3230  InputLeafNodeType const * const * const mInputNodes;
3231 
3232  BoolTreeType mIntersectionTree;
3233  tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3234 
3235  const LeafNodeVoxelOffsets& mOffsets;
3236  const InputValueType mIsovalue;
3237 }; // struct IdentifyIntersectingVoxels
3238 
3239 
3240 template<typename InputTreeType>
3241 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3242  const InputTreeType& inputTree,
3243  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3244  const LeafNodeVoxelOffsets& offsets,
3245  BoolTreeType& intersectionTree,
3246  InputValueType iso)
3247  : mInputAccessor(inputTree)
3248  , mInputNodes(inputLeafNodes.data())
3249  , mIntersectionTree(false)
3250  , mIntersectionAccessor(intersectionTree)
3251  , mOffsets(offsets)
3252  , mIsovalue(iso)
3253 {
3254 }
3255 
3256 
3257 template<typename InputTreeType>
3258 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3259  IdentifyIntersectingVoxels& rhs, tbb::split)
3260  : mInputAccessor(rhs.mInputAccessor.tree())
3261  , mInputNodes(rhs.mInputNodes)
3262  , mIntersectionTree(false)
3263  , mIntersectionAccessor(mIntersectionTree) // use local tree.
3264  , mOffsets(rhs.mOffsets)
3265  , mIsovalue(rhs.mIsovalue)
3266 {
3267 }
3268 
3269 
3270 template<typename InputTreeType>
3271 void
3272 IdentifyIntersectingVoxels<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3273 {
3274  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3275  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3276  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3277 
3278  for (size_t n = range.begin(); n != range.end(); ++n) {
3279 
3280  const InputLeafNodeType& node = *mInputNodes[n];
3281 
3282  // internal x + 1 voxel edges
3283  evalInternalVoxelEdges(xEdgeAcc, node, mOffsets, mIsovalue);
3284  // internal y + 1 voxel edges
3285  evalInternalVoxelEdges(yEdgeAcc, node, mOffsets, mIsovalue);
3286  // internal z + 1 voxel edges
3287  evalInternalVoxelEdges(zEdgeAcc, node, mOffsets, mIsovalue);
3288 
3289  // external x + 1 voxels edges (back face)
3290  evalExternalVoxelEdges(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3291  // external y + 1 voxels edges (top face)
3292  evalExternalVoxelEdges(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3293  // external z + 1 voxels edges (right face)
3294  evalExternalVoxelEdges(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3295 
3296  // The remaining edges are only checked if the leafnode neighbour, in the
3297  // corresponding direction, is an inactive tile.
3298 
3299  // external x - 1 voxels edges (front face)
3300  evalExternalVoxelEdgesInv(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3301  // external y - 1 voxels edges (bottom face)
3302  evalExternalVoxelEdgesInv(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3303  // external z - 1 voxels edges (left face)
3304  evalExternalVoxelEdgesInv(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3305  }
3306 } // IdentifyIntersectingVoxels::operator()
3307 
3308 
3309 template<typename InputTreeType>
3310 inline void
3311 identifySurfaceIntersectingVoxels(
3312  typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3313  const InputTreeType& inputTree,
3314  typename InputTreeType::ValueType isovalue)
3315 {
3316  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3317 
3318  std::vector<const InputLeafNodeType*> inputLeafNodes;
3319  inputTree.getNodes(inputLeafNodes);
3320 
3321  LeafNodeVoxelOffsets offsets;
3322  offsets.constructOffsetList<InputLeafNodeType>();
3323 
3324  IdentifyIntersectingVoxels<InputTreeType> op(
3325  inputTree, inputLeafNodes, offsets, intersectionTree, isovalue);
3326 
3327  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3328 
3329  maskActiveTileBorders(inputTree, isovalue, intersectionTree);
3330 }
3331 
3332 
3333 ////////////////////////////////////////
3334 
3335 
3336 template<typename InputTreeType>
3337 struct MaskIntersectingVoxels
3338 {
3339  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3340  using InputValueType = typename InputLeafNodeType::ValueType;
3341 
3342  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3343  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3344 
3345  MaskIntersectingVoxels(
3346  const InputTreeType& inputTree,
3347  const std::vector<BoolLeafNodeType*>& nodes,
3348  BoolTreeType& intersectionTree,
3349  InputValueType iso);
3350 
3351  MaskIntersectingVoxels(MaskIntersectingVoxels&, tbb::split);
3352  void operator()(const tbb::blocked_range<size_t>&);
3353  void join(const MaskIntersectingVoxels& rhs) {
3354  mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3355  }
3356 
3357 private:
3359  BoolLeafNodeType const * const * const mNodes;
3360 
3361  BoolTreeType mIntersectionTree;
3362  tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3363 
3364  const InputValueType mIsovalue;
3365 }; // struct MaskIntersectingVoxels
3366 
3367 
3368 template<typename InputTreeType>
3369 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3370  const InputTreeType& inputTree,
3371  const std::vector<BoolLeafNodeType*>& nodes,
3372  BoolTreeType& intersectionTree,
3373  InputValueType iso)
3374  : mInputAccessor(inputTree)
3375  , mNodes(nodes.data())
3376  , mIntersectionTree(false)
3377  , mIntersectionAccessor(intersectionTree)
3378  , mIsovalue(iso)
3379 {
3380 }
3381 
3382 
3383 template<typename InputTreeType>
3384 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3385  MaskIntersectingVoxels& rhs, tbb::split)
3386  : mInputAccessor(rhs.mInputAccessor.tree())
3387  , mNodes(rhs.mNodes)
3388  , mIntersectionTree(false)
3389  , mIntersectionAccessor(mIntersectionTree) // use local tree.
3390  , mIsovalue(rhs.mIsovalue)
3391 {
3392 }
3393 
3394 
3395 template<typename InputTreeType>
3396 void
3397 MaskIntersectingVoxels<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3398 {
3399  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3400  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3401  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3402 
3403  Coord ijk;
3404 
3405  for (size_t n = range.begin(); n != range.end(); ++n) {
3406 
3407  const BoolLeafNodeType& node = *mNodes[n];
3408 
3409  for (auto it = node.cbeginValueOn(); it; ++it) {
3410 
3411  if (!it.getValue()) {
3412 
3413  ijk = it.getCoord();
3414 
3415  const bool inside = isInsideValue(mInputAccessor.getValue(ijk), mIsovalue);
3416 
3417  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(1, 0, 0)), mIsovalue)) {
3418  xEdgeAcc.set(ijk);
3419  }
3420 
3421  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 1, 0)), mIsovalue)) {
3422  yEdgeAcc.set(ijk);
3423  }
3424 
3425  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 0, 1)), mIsovalue)) {
3426  zEdgeAcc.set(ijk);
3427  }
3428  }
3429  }
3430  }
3431 } // MaskIntersectingVoxels::operator()
3432 
3433 
3434 template<typename BoolTreeType>
3435 struct MaskBorderVoxels
3436 {
3437  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3438 
3439  MaskBorderVoxels(const BoolTreeType& maskTree,
3440  const std::vector<BoolLeafNodeType*>& maskNodes,
3441  BoolTreeType& borderTree)
3442  : mMaskTree(&maskTree)
3443  , mMaskNodes(maskNodes.data())
3444  , mTmpBorderTree(false)
3445  , mBorderTree(&borderTree) {}
3446 
3447  MaskBorderVoxels(MaskBorderVoxels& rhs, tbb::split)
3448  : mMaskTree(rhs.mMaskTree)
3449  , mMaskNodes(rhs.mMaskNodes)
3450  , mTmpBorderTree(false)
3451  , mBorderTree(&mTmpBorderTree) {}
3452 
3453  void join(MaskBorderVoxels& rhs) { mBorderTree->merge(*rhs.mBorderTree); }
3454 
3455  void operator()(const tbb::blocked_range<size_t>& range)
3456  {
3457  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
3458  tree::ValueAccessor<BoolTreeType> borderAcc(*mBorderTree);
3459  Coord ijk;
3460 
3461  for (size_t n = range.begin(); n != range.end(); ++n) {
3462 
3463  const BoolLeafNodeType& node = *mMaskNodes[n];
3464 
3465  for (auto it = node.cbeginValueOn(); it; ++it) {
3466 
3467  ijk = it.getCoord();
3468 
3469  const bool lhs = it.getValue();
3470  bool rhs = lhs;
3471 
3472  bool isEdgeVoxel = false;
3473 
3474  ijk[2] += 1; // i, j, k+1
3475  isEdgeVoxel = (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3476 
3477  ijk[1] += 1; // i, j+1, k+1
3478  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3479 
3480  ijk[0] += 1; // i+1, j+1, k+1
3481  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3482 
3483  ijk[1] -= 1; // i+1, j, k+1
3484  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3485 
3486 
3487  ijk[2] -= 1; // i+1, j, k
3488  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3489 
3490  ijk[1] += 1; // i+1, j+1, k
3491  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3492 
3493  ijk[0] -= 1; // i, j+1, k
3494  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3495 
3496  if (isEdgeVoxel) {
3497  ijk[1] -= 1; // i, j, k
3498  borderAcc.setValue(ijk, true);
3499  }
3500  }
3501  }
3502  }
3503 
3504 private:
3505  BoolTreeType const * const mMaskTree;
3506  BoolLeafNodeType const * const * const mMaskNodes;
3507 
3508  BoolTreeType mTmpBorderTree;
3509  BoolTreeType * const mBorderTree;
3510 }; // struct MaskBorderVoxels
3511 
3512 
3513 template<typename BoolTreeType>
3514 struct SyncMaskValues
3515 {
3516  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3517 
3518  SyncMaskValues(const std::vector<BoolLeafNodeType*>& nodes, const BoolTreeType& mask)
3519  : mNodes(nodes.data())
3520  , mMaskTree(&mask) {}
3521 
3522  void operator()(const tbb::blocked_range<size_t>& range) const
3523  {
3524  using ValueOnIter = typename BoolLeafNodeType::ValueOnIter;
3525 
3526  tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3527 
3528  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3529 
3530  BoolLeafNodeType& node = *mNodes[n];
3531 
3532  const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3533 
3534  if (maskNode) {
3535  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3536  const Index pos = it.pos();
3537  if (maskNode->getValue(pos)) {
3538  node.setValueOnly(pos, true);
3539  }
3540  }
3541  }
3542  }
3543  }
3544 
3545 private:
3546  BoolLeafNodeType * const * const mNodes;
3547  BoolTreeType const * const mMaskTree;
3548 }; // struct SyncMaskValues
3549 
3550 
3551 ////////////////////////////////////////
3552 
3553 
3554 template<typename BoolTreeType>
3555 struct MaskSurface
3556 {
3557  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3558 
3559  MaskSurface(const std::vector<BoolLeafNodeType*>& nodes,
3560  const BoolTreeType& mask,
3561  const math::Transform& inputTransform,
3562  const math::Transform& maskTransform,
3563  const bool invert)
3564  : mNodes(nodes.data())
3565  , mMaskTree(&mask)
3566  , mInputTransform(inputTransform)
3567  , mMaskTransform(maskTransform)
3568  , mMatchingTransforms(mInputTransform == mMaskTransform)
3569  , mInvertMask(invert) {}
3570 
3571  void operator()(const tbb::blocked_range<size_t>& range) const
3572  {
3573  using ValueOnIter = typename BoolLeafNodeType::ValueOnIter;
3574 
3575  tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3576 
3577  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3578 
3579  BoolLeafNodeType& node = *mNodes[n];
3580 
3581  if (mMatchingTransforms) {
3582 
3583  const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3584 
3585  if (maskNode) {
3586 
3587  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3588  const Index pos = it.pos();
3589  if (maskNode->isValueOn(pos) == mInvertMask) {
3590  node.setValueOnly(pos, true);
3591  }
3592  }
3593 
3594  } else {
3595 
3596  if (maskTreeAcc.isValueOn(node.origin()) == mInvertMask) {
3597  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3598  node.setValueOnly(it.pos(), true);
3599  }
3600  }
3601 
3602  }
3603 
3604  } else {
3605 
3606  Coord ijk;
3607 
3608  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3609 
3610  ijk = mMaskTransform.worldToIndexCellCentered(
3611  mInputTransform.indexToWorld(it.getCoord()));
3612 
3613  if (maskTreeAcc.isValueOn(ijk) == mInvertMask) {
3614  node.setValueOnly(it.pos(), true);
3615  }
3616  }
3617 
3618  }
3619  }
3620  }
3621 
3622 private:
3623  BoolLeafNodeType * const * const mNodes;
3624  BoolTreeType const * const mMaskTree;
3625  const math::Transform& mInputTransform;
3626  const math::Transform& mMaskTransform;
3627  const bool mMatchingTransforms;
3628  const bool mInvertMask;
3629 }; // struct MaskSurface
3630 
3631 
3632 template<typename InputGridType>
3633 inline void
3634 applySurfaceMask(
3635  typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3636  typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3637  const InputGridType& inputGrid,
3638  const GridBase::ConstPtr& maskGrid,
3639  const bool invertMask,
3640  const typename InputGridType::ValueType isovalue)
3641 {
3642  using InputTreeType = typename InputGridType::TreeType;
3643  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3644  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3645  using BoolGridType = Grid<BoolTreeType>;
3646 
3647  if (!maskGrid) return;
3648  if (maskGrid->type() != BoolGridType::gridType()) return;
3649 
3650  const math::Transform& transform = inputGrid.transform();
3651  const InputTreeType& inputTree = inputGrid.tree();
3652 
3653  const BoolGridType * surfaceMask = static_cast<const BoolGridType*>(maskGrid.get());
3654 
3655  const BoolTreeType& maskTree = surfaceMask->tree();
3656  const math::Transform& maskTransform = surfaceMask->transform();
3657 
3658  // mark masked voxels
3659 
3660  std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3661  intersectionTree.getNodes(intersectionLeafNodes);
3662 
3663  const tbb::blocked_range<size_t> intersectionRange(0, intersectionLeafNodes.size());
3664 
3665  tbb::parallel_for(intersectionRange,
3666  MaskSurface<BoolTreeType>(
3667  intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3668 
3669 
3670  // mask surface-mask border
3671 
3672  MaskBorderVoxels<BoolTreeType> borderOp(
3673  intersectionTree, intersectionLeafNodes, borderTree);
3674  tbb::parallel_reduce(intersectionRange, borderOp);
3675 
3676 
3677  // recompute isosurface intersection mask
3678 
3679  BoolTreeType tmpIntersectionTree(false);
3680 
3681  MaskIntersectingVoxels<InputTreeType> op(
3682  inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3683 
3684  tbb::parallel_reduce(intersectionRange, op);
3685 
3686  std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3687  tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3688 
3689  tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3690  SyncMaskValues<BoolTreeType>(tmpIntersectionLeafNodes, intersectionTree));
3691 
3692  intersectionTree.clear();
3693  intersectionTree.merge(tmpIntersectionTree);
3694 }
3695 
3696 
3697 ////////////////////////////////////////
3698 
3699 
3700 template<typename InputTreeType>
3701 struct ComputeAuxiliaryData
3702 {
3703  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3704  using InputValueType = typename InputLeafNodeType::ValueType;
3705 
3706  using BoolLeafNodeType = tree::LeafNode<bool, InputLeafNodeType::LOG2DIM>;
3707 
3708  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
3709  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
3710 
3711 
3712  ComputeAuxiliaryData(const InputTreeType& inputTree,
3713  const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3714  Int16TreeType& signFlagsTree,
3715  Index32TreeType& pointIndexTree,
3716  InputValueType iso);
3717 
3718  ComputeAuxiliaryData(ComputeAuxiliaryData&, tbb::split);
3719  void operator()(const tbb::blocked_range<size_t>&);
3720  void join(const ComputeAuxiliaryData& rhs) {
3721  mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.tree());
3722  mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.tree());
3723  }
3724 
3725 private:
3727  BoolLeafNodeType const * const * const mIntersectionNodes;
3728 
3729  Int16TreeType mSignFlagsTree;
3730  tree::ValueAccessor<Int16TreeType> mSignFlagsAccessor;
3731  Index32TreeType mPointIndexTree;
3732  tree::ValueAccessor<Index32TreeType> mPointIndexAccessor;
3733 
3734  const InputValueType mIsovalue;
3735 };
3736 
3737 
3738 template<typename InputTreeType>
3739 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(
3740  const InputTreeType& inputTree,
3741  const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3742  Int16TreeType& signFlagsTree,
3743  Index32TreeType& pointIndexTree,
3744  InputValueType iso)
3745  : mInputAccessor(inputTree)
3746  , mIntersectionNodes(intersectionLeafNodes.data())
3747  , mSignFlagsTree(0)
3748  , mSignFlagsAccessor(signFlagsTree)
3749  , mPointIndexTree(std::numeric_limits<Index32>::max())
3750  , mPointIndexAccessor(pointIndexTree)
3751  , mIsovalue(iso)
3752 {
3753  pointIndexTree.root().setBackground(std::numeric_limits<Index32>::max(), false);
3754 }
3755 
3756 
3757 template<typename InputTreeType>
3758 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(ComputeAuxiliaryData& rhs, tbb::split)
3759  : mInputAccessor(rhs.mInputAccessor.tree())
3760  , mIntersectionNodes(rhs.mIntersectionNodes)
3761  , mSignFlagsTree(0)
3762  , mSignFlagsAccessor(mSignFlagsTree)
3763  , mPointIndexTree(std::numeric_limits<Index32>::max())
3764  , mPointIndexAccessor(mPointIndexTree)
3765  , mIsovalue(rhs.mIsovalue)
3766 {
3767 }
3768 
3769 
3770 template<typename InputTreeType>
3771 void
3772 ComputeAuxiliaryData<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3773 {
3774  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
3775 
3776  Coord ijk;
3777  std::array<InputValueType, 8> cellVertexValues;
3778  std::unique_ptr<Int16LeafNodeType> signsNodePt(nullptr);
3779 
3780  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3781 
3782  const BoolLeafNodeType& maskNode = *mIntersectionNodes[n];
3783  const Coord& origin = maskNode.origin();
3784 
3785  const InputLeafNodeType* leafPt = mInputAccessor.probeConstLeaf(origin);
3786 
3787  if (!signsNodePt.get()) signsNodePt.reset(new Int16LeafNodeType(origin, 0));
3788  else signsNodePt->setOrigin(origin);
3789 
3790  bool updatedNode = false;
3791 
3792  for (auto it = maskNode.cbeginValueOn(); it; ++it) {
3793 
3794  const Index pos = it.pos();
3795  ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3796 
3797  const bool inclusiveCell = leafPt && isInternalLeafCoord<InputLeafNodeType>(ijk);
3798 
3799  if (inclusiveCell) getCellVertexValues(*leafPt, pos, cellVertexValues);
3800  else getCellVertexValues(mInputAccessor, origin + ijk, cellVertexValues);
3801 
3802  uint8_t signFlags = computeSignFlags(cellVertexValues, mIsovalue);
3803 
3804  if (signFlags != 0 && signFlags != 0xFF) {
3805 
3806  const bool inside = signFlags & 0x1;
3807 
3808  int edgeFlags = inside ? INSIDE : 0;
3809 
3810  if (!it.getValue()) {
3811  edgeFlags |= inside != ((signFlags & 0x02) != 0) ? XEDGE : 0;
3812  edgeFlags |= inside != ((signFlags & 0x10) != 0) ? YEDGE : 0;
3813  edgeFlags |= inside != ((signFlags & 0x08) != 0) ? ZEDGE : 0;
3814  }
3815 
3816  const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3817  if (ambiguousCellFlags != 0) {
3818  correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor,
3819  origin + ijk, mIsovalue);
3820  }
3821 
3822  edgeFlags |= int(signFlags);
3823 
3824  signsNodePt->setValueOn(pos, Int16(edgeFlags));
3825  updatedNode = true;
3826  }
3827  }
3828 
3829  if (updatedNode) {
3830  typename Index32TreeType::LeafNodeType* idxNode =
3831  mPointIndexAccessor.touchLeaf(origin);
3832  idxNode->topologyUnion(*signsNodePt);
3833 
3834  // zero fill
3835  auto* const idxData = idxNode->buffer().data();
3836  for (auto it = idxNode->beginValueOn(); it; ++it) {
3837  idxData[it.pos()] = 0;
3838  }
3839 
3840  mSignFlagsAccessor.addLeaf(signsNodePt.release());
3841  }
3842  }
3843 } // ComputeAuxiliaryData::operator()
3844 
3845 
3846 template<typename InputTreeType>
3847 inline void
3848 computeAuxiliaryData(
3849  typename InputTreeType::template ValueConverter<Int16>::Type& signFlagsTree,
3850  typename InputTreeType::template ValueConverter<Index32>::Type& pointIndexTree,
3851  const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3852  const InputTreeType& inputTree,
3853  typename InputTreeType::ValueType isovalue)
3854 {
3855  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3856  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3857 
3858  std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3859  intersectionTree.getNodes(intersectionLeafNodes);
3860 
3861  ComputeAuxiliaryData<InputTreeType> op(
3862  inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3863 
3864  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3865 }
3866 
3867 
3868 ////////////////////////////////////////
3869 
3870 
3871 template<Index32 LeafNodeLog2Dim>
3872 struct LeafNodePointCount
3873 {
3874  using Int16LeafNodeType = tree::LeafNode<Int16, LeafNodeLog2Dim>;
3875 
3876  LeafNodePointCount(const std::vector<Int16LeafNodeType*>& leafNodes,
3877  std::unique_ptr<Index32[]>& leafNodeCount)
3878  : mLeafNodes(leafNodes.data())
3879  , mData(leafNodeCount.get()) {}
3880 
3881  void operator()(const tbb::blocked_range<size_t>& range) const {
3882 
3883  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3884 
3885  Index32 count = 0;
3886 
3887  Int16 const * p = mLeafNodes[n]->buffer().data();
3888  Int16 const * const endP = p + Int16LeafNodeType::SIZE;
3889 
3890  while (p < endP) {
3891  count += Index32(sEdgeGroupTable[(SIGNS & *p)][0]);
3892  ++p;
3893  }
3894 
3895  mData[n] = count;
3896  }
3897  }
3898 
3899 private:
3900  Int16LeafNodeType * const * const mLeafNodes;
3901  Index32 *mData;
3902 }; // struct LeafNodePointCount
3903 
3904 
3905 template<typename PointIndexLeafNode>
3906 struct AdaptiveLeafNodePointCount
3907 {
3908  using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3909 
3910  AdaptiveLeafNodePointCount(const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3911  const std::vector<Int16LeafNodeType*>& signDataNodes,
3912  std::unique_ptr<Index32[]>& leafNodeCount)
3913  : mPointIndexNodes(pointIndexNodes.data())
3914  , mSignDataNodes(signDataNodes.data())
3915  , mData(leafNodeCount.get()) {}
3916 
3917  void operator()(const tbb::blocked_range<size_t>& range) const
3918  {
3919  using IndexType = typename PointIndexLeafNode::ValueType;
3920 
3921  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3922 
3923  const PointIndexLeafNode& node = *mPointIndexNodes[n];
3924  const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3925 
3926  size_t count = 0;
3927 
3928  std::set<IndexType> uniqueRegions;
3929 
3930  for (typename PointIndexLeafNode::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3931 
3932  IndexType id = it.getValue();
3933 
3934  if (id == 0) {
3935  count += size_t(sEdgeGroupTable[(SIGNS & signNode.getValue(it.pos()))][0]);
3936  } else if (id != IndexType(util::INVALID_IDX)) {
3937  uniqueRegions.insert(id);
3938  }
3939  }
3940 
3941  mData[n] = Index32(count + uniqueRegions.size());
3942  }
3943  }
3944 
3945 private:
3946  PointIndexLeafNode const * const * const mPointIndexNodes;
3947  Int16LeafNodeType const * const * const mSignDataNodes;
3948  Index32 *mData;
3949 }; // struct AdaptiveLeafNodePointCount
3950 
3951 
3952 template<typename PointIndexLeafNode>
3953 struct MapPoints
3954 {
3955  using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3956 
3957  MapPoints(const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3958  const std::vector<Int16LeafNodeType*>& signDataNodes,
3959  std::unique_ptr<Index32[]>& leafNodeCount)
3960  : mPointIndexNodes(pointIndexNodes.data())
3961  , mSignDataNodes(signDataNodes.data())
3962  , mData(leafNodeCount.get()) {}
3963 
3964  void operator()(const tbb::blocked_range<size_t>& range) const {
3965 
3966  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3967 
3968  const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3969  PointIndexLeafNode& indexNode = *mPointIndexNodes[n];
3970 
3971  Index32 pointOffset = mData[n];
3972 
3973  for (auto it = indexNode.beginValueOn(); it; ++it) {
3974  const Index pos = it.pos();
3975  indexNode.setValueOnly(pos, pointOffset);
3976  const int signs = SIGNS & int(signNode.getValue(pos));
3977  pointOffset += Index32(sEdgeGroupTable[signs][0]);
3978  }
3979  }
3980  }
3981 
3982 private:
3983  PointIndexLeafNode * const * const mPointIndexNodes;
3984  Int16LeafNodeType const * const * const mSignDataNodes;
3985  Index32 * const mData;
3986 }; // struct MapPoints
3987 
3988 
3989 template<typename TreeType, typename PrimBuilder>
3990 struct ComputePolygons
3991 {
3992  using Int16TreeType = typename TreeType::template ValueConverter<Int16>::Type;
3993  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
3994 
3995  using Index32TreeType = typename TreeType::template ValueConverter<Index32>::Type;
3996  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
3997 
3998  ComputePolygons(
3999  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4000  const Int16TreeType& signFlagsTree,
4001  const Index32TreeType& idxTree,
4002  PolygonPoolList& polygons,
4003  bool invertSurfaceOrientation);
4004 
4005  void setRefSignTree(const Int16TreeType * r) { mRefSignFlagsTree = r; }
4006 
4007  void operator()(const tbb::blocked_range<size_t>&) const;
4008 
4009 private:
4010  Int16LeafNodeType * const * const mSignFlagsLeafNodes;
4011  Int16TreeType const * const mSignFlagsTree;
4012  Int16TreeType const * mRefSignFlagsTree;
4013  Index32TreeType const * const mIndexTree;
4014  PolygonPoolList * const mPolygonPoolList;
4015  bool const mInvertSurfaceOrientation;
4016 }; // struct ComputePolygons
4017 
4018 
4019 template<typename TreeType, typename PrimBuilder>
4020 ComputePolygons<TreeType, PrimBuilder>::ComputePolygons(
4021  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4022  const Int16TreeType& signFlagsTree,
4023  const Index32TreeType& idxTree,
4024  PolygonPoolList& polygons,
4025  bool invertSurfaceOrientation)
4026  : mSignFlagsLeafNodes(signFlagsLeafNodes.data())
4027  , mSignFlagsTree(&signFlagsTree)
4028  , mRefSignFlagsTree(nullptr)
4029  , mIndexTree(&idxTree)
4030  , mPolygonPoolList(&polygons)
4031  , mInvertSurfaceOrientation(invertSurfaceOrientation)
4032 {
4033 }
4034 
4035 template<typename InputTreeType, typename PrimBuilder>
4036 void
4037 ComputePolygons<InputTreeType, PrimBuilder>::operator()(const tbb::blocked_range<size_t>& range) const
4038 {
4039  using Int16ValueAccessor = tree::ValueAccessor<const Int16TreeType>;
4040  Int16ValueAccessor signAcc(*mSignFlagsTree);
4041 
4042  tree::ValueAccessor<const Index32TreeType> idxAcc(*mIndexTree);
4043 
4044  const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4045 
4046  PrimBuilder mesher;
4047  size_t edgeCount;
4048  Coord ijk, origin;
4049 
4050  // reference data
4051  std::unique_ptr<Int16ValueAccessor> refSignAcc;
4052  if (mRefSignFlagsTree) refSignAcc.reset(new Int16ValueAccessor(*mRefSignFlagsTree));
4053 
4054  for (size_t n = range.begin(); n != range.end(); ++n) {
4055 
4056  const Int16LeafNodeType& node = *mSignFlagsLeafNodes[n];
4057  origin = node.origin();
4058 
4059  // Get an upper bound on the number of primitives.
4060  edgeCount = 0;
4061  typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4062  for (; iter; ++iter) {
4063  if (iter.getValue() & XEDGE) ++edgeCount;
4064  if (iter.getValue() & YEDGE) ++edgeCount;
4065  if (iter.getValue() & ZEDGE) ++edgeCount;
4066  }
4067 
4068  if (edgeCount == 0) continue;
4069 
4070  mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4071 
4072  const Int16LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
4073  const Index32LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
4074 
4075  if (!signleafPt || !idxLeafPt) continue;
4076 
4077  const Int16LeafNodeType *refSignLeafPt = nullptr;
4078  if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4079 
4080  Vec3i offsets;
4081 
4082  for (iter = node.cbeginValueOn(); iter; ++iter) {
4083  ijk = iter.getCoord();
4084 
4085  const Int16 flags = iter.getValue();
4086  if (!(flags & 0xE00)) continue;
4087 
4088  Int16 refFlags = 0;
4089  if (refSignLeafPt) {
4090  refFlags = refSignLeafPt->getValue(iter.pos());
4091  }
4092 
4093  const uint8_t cell = uint8_t(SIGNS & flags);
4094 
4095  if (sEdgeGroupTable[cell][0] > 1) {
4096  offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4097  offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4098  offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4099  }
4100  else {
4101  offsets.setZero();
4102  }
4103 
4104  if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4105  constructPolygons(invertSurfaceOrientation,
4106  flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4107  } else {
4108  constructPolygons(invertSurfaceOrientation,
4109  flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4110  }
4111  }
4112 
4113  mesher.done();
4114  }
4115 
4116 } // ComputePolygons::operator()
4117 
4118 
4119 ////////////////////////////////////////
4120 
4121 
4122 template<typename T>
4123 struct CopyArray
4124 {
4125  CopyArray(T * outputArray, const T * inputArray, size_t outputOffset = 0)
4126  : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4127  {
4128  }
4129 
4130  void operator()(const tbb::blocked_range<size_t>& inputArrayRange) const
4131  {
4132  const size_t offset = mOutputOffset;
4133  for (size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n < N; ++n) {
4134  mOutputArray[offset + n] = mInputArray[n];
4135  }
4136  }
4137 
4138 private:
4139  T * const mOutputArray;
4140  T const * const mInputArray;
4141  size_t const mOutputOffset;
4142 }; // struct CopyArray
4143 
4144 
4145 struct FlagAndCountQuadsToSubdivide
4146 {
4147  FlagAndCountQuadsToSubdivide(PolygonPoolList& polygons,
4148  const std::vector<uint8_t>& pointFlags,
4149  std::unique_ptr<openvdb::Vec3s[]>& points,
4150  std::unique_ptr<unsigned[]>& numQuadsToDivide)
4151  : mPolygonPoolList(&polygons)
4152  , mPointFlags(pointFlags.data())
4153  , mPoints(points.get())
4154  , mNumQuadsToDivide(numQuadsToDivide.get())
4155  {
4156  }
4157 
4158  void operator()(const tbb::blocked_range<size_t>& range) const
4159  {
4160  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4161 
4162  PolygonPool& polygons = (*mPolygonPoolList)[n];
4163 
4164  unsigned count = 0;
4165 
4166  // count and tag nonplanar seam line quads.
4167  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4168 
4169  char& flags = polygons.quadFlags(i);
4170 
4171  if ((flags & POLYFLAG_FRACTURE_SEAM) && !(flags & POLYFLAG_EXTERIOR)) {
4172 
4173  Vec4I& quad = polygons.quad(i);
4174 
4175  const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4176  || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4177 
4178  if (!edgePoly) continue;
4179 
4180  const Vec3s& p0 = mPoints[quad[0]];
4181  const Vec3s& p1 = mPoints[quad[1]];
4182  const Vec3s& p2 = mPoints[quad[2]];
4183  const Vec3s& p3 = mPoints[quad[3]];
4184 
4185  if (!isPlanarQuad(p0, p1, p2, p3, 1e-6f)) {
4186  flags |= POLYFLAG_SUBDIVIDED;
4187  count++;
4188  }
4189  }
4190  }
4191 
4192  mNumQuadsToDivide[n] = count;
4193  }
4194  }
4195 
4196 private:
4197  PolygonPoolList * const mPolygonPoolList;
4198  uint8_t const * const mPointFlags;
4199  Vec3s const * const mPoints;
4200  unsigned * const mNumQuadsToDivide;
4201 }; // struct FlagAndCountQuadsToSubdivide
4202 
4203 
4204 struct SubdivideQuads
4205 {
4206  SubdivideQuads(PolygonPoolList& polygons,
4207  const std::unique_ptr<openvdb::Vec3s[]>& points,
4208  size_t pointCount,
4209  std::unique_ptr<openvdb::Vec3s[]>& centroids,
4210  std::unique_ptr<unsigned[]>& numQuadsToDivide,
4211  std::unique_ptr<unsigned[]>& centroidOffsets)
4212  : mPolygonPoolList(&polygons)
4213  , mPoints(points.get())
4214  , mCentroids(centroids.get())
4215  , mNumQuadsToDivide(numQuadsToDivide.get())
4216  , mCentroidOffsets(centroidOffsets.get())
4217  , mPointCount(pointCount)
4218  {
4219  }
4220 
4221  void operator()(const tbb::blocked_range<size_t>& range) const
4222  {
4223  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4224 
4225  PolygonPool& polygons = (*mPolygonPoolList)[n];
4226 
4227  const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4228 
4229  if (nonplanarCount > 0) {
4230 
4231  PolygonPool tmpPolygons;
4232  tmpPolygons.resetQuads(polygons.numQuads() - nonplanarCount);
4233  tmpPolygons.resetTriangles(polygons.numTriangles() + size_t(4) * nonplanarCount);
4234 
4235  size_t offset = mCentroidOffsets[n];
4236 
4237  size_t triangleIdx = 0;
4238 
4239  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4240 
4241  const char quadFlags = polygons.quadFlags(i);
4242  if (!(quadFlags & POLYFLAG_SUBDIVIDED)) continue;
4243 
4244  unsigned newPointIdx = unsigned(offset + mPointCount);
4245 
4246  openvdb::Vec4I& quad = polygons.quad(i);
4247 
4248  mCentroids[offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4249  mPoints[quad[2]] + mPoints[quad[3]]) * 0.25f;
4250 
4251  ++offset;
4252 
4253  {
4254  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4255 
4256  triangle[0] = quad[0];
4257  triangle[1] = newPointIdx;
4258  triangle[2] = quad[3];
4259 
4260  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4261  }
4262 
4263  ++triangleIdx;
4264 
4265  {
4266  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4267 
4268  triangle[0] = quad[0];
4269  triangle[1] = quad[1];
4270  triangle[2] = newPointIdx;
4271 
4272  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4273  }
4274 
4275  ++triangleIdx;
4276 
4277  {
4278  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4279 
4280  triangle[0] = quad[1];
4281  triangle[1] = quad[2];
4282  triangle[2] = newPointIdx;
4283 
4284  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4285  }
4286 
4287 
4288  ++triangleIdx;
4289 
4290  {
4291  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4292 
4293  triangle[0] = quad[2];
4294  triangle[1] = quad[3];
4295  triangle[2] = newPointIdx;
4296 
4297  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4298  }
4299 
4300  ++triangleIdx;
4301 
4302  quad[0] = util::INVALID_IDX; // mark for deletion
4303  }
4304 
4305 
4306  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4307  tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
4308  tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
4309  ++triangleIdx;
4310  }
4311 
4312  size_t quadIdx = 0;
4313  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4314  openvdb::Vec4I& quad = polygons.quad(i);
4315 
4316  if (quad[0] != util::INVALID_IDX) { // ignore invalid quads
4317  tmpPolygons.quad(quadIdx) = quad;
4318  tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
4319  ++quadIdx;
4320  }
4321  }
4322 
4323  polygons.copy(tmpPolygons);
4324  }
4325  }
4326  }
4327 
4328 private:
4329  PolygonPoolList * const mPolygonPoolList;
4330  Vec3s const * const mPoints;
4331  Vec3s * const mCentroids;
4332  unsigned * const mNumQuadsToDivide;
4333  unsigned * const mCentroidOffsets;
4334  size_t const mPointCount;
4335 }; // struct SubdivideQuads
4336 
4337 
4338 inline void
4339 subdivideNonplanarSeamLineQuads(
4340  PolygonPoolList& polygonPoolList,
4341  size_t polygonPoolListSize,
4342  PointList& pointList,
4343  size_t& pointListSize,
4344  std::vector<uint8_t>& pointFlags)
4345 {
4346  const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4347 
4348  std::unique_ptr<unsigned[]> numQuadsToDivide(new unsigned[polygonPoolListSize]);
4349 
4350  tbb::parallel_for(polygonPoolListRange,
4351  FlagAndCountQuadsToSubdivide(polygonPoolList, pointFlags, pointList, numQuadsToDivide));
4352 
4353  std::unique_ptr<unsigned[]> centroidOffsets(new unsigned[polygonPoolListSize]);
4354 
4355  size_t centroidCount = 0;
4356 
4357  {
4358  unsigned sum = 0;
4359  for (size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4360  centroidOffsets[n] = sum;
4361  sum += numQuadsToDivide[n];
4362  }
4363  centroidCount = size_t(sum);
4364  }
4365 
4366  std::unique_ptr<Vec3s[]> centroidList(new Vec3s[centroidCount]);
4367 
4368  tbb::parallel_for(polygonPoolListRange,
4369  SubdivideQuads(polygonPoolList, pointList, pointListSize,
4370  centroidList, numQuadsToDivide, centroidOffsets));
4371 
4372  if (centroidCount > 0) {
4373 
4374  const size_t newPointListSize = centroidCount + pointListSize;
4375 
4376  std::unique_ptr<openvdb::Vec3s[]> newPointList(new openvdb::Vec3s[newPointListSize]);
4377 
4378  tbb::parallel_for(tbb::blocked_range<size_t>(0, pointListSize),
4379  CopyArray<Vec3s>(newPointList.get(), pointList.get()));
4380 
4381  tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4382  CopyArray<Vec3s>(newPointList.get(), centroidList.get(), pointListSize));
4383 
4384  pointListSize = newPointListSize;
4385  pointList.swap(newPointList);
4386  pointFlags.resize(pointListSize, 0);
4387  }
4388 }
4389 
4390 
4391 struct ReviseSeamLineFlags
4392 {
4393  ReviseSeamLineFlags(PolygonPoolList& polygons,
4394  const std::vector<uint8_t>& pointFlags)
4395  : mPolygonPoolList(&polygons)
4396  , mPointFlags(pointFlags.data())
4397  {
4398  }
4399 
4400  void operator()(const tbb::blocked_range<size_t>& range) const
4401  {
4402  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4403 
4404  PolygonPool& polygons = (*mPolygonPoolList)[n];
4405 
4406  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4407 
4408  char& flags = polygons.quadFlags(i);
4409 
4410  if (flags & POLYFLAG_FRACTURE_SEAM) {
4411 
4412  openvdb::Vec4I& verts = polygons.quad(i);
4413 
4414  const bool hasSeamLinePoint =
4415  mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4416  mPointFlags[verts[2]] || mPointFlags[verts[3]];
4417 
4418  if (!hasSeamLinePoint) {
4419  flags &= ~POLYFLAG_FRACTURE_SEAM;
4420  }
4421  }
4422  } // end quad loop
4423 
4424  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4425 
4426  char& flags = polygons.triangleFlags(i);
4427 
4428  if (flags & POLYFLAG_FRACTURE_SEAM) {
4429 
4430  openvdb::Vec3I& verts = polygons.triangle(i);
4431 
4432  const bool hasSeamLinePoint =
4433  mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4434 
4435  if (!hasSeamLinePoint) {
4436  flags &= ~POLYFLAG_FRACTURE_SEAM;
4437  }
4438 
4439  }
4440  } // end triangle loop
4441 
4442  } // end polygon pool loop
4443  }
4444 
4445 private:
4446  PolygonPoolList * const mPolygonPoolList;
4447  uint8_t const * const mPointFlags;
4448 }; // struct ReviseSeamLineFlags
4449 
4450 
4451 inline void
4452 reviseSeamLineFlags(PolygonPoolList& polygonPoolList, size_t polygonPoolListSize,
4453  std::vector<uint8_t>& pointFlags)
4454 {
4455  tbb::parallel_for(tbb::blocked_range<size_t>(0, polygonPoolListSize),
4456  ReviseSeamLineFlags(polygonPoolList, pointFlags));
4457 }
4458 
4459 
4460 ////////////////////////////////////////
4461 
4462 
4463 template<typename InputTreeType>
4464 struct MaskDisorientedTrianglePoints
4465 {
4466  MaskDisorientedTrianglePoints(const InputTreeType& inputTree, const PolygonPoolList& polygons,
4467  const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4468  const math::Transform& transform, bool invertSurfaceOrientation)
4469  : mInputTree(&inputTree)
4470  , mPolygonPoolList(&polygons)
4471  , mPointList(&pointList)
4472  , mPointMask(pointMask.get())
4473  , mTransform(transform)
4474  , mInvertSurfaceOrientation(invertSurfaceOrientation)
4475  {
4476  }
4477 
4478  void operator()(const tbb::blocked_range<size_t>& range) const
4479  {
4480  using ValueType = typename InputTreeType::LeafNodeType::ValueType;
4481 
4482  tree::ValueAccessor<const InputTreeType> inputAcc(*mInputTree);
4483  Vec3s centroid, normal;
4484  Coord ijk;
4485 
4486  const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4487 
4488  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4489 
4490  const PolygonPool& polygons = (*mPolygonPoolList)[n];
4491 
4492  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4493 
4494  const Vec3I& verts = polygons.triangle(i);
4495 
4496  const Vec3s& v0 = (*mPointList)[verts[0]];
4497  const Vec3s& v1 = (*mPointList)[verts[1]];
4498  const Vec3s& v2 = (*mPointList)[verts[2]];
4499 
4500  normal = (v2 - v0).cross((v1 - v0));
4501  normal.normalize();
4502 
4503  centroid = (v0 + v1 + v2) * (1.0f / 3.0f);
4504  ijk = mTransform.worldToIndexCellCentered(centroid);
4505 
4506  Vec3s dir( math::ISGradient<math::CD_2ND>::result(inputAcc, ijk) );
4507  dir.normalize();
4508 
4509  if (invertGradientDir) {
4510  dir = -dir;
4511  }
4512 
4513  // check if the angle is obtuse
4514  if (dir.dot(normal) < -0.5f) {
4515  // Concurrent writes to same memory address can occur, but
4516  // all threads are writing the same value and char is atomic.
4517  // (It is extremely rare that disoriented triangles share points,
4518  // false sharing related performance impacts are not a concern.)
4519  mPointMask[verts[0]] = 1;
4520  mPointMask[verts[1]] = 1;
4521  mPointMask[verts[2]] = 1;
4522  }
4523 
4524  } // end triangle loop
4525 
4526  } // end polygon pool loop
4527  }
4528 
4529 private:
4530  InputTreeType const * const mInputTree;
4531  PolygonPoolList const * const mPolygonPoolList;
4532  PointList const * const mPointList;
4533  uint8_t * const mPointMask;
4534  math::Transform const mTransform;
4535  bool const mInvertSurfaceOrientation;
4536 }; // struct MaskDisorientedTrianglePoints
4537 
4538 
4539 template<typename InputTree>
4540 inline void
4541 relaxDisorientedTriangles(
4542  bool invertSurfaceOrientation,
4543  const InputTree& inputTree,
4544  const math::Transform& transform,
4545  PolygonPoolList& polygonPoolList,
4546  size_t polygonPoolListSize,
4547  PointList& pointList,
4548  const size_t pointListSize)
4549 {
4550  const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4551 
4552  std::unique_ptr<uint8_t[]> pointMask(new uint8_t[pointListSize]);
4553  fillArray(pointMask.get(), uint8_t(0), pointListSize);
4554 
4555  tbb::parallel_for(polygonPoolListRange,
4556  MaskDisorientedTrianglePoints<InputTree>(
4557  inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4558 
4559  std::unique_ptr<uint8_t[]> pointUpdates(new uint8_t[pointListSize]);
4560  fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4561 
4562  std::unique_ptr<Vec3s[]> newPoints(new Vec3s[pointListSize]);
4563  fillArray(newPoints.get(), Vec3s(0.0f, 0.0f, 0.0f), pointListSize);
4564 
4565  for (size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4566 
4567  PolygonPool& polygons = polygonPoolList[n];
4568 
4569  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4570  openvdb::Vec4I& verts = polygons.quad(i);
4571 
4572  for (int v = 0; v < 4; ++v) {
4573 
4574  const unsigned pointIdx = verts[v];
4575 
4576  if (pointMask[pointIdx] == 1) {
4577 
4578  newPoints[pointIdx] +=
4579  pointList[verts[0]] + pointList[verts[1]] +
4580  pointList[verts[2]] + pointList[verts[3]];
4581 
4582  pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4583  }
4584  }
4585  }
4586 
4587  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4588  openvdb::Vec3I& verts = polygons.triangle(i);
4589 
4590  for (int v = 0; v < 3; ++v) {
4591 
4592  const unsigned pointIdx = verts[v];
4593 
4594  if (pointMask[pointIdx] == 1) {
4595  newPoints[pointIdx] +=
4596  pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4597 
4598  pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4599  }
4600  }
4601  }
4602  }
4603 
4604  for (size_t n = 0, N = pointListSize; n < N; ++n) {
4605  if (pointUpdates[n] > 0) {
4606  const double weight = 1.0 / double(pointUpdates[n]);
4607  pointList[n] = newPoints[n] * float(weight);
4608  }
4609  }
4610 }
4611 
4612 
4613 template<typename GridType>
4614 void
4615 doVolumeToMesh(
4616  const GridType& grid,
4617  std::vector<Vec3s>& points,
4618  std::vector<Vec3I>& triangles,
4619  std::vector<Vec4I>& quads,
4620  double isovalue,
4621  double adaptivity,
4622  bool relaxDisorientedTriangles)
4623 {
4624  static_assert(std::is_scalar<typename GridType::ValueType>::value,
4625  "volume to mesh conversion is supported only for scalar grids");
4626 
4627  VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
4628  mesher(grid);
4629 
4630  // Preallocate the point list
4631  points.clear();
4632  points.resize(mesher.pointListSize());
4633 
4634  { // Copy points
4635  volume_to_mesh_internal::PointListCopy ptnCpy(mesher.pointList(), points);
4636  tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
4637  mesher.pointList().reset(nullptr);
4638  }
4639 
4640  PolygonPoolList& polygonPoolList = mesher.polygonPoolList();
4641 
4642  { // Preallocate primitive lists
4643  size_t numQuads = 0, numTriangles = 0;
4644  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4645  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4646  numTriangles += polygons.numTriangles();
4647  numQuads += polygons.numQuads();
4648  }
4649 
4650  triangles.clear();
4651  triangles.resize(numTriangles);
4652  quads.clear();
4653  quads.resize(numQuads);
4654  }
4655 
4656  // Copy primitives
4657  size_t qIdx = 0, tIdx = 0;
4658  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4659  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4660 
4661  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4662  quads[qIdx++] = polygons.quad(i);
4663  }
4664 
4665  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4666  triangles[tIdx++] = polygons.triangle(i);
4667  }
4668  }
4669 }
4670 
4671 
4672 } // volume_to_mesh_internal namespace
4673 
4674 /// @endcond
4675 
4676 ////////////////////////////////////////
4677 
4678 
4679 inline
4680 PolygonPool::PolygonPool()
4681  : mNumQuads(0)
4682  , mNumTriangles(0)
4683  , mQuads(nullptr)
4684  , mTriangles(nullptr)
4685  , mQuadFlags(nullptr)
4686  , mTriangleFlags(nullptr)
4687 {
4688 }
4689 
4690 
4691 inline
4693  : mNumQuads(numQuads)
4694  , mNumTriangles(numTriangles)
4695  , mQuads(new openvdb::Vec4I[mNumQuads])
4696  , mTriangles(new openvdb::Vec3I[mNumTriangles])
4697  , mQuadFlags(new char[mNumQuads])
4698  , mTriangleFlags(new char[mNumTriangles])
4699 {
4700 }
4701 
4702 
4703 inline void
4705 {
4706  resetQuads(rhs.numQuads());
4708 
4709  for (size_t i = 0; i < mNumQuads; ++i) {
4710  mQuads[i] = rhs.mQuads[i];
4711  mQuadFlags[i] = rhs.mQuadFlags[i];
4712  }
4713 
4714  for (size_t i = 0; i < mNumTriangles; ++i) {
4715  mTriangles[i] = rhs.mTriangles[i];
4716  mTriangleFlags[i] = rhs.mTriangleFlags[i];
4717  }
4718 }
4719 
4720 
4721 inline void
4723 {
4724  mNumQuads = size;
4725  mQuads.reset(new openvdb::Vec4I[mNumQuads]);
4726  mQuadFlags.reset(new char[mNumQuads]);
4727 }
4728 
4729 
4730 inline void
4732 {
4733  mNumQuads = 0;
4734  mQuads.reset(nullptr);
4735  mQuadFlags.reset(nullptr);
4736 }
4737 
4738 
4739 inline void
4741 {
4742  mNumTriangles = size;
4743  mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
4744  mTriangleFlags.reset(new char[mNumTriangles]);
4745 }
4746 
4747 
4748 inline void
4750 {
4751  mNumTriangles = 0;
4752  mTriangles.reset(nullptr);
4753  mTriangleFlags.reset(nullptr);
4754 }
4755 
4756 
4757 inline bool
4758 PolygonPool::trimQuads(const size_t n, bool reallocate)
4759 {
4760  if (!(n < mNumQuads)) return false;
4761 
4762  if (reallocate) {
4763 
4764  if (n == 0) {
4765  mQuads.reset(nullptr);
4766  } else {
4767 
4768  std::unique_ptr<openvdb::Vec4I[]> quads(new openvdb::Vec4I[n]);
4769  std::unique_ptr<char[]> flags(new char[n]);
4770 
4771  for (size_t i = 0; i < n; ++i) {
4772  quads[i] = mQuads[i];
4773  flags[i] = mQuadFlags[i];
4774  }
4775 
4776  mQuads.swap(quads);
4777  mQuadFlags.swap(flags);
4778  }
4779  }
4780 
4781  mNumQuads = n;
4782  return true;
4783 }
4784 
4785 
4786 inline bool
4787 PolygonPool::trimTrinagles(const size_t n, bool reallocate)
4788 {
4789  if (!(n < mNumTriangles)) return false;
4790 
4791  if (reallocate) {
4792 
4793  if (n == 0) {
4794  mTriangles.reset(nullptr);
4795  } else {
4796 
4797  std::unique_ptr<openvdb::Vec3I[]> triangles(new openvdb::Vec3I[n]);
4798  std::unique_ptr<char[]> flags(new char[n]);
4799 
4800  for (size_t i = 0; i < n; ++i) {
4801  triangles[i] = mTriangles[i];
4802  flags[i] = mTriangleFlags[i];
4803  }
4804 
4805  mTriangles.swap(triangles);
4806  mTriangleFlags.swap(flags);
4807  }
4808  }
4809 
4810  mNumTriangles = n;
4811  return true;
4812 }
4813 
4814 
4815 ////////////////////////////////////////
4816 
4817 
4818 inline
4819 VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity, bool relaxDisorientedTriangles)
4820  : mPoints(nullptr)
4821  , mPolygons()
4822  , mPointListSize(0)
4823  , mSeamPointListSize(0)
4824  , mPolygonPoolListSize(0)
4825  , mIsovalue(isovalue)
4826  , mPrimAdaptivity(adaptivity)
4827  , mSecAdaptivity(0.0)
4828  , mRefGrid(GridBase::ConstPtr())
4829  , mSurfaceMaskGrid(GridBase::ConstPtr())
4830  , mAdaptivityGrid(GridBase::ConstPtr())
4831  , mAdaptivityMaskTree(TreeBase::ConstPtr())
4832  , mRefSignTree(TreeBase::Ptr())
4833  , mRefIdxTree(TreeBase::Ptr())
4834  , mInvertSurfaceMask(false)
4835  , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4836  , mQuantizedSeamPoints(nullptr)
4837  , mPointFlags(0)
4838 {
4839 }
4840 
4841 
4842 inline void
4843 VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity)
4844 {
4845  mRefGrid = grid;
4846  mSecAdaptivity = secAdaptivity;
4847 
4848  // Clear out old auxiliary data
4849  mRefSignTree = TreeBase::Ptr();
4850  mRefIdxTree = TreeBase::Ptr();
4851  mSeamPointListSize = 0;
4852  mQuantizedSeamPoints.reset(nullptr);
4853 }
4854 
4855 
4856 inline void
4858 {
4859  mSurfaceMaskGrid = mask;
4860  mInvertSurfaceMask = invertMask;
4861 }
4862 
4863 
4864 inline void
4866 {
4867  mAdaptivityGrid = grid;
4868 }
4869 
4870 
4871 inline void
4873 {
4874  mAdaptivityMaskTree = tree;
4875 }
4876 
4877 
4878 template<typename InputGridType>
4879 inline void
4880 VolumeToMesh::operator()(const InputGridType& inputGrid)
4881 {
4882  // input data types
4883 
4884  using InputTreeType = typename InputGridType::TreeType;
4885  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
4886  using InputValueType = typename InputLeafNodeType::ValueType;
4887 
4888  // auxiliary data types
4889 
4890  using FloatTreeType = typename InputTreeType::template ValueConverter<float>::Type;
4891  using FloatGridType = Grid<FloatTreeType>;
4892  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
4893  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
4894  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
4895  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
4896  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
4897 
4898  // clear old data
4899  mPointListSize = 0;
4900  mPoints.reset();
4901  mPolygonPoolListSize = 0;
4902  mPolygons.reset();
4903  mPointFlags.clear();
4904 
4905  // settings
4906 
4907  const math::Transform& transform = inputGrid.transform();
4908  const InputValueType isovalue = InputValueType(mIsovalue);
4909  const float adaptivityThreshold = float(mPrimAdaptivity);
4910  const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4911 
4912  // The default surface orientation is setup for level set and bool/mask grids.
4913  // Boolean grids are handled correctly by their value type. Signed distance fields,
4914  // unsigned distance fields and fog volumes have the same value type but use different
4915  // inside value classifications.
4916  const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4917  && (inputGrid.getGridClass() != openvdb::GRID_LEVEL_SET));
4918 
4919  // references, masks and auxiliary data
4920 
4921  const InputTreeType& inputTree = inputGrid.tree();
4922 
4923  BoolTreeType intersectionTree(false), adaptivityMask(false);
4924 
4925  if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4926  const BoolTreeType *refAdaptivityMask=
4927  static_cast<const BoolTreeType*>(mAdaptivityMaskTree.get());
4928  adaptivityMask.topologyUnion(*refAdaptivityMask);
4929  }
4930 
4931  Int16TreeType signFlagsTree(0);
4932  Index32TreeType pointIndexTree(std::numeric_limits<Index32>::max());
4933 
4934 
4935  // collect auxiliary data
4936 
4937  volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4938  intersectionTree, inputTree, isovalue);
4939 
4940  volume_to_mesh_internal::applySurfaceMask(intersectionTree, adaptivityMask,
4941  inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4942 
4943  if (intersectionTree.empty()) return;
4944 
4945  volume_to_mesh_internal::computeAuxiliaryData(
4946  signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4947 
4948  intersectionTree.clear();
4949 
4950  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4951  pointIndexTree.getNodes(pointIndexLeafNodes);
4952 
4953  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4954  signFlagsTree.getNodes(signFlagsLeafNodes);
4955 
4956  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4957 
4958 
4959  // optionally collect auxiliary data from a reference volume.
4960 
4961  Int16TreeType* refSignFlagsTree = nullptr;
4962  Index32TreeType* refPointIndexTree = nullptr;
4963  InputTreeType const* refInputTree = nullptr;
4964 
4965  if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4966 
4967  const InputGridType* refGrid = static_cast<const InputGridType*>(mRefGrid.get());
4968  refInputTree = &refGrid->tree();
4969 
4970  if (!mRefSignTree && !mRefIdxTree) {
4971 
4972  // first time, collect and cache auxiliary data.
4973 
4974  typename Int16TreeType::Ptr refSignFlagsTreePt(new Int16TreeType(0));
4975  typename Index32TreeType::Ptr refPointIndexTreePt(
4976  new Index32TreeType(std::numeric_limits<Index32>::max()));
4977 
4978  BoolTreeType refIntersectionTree(false);
4979 
4980  volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4981  refIntersectionTree, *refInputTree, isovalue);
4982 
4983  volume_to_mesh_internal::computeAuxiliaryData(*refSignFlagsTreePt,
4984  *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
4985 
4986  mRefSignTree = refSignFlagsTreePt;
4987  mRefIdxTree = refPointIndexTreePt;
4988  }
4989 
4990  if (mRefSignTree && mRefIdxTree) {
4991 
4992  // get cached auxiliary data
4993 
4994  refSignFlagsTree = static_cast<Int16TreeType*>(mRefSignTree.get());
4995  refPointIndexTree = static_cast<Index32TreeType*>(mRefIdxTree.get());
4996  }
4997 
4998 
4999  if (refSignFlagsTree && refPointIndexTree) {
5000 
5001  // generate seam line sample points
5002 
5003  volume_to_mesh_internal::markSeamLineData(signFlagsTree, *refSignFlagsTree);
5004 
5005  if (mSeamPointListSize == 0) {
5006 
5007  // count unique points on reference surface
5008 
5009  std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
5010  refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
5011 
5012  std::unique_ptr<Index32[]> leafNodeOffsets(
5013  new Index32[refSignFlagsLeafNodes.size()]);
5014 
5015  tbb::parallel_for(tbb::blocked_range<size_t>(0, refSignFlagsLeafNodes.size()),
5016  volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>(
5017  refSignFlagsLeafNodes, leafNodeOffsets));
5018 
5019  {
5020  Index32 count = 0;
5021  for (size_t n = 0, N = refSignFlagsLeafNodes.size(); n < N; ++n) {
5022  const Index32 tmp = leafNodeOffsets[n];
5023  leafNodeOffsets[n] = count;
5024  count += tmp;
5025  }
5026  mSeamPointListSize = size_t(count);
5027  }
5028 
5029  if (mSeamPointListSize != 0) {
5030 
5031  mQuantizedSeamPoints.reset(new uint32_t[mSeamPointListSize]);
5032 
5033  std::memset(mQuantizedSeamPoints.get(), 0, sizeof(uint32_t) * mSeamPointListSize);
5034 
5035  std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5036  refPointIndexTree->getNodes(refPointIndexLeafNodes);
5037 
5038  tbb::parallel_for(tbb::blocked_range<size_t>(0, refPointIndexLeafNodes.size()),
5039  volume_to_mesh_internal::MapPoints<Index32LeafNodeType>(
5040  refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5041  }
5042  }
5043 
5044  if (mSeamPointListSize != 0) {
5045 
5046  tbb::parallel_for(auxiliaryLeafNodeRange,
5047  volume_to_mesh_internal::SeamLineWeights<InputTreeType>(
5048  signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5049  mQuantizedSeamPoints.get(), isovalue));
5050  }
5051  }
5052  }
5053 
5054  const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5055 
5056 
5057  // adapt and count unique points
5058 
5059  std::unique_ptr<Index32[]> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
5060 
5061  if (adaptive) {
5062  volume_to_mesh_internal::MergeVoxelRegions<InputGridType> mergeOp(
5063  inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5064  isovalue, adaptivityThreshold, invertSurfaceOrientation);
5065 
5066  if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5067  const FloatGridType* adaptivityGrid =
5068  static_cast<const FloatGridType*>(mAdaptivityGrid.get());
5069  mergeOp.setSpatialAdaptivity(*adaptivityGrid);
5070  }
5071 
5072  if (!adaptivityMask.empty()) {
5073  mergeOp.setAdaptivityMask(adaptivityMask);
5074  }
5075 
5076  if (referenceMeshing) {
5077  mergeOp.setRefSignFlagsData(*refSignFlagsTree, float(mSecAdaptivity));
5078  }
5079 
5080  tbb::parallel_for(auxiliaryLeafNodeRange, mergeOp);
5081 
5082  volume_to_mesh_internal::AdaptiveLeafNodePointCount<Index32LeafNodeType>
5083  op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5084 
5085  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5086 
5087  } else {
5088 
5089  volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
5090  op(signFlagsLeafNodes, leafNodeOffsets);
5091 
5092  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5093  }
5094 
5095 
5096  {
5097  Index32 pointCount = 0;
5098  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
5099  const Index32 tmp = leafNodeOffsets[n];
5100  leafNodeOffsets[n] = pointCount;
5101  pointCount += tmp;
5102  }
5103 
5104  mPointListSize = size_t(pointCount);
5105  mPoints.reset(new openvdb::Vec3s[mPointListSize]);
5106  mPointFlags.clear();
5107  }
5108 
5109 
5110  // compute points
5111 
5112  {
5113  volume_to_mesh_internal::ComputePoints<InputTreeType>
5114  op(mPoints.get(), inputTree, pointIndexLeafNodes,
5115  signFlagsLeafNodes, leafNodeOffsets, transform, mIsovalue);
5116 
5117  if (referenceMeshing) {
5118  mPointFlags.resize(mPointListSize);
5119  op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5120  mQuantizedSeamPoints.get(), mPointFlags.data());
5121  }
5122 
5123  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5124  }
5125 
5126 
5127  // compute polygons
5128 
5129  mPolygonPoolListSize = signFlagsLeafNodes.size();
5130  mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
5131 
5132  if (adaptive) {
5133 
5134  using PrimBuilder = volume_to_mesh_internal::AdaptivePrimBuilder;
5135 
5136  volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5137  op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5138  mPolygons, invertSurfaceOrientation);
5139 
5140  if (referenceMeshing) {
5141  op.setRefSignTree(refSignFlagsTree);
5142  }
5143 
5144  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5145 
5146  } else {
5147 
5148  using PrimBuilder = volume_to_mesh_internal::UniformPrimBuilder;
5149 
5150  volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5151  op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5152  mPolygons, invertSurfaceOrientation);
5153 
5154  if (referenceMeshing) {
5155  op.setRefSignTree(refSignFlagsTree);
5156  }
5157 
5158  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5159  }
5160 
5161 
5162  signFlagsTree.clear();
5163  pointIndexTree.clear();
5164 
5165 
5166  if (adaptive && mRelaxDisorientedTriangles) {
5167  volume_to_mesh_internal::relaxDisorientedTriangles(invertSurfaceOrientation,
5168  inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5169  }
5170 
5171 
5172  if (referenceMeshing) {
5173  volume_to_mesh_internal::subdivideNonplanarSeamLineQuads(
5174  mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5175 
5176  volume_to_mesh_internal::reviseSeamLineFlags(mPolygons, mPolygonPoolListSize, mPointFlags);
5177  }
5178 
5179 }
5180 
5181 
5182 ////////////////////////////////////////
5183 
5184 
5185 template<typename GridType>
5187  const GridType& grid,
5188  std::vector<Vec3s>& points,
5189  std::vector<Vec3I>& triangles,
5190  std::vector<Vec4I>& quads,
5191  double isovalue,
5192  double adaptivity,
5193  bool relaxDisorientedTriangles)
5194 {
5195  volume_to_mesh_internal::doVolumeToMesh(grid, points, triangles, quads,
5196  isovalue, adaptivity, relaxDisorientedTriangles);
5197 }
5198 
5199 template<typename GridType>
5201  const GridType& grid,
5202  std::vector<Vec3s>& points,
5203  std::vector<Vec4I>& quads,
5204  double isovalue)
5205 {
5206  std::vector<Vec3I> triangles;
5207  volumeToMesh(grid, points, triangles, quads, isovalue, 0.0, true);
5208 }
5209 
5210 
5211 ////////////////////////////////////////
5212 
5213 
5214 // Explicit Template Instantiation
5215 
5216 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
5217 
5218 #ifdef OPENVDB_INSTANTIATE_VOLUMETOMESH
5220 #endif
5221 
5222 #define _FUNCTION(TreeT) \
5223  void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec4I>&, double)
5225 #undef _FUNCTION
5226 
5227 #define _FUNCTION(TreeT) \
5228  void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec3I>&, std::vector<Vec4I>&, double, double, bool)
5230 #undef _FUNCTION
5231 
5232 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
5233 
5234 
5235 } // namespace tools
5236 } // namespace OPENVDB_VERSION_NAME
5237 } // namespace openvdb
5238 
5239 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
ValueOnCIter cbeginValueOn() const
Definition: PointIndexGrid.h:1612
void setRefGrid(const GridBase::ConstPtr &grid, double secAdaptivity=0)
When surfacing fractured SDF fragments, the original unfractured SDF grid can be used to eliminate se...
Definition: VolumeToMesh.h:4843
void setValueOnly(const Coord &, const ValueType &)
Definition: PointIndexGrid.h:1509
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:38
Mat3< double > Mat3d
Definition: Mat3.h:834
Definition: VolumeToMesh.h:99
const size_t & numQuads() const
Definition: VolumeToMesh.h:121
void setAdaptivityMask(const TreeBase::ConstPtr &tree)
Definition: VolumeToMesh.h:4872
const PolygonPoolList & polygonPoolList() const
Definition: VolumeToMesh.h:190
Definition: PointIndexGrid.h:54
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
Mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:170
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
std::unique_ptr< openvdb::Vec3s[]> PointList
Point and primitive list types.
Definition: VolumeToMesh.h:161
int getValueDepth(const Coord &xyz) const
Return the tree depth (0 = root) at which the value of voxel (x, y, z) resides, or -1 if (x...
Definition: ValueAccessor.h:516
openvdb::Vec4I & quad(size_t n)
Definition: VolumeToMesh.h:123
#define OPENVDB_NUMERIC_TREE_INSTANTIATE(Function)
Definition: version.h.in:163
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: LeafNode.h:1075
void setSpatialAdaptivity(const GridBase::ConstPtr &grid)
Definition: VolumeToMesh.h:4865
int16_t Int16
Definition: Types.h:55
char & triangleFlags(size_t n)
Definition: VolumeToMesh.h:138
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:454
T & z()
Definition: Vec3.h:88
std::unique_ptr< PolygonPool[]> PolygonPoolList
Point and primitive list types.
Definition: VolumeToMesh.h:162
const openvdb::Vec3I & triangle(size_t n) const
Definition: VolumeToMesh.h:130
Index32 Index
Definition: Types.h:54
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
OutGridT XformOp & op
Definition: ValueTransformer.h:139
PolygonPool()
Definition: VolumeToMesh.h:4680
void copy(const PolygonPool &rhs)
Definition: VolumeToMesh.h:4704
Vec3d findFeaturePoint(const std::vector< Vec3d > &points, const std::vector< Vec3d > &normals)
Given a set of tangent elements, points with corresponding normals, this method returns the intersect...
Definition: VolumeToMesh.h:280
constexpr Index32 INVALID_IDX
Definition: Util.h:19
void clearQuads()
Definition: VolumeToMesh.h:4731
Definition: Exceptions.h:65
void volumeToMesh(const GridType &grid, std::vector< Vec3s > &points, std::vector< Vec3I > &triangles, std::vector< Vec4I > &quads, double isovalue=0.0, double adaptivity=0.0, bool relaxDisorientedTriangles=true)
Adaptively mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:5186
typename BaseLeaf::template ValueIter< MaskOnIterator, const PointIndexLeafNode, const ValueType, ValueOn > ValueOnCIter
Definition: PointIndexGrid.h:1589
char & quadFlags(size_t n)
Definition: VolumeToMesh.h:135
ValueOnCIter beginValueOn() const
Definition: PointIndexGrid.h:1613
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:363
VolumeToMesh(double isovalue=0, double adaptivity=0, bool relaxDisorientedTriangles=true)
Definition: VolumeToMesh.h:4819
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
openvdb::Vec3I & triangle(size_t n)
Definition: VolumeToMesh.h:129
Vec3< float > Vec3s
Definition: Vec3.h:664
bool trimTrinagles(const size_t n, bool reallocate=false)
Definition: VolumeToMesh.h:4787
Definition: Types.h:455
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
Definition: ValueAccessor.h:492
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
size_t polygonPoolListSize() const
Definition: VolumeToMesh.h:188
Base class for typed trees.
Definition: Tree.h:37
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Definition: PointCountImpl.h:18
SharedPtr< const TreeBase > ConstPtr
Definition: Tree.h:41
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
Definition: ValueAccessor.h:838
Definition: Exceptions.h:13
Collection of quads and triangles.
Definition: VolumeToMesh.h:103
void operator()(const InputGridType &)
Main call.
Definition: VolumeToMesh.h:4880
const std::vector< uint8_t > & pointFlags() const
Definition: VolumeToMesh.h:193
SharedPtr< const GridBase > ConstPtr
Definition: Grid.h:81
T ValueType
Definition: PointIndexGrid.h:1363
Definition: Mat.h:165
OutGridT const XformOp bool bool
Definition: ValueTransformer.h:609
void volumeToMesh(const GridType &grid, std::vector< Vec3s > &points, std::vector< Vec4I > &quads, double isovalue=0.0)
Uniformly mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:5200
Definition: VolumeToMesh.h:99
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
void clearTriangles()
Definition: VolumeToMesh.h:4749
const char & quadFlags(size_t n) const
Definition: VolumeToMesh.h:136
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:219
uint32_t Index32
Definition: Types.h:52
void resetTriangles(size_t size)
Definition: VolumeToMesh.h:4740
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:455
Definition: Transform.h:39
PolygonPoolList & polygonPoolList()
Definition: VolumeToMesh.h:189
const PointList & pointList() const
Definition: VolumeToMesh.h:186
void resetQuads(size_t size)
Definition: VolumeToMesh.h:4722
bool diagonalizeSymmetricMatrix(const Mat3< T > &input, Mat3< T > &Q, Vec3< T > &D, unsigned int MAX_ITERATIONS=250)
Use Jacobi iterations to decompose a symmetric 3x3 matrix (diagonalize and compute eigenvectors) ...
Definition: Mat3.h:737
Definition: Mat4.h:24
const size_t & numTriangles() const
Definition: VolumeToMesh.h:127
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
Vec4< int32_t > Vec4i
Definition: Vec4.h:559
const char & triangleFlags(size_t n) const
Definition: VolumeToMesh.h:139
SharedPtr< TreeBase > Ptr
Definition: Tree.h:40
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
const openvdb::Vec4I & quad(size_t n) const
Definition: VolumeToMesh.h:124
Vec3< int32_t > Vec3i
Definition: Vec3.h:662
std::vector< uint8_t > & pointFlags()
Definition: VolumeToMesh.h:192
bool trimQuads(const size_t n, bool reallocate=false)
Definition: VolumeToMesh.h:4758
PointList & pointList()
Definition: VolumeToMesh.h:185
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:192
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
void setValue(const Coord &xyz, const ValueType &value)
Set a particular value at the given coordinate and mark the coordinate as active. ...
Definition: ValueAccessor.h:550
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
Definition: ValueAccessor.h:480
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
size_t pointListSize() const
Definition: VolumeToMesh.h:184
Gradient operators defined in index space of various orders.
Definition: Operators.h:99
Abstract base class for typed grids.
Definition: Grid.h:77
void setSurfaceMask(const GridBase::ConstPtr &mask, bool invertMask=false)
Definition: VolumeToMesh.h:4857