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