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