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