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