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