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