HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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, weightcount;
1618 
1619  if (refSigns == 0) {
1620  count = computeCellPoints(points, values, signs, iso);
1621  weightcount = 0;
1622  } else {
1623  if (inclusiveCell && refInputNode) {
1624  getCellVertexValues(*refInputNode, offset, refValues);
1625  } else {
1626  getCellVertexValues(*refInputAcc, ijk, refValues);
1627  }
1628  count = computeCellPoints(points, weightedPointMask, values, refValues, signs, refSigns,
1629  iso, refPointIndexNode->getValue(offset), mQuantizedSeamLinePoints);
1630  weightcount = count;
1631  }
1632 
1633  xyz = ijk.asVec3d();
1634 
1635  for (size_t i = 0; i < count; ++i) {
1636 
1637  Vec3d& point = points[i];
1638 
1639  // Checks for both NaN and inf vertex positions, i.e. any value that is not finite.
1640  if (!std::isfinite(point[0]) ||
1641  !std::isfinite(point[1]) ||
1642  !std::isfinite(point[2]))
1643  {
1644  OPENVDB_THROW(ValueError,
1645  "VolumeToMesh encountered NaNs or infs in the input VDB!"
1646  " Hint: Check the input and consider using the \"Diagnostics\" tool "
1647  "to detect and resolve the NaNs.");
1648  }
1649 
1650  point += xyz;
1651  point = mTransform.indexToWorld(point);
1652 
1653  Vec3s& pos = mPoints[pointOffset];
1654  pos[0] = float(point[0]);
1655  pos[1] = float(point[1]);
1656  pos[2] = float(point[2]);
1657 
1658  if (mSeamLinePointsFlags && weightcount && weightedPointMask[i]) {
1659  mSeamLinePointsFlags[pointOffset] = uint8_t(1);
1660  }
1661 
1662  ++pointOffset;
1663  }
1664  }
1665 
1666  // generate collapsed region points
1667  for (auto it = regions.begin(); it != regions.end(); ++it)
1668  {
1669  Vec3d avg(0.0);
1670 
1671  const IndexArray& voxels = it->second;
1672  for (size_t i = 0, I = voxels.size(); i < I; ++i) {
1673 
1674  const Index offset = voxels[i];
1675  ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1676 
1677  const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1678 
1679  ijk += origin;
1680 
1681  pidxData[offset] = pointOffset;
1682 
1683  const uint8_t signs = uint8_t(SIGNS & sfData[offset]);
1684 
1685  if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1686  else getCellVertexValues(inputAcc, ijk, values);
1687 
1688  computeCellPoints(points, values, signs, iso);
1689 
1690  avg[0] += double(ijk[0]) + points[0][0];
1691  avg[1] += double(ijk[1]) + points[0][1];
1692  avg[2] += double(ijk[2]) + points[0][2];
1693  }
1694 
1695  if (voxels.size() > 1) {
1696  const double w = 1.0 / double(voxels.size());
1697  avg *= w;
1698  }
1699 
1700  avg = mTransform.indexToWorld(avg);
1701 
1702  Vec3s& pos = mPoints[pointOffset];
1703  pos[0] = float(avg[0]);
1704  pos[1] = float(avg[1]);
1705  pos[2] = float(avg[2]);
1706 
1707  ++pointOffset;
1708  }
1709  }
1710 } // ComputePoints::operator()
1711 
1712 
1713 ////////////////////////////////////////
1714 
1715 
1716 template <typename InputTreeType>
1717 struct SeamLineWeights
1718 {
1719  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
1720  using InputValueType = typename InputLeafNodeType::ValueType;
1721 
1722  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
1723  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
1724 
1725  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
1726  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
1727 
1728  SeamLineWeights(const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
1729  const InputTreeType& inputTree,
1730  const Index32TreeType& refPointIndexTree,
1731  const Int16TreeType& refSignFlagsTree,
1732  uint32_t * quantizedPoints,
1733  InputValueType iso)
1734  : mSignFlagsNodes(signFlagsLeafNodes.data())
1735  , mInputTree(&inputTree)
1736  , mRefPointIndexTree(&refPointIndexTree)
1737  , mRefSignFlagsTree(&refSignFlagsTree)
1738  , mQuantizedPoints(quantizedPoints)
1739  , mIsovalue(iso)
1740  {
1741  }
1742 
1743  void operator()(const tbb::blocked_range<size_t>& range) const
1744  {
1745  tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
1746  tree::ValueAccessor<const Index32TreeType> pointIndexTreeAcc(*mRefPointIndexTree);
1747  tree::ValueAccessor<const Int16TreeType> signFlagsTreeAcc(*mRefSignFlagsTree);
1748 
1749  std::array<double, 8> values;
1750  const double iso = double(mIsovalue);
1751  Coord ijk;
1752  Vec3d pos;
1753 
1754  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1755 
1756  const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1757  const Coord& origin = signFlagsNode.origin();
1758 
1759  const Int16LeafNodeType * refSignNode = signFlagsTreeAcc.probeConstLeaf(origin);
1760  if (!refSignNode) continue;
1761 
1762  const Index32LeafNodeType* refPointIndexNode =
1763  pointIndexTreeAcc.probeConstLeaf(origin);
1764  if (!refPointIndexNode) continue;
1765 
1766  const InputLeafNodeType* inputNode = inputTreeAcc.probeConstLeaf(origin);
1767 
1768  const auto* const sfData = signFlagsNode.buffer().data();
1769  const auto* const rfIdxData = refPointIndexNode->buffer().data();
1770  const auto* const rsData = refSignNode->buffer().data();
1771 
1772  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it)
1773  {
1774  const Index offset = it.pos();
1775  const Int16 flags = sfData[offset];
1776 
1777  ijk = Index32LeafNodeType::offsetToLocalCoord(offset);
1778 
1779  const bool inclusiveCell = inputNode && isInternalLeafCoord<InputLeafNodeType>(ijk);
1780 
1781  ijk += origin;
1782 
1783  if ((flags & SEAM) && refSignNode->isValueOn(offset)) {
1784 
1785  const uint8_t lhsSigns = uint8_t(SIGNS & flags);
1786  const uint8_t rhsSigns = uint8_t(SIGNS & rsData[offset]);
1787 
1788  if (inclusiveCell) getCellVertexValues(*inputNode, offset, values);
1789  else getCellVertexValues(inputTreeAcc, ijk, values);
1790 
1791  for (unsigned i = 1, I = sEdgeGroupTable[lhsSigns][0] + 1; i < I; ++i) {
1792 
1793  const int id = matchEdgeGroup(uint8_t(i), lhsSigns, rhsSigns);
1794 
1795  if (id != -1) {
1796 
1797  uint32_t& data = mQuantizedPoints[rfIdxData[offset] + (id - 1)];
1798 
1799  if (!(data & MASK_DIRTY_BIT)) {
1800 
1801  const int samples = computeMaskedPoint(
1802  pos, values, lhsSigns, rhsSigns, uint8_t(i), iso);
1803 
1804  if (samples > 0) data = packPoint(pos);
1805  else data = MASK_INVALID_BIT;
1806 
1807  data |= MASK_DIRTY_BIT;
1808  }
1809  }
1810  } // end point group loop
1811  }
1812  } // end value on loop
1813  } // end leaf node loop
1814  }
1815 
1816 private:
1817  Int16LeafNodeType const * const * const mSignFlagsNodes;
1818  InputTreeType const * const mInputTree;
1819  Index32TreeType const * const mRefPointIndexTree;
1820  Int16TreeType const * const mRefSignFlagsTree;
1821  uint32_t * const mQuantizedPoints;
1822  InputValueType const mIsovalue;
1823 }; // struct SeamLineWeights
1824 
1825 
1826 template <typename TreeType>
1827 struct SetSeamLineFlags
1828 {
1829  using LeafNodeType = typename TreeType::LeafNodeType;
1830 
1831  SetSeamLineFlags(const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1832  const TreeType& refSignFlagsTree)
1833  : mSignFlagsNodes(signFlagsLeafNodes.data())
1834  , mRefSignFlagsTree(&refSignFlagsTree)
1835  {
1836  }
1837 
1838  void operator()(const tbb::blocked_range<size_t>& range) const
1839  {
1840  tree::ValueAccessor<const TreeType> refSignFlagsTreeAcc(*mRefSignFlagsTree);
1841 
1842  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1843 
1844  LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1845  const Coord& origin = signFlagsNode.origin();
1846 
1847  const LeafNodeType* refSignNode = refSignFlagsTreeAcc.probeConstLeaf(origin);
1848  if (!refSignNode) continue;
1849 
1850  const auto* const rsData = refSignNode->buffer().data();
1851  auto* const sfData = signFlagsNode.buffer().data();
1852 
1853  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1854  const Index offset = it.pos();
1855 
1856  const uint8_t rhsSigns = uint8_t(rsData[offset] & SIGNS);
1857 
1858  if (sEdgeGroupTable[rhsSigns][0] > 0) {
1859 
1860  typename LeafNodeType::ValueType& value = sfData[offset];
1861  const uint8_t lhsSigns = uint8_t(value & SIGNS);
1862 
1863  if (rhsSigns != lhsSigns) {
1864  value |= SEAM;
1865  }
1866  }
1867 
1868  } // end value on loop
1869 
1870  } // end leaf node loop
1871  }
1872 
1873 private:
1874  LeafNodeType * const * const mSignFlagsNodes;
1875  TreeType const * const mRefSignFlagsTree;
1876 }; // struct SetSeamLineFlags
1877 
1878 
1879 template <typename BoolTreeType, typename SignDataType>
1880 struct TransferSeamLineFlags
1881 {
1882  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1883 
1884  using SignDataTreeType = typename BoolTreeType::template ValueConverter<SignDataType>::Type;
1885  using SignDataLeafNodeType = typename SignDataTreeType::LeafNodeType;
1886 
1887  TransferSeamLineFlags(const std::vector<SignDataLeafNodeType*>& signFlagsLeafNodes,
1888  const BoolTreeType& maskTree)
1889  : mSignFlagsNodes(signFlagsLeafNodes.data())
1890  , mMaskTree(&maskTree)
1891  {
1892  }
1893 
1894  void operator()(const tbb::blocked_range<size_t>& range) const
1895  {
1896  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1897 
1898  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1899 
1900  SignDataLeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1901  const Coord& origin = signFlagsNode.origin();
1902 
1903  const BoolLeafNodeType * maskNode = maskAcc.probeConstLeaf(origin);
1904  if (!maskNode) continue;
1905 
1906  auto* const sfData = signFlagsNode.buffer().data();
1907 
1908  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
1909  const Index offset = it.pos();
1910 
1911  if (maskNode->isValueOn(offset)) {
1912  sfData[offset] |= SEAM;
1913  }
1914  } // end value on loop
1915  } // end leaf node loop
1916  }
1917 
1918 private:
1919  SignDataLeafNodeType * const * const mSignFlagsNodes;
1920  BoolTreeType const * const mMaskTree;
1921 }; // struct TransferSeamLineFlags
1922 
1923 
1924 template <typename TreeType>
1925 struct MaskSeamLineVoxels
1926 {
1927  using LeafNodeType = typename TreeType::LeafNodeType;
1928  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1929 
1930  MaskSeamLineVoxels(const std::vector<LeafNodeType*>& signFlagsLeafNodes,
1931  const TreeType& signFlagsTree,
1932  BoolTreeType& mask)
1933  : mSignFlagsNodes(signFlagsLeafNodes.data())
1934  , mSignFlagsTree(&signFlagsTree)
1935  , mTempMask(false)
1936  , mMask(&mask)
1937  {
1938  }
1939 
1940  MaskSeamLineVoxels(MaskSeamLineVoxels& rhs, tbb::split)
1941  : mSignFlagsNodes(rhs.mSignFlagsNodes)
1942  , mSignFlagsTree(rhs.mSignFlagsTree)
1943  , mTempMask(false)
1944  , mMask(&mTempMask)
1945  {
1946  }
1947 
1948  void join(MaskSeamLineVoxels& rhs) { mMask->merge(*rhs.mMask); }
1949 
1950  void operator()(const tbb::blocked_range<size_t>& range)
1951  {
1952  using ValueOnCIter = typename LeafNodeType::ValueOnCIter;
1953  using ValueType = typename LeafNodeType::ValueType;
1954 
1955  tree::ValueAccessor<const TreeType> signFlagsAcc(*mSignFlagsTree);
1956  tree::ValueAccessor<BoolTreeType> maskAcc(*mMask);
1957  Coord ijk;
1958 
1959  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
1960 
1961  LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
1962  auto* const sfData = signFlagsNode.buffer().data();
1963 
1964  for (ValueOnCIter it = signFlagsNode.cbeginValueOn(); it; ++it) {
1965 
1966  const ValueType flags = sfData[it.pos()];
1967 
1968  if (!(flags & SEAM) && (flags & EDGES)) {
1969 
1970  ijk = it.getCoord();
1971 
1972  bool isSeamLineVoxel = false;
1973 
1974  if (flags & XEDGE) {
1975  ijk[1] -= 1;
1976  isSeamLineVoxel = (signFlagsAcc.getValue(ijk) & SEAM);
1977  ijk[2] -= 1;
1978  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1979  ijk[1] += 1;
1980  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1981  ijk[2] += 1;
1982  }
1983 
1984  if (!isSeamLineVoxel && flags & YEDGE) {
1985  ijk[2] -= 1;
1986  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1987  ijk[0] -= 1;
1988  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1989  ijk[2] += 1;
1990  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1991  ijk[0] += 1;
1992  }
1993 
1994  if (!isSeamLineVoxel && flags & ZEDGE) {
1995  ijk[1] -= 1;
1996  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1997  ijk[0] -= 1;
1998  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
1999  ijk[1] += 1;
2000  isSeamLineVoxel = isSeamLineVoxel || (signFlagsAcc.getValue(ijk) & SEAM);
2001  ijk[0] += 1;
2002  }
2003 
2004  if (isSeamLineVoxel) {
2005  maskAcc.setValue(it.getCoord(), true);
2006  }
2007  }
2008  } // end value on loop
2009 
2010  } // end leaf node loop
2011  }
2012 
2013 private:
2014  LeafNodeType * const * const mSignFlagsNodes;
2015  TreeType const * const mSignFlagsTree;
2016  BoolTreeType mTempMask;
2017  BoolTreeType * const mMask;
2018 }; // struct MaskSeamLineVoxels
2019 
2020 
2021 template<typename SignDataTreeType>
2022 inline void
2023 markSeamLineData(SignDataTreeType& signFlagsTree, const SignDataTreeType& refSignFlagsTree)
2024 {
2025  using SignDataType = typename SignDataTreeType::ValueType;
2026  using SignDataLeafNodeType = typename SignDataTreeType::LeafNodeType;
2027  using BoolTreeType = typename SignDataTreeType::template ValueConverter<bool>::Type;
2028 
2029  std::vector<SignDataLeafNodeType*> signFlagsLeafNodes;
2030  signFlagsTree.getNodes(signFlagsLeafNodes);
2031 
2032  const tbb::blocked_range<size_t> nodeRange(0, signFlagsLeafNodes.size());
2033 
2034  tbb::parallel_for(nodeRange,
2035  SetSeamLineFlags<SignDataTreeType>(signFlagsLeafNodes, refSignFlagsTree));
2036 
2037  BoolTreeType seamLineMaskTree(false);
2038 
2039  MaskSeamLineVoxels<SignDataTreeType>
2040  maskSeamLine(signFlagsLeafNodes, signFlagsTree, seamLineMaskTree);
2041 
2042  tbb::parallel_reduce(nodeRange, maskSeamLine);
2043 
2044  tbb::parallel_for(nodeRange,
2045  TransferSeamLineFlags<BoolTreeType, SignDataType>(signFlagsLeafNodes, seamLineMaskTree));
2046 }
2047 
2048 
2049 ////////////////////////////////////////
2050 
2051 
2052 template <typename InputGridType>
2053 struct MergeVoxelRegions
2054 {
2055  using InputTreeType = typename InputGridType::TreeType;
2056  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
2057  using InputValueType = typename InputLeafNodeType::ValueType;
2058 
2059  using FloatTreeType = typename InputTreeType::template ValueConverter<float>::Type;
2060  using FloatLeafNodeType = typename FloatTreeType::LeafNodeType;
2061  using FloatGridType = Grid<FloatTreeType>;
2062 
2063  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
2064  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
2065 
2066  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
2067  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
2068 
2069  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
2070  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2071 
2072  using MaskTreeType = typename InputTreeType::template ValueConverter<ValueMask>::Type;
2073  using MaskLeafNodeType = typename MaskTreeType::LeafNodeType;
2074 
2075  MergeVoxelRegions(const InputGridType& inputGrid,
2076  const Index32TreeType& pointIndexTree,
2077  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2078  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2079  InputValueType iso,
2080  float adaptivity,
2081  bool invertSurfaceOrientation);
2082 
2083  void setSpatialAdaptivity(const FloatGridType& grid)
2084  {
2085  mSpatialAdaptivityTree = &grid.tree();
2086  mSpatialAdaptivityTransform = &grid.transform();
2087  }
2088 
2089  void setAdaptivityMask(const BoolTreeType& mask)
2090  {
2091  mMaskTree = &mask;
2092  }
2093 
2094  void setRefSignFlagsData(const Int16TreeType& signFlagsData, float internalAdaptivity)
2095  {
2096  mRefSignFlagsTree = &signFlagsData;
2097  mInternalAdaptivity = internalAdaptivity;
2098  }
2099 
2100  void operator()(const tbb::blocked_range<size_t>&) const;
2101 
2102 private:
2103  InputTreeType const * const mInputTree;
2104  math::Transform const * const mInputTransform;
2105 
2106  Index32TreeType const * const mPointIndexTree;
2107  Index32LeafNodeType * const * const mPointIndexNodes;
2108  Int16LeafNodeType const * const * const mSignFlagsNodes;
2109 
2110  InputValueType mIsovalue;
2111  float mSurfaceAdaptivity, mInternalAdaptivity;
2112  bool mInvertSurfaceOrientation;
2113 
2114  FloatTreeType const * mSpatialAdaptivityTree;
2115  BoolTreeType const * mMaskTree;
2116  Int16TreeType const * mRefSignFlagsTree;
2117  math::Transform const * mSpatialAdaptivityTransform;
2118 }; // struct MergeVoxelRegions
2119 
2120 
2121 template <typename InputGridType>
2122 MergeVoxelRegions<InputGridType>::MergeVoxelRegions(
2123  const InputGridType& inputGrid,
2124  const Index32TreeType& pointIndexTree,
2125  const std::vector<Index32LeafNodeType*>& pointIndexLeafNodes,
2126  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
2127  InputValueType iso,
2128  float adaptivity,
2129  bool invertSurfaceOrientation)
2130  : mInputTree(&inputGrid.tree())
2131  , mInputTransform(&inputGrid.transform())
2132  , mPointIndexTree(&pointIndexTree)
2133  , mPointIndexNodes(pointIndexLeafNodes.data())
2134  , mSignFlagsNodes(signFlagsLeafNodes.data())
2135  , mIsovalue(iso)
2136  , mSurfaceAdaptivity(adaptivity)
2137  , mInternalAdaptivity(adaptivity)
2138  , mInvertSurfaceOrientation(invertSurfaceOrientation)
2139  , mSpatialAdaptivityTree(nullptr)
2140  , mMaskTree(nullptr)
2141  , mRefSignFlagsTree(nullptr)
2142  , mSpatialAdaptivityTransform(nullptr)
2143 {
2144 }
2145 
2146 
2147 template <typename InputGridType>
2148 void
2149 MergeVoxelRegions<InputGridType>::operator()(const tbb::blocked_range<size_t>& range) const
2150 {
2151  using Vec3sType = math::Vec3<float>;
2152  using Vec3sLeafNodeType = typename InputLeafNodeType::template ValueConverter<Vec3sType>::Type;
2153 
2154  using InputTreeAccessor = tree::ValueAccessor<const InputTreeType>;
2155  using FloatTreeAccessor = tree::ValueAccessor<const FloatTreeType>;
2156  using Index32TreeAccessor = tree::ValueAccessor<const Index32TreeType>;
2157  using Int16TreeAccessor = tree::ValueAccessor<const Int16TreeType>;
2158  using BoolTreeAccessor = tree::ValueAccessor<const BoolTreeType>;
2159 
2160  std::unique_ptr<FloatTreeAccessor> spatialAdaptivityAcc;
2161  if (mSpatialAdaptivityTree && mSpatialAdaptivityTransform) {
2162  spatialAdaptivityAcc.reset(new FloatTreeAccessor(*mSpatialAdaptivityTree));
2163  }
2164 
2165  std::unique_ptr<BoolTreeAccessor> maskAcc;
2166  if (mMaskTree) {
2167  maskAcc.reset(new BoolTreeAccessor(*mMaskTree));
2168  }
2169 
2170  std::unique_ptr<Int16TreeAccessor> refSignFlagsAcc;
2171  if (mRefSignFlagsTree) {
2172  refSignFlagsAcc.reset(new Int16TreeAccessor(*mRefSignFlagsTree));
2173  }
2174 
2175  InputTreeAccessor inputAcc(*mInputTree);
2176  Index32TreeAccessor pointIndexAcc(*mPointIndexTree);
2177 
2178  MaskLeafNodeType mask;
2179 
2180  const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<InputValueType>();
2181  std::unique_ptr<Vec3sLeafNodeType> gradientNode;
2182 
2183  Coord ijk, end;
2184  const int LeafDim = InputLeafNodeType::DIM;
2185 
2186  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
2187 
2188  mask.setValuesOff();
2189 
2190  const Int16LeafNodeType& signFlagsNode = *mSignFlagsNodes[n];
2191  Index32LeafNodeType& pointIndexNode = *mPointIndexNodes[n];
2192 
2193  const Coord& origin = pointIndexNode.origin();
2194  end = origin.offsetBy(LeafDim);
2195 
2196  // Mask off seam line adjacent voxels
2197  if (maskAcc) {
2198  const BoolLeafNodeType* maskLeaf = maskAcc->probeConstLeaf(origin);
2199  if (maskLeaf != nullptr) {
2200  for (auto it = maskLeaf->cbeginValueOn(); it; ++it)
2201  {
2202  mask.setActiveState(it.getCoord() & ~1u, true);
2203  }
2204  }
2205  }
2206 
2207  float adaptivity = (refSignFlagsAcc && !refSignFlagsAcc->probeConstLeaf(origin)) ?
2208  mInternalAdaptivity : mSurfaceAdaptivity;
2209 
2210  bool useGradients = adaptivity < 1.0f;
2211 
2212  // Set region adaptivity
2213  FloatLeafNodeType adaptivityLeaf(origin, adaptivity);
2214 
2215  if (spatialAdaptivityAcc) {
2216  useGradients = false;
2217  for (Index offset = 0; offset < FloatLeafNodeType::NUM_VALUES; ++offset) {
2218  ijk = adaptivityLeaf.offsetToGlobalCoord(offset);
2219  ijk = mSpatialAdaptivityTransform->worldToIndexCellCentered(
2220  mInputTransform->indexToWorld(ijk));
2221  float weight = spatialAdaptivityAcc->getValue(ijk);
2222  float adaptivityValue = weight * adaptivity;
2223  if (adaptivityValue < 1.0f) useGradients = true;
2224  adaptivityLeaf.setValueOnly(offset, adaptivityValue);
2225  }
2226  }
2227 
2228  // Mask off ambiguous voxels
2229  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it) {
2230  const Int16 flags = it.getValue();
2231  const unsigned char signs = static_cast<unsigned char>(SIGNS & int(flags));
2232 
2233  if ((flags & SEAM) || !sAdaptable[signs] || sEdgeGroupTable[signs][0] > 1) {
2234 
2235  mask.setActiveState(it.getCoord() & ~1u, true);
2236 
2237  } else if (flags & EDGES) {
2238 
2239  bool maskRegion = false;
2240 
2241  ijk = it.getCoord();
2242  if (!pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2243 
2244  if (!maskRegion && flags & XEDGE) {
2245  ijk[1] -= 1;
2246  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2247  ijk[2] -= 1;
2248  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2249  ijk[1] += 1;
2250  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2251  ijk[2] += 1;
2252  }
2253 
2254  if (!maskRegion && flags & YEDGE) {
2255  ijk[2] -= 1;
2256  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2257  ijk[0] -= 1;
2258  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2259  ijk[2] += 1;
2260  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2261  ijk[0] += 1;
2262  }
2263 
2264  if (!maskRegion && flags & ZEDGE) {
2265  ijk[1] -= 1;
2266  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2267  ijk[0] -= 1;
2268  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2269  ijk[1] += 1;
2270  if (!maskRegion && !pointIndexAcc.isValueOn(ijk)) maskRegion = true;
2271  ijk[0] += 1;
2272  }
2273 
2274  if (maskRegion) {
2275  mask.setActiveState(it.getCoord() & ~1u, true);
2276  }
2277  }
2278  }
2279 
2280  // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
2281  int dim = 2;
2282  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2283  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2284  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2285  if (!mask.isValueOn(ijk) && isNonManifold(inputAcc, ijk, mIsovalue, dim)) {
2286  mask.setActiveState(ijk, true);
2287  }
2288  }
2289  }
2290  }
2291 
2292  // Compute the gradient for the remaining voxels
2293 
2294  if (useGradients) {
2295 
2296  if (gradientNode) {
2297  gradientNode->setValuesOff();
2298  } else {
2299  gradientNode.reset(new Vec3sLeafNodeType());
2300  }
2301 
2302  for (auto it = signFlagsNode.cbeginValueOn(); it; ++it)
2303  {
2304  ijk = it.getCoord();
2305  if (!mask.isValueOn(ijk & ~1u))
2306  {
2307  Vec3sType dir(math::ISGradient<math::CD_2ND>::result(inputAcc, ijk));
2308  dir.normalize();
2309 
2310  if (invertGradientDir) {
2311  dir = -dir;
2312  }
2313 
2314  gradientNode->setValueOn(it.pos(), dir);
2315  }
2316  }
2317  }
2318 
2319  // Merge regions
2320  int regionId = 1;
2321  for ( ; dim <= LeafDim; dim = dim << 1) {
2322  const unsigned coordMask = ~((dim << 1) - 1);
2323  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
2324  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
2325  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
2326 
2327  adaptivity = adaptivityLeaf.getValue(ijk);
2328 
2329  if (mask.isValueOn(ijk)
2330  || isNonManifold(inputAcc, ijk, mIsovalue, dim)
2331  || (useGradients && !isMergeable(*gradientNode, ijk, dim, adaptivity)))
2332  {
2333  mask.setActiveState(ijk & coordMask, true);
2334  } else {
2335  mergeVoxels(pointIndexNode, ijk, dim, regionId++);
2336  }
2337  }
2338  }
2339  }
2340  }
2341  }
2342 } // MergeVoxelRegions::operator()
2343 
2344 
2345 ////////////////////////////////////////
2346 
2347 
2348 // Constructs qudas
2349 struct UniformPrimBuilder
2350 {
2351  UniformPrimBuilder(): mIdx(0), mPolygonPool(nullptr) {}
2352 
2353  void init(const size_t upperBound, PolygonPool& quadPool)
2354  {
2355  mPolygonPool = &quadPool;
2356  mPolygonPool->resetQuads(upperBound);
2357  mIdx = 0;
2358  }
2359 
2360  template<typename IndexType>
2361  void addPrim(const math::Vec4<IndexType>& verts, bool reverse, char flags = 0)
2362  {
2363  if (!reverse) {
2364  mPolygonPool->quad(mIdx) = verts;
2365  } else {
2366  Vec4I& quad = mPolygonPool->quad(mIdx);
2367  quad[0] = verts[3];
2368  quad[1] = verts[2];
2369  quad[2] = verts[1];
2370  quad[3] = verts[0];
2371  }
2372  mPolygonPool->quadFlags(mIdx) = flags;
2373  ++mIdx;
2374  }
2375 
2376  void done()
2377  {
2378  mPolygonPool->trimQuads(mIdx);
2379  }
2380 
2381 private:
2382  size_t mIdx;
2383  PolygonPool* mPolygonPool;
2384 };
2385 
2386 
2387 // Constructs qudas and triangles
2388 struct AdaptivePrimBuilder
2389 {
2390  AdaptivePrimBuilder() : mQuadIdx(0), mTriangleIdx(0), mPolygonPool(nullptr) {}
2391 
2392  void init(const size_t upperBound, PolygonPool& polygonPool)
2393  {
2394  mPolygonPool = &polygonPool;
2395  mPolygonPool->resetQuads(upperBound);
2396  mPolygonPool->resetTriangles(upperBound);
2397 
2398  mQuadIdx = 0;
2399  mTriangleIdx = 0;
2400  }
2401 
2402  template<typename IndexType>
2403  void addPrim(const math::Vec4<IndexType>& verts, bool reverse, char flags = 0)
2404  {
2405  if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
2406  && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
2407  mPolygonPool->quadFlags(mQuadIdx) = flags;
2408  addQuad(verts, reverse);
2409  } else if (
2410  verts[0] == verts[3] &&
2411  verts[1] != verts[2] &&
2412  verts[1] != verts[0] &&
2413  verts[2] != verts[0]) {
2414  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2415  addTriangle(verts[0], verts[1], verts[2], reverse);
2416  } else if (
2417  verts[1] == verts[2] &&
2418  verts[0] != verts[3] &&
2419  verts[0] != verts[1] &&
2420  verts[3] != verts[1]) {
2421  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2422  addTriangle(verts[0], verts[1], verts[3], reverse);
2423  } else if (
2424  verts[0] == verts[1] &&
2425  verts[2] != verts[3] &&
2426  verts[2] != verts[0] &&
2427  verts[3] != verts[0]) {
2428  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2429  addTriangle(verts[0], verts[2], verts[3], reverse);
2430  } else if (
2431  verts[2] == verts[3] &&
2432  verts[0] != verts[1] &&
2433  verts[0] != verts[2] &&
2434  verts[1] != verts[2]) {
2435  mPolygonPool->triangleFlags(mTriangleIdx) = flags;
2436  addTriangle(verts[0], verts[1], verts[2], reverse);
2437  }
2438  }
2439 
2440 
2441  void done()
2442  {
2443  mPolygonPool->trimQuads(mQuadIdx, /*reallocate=*/true);
2444  mPolygonPool->trimTrinagles(mTriangleIdx, /*reallocate=*/true);
2445  }
2446 
2447 private:
2448 
2449  template<typename IndexType>
2450  void addQuad(const math::Vec4<IndexType>& verts, bool reverse)
2451  {
2452  if (!reverse) {
2453  mPolygonPool->quad(mQuadIdx) = verts;
2454  } else {
2455  Vec4I& quad = mPolygonPool->quad(mQuadIdx);
2456  quad[0] = verts[3];
2457  quad[1] = verts[2];
2458  quad[2] = verts[1];
2459  quad[3] = verts[0];
2460  }
2461  ++mQuadIdx;
2462  }
2463 
2464  void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
2465  {
2466  Vec3I& prim = mPolygonPool->triangle(mTriangleIdx);
2467 
2468  prim[1] = v1;
2469 
2470  if (!reverse) {
2471  prim[0] = v0;
2472  prim[2] = v2;
2473  } else {
2474  prim[0] = v2;
2475  prim[2] = v0;
2476  }
2477  ++mTriangleIdx;
2478  }
2479 
2480  size_t mQuadIdx, mTriangleIdx;
2481  PolygonPool *mPolygonPool;
2482 };
2483 
2484 
2485 template<typename SignAccT, typename IdxAccT, typename PrimBuilder>
2486 inline void
2487 constructPolygons(
2488  bool invertSurfaceOrientation,
2489  Int16 flags,
2490  Int16 refFlags,
2491  const Vec3i& offsets,
2492  const Coord& ijk,
2493  const SignAccT& signAcc,
2494  const IdxAccT& idxAcc,
2495  PrimBuilder& mesher)
2496 {
2497  using IndexType = typename IdxAccT::ValueType;
2498 
2500  const bool isActive = idxAcc.probeValue(ijk, v0);
2501  if (isActive == false || v0 == IndexType(util::INVALID_IDX)) return;
2502 
2503  char tag[2];
2504  tag[0] = (flags & SEAM) ? POLYFLAG_FRACTURE_SEAM : 0;
2505  tag[1] = tag[0] | char(POLYFLAG_EXTERIOR);
2506 
2507  bool isInside = flags & INSIDE;
2508 
2509  isInside = invertSurfaceOrientation ? !isInside : isInside;
2510 
2511  Coord coord = ijk;
2512  math::Vec4<IndexType> quad(0,0,0,0);
2513 
2514  if (flags & XEDGE) {
2515 
2516  quad[0] = v0 + offsets[0];
2517 
2518  // i, j-1, k
2519  coord[1]--;
2520  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2521  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2522  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][5] - 1 : 0;
2523 
2524  // i, j-1, k-1
2525  coord[2]--;
2526  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2527  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2528  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][7] - 1 : 0;
2529 
2530  // i, j, k-1
2531  coord[1]++;
2532  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2533  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2534  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][3] - 1 : 0;
2535 
2536  if (activeValues) {
2537  mesher.addPrim(quad, isInside, tag[bool(refFlags & XEDGE)]);
2538  }
2539 
2540  coord[2]++; // i, j, k
2541  }
2542 
2543 
2544  if (flags & YEDGE) {
2545 
2546  quad[0] = v0 + offsets[1];
2547 
2548  // i, j, k-1
2549  coord[2]--;
2550  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2551  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2552  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][12] - 1 : 0;
2553 
2554  // i-1, j, k-1
2555  coord[0]--;
2556  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2557  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2558  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][11] - 1 : 0;
2559 
2560  // i-1, j, k
2561  coord[2]++;
2562  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2563  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2564  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][10] - 1 : 0;
2565 
2566  if (activeValues) {
2567  mesher.addPrim(quad, isInside, tag[bool(refFlags & YEDGE)]);
2568  }
2569 
2570  coord[0]++; // i, j, k
2571  }
2572 
2573 
2574  if (flags & ZEDGE) {
2575 
2576  quad[0] = v0 + offsets[2];
2577 
2578  // i, j-1, k
2579  coord[1]--;
2580  bool activeValues = idxAcc.probeValue(coord, quad[1]);
2581  uint8_t cell = uint8_t(SIGNS & signAcc.getValue(coord));
2582  quad[1] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][8] - 1 : 0;
2583 
2584  // i-1, j-1, k
2585  coord[0]--;
2586  activeValues = activeValues && idxAcc.probeValue(coord, quad[2]);
2587  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2588  quad[2] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][6] - 1 : 0;
2589 
2590  // i-1, j, k
2591  coord[1]++;
2592  activeValues = activeValues && idxAcc.probeValue(coord, quad[3]);
2593  cell = uint8_t(SIGNS & signAcc.getValue(coord));
2594  quad[3] += sEdgeGroupTable[cell][0] > 1 ? sEdgeGroupTable[cell][2] - 1 : 0;
2595 
2596  if (activeValues) {
2597  mesher.addPrim(quad, !isInside, tag[bool(refFlags & ZEDGE)]);
2598  }
2599  }
2600 }
2601 
2602 
2603 ////////////////////////////////////////
2604 
2605 
2606 template<typename InputTreeType>
2607 struct MaskTileBorders
2608 {
2609  using InputValueType = typename InputTreeType::ValueType;
2610  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
2611 
2612 
2613  MaskTileBorders(const InputTreeType& inputTree, InputValueType iso,
2614  BoolTreeType& mask, const Vec4i* tileArray)
2615  : mInputTree(&inputTree)
2616  , mIsovalue(iso)
2617  , mTempMask(false)
2618  , mMask(&mask)
2619  , mTileArray(tileArray)
2620  {
2621  }
2622 
2623  MaskTileBorders(MaskTileBorders& rhs, tbb::split)
2624  : mInputTree(rhs.mInputTree)
2625  , mIsovalue(rhs.mIsovalue)
2626  , mTempMask(false)
2627  , mMask(&mTempMask)
2628  , mTileArray(rhs.mTileArray)
2629  {
2630  }
2631 
2632  void join(MaskTileBorders& rhs) { mMask->merge(*rhs.mMask); }
2633 
2634  void operator()(const tbb::blocked_range<size_t>&);
2635 
2636 private:
2637  InputTreeType const * const mInputTree;
2638  InputValueType const mIsovalue;
2639  BoolTreeType mTempMask;
2640  BoolTreeType * const mMask;
2641  Vec4i const * const mTileArray;
2642 }; // MaskTileBorders
2643 
2644 
2645 template<typename InputTreeType>
2646 void
2647 MaskTileBorders<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
2648 {
2649  tree::ValueAccessor<const InputTreeType> inputTreeAcc(*mInputTree);
2650 
2651  CoordBBox region, bbox;
2652  Coord ijk, nijk;
2653 
2654  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
2655 
2656  const Vec4i& tile = mTileArray[n];
2657 
2658  bbox.min()[0] = tile[0];
2659  bbox.min()[1] = tile[1];
2660  bbox.min()[2] = tile[2];
2661  bbox.max() = bbox.min();
2662  bbox.max().offset(tile[3]);
2663 
2664  InputValueType value = mInputTree->background();
2665 
2666  const bool isInside = isInsideValue(inputTreeAcc.getValue(bbox.min()), mIsovalue);
2667  const int valueDepth = inputTreeAcc.getValueDepth(bbox.min());
2668 
2669  // eval x-edges
2670 
2671  ijk = bbox.max();
2672  nijk = ijk;
2673  ++nijk[0];
2674 
2675  bool processRegion = true;
2676  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2677  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2678  }
2679 
2680  if (processRegion) {
2681  region = bbox;
2682  region.expand(1);
2683  region.min()[0] = region.max()[0] = ijk[0];
2684  mMask->fill(region, false);
2685  }
2686 
2687 
2688  ijk = bbox.min();
2689  --ijk[0];
2690 
2691  processRegion = true;
2692  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2693  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2694  && isInside != isInsideValue(value, mIsovalue));
2695  }
2696 
2697  if (processRegion) {
2698  region = bbox;
2699  region.expand(1);
2700  region.min()[0] = region.max()[0] = ijk[0];
2701  mMask->fill(region, false);
2702  }
2703 
2704 
2705  // eval y-edges
2706 
2707  ijk = bbox.max();
2708  nijk = ijk;
2709  ++nijk[1];
2710 
2711  processRegion = true;
2712  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2713  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2714  }
2715 
2716  if (processRegion) {
2717  region = bbox;
2718  region.expand(1);
2719  region.min()[1] = region.max()[1] = ijk[1];
2720  mMask->fill(region, false);
2721  }
2722 
2723 
2724  ijk = bbox.min();
2725  --ijk[1];
2726 
2727  processRegion = true;
2728  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2729  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2730  && isInside != isInsideValue(value, mIsovalue));
2731  }
2732 
2733  if (processRegion) {
2734  region = bbox;
2735  region.expand(1);
2736  region.min()[1] = region.max()[1] = ijk[1];
2737  mMask->fill(region, false);
2738  }
2739 
2740 
2741  // eval z-edges
2742 
2743  ijk = bbox.max();
2744  nijk = ijk;
2745  ++nijk[2];
2746 
2747  processRegion = true;
2748  if (valueDepth >= inputTreeAcc.getValueDepth(nijk)) {
2749  processRegion = isInside != isInsideValue(inputTreeAcc.getValue(nijk), mIsovalue);
2750  }
2751 
2752  if (processRegion) {
2753  region = bbox;
2754  region.expand(1);
2755  region.min()[2] = region.max()[2] = ijk[2];
2756  mMask->fill(region, false);
2757  }
2758 
2759  ijk = bbox.min();
2760  --ijk[2];
2761 
2762  processRegion = true;
2763  if (valueDepth >= inputTreeAcc.getValueDepth(ijk)) {
2764  processRegion = (!inputTreeAcc.probeValue(ijk, value)
2765  && isInside != isInsideValue(value, mIsovalue));
2766  }
2767 
2768  if (processRegion) {
2769  region = bbox;
2770  region.expand(1);
2771  region.min()[2] = region.max()[2] = ijk[2];
2772  mMask->fill(region, false);
2773  }
2774  }
2775 } // MaskTileBorders::operator()
2776 
2777 
2778 template<typename InputTreeType>
2779 inline void
2780 maskActiveTileBorders(const InputTreeType& inputTree,
2781  const typename InputTreeType::ValueType iso,
2782  typename InputTreeType::template ValueConverter<bool>::Type& mask)
2783 {
2784  typename InputTreeType::ValueOnCIter tileIter(inputTree);
2785  tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2786 
2787  size_t tileCount = 0;
2788  for ( ; tileIter; ++tileIter) {
2789  ++tileCount;
2790  }
2791 
2792  if (tileCount > 0) {
2793  std::unique_ptr<Vec4i[]> tiles(new Vec4i[tileCount]);
2794 
2795  CoordBBox bbox;
2796  size_t index = 0;
2797 
2798  tileIter = inputTree.cbeginValueOn();
2799  tileIter.setMaxDepth(InputTreeType::ValueOnCIter::LEAF_DEPTH - 1);
2800 
2801  for (; tileIter; ++tileIter) {
2802  Vec4i& tile = tiles[index++];
2803  tileIter.getBoundingBox(bbox);
2804  tile[0] = bbox.min()[0];
2805  tile[1] = bbox.min()[1];
2806  tile[2] = bbox.min()[2];
2807  tile[3] = bbox.max()[0] - bbox.min()[0];
2808  }
2809 
2810  MaskTileBorders<InputTreeType> op(inputTree, iso, mask, tiles.get());
2811  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, tileCount), op);
2812  }
2813 }
2814 
2815 
2816 ////////////////////////////////////////
2817 
2818 
2819 // Utility class for the volumeToMesh wrapper
2820 class PointListCopy
2821 {
2822 public:
2823  PointListCopy(const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2824  : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2825  {
2826  }
2827 
2828  void operator()(const tbb::blocked_range<size_t>& range) const
2829  {
2830  for (size_t n = range.begin(); n < range.end(); ++n) {
2831  mPointsOut[n] = mPointsIn[n];
2832  }
2833  }
2834 
2835 private:
2836  const PointList& mPointsIn;
2837  std::vector<Vec3s>& mPointsOut;
2838 };
2839 
2840 
2841 ////////////////////////////////////////////////////////////////////////////////
2842 ////////////////////////////////////////////////////////////////////////////////
2843 
2844 
2845 struct LeafNodeVoxelOffsets
2846 {
2847  using IndexVector = std::vector<Index>;
2848 
2849  template<typename LeafNodeType>
2850  void constructOffsetList();
2851 
2852  /// Return internal core voxel offsets.
2853  const IndexVector& core() const { return mCore; }
2854 
2855 
2856  /// Return front face voxel offsets.
2857  const IndexVector& minX() const { return mMinX; }
2858 
2859  /// Return back face voxel offsets.
2860  const IndexVector& maxX() const { return mMaxX; }
2861 
2862 
2863  /// Return bottom face voxel offsets.
2864  const IndexVector& minY() const { return mMinY; }
2865 
2866  /// Return top face voxel offsets.
2867  const IndexVector& maxY() const { return mMaxY; }
2868 
2869 
2870  /// Return left face voxel offsets.
2871  const IndexVector& minZ() const { return mMinZ; }
2872 
2873  /// Return right face voxel offsets.
2874  const IndexVector& maxZ() const { return mMaxZ; }
2875 
2876 
2877  /// Return voxel offsets with internal neighbours in x + 1.
2878  const IndexVector& internalNeighborsX() const { return mInternalNeighborsX; }
2879 
2880  /// Return voxel offsets with internal neighbours in y + 1.
2881  const IndexVector& internalNeighborsY() const { return mInternalNeighborsY; }
2882 
2883  /// Return voxel offsets with internal neighbours in z + 1.
2884  const IndexVector& internalNeighborsZ() const { return mInternalNeighborsZ; }
2885 
2886 
2887 private:
2888  IndexVector mCore, mMinX, mMaxX, mMinY, mMaxY, mMinZ, mMaxZ,
2889  mInternalNeighborsX, mInternalNeighborsY, mInternalNeighborsZ;
2890 }; // struct LeafNodeOffsets
2891 
2892 
2893 template<typename LeafNodeType>
2894 inline void
2895 LeafNodeVoxelOffsets::constructOffsetList()
2896 {
2897  // internal core voxels
2898  mCore.clear();
2899  mCore.reserve((LeafNodeType::DIM - 2) * (LeafNodeType::DIM - 2));
2900 
2901  for (Index x = 1; x < (LeafNodeType::DIM - 1); ++x) {
2902  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2903  for (Index y = 1; y < (LeafNodeType::DIM - 1); ++y) {
2904  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2905  for (Index z = 1; z < (LeafNodeType::DIM - 1); ++z) {
2906  mCore.push_back(offsetXY + z);
2907  }
2908  }
2909  }
2910 
2911  // internal neighbors in x + 1
2912  mInternalNeighborsX.clear();
2913  mInternalNeighborsX.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2914 
2915  for (Index x = 0; x < (LeafNodeType::DIM - 1); ++x) {
2916  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2917  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2918  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2919  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2920  mInternalNeighborsX.push_back(offsetXY + z);
2921  }
2922  }
2923  }
2924 
2925  // internal neighbors in y + 1
2926  mInternalNeighborsY.clear();
2927  mInternalNeighborsY.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2928 
2929  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2930  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2931  for (Index y = 0; y < (LeafNodeType::DIM - 1); ++y) {
2932  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2933  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2934  mInternalNeighborsY.push_back(offsetXY + z);
2935  }
2936  }
2937  }
2938 
2939  // internal neighbors in z + 1
2940  mInternalNeighborsZ.clear();
2941  mInternalNeighborsZ.reserve(LeafNodeType::SIZE - (LeafNodeType::DIM * LeafNodeType::DIM));
2942 
2943  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2944  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2945  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2946  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2947  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2948  mInternalNeighborsZ.push_back(offsetXY + z);
2949  }
2950  }
2951  }
2952 
2953  // min x
2954  mMinX.clear();
2955  mMinX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2956  {
2957  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2958  const Index offsetXY = (y << LeafNodeType::LOG2DIM);
2959  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2960  mMinX.push_back(offsetXY + z);
2961  }
2962  }
2963  }
2964 
2965  // max x
2966  mMaxX.clear();
2967  mMaxX.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2968  {
2969  const Index offsetX = (LeafNodeType::DIM - 1) << (2 * LeafNodeType::LOG2DIM);
2970  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
2971  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
2972  for (Index z = 0; z < LeafNodeType::DIM; ++z) {
2973  mMaxX.push_back(offsetXY + z);
2974  }
2975  }
2976  }
2977 
2978  // min y
2979  mMinY.clear();
2980  mMinY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2981  {
2982  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2983  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2984  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2985  mMinY.push_back(offsetX + z);
2986  }
2987  }
2988  }
2989 
2990  // max y
2991  mMaxY.clear();
2992  mMaxY.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
2993  {
2994  const Index offsetY = (LeafNodeType::DIM - 1) << LeafNodeType::LOG2DIM;
2995  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
2996  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
2997  for (Index z = 0; z < (LeafNodeType::DIM - 1); ++z) {
2998  mMaxY.push_back(offsetX + offsetY + z);
2999  }
3000  }
3001  }
3002 
3003  // min z
3004  mMinZ.clear();
3005  mMinZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3006  {
3007  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
3008  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3009  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
3010  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3011  mMinZ.push_back(offsetXY);
3012  }
3013  }
3014  }
3015 
3016  // max z
3017  mMaxZ.clear();
3018  mMaxZ.reserve(LeafNodeType::DIM * LeafNodeType::DIM);
3019  {
3020  for (Index x = 0; x < LeafNodeType::DIM; ++x) {
3021  const Index offsetX = x << (2 * LeafNodeType::LOG2DIM);
3022  for (Index y = 0; y < LeafNodeType::DIM; ++y) {
3023  const Index offsetXY = offsetX + (y << LeafNodeType::LOG2DIM);
3024  mMaxZ.push_back(offsetXY + (LeafNodeType::DIM - 1));
3025  }
3026  }
3027  }
3028 }
3029 
3030 
3031 ////////////////////////////////////////
3032 
3033 
3034 /// Utility method to marks all voxels that share an edge.
3035 template<typename AccessorT, int _AXIS>
3036 struct VoxelEdgeAccessor
3037 {
3038  enum { AXIS = _AXIS };
3039  AccessorT& acc;
3040 
3041  VoxelEdgeAccessor(AccessorT& _acc) : acc(_acc) {}
3042 
3043  void set(Coord ijk) {
3044  if (_AXIS == 0) { // x + 1 edge
3045  acc.setActiveState(ijk);
3046  --ijk[1]; // set i, j-1, k
3047  acc.setActiveState(ijk);
3048  --ijk[2]; // set i, j-1, k-1
3049  acc.setActiveState(ijk);
3050  ++ijk[1]; // set i, j, k-1
3051  acc.setActiveState(ijk);
3052  } else if (_AXIS == 1) { // y + 1 edge
3053  acc.setActiveState(ijk);
3054  --ijk[2]; // set i, j, k-1
3055  acc.setActiveState(ijk);
3056  --ijk[0]; // set i-1, j, k-1
3057  acc.setActiveState(ijk);
3058  ++ijk[2]; // set i-1, j, k
3059  acc.setActiveState(ijk);
3060  } else { // z + 1 edge
3061  acc.setActiveState(ijk);
3062  --ijk[1]; // set i, j-1, k
3063  acc.setActiveState(ijk);
3064  --ijk[0]; // set i-1, j-1, k
3065  acc.setActiveState(ijk);
3066  ++ijk[1]; // set i-1, j, k
3067  acc.setActiveState(ijk);
3068  }
3069  }
3070 };
3071 
3072 
3073 /// Utility method to check for sign changes along the x + 1, y + 1 or z + 1 directions.
3074 /// The direction is determined by the @a edgeAcc parameter. Only voxels that have internal
3075 /// neighbours are evaluated.
3076 template<typename VoxelEdgeAcc, typename LeafNodeT>
3077 void
3078 evalInternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3079  const LeafNodeT& leafnode,
3080  const LeafNodeVoxelOffsets& voxels,
3081  const typename LeafNodeT::ValueType iso)
3082 {
3083  Index nvo = 1; // neighbour voxel offset, z + 1 direction assumed initially.
3084  const std::vector<Index>* offsets = &voxels.internalNeighborsZ();
3085 
3086  if (VoxelEdgeAcc::AXIS == 0) { // x + 1 direction
3087  nvo = LeafNodeT::DIM * LeafNodeT::DIM;
3088  offsets = &voxels.internalNeighborsX();
3089  } else if (VoxelEdgeAcc::AXIS == 1) { // y + 1 direction
3090  nvo = LeafNodeT::DIM;
3091  offsets = &voxels.internalNeighborsY();
3092  }
3093 
3094  const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3095 
3096  for (size_t n = 0, N = offsets->size(); n < N; ++n) {
3097  const Index& pos = (*offsets)[n];
3098  const bool isActive = leafnode.isValueOn(pos) || leafnode.isValueOn(pos + nvo);
3099  if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3100  isInsideValue(lhsAcc.get((pos + nvo)), iso))) {
3101  edgeAcc.set(leafnode.offsetToGlobalCoord(pos));
3102  }
3103  }
3104 }
3105 
3106 
3107 /// Utility method to check for sign changes along the x + 1, y + 1 or z + 1 directions.
3108 /// The direction is determined by the @a edgeAcc parameter. All voxels that reside in the
3109 /// specified leafnode face: back, top or right are evaluated.
3110 template<typename LeafNodeT, typename TreeAcc, typename VoxelEdgeAcc>
3111 void
3112 evalExternalVoxelEdges(VoxelEdgeAcc& edgeAcc,
3113  const TreeAcc& acc,
3114  const LeafNodeT& lhsNode,
3115  const LeafNodeVoxelOffsets& voxels,
3116  const typename LeafNodeT::ValueType iso)
3117 {
3118  const std::vector<Index>* lhsOffsets = &voxels.maxX();
3119  const std::vector<Index>* rhsOffsets = &voxels.minX();
3120  Coord ijk = lhsNode.origin();
3121 
3122  if (VoxelEdgeAcc::AXIS == 0) { // back leafnode face
3123  ijk[0] += LeafNodeT::DIM;
3124  } else if (VoxelEdgeAcc::AXIS == 1) { // top leafnode face
3125  ijk[1] += LeafNodeT::DIM;
3126  lhsOffsets = &voxels.maxY();
3127  rhsOffsets = &voxels.minY();
3128  } else if (VoxelEdgeAcc::AXIS == 2) { // right leafnode face
3129  ijk[2] += LeafNodeT::DIM;
3130  lhsOffsets = &voxels.maxZ();
3131  rhsOffsets = &voxels.minZ();
3132  }
3133 
3134  typename LeafNodeT::ValueType value;
3135  const LeafNodeT* rhsNodePt = acc.probeConstLeaf(ijk);
3136 
3137  const LeafBufferAccessor<LeafNodeT> lhsAcc(lhsNode);
3138 
3139  if (rhsNodePt) {
3140  const LeafBufferAccessor<LeafNodeT> rhsAcc(*rhsNodePt);
3141 
3142  for (size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3143  const Index& pos = (*lhsOffsets)[n];
3144  bool isActive = lhsNode.isValueOn(pos) || rhsNodePt->isValueOn((*rhsOffsets)[n]);
3145  if (isActive && (isInsideValue(lhsAcc.get(pos), iso) !=
3146  isInsideValue(rhsAcc.get((*rhsOffsets)[n]), iso))) {
3147  edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3148  }
3149  }
3150  } else if (!acc.probeValue(ijk, value)) {
3151  const bool inside = isInsideValue(value, iso);
3152  for (size_t n = 0, N = lhsOffsets->size(); n < N; ++n) {
3153  const Index& pos = (*lhsOffsets)[n];
3154  if (lhsNode.isValueOn(pos) && (inside != isInsideValue(lhsAcc.get(pos), iso))) {
3155  edgeAcc.set(lhsNode.offsetToGlobalCoord(pos));
3156  }
3157  }
3158  }
3159 }
3160 
3161 
3162 /// Utility method to check for sign changes along the x - 1, y - 1 or z - 1 directions.
3163 /// The direction is determined by the @a edgeAcc parameter. All voxels that reside in the
3164 /// specified leafnode face: front, bottom or left are evaluated.
3165 template<typename LeafNodeT, typename TreeAcc, typename VoxelEdgeAcc>
3166 void
3167 evalExternalVoxelEdgesInv(VoxelEdgeAcc& edgeAcc,
3168  const TreeAcc& acc,
3169  const LeafNodeT& leafnode,
3170  const LeafNodeVoxelOffsets& voxels,
3171  const typename LeafNodeT::ValueType iso)
3172 {
3173  Coord ijk = leafnode.origin();
3174  if (VoxelEdgeAcc::AXIS == 0) --ijk[0]; // front leafnode face
3175  else if (VoxelEdgeAcc::AXIS == 1) --ijk[1]; // bottom leafnode face
3176  else if (VoxelEdgeAcc::AXIS == 2) --ijk[2]; // left leafnode face
3177 
3178  typename LeafNodeT::ValueType value;
3179  if (!acc.probeConstLeaf(ijk) && !acc.probeValue(ijk, value)) {
3180 
3181  const std::vector<Index>* offsets = &voxels.internalNeighborsX();
3182  if (VoxelEdgeAcc::AXIS == 1) offsets = &voxels.internalNeighborsY();
3183  else if (VoxelEdgeAcc::AXIS == 2) offsets = &voxels.internalNeighborsZ();
3184 
3185  const LeafBufferAccessor<LeafNodeT> lhsAcc(leafnode);
3186 
3187  const bool inside = isInsideValue(value, iso);
3188  for (size_t n = 0, N = offsets->size(); n < N; ++n) {
3189 
3190  const Index& pos = (*offsets)[n];
3191  if (leafnode.isValueOn(pos)
3192  && (inside != isInsideValue(lhsAcc.get(pos), iso)))
3193  {
3194  ijk = leafnode.offsetToGlobalCoord(pos);
3195  if (VoxelEdgeAcc::AXIS == 0) --ijk[0];
3196  else if (VoxelEdgeAcc::AXIS == 1) --ijk[1];
3197  else if (VoxelEdgeAcc::AXIS == 2) --ijk[2];
3198 
3199  edgeAcc.set(ijk);
3200  }
3201  }
3202  }
3203 }
3204 
3205 
3206 template<typename InputTreeType>
3207 struct IdentifyIntersectingVoxels
3208 {
3209  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3210  using InputValueType = typename InputLeafNodeType::ValueType;
3211 
3212  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3213 
3214  IdentifyIntersectingVoxels(
3215  const InputTreeType& inputTree,
3216  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3217  const LeafNodeVoxelOffsets& offsets,
3218  BoolTreeType& intersectionTree,
3219  InputValueType iso);
3220 
3221  IdentifyIntersectingVoxels(IdentifyIntersectingVoxels&, tbb::split);
3222  void operator()(const tbb::blocked_range<size_t>&);
3223  void join(const IdentifyIntersectingVoxels& rhs) {
3224  mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3225  }
3226 
3227 private:
3228  tree::ValueAccessor<const InputTreeType> mInputAccessor;
3229  InputLeafNodeType const * const * const mInputNodes;
3230 
3231  BoolTreeType mIntersectionTree;
3232  tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3233 
3234  const LeafNodeVoxelOffsets& mOffsets;
3235  const InputValueType mIsovalue;
3236 }; // struct IdentifyIntersectingVoxels
3237 
3238 
3239 template<typename InputTreeType>
3240 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3241  const InputTreeType& inputTree,
3242  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
3243  const LeafNodeVoxelOffsets& offsets,
3244  BoolTreeType& intersectionTree,
3245  InputValueType iso)
3246  : mInputAccessor(inputTree)
3247  , mInputNodes(inputLeafNodes.data())
3248  , mIntersectionTree(false)
3249  , mIntersectionAccessor(intersectionTree)
3250  , mOffsets(offsets)
3251  , mIsovalue(iso)
3252 {
3253 }
3254 
3255 
3256 template<typename InputTreeType>
3257 IdentifyIntersectingVoxels<InputTreeType>::IdentifyIntersectingVoxels(
3258  IdentifyIntersectingVoxels& rhs, tbb::split)
3259  : mInputAccessor(rhs.mInputAccessor.tree())
3260  , mInputNodes(rhs.mInputNodes)
3261  , mIntersectionTree(false)
3262  , mIntersectionAccessor(mIntersectionTree) // use local tree.
3263  , mOffsets(rhs.mOffsets)
3264  , mIsovalue(rhs.mIsovalue)
3265 {
3266 }
3267 
3268 
3269 template<typename InputTreeType>
3270 void
3271 IdentifyIntersectingVoxels<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3272 {
3273  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3274  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3275  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3276 
3277  for (size_t n = range.begin(); n != range.end(); ++n) {
3278 
3279  const InputLeafNodeType& node = *mInputNodes[n];
3280 
3281  // internal x + 1 voxel edges
3282  evalInternalVoxelEdges(xEdgeAcc, node, mOffsets, mIsovalue);
3283  // internal y + 1 voxel edges
3284  evalInternalVoxelEdges(yEdgeAcc, node, mOffsets, mIsovalue);
3285  // internal z + 1 voxel edges
3286  evalInternalVoxelEdges(zEdgeAcc, node, mOffsets, mIsovalue);
3287 
3288  // external x + 1 voxels edges (back face)
3289  evalExternalVoxelEdges(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3290  // external y + 1 voxels edges (top face)
3291  evalExternalVoxelEdges(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3292  // external z + 1 voxels edges (right face)
3293  evalExternalVoxelEdges(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3294 
3295  // The remaining edges are only checked if the leafnode neighbour, in the
3296  // corresponding direction, is an inactive tile.
3297 
3298  // external x - 1 voxels edges (front face)
3299  evalExternalVoxelEdgesInv(xEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3300  // external y - 1 voxels edges (bottom face)
3301  evalExternalVoxelEdgesInv(yEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3302  // external z - 1 voxels edges (left face)
3303  evalExternalVoxelEdgesInv(zEdgeAcc, mInputAccessor, node, mOffsets, mIsovalue);
3304  }
3305 } // IdentifyIntersectingVoxels::operator()
3306 
3307 
3308 template<typename InputTreeType>
3309 inline void
3310 identifySurfaceIntersectingVoxels(
3311  typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3312  const InputTreeType& inputTree,
3313  typename InputTreeType::ValueType isovalue)
3314 {
3315  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3316 
3317  std::vector<const InputLeafNodeType*> inputLeafNodes;
3318  inputTree.getNodes(inputLeafNodes);
3319 
3320  LeafNodeVoxelOffsets offsets;
3321  offsets.constructOffsetList<InputLeafNodeType>();
3322 
3323  IdentifyIntersectingVoxels<InputTreeType> op(
3324  inputTree, inputLeafNodes, offsets, intersectionTree, isovalue);
3325 
3326  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, inputLeafNodes.size()), op);
3327 
3328  maskActiveTileBorders(inputTree, isovalue, intersectionTree);
3329 }
3330 
3331 
3332 ////////////////////////////////////////
3333 
3334 
3335 template<typename InputTreeType>
3336 struct MaskIntersectingVoxels
3337 {
3338  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3339  using InputValueType = typename InputLeafNodeType::ValueType;
3340 
3341  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3342  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3343 
3344  MaskIntersectingVoxels(
3345  const InputTreeType& inputTree,
3346  const std::vector<BoolLeafNodeType*>& nodes,
3347  BoolTreeType& intersectionTree,
3348  InputValueType iso);
3349 
3350  MaskIntersectingVoxels(MaskIntersectingVoxels&, tbb::split);
3351  void operator()(const tbb::blocked_range<size_t>&);
3352  void join(const MaskIntersectingVoxels& rhs) {
3353  mIntersectionAccessor.tree().merge(rhs.mIntersectionAccessor.tree());
3354  }
3355 
3356 private:
3357  tree::ValueAccessor<const InputTreeType> mInputAccessor;
3358  BoolLeafNodeType const * const * const mNodes;
3359 
3360  BoolTreeType mIntersectionTree;
3361  tree::ValueAccessor<BoolTreeType> mIntersectionAccessor;
3362 
3363  const InputValueType mIsovalue;
3364 }; // struct MaskIntersectingVoxels
3365 
3366 
3367 template<typename InputTreeType>
3368 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3369  const InputTreeType& inputTree,
3370  const std::vector<BoolLeafNodeType*>& nodes,
3371  BoolTreeType& intersectionTree,
3372  InputValueType iso)
3373  : mInputAccessor(inputTree)
3374  , mNodes(nodes.data())
3375  , mIntersectionTree(false)
3376  , mIntersectionAccessor(intersectionTree)
3377  , mIsovalue(iso)
3378 {
3379 }
3380 
3381 
3382 template<typename InputTreeType>
3383 MaskIntersectingVoxels<InputTreeType>::MaskIntersectingVoxels(
3384  MaskIntersectingVoxels& rhs, tbb::split)
3385  : mInputAccessor(rhs.mInputAccessor.tree())
3386  , mNodes(rhs.mNodes)
3387  , mIntersectionTree(false)
3388  , mIntersectionAccessor(mIntersectionTree) // use local tree.
3389  , mIsovalue(rhs.mIsovalue)
3390 {
3391 }
3392 
3393 
3394 template<typename InputTreeType>
3395 void
3396 MaskIntersectingVoxels<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3397 {
3398  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 0> xEdgeAcc(mIntersectionAccessor);
3399  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 1> yEdgeAcc(mIntersectionAccessor);
3400  VoxelEdgeAccessor<tree::ValueAccessor<BoolTreeType>, 2> zEdgeAcc(mIntersectionAccessor);
3401 
3402  Coord ijk;
3403 
3404  for (size_t n = range.begin(); n != range.end(); ++n) {
3405 
3406  const BoolLeafNodeType& node = *mNodes[n];
3407 
3408  for (auto it = node.cbeginValueOn(); it; ++it) {
3409 
3410  if (!it.getValue()) {
3411 
3412  ijk = it.getCoord();
3413 
3414  const bool inside = isInsideValue(mInputAccessor.getValue(ijk), mIsovalue);
3415 
3416  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(1, 0, 0)), mIsovalue)) {
3417  xEdgeAcc.set(ijk);
3418  }
3419 
3420  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 1, 0)), mIsovalue)) {
3421  yEdgeAcc.set(ijk);
3422  }
3423 
3424  if (inside != isInsideValue(mInputAccessor.getValue(ijk.offsetBy(0, 0, 1)), mIsovalue)) {
3425  zEdgeAcc.set(ijk);
3426  }
3427  }
3428  }
3429  }
3430 } // MaskIntersectingVoxels::operator()
3431 
3432 
3433 template<typename BoolTreeType>
3434 struct MaskBorderVoxels
3435 {
3436  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3437 
3438  MaskBorderVoxels(const BoolTreeType& maskTree,
3439  const std::vector<BoolLeafNodeType*>& maskNodes,
3440  BoolTreeType& borderTree)
3441  : mMaskTree(&maskTree)
3442  , mMaskNodes(maskNodes.data())
3443  , mTmpBorderTree(false)
3444  , mBorderTree(&borderTree) {}
3445 
3446  MaskBorderVoxels(MaskBorderVoxels& rhs, tbb::split)
3447  : mMaskTree(rhs.mMaskTree)
3448  , mMaskNodes(rhs.mMaskNodes)
3449  , mTmpBorderTree(false)
3450  , mBorderTree(&mTmpBorderTree) {}
3451 
3452  void join(MaskBorderVoxels& rhs) { mBorderTree->merge(*rhs.mBorderTree); }
3453 
3454  void operator()(const tbb::blocked_range<size_t>& range)
3455  {
3456  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
3457  tree::ValueAccessor<BoolTreeType> borderAcc(*mBorderTree);
3458  Coord ijk;
3459 
3460  for (size_t n = range.begin(); n != range.end(); ++n) {
3461 
3462  const BoolLeafNodeType& node = *mMaskNodes[n];
3463 
3464  for (auto it = node.cbeginValueOn(); it; ++it) {
3465 
3466  ijk = it.getCoord();
3467 
3468  const bool lhs = it.getValue();
3469  bool rhs = lhs;
3470 
3471  bool isEdgeVoxel = false;
3472 
3473  ijk[2] += 1; // i, j, k+1
3474  isEdgeVoxel = (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3475 
3476  ijk[1] += 1; // i, j+1, k+1
3477  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3478 
3479  ijk[0] += 1; // i+1, j+1, k+1
3480  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3481 
3482  ijk[1] -= 1; // i+1, j, k+1
3483  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3484 
3485 
3486  ijk[2] -= 1; // i+1, j, k
3487  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3488 
3489  ijk[1] += 1; // i+1, j+1, k
3490  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3491 
3492  ijk[0] -= 1; // i, j+1, k
3493  isEdgeVoxel = isEdgeVoxel || (maskAcc.probeValue(ijk, rhs) && lhs != rhs);
3494 
3495  if (isEdgeVoxel) {
3496  ijk[1] -= 1; // i, j, k
3497  borderAcc.setValue(ijk, true);
3498  }
3499  }
3500  }
3501  }
3502 
3503 private:
3504  BoolTreeType const * const mMaskTree;
3505  BoolLeafNodeType const * const * const mMaskNodes;
3506 
3507  BoolTreeType mTmpBorderTree;
3508  BoolTreeType * const mBorderTree;
3509 }; // struct MaskBorderVoxels
3510 
3511 
3512 template<typename BoolTreeType>
3513 struct SyncMaskValues
3514 {
3515  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3516 
3517  SyncMaskValues(const std::vector<BoolLeafNodeType*>& nodes, const BoolTreeType& mask)
3518  : mNodes(nodes.data())
3519  , mMaskTree(&mask) {}
3520 
3521  void operator()(const tbb::blocked_range<size_t>& range) const
3522  {
3523  using ValueOnIter = typename BoolLeafNodeType::ValueOnIter;
3524 
3525  tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3526 
3527  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3528 
3529  BoolLeafNodeType& node = *mNodes[n];
3530 
3531  const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3532 
3533  if (maskNode) {
3534  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3535  const Index pos = it.pos();
3536  if (maskNode->getValue(pos)) {
3537  node.setValueOnly(pos, true);
3538  }
3539  }
3540  }
3541  }
3542  }
3543 
3544 private:
3545  BoolLeafNodeType * const * const mNodes;
3546  BoolTreeType const * const mMaskTree;
3547 }; // struct SyncMaskValues
3548 
3549 
3550 ////////////////////////////////////////
3551 
3552 
3553 template<typename BoolTreeType>
3554 struct MaskSurface
3555 {
3556  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3557 
3558  MaskSurface(const std::vector<BoolLeafNodeType*>& nodes,
3559  const BoolTreeType& mask,
3560  const math::Transform& inputTransform,
3561  const math::Transform& maskTransform,
3562  const bool invert)
3563  : mNodes(nodes.data())
3564  , mMaskTree(&mask)
3565  , mInputTransform(inputTransform)
3566  , mMaskTransform(maskTransform)
3567  , mMatchingTransforms(mInputTransform == mMaskTransform)
3568  , mInvertMask(invert) {}
3569 
3570  void operator()(const tbb::blocked_range<size_t>& range) const
3571  {
3572  using ValueOnIter = typename BoolLeafNodeType::ValueOnIter;
3573 
3574  tree::ValueAccessor<const BoolTreeType> maskTreeAcc(*mMaskTree);
3575 
3576  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3577 
3578  BoolLeafNodeType& node = *mNodes[n];
3579 
3580  if (mMatchingTransforms) {
3581 
3582  const BoolLeafNodeType * maskNode = maskTreeAcc.probeConstLeaf(node.origin());
3583 
3584  if (maskNode) {
3585 
3586  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3587  const Index pos = it.pos();
3588  if (maskNode->isValueOn(pos) == mInvertMask) {
3589  node.setValueOnly(pos, true);
3590  }
3591  }
3592 
3593  } else {
3594 
3595  if (maskTreeAcc.isValueOn(node.origin()) == mInvertMask) {
3596  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3597  node.setValueOnly(it.pos(), true);
3598  }
3599  }
3600 
3601  }
3602 
3603  } else {
3604 
3605  Coord ijk;
3606 
3607  for (ValueOnIter it = node.beginValueOn(); it; ++it) {
3608 
3609  ijk = mMaskTransform.worldToIndexCellCentered(
3610  mInputTransform.indexToWorld(it.getCoord()));
3611 
3612  if (maskTreeAcc.isValueOn(ijk) == mInvertMask) {
3613  node.setValueOnly(it.pos(), true);
3614  }
3615  }
3616 
3617  }
3618  }
3619  }
3620 
3621 private:
3622  BoolLeafNodeType * const * const mNodes;
3623  BoolTreeType const * const mMaskTree;
3624  const math::Transform& mInputTransform;
3625  const math::Transform& mMaskTransform;
3626  const bool mMatchingTransforms;
3627  const bool mInvertMask;
3628 }; // struct MaskSurface
3629 
3630 
3631 template<typename InputGridType>
3632 inline void
3633 applySurfaceMask(
3634  typename InputGridType::TreeType::template ValueConverter<bool>::Type& intersectionTree,
3635  typename InputGridType::TreeType::template ValueConverter<bool>::Type& borderTree,
3636  const InputGridType& inputGrid,
3637  const GridBase::ConstPtr& maskGrid,
3638  const bool invertMask,
3639  const typename InputGridType::ValueType isovalue)
3640 {
3641  using InputTreeType = typename InputGridType::TreeType;
3642  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3643  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3644  using BoolGridType = Grid<BoolTreeType>;
3645 
3646  if (!maskGrid) return;
3647  if (maskGrid->type() != BoolGridType::gridType()) return;
3648 
3649  const math::Transform& transform = inputGrid.transform();
3650  const InputTreeType& inputTree = inputGrid.tree();
3651 
3652  const BoolGridType * surfaceMask = static_cast<const BoolGridType*>(maskGrid.get());
3653 
3654  const BoolTreeType& maskTree = surfaceMask->tree();
3655  const math::Transform& maskTransform = surfaceMask->transform();
3656 
3657  // mark masked voxels
3658 
3659  std::vector<BoolLeafNodeType*> intersectionLeafNodes;
3660  intersectionTree.getNodes(intersectionLeafNodes);
3661 
3662  const tbb::blocked_range<size_t> intersectionRange(0, intersectionLeafNodes.size());
3663 
3664  tbb::parallel_for(intersectionRange,
3665  MaskSurface<BoolTreeType>(
3666  intersectionLeafNodes, maskTree, transform, maskTransform, invertMask));
3667 
3668 
3669  // mask surface-mask border
3670 
3671  MaskBorderVoxels<BoolTreeType> borderOp(
3672  intersectionTree, intersectionLeafNodes, borderTree);
3673  tbb::parallel_reduce(intersectionRange, borderOp);
3674 
3675 
3676  // recompute isosurface intersection mask
3677 
3678  BoolTreeType tmpIntersectionTree(false);
3679 
3680  MaskIntersectingVoxels<InputTreeType> op(
3681  inputTree, intersectionLeafNodes, tmpIntersectionTree, isovalue);
3682 
3683  tbb::parallel_reduce(intersectionRange, op);
3684 
3685  std::vector<BoolLeafNodeType*> tmpIntersectionLeafNodes;
3686  tmpIntersectionTree.getNodes(tmpIntersectionLeafNodes);
3687 
3688  tbb::parallel_for(tbb::blocked_range<size_t>(0, tmpIntersectionLeafNodes.size()),
3689  SyncMaskValues<BoolTreeType>(tmpIntersectionLeafNodes, intersectionTree));
3690 
3691  intersectionTree.clear();
3692  intersectionTree.merge(tmpIntersectionTree);
3693 }
3694 
3695 
3696 ////////////////////////////////////////
3697 
3698 
3699 template<typename InputTreeType>
3700 struct ComputeAuxiliaryData
3701 {
3702  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
3703  using InputValueType = typename InputLeafNodeType::ValueType;
3704 
3705  using BoolLeafNodeType = tree::LeafNode<bool, InputLeafNodeType::LOG2DIM>;
3706 
3707  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
3708  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
3709 
3710 
3711  ComputeAuxiliaryData(const InputTreeType& inputTree,
3712  const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3713  Int16TreeType& signFlagsTree,
3714  Index32TreeType& pointIndexTree,
3715  InputValueType iso);
3716 
3717  ComputeAuxiliaryData(ComputeAuxiliaryData&, tbb::split);
3718  void operator()(const tbb::blocked_range<size_t>&);
3719  void join(const ComputeAuxiliaryData& rhs) {
3720  mSignFlagsAccessor.tree().merge(rhs.mSignFlagsAccessor.tree());
3721  mPointIndexAccessor.tree().merge(rhs.mPointIndexAccessor.tree());
3722  }
3723 
3724 private:
3725  tree::ValueAccessor<const InputTreeType> mInputAccessor;
3726  BoolLeafNodeType const * const * const mIntersectionNodes;
3727 
3728  Int16TreeType mSignFlagsTree;
3729  tree::ValueAccessor<Int16TreeType> mSignFlagsAccessor;
3730  Index32TreeType mPointIndexTree;
3731  tree::ValueAccessor<Index32TreeType> mPointIndexAccessor;
3732 
3733  const InputValueType mIsovalue;
3734 };
3735 
3736 
3737 template<typename InputTreeType>
3738 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(
3739  const InputTreeType& inputTree,
3740  const std::vector<const BoolLeafNodeType*>& intersectionLeafNodes,
3741  Int16TreeType& signFlagsTree,
3742  Index32TreeType& pointIndexTree,
3743  InputValueType iso)
3744  : mInputAccessor(inputTree)
3745  , mIntersectionNodes(intersectionLeafNodes.data())
3746  , mSignFlagsTree(0)
3747  , mSignFlagsAccessor(signFlagsTree)
3748  , mPointIndexTree(std::numeric_limits<Index32>::max())
3749  , mPointIndexAccessor(pointIndexTree)
3750  , mIsovalue(iso)
3751 {
3752  pointIndexTree.root().setBackground(std::numeric_limits<Index32>::max(), false);
3753 }
3754 
3755 
3756 template<typename InputTreeType>
3757 ComputeAuxiliaryData<InputTreeType>::ComputeAuxiliaryData(ComputeAuxiliaryData& rhs, tbb::split)
3758  : mInputAccessor(rhs.mInputAccessor.tree())
3759  , mIntersectionNodes(rhs.mIntersectionNodes)
3760  , mSignFlagsTree(0)
3761  , mSignFlagsAccessor(mSignFlagsTree)
3762  , mPointIndexTree(std::numeric_limits<Index32>::max())
3763  , mPointIndexAccessor(mPointIndexTree)
3764  , mIsovalue(rhs.mIsovalue)
3765 {
3766 }
3767 
3768 
3769 template<typename InputTreeType>
3770 void
3771 ComputeAuxiliaryData<InputTreeType>::operator()(const tbb::blocked_range<size_t>& range)
3772 {
3773  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
3774 
3775  Coord ijk;
3776  std::array<InputValueType, 8> cellVertexValues;
3777  std::unique_ptr<Int16LeafNodeType> signsNodePt(nullptr);
3778 
3779  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3780 
3781  const BoolLeafNodeType& maskNode = *mIntersectionNodes[n];
3782  const Coord& origin = maskNode.origin();
3783 
3784  const InputLeafNodeType* leafPt = mInputAccessor.probeConstLeaf(origin);
3785 
3786  if (!signsNodePt.get()) signsNodePt.reset(new Int16LeafNodeType(origin, 0));
3787  else signsNodePt->setOrigin(origin);
3788 
3789  bool updatedNode = false;
3790 
3791  for (auto it = maskNode.cbeginValueOn(); it; ++it) {
3792 
3793  const Index pos = it.pos();
3794  ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
3795 
3796  const bool inclusiveCell = leafPt && isInternalLeafCoord<InputLeafNodeType>(ijk);
3797 
3798  if (inclusiveCell) getCellVertexValues(*leafPt, pos, cellVertexValues);
3799  else getCellVertexValues(mInputAccessor, origin + ijk, cellVertexValues);
3800 
3801  uint8_t signFlags = computeSignFlags(cellVertexValues, mIsovalue);
3802 
3803  if (signFlags != 0 && signFlags != 0xFF) {
3804 
3805  const bool inside = signFlags & 0x1;
3806 
3807  int edgeFlags = inside ? INSIDE : 0;
3808 
3809  if (!it.getValue()) {
3810  edgeFlags |= inside != ((signFlags & 0x02) != 0) ? XEDGE : 0;
3811  edgeFlags |= inside != ((signFlags & 0x10) != 0) ? YEDGE : 0;
3812  edgeFlags |= inside != ((signFlags & 0x08) != 0) ? ZEDGE : 0;
3813  }
3814 
3815  const uint8_t ambiguousCellFlags = sAmbiguousFace[signFlags];
3816  if (ambiguousCellFlags != 0) {
3817  correctCellSigns(signFlags, ambiguousCellFlags, mInputAccessor,
3818  origin + ijk, mIsovalue);
3819  }
3820 
3821  edgeFlags |= int(signFlags);
3822 
3823  signsNodePt->setValueOn(pos, Int16(edgeFlags));
3824  updatedNode = true;
3825  }
3826  }
3827 
3828  if (updatedNode) {
3829  typename Index32TreeType::LeafNodeType* idxNode =
3830  mPointIndexAccessor.touchLeaf(origin);
3831  idxNode->topologyUnion(*signsNodePt);
3832 
3833  // zero fill
3834  auto* const idxData = idxNode->buffer().data();
3835  for (auto it = idxNode->beginValueOn(); it; ++it) {
3836  idxData[it.pos()] = 0;
3837  }
3838 
3839  mSignFlagsAccessor.addLeaf(signsNodePt.release());
3840  }
3841  }
3842 } // ComputeAuxiliaryData::operator()
3843 
3844 
3845 template<typename InputTreeType>
3846 inline void
3847 computeAuxiliaryData(
3848  typename InputTreeType::template ValueConverter<Int16>::Type& signFlagsTree,
3849  typename InputTreeType::template ValueConverter<Index32>::Type& pointIndexTree,
3850  const typename InputTreeType::template ValueConverter<bool>::Type& intersectionTree,
3851  const InputTreeType& inputTree,
3852  typename InputTreeType::ValueType isovalue)
3853 {
3854  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
3855  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
3856 
3857  std::vector<const BoolLeafNodeType*> intersectionLeafNodes;
3858  intersectionTree.getNodes(intersectionLeafNodes);
3859 
3860  ComputeAuxiliaryData<InputTreeType> op(
3861  inputTree, intersectionLeafNodes, signFlagsTree, pointIndexTree, isovalue);
3862 
3863  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, intersectionLeafNodes.size()), op);
3864 }
3865 
3866 
3867 ////////////////////////////////////////
3868 
3869 
3870 template<Index32 LeafNodeLog2Dim>
3871 struct LeafNodePointCount
3872 {
3873  using Int16LeafNodeType = tree::LeafNode<Int16, LeafNodeLog2Dim>;
3874 
3875  LeafNodePointCount(const std::vector<Int16LeafNodeType*>& leafNodes,
3876  std::unique_ptr<Index32[]>& leafNodeCount)
3877  : mLeafNodes(leafNodes.data())
3878  , mData(leafNodeCount.get()) {}
3879 
3880  void operator()(const tbb::blocked_range<size_t>& range) const {
3881 
3882  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3883 
3884  Index32 count = 0;
3885 
3886  Int16 const * p = mLeafNodes[n]->buffer().data();
3887  Int16 const * const endP = p + Int16LeafNodeType::SIZE;
3888 
3889  while (p < endP) {
3890  count += Index32(sEdgeGroupTable[(SIGNS & *p)][0]);
3891  ++p;
3892  }
3893 
3894  mData[n] = count;
3895  }
3896  }
3897 
3898 private:
3899  Int16LeafNodeType * const * const mLeafNodes;
3900  Index32 *mData;
3901 }; // struct LeafNodePointCount
3902 
3903 
3904 template<typename PointIndexLeafNode>
3905 struct AdaptiveLeafNodePointCount
3906 {
3907  using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3908 
3909  AdaptiveLeafNodePointCount(const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3910  const std::vector<Int16LeafNodeType*>& signDataNodes,
3911  std::unique_ptr<Index32[]>& leafNodeCount)
3912  : mPointIndexNodes(pointIndexNodes.data())
3913  , mSignDataNodes(signDataNodes.data())
3914  , mData(leafNodeCount.get()) {}
3915 
3916  void operator()(const tbb::blocked_range<size_t>& range) const
3917  {
3918  using IndexType = typename PointIndexLeafNode::ValueType;
3919 
3920  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3921 
3922  const PointIndexLeafNode& node = *mPointIndexNodes[n];
3923  const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3924 
3925  size_t count = 0;
3926 
3927  std::set<IndexType> uniqueRegions;
3928 
3929  for (typename PointIndexLeafNode::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
3930 
3931  IndexType id = it.getValue();
3932 
3933  if (id == 0) {
3934  count += size_t(sEdgeGroupTable[(SIGNS & signNode.getValue(it.pos()))][0]);
3935  } else if (id != IndexType(util::INVALID_IDX)) {
3936  uniqueRegions.insert(id);
3937  }
3938  }
3939 
3940  mData[n] = Index32(count + uniqueRegions.size());
3941  }
3942  }
3943 
3944 private:
3945  PointIndexLeafNode const * const * const mPointIndexNodes;
3946  Int16LeafNodeType const * const * const mSignDataNodes;
3947  Index32 *mData;
3948 }; // struct AdaptiveLeafNodePointCount
3949 
3950 
3951 template<typename PointIndexLeafNode>
3952 struct MapPoints
3953 {
3954  using Int16LeafNodeType = tree::LeafNode<Int16, PointIndexLeafNode::LOG2DIM>;
3955 
3956  MapPoints(const std::vector<PointIndexLeafNode*>& pointIndexNodes,
3957  const std::vector<Int16LeafNodeType*>& signDataNodes,
3958  std::unique_ptr<Index32[]>& leafNodeCount)
3959  : mPointIndexNodes(pointIndexNodes.data())
3960  , mSignDataNodes(signDataNodes.data())
3961  , mData(leafNodeCount.get()) {}
3962 
3963  void operator()(const tbb::blocked_range<size_t>& range) const {
3964 
3965  for (size_t n = range.begin(), N = range.end(); n != N; ++n) {
3966 
3967  const Int16LeafNodeType& signNode = *mSignDataNodes[n];
3968  PointIndexLeafNode& indexNode = *mPointIndexNodes[n];
3969 
3970  Index32 pointOffset = mData[n];
3971 
3972  for (auto it = indexNode.beginValueOn(); it; ++it) {
3973  const Index pos = it.pos();
3974  indexNode.setValueOnly(pos, pointOffset);
3975  const int signs = SIGNS & int(signNode.getValue(pos));
3976  pointOffset += Index32(sEdgeGroupTable[signs][0]);
3977  }
3978  }
3979  }
3980 
3981 private:
3982  PointIndexLeafNode * const * const mPointIndexNodes;
3983  Int16LeafNodeType const * const * const mSignDataNodes;
3984  Index32 * const mData;
3985 }; // struct MapPoints
3986 
3987 
3988 template<typename TreeType, typename PrimBuilder>
3989 struct ComputePolygons
3990 {
3991  using Int16TreeType = typename TreeType::template ValueConverter<Int16>::Type;
3992  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
3993 
3994  using Index32TreeType = typename TreeType::template ValueConverter<Index32>::Type;
3995  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
3996 
3997  ComputePolygons(
3998  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
3999  const Int16TreeType& signFlagsTree,
4000  const Index32TreeType& idxTree,
4001  PolygonPoolList& polygons,
4002  bool invertSurfaceOrientation);
4003 
4004  void setRefSignTree(const Int16TreeType * r) { mRefSignFlagsTree = r; }
4005 
4006  void operator()(const tbb::blocked_range<size_t>&) const;
4007 
4008 private:
4009  Int16LeafNodeType * const * const mSignFlagsLeafNodes;
4010  Int16TreeType const * const mSignFlagsTree;
4011  Int16TreeType const * mRefSignFlagsTree;
4012  Index32TreeType const * const mIndexTree;
4013  PolygonPoolList * const mPolygonPoolList;
4014  bool const mInvertSurfaceOrientation;
4015 }; // struct ComputePolygons
4016 
4017 
4018 template<typename TreeType, typename PrimBuilder>
4019 ComputePolygons<TreeType, PrimBuilder>::ComputePolygons(
4020  const std::vector<Int16LeafNodeType*>& signFlagsLeafNodes,
4021  const Int16TreeType& signFlagsTree,
4022  const Index32TreeType& idxTree,
4023  PolygonPoolList& polygons,
4024  bool invertSurfaceOrientation)
4025  : mSignFlagsLeafNodes(signFlagsLeafNodes.data())
4026  , mSignFlagsTree(&signFlagsTree)
4027  , mRefSignFlagsTree(nullptr)
4028  , mIndexTree(&idxTree)
4029  , mPolygonPoolList(&polygons)
4030  , mInvertSurfaceOrientation(invertSurfaceOrientation)
4031 {
4032 }
4033 
4034 template<typename InputTreeType, typename PrimBuilder>
4035 void
4036 ComputePolygons<InputTreeType, PrimBuilder>::operator()(const tbb::blocked_range<size_t>& range) const
4037 {
4038  using Int16ValueAccessor = tree::ValueAccessor<const Int16TreeType>;
4039  Int16ValueAccessor signAcc(*mSignFlagsTree);
4040 
4041  tree::ValueAccessor<const Index32TreeType> idxAcc(*mIndexTree);
4042 
4043  const bool invertSurfaceOrientation = mInvertSurfaceOrientation;
4044 
4045  PrimBuilder mesher;
4046  size_t edgeCount;
4047  Coord ijk, origin;
4048 
4049  // reference data
4050  std::unique_ptr<Int16ValueAccessor> refSignAcc;
4051  if (mRefSignFlagsTree) refSignAcc.reset(new Int16ValueAccessor(*mRefSignFlagsTree));
4052 
4053  for (size_t n = range.begin(); n != range.end(); ++n) {
4054 
4055  const Int16LeafNodeType& node = *mSignFlagsLeafNodes[n];
4056  origin = node.origin();
4057 
4058  // Get an upper bound on the number of primitives.
4059  edgeCount = 0;
4060  typename Int16LeafNodeType::ValueOnCIter iter = node.cbeginValueOn();
4061  for (; iter; ++iter) {
4062  if (iter.getValue() & XEDGE) ++edgeCount;
4063  if (iter.getValue() & YEDGE) ++edgeCount;
4064  if (iter.getValue() & ZEDGE) ++edgeCount;
4065  }
4066 
4067  if (edgeCount == 0) continue;
4068 
4069  mesher.init(edgeCount, (*mPolygonPoolList)[n]);
4070 
4071  const Int16LeafNodeType *signleafPt = signAcc.probeConstLeaf(origin);
4072  const Index32LeafNodeType *idxLeafPt = idxAcc.probeConstLeaf(origin);
4073 
4074  if (!signleafPt || !idxLeafPt) continue;
4075 
4076  const Int16LeafNodeType *refSignLeafPt = nullptr;
4077  if (refSignAcc) refSignLeafPt = refSignAcc->probeConstLeaf(origin);
4078 
4079  Vec3i offsets;
4080 
4081  for (iter = node.cbeginValueOn(); iter; ++iter) {
4082  ijk = iter.getCoord();
4083 
4084  const Int16 flags = iter.getValue();
4085  if (!(flags & 0xE00)) continue;
4086 
4087  Int16 refFlags = 0;
4088  if (refSignLeafPt) {
4089  refFlags = refSignLeafPt->getValue(iter.pos());
4090  }
4091 
4092  const uint8_t cell = uint8_t(SIGNS & flags);
4093 
4094  if (sEdgeGroupTable[cell][0] > 1) {
4095  offsets[0] = (sEdgeGroupTable[cell][1] - 1);
4096  offsets[1] = (sEdgeGroupTable[cell][9] - 1);
4097  offsets[2] = (sEdgeGroupTable[cell][4] - 1);
4098  }
4099  else {
4100  offsets.setZero();
4101  }
4102 
4103  if (ijk[0] > origin[0] && ijk[1] > origin[1] && ijk[2] > origin[2]) {
4104  constructPolygons(invertSurfaceOrientation,
4105  flags, refFlags, offsets, ijk, *signleafPt, *idxLeafPt, mesher);
4106  } else {
4107  constructPolygons(invertSurfaceOrientation,
4108  flags, refFlags, offsets, ijk, signAcc, idxAcc, mesher);
4109  }
4110  }
4111 
4112  mesher.done();
4113  }
4114 
4115 } // ComputePolygons::operator()
4116 
4117 
4118 ////////////////////////////////////////
4119 
4120 
4121 template<typename T>
4122 struct CopyArray
4123 {
4124  CopyArray(T * outputArray, const T * inputArray, size_t outputOffset = 0)
4125  : mOutputArray(outputArray), mInputArray(inputArray), mOutputOffset(outputOffset)
4126  {
4127  }
4128 
4129  void operator()(const tbb::blocked_range<size_t>& inputArrayRange) const
4130  {
4131  const size_t offset = mOutputOffset;
4132  for (size_t n = inputArrayRange.begin(), N = inputArrayRange.end(); n < N; ++n) {
4133  mOutputArray[offset + n] = mInputArray[n];
4134  }
4135  }
4136 
4137 private:
4138  T * const mOutputArray;
4139  T const * const mInputArray;
4140  size_t const mOutputOffset;
4141 }; // struct CopyArray
4142 
4143 
4144 struct FlagAndCountQuadsToSubdivide
4145 {
4146  FlagAndCountQuadsToSubdivide(PolygonPoolList& polygons,
4147  const std::vector<uint8_t>& pointFlags,
4148  std::unique_ptr<openvdb::Vec3s[]>& points,
4149  std::unique_ptr<unsigned[]>& numQuadsToDivide)
4150  : mPolygonPoolList(&polygons)
4151  , mPointFlags(pointFlags.data())
4152  , mPoints(points.get())
4153  , mNumQuadsToDivide(numQuadsToDivide.get())
4154  {
4155  }
4156 
4157  void operator()(const tbb::blocked_range<size_t>& range) const
4158  {
4159  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4160 
4161  PolygonPool& polygons = (*mPolygonPoolList)[n];
4162 
4163  unsigned count = 0;
4164 
4165  // count and tag nonplanar seam line quads.
4166  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4167 
4168  char& flags = polygons.quadFlags(i);
4169 
4170  if ((flags & POLYFLAG_FRACTURE_SEAM) && !(flags & POLYFLAG_EXTERIOR)) {
4171 
4172  Vec4I& quad = polygons.quad(i);
4173 
4174  const bool edgePoly = mPointFlags[quad[0]] || mPointFlags[quad[1]]
4175  || mPointFlags[quad[2]] || mPointFlags[quad[3]];
4176 
4177  if (!edgePoly) continue;
4178 
4179  const Vec3s& p0 = mPoints[quad[0]];
4180  const Vec3s& p1 = mPoints[quad[1]];
4181  const Vec3s& p2 = mPoints[quad[2]];
4182  const Vec3s& p3 = mPoints[quad[3]];
4183 
4184  if (!isPlanarQuad(p0, p1, p2, p3, 1e-6f)) {
4185  flags |= POLYFLAG_SUBDIVIDED;
4186  count++;
4187  }
4188  }
4189  }
4190 
4191  mNumQuadsToDivide[n] = count;
4192  }
4193  }
4194 
4195 private:
4196  PolygonPoolList * const mPolygonPoolList;
4197  uint8_t const * const mPointFlags;
4198  Vec3s const * const mPoints;
4199  unsigned * const mNumQuadsToDivide;
4200 }; // struct FlagAndCountQuadsToSubdivide
4201 
4202 
4203 struct SubdivideQuads
4204 {
4205  SubdivideQuads(PolygonPoolList& polygons,
4206  const std::unique_ptr<openvdb::Vec3s[]>& points,
4207  size_t pointCount,
4208  std::unique_ptr<openvdb::Vec3s[]>& centroids,
4209  std::unique_ptr<unsigned[]>& numQuadsToDivide,
4210  std::unique_ptr<unsigned[]>& centroidOffsets)
4211  : mPolygonPoolList(&polygons)
4212  , mPoints(points.get())
4213  , mCentroids(centroids.get())
4214  , mNumQuadsToDivide(numQuadsToDivide.get())
4215  , mCentroidOffsets(centroidOffsets.get())
4216  , mPointCount(pointCount)
4217  {
4218  }
4219 
4220  void operator()(const tbb::blocked_range<size_t>& range) const
4221  {
4222  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4223 
4224  PolygonPool& polygons = (*mPolygonPoolList)[n];
4225 
4226  const size_t nonplanarCount = size_t(mNumQuadsToDivide[n]);
4227 
4228  if (nonplanarCount > 0) {
4229 
4230  PolygonPool tmpPolygons;
4231  tmpPolygons.resetQuads(polygons.numQuads() - nonplanarCount);
4232  tmpPolygons.resetTriangles(polygons.numTriangles() + size_t(4) * nonplanarCount);
4233 
4234  size_t offset = mCentroidOffsets[n];
4235 
4236  size_t triangleIdx = 0;
4237 
4238  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4239 
4240  const char quadFlags = polygons.quadFlags(i);
4241  if (!(quadFlags & POLYFLAG_SUBDIVIDED)) continue;
4242 
4243  unsigned newPointIdx = unsigned(offset + mPointCount);
4244 
4245  openvdb::Vec4I& quad = polygons.quad(i);
4246 
4247  mCentroids[offset] = (mPoints[quad[0]] + mPoints[quad[1]] +
4248  mPoints[quad[2]] + mPoints[quad[3]]) * 0.25f;
4249 
4250  ++offset;
4251 
4252  {
4253  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4254 
4255  triangle[0] = quad[0];
4256  triangle[1] = newPointIdx;
4257  triangle[2] = quad[3];
4258 
4259  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4260  }
4261 
4262  ++triangleIdx;
4263 
4264  {
4265  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4266 
4267  triangle[0] = quad[0];
4268  triangle[1] = quad[1];
4269  triangle[2] = newPointIdx;
4270 
4271  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4272  }
4273 
4274  ++triangleIdx;
4275 
4276  {
4277  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4278 
4279  triangle[0] = quad[1];
4280  triangle[1] = quad[2];
4281  triangle[2] = newPointIdx;
4282 
4283  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4284  }
4285 
4286 
4287  ++triangleIdx;
4288 
4289  {
4290  Vec3I& triangle = tmpPolygons.triangle(triangleIdx);
4291 
4292  triangle[0] = quad[2];
4293  triangle[1] = quad[3];
4294  triangle[2] = newPointIdx;
4295 
4296  tmpPolygons.triangleFlags(triangleIdx) = quadFlags;
4297  }
4298 
4299  ++triangleIdx;
4300 
4301  quad[0] = util::INVALID_IDX; // mark for deletion
4302  }
4303 
4304 
4305  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4306  tmpPolygons.triangle(triangleIdx) = polygons.triangle(i);
4307  tmpPolygons.triangleFlags(triangleIdx) = polygons.triangleFlags(i);
4308  ++triangleIdx;
4309  }
4310 
4311  size_t quadIdx = 0;
4312  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4313  openvdb::Vec4I& quad = polygons.quad(i);
4314 
4315  if (quad[0] != util::INVALID_IDX) { // ignore invalid quads
4316  tmpPolygons.quad(quadIdx) = quad;
4317  tmpPolygons.quadFlags(quadIdx) = polygons.quadFlags(i);
4318  ++quadIdx;
4319  }
4320  }
4321 
4322  polygons.copy(tmpPolygons);
4323  }
4324  }
4325  }
4326 
4327 private:
4328  PolygonPoolList * const mPolygonPoolList;
4329  Vec3s const * const mPoints;
4330  Vec3s * const mCentroids;
4331  unsigned * const mNumQuadsToDivide;
4332  unsigned * const mCentroidOffsets;
4333  size_t const mPointCount;
4334 }; // struct SubdivideQuads
4335 
4336 
4337 inline void
4338 subdivideNonplanarSeamLineQuads(
4339  PolygonPoolList& polygonPoolList,
4340  size_t polygonPoolListSize,
4341  PointList& pointList,
4342  size_t& pointListSize,
4343  std::vector<uint8_t>& pointFlags)
4344 {
4345  const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4346 
4347  std::unique_ptr<unsigned[]> numQuadsToDivide(new unsigned[polygonPoolListSize]);
4348 
4349  tbb::parallel_for(polygonPoolListRange,
4350  FlagAndCountQuadsToSubdivide(polygonPoolList, pointFlags, pointList, numQuadsToDivide));
4351 
4352  std::unique_ptr<unsigned[]> centroidOffsets(new unsigned[polygonPoolListSize]);
4353 
4354  size_t centroidCount = 0;
4355 
4356  {
4357  unsigned sum = 0;
4358  for (size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4359  centroidOffsets[n] = sum;
4360  sum += numQuadsToDivide[n];
4361  }
4362  centroidCount = size_t(sum);
4363  }
4364 
4365  std::unique_ptr<Vec3s[]> centroidList(new Vec3s[centroidCount]);
4366 
4367  tbb::parallel_for(polygonPoolListRange,
4368  SubdivideQuads(polygonPoolList, pointList, pointListSize,
4369  centroidList, numQuadsToDivide, centroidOffsets));
4370 
4371  if (centroidCount > 0) {
4372 
4373  const size_t newPointListSize = centroidCount + pointListSize;
4374 
4375  std::unique_ptr<openvdb::Vec3s[]> newPointList(new openvdb::Vec3s[newPointListSize]);
4376 
4377  tbb::parallel_for(tbb::blocked_range<size_t>(0, pointListSize),
4378  CopyArray<Vec3s>(newPointList.get(), pointList.get()));
4379 
4380  tbb::parallel_for(tbb::blocked_range<size_t>(0, newPointListSize - pointListSize),
4381  CopyArray<Vec3s>(newPointList.get(), centroidList.get(), pointListSize));
4382 
4383  pointListSize = newPointListSize;
4384  pointList.swap(newPointList);
4385  pointFlags.resize(pointListSize, 0);
4386  }
4387 }
4388 
4389 
4390 struct ReviseSeamLineFlags
4391 {
4392  ReviseSeamLineFlags(PolygonPoolList& polygons,
4393  const std::vector<uint8_t>& pointFlags)
4394  : mPolygonPoolList(&polygons)
4395  , mPointFlags(pointFlags.data())
4396  {
4397  }
4398 
4399  void operator()(const tbb::blocked_range<size_t>& range) const
4400  {
4401  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4402 
4403  PolygonPool& polygons = (*mPolygonPoolList)[n];
4404 
4405  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4406 
4407  char& flags = polygons.quadFlags(i);
4408 
4409  if (flags & POLYFLAG_FRACTURE_SEAM) {
4410 
4411  openvdb::Vec4I& verts = polygons.quad(i);
4412 
4413  const bool hasSeamLinePoint =
4414  mPointFlags[verts[0]] || mPointFlags[verts[1]] ||
4415  mPointFlags[verts[2]] || mPointFlags[verts[3]];
4416 
4417  if (!hasSeamLinePoint) {
4418  flags &= ~POLYFLAG_FRACTURE_SEAM;
4419  }
4420  }
4421  } // end quad loop
4422 
4423  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4424 
4425  char& flags = polygons.triangleFlags(i);
4426 
4427  if (flags & POLYFLAG_FRACTURE_SEAM) {
4428 
4429  openvdb::Vec3I& verts = polygons.triangle(i);
4430 
4431  const bool hasSeamLinePoint =
4432  mPointFlags[verts[0]] || mPointFlags[verts[1]] || mPointFlags[verts[2]];
4433 
4434  if (!hasSeamLinePoint) {
4435  flags &= ~POLYFLAG_FRACTURE_SEAM;
4436  }
4437 
4438  }
4439  } // end triangle loop
4440 
4441  } // end polygon pool loop
4442  }
4443 
4444 private:
4445  PolygonPoolList * const mPolygonPoolList;
4446  uint8_t const * const mPointFlags;
4447 }; // struct ReviseSeamLineFlags
4448 
4449 
4450 inline void
4451 reviseSeamLineFlags(PolygonPoolList& polygonPoolList, size_t polygonPoolListSize,
4452  std::vector<uint8_t>& pointFlags)
4453 {
4454  tbb::parallel_for(tbb::blocked_range<size_t>(0, polygonPoolListSize),
4455  ReviseSeamLineFlags(polygonPoolList, pointFlags));
4456 }
4457 
4458 
4459 ////////////////////////////////////////
4460 
4461 
4462 template<typename InputTreeType>
4463 struct MaskDisorientedTrianglePoints
4464 {
4465  MaskDisorientedTrianglePoints(const InputTreeType& inputTree, const PolygonPoolList& polygons,
4466  const PointList& pointList, std::unique_ptr<uint8_t[]>& pointMask,
4467  const math::Transform& transform, bool invertSurfaceOrientation)
4468  : mInputTree(&inputTree)
4469  , mPolygonPoolList(&polygons)
4470  , mPointList(&pointList)
4471  , mPointMask(pointMask.get())
4472  , mTransform(transform)
4473  , mInvertSurfaceOrientation(invertSurfaceOrientation)
4474  {
4475  }
4476 
4477  void operator()(const tbb::blocked_range<size_t>& range) const
4478  {
4480 
4481  tree::ValueAccessor<const InputTreeType> inputAcc(*mInputTree);
4482  Vec3s centroid, normal;
4483  Coord ijk;
4484 
4485  const bool invertGradientDir = mInvertSurfaceOrientation || isBoolValue<ValueType>();
4486 
4487  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
4488 
4489  const PolygonPool& polygons = (*mPolygonPoolList)[n];
4490 
4491  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4492 
4493  const Vec3I& verts = polygons.triangle(i);
4494 
4495  const Vec3s& v0 = (*mPointList)[verts[0]];
4496  const Vec3s& v1 = (*mPointList)[verts[1]];
4497  const Vec3s& v2 = (*mPointList)[verts[2]];
4498 
4499  normal = (v2 - v0).cross((v1 - v0));
4500  normal.normalize();
4501 
4502  centroid = (v0 + v1 + v2) * (1.0f / 3.0f);
4503  ijk = mTransform.worldToIndexCellCentered(centroid);
4504 
4505  Vec3s dir( math::ISGradient<math::CD_2ND>::result(inputAcc, ijk) );
4506  dir.normalize();
4507 
4508  if (invertGradientDir) {
4509  dir = -dir;
4510  }
4511 
4512  // check if the angle is obtuse
4513  if (dir.dot(normal) < -0.5f) {
4514  // Concurrent writes to same memory address can occur, but
4515  // all threads are writing the same value and char is atomic.
4516  // (It is extremely rare that disoriented triangles share points,
4517  // false sharing related performance impacts are not a concern.)
4518  mPointMask[verts[0]] = 1;
4519  mPointMask[verts[1]] = 1;
4520  mPointMask[verts[2]] = 1;
4521  }
4522 
4523  } // end triangle loop
4524 
4525  } // end polygon pool loop
4526  }
4527 
4528 private:
4529  InputTreeType const * const mInputTree;
4530  PolygonPoolList const * const mPolygonPoolList;
4531  PointList const * const mPointList;
4532  uint8_t * const mPointMask;
4533  math::Transform const mTransform;
4534  bool const mInvertSurfaceOrientation;
4535 }; // struct MaskDisorientedTrianglePoints
4536 
4537 
4538 template<typename InputTree>
4539 inline void
4540 relaxDisorientedTriangles(
4541  bool invertSurfaceOrientation,
4542  const InputTree& inputTree,
4543  const math::Transform& transform,
4544  PolygonPoolList& polygonPoolList,
4545  size_t polygonPoolListSize,
4546  PointList& pointList,
4547  const size_t pointListSize)
4548 {
4549  const tbb::blocked_range<size_t> polygonPoolListRange(0, polygonPoolListSize);
4550 
4551  std::unique_ptr<uint8_t[]> pointMask(new uint8_t[pointListSize]);
4552  fillArray(pointMask.get(), uint8_t(0), pointListSize);
4553 
4554  tbb::parallel_for(polygonPoolListRange,
4555  MaskDisorientedTrianglePoints<InputTree>(
4556  inputTree, polygonPoolList, pointList, pointMask, transform, invertSurfaceOrientation));
4557 
4558  std::unique_ptr<uint8_t[]> pointUpdates(new uint8_t[pointListSize]);
4559  fillArray(pointUpdates.get(), uint8_t(0), pointListSize);
4560 
4561  std::unique_ptr<Vec3s[]> newPoints(new Vec3s[pointListSize]);
4562  fillArray(newPoints.get(), Vec3s(0.0f, 0.0f, 0.0f), pointListSize);
4563 
4564  for (size_t n = 0, N = polygonPoolListSize; n < N; ++n) {
4565 
4566  PolygonPool& polygons = polygonPoolList[n];
4567 
4568  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4569  openvdb::Vec4I& verts = polygons.quad(i);
4570 
4571  for (int v = 0; v < 4; ++v) {
4572 
4573  const unsigned pointIdx = verts[v];
4574 
4575  if (pointMask[pointIdx] == 1) {
4576 
4577  newPoints[pointIdx] +=
4578  pointList[verts[0]] + pointList[verts[1]] +
4579  pointList[verts[2]] + pointList[verts[3]];
4580 
4581  pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 4);
4582  }
4583  }
4584  }
4585 
4586  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4587  openvdb::Vec3I& verts = polygons.triangle(i);
4588 
4589  for (int v = 0; v < 3; ++v) {
4590 
4591  const unsigned pointIdx = verts[v];
4592 
4593  if (pointMask[pointIdx] == 1) {
4594  newPoints[pointIdx] +=
4595  pointList[verts[0]] + pointList[verts[1]] + pointList[verts[2]];
4596 
4597  pointUpdates[pointIdx] = uint8_t(pointUpdates[pointIdx] + 3);
4598  }
4599  }
4600  }
4601  }
4602 
4603  for (size_t n = 0, N = pointListSize; n < N; ++n) {
4604  if (pointUpdates[n] > 0) {
4605  const double weight = 1.0 / double(pointUpdates[n]);
4606  pointList[n] = newPoints[n] * float(weight);
4607  }
4608  }
4609 }
4610 
4611 
4612 template<typename GridType>
4613 void
4614 doVolumeToMesh(
4615  const GridType& grid,
4616  std::vector<Vec3s>& points,
4617  std::vector<Vec3I>& triangles,
4618  std::vector<Vec4I>& quads,
4619  double isovalue,
4620  double adaptivity,
4621  bool relaxDisorientedTriangles)
4622 {
4624  "volume to mesh conversion is supported only for scalar grids");
4625 
4626  VolumeToMesh mesher(isovalue, adaptivity, relaxDisorientedTriangles);
4627  mesher(grid);
4628 
4629  // Preallocate the point list
4630  points.clear();
4631  points.resize(mesher.pointListSize());
4632 
4633  { // Copy points
4634  volume_to_mesh_internal::PointListCopy ptnCpy(mesher.pointList(), points);
4635  tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
4636  mesher.pointList().reset(nullptr);
4637  }
4638 
4639  PolygonPoolList& polygonPoolList = mesher.polygonPoolList();
4640 
4641  { // Preallocate primitive lists
4642  size_t numQuads = 0, numTriangles = 0;
4643  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4644  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4645  numTriangles += polygons.numTriangles();
4646  numQuads += polygons.numQuads();
4647  }
4648 
4649  triangles.clear();
4650  triangles.resize(numTriangles);
4651  quads.clear();
4652  quads.resize(numQuads);
4653  }
4654 
4655  // Copy primitives
4656  size_t qIdx = 0, tIdx = 0;
4657  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
4658  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
4659 
4660  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
4661  quads[qIdx++] = polygons.quad(i);
4662  }
4663 
4664  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
4665  triangles[tIdx++] = polygons.triangle(i);
4666  }
4667  }
4668 }
4669 
4670 
4671 } // volume_to_mesh_internal namespace
4672 
4673 /// @endcond
4674 
4675 ////////////////////////////////////////
4676 
4677 
4678 inline
4679 PolygonPool::PolygonPool()
4680  : mNumQuads(0)
4681  , mNumTriangles(0)
4682  , mQuads(nullptr)
4683  , mTriangles(nullptr)
4684  , mQuadFlags(nullptr)
4685  , mTriangleFlags(nullptr)
4686 {
4687 }
4688 
4689 
4690 inline
4691 PolygonPool::PolygonPool(const size_t numQuads, const size_t numTriangles)
4692  : mNumQuads(numQuads)
4693  , mNumTriangles(numTriangles)
4694  , mQuads(new openvdb::Vec4I[mNumQuads])
4695  , mTriangles(new openvdb::Vec3I[mNumTriangles])
4696  , mQuadFlags(new char[mNumQuads])
4697  , mTriangleFlags(new char[mNumTriangles])
4698 {
4699 }
4700 
4701 
4702 inline void
4704 {
4705  resetQuads(rhs.numQuads());
4707 
4708  for (size_t i = 0; i < mNumQuads; ++i) {
4709  mQuads[i] = rhs.mQuads[i];
4710  mQuadFlags[i] = rhs.mQuadFlags[i];
4711  }
4712 
4713  for (size_t i = 0; i < mNumTriangles; ++i) {
4714  mTriangles[i] = rhs.mTriangles[i];
4715  mTriangleFlags[i] = rhs.mTriangleFlags[i];
4716  }
4717 }
4718 
4719 
4720 inline void
4722 {
4723  mNumQuads = size;
4724  mQuads.reset(new openvdb::Vec4I[mNumQuads]);
4725  mQuadFlags.reset(new char[mNumQuads]);
4726 }
4727 
4728 
4729 inline void
4731 {
4732  mNumQuads = 0;
4733  mQuads.reset(nullptr);
4734  mQuadFlags.reset(nullptr);
4735 }
4736 
4737 
4738 inline void
4740 {
4741  mNumTriangles = size;
4742  mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
4743  mTriangleFlags.reset(new char[mNumTriangles]);
4744 }
4745 
4746 
4747 inline void
4749 {
4750  mNumTriangles = 0;
4751  mTriangles.reset(nullptr);
4752  mTriangleFlags.reset(nullptr);
4753 }
4754 
4755 
4756 inline bool
4757 PolygonPool::trimQuads(const size_t n, bool reallocate)
4758 {
4759  if (!(n < mNumQuads)) return false;
4760 
4761  if (reallocate) {
4762 
4763  if (n == 0) {
4764  mQuads.reset(nullptr);
4765  } else {
4766 
4767  std::unique_ptr<openvdb::Vec4I[]> quads(new openvdb::Vec4I[n]);
4768  std::unique_ptr<char[]> flags(new char[n]);
4769 
4770  for (size_t i = 0; i < n; ++i) {
4771  quads[i] = mQuads[i];
4772  flags[i] = mQuadFlags[i];
4773  }
4774 
4775  mQuads.swap(quads);
4776  mQuadFlags.swap(flags);
4777  }
4778  }
4779 
4780  mNumQuads = n;
4781  return true;
4782 }
4783 
4784 
4785 inline bool
4786 PolygonPool::trimTrinagles(const size_t n, bool reallocate)
4787 {
4788  if (!(n < mNumTriangles)) return false;
4789 
4790  if (reallocate) {
4791 
4792  if (n == 0) {
4793  mTriangles.reset(nullptr);
4794  } else {
4795 
4796  std::unique_ptr<openvdb::Vec3I[]> triangles(new openvdb::Vec3I[n]);
4797  std::unique_ptr<char[]> flags(new char[n]);
4798 
4799  for (size_t i = 0; i < n; ++i) {
4800  triangles[i] = mTriangles[i];
4801  flags[i] = mTriangleFlags[i];
4802  }
4803 
4804  mTriangles.swap(triangles);
4805  mTriangleFlags.swap(flags);
4806  }
4807  }
4808 
4809  mNumTriangles = n;
4810  return true;
4811 }
4812 
4813 
4814 ////////////////////////////////////////
4815 
4816 
4817 inline
4818 VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity, bool relaxDisorientedTriangles)
4819  : mPoints(nullptr)
4820  , mPolygons()
4821  , mPointListSize(0)
4822  , mSeamPointListSize(0)
4823  , mPolygonPoolListSize(0)
4824  , mIsovalue(isovalue)
4825  , mPrimAdaptivity(adaptivity)
4826  , mSecAdaptivity(0.0)
4827  , mRefGrid(GridBase::ConstPtr())
4828  , mSurfaceMaskGrid(GridBase::ConstPtr())
4829  , mAdaptivityGrid(GridBase::ConstPtr())
4830  , mAdaptivityMaskTree(TreeBase::ConstPtr())
4831  , mRefSignTree(TreeBase::Ptr())
4832  , mRefIdxTree(TreeBase::Ptr())
4833  , mInvertSurfaceMask(false)
4834  , mRelaxDisorientedTriangles(relaxDisorientedTriangles)
4835  , mQuantizedSeamPoints(nullptr)
4836  , mPointFlags(0)
4837 {
4838 }
4839 
4840 
4841 inline void
4842 VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity)
4843 {
4844  mRefGrid = grid;
4845  mSecAdaptivity = secAdaptivity;
4846 
4847  // Clear out old auxiliary data
4848  mRefSignTree = TreeBase::Ptr();
4849  mRefIdxTree = TreeBase::Ptr();
4850  mSeamPointListSize = 0;
4851  mQuantizedSeamPoints.reset(nullptr);
4852 }
4853 
4854 
4855 inline void
4857 {
4858  mSurfaceMaskGrid = mask;
4859  mInvertSurfaceMask = invertMask;
4860 }
4861 
4862 
4863 inline void
4865 {
4866  mAdaptivityGrid = grid;
4867 }
4868 
4869 
4870 inline void
4872 {
4873  mAdaptivityMaskTree = tree;
4874 }
4875 
4876 
4877 template<typename InputGridType>
4878 inline void
4879 VolumeToMesh::operator()(const InputGridType& inputGrid)
4880 {
4881  // input data types
4882 
4883  using InputTreeType = typename InputGridType::TreeType;
4884  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
4885  using InputValueType = typename InputLeafNodeType::ValueType;
4886 
4887  // auxiliary data types
4888 
4889  using FloatTreeType = typename InputTreeType::template ValueConverter<float>::Type;
4890  using FloatGridType = Grid<FloatTreeType>;
4891  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
4892  using Int16TreeType = typename InputTreeType::template ValueConverter<Int16>::Type;
4893  using Int16LeafNodeType = typename Int16TreeType::LeafNodeType;
4894  using Index32TreeType = typename InputTreeType::template ValueConverter<Index32>::Type;
4895  using Index32LeafNodeType = typename Index32TreeType::LeafNodeType;
4896 
4897  // clear old data
4898  mPointListSize = 0;
4899  mPoints.reset();
4900  mPolygonPoolListSize = 0;
4901  mPolygons.reset();
4902  mPointFlags.clear();
4903 
4904  // settings
4905 
4906  const math::Transform& transform = inputGrid.transform();
4907  const InputValueType isovalue = InputValueType(mIsovalue);
4908  const float adaptivityThreshold = float(mPrimAdaptivity);
4909  const bool adaptive = mPrimAdaptivity > 1e-7 || mSecAdaptivity > 1e-7;
4910 
4911  // The default surface orientation is setup for level set and bool/mask grids.
4912  // Boolean grids are handled correctly by their value type. Signed distance fields,
4913  // unsigned distance fields and fog volumes have the same value type but use different
4914  // inside value classifications.
4915  const bool invertSurfaceOrientation = (!volume_to_mesh_internal::isBoolValue<InputValueType>()
4916  && (inputGrid.getGridClass() != openvdb::GRID_LEVEL_SET));
4917 
4918  // references, masks and auxiliary data
4919 
4920  const InputTreeType& inputTree = inputGrid.tree();
4921 
4922  BoolTreeType intersectionTree(false), adaptivityMask(false);
4923 
4924  if (mAdaptivityMaskTree && mAdaptivityMaskTree->type() == BoolTreeType::treeType()) {
4925  const BoolTreeType *refAdaptivityMask=
4926  static_cast<const BoolTreeType*>(mAdaptivityMaskTree.get());
4927  adaptivityMask.topologyUnion(*refAdaptivityMask);
4928  }
4929 
4930  Int16TreeType signFlagsTree(0);
4931  Index32TreeType pointIndexTree(std::numeric_limits<Index32>::max());
4932 
4933 
4934  // collect auxiliary data
4935 
4936  volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4937  intersectionTree, inputTree, isovalue);
4938 
4939  volume_to_mesh_internal::applySurfaceMask(intersectionTree, adaptivityMask,
4940  inputGrid, mSurfaceMaskGrid, mInvertSurfaceMask, isovalue);
4941 
4942  if (intersectionTree.empty()) return;
4943 
4944  volume_to_mesh_internal::computeAuxiliaryData(
4945  signFlagsTree, pointIndexTree, intersectionTree, inputTree, isovalue);
4946 
4947  intersectionTree.clear();
4948 
4949  std::vector<Index32LeafNodeType*> pointIndexLeafNodes;
4950  pointIndexTree.getNodes(pointIndexLeafNodes);
4951 
4952  std::vector<Int16LeafNodeType*> signFlagsLeafNodes;
4953  signFlagsTree.getNodes(signFlagsLeafNodes);
4954 
4955  const tbb::blocked_range<size_t> auxiliaryLeafNodeRange(0, signFlagsLeafNodes.size());
4956 
4957 
4958  // optionally collect auxiliary data from a reference volume.
4959 
4960  Int16TreeType* refSignFlagsTree = nullptr;
4961  Index32TreeType* refPointIndexTree = nullptr;
4962  InputTreeType const* refInputTree = nullptr;
4963 
4964  if (mRefGrid && mRefGrid->type() == InputGridType::gridType()) {
4965 
4966  const InputGridType* refGrid = static_cast<const InputGridType*>(mRefGrid.get());
4967  refInputTree = &refGrid->tree();
4968 
4969  if (!mRefSignTree && !mRefIdxTree) {
4970 
4971  // first time, collect and cache auxiliary data.
4972 
4973  typename Int16TreeType::Ptr refSignFlagsTreePt(new Int16TreeType(0));
4974  typename Index32TreeType::Ptr refPointIndexTreePt(
4975  new Index32TreeType(std::numeric_limits<Index32>::max()));
4976 
4977  BoolTreeType refIntersectionTree(false);
4978 
4979  volume_to_mesh_internal::identifySurfaceIntersectingVoxels(
4980  refIntersectionTree, *refInputTree, isovalue);
4981 
4982  volume_to_mesh_internal::computeAuxiliaryData(*refSignFlagsTreePt,
4983  *refPointIndexTreePt, refIntersectionTree, *refInputTree, isovalue);
4984 
4985  mRefSignTree = refSignFlagsTreePt;
4986  mRefIdxTree = refPointIndexTreePt;
4987  }
4988 
4989  if (mRefSignTree && mRefIdxTree) {
4990 
4991  // get cached auxiliary data
4992 
4993  refSignFlagsTree = static_cast<Int16TreeType*>(mRefSignTree.get());
4994  refPointIndexTree = static_cast<Index32TreeType*>(mRefIdxTree.get());
4995  }
4996 
4997 
4998  if (refSignFlagsTree && refPointIndexTree) {
4999 
5000  // generate seam line sample points
5001 
5002  volume_to_mesh_internal::markSeamLineData(signFlagsTree, *refSignFlagsTree);
5003 
5004  if (mSeamPointListSize == 0) {
5005 
5006  // count unique points on reference surface
5007 
5008  std::vector<Int16LeafNodeType*> refSignFlagsLeafNodes;
5009  refSignFlagsTree->getNodes(refSignFlagsLeafNodes);
5010 
5011  std::unique_ptr<Index32[]> leafNodeOffsets(
5012  new Index32[refSignFlagsLeafNodes.size()]);
5013 
5014  tbb::parallel_for(tbb::blocked_range<size_t>(0, refSignFlagsLeafNodes.size()),
5015  volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>(
5016  refSignFlagsLeafNodes, leafNodeOffsets));
5017 
5018  {
5019  Index32 count = 0;
5020  for (size_t n = 0, N = refSignFlagsLeafNodes.size(); n < N; ++n) {
5021  const Index32 tmp = leafNodeOffsets[n];
5022  leafNodeOffsets[n] = count;
5023  count += tmp;
5024  }
5025  mSeamPointListSize = size_t(count);
5026  }
5027 
5028  if (mSeamPointListSize != 0) {
5029 
5030  mQuantizedSeamPoints.reset(new uint32_t[mSeamPointListSize]);
5031 
5032  std::memset(mQuantizedSeamPoints.get(), 0, sizeof(uint32_t) * mSeamPointListSize);
5033 
5034  std::vector<Index32LeafNodeType*> refPointIndexLeafNodes;
5035  refPointIndexTree->getNodes(refPointIndexLeafNodes);
5036 
5037  tbb::parallel_for(tbb::blocked_range<size_t>(0, refPointIndexLeafNodes.size()),
5038  volume_to_mesh_internal::MapPoints<Index32LeafNodeType>(
5039  refPointIndexLeafNodes, refSignFlagsLeafNodes, leafNodeOffsets));
5040  }
5041  }
5042 
5043  if (mSeamPointListSize != 0) {
5044 
5045  tbb::parallel_for(auxiliaryLeafNodeRange,
5046  volume_to_mesh_internal::SeamLineWeights<InputTreeType>(
5047  signFlagsLeafNodes, inputTree, *refPointIndexTree, *refSignFlagsTree,
5048  mQuantizedSeamPoints.get(), isovalue));
5049  }
5050  }
5051  }
5052 
5053  const bool referenceMeshing = refSignFlagsTree && refPointIndexTree && refInputTree;
5054 
5055 
5056  // adapt and count unique points
5057 
5058  std::unique_ptr<Index32[]> leafNodeOffsets(new Index32[signFlagsLeafNodes.size()]);
5059 
5060  if (adaptive) {
5061  volume_to_mesh_internal::MergeVoxelRegions<InputGridType> mergeOp(
5062  inputGrid, pointIndexTree, pointIndexLeafNodes, signFlagsLeafNodes,
5063  isovalue, adaptivityThreshold, invertSurfaceOrientation);
5064 
5065  if (mAdaptivityGrid && mAdaptivityGrid->type() == FloatGridType::gridType()) {
5066  const FloatGridType* adaptivityGrid =
5067  static_cast<const FloatGridType*>(mAdaptivityGrid.get());
5068  mergeOp.setSpatialAdaptivity(*adaptivityGrid);
5069  }
5070 
5071  if (!adaptivityMask.empty()) {
5072  mergeOp.setAdaptivityMask(adaptivityMask);
5073  }
5074 
5075  if (referenceMeshing) {
5076  mergeOp.setRefSignFlagsData(*refSignFlagsTree, float(mSecAdaptivity));
5077  }
5078 
5079  tbb::parallel_for(auxiliaryLeafNodeRange, mergeOp);
5080 
5081  volume_to_mesh_internal::AdaptiveLeafNodePointCount<Index32LeafNodeType>
5082  op(pointIndexLeafNodes, signFlagsLeafNodes, leafNodeOffsets);
5083 
5084  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5085 
5086  } else {
5087 
5088  volume_to_mesh_internal::LeafNodePointCount<Int16LeafNodeType::LOG2DIM>
5089  op(signFlagsLeafNodes, leafNodeOffsets);
5090 
5091  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5092  }
5093 
5094 
5095  {
5096  Index32 pointCount = 0;
5097  for (size_t n = 0, N = signFlagsLeafNodes.size(); n < N; ++n) {
5098  const Index32 tmp = leafNodeOffsets[n];
5099  leafNodeOffsets[n] = pointCount;
5100  pointCount += tmp;
5101  }
5102 
5103  mPointListSize = size_t(pointCount);
5104  mPoints.reset(new openvdb::Vec3s[mPointListSize]);
5105  mPointFlags.clear();
5106  }
5107 
5108 
5109  // compute points
5110 
5111  {
5112  volume_to_mesh_internal::ComputePoints<InputTreeType>
5113  op(mPoints.get(), inputTree, pointIndexLeafNodes,
5114  signFlagsLeafNodes, leafNodeOffsets, transform, mIsovalue);
5115 
5116  if (referenceMeshing) {
5117  mPointFlags.resize(mPointListSize);
5118  op.setRefData(*refInputTree, *refPointIndexTree, *refSignFlagsTree,
5119  mQuantizedSeamPoints.get(), mPointFlags.data());
5120  }
5121 
5122  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5123  }
5124 
5125 
5126  // compute polygons
5127 
5128  mPolygonPoolListSize = signFlagsLeafNodes.size();
5129  mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
5130 
5131  if (adaptive) {
5132 
5133  using PrimBuilder = volume_to_mesh_internal::AdaptivePrimBuilder;
5134 
5135  volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5136  op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5137  mPolygons, invertSurfaceOrientation);
5138 
5139  if (referenceMeshing) {
5140  op.setRefSignTree(refSignFlagsTree);
5141  }
5142 
5143  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5144 
5145  } else {
5146 
5147  using PrimBuilder = volume_to_mesh_internal::UniformPrimBuilder;
5148 
5149  volume_to_mesh_internal::ComputePolygons<Int16TreeType, PrimBuilder>
5150  op(signFlagsLeafNodes, signFlagsTree, pointIndexTree,
5151  mPolygons, invertSurfaceOrientation);
5152 
5153  if (referenceMeshing) {
5154  op.setRefSignTree(refSignFlagsTree);
5155  }
5156 
5157  tbb::parallel_for(auxiliaryLeafNodeRange, op);
5158  }
5159 
5160 
5161  signFlagsTree.clear();
5162  pointIndexTree.clear();
5163 
5164 
5165  if (adaptive && mRelaxDisorientedTriangles) {
5166  volume_to_mesh_internal::relaxDisorientedTriangles(invertSurfaceOrientation,
5167  inputTree, transform, mPolygons, mPolygonPoolListSize, mPoints, mPointListSize);
5168  }
5169 
5170 
5171  if (referenceMeshing) {
5172  volume_to_mesh_internal::subdivideNonplanarSeamLineQuads(
5173  mPolygons, mPolygonPoolListSize, mPoints, mPointListSize, mPointFlags);
5174 
5175  volume_to_mesh_internal::reviseSeamLineFlags(mPolygons, mPolygonPoolListSize, mPointFlags);
5176  }
5177 
5178 }
5179 
5180 
5181 ////////////////////////////////////////
5182 
5183 
5184 template<typename GridType>
5186  const GridType& grid,
5187  std::vector<Vec3s>& points,
5188  std::vector<Vec3I>& triangles,
5189  std::vector<Vec4I>& quads,
5190  double isovalue,
5191  double adaptivity,
5192  bool relaxDisorientedTriangles)
5193 {
5194  volume_to_mesh_internal::doVolumeToMesh(grid, points, triangles, quads,
5195  isovalue, adaptivity, relaxDisorientedTriangles);
5196 }
5197 
5198 template<typename GridType>
5200  const GridType& grid,
5201  std::vector<Vec3s>& points,
5202  std::vector<Vec4I>& quads,
5203  double isovalue)
5204 {
5205  std::vector<Vec3I> triangles;
5206  volumeToMesh(grid, points, triangles, quads, isovalue, 0.0, true);
5207 }
5208 
5209 
5210 ////////////////////////////////////////
5211 
5212 
5213 // Explicit Template Instantiation
5214 
5215 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
5216 
5217 #ifdef OPENVDB_INSTANTIATE_VOLUMETOMESH
5219 #endif
5220 
5221 #define _FUNCTION(TreeT) \
5222  void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec4I>&, double)
5224 #undef _FUNCTION
5225 
5226 #define _FUNCTION(TreeT) \
5227  void volumeToMesh(const Grid<TreeT>&, std::vector<Vec3s>&, std::vector<Vec3I>&, std::vector<Vec4I>&, double, double, bool)
5229 #undef _FUNCTION
5230 
5231 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
5232 
5233 
5234 } // namespace tools
5235 } // namespace OPENVDB_VERSION_NAME
5236 } // namespace openvdb
5237 
5238 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
std::unique_ptr< openvdb::Vec3s[]> PointList
Point and primitive list types.
Definition: VolumeToMesh.h:160
void operator()(const InputGridType &)
Main call.
const char & triangleFlags(size_t n) const
Definition: VolumeToMesh.h:138
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
GU_API exint quads(GU_Detail *gdp, UT_Array< GA_OffsetArray > &rings, UT_Array< GA_OffsetArray > &originalRings, GA_PrimitiveGroup *patchgroup, GA_PrimitiveGroup *loopgroup, bool smooth, fpreal smoothstrength, bool edgeloop, fpreal edgeloopPercentage)
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
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
__hostdev__ uint64_t pointCount() const
Definition: NanoVDB.h:5974
GLbitfield flags
Definition: glcorearb.h:1596
GLenum GLint * range
Definition: glcorearb.h:1925
GLdouble GLdouble GLint GLint const GLdouble * points
Definition: glad.h:2676
constexpr Index32 INVALID_IDX
Definition: Util.h:19
const PolygonPoolList & polygonPoolList() const
Definition: VolumeToMesh.h:189
GridType
List of types that are currently supported by NanoVDB.
Definition: NanoVDB.h:294
GLboolean invert
Definition: glcorearb.h:549
GLboolean * data
Definition: glcorearb.h:131
const GLdouble * v
Definition: glcorearb.h:837
GLuint start
Definition: glcorearb.h:475
GLsizei const GLfloat * value
Definition: glcorearb.h:824
bool trimQuads(const size_t n, bool reallocate=false)
**And then you can **find out if it s done
Definition: thread.h:622
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
void reverse(I begin, I end)
Definition: pugixml.cpp:7190
GLdouble s
Definition: glad.h:3009
const std::vector< uint8_t > & pointFlags() const
Definition: VolumeToMesh.h:192
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GLint y
Definition: glcorearb.h:103
MaskT< LOG2DIM > mMask
Definition: NanoVDB.h:5924
**But if you need a result
Definition: thread.h:613
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:818
uint64 value_type
Definition: GA_PrimCompat.h:29
void setSurfaceMask(const GridBase::ConstPtr &mask, bool invertMask=false)
std::unique_ptr< PolygonPool[]> PolygonPoolList
Point and primitive list types.
Definition: VolumeToMesh.h:161
GLuint GLsizei const GLuint const GLintptr * offsets
Definition: glcorearb.h:2621
SIM_DerVector3 normalize() const
Abstract base class for typed grids.
Definition: Grid.h:77
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:454
SYS_FORCE_INLINE const_iterator end() const
Vec3< double > Vec3d
Definition: NanoVDB.h:1685
GLdouble n
Definition: glcorearb.h:2008
GLfloat f
Definition: glcorearb.h:1926
GLintptr offset
Definition: glcorearb.h:665
Definition: core.h:760
IMATH_NAMESPACE::V2f float
math::Vec4< Index32 > Vec4I
Definition: Types.h:88
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...
bool diagonalizeSymmetricMatrix(const Mat3< T > &input, Mat3< T > &Q, Vec3< T > &D, unsigned int MAX_ITERATIONS=250)
Use Jacobi iterations to decompose a symmetric 3x3 matrix (diagonalize and compute eigenvectors) ...
Definition: Mat3.h:737
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
BBox< Coord > CoordBBox
Definition: NanoVDB.h:2516
uint64_t mPointCount
Definition: NanoVDB.h:5960
GLuint GLuint end
Definition: glcorearb.h:475
GLint GLuint mask
Definition: glcorearb.h:124
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
Definition: UT_ARTMap.h:1073
GLuint id
Definition: glcorearb.h:655
void setAdaptivityMask(const TreeBase::ConstPtr &tree)
Mesh any scalar grid that has a continuous isosurface.
Definition: VolumeToMesh.h:169
const openvdb::Vec4I & quad(size_t n) const
Definition: VolumeToMesh.h:123
GA_API const UT_StringHolder transform
GLint GLenum GLint x
Definition: glcorearb.h:409
bool trimTrinagles(const size_t n, bool reallocate=false)
Vec3< int32_t > Vec3i
Definition: NanoVDB.h:1687
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
GLsizei samples
Definition: glcorearb.h:1298
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.
GLfloat v0
Definition: glcorearb.h:816
#define OPENVDB_NUMERIC_TREE_INSTANTIATE(Function)
Definition: version.h:160
Base class for typed trees.
Definition: Tree.h:36
GLint j
Definition: glad.h:2733
GLsizeiptr size
Definition: glcorearb.h:664
__hostdev__ T dot(const Vec3T &v) const
Definition: NanoVDB.h:1556
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1602
GLuint index
Definition: glcorearb.h:786
SharedPtr< const GridBase > ConstPtr
Definition: Grid.h:81
math::Vec3< Index32 > Vec3I
Definition: Types.h:73
void setSpatialAdaptivity(const GridBase::ConstPtr &grid)
GLfloat GLfloat v1
Definition: glcorearb.h:817
GLuint GLfloat * val
Definition: glcorearb.h:1608
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GA_API const UT_StringHolder N
SharedPtr< const TreeBase > ConstPtr
Definition: Tree.h:40
#define SIZE
Definition: simple.C:41
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:857
Definition: core.h:1131
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE constexpr T abs(T a) IMATH_NOEXCEPT
Definition: ImathFun.h:26
VolumeToMesh(double isovalue=0, double adaptivity=0, bool relaxDisorientedTriangles=true)
GLboolean r
Definition: glcorearb.h:1222
void OIIO_UTIL_API split(string_view str, std::vector< string_view > &result, string_view sep=string_view(), int maxsplit=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
Vec4< int > Vec4i
Definition: NanoVDB.h:1834
bool ValueType
Definition: NanoVDB.h:5729
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
GLint GLsizei count
Definition: glcorearb.h:405
Definition: format.h:895
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
auto join(It begin, Sentinel end, string_view sep) -> join_view< It, Sentinel >
Definition: format.h:2559
Collection of quads and triangles.
Definition: VolumeToMesh.h:102
const openvdb::Vec3I & triangle(size_t n) const
Definition: VolumeToMesh.h:129